nextnano.com
nextnano³
Download | Search | Copyright | Publications
* password protected

nextnano³ software

numeric-control

Numeric control

With this flag all of the numerical methods and input parameters can be controlled when
- solving strain
- solving the Poisson equation
- solving the nonlinear Poisson equation
- solving the current equation
- solving the coupled Poisson-current equation
- solving the Schrödinger equation (1-band, 6-band k.p, 8-band k.p)
- solving the coupled Schrödinger-Poisson equation
- solving the coupled Schroedinger-Poisson-current equation
- ...

Other options:
- Schrödinger-Poisson predictor-corrector method
- scale Poisson equation
- scale current equation
- factor for separation energy
- evaluation of Fermi functions
- starting value for the electrostatic potential - initial guess
- k.p options
- ...

!---------------------------------------------------------------! $numeric-control optional ! simulation-dimension integer required ! ! strain-lin-eq-solv character optional ! strain-symm-sparse-matrix character optional ! strain-iterations integer optional ! strain-residual double optional ! strain-volume-correction-residual double optional ! strain-volume-correction-iterations integer optional ! filling-factor-for-ILU-decomposition integer optional ! ! discretize-only-once character optional ! ! newton-method character optional ! 1D/2D/3D newton-max-linesearch-steps character optional ! 1D/2D/3D ! nonlinear-poisson-iterations integer optional ! 1D/2D/3D nonlinear-poisson-residual double optional ! 1D/2D/3D nonlinear-poisson-stepmax double optional ! 1D/2D/3D nonlinear-poisson-cg-lin-eq-solv character optional ! 1D/2D/3D nonlinear-poisson-cg-iterations integer optional ! 1D/2D/3D nonlinear-poisson-cg-residual double optional ! 1D/2D/3D ! current-poisson-method character optional ! 1D/2D/3D current-poisson-lin-eq-solv character optional ! 1D/2D/3D current-poisson-lin-eq-solv-iterations integer optional ! 1D/2D/3D current-poisson-lin-eq-solv-residual double optional ! 1D/2D/3D current-poisson-iterations integer optional ! 1D/2D/3D current-poisson-residual double optional ! 1D/2D/3D ! current-block-iterations-Fermi integer optional ! 1D/2D/3D current-block-relaxation-Fermi double optional ! 1D/2D/3D ! coupled-current-poisson-iterations integer optional ! 1D/2D/3D coupled-current-poisson-residual double optional ! 1D/2D/3D coupled-current-poisson-stepmax double optional ! 1D/2D/3D ! current-problem character optional ! 1D/2D/3D current-problem-iterations integer optional ! 1D/2D/3D current-problem-cg-iterations integer optional ! 1D/2D/3D current-problem-cg-residual double optional ! 1D/2D/3D ! schroedinger-1band-ev-solv character optional ! schroedinger-1band-iterations integer optional ! schroedinger-1band-residual double optional ! schroedinger-1band-more-ev integer optional ! ! schroedinger-chearn-el-cutoff double optional ! schroedinger-chearn-hl-cutoff double optional ! ! schroedinger-kp-ev-solv character optional ! schroedinger-kp-iterations integer optional ! schroedinger-kp-residual double optional ! schroedinger-kp-more-ev integer optional ! schroedinger-kp-discretization character optional ! schroedinger-kp-basis character optional ! ! schroedinger-poisson-problem character optional ! 1D/2D/3D schroedinger-poisson-precor-iterations integer optional ! 1D/2D/3D schroedinger-poisson-precor-residual double optional ! 1D/2D/3D ! poisson-boundary-condition-along-x character optional ! 1D/2D/3D poisson-boundary-condition-along-y character optional ! 1D/2D/3D poisson-boundary-condition-along-z character optional ! 1D/2D/3D scale-poisson double optional ! 1D/2D/3D scale-current double optional ! 1D/2D/3D ! fermi-function-mode character optional ! 1D/2D/3D fermi-function-precision double optional ! 1D/2D/3D carrier-statistics character optional ! 1D/2D/3D ! potential-from-function character optional ! new initial-potential double optional ! zero-potential character optional ! built-in-potential-qm character optional ! ! separation-energy-shift double optional ! 1D/2D/3D separation-energy-shift-eV double optional ! 1D/2D/3D ! schroedinger-masses-anisotropic character optional ! 1D/2D/3D ! use-band-gaps character optional ! varshni-parameters-on character optional ! lattice-constants-temp-coeff-on character optional ! ! Luttinger-parameters character optional ! Kane-momentum-matrix-element character optional ! ! valence-band-masses-from-kp character optional ! 1D/2D/3D ! kp-cv-term-symmetrization character optional ! kp-vv-term-symmetrization character optional ! ! 8x8kp-params-from-6x6kp-params character optional ! 8x8kp-params-rescale-S-to character optional ! ! quantization-axis-of-spin character optional ! ! valence-bandedge-energies-zb character optional ! valence-bandedge-energies-wz character optional ! ! strain-transforms-k-vectors character optional ! 1D/2D/3D ! broken-gap character optional ! ! piezo-charge-at-boundaries character optional ! 1D/2D/3D pyro-charge-at-boundaries character optional ! 1D/2D/3D ! piezo-second-order character optional ! piezo-constants-zero character optional ! pyro-constants-zero character optional ! ! discontinous-charge-density-1band character optional ! 1D/2D ! superlattice-option character optional ! superlattice-save-wavefunctions character optional ! ! exchange-correlation character optional ! ! coulomb-matrix-element character optional ! 1D/2D/3D calculate-exciton character optional ! 1D/2D/3D exciton-lambda double optional ! 1D exciton-lambda-step-length double optional ! 1D exciton-iterations integer optional ! 1D/2D/3D exciton-residual double optional ! 1D/2D/3D exciton-electron-state-number integer_array optional ! 1D/2D/3D exciton-hole-state-number integer_array optional ! 1D/2D/3D number-of-electron-states-for-exciton integer optional ! 1D/2D/3D number-of-hole-states-for-exciton integer optional ! 1D/2D/3D ! tight-binding-method character optional ! 1D/2D/3D tight-binding-parameters double_array optional ! 1D/2D/3D tight-binding-calculate-parameter character optional ! 1D/2D/3D ! graphene-potential-fluctuation double_array optional ! ! get-k-vector-dispersion-for-lead-modes character optional ! 1D/2D/3D ! $end_numeric-control optional ! !---------------------------------------------------------------!

Syntax

simulation-dimension = 1 !for 1D simulation = 2 !for 2D simulation = 3 !for 3D simulation

Strain (strain-minimization, 2D/3D only)

Minimization of the elastic energy (i.e. minimization of the integral over the elastic energy density)

strain-lin-eq-solv = BiCGSTAB !
= PCG_Kent !
Switch between
- BiCGSTAB (BiConjugate Gradient Stabilized) algorithm and
- Preconditioned Conjugate Gradient (written by Kent Smith, Agere Systems), ILU as preconditioner
for solving strain equation: solve linear system Ax=b for x.
1D: not implemented
2D: BiCGSTAB/PCG_Kent
3D: BiCGSTAB/PCG_Kent

filling-factor-for-ILU-decomposition = 10 !(default)
This is an optimization parameter for thePCG_Kent algorithm.
Number of filling levels [0,n]
ILU: Incomplete LU factorization of a sparse matrix A
It is used as a preconditioner for solving the system of linear equations A.x = b, where the matrix A is sparse, and A.x=b is solved with BiCGSTAB or another conjugate gradient method.
For large matrices a smaller value than 10 is recommended.
ILU(0) has the same number of offdiagonal matrix elements than the original sparse matrix A.

strain-symm-sparse-matrix = yes ! strain-symm-sparse-matrix = no !(default)Sets up symmetric (upper half of matrix) or nonsymmetric matrix (full matrix) assuming that the matrix is symmetric. The latter option also allows that the matrix is not symmetric.
Obviously, the required storage is twice as much but the matrix-vector products are faster.
If periodic boundary conditions are used for the strain matrix, then in any case a nonsymmetric matrix is assumed.
Also, for free-standing structures, the strain matrix is nonsymmetric.

strain-iterations = 2000 ! -->number of iterations for BiCGSTAB when calculating strain

strain-residual = 1d-12 ! -->residual for strain calculation (BiCGSTAB)

These refer to the following variables incontrol_numeric: -> INTEGER: StrainIterations2D -> REAL: StrainResidual2D -> INTEGER: StrainIterations3D -> REAL: StrainResidual3D

Example:
Dot calculations calculations for strain (4nm_dot.in with finest grid).

Results:

Input file:
strain-iterations = 4000 4000 4000 4000 4000 4000 strain-residual = 1.00d-13 1.00d-12 1.00d-10 1.00d-09 1.00d-08 1.00d-06
Output:
No. of iterations 1356 1221 615 411 67 3 final residual 2.23E-08 ? 4.31E-05 4.18E-04 4.40E-03 0.31E-00 maximum value f. displ. 3.79E-11 3.79E-11 3.77E-11 4.38E-11 6.16E-11 1.76E-11(units m)minimum value f. displ. -8.73E-11 -8.73E-11 -8.75E-11 -8.06E-11 -6.16E-11 -1.55E-11(units m)

So what is interesting is the fact that the final residual is always less than the one specified as the actual residual. But the BiCGSTAB routine leaves without errors of convergence.

As a conclusion the residual should be (at least 1.00d-10) in the input file.

Strain algorithm which also works for free-standing structures

strain-volume-correction-residual = 1d-6 !residual criteria for the volume correction iteration (only relevant forstrain-calculation = strain-minimization-new) strain-volume-correction-iterations = 50 !maximum number of iterations for the volume correction iteration (only relevant forstrain-calculation = strain-minimization-new)

Discretization of equations

Discretize equation twice or use linked list.
The thing is that if one uses lists it takes a lot of memory. So now one can choose between CPU time or RAM.

discretize-only-once = yes ! (default) = no

If 'yes', the matrix elements of the discretized equations are saved into a linked list,
and after that the elements from the list are written into a (sparse) matrix.
So the equation has to be discretized only once.
But this requires more memory (RAM) and can be critical for large simulations.
That's why the user can say 'no' and then the discretization will be done twice,
first to count the number of matrix elements to be stored and second to store them.
Of course, this will require the double time necessary for discretization.

This is related to the strain routine and to the box discretization method of the Schrödinger equations.

Newton's method

newton-method = Newton-1 ! 1D/2D/3Dnewton-method = Newton-2 !1D/2D/3D
newton-method = Newton-3 ! 1D/2D/3D (default)

There are three different Newton method's implemented. 'Newton-1', 'Newton-2', and 'Newton-3'. 'Newton-3' is based on the C++ version of nextnano³.

This will affect the following equations:

Classical nonlinear Poisson equation (nonlinear-poisson-...)
Nonlinear Schrödinger-Poisson equation based on predictor-corrector approach (nonlinear-poisson-...)
schroedinger-poisson-problem = precor ! (predictor-corrector)
Coupled Current-Poisson equation in 1D (for Poisson equation) (nonlinear-poisson-...)
current-poisson-method = couple-all-true = couple-all-false
Coupled Current-Poisson equation in 1D/2D (for current equation) (coupled-current-poisson-)
current-poisson-method = couple-all-true = couple-all-false

Note: For 'Newton-2' and 'Newton-3' there is an additional specifier available which determines the maximum number of steps for the linesearch:

newton-max-linesearch-steps = 20 !1D/2D/3D

Nonlinear Poisson equation (classical or quantum mechanical):

nonlinear-poisson-iterations = 100 !1Dnonlinear-poisson-iterations = 50 !2Dnonlinear-poisson-iterations = 50 !3Dnonlinear-poisson-residual = 1.0d-8 !1Dnonlinear-poisson-residual = 1.0d-12 !2Dnonlinear-poisson-residual = 1.0d-8 !3D

nonlinear-poisson-stepmax = 1.0d-2 !1Dnonlinear-poisson-stepmax = 1.0d-2 !2Dnonlinear-poisson-stepmax = 1.0d-2 !3D

-> NonLinPoissonIterationsXD ... maximum number of iterations
-> NonLinPoissonResidualXD   ... upper limit for residual of functional
-> NonLinPoissonStepMaxXD    ... maximum step length

NonLinPoissonStepMax1D/2D/3D is an input quantity that limits the length of the steps so that you do not try to evaluate the function in regions where it is undefined or subject to overflow (used in subroutine linesearch in subroutinenewton).dtolmin                = 1.0d-2 --> tolmin(in Newton) = tolf*dtolmin ! These values should not be altered.
NonLinPoissonStepMaxXD= 1.0d-2 -->in Newton                        !These values should not be altered.
alf                  = 1.0d-4 -->in Newton                        !These values should not be altered.

not implemented yet (maybe later):
NonLinPoissonIterations_actXD ... Values that are actually used,
NonLinPoissonResidual_actXD   ... can be varied in self-consistent calculation.

stage, delta, Arnoldi

epsilon_ev_rel_to_precision

Conjugate gradient method for nonlinear Poisson equation
Used for calculating one Newton correction step.

nonlinear-poisson-cg-lin-eq-solv = LAPACK             !1D           (default for 1D)                                 = LAPACK-full        !1D                                 = LAPACK-tridiagonal !1D                                 = cg                 !1D/2D/3D                                 = cg-MICCG           !1D/2D/3D                                 = cg-Jacobi          !1D/2D/3D                                 = PCG_Kent           !1D/2D/3D (default for 2D/3D)
- LAPACK             LAPACK routine (banded matrix)
- LAPACK-full        LAPACK routine (full matrix)
- LAPACK-tridiagonal LAPACK routine (tridiagonal matrix) (not implemented yet)
- cg              (conjugate gradient with default preconditioning (Jacobi))
- cg-MICCG        (conjugate gradient with modified incomplete Choleski preconditioning)
- cg-Jacobi       (conjugate gradient with Jacobi preconditioning)
- PCG_Kent           (Preconditioned conjugate gradient written by Kent Smith, faster thancg)

nonlinear-poisson-cg-iterations = 1000 !1Dnonlinear-poisson-cg-iterations = 1000 !2Dnonlinear-poisson-cg-iterations = 1000 !3Dnonlinear-poisson-cg-residual = 1.0d-6 !1D 1.0d-15(classical ?)nonlinear-poisson-cg-residual = 1.0d-6 !2Dnonlinear-poisson-cg-residual = 1.0d-6 !3D

NonLinPoissonCGLinEqSolvXD      ... linear equation solver for nonlinear Poisson equation
NonLinPoissonCGIterationsXD     ... upper limit for number of iterations for conjugate gradient
NonLinPoissonCGResidualXD       ... value of residual

not implemented yet:NonLinPoissonCGIterationsXD_act ...Values that are actually used.
NonLinPoissonCGResidualXD_axt   ... Can be varied in self-consistent solution.

More information can be found here: nonlinear Poisson equation.

Current-Poisson equation

current-poisson-method = block-iterative        !Block-iterative                       = couple-all-true        !Current-Poisson coupled
                                               !(i.e. electrostatic potential, Fermi level electrons, Fermi level holes are coupled)                       = couple-all-false       !Current-Poisson coupled
                                                !(i.e. electrostatic potential, Fermi level electrons, Fermi level holes are not coupled)                       = couple-all-true-aniso !Current-Poisson coupled ('anisotropic' setup of matrix)                                               ! (i.e. electrostatic potential, Fermi level electrons, Fermi level holes are coupled)                       = couple-all-false-aniso !Current-Poisson coupled ('anisotropic' setup of matrix)                                                ! (i.e. electrostatic potential, Fermi level electrons, Fermi level holes are not coupled)1D: couple-all-true  (default)2D: couple-all-true  (default)
3D: block-iterative  (default)

CurrentPoissonMethod1D/2D/3D

It seems that for a classical Double Gate MOSFET simulation, at least in 2D, block-iterativeis much faster than couple-all-false.

current-poisson-lin-eq-solv = BiCGSTAB !(default) = QMRducpl ! = LAPACK !1D = PCG_KENT !2DSwitch between
- BiCGSTAB (BiConjugate Gradient Stabilized) algorithm and
- QMRDUCPL (Quasi-Minimal Residual for Double, Unsymmetric based on CouPLed two-term look-ahead Lanczos)
- PCG_KENT (Preconditioned conjugate gradient written by Kent Smith)
- LAPACK (Lapack routine DGBSV, Lapack only possible for noncoupled equations.)
for solving current-Poisson equation: solve linear system Ax=b for x.
1D: QMRducpl
2D: not implemented
3D: not implemented

Number of iterations and residual for Current-Poisson linear equation solver:

current-poisson-lin-eq-solv-iterations = 4000 !1Dcurrent-poisson-lin-eq-solv-iterations = 4000 !2D

current-poisson-lin-eq-solv-residual = 1.0d-10 !1Dcurrent-poisson-lin-eq-solv-residual = 1.0d-10 !2D

Maximum number of iterations (should be used with care, for some devices 1000 iterations and a residual of 1.0e-16 might be required)

current-poisson-iterations = 25     !1D
current-poisson-iterations = 25     !2D
current-poisson-iterations = 25    !3D

Precision for solution of current and Poisson equation

current-poisson-residual = 1.0e-10 !1Dcurrent-poisson-residual = 1.0e-10 !2Dcurrent-poisson-residual = 1.0e-10 !3D

Block-iterative current-Poisson equation:

In the current equation, for a fixed electrostatic potential the Fermi levels for electrons and holes are solved.
With these new Fermi levels, the mobility, the recombination and generation rates, and the electron and hole densities are calculated anew.
Then the new Fermi levels are calculated until self-consistency is achieved. The maximum number of iterations for this self-consistency loop can be adjusted as well as the relaxation parameter w that relaxes the quasi-Fermi levels in each iteration step k.

E_F,n^(k+1) = w E_F,n + (1 - w) E_F,n^(k) E_F,p^(k+1) = w E_F,p + (1 - w) E_F,p^(k)

Maximum no. of iterations for Fermi levels in block-iterative current equation: current-block-iterations-Fermi = 10 ! 1D current-block-iterations-Fermi = 10 ! 2D current-block-iterations-Fermi = 10 ! 3D current-block-residual-Fermi = 1d-10 ! 1D current-block-residual-Fermi = 1d-10 ! 2D current-block-residual-Fermi = 1d-10 ! 3DRelaxation parameter w: current-block-relaxation-Fermi = 0.2d0 ! 1D current-block-relaxation-Fermi = 0.2d0 ! 2D current-block-relaxation-Fermi = 0.2d0 ! 3D

Note: Choosing current-block-relaxation-Fermi = 1d0 means no relaxation at all.,
choosingcurrent-block-relaxation-Fermi = 0d0 is not really a good idea.

Coupled current-Poisson equation:

coupled-current-poisson-iterations = 100 !1Dcoupled-current-poisson-iterations = 100 !2Dcoupled-current-poisson-iterations = 100 !3Dcoupled-current-poisson-residual = 1.0d-6 !1Dcoupled-current-poisson-residual = 1.0d-6 !2Dcoupled-current-poisson-residual = 1.0d-10 !3D

coupled-current-poisson-stepmax = 1.0d-2 !1Dcoupled-current-poisson-stepmax = 1.0d-2 !2Dcoupled-current-poisson-stepmax = 1.0d-2 !3D

-> CouplCurrentPoissonIter1D/2D/3D    ... maximum number of iterations
-> CouplCurrentPoissonResid1D/2D/3D   ...upper limit for residual of functional
-> CouplCurrentPoissonStepMax1D/2D/3D ... maximum step length

CouplCurrentPoissonStepMax1D/2D/3D is an input quantity that limits the length of the steps so that you do not try to evaluate the function in regions where it is undefined or subject to overflow (used in subroutine linesearch in subroutinenewton).dtolmin                = 1.0d-2 --> tolmin(in Newton) = tolf*dtolmin ! These values should not be altered.
NonLinPoissonStepMaxXD= 1.0d-2 -->in Newton                        !These values should not be altered.
alf                  = 1.0d-4 -->in Newton                        !These values should not be altered.

not implemented yet (maybe later):
CouplCurrentPoissonIterations_act1D/2D/3D ... Values that are actually used,
CouplCurrentPoissonResidual_act1D/2D/3D   ... can be varied in self-consistent calculation.

Current problem (Drift diffusion)

Here one can specify the method of the current calculation.

Note: When solving the current equations one takes the Fermi levels fromMODULE fermi_leveland calculates the relevant densities on that basis. The result is a new Fermi level. The potential is kept fixed throughout the whole procedure.

current-problem               = solve-for-Fermi-integrate-current   !1D only (default)                              = solve-for-Fermi-LAPACK              !tridiagonal LAPACK solver (1D only)                              = solve-for-Fermi-cg                  !conjugate gradient algorithm (default in 2D/3D)Solves drift-diffusion equation by
a) 1D:           integrating the current equation (solve-for-Fermi-integrate-current) or by
(Comment S. Birner: I am not sure ifsolve-for-Fermi-integrate-current works if recombination terms are present, e.g. in the case of pn junctions.)
b) 1D/2D/3D: using a LAPACK (solve-for-Fermi-LAPACK, 1D only) or conjugate gradient algorithm (solve-for-Fermi-cg). Here, the drift-diffusion equation is solved by discretizing for the Fermi levels, while considering the electron and hole densities to be constant in each iteration step.

The electrostatic potential is kept constant until convergence of the current equation is obtained.

current-problem-iterations    = 5 !1D                              = 1 !2D/3D current-problem-residual     = 1d-10This loop is iterated CurrentProblemIterations1D/2D/3D times, its default value is CurrentProblemIterations1D/2D/3D (default: 5 (1D) / 1 (2D/3D)).
current-problem-residualchecks for convergence.
This flag is used in the current routine (SUBROUTINE solve_current_problem) where the quasi-Fermi levels are kept fixed (usually called 'block-iterative').

current-problem-cg-iterations = 1000   !1D                              = 2000   !2D                              = 500   !3DIteration steps for conjugate gradient method for solution for Fermi level.
(Value that is actually used: CurrentProblemCGIterations1D/2D/3D_actual)
Old name was: iter_cg_fermi_default

current-problem-cg-residual   = 1.0d-14Default value for residual is CurrentProblemCGResidual1D/2D/3D = 1.0d-14.
(Value that is actually used: CurrentProblemCGIterations1D/2D/3D_actual)
Old name was: resid_cg_fermiFor programmers: current-problem               = CurrentProblem1D/2D/3D               current-problem-iterations    = CurrentProblemIterations1D/2D/3D                                               CurrentProblemIterations1D/2D/3D_actual: not implemented yet
              current-problem-residual     = CurrentProblemResidual1D/2D/3D               current-problem-cg-iterations = CurrentProblemCGIterations1D/2D/3D                                               CurrentProblemCGIterations1D/2D/3D_actual: not implemented yet              current-problem-cg-residual   = CurrentProblemCGResidual1D/2D/3D                                               CurrentProblemCGResidual1D/2D/3D_actual: not implemented yet
              solve-for-Fermi-integrate-current               solve-for-Fermi-LAPACK               solve-for-Fermi-cg

1-band Schrödinger equation

Eigenvalue solver option for single-band Schrödinger equation (e.g. nonlinear Poisson equation):

schroedinger-1band-ev-solv = LAPACK              !1D                           = LAPACK-DSYEV         !1D                           = LAPACK-DSYEVX        !1D                           = LAPACK-DSTEV         !1D (tridiagonal)                           = LAPACK-DSTEVD        ! 1D (tridiagonal)                           = LAPACK-ZHEEV         ! 1D                           = LAPACK-ZHEEVX        ! 1D                           = LAPACK-ZHEEVR        !1D (should be faster than LAPACK-ZHEEV, and LAPACK-ZHEEVX)                           = LAPACK-ZHBGV         !1D/2D/3D                           = LAPACK-ZHBGVX        ! 1D/2D/3D                           = laband               ! 1D/2D/3D                           = ARPACK               ! 1D/2D/3D                           = ARPACK-shift-invert !1D/2D/3D                           = chearn               ! 1D/2D/3D                           = it_jam               !1D/2D/3D                           = FEASTd               !1D (FEAST-dense)                           = FEAST-CSR            !1D (FEAST-CSR)                           = FEAST-RCI            !1D (FEAST-RCI)LAPACK        = DSYEVX(default if real matrix) or ZHEEVR (default if complex matrix) from LAPACK package (calculates all eigenvalues)
LAPACK-DSYEVX = DSYEVX from LAPACK package (calculates selected eigenvalues thus faster thanLAPACK-DSYEV)
                1D default
                2D/3D not implemented
LAPACK-DSYEV = DSYEV (if real matrix) or ZHEEV (if complex matrix) from LAPACK package (calculates all eigenvalues)
LAPACK-ZHEEV = DSYEV (if real matrix) or ZHEEV (if complex matrix) from LAPACK package (calculates all eigenvalues)
                1D
              2D/3D not implemented
LAPACK-ZHEEVX = ZHEEVX from LAPACK package (calculates selected eigenvalues thus faster thanLAPACK-ZHEEV)
                1D
                2D/3D not implemented
LAPACK-ZHEEVR = ZHEEVR from LAPACK package (calculates selected eigenvalues faster thanLAPACK-ZHEEVX).
                1D (default)
                2D/3D not implemented
LAPACK-DSTEV = DSTEV (if real matrix) from LAPACK package (calculates all eigenvalues of a tridiagonal matrix)
LAPACK-DSTEVD = DSTEVD (if real matrix) from LAPACK package (calculates all eigenvalues of a tridiagonal matrix - divide and conquer algorithm)
laband      = ZHBGV from LAPACK package for banded matrix (1D)
LAPACK-ZHBGV = ZHBGV from LAPACK package for banded matrix (1D)
                ZHBGVX from LAPACK package for banded matrix (2D/3D) forschroedinger-masses-anisotropic = boxonly
LAPACK-ZHBGVX = ZHBGVX from LAPACK package for banded matrix
                2D/3D forschroedinger-masses-anisotropic = boxonly
ARPACK       = Arnoldi from ARPACK package with option              'SA' !smallest (algebraic) eigenvalues (for electrons) (real version)                'LA' ! largest (algebraic) eigenvalues (for holes)       (real version)                'SR' !smallest real parts                     (for electrons) (complex version)                'LR' ! largest real parts                     (for holes)       (complex version)               1D not implemented
               2D default
               3D default
ARPACK-shift-invert =
               Arnoldi from ARPACK package in shift-invert mode with option               1D not implemented
                2D/3D not much tested yet
CheArn        =Chebychev-Arnoldi (Arnoldi form with Chebychev preconditioning)
               (eigenvalue solver for real, symmetric matrices and for hermitian matrices)                ~10 times faster thanARPACKbut this depends a lot on the size of the matrix and on the number of requested eigenvalues, and for setting the correct energy range (cutoff).               You will have to specify a cutoff.
it_jam        = iterative method by Jacek Majewski (not tested sufficiently)
               Only for extreme eigenvalues if inner eigenvalue H² is diagonalized.              1D not implemented
               2D not tested
              3D not tested
FEASTd        = FEAST dense solver
FEAST-CSR     = FEAST CSR solver
FEAST-RCI     = FEAST RCI solver
               FEAST solver package:               E. Polizzi, Density-Matrix-Based Algorithms for Solving Eigenvalue Problems, Phys. Rev. B 79, 115112 (2009)
             Note: FEAST requires an upper limit where to search for eigenvalues, similar to 'chearn'.

In 1D, LAPACK is default (which is fine for small matrices dim < 1000), in 2D/3DARPACK is default.
The fastest results are obtained with chearn provided that the user sets the cutoff energy appropriately.

Eigensolver for Schrödinger equation: ARPACK/Chebychev-Arnoldi

schroedinger-1band-iterations = 1000 !1Dschroedinger-1band-iterations = 1000 !2Dschroedinger-1band-iterations = 1000 !3Dmaximum number of iterations
schroedinger-1band-residual = 1.0d-12 !1Dschroedinger-1band-residual = 1.0d-12 !2Dschroedinger-1band-residual = 1.0d-12 !3D It seems that 1.0d0-12 is better than 1.0d-10 in order to resolve degenerate eigenvalues.Residual

This currently applies only to: ARPACK schroedinger-1band-more-ev = 6 ! 1D schroedinger-1band-more-ev = 6 ! 2D schroedinger-1band-more-ev = 6 ! 3D If n eigenvalues are required, then ARPACK is asked to calculate n + delta_num_ev eigenvalues. This helps to improve convergence of eigenvalues. (Very often some eigenvalues are zero.) This does not apply to any any other solver. epsilon_ev_rel_to_precision --> From outer iteration the actual relative precision for the potential is known. The precision for the eigensolver relative to that value is given by epsilon_ev_relative_to_precision = = 0 --> equal = -1 --> eigenvalues have epsilon precision*1d-10 = 1 --> eigenvalues have epsilon precision*1d10 epsilon_it_arnoldi --> precision which is actually used



Hint: When magnetic field is switched on (only 2D and 3D), the 1-band 
Schrödinger equation is solved with a complex eigenvalue solver.

However, 
if the magnetic field strength is zero, all imaginary entries are zero.

In that case, a real eigenvalue solver would be faster.

Exactly, the same applies for the superlattice option.

You could specify "magnetic-field-on = yes" and "magnetic-field-strength 
= 0.0d0" (keyword $magnetic-field) 
in order to solve the 1-band Schrödinger equation with the complex 
method rather than using the real solvers.
 

Cutoff for Schrödinger equation eigenvalue solver
schroedinger-chearn-el-cutoff = 5d0   
! [eV] 2D/3D

schroedinger-chearn-hl-cutoff = 5d0   ! [eV] 2D/3D

Relative cutoff for Chebychev-Arnoldi eigenvalue solver in units of [eV].

Reference is the 
minimum band edge for electrons and the maximum band edge for holes.

Absolute band edge for electrons: Minimum band edge for electrons + 
schroedinger-chearn-el-cutoff

Absolute band edge for holes:      Maximum band edge 
for holes       - 
schroedinger-chearn-hl-cutoff

Note: This parameter is also necessary for the 
FEAST eigenvalue solver.



Essentially, one has to tell the eigenvalue solver to in which energy interval 
[E_min , E_max] the eigenvalues are sought.

For electrons: [E_min , E_max] = [ E_c,min , E_c,min+ schroedinger-chearn-el-cutoff ]



For holes:   [E_min , E_max] = [ E_v,max- schroedinger-chearn-hl-cutoff , E_v,max ]
 

6-band and 8-band k.p Schrödinger equation

Eigenvalue solver option for 6-band and 8-band k.p Schrödinger equation:

schroedinger-kp-ev-solv = LAPACK ! 1D = LAPACK-ZHEEVR ! 1D = LAPACK-ZHEEV ! 1D = LAPACK-ZHEEVX ! 1D = laband ! 1D/2D/3D = LAPACK-ZHBGV ! 1D/2D/3D = LAPACK-ZHBGVX ! 1D/2D/3D = ARPACK ! 1D/2D/3D = ARPACK-shift-invert ! 2D/3D = chearn ! 1D/2D/3D = it_jam ! 1D/2D/3D = FEASTd ! 1D (FEAST-dense) = FEAST-CSR ! 1D (FEAST-CSR) = FEAST-RCI ! 1D (FEAST-RCI)

LAPACK =use default solver from LAPACK packageLAPACK-ZHEEVR= ZHEEVR from LAPACK package (calculates selected eigenvalues faster than LAPACK-ZHEEV or LAPACK-ZHEEVX) 1D default 2D/3D not implemented LAPACK-ZHEEV = ZHEEV from LAPACK package (calculates all eigenvalues) 1D 2D not implemented 3D not implemented LAPACK-ZHEEVX = ZHEEVX from LAPACK package (calculates selected eigenvalues thus faster than LAPACK-ZHEEV) 1D 2D/3D not implemented laband = ZHBGV from LAPACK package for banded matrix (1D) LAPACK-ZHBGV = ZHBGV from LAPACK package for banded matrix (1D) ZHBGVX from LAPACK package for banded matrix (2D/3D) for schroedinger-kp-discretization = box-integration only LAPACK-ZHBGVX = ZHBGVX from LAPACK package for banded matrix 2D/3D for schroedinger-kp-discretization = box-integration only ARPACK = Arnoldi from ARPACK package with option 'SM' ! smallest magnitude (8-band k.p, for holes and electrons) (complex version) 'LR' ! largest real parts (6-band k.p, for holes) (complex version) ARPACK-shift-invert = Arnoldi from ARPACK package with option 1D not tested yet 2D/3D not tested yet chearn = Chebychev-Arnoldi (Arnoldi form with Chebychev preconditioning) (eigenvalue solver for hermitian matrices) 1D 2D 3D it_jam = iterative method by Jacek Majewski (not tested sufficiently) Only for extreme eigenvalues if inner eigenvalue H² is diagonalized. 1D not tested yet 2D not tested yet 3D not tested yet FEASTd = FEAST dense solver FEAST-CSR = FEAST CSR solver FEAST-RCI = FEAST RCI solver FEAST solver package: E. Polizzi, Density-Matrix-Based Algorithms for Solving Eigenvalue Problems, Phys. Rev. B 79, 115112 (2009) Note: FEAST requires an upper limit where to search for eigenvalues, similar to 'chearn'.



Note, for 6-band k.p in 2D and 3D it is recommended to use 
ARPACK. In 1D, 
LAPACK 
should be used in all cases unless the matrices are larger than dimension = 1000.

The fastest results are obtained with chearn 
provided that the user sets the cutoff energy appropriately. Note that
chearn does not work for 8-band k.p 
where interior eigenvalues have to be calculated.
 

Some suggestions for the parameters:


  it_jam



schroedinger-kp-iterations   =   
      ! 1D

schroedinger-kp-iterations   =  30     
! 2D

schroedinger-kp-iterations   =         
! 3D

 ---> SchroedingerkpIterations1D (2D/3D)

schroedinger-kp-residual     =          ! 1D (value is a guess, needs 
some testing)

schroedinger-kp-residual     =  1.0d-10
! 2D

schroedinger-kp-residual     =         
! 3D

 ---> SchroedingerkpResidual1D (2D/3D)

! itere_max = 30

! eps_itere = 1.0d-10

! itsub_max = 200

! eps_itsub = 1.0d-10



  Arnoldi
  Chebychev-Arnoldi

  => You will need a cutoff for the Schrödinger equation eigenvalue solver 
  (see there for details):

     schroedinger-chearn-el-cutoff = 5d0

  
schroedinger-chearn-hl-cutoff = 5d0


 

 

Eigensolver: arnoldi-complex

for 6-band k.p or 8-band k.p, i.e. Hermitian version

 

SchroedingerkpIterations1D (2D/3D) --> number of iterations

Set to a default value, does not change during simulation.

SchroedingerkpResidual1D (2D/3D)   --> precision for arnoldi method; is 
defined according to precision of outer loop



Note: For extreme eigenvalues --> arnoldi with option 
'LM' or 'SM'

         For inner eigenvalues      
--> arnoldi with option 'SM'



NOTE: All eigenvalues around E_separate are calculated due 
to H --> H-E_separate.

For electrons, however, one only needs eigenvalues above this edge (holes: 
below).

--> Number of eigenvalues is multiplied with a factor 
mult_num_ev_arnoldi and only the relevant eigenvalues are taken. 
mult_num_ev_arnoldi here is just a default value, it is specified in each
k.p region.

If the calculated eigenvalues are not sufficient --> mult_num_ev_arnoldi 
(in the k.p region) is increased by add_num_ev_arnoldi, 
otherwise decreased by dec_num_ev_arnoldi.

If the calculated eigenvalues are not sufficient, H --> 
H-highest_eigenvalue and again the specified number of eigenvalues is 
calculated...



Note: arnoldi seems to have problems with degeneracies:

--> disturb_arnoldi is added to diagonal elements of the upper 
block in order to destroy the symmetry.

 
This currently applies only to: ARPACK

schroedinger-kp-more-ev    =
 6     ! 1D

schroedinger-kp-more-ev    =
 6     ! 2D

schroedinger-kp-more-ev    =
 6     ! 3D



If n eigenvalues are required, then 
ARPACK is asked to
calculate n + delta_num_ev eigenvalues.

This helps to improve convergence of eigenvalues. (Very often some eigenvalues are zero.)

schroedinger-kp-discretization = box-integration       ! 
1D/2D/3D    (default) Here, the output is with respect to heavy 
hole, light hole and split-off hole basis (HH LH SO).

                               
= box-integration-XYZ   ! 1D/2D/3D    
Here, the output is with respect to XYZ basis.

                               
= finite-differences    ! 1D/2D/3D

                               
= finite-differences1D  ! 1D (for testing purposes)

Switch between box integration and finite differences discretization of k.p 
Schrödinger equation (similar to 1-band Schrödinger equation).

We strongly recommend box-integration. Finite-differences seems to be not fully reliable.
For box-integration, the k.p Hamiltonian is set up with respect to the XYZ crystal coordinate system. It is then rotated into the heavy hole, light hole and split-off hole basis, and then the eigenvalues are calculated.
For box-integration-XYZ, the k.p Hamiltonian is set up with respect to the XYZ crystal coordinate system. Then the eigenvalues are calculated.
The eigenvalues and square of the wavefunctions are identical for both methods. The spinor components are different, however.

Schrödinger-Poisson problem

For coupled Schrödinger-Poisson equation.

schroedinger-poisson-problem = precor !(predictor-corrector) = newton !not implemented yet'precor' -->Predictor-corrector approach (default, only option so far)'newton' -->Newton method is used (not implemented yet)

For predictor-corrector method

Number of iterations for Schrödinger-Poisson predictor-corrector cycle:

schroedinger-poisson-precor-iterations = 12 ! 1D schroedinger-poisson-precor-iterations = 12 ! 2D schroedinger-poisson-precor-iterations = 12 ! 3D

Final precision for Schrödinger-Poisson predictor-corrector cycle:

schroedinger-poisson-precor-residual = 1.0d-9 !1D schroedinger-poisson-precor-residual = 1.0d-9 !2D schroedinger-poisson-precor-residual = 1.0d-9 !3D = 1.0d-16 !3D ???

SchroedingerPoissonPreCorIterXD --> default value (default: 12)
SchroedingerPoissonPreCorResidXD --> default value (default: 1.0d-9)
Maybe introduce later a variable for these that change during execution.

Boundary conditions for Poisson equation along x, y and z directions

poisson-boundary-condition-along-x = periodic ! = Neumann !Neumann boundary conditions (default) = Dirichlet !not implemented yetpoisson-boundary-condition-along-y = periodic ! = Neumann !Neumann boundary conditions (default) = Dirichlet !not implemented yetpoisson-boundary-condition-along-z = periodic ! = Neumann !Neumann boundary conditions (default) = Dirichlet !not implemented yet

This feature is useful if one has a superlattice and wants to have periodic boundary condtions for the Poisson equation as well.
It currently only works in 1D and 3D if the Poisson equation is solved over the whole device, i.e. contacts ($poisson-boundary-conditions) are not allowed.

Scale Poisson

scale-poisson = 1.0d0 !1D (default:1.0d0)scale-poisson = 1.0d9 !2D (default: 1.0d9)scale-poisson = 1.0d18 !3D (default: 1.0d18)(Variable in code: ScalePoissonXD = 1.0d0)

In 1D the Poisson problem is scaled by 1.
In 3D and 2D the Poisson problem is scaled, as described in MODULE control_numeric:

Linear Poisson:
------------ Note: Using SI system has the effect that in equation

PoissonM * xV = bV
both sides (i.e. PoissonM and bV) are scaled by a factor about 10^-22 in 3D (2D: 10^-12).
We have noticed that the conjugate gradient method in templates does not work well within that regime.
Therefore both sides are scaled by a factorScalePoisson3D.
To be specific, all elements of PoissonM and bV are scaled AFTER definition of Dirichlet boundary conditions in subroutine scale_poisson.
ScalePoisson3D=1d18 has been proved to be adequate in 3D (2D: 1d9).
This value should not be altered.

Scale current

scale-current = 1.0d0 !1D (default:1.0d0)scale-current = 1.0d0 !2D (default: 1.0d0)scale-current = 1.0d0 !3D (default: 1.0d0)(Variable in code: ScaleCurrentXD = 1.0d0)

Evaluation of Fermi functions (Fermi-Dirac integrals)

Note:Aymerich-Humet andvanHalen-Pulfrey are not implemented any more, soTrellakisis the default.Goanois possible, however.

fermi-function-mode = Trellakis ! ='4' Trellakis = Goano ! ='3' Goano: Uses algorithm of Goano = vanHalen-Pulfrey ! ='2' van Halen, Pulfrey: Uses approximations as documented inMODULE fermi_functions. = Aymerich-Humet ! ='1' Aymerich-Humet FermiFunctionMode1D = '4' ! 'Trellakis' !1D (default)FermiFunctionMode2D = '4' ! 'Trellakis' !2D (default)FermiFunctionMode3D = '4' ! 'Trellakis' !3D (default)fermi-function-precision = 1.0d-6 ! 1Dfermi-function-precision = 1.0d-6 !2D
fermi-function-precision = 1.0d-6 ! 3D

FermiFunctionPrecision1D/2D/3D --> only valid for 'Goano'
In case no approximation exists, 'Goano' is taken with FermiFunctionPrecision1D/2D/3D.

Not implemented yet:
Given are default values, they may be overwritten during execution. Defaults which are not changed during execution:
- FermiFunctionMode1D/2D/3D_default -FermiFunctionPrecision1D/2D/3D_default

The following Fermi-Dirac integrals are used inside the code:
- F_3/2:- Goano(if(ABS(x) > 353.6d0) .AND. (ABS(x) < 354.2d0) Trellakisis used) (The calculation forABS(x)from 353.6 up to 354.2 fails withGoanoon SUN.) - Trellakis- F_1/2:- Goano(if(ABS(x) > 353.6d0) .AND. (ABS(x) < 354.2d0) Trellakisis used) (The calculation forABS(x)from 353.6 up to 354.2 fails withGoanoon SUN.) - Trellakis- F_-1/2:- Goano(if(ABS(x) > 353.6d0) .AND. (ABS(x) < 354.2d0) Trellakisis used) (The calculation forABS(x)from 353.6 up to 354.2 fails withGoanoon SUN.) - Trellakis- F_-3/2:- Goano(if(ABS(x) > 353.6d0) .AND. (ABS(x) < 354.2d0) Trellakisis used) (The calculation forABS(x)from 353.6 up to 354.2 fails withGoanoon SUN.) - Trellakis

These Fermi-Dirac integrals include the Gamma prefactors of the Fermi-Dirac integral.

Carrier statistics for classical densities

carrier-statistics = Fermi-Dirac !Classical densities are calculated using Fermi-Dirac statistics (default) = Maxwell-Boltzmann !Classical densities are calculated using Maxwell-Boltzmann statistics (not recommended as this is only an approximation which is only applicable in certain cases)

n = N_c * FermiIntegral_1/2( (E_F - E_c) / (k_BT) ) !electron density forFermi-Dirac p = N_c * FermiIntegral_1/2( (E_v - E_F) / (k_BT) ) ! hole density for Fermi-Dirac

n = N_c * exp( (E_F - E_c) / (k_BT) ) !electron density for Maxwell-Boltzmann p = N_c * exp( (E_v - E_F) / (k_BT) ) ! hole density for Maxwell-Boltzmann

User-defined potential profile (new)

Using this specifier, one can define a customized function that defines the electrostatic potential phi.
This electrostatic potential is subtracted from the conduction and valence band edge profiles.
E_c' = E_c - e phi
E_v' = E_v - e phi

potential-from-function = "2 * x + y" ! phi(x,y,z) = ... [V] = no !do not use user-defined electrostatic potential

The variables x, y, z refer to the grid point coordinates of the simulation area in units of [nm].

Examples:

potential-from-function = " - 1/2 (x-15)^2 " !phi(x,y,z) = ...: parabolic potential shifted by 15 nm along the positive x axis
%a = 0.5 potential-from-function = "%a * sin(x)" !phi(x,y,z) = ...: sinus like potential using%aas a variable

The following operators and functions are supported:
+ , - , * , / , ^ abs , exp , sqrt , log , log10 , sin , cos , tan , sinh , cosh , tanh , asin , acos , atan

If you use potential-from-function = "...", you typically want to calculate the wave functions in a particular potential landscape.
For this, additionally, you need.
$simulation-flow-control flow-scheme = 3 !solve Schrödinger equation only

Starting value for the potential - initial guess

This might be helpful in rare cases to improve the convergence of the initial Poisson equation.

initial-potential = 0.0 ! [V]

By default, this flag is ignored. An initial guess for the electrostatic potential phi(x) is calculated automatically taking into account doping properties.
If this flag is present, then instead, the initial guess for the electrostatic potential is set to phi(x) = initial-potential.

If you read in raw data (e.g. the potential) from a sweep-step or not, then initial-potential is ignored in any case because you read in a previously calculated potential rather than calculating it from scratch.

If you start a simulation, first the built-in potential is calculated for the equilibrium, called built-in potential, then voltage is applied with the built-in potential as starting value to determine the potential in nonequilibrium.
So the initial-potential is just relevant (as a starting guess) to calculate the built-in potential, after that it doesn't have any influence.
The units are in [V] and should refer somehow to the values (with respect to magnitude) in the database (conduction-band-energies, valence-band-energies in units of [eV])
because it holds for the resulting conduction band edgeE_c(x) = E_c,0(x) - phi(x). The band edge is calculated taking into account that the Fermi level is located at 0 eV.

Skip calculation of (classical) built-in potential

zero-potential = yes                       ! set electrostatic potential phi to 0 [V]
                = no                       !(default is 'no'), initial guess of phi for each grid point based on doping (if present) and local charge neutrality in undoped regions
                = use-mid-gap-initial-guess !             initial guess of phi for each grid point based on doping (if present) and mid band-gap position of Fermi level in undoped regions

Using zero-potential = yes omits the first calculation of the Poisson equation (which is solved without quantum mechanics, i.e. classically) which is used to obtain an initial guess for the electrostatic potential.
This is useful for e.g.
- to improve speed (e.g. if no doping or charge carriers are present)
- One can set the electrostatic potential to zero and then calculate the eigenstates withflow-scheme =3, i.e. if one is interested in solving the Schrödinger equation only once and without self-consistency.
- For testing and debugging purposes, e.g. testing the eigenvalue solvers.

The potential of this initial guess and the built-in potential can be written out in$output-bandstructure (built-in-potential = yes).
potential_built_in.dat potential_built_in_cl.dat potential_built_in_qm.dat potential_initial_guess.dat

Calculation of built-in potential quantum mechanically

Calculate built-in potential quantum mechanically directly after it has been calculated classically.

built-in-potential-qm = yes !
= no !(default is 'no')

Calculation of interior eigenvalues

The following only applies to 8-band k.p where the interior eigenvalues of the energy spectrum are sought.
(For single-band and 6-band k.p one is calculating the eigenvalues at the end of the energy spectrum which is much easier.)

In almost all cases, one is only interested in the energies around the band gap, i.e. only a few electron levels and only a few hole eigenvalues are needed.
Therefore, one is using fast eigenvalue solvers such as ARPACK that pick out the relevant eigenvalues around the band gap.

The eigenvalue solver then needs the information at which energy the eigenvalues shall be calculated.

In most cases, we have the information where in the device the conduction band edge E_c(x) has a minimum and where the valence band edge E_v(x) has a maximum.
Note that the position x might not be identical for MIN(E_c) and MAX(E_v).
So we can tell the eigenvalue solver to search for energies around the middle of the band gap, i.e. E_separate = (E_c - E_v) / 2.
Then all eigenvalues that are larger than E_separate are electron eigenvalues and those that are smaller in energy are hole eigenvalues.

The actual value of E_separate influences the CPU time for finding the eigenvalues.
Note that the Hamiltonian is solved twice. First, one determines the electron eigenvalues, and then the hole eigenvalues.
(In contrast, the nextnano++ software solves the Hamiltonian only once.)

1. Factor for separation energy

separation-energy-shift = 0.3d0 ! must be within range [0,1]

Specifies where between conduction and valence bands, i.e. typically within the band gap, the separation energy E_separate should be located.
It determines the distance of the separation energy relative to the conduction and valence band,
e.g. SeparationEnergyShift = 0.1 means thatE_separateis (0.1 * band gap) below the lowest conduction band energy (or above the highest valence band for holes).

CASE('electrons') E_separate = E_cond - SeparationEnergyShift * (E_cond - E_val) !for electrons

CASE('holes') E_separate = E_val + SeparationEnergyShift * (E_cond - E_val) !hor holes

This energy is used for the shift-inverse ARPACK eigenvalue solver.

(It could happen that the hole states are assigned to electron states if the separation energy shift is wrong, or vice versa.)

Note: The implementation is different for type-II heterostructures.

2. Separation energy

separation-energy-shift-eV = 0.0d0 ! [eV]

SeparationEnergyShift_eV ==> E_separate in units of[eV]

Separation energy between conduction and valence band states on an absolute energy scale,
e.g. separation-energy-shift-eV = 0.1d0 means that E_separate is at 0.1 eV.

This energy is used for the shift-inverse ARPACK eigenvalue solver.

(It could happen that the hole states are assigned to electron states if the separation energy is wrong, or vice versa.)

Essentially, this parameter tells the eigenvalue solver where to look for eigenvalues, e.g.
- electron eigenvalues (if$quantum-model-electronsis used) are sought aboveseparation-energy-shift-eV [eV]- hole eigenvalues (if$quantum-model-holesis used) are sought below separation-energy-shift-eV [eV]

Note: The flag separation-energy-shift-eV might be useful for type-II heterostructures where electron and hole eigenvalues might overlap.
Note: If separation-energy-shift-eV is not equal to zero, then this separation energy will be used.
If it is equal to zero, then the other model (described above) using a factor for the separation energy is used (separation-energy-shift).

Schrödinger masses isotropic/anisotropic

schroedinger-masses-anisotropic = yes ! 1D/2D/3D
= box ! 1D/2D/3D = no ! 1D/2D/3D (CHECK: Are the masses treated correctly at material interfaces?) (default for 2D/3D) = 1D ! 1D discretization of the Schrödinger equation (default for 1D)

 
The effective mass tensor can be written out with:

$output-1-band-schroedinger

 effective-mass-tensor = yes

If effective mass is isotropic (at Gamma point: conduction band and 
heavy, light and split-off hole band) then the mass tensor is diagonal in any 
coordinate system. So one needs only pure derivatives for discretization (dx²,dy²,dz²). 
If mass is anisotropic (L point, X point), one needs also mixed 
derivatives (dx dy, dx dz, dy dz).
The mass tensor can become anisotropic when strain is applied. In the case of 
[001] quasi-homogeneous strain, mass tensor will be anisotropic, but still 
diagonal in the calculation system, so one doesn't need schroedinger-masses-anisotropic = 
yes.

In the case of arbitrary grown strained structures hole effective masses 
are nondiagonal in the calculation system. For this we need schroedinger-masses-anisotropic = 
yes.
box 
means that it is anisotropic but box integration (finite volume) 
is carried out instead of finite differences. Note that 
Neumann boundary conditions are not included yet for option "box", 
only Dirichlet is possible. (The 
boundary conditions for the Schrödinger equation are set within the keywords
$quantum-model-electrons/$quantum-model-holes.)

Restrictions:

schroedinger-masses-anisotropic = 
yes does 
not support periodic boundary conditions so far.


 
		 
Use band gaps or conduction band energies


use-band-gaps = yes  ! Use band-gaps and bow-band-gaps.

              = no  
 ! Use conduction-band-energies 
and bow-conduction-band-energies. (default)

In the database file or in the input file one can either specify

  - conduction-band-energies (and bow-conduction-band-energies) 
(default) or

  - band-gaps (and bow-band-gaps).

This flag specifies which option should be used.
 
 

Varshni parameters


varshni-parameters-on = yes  ! 
Temperature dependent energy gaps.

                      
= no   ! 
Band gaps independent of temperature. 
Absolute values from database are taken.

Varshni parameters are used to determine temperature dependent energy gaps 
(i.e. temperature dependent conduction-band-energies).


More information...

 

 

Temperature dependent lattice constants


lattice-constants-temp-coeff-on = yes  ! 
Temperature dependent lattice constants.

                                
= no   ! 
Lattice constants independent of temperature. 
Absolute values from database are taken.

The coefficients are used to determine temperature dependent 
lattice-constants.

More information...

 
 
Choice of k.p parameters
There are two sets of k.p parameters that are used in the literature 
to describe the bulk k.p energy dispersion E(k_x,k_y,k_z) 
in a semiconductor:

	Luttinger parameters:                              gamma₁, gamma₂, gamma₃ 
	and kappa  (kappa is often not specified, it can be approximated from 
	the other three parameters.)
	Dresselhaus-Kip-Kittel (DKK) parameters: L, M, N                               
	and kappa  (kappa is often not specified, it can be approximated from 
	the other three parameters and is related to N⁺, N^- 
	where N = N⁺ + N^-.)

 
For 6-band k.p for holes, the following parameters are relevant where 
kappa is optional:

	Luttinger parameters: gamma₁, gamma₂, gamma₃, 
	Delta_so, (kappa)
	DKK parameters:       L, M, N,                               Delta_so, (kappa)

For 8-band k.p for electrons and holes, the following parameters are 
relevant:

	Luttinger parameters: gamma₁', gamma₂', gamma₃', 
	Delta_so, (kappa'), E_gap, E_p, S ( = 1 + 2 
	F), B
	DKK parameters:       L', M'=M, N',                          Delta_so, (kappa'), E_gap, E_p, S ( = 1 + 2 
	F), B

Note the primes on the 8-band k.p parameters which distinguish them 
from the unprimed 6-band k.p parameters.
 
nextnano³ provides the user full flexibility to input the k.p 
parameters.

	For 6-band k.p, the user can either use the the LMN parameters 
	(default) or the Luttinger parameters gamma₁, gamma₂, gamma₃ 
	(Luttinger-parameters = 
	6x6kp).

	For both choices, the kappa parameter can be used or not. In the latter case 
	(default), it is calculated automatically.
	The 8-band k.p parameters depend on the band gap E_gap.

	For 8-band k.p, the user can either use the values listed in 
	8x8kp-parameters (default), or calculate the 8-band k.p 
	parameters from the 6-band k.p parameters (recommended).

	In this way, the temperature dependent band gap E_gap is 
	automatically taken into account.

	This is important, as L', N', kappa' (or gamma₁', gamma₂', gamma₃', 
	kappa') and S depend on the band gap.

	The B parameter is used in any case.
	The user can either use the S parameter specified in 
	8x8kp-parameters (default) or calculate this parameter from the 
	conduction band electron mass (recommended).
	Kane's momentum matrix element E_p is usually given in units 
	of [eV] (default) or [eV Angstrom].

	The user can use his preferred choice using

	Kane-momentum-matrix-element = E_P  ! E_P 
in units of [eV]         
	 (default)
                              
	P    ! P in units of [eV Angstrom]

	This choice affects the value of E_p in 8x8kp-parameters.

The actually used k.p parameters can be found in these files:

- material_parameters/8x8kp_Luttinger_parameters_used.dat

- material_parameters/8x8kp_parameters_used.dat

Use Luttinger parameters instead of Dresselhaus-Kip-Kittel (DKK) k.p parameters


The Luttinger parameters are called gamma₁, gamma₂ and 
gamma₃ and are typically given for a 6-band k.p 
Hamiltonian.

There are additional Luttinger parameters called kappa and q.
q is not implemented yet.

Luttinger-parameters = 6x6kp (or) yes ! Use database entries Luttinger-parameters = gamma₁ gamma₂ gamma₃ ---.
= 6x6kp-kappa ! Use database entries Luttinger-parameters = gamma₁ gamma₂ gamma₃ kappa. = 6x6kp-kappa-only ! Use database entries Luttinger-parameters = --- --- --- kappa, and L, M, N (or L', M', N') parameters. = 8x8kp ! Use database entriesLuttinger-parameters = gamma₁' gamma₂' gamma₃' ---. Note: Here, modified (i.e. 8-band k.p) Luttinger parameters have to be specified. = 8x8kp-kappa ! Use database entries Luttinger-parameters = gamma₁' gamma₂' gamma₃' kappa'. Note: Here, modified (i.e. 8-band k.p) Luttinger parameters have to be specified. = 8x8kp-kappa-only ! Use database entriesLuttinger-parameters = --- --- --- kappa', and L, M, N (or L', M', N') parameters. Note: Here, the modified (i.e. 8-band k.p) kappa parameter has to be specified. = no ! default, i.e. use L, M, N (or L', M', N') parameters.

This flag can be used to take from the database (or input file) the Luttinger parameters (gamma₁, gamma₂, gamma₃, i.e.Luttinger-parameters, and optionally also kappa) instead of the DKK notation (L,M,N, i.e.6x6kp-parametersand 8x8kp-parameters). This flag currently only works for zinc blende.

Example: For 6-band k.p, the following options are possible:
- DKK 6-band k.p parameters L,M,N (default)                                       Luttinger-parameters = no
- DKK 6-band k.p parameters L,M,N including kappa                 Luttinger-parameters = 6x6kp-kappa-only
- Luttinger 6-band k.p parameters gamma1,gamma2,gamma3             Luttinger-parameters = 6x6kp
- Luttinger 6-band k.p parameters gamma1,gamma2,gamma3,kappa   Luttinger-parameters = 6x6kp-kappa

In the default database, the Luttinger parameters are defined for 6-band k.p. i.e. not for 8-band k.p.
The 8-band k.p DKK parameters L', M', N' can be calculated from the Luttinger parameters if the following flags are chosen:
8x8kp-params-from-6x6kp-params = LMNS(or)yes ! Calculate L',M'=M,N',S. Note: S is calculated from E_gap, Delta_so, E_P and m_e.                                = LMN          ! Calculate L',M'=M,N'. Note: S is taken from database entry.For 8-band k.p one can alternatively specify the modified Luttinger parameters gamma₁', gamma₂', gamma₃' directly ifLuttinger-parameters = 8x8kpor =8x8kp-kappa is used.

In wurtzite, the "Luttinger" parameters are called Rashba-Sheka-Pikus parameters. This is used by default in wurtzite.

Summary of the options:

!----------------------------------------------------------------------------------------------- !We will now treat the different options that could be present in the input file, i.e.! ==> 8x8kp-params-from-6x6kp-params = no ! o use 8-band DKK k.p parameters L',M',N' (default) ==> Luttinger-parameters = no ! o use 8-band DKK k.p parameters L',M',N' and kappa' ==> Luttinger-parameters = 8x8kp-kappa-only ! o use 8-band Luttinger k.p parameters gamma1',gamma2',gamma3' ==> Luttinger-parameters = 8x8kp ! o use 8-band Luttinger k.p parameters gamma1',gamma2',gamma3',kappa' ==> Luttinger-parameters = 8x8kp-kappa ! ! ==> 8x8kp-params-from-6x6kp-params = LMN (or LMNS) ! o calculate 8-band DKK k.p parameters from 6-band parameters L,M,N ==> Luttinger-parameters = no ! o calculate 8-band DKK k.p parameters from 6-band parameters L,M,N including kappa ==> Luttinger-parameters = 6x6kp-kappa-only ! o calculate 8-band DKK k.p parameters from 6-band parameters gamma1,gamma2,gamma3 ==> Luttinger-parameters = 6x6kp ! o calculate 8-band DKK k.p parameters from 6-band parameters gamma1,gamma2,gamma3,kappa ==> Luttinger-parameters = 6x6kp-kappa !-----------------------------------------------------------------------------------------------

Calculate `8x8kp-parameters` from `6x6kp-parameters`

8x8kp-params-from-6x6kp-params = no             ! (default)                                = LMNS(or)yes    !Calculate L',M'=M,N',S: Note: S is calculated from E_gap, Delta_so, E_P and m_e.
                              = LMN           ! Calculate L',M'=M,N'. Note: S is taken from database entry.
                              = S              !S is calculated from E_gap, Delta_so, E_P and m_e.
no: Don't calculate 8-band k.p parameters from 6-band k.p parameters. Take 8x8kp-parameters from database or input file (L',M'=M,N',S).

LMNS: Calculate 8-band k.p parameters from 6-band k.p parameters (L',M'=M,N',S). Don't take8x8kp-parametersfrom database or input file (apart from B, E_P).
Note: S is calculated from E_gap, Delta_so, E_P and m_e.

LMN: Calculate 8-band k.p parameters from 6-band k.p parameters (L',M'=M,N'). Don't take 8x8kp-parameters from database or input file (apart from B, E_P, S).
Note: S is taken from database and is not calculated.

S: S is calculated from E_gap, Delta_so, E_P and m_e.
Don't take8x8kp-parametersfrom database or input file (apart from L',M'=M,N',B, E_P).

These flags will overwrite 8x8kp-parameters (L',M',N',S or L',M',N') from database or input file by calculating them from 6x6kp-parameters (L,M,N,E_gap,Delta_so,E_P,m_e).

Note: This flag currently only applies to zinc blende and wurtzite with quantum-model "8x8kp".

Rescale 8-band k.p parameters in order to avoid spurious solutions

Sometimes it might be necessary to rescale the 8-band k.p parameters in order to avoid spurious solutions, see eq. (3.158) and eq. (3.159) in the PhD thesis of S. Birner.

8x8kp-params-rescale-S-to = ONE ! 1D/2D/3D = ZERO ! 1D/2D/3D
= no ! 1D/2D/3D(default)

no: Don't rescale 8-band k.p parameters (default).
ONE: Rescale 8-band k.p parameters so that S=1. This affects S,E_P,L',N',N⁺',N^-' that are calculated anew.
ZERO: Rescale 8-band k.p parameters so that S=0. This affects S,E_P,L',N',N⁺',N^-' that are calculated anew.

This flag will overwrite the 8-band k.p parameters (L',N',N⁺',N^-',E_P,S) by rescaling them to fit S=1 or S=0. This is sometimes necessary to get rid of spurious 8-band k.p eigenfunctions.

Obviously, rescaling so that S=1/S=0 affects the eigenvalues of the electrons more than the eigenvalues of the holes.

Option ZERO: In this case, spurious k.p solutions are avoided in any case (but this argument only holds for the k space and not for a real space discretization with grid spacing Delta_xwhere k_max = 2pi /Delta_x.
                     This corresponds to the situations that remote-band contributions cancel the free-electron term.
                      (S=0: suggestion by B. Foreman)
                      It is likely that one still gets spurious solutions in the conduction band.
Option ONE:   In this case, spurious k.p solutions are avoided in most cases. This corresponds to entirely neglecting remote bands.
                     This is the default setting in the nextnano++ software.

Kane's momentum matrix element

Kane-momentum-matrix-element = E_P ! E_P in units of [eV] (default) P ! P in units of [eV Angstrom]

Refers to E_P parameter specified in 8x8kp-parameters in input file or database.
All parameters in the database are E_P. If the users prefers P instead, he is welcome to do so by using this flag.
This works for both zinc blende and wurtzite.

Valence band effective masses calculated from k.p parameters

valence-band-masses-from-kp = yes ! 1D/2D/3D = no ! 1D/2D/3D

The valence band effective masses for heavy hole, light hole and split-off hole are calculated from the L, M and N parameters (Dresselhaus notation) given in the database for each material (6x6kp-parameters) considering strain effects.
Through this, one has the option to calculate effective masses of holes, depending on k.p parameters and strain and then compare single-band Schrödinger results ('quantum-model-holes = effective-mass' and using effective masses derived from k.p and strain) with "6x6kp" calculation for instance.

Nonsymmetrized or symmetrized k.p Hamiltonian

kp-cv-term-symmetrization = no ! 1D/2D/3D(default) = yes ! 1D/2D/3D kp-vv-term-symmetrization = no ! 1D/2D/3D(default) = yes ! 1D/2D/3D

cv is related to the H_cv part of the k.p Hamiltonian.vvis related to the H_vv part of the k.p Hamiltonian.

Here one can choose between two methods of discretizing the k.p Hamiltonian, a symmetrized and a nonsymmetrized k.p Hamiltonian. The symmetrized Hamiltonian is based on the bulk k.p Hamiltonian, whereas the nonsymmetrized version is based on Burt's k.p formulism for heterostructures (see also Foreman, Elimination of spurious solutions from eight-band kp theory, PRB 56, R12748 (1997) or M.G. Burt, Fundamentals of envelope function theory for electronic states and photonic modes in nanostructures, J. Phys.: Condens. Matter 11, R53 (1999)).
Note: These two Hamiltonians are different only for heterostructures. For bulk structures with infinite barriers (e.g. GaN or Si nanowires), they are identical. Only the nonsymmetrized Hamiltonian corresponds to the correct physics for heterostructures.

Quantization axis of spin

quantization-axis-of-spin = x ! = y ! = z ! = default !Typically, in the physics literature, the spin is quantized with respect to thezdirection.
In a numerical calculation using the nextnano software, the user might want to define a different direction, e.g. if the quantum well is grown along the x, y or z direction.
By default, the program uses adefaultdirection for the spin quantization depending on symmetry and/or growth direction.
This flag is relevant for the k.p Hamiltonian where the basis is defined with respect to a particular axis of spin quantization.
Depending on the choice of quantization-axis-of-spin, the respective unitary transformation is performed.

Currently, this flag is only used within box-integration routine but not within finite-differences.
For finite differences, spinors are rotated afterwards (post-processing), seeSUBROUTINE rotate_spinor_calculation2sim.
In 1D, a simulation can be along the x, y or z axis of the simulation coordinate system.
We consider the quantization axis of spin to be parallel to the simulation axis for a 1D simulation.
This makes sense as the labeling of heavy hole and light hole for a quantum well only makes sense
if the quantization axis of spin is parallel to the growth direction which is the simulation direction.
This is here taken into account automatically.

However, the situation is different for a quantum well simulated in a 2D or 3D domain.
In this case, the quantum well can be oriented along the x, y or z axis.
The program does not know this.
A possible solution could be to use the growth direction, 'growth-coordinate-axis', as a means to determine the quantization axis or angular momentum.
So far, it is the z axis for a 2D and 3D simulation.
In a 1D simulation, this z axis is adjusted automatically to x or y, if the simulation is along x or along y.

(CHECK: How about rotations? Crystal vs. simulation vs. calculation coordinate system?)

Basis used for k.p Hamiltonian

(Preliminary feature for testing. Use with care. If you need more information, please contact us!)

For the k.p Hamiltonian, different bases can be used. For details, see e.g. the book of L. C. Lew Yan Voon, M. Willatzen, The k.p Method (2009)[VoonWillatzen](e.g. Section 3.3.1.1).
The differences are due to zinc blende vs. wurtzite symmetry and due to the phase factors which can be defined differently (see Table C.2).

This flag only makes sense for schroedinger-kp-discretization = box-integration.

schroedinger-kp-basis = default                               ! Cho--Chuang (for zinc blende and box-integration) / Chuang-Chang_u1u2u4u5u3u6 (for wurtzite and box-integration)
(for zinc blende)    = Voon-Willatzen--Bastard--Foreman                          ! L. C. Lew Yan Voon, M. Willatzen, The k.p method (2009);
                                                                 G. Bastard, Wave Mechanics Applied to Semiconductor Heterostructures (1988); B. A. Foreman, PRB 48, 4964 (1993)
   = Cho--Chuang                ! default (for box-integration) ! K. Cho, in 'Excitons, Topics in Current Physics', Vol. 14, Chap. 2 (1979); S. L. Chuang, Physics of Optoelectronic Devices, p. 138
   = Elliot                                                    ! R. J. Elliott, PR 96, 266 (1954)
   = Luttinger-Kohn--Pidgeon-Brown--Bir-Pikus                   ! J. M. Luttinger, W. Kohn, PR 97, 869 (1955); C. R. Pidgeon, R.N. Brown PR 146, 575 (1966); G. L. Bir, G.E. Pikus (1975)
   = Dresselhaus-Kip-Kittel                                  ! G. Dresselhaus, A. F. Kip, C. Kittel, PR 98, 368 (1955)
   = Roth-Lax-Zwerdling                                       ! L. M. Roth, B. Lax, S. Zwerdling, PR 114, 90 (1959)
   = Kane_1956                                                ! E. O. Kane, J. Phys. Chem. Sol. 1, 82 (1956)
   = Kane_1966--Pollak                                        ! E. O. Kane, in 'Physics of III-V Compounds, Semiconductors and Semimetals', Vol. 1, Chapt. 3 (1966);
                                                                 F. H. Pollak, in 'Strained-layer Superlattices: Physics, Semiconductors, and Semimetals', Vol. 32, Chap. 2 (1990)
   = Kane_1982                                                  ! E. O. Kane, in 'Handbook of Semiconductors', Vol. 1 (1982)
   = Koster-Dimmock-Wheeler-Statz--Weiler-Aggrawal-Lax--Fishman ! G. F. Koster, J. O. Dimmock, R. G. Wheeler, H. Statz, Properties of the Thirty-Two Point Groups (1963);
                                                                 M. H. Weiler, R. L. Aggrawal, B. Lax, PRB 17, 3269 (1978); G. Fishman, Energie et Fonction D'onde des Semiconducteurs (1988)
   = Broido-Sham--Winkler                                       ! D. Broido, L. J. Sham, PRB 31, 888 (1985); R. Winkler, Spin-Orbit Coupling in Two-Dimensional Electron and Hole Systems (2003)
   = Sercel-Vahala                                              ! P. C. Sercel, K. J. Vahala, PRB 42, 3690 (1990)
   = deSterke                                                   ! C. M. de Sterke, PRB 36, 6574 (1987) to do: value for phi is not clear
   = Los-Fasolino-Catellani                                     ! J. Los, A. Fasolino, A. Catellani, PRB 53, 4630 (1996), PhD thesis Birner eq. (3.85), PhD thesis Hackenbuchner eq. (2.31)-(2.36)) (for wurtzite)   = Chuang-Chang_u1u2u4u5u3u6 ! default (for box-integration) ! S. L. Chuang, C. S. Chang, PRB 54, 2491 (1996) (eq. (4))   = Chuang-Chang                                               ! S. L. Chuang, C. S. Chang, PRB 54, 2491 (1996) (eq. (4))   = Identity                                                   ! The identity matrix is used, i.e. keep SXYZ basis. (for testing purposes, has the same effect as box-integration-XYZ)   = TEST                                                       ! (user-defined)

The cartesian (S,X,Y,Z) and the angular momentum states (S,|J,m>) basis sets are related by a unitary transformation, see also Section 3.3.1.1 in [VoonWillatzen].
- One can define a 6x6 matrix that transforms from |X+>,|Y+>,|Z+>,|X->,|Y->,|Z-> basis to the angular momentum |J,m> basis.
- One can define an 8x8 matrix that transforms from |S+>,|S->,|X+>,|Y+>,|Z+>,|X->,|Y->,|Z-> basis to the angular momentum |J,m> basis.

We assume the z axis as the quantization axis for the angular momentum. (This can be changed, seequantization-axis-of-spin = x/y/z/default.)

When is this basis used?

The nextnano software defines the k.p Hamiltonian and its parameters always with respect to the SXYZ basis where the spin quantization axis is assumed to be along the cartesian z direction.
There are several rotations involved, e.g. if the crystal is rotated with respect to the simulation coordinate system or if the simulation is performed along the x, y or z direction.

The k.p matrix in the SXYZ basis can be transformed to another angular momentum basis before the eigenvalues are calculated.

schroedinger-kp-discretization = box-integration !Transform SXYZ Hamiltonian into angular momentum basis before the eigenvalues are calculated. Eigenfunctions are then defined with respect to basis specified in schroedinger-kp-basis. = box-integration-XYZ !Use SXYZ Hamiltonian. Eigenfunctions are defined with respect to SXYZ basis.The numercial speed of the iterative eigenvalue solver depends on the basis used.Once the eigenfunctions are known, they can be transformed to any other basis of interest.

to do: CHECK subroutine Examine_kp_State_BasisXYZ



Calculate valence band edge energies

There are several algorithms implemented to calculate the valence band edge 
energies that are used as input for the

- classical densities for the holes and

- single-band Schrödinger equations for 
the holes.

Here, for each grid point, the strain-dependent 6-band k.p 
Hamiltonian is solved at k = (k_x,k_y,k_z) 
= 0 (Gamma point).

The resulting eigenvalues are known as heavy hole, light hole and splif-off 
hole. They represent the valence band edge energies for the single-band 
Schrödinger equation for each hole band.

The challenge is to determine which eigenstate is light hole and which eigenstate is crystal-field hole for wurtzite. Therefore, we implemented several 
methods to test each method.
Zinc blende

valence-bandedge-energies-zb = default    ! (default) Method implemented by M. Povolotskyi for zincblende.
                             = old        !              Method implemented by M. Sabathil for zincblende.
                           = nextnano++ ! (not implemented yet)



Wurtzite

valence-bandedge-energies-wz = default    ! (default) Method implemented by M. Povolotskyi for zincblende, generalized by S. Birner for wurtzite.
                             = old        !              Method implemented by M. Sabathil for zincblende, generalized by S. Birner for wurtzite.
                             = old-zb     !              Method implemented by M. Sabathil for zincblende and using the basis for zincblende (even in the wurtzite case) - not recommended
                           = nextnano++ ! (not implemented yet)



 
Transform crystal momentum because of strain

 strain-transforms-k-vectors = yes 
! 1D/2D/3D

                             
= no   ! 1D/2D/3D
Ref. P. Enders et al., PRB 51, 16695

"k.p theory of energy bands, wave functions, and optical selection rules 
in strained tetrahedral semiconductors"



The idea is the following:

If the elastic deformation is given by a relation r' = (1 + eps)r

then one has to perform a transformation k' = (1 - eps)k.

eps is du/dr and not just the symmetric part of it.

The flag is to decide if we want to do this transformation or not.



LOGICAL :: StrainTransforms_k_Vectors = .FALSE.  ! (default)

                                      
= .TRUE.
Example by Michael Povolotskyi, University of Rome "Tor Vergata" 
(PhD Thesis):

The effect is small. And this is good, because if it were big that would 
mean that the method cannot be applied.

(5 nm 211 InAs/GaAs quantum well, 8-band k.p calculation)




 
Broken gap type-II band alignment (for 8-band k.p calculations)

broken-gap = no ! (default) = yes ! = full-band-density ! FB-EFA: Either$quantum-model-electronsor$quantum-model-holeshas to be specified. = full-band-density-electrons ! FB-EFA: The k.p Hamiltonian is solved only once. $quantum-model-electronshas to be specified. = full-band-density-holes ! FB-EFA: The k.p Hamiltonian is solved only once. $quantum-model-holes has to be specified.
= full-band-density-electrons-only ! FB-EFA: The k.p Hamiltonian is solved only once. Both$quantum-model-electronsand$quantum-model-holesshould be specified. = full-band-density-holes-only ! FB-EFA: The k.p Hamiltonian is solved only once. Both$quantum-model-electronsand$quantum-model-holesshould be specified.



FB-EFA = full-band envelope-function approach



The following error for 8-band k.p

Error update_kp: Cannot handle conduction band below valence band (E_cond_min 
< E_val_max).

can be avoided when using broken-gap =
yes or broken-gap =
full-band-density.

This error appears when the lowest conduction band edge value is lower than the 
highest valence band edge, i.e. a negative band gap is present.

In this case the code does not know how to distinguish between electron and hole 
states in an 8-band k.p approach.

(Don't try to calculate the density when using broken-gap =
yes.)

The density can be calculated using broken-gap = 
full-band-density.

This approach is termed FB-EFA (full-band envelope-function approach).

It is documented in detail in the following two papers:

	T. Andlauer, T. Zibold, P. Vogl

	Proc. SPIE 7222, 722211 (2009)
	Full-band envelope function approach for type-II broken-gap 
	superlattices

	T. Andlauer, P. Vogl

	Physical Review B 80, 035304 (2009)

 
Piezo and pyroelectric charges at the boundaries
 piezo-charge-at-boundaries =
no    ! 1D/2D/3D (default)

                            
=
yes   ! 1D/2D/3D

 pyro-charge-at-boundaries  =
no    ! 1D/2D/3D (default)

                            
=
yes   ! 1D/2D/3D



To allow for piezo and pyroelectric charges at the boundaries. Usually 
nextnano³ sets the piezo and pyroelectric charges to zero at the 
boundaries of the device.



If one assumes that one only simulates a small part of the device (as is the 
case for most situations) one does not have an interface at the boundary (and 
therefore no charge).

If one wants to have an explicit interface at the boundary (because the device 
has a boundary there), one has to define a material-air interface in the input 
file. Alternatively one can use this option.

The corresponding charges can be looked up in the output file 
interface_densitiesD.txt ($output-densities).

This option might be useful for wurtzite (pyroelectric charges) and strained 
zinc blende structures grown along other directions than [100] or [011] where 
there is a nonzero piezoelectric polarization P at the boundary. In vacuum 
P = 
0. Therefore, there is div P nonequal to zero at the surface.



 
Set piezo and pyroelectric constants to zero
 piezo-constants-zero =
no    ! (default)

                      
=
yes   !

 pyro-constants-zero  =
no    ! (default)

                      
=
yes   !



To set the piezo any pyroelectric constants (piezo-electric-constants, 
pyro-polarization) from the database (database_nn3.in) 
or input file to zero. This option is useful to study the effects of piezo and 
pyroelectric charges in wurtzite independently. It is also useful for strained 
structures in order to understand the influence of strain much better.
 
Piezoelectric electric constant: Second order
 piezo-second-order 
=
no                  ! (default)

                    
= 2nd-order           !

                    
= 2nd-order-Tse-Pal   !

identical to 2nd-order for 
zincblende but different for wurtzite. B₁₅₆ for zinc 
blende is zero in Tse-Pal paper.

                    
= 4th-order-Tse-Pal   ! Tse-Pal model and Tse-Pal material parameters (4^th order for 
zinc blende but only 2^nd order for wurtzite because in wurtzite 4^th 
order is not implemented.)
 

See e.g. - First- and second-order piezoelectricity in III-V 
semiconductors

         [BeyaWakata2011] A. Beya-Wakata, P.-Y. Prodhomme, G. Bester, Phys. Rev. B 84, 195207 (2011)

       - [Bester2006]     G. Bester et al., Phys. Rev. Lett. 96, 187602 (2006) 
(zinc blende)

       - [Patra2017]      S. K. Patra, S. Schulz, Phys. 
Rev. B 96, 
155307 (2017) (wurtzite)

       - [Pedesseau2012]  L. Pedesseau et al., Appl. Phys. Lett.
100, 
031903 (2012) (wurtzite)
Here, one can switch on the second order effect of piezoelectricity. 
Obviously several piezoelectric parameters are needed.

  - zinc blende: e₁₄              
(1^st order effect); B₁₁₄  B₁₂₄  B₁₅₆ 
(2^nd order effect)

  - wurtzite:      e₃₃  e₃₁  
e₁₅ (1^st order effect); B₃₁₁  B₃₁₂  B₃₁₃  B₃₃₃  B₁₁₅  B₁₂₅  B₁₃₅  B₃₄₄ (2^nd order effect)
Note: In the simulated structure, a mixture of zinc blende and wurtzite 
materials is allowed.

	2nd-order

	This is the model supported by the database_nn3.in 
	file, i.e. these parameters can be specified in the database or in the input 
	file.

	- For zinc blende materials 
	the 2^nd order model of [BeyaWakata2011] is used.

	  4 parameters are needed.

	  e₁₄  B₁₁₄  B₁₂₄  B₁₅₆   
	(default)

	  The equations for the second order terms read:

	  P_x = 2 * B₁₁₄ * eps_xx * eps_yz + 2 * 
B₁₂₄ * eps_yz * (eps_yy + eps_zz) + 4 * B₁₅₆ 
* eps_xz * eps_xy
  P_x = 2 * B₁₁₄ * eps_yy * eps_xz + 2 * 
B₁₂₄ * eps_xz * (eps_zz + eps_xx) + 4 * B₁₅₆ 
* eps_yz * eps_xy
  P_z = 2 * B₁₁₄ * eps_zz * eps_xy + 2 * 
B₁₂₄ * eps_xy * (eps_xx + eps_yy) + 4 * B₁₅₆ 
* eps_yz * eps_xz


- For wurtzite 
	materials, the 2^nd order model of [Patra2017] is used.

	  11 parameters are needed.

	  e₃₃  e₃₁  
e₁₅  B₃₁₁  B₃₁₂  B₃₁₃  B₃₃₃  B₁₁₅  B₁₂₅  B₁₃₅  B₃₄₄   
	(default)  The equations for the second order terms 
	are given in eq. (3) in [Patra2017].
	2nd-order-Tse-Pal

	- For 
	zinc blende materials it is identical to 2nd-order.

	  Additionally, it is 
	identical to 4th-order-Tse-Pal if one neglects 3^rd and 4^th order terms.

	  4 parameters are needed.

	  e₁₄  B₁₁₄  B₁₂₄  B₁₅₆   
	(default)

	  Note that the Tse-Pal model 
	assumes B₁₁₄ = e₁₁₄, B₁₂₄ = e₁₂₄ 
	and B₁₅₆ = e₁₅₆ = 0.

	- For wurtzite 
	materials, the 2^nd order model of J. Pal, G. Tse et al., Phys. 
	Rev. B 84, 085211 (2011) is used.

	  6 parameters are needed. 
	See also H. Y. S. Al-Zahrani et al., Nano Energy 2, 
1214 (2013).

	  e₃₃ e₃₁ 
	e₁₅  e₃₁₁ e₃₁₃ e₃₃₃  The wurtzite material 
	parameters for this model can only be specified in the input file but not in 
	the database file database_nn3.in.

	    The implemented equations 
	are only valid for growth along the [0001] direction.

 
	4th-order-Tse-Pal

	The material parameters for this model can only be specified 
	in the input file but not in the database file database_nn3.in.

	- For zinc blende materials 
	the model of G. Tse, J. Pal et al., J. Appl. Phys. 114, 073515 (2013) 
	is used.

	  It includes linear terms, 2^nd, 
	3^rd and 4^th order terms for zinc blende.

	  13 parameters are needed.

	  e₁₄  B₁₁₄  B₁₂₄  B₁₅₆  e₁₁₁₄ e₁₁₂₄ e₁₂₂₄ e₁₂₃₄  
	e₁₁₁₁₄ e₁₁₁₂₄ e₁₂₂₂₄ e₁₂₂₃₄ e₁₁₂₃₄  Note that the Tse-Pal 
	model assumes B₁₁₄ = e₁₁₄, B₁₂₄ = e₁₂₄ 
	and B₁₅₆ = e₁₅₆ = 0.

	- For wurtzite 
	materials it is identical to 2nd-order-Tse-Pal.



 
Discontinous quantum charge density at material interfaces in 1D and 2D
For the single-band Schrödinger equation in 1D and 2D, we need a parallel 
effective mass to calculate the quantum mechanical charge density. At material 
interfaces, this parallel effective mass is not continous (because the materials 
have different effective masses) which will lead to a discontinuity in the 
quantum mechanical charge density.

By default, we use an averaged parallel effective mass value which is weighted 
by the parallel effective masses at each grid point and the square of the 
wave function psi² at each grid point (FUNCTION

get_mass_par_averaged). This will lead to a contious quantum 
mechanical charge density.

However, especially in 2D this could lead to an enormous increase in CPU time, 
especially if lots of eigenvalues are specified. Therefore the user can allow 
for discontinous parallel effective masses (and thus a discontinous quantum 
mechanical charge density) to speed up the calculations.
 discontinous-charge-density-1band =
no    ! 1D/2D (default)

                                   
=
yes   ! 1D/2D
The figure on the left shows the quantum mechanical charge density of a 1D 
quantum well for the case of discontinous-charge-density-1band =
yes.

The discontinuity in the parallel effective masses m(z) at material interfaces 
leads to a discontinuity in the quantum charge density at interfaces.



The right figure shows (for a different quantum well!) a comparison of the two 
approaches. The red line uses an averaged 
parallel effective m_o mass which leads to a continous quantum charge 
density (discontinous-charge-density-1band =
no).

The x axis corresponds to the position in nm, the y axis to the charge density 
in |e| * 1 * 10¹⁸ cm^-3.




 
Superlattice
 superlattice-option             
= whole                 
 ! (default)

                                 
= 100-only              
 ! for dispersion along [100] direction, i.e. from Gamma to X

                                 
= 010-only              
 ! for dispersion along [010] direction, i.e. from Gamma to X

                                 
= 001-only              
 ! for dispersion along [001] direction, i.e. from Gamma to X

                                 
= 110-only              
 ! for dispersion along [110] direction, i.e. from Gamma to K

                                 
= 101-only              
 ! for dispersion along [101] direction, i.e. from Gamma to K

                                 
= 001-only              
 ! for dispersion along [001] direction, i.e. from Gamma to K

                                 
= 111-only              
 ! for dispersion along [111] direction, i.e. from Gamma to L

 Here one can specify if one wants to calculate the superlattice 
dispersion for the whole superlattice Brillouin zone,

  or along special directions in the superlattice Brillouin zone.

 Note that the directions refer to the simulation coordinate 
system and not to the crystal coordinate system.

  Internally, only positive K_SL values are 
calculated.



 superlattice-save-wavefunctions = yes                   
 ! (default)

                                 
= no                    
 ! 

 no: Program will not save 
wave functions for K_SL vectors.

 Usually one is interested only in the eigenvalues and not so 
much in the wave functions itself.

 This option is necessary to avoid huge memory allocation, 
e.g. especially in 3D for option whole.
 
Exchange-Correlation (XC)
 exchange-correlation =
no                      ! (default)

                      =
LDA                     ! LDA   = local        density 
approximation

                      =
LSDA                    ! LSDA = local spin density approximation

                      =
LDA-exchange            !             
local        density approximation (only 
exchange,  for testing)

                      =
LSDA-exchange           !             
local spin density approximation (only exchange,  for testing)

                      =
LDA-correlation         !             
local        density approximation (only 
correlation, for testing)

                      =
LSDA-correlation         !             
local spin density approximation (only correlation, for testing)

 
                     =
LDA-not-selfconsistent  !             
local        density approximation (for testing) 
- The exchange-correlation potential energies are written out but are not used 
in the calculation.

                      =
LSDA-not-selfconsistent !             
local spin density approximation (for testing) - The exchange-correlation 
potential energies are written out but are not used in the calculation.

 
Calculates exchange and correlation potentials as a functional of the density 
in LDA or LSDA.



Parameters for Ceperley-Alder model from J.P. Perdew and A. Zunger, Phys. Rev. B
23, 5048 (1981)

Depending on the type of polarization ([ unpolarized (LDA) , polarized (LSDA) 
]), different sets of parameters have to be used.



        [ 
unpolarized , polarized
]

gamma = [  -0.1423d0  , -0.0843d0 ]

beta1 = [   1.0529d0  ,  1.3981d0 ]

beta2 = [   0.3334d0  ,  0.2611d0 ]

A     = [   0.0311d0  ,  0.0155d0 ]

B     = [  -0.048d0   , -0.0269d0 ]

C     = [   0.0020d0  ,  0.0007d0 ]

D     = [  -0.0116d0  , -0.0048d0 ]
 
Exciton

 coulomb-matrix-element                 
= yes

                                        
= no

 To output Coulomb energy (= Hartree 
energy) only (does not work yet).
 calculate-exciton                      
= yes  ! 1D/3D

                                        
= no   ! 1D/3D

 3D: To calculate exciton energy of electron and heavy hole state 
self-consistently (Hartree approximation only).

        It is also possible to calculate the 
exciton energy of electron and light hole, or electron and split-off hole.

 
 exciton-iterations                     
= 5    ! 3D for self-consistency loop (Hartree only)

 Note: For self-consistent Kohn-Sham eigenvalue problem, 
exiton-iterations must be set too a much higher value (!).

 There are two situations where exciton-iterations
is used for self-consistent Kohn-Sham eigenvalue problem:

  - SUBROUTINE calc_exciton             
(=> exciton-iterations =  50)

  - SUBROUTINE calc_many_body_state_LSD (=> exciton-iterations = 
150)



 exciton-residual                       
= 1d-5 ! 3D for self-consistency loop

 There are two further situations where exciton-residual
is used, namely for self-consistent Kohn-Sham eigenvalue problem:

  - SUBROUTINE calc_exciton             
(=> exciton-residual = 1d-5)

  - SUBROUTINE calc_many_body_state_LSD (=> exciton-residual = 1d-6)

exciton-electron-state-number = 1 ! electron ground state = 2 ! 2^nd electron state = 1 3 ! electron ground state and 3^rd state exciton-hole-state-number = 1 ! hole ground state = 2 ! 2^nd hole state = 2 3 ! 2^nd and 3^rd hole state = ... For example with these integer numbers it is possible to specify: - Calculate exciton energy of the electron ground state with the hole ground state. - Calculate exciton energy of the 2^ndelectron state with the hole ground state. - Calculate exciton energy of the electron ground state with the 1^st hole state, and then with the 2^ndhole state. - Calculate exciton energy of the 2^ndelectron state with the 2^ndhole state. Note: The numbering of k.p states is twice as many as single-band because spin up and spin down states are included. Example: The 1^st k.p hole state corresponds to the 1^st single-band holes state. The 3^rd k.p hole state corresponds to the 2^nd single-band holes state. The 5^rd k.p hole state corresponds to the 3^nd single-band holes state.

number-of-electron-states-for-exciton = 3 ! 2D number-of-hole-states-for-exciton = 3 ! 2D

In 1D, these two specifiers are not used.In 2D, these two specifiers specify the number of electron and hole states that are considered for the wave function expansion, i.e. 3 means 1 electron ground state and 2 excited electron states and similar for the holes.In 3D, these two specifiers are used to distiguish between exciton (1 electron + 1 hole), biexciton (2 electrons + 2 holes) and charged excitons (X+: 1 electron + 2 holes, or X-: 2 electrons + 1 hole).

3D: exciton: X0 ----------- number-of-electron-states-for-exciton = 1 ! electron ground state number-of-hole-states-for-exciton = 1 ! hole ground state biexciton: XX0 -------------- number-of-electron-states-for-exciton = 2 ! 2 lowest electron states number-of-hole-states-for-exciton = 2 ! 2 highest hole states Not implemented yet: charged exciton: X- ------------------- number-of-electron-states-for-exciton = 2 ! 2 lowest electron states number-of-hole-states-for-exciton = 1 ! hole ground state charged exciton: X+ ------------------- number-of-electron-states-for-exciton = 1 ! electron ground state number-of-hole-states-for-exciton = 2 ! 2 highest hole states Note: The numbering of k.p states is twice as many as single-band because spin up and spin down states are included. Example: The 1^st k.p hole state corresponds to the 1^st single-band holes state. The 3^rd k.p hole state corresponds to the 2^nd single-band holes state. The 5^rd k.p hole state corresponds to the 3^nd single-band holes state.


In 1D we currently have a routine that calculates the 
exciton correction in type-I quantum wells (single-band only).

See tutorial "Exciton 
energy in quantum wells":

The 1D output can be found here: Schroedinger_1band/exciton_energy1D.dat

It contains information on the exciton binding energy (in units of [meV] 
and [Rydberg]) and on the exciton Bohr radius (which is the variational 
parameter) in units of [nm].



To speed up the calculations, the user can provide an educated guess for the 
start value of the variational parameter lambda.

Both of the following inputs are optional;

 exciton-lambda             
= 2.0   ! [nm] 1D only (optional)

 exciton-lambda-step-length = 0.05  ! [nm] 1D only 
(optional)

If they are not present, then internally calculated values are taken:

 ==> exciton-lambda             
=  BohrRadius2DLimit

 ==> exciton-lambda-step-length = (BohrRadiusBulkExciton - 
BohrRadius2DLimit) / (lambda_steps - 1)






In 3D we can calculate the exciton correction in quantum dots taking into 
account the Coulomb potential (Hartree approximation).

This works for single-band Schrödinger equation for both electrons and holes, 
as well as for 8-band k.p for both electrons and holes. One can 
also treat the electrons within single-band and the holes within 6-band k.p.
The following output files are provided:



Exciton ground state



    Schroedinger_1band/exciton_energy3D_cb001_vb001.dat 
(At present, also for 
k.p calculations the output directory is the single-band directory.)

                                     
vb001 (= heavy hole)

                                     
vb002 (= light hole)

                                     
vb003 (= split-off hole)



    E field [V/m],   E_ex [eV],      
E_el - E_hl,    E_el0 - E_hl0,   Delta_Ex,       REAL(inter_matV(1))

    ...

    0.000000E+000    0.788262E+000   
0.778266E+000   0.797671E+000    0.940847E-002   
0.992412E+000 



E field [V/m]:      
Electric field value (not fully implemented and tested yet)

                     
E field = voltage-offset / lever-arm-length

                    
For details, see 
$output-1-band-schroedinger.

E_ex [eV]:           
Energy of the exciton state

                     
E_ex = E_el - E_hl - 0.5d0 * (mat_el - mat_hl)  (Eq. (6.37), p. 79 
in diploma thesis of Matthias Sabathil)

                    
matrix element mat_el = < psi_hl | V_psi_el | psi_hl >

                    
matrix element mat_hl = < psi_el | V_psi_hl | psi_el >

E_el  - E_hl  [eV]:  Energy of the electron state minus 
the hole state

E_el0 - E_hl0 [eV]: 
Transition energy of the original electron state minus the original hole 
state (without exciton correction), i.e.

                    
energy difference between non-interacting electron and hole 
one-particle states.

Delta_Ex [eV]:      
Energy difference due to exciton correction, i.e. between the energy 
difference of the eigenvalues of the

                    
non-interacting electron and hole one-particle states (E_el0 - 
E_hl0) and the energy of the exciton state E_ex.

                     
Delta_Ex = E_el0 - E_hl0 - E_ex 
(Eq. (6.38), p. 79 in diploma thesis of Matthias Sabathil)

REAL(inter_matV(1)):
single-band case: | SUM(psi_el * psi_hl) |²

                    
k.p: optical matrix element





Biexciton ground state



    Schroedinger_1band/biexciton_energy3D_cb001_vb001.dat 
(At present, also for 
k.p calculations the output directory is the single-band directory.)

                                       
vb001 (= heavy hole)

                                       
vb002 (= light hole)

                                       
vb003 (= split-off hole)



    E field [V/m],   E_biex [eV],    E_el - E_hl,    E_el0 - E_hl0,   Delta_Ex,       REAL(inter_matV(1))

    ...

    0.000000E+000    0.254603E+001   
0.128976E+001   0.135468E+001    0.163326E+000   
0.000000E+000 



E field [V/m]:      
Electric field value (not fully implemented and tested yet)

                     
E field = voltage-offset / lever-arm-length

                    
For details, see 
$output-1-band-schroedinger.

E_biex [eV]:         
Energy of the biexciton state

                     
E_biex = 2 (E_el - E_hl) - 4 * 0.5d0 * (mat_el - mat_hl) - mat_elel - 
mat_hlhl (Eq. on p. 65 
in PhD thesis of O. Stier, TU Berlin)

                    
matrix element mat_el   =   < psi_hl | V_psi_el | psi_hl >

                    
matrix element mat_hl   =   < psi_el | V_psi_hl | psi_el >

                    
matrix element mat_elel =   < psi_el | V_psi_el | psi_el >

                    
matrix element mat_hlhl = - < psi_hl | V_psi_hl | psi_hl >

E_el  - E_hl  [eV]:  Energy of the electron state minus 
the hole state

E_el0 - E_hl0 [eV]: 
Transition energy of the original electron state minus the original hole 
state (without biexciton correction), i.e.

                    
energy difference between non-interacting electron and hole 
one-particle states.

Delta_Ex [eV]:      
Energy difference due to biexciton correction, i.e. between the energy 
difference of the eigenvalues of the

                    
non-interacting electron and hole one-particle states 2*(E_el0 - 
E_hl0) and the energy of the biexciton state E_biex.

                     
Delta_Ex = 2*(E_el0 - E_hl0) - E_biex

REAL(inter_matV(1)):
single-band case: | SUM(psi_el * psi_hl) |²

                    

k.p: optical matrix element





Note: There might be problems with the current nextnano³ implementation 
of the exciton routines for electrons at other valleys than Gamma (due to 
degeneracies, splitting of bands) or for the light holes (single-band 
approximation).
 
Tight-binding: Method and parameters

tight-binding-method = bulk-zincblende-sp3s*-nn     
!     nn =        
nearest-neighbor (20x20 matrix, including spin)

                     
= bulk-zincblende-sp3-nn       
!     nn = 

       nearest-neighbor

(16x16 matrix, including spin)

                     
= bulk-graphene-Saito-nn       
!     nn = 

       nearest-neighbor

                     
= bulk-graphene-Scholz-nn      !
    nn = 

       nearest-neighbor

                     
= bulk-graphene-Scholz-3rd-nn  !
3rd-nn = 

third-nearest-neighbor

                     
= bulk-graphene-Reich-3rd-nn   !
3rd-nn = 

third-nearest-neighbor

The units are [eV]. 16 parameters are required for zinc blende 
materials, even in the case of bulk-zincblende-sp3-nn. 
In this case, the s* related parameters must be zero.


tight-binding-parameters = -3.53284 ! 
Esa         (GaAs)

                           
0.27772 ! Epa

                          
-8.11499 ! Esc

                          
 4.57341 ! Epc

                         
 12.33930 ! Es_a

                           
4.31241 ! Es_c

                        
  -6.87653 ! Vss

                          
 1.33572 ! Vxx

                         
  5.07596 ! Vxy

                          
 0.0     ! Vs_s_

                           
2.85929 ! Vsa_pc

                          
11.09774 ! Vsc_pa

                           
6.31619 ! Vs_a_pc

                           
5.02335 ! Vs_c_pa

                           
0.32703 ! Delta_so_a

                           
0.12000 ! Delta_so_c
Note: Special case for graphene (3 parameters only, nearest-neighbor 
tight-binding approximation):
 tight-binding-parameters =  0.0   
 ! [eV] E_2p: site energy of the 2p_z 
atomic orbital (orbital energy)

                           
-3.013   ! [eV] gamma₀: C-C 
transfer energy (nearest-neighbor, nn)

                           
 0.129   ! []   s₀: 
denotes the overlap of the electronic wave function on adjacent sites (nn)
Note: Special case for graphene (7 parameters only, 
third-nearest-neighbor tight-binding approximation):
 tight-binding-parameters = -0.28    ! 
[eV] E_2p: site energy of the 2p_z atomic orbital 
(orbital energy)

                           
-2.97    ! [eV] gamma₀: C-C 
transfer energy (nearest-neighbor, nn)

                           
 0.073   ! []   s₀: 
denotes the overlap of the electronic wave function on adjacent sites (nn)

                           
-0.073   ! [eV] gamma₁: (2^nd-nn)

                           
 0.018   ! []   s₁:          
(2^nd-nn)

                           
-0.33    ! [eV] gamma₂: (3^rd-nn)

                           
 0.026   ! []   s₂:          
(3^rd-nn)
Special option for third-nearest-neighbor tight-binding approximation in 
graphene:
tight-binding-calculate-parameter = no   
 ! (default) use 7 
parameters (E_2p,gamma₀,gamma₁,gamma₂,s₀,s₁,s₂)

                                  
= E2p    !             
use 6 parameters (       gamma₀,gamma₁,gamma₂,s₀,s₁,s₂)

                                  
= gamma1 !             
use 6 parameters (E_2p,gamma₀,             gamma₂,s₀,s₁,s₂)

E_2p or gamma₁ can be calculated internally in order to 
force E(K) = 0 eV where K is the K point in the Brillouin zone. In that case, 
only 6 parameters are used (although 7 parameters have to be present in the 
input file). The parameter to be calculated is simply ignored inside the code.
 
Graphene potential fluctuation
graphene-potential-fluctuation = 0.0                   
0.0001                 
4.0       
 ! [eV] [eV] [] (default)

                               
= 0.0                   
0.0001        
                    ! 
[eV] [eV]

                               
= 0.0                                                     
! [eV] (default, i.e. "ideal" graphene layer)

                               
= 0.100                 
                                  ! 
[eV]

                               
= 0.100                 
0.001                  
3.0
        ! [eV] [eV] [] These 
are useful values, especially the latter two for CPU efficiency.

                               
= PotentialFluctuation  EnergyGridResolution  
sigma_factor ! [eV] [eV] []
For Dirac point potential fluctuation see H. Xu et al,. APL 98, 133122 
(2011).

We assume a Gaussian distribution of the potential fluctuation.

The 1^st value of these array entries corresponds to the fluctuation 
in potential energy of the Dirac point in units of [eV]. Useful 
values lie in the range from 0 eV to 0.200 eV.

The 2^nd value of these array entries (optional) corresponds to the 
energy grid resolution of the Gaussian integration in units of [eV].

The 3^rd value of these array entries (optional) is related to the 
integration range of the Gaussian integration in units of [].

The integration range is from [ - sigma_factor
*  PotentialFluctuation 
, + sigma_factor * 
PotentialFluctuation ] [eV].

Adjusting the 2^nd and 3^rd value has implications on CPU 
time. 

3.0 
is a very reasonable value for 
sigma_factor 

and also an energy grid resolution of 1 meV (0.001) 
is sufficient.


 
CBR method
 get-k-vector-dispersion-for-lead-modes = 
arccos ! (default)

                                        
= use_ln ! 

                                        
= acosIN ! use intrinisc arccosine function 
DACOS (only for testing)

 
In principle, all methods for arcos(x) = cos^-1(x) should be 
identical. (acosIN  only for certain 
energy range where k is real.)
At present, it seems that use_ln 
works for electrons only.                 
(for holes: 'mass' is negative)

                 
It seems that arccos works for both, 
electrons and holes. (for holes: 'mass' is negative)

	k = 1 / dx * cos^-1 ( ( E - E_l,m ) / 2 t - 1 )

For more details, see Matthias Sabathil's documentation of CBR_example, 
page 4 ("Lead modes").

Numeric control

Syntax

Strain (strain-minimization, 2D/3D only)

Discretization of equations

Newton's method

Nonlinear Poisson equation (classical or quantum mechanical):

Current-Poisson equation

Block-iterative current-Poisson equation:

Coupled current-Poisson equation:

Current problem (Drift diffusion)

1-band Schrödinger equation

6-band and 8-band k.p Schrödinger equation

Schrödinger-Poisson problem

For predictor-corrector method

Boundary conditions for Poisson equation along x, y and z directions

Scale Poisson

Scale current

Evaluation of Fermi functions (Fermi-Dirac integrals)

Carrier statistics for classical densities

User-defined potential profile (new)

Starting value for the potential - initial guess

Skip calculation of (classical) built-in potential

Calculation of built-in potential quantum mechanically

Calculation of interior eigenvalues

1. Factor for separation energy

2. Separation energy

Schrödinger masses isotropic/anisotropic

Use band gaps or conduction band energies

Varshni parameters

Temperature dependent lattice constants

Choice of k.p parameters

Use Luttinger parameters instead of Dresselhaus-Kip-Kittel (DKK) k.p parameters

Calculate 8x8kp-parameters from 6x6kp-parameters

Rescale 8-band k.p parameters in order to avoid spurious solutions

Kane's momentum matrix element

Valence band effective masses calculated from k.p parameters

Nonsymmetrized or symmetrized k.p Hamiltonian

Quantization axis of spin

Basis used for k.p Hamiltonian

Calculate valence band edge energies

Transform crystal momentum because of strain

Broken gap type-II band alignment (for 8-band k.p calculations)

Piezo and pyroelectric charges at the boundaries

Set piezo and pyroelectric constants to zero

Piezoelectric electric constant: Second order

Discontinous quantum charge density at material interfaces in 1D and 2D

Superlattice

Exchange-Correlation (XC)

Exciton

Tight-binding: Method and parameters

Graphene potential fluctuation

CBR method

Calculate `8x8kp-parameters` from `6x6kp-parameters`