[ Home ] [ News ] [ Contact ] [ Search ] [ nextnano.com ] numeric-control

 Download Software
 Copyright notice
 Publications

 GUI: nextnanomat
 Tool: nextnano³
 Tool: nextnano++
 Tool: nextnano.QCL
 * password protected

 

 
Up
 

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-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                    integer      optional ! 1D/2D/3D
 fermi-function-precision               double       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 !
                                                              !
 calculate-bandedge-energies            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      optional ! 1D/2D/3D
 exciton-hole-state-number              integer      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 the PCG_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 in control_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 for strain-calculation = strain-minimization-new)
 strain-volume-correction-iterations = 50   !
maximum number of iterations for the volume correction iteration (only relevant for strain-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/3D
newton-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      ! 1D
nonlinear-poisson-iterations     =   50      !
2D
nonlinear-poisson-iterations     =   50      !
3D

nonlinear-poisson-residual       =  1.0d-8   !
1D
nonlinear-poisson-residual       =  1.0d-12  !
2D
nonlinear-poisson-residual       =  1.0d-8   !
3D

nonlinear-poisson-stepmax        =  1.0d-2   ! 1D
nonlinear-poisson-stepmax        =  1.0d-2   !
2D
nonlinear-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 subroutine newton).

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 than cg)

nonlinear-poisson-cg-iterations =  1000     ! 1D
nonlinear-poisson-cg-iterations =  1000     !
2D
nonlinear-poisson-cg-iterations =  1000     !
3D

nonlinear-poisson-cg-residual   =  1.0d-6   !
1D   1.0d-15 (classical ?)
nonlinear-poisson-cg-residual   =  1.0d-6   !
2D
nonlinear-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-iterative is much faster than couple-all-false.

 

current-poisson-lin-eq-solv = BiCGSTAB  ! (default)
                            = QMRducpl  !
                            = LAPACK    !
1D
                            = PCG_KENT  !
2D
Switch 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:

CurrentPoissonLinEqSolvIter1D/2D/3D
CurrentPoissonLinEqSolvResid1D/2D/3D

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

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

current-poisson-iterations =  10       ! 1D  (15?)
current-poisson-iterations =  100      ! 2D  (25?)
current-poisson-iterations =  5        ! 3D

current-poisson-residual   =  1.0d-8   ! 1D
current-poisson-residual   =  1.0d-10  !
2D
current-poisson-residual   =  ???      !
3D

CurrentPoissonIterationsXD --> Maximum number of iterations.
CurrentPoissonResidualXD   -->
Precision for solution of current and Poisson equation.

 

 

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.

  EF,n(k+1) = w EF,n + (1 - w) EF,n(k)
  EF,p(k+1) = w EF,p + (1 - w) EF,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   ! 3D

Relaxation 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.,
         choosing current-block-relaxation-Fermi = 0d0  is not really a good idea.

 

 

Coupled current-Poisson equation:

coupled-current-poisson-iterations =  100      ! 1D
coupled-current-poisson-iterations =  100      !
2D
coupled-current-poisson-iterations =  100      !
3D

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

coupled-current-poisson-stepmax    =  1.0d-2   ! 1D
coupled-current-poisson-stepmax    =  1.0d-2   !
2D
coupled-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 subroutine newton).

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 from MODULE fermi_level and 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 if solve-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    = ! 1D
                                  = !
    2D/3D
    current-problem-residual      = 1d-10
    This loop is iterated CurrentProblemIterations1D/2D/3D times, its default value is CurrentProblemIterations1D/2D/3D (default: 5 (1D) / 1 (2D/3D)).
    current-problem-residual checks 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   !
    3D
    Iteration 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-14
    Default value for residual is CurrentProblemCGResidual1D/2D/3D = 1.0d-14.
    (Value that is actually used: CurrentProblemCGIterations1D/2D/3D_actual)
    Old name was: resid_cg_fermi

    For 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 than LAPACK-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 than LAPACK-ZHEEV)
                1D
                2D/3D not implemented
LAPACK-ZHEEVR = ZHEEVR from LAPACK package (calculates selected eigenvalues faster than LAPACK-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) for schroedinger-masses-anisotropic = box only
LAPACK-ZHBGVX = ZHBGVX from LAPACK package for banded matrix
      
          2D/3D for schroedinger-masses-anisotropic = box only
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 than ARPACK but 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 H2 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/3D ARPACK 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     ! 1D
schroedinger-1band-iterations =  1000     !
2D
schroedinger-1band-iterations =  1000     !
3D
maximum number of iterations

schroedinger-1band-residual   =  1.0d-12  !
1D
schroedinger-1band-residual   =  1.0d-12  !
2D
schroedinger-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 [Emin , Emax] the eigenvalues are sought.
For electrons: [Emin , Emax] = [ Ec,min , Ec,min + schroedinger-chearn-el-cutoff ]
For holes:   [Emin , Emax] = [ Ev,max - schroedinger-chearn-hl-cutoff , Ev,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 package
LAPACK-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 H2 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    =      !
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 yet
poisson-boundary-condition-along-y = periodic  !
                                   = Neumann   !
Neumann boundary conditions (default)
                                   = Dirichlet !
not implemented yet
poisson-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 factor ScalePoisson3D.
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 and vanHalen-Pulfrey are not implemented any more, so Trellakis is the default. Goano is 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 in MODULE 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  !
1D
fermi-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:
- F3/2: - Goano (if (ABS(x) > 353.6d0) .AND. (ABS(x) < 354.2d0) Trellakis is used)
         (The calculation for ABS(x) from 353.6 up to 354.2 fails with Goano on SUN.)
       - Trellakis
- F1/2: - Goano (if (ABS(x) > 353.6d0) .AND. (ABS(x) < 354.2d0) Trellakis is used)
         (The calculation for ABS(x) from 353.6 up to 354.2 fails with Goano on SUN.)
       - Trellakis
- F-1/2: - Goano (if (ABS(x) > 353.6d0) .AND. (ABS(x) < 354.2d0) Trellakis is used)
         (The calculation for ABS(x) from 353.6 up to 354.2 fails with Goano on SUN.)
       - Trellakis
- F-3/2: - Goano (if (ABS(x) > 353.6d0) .AND. (ABS(x) < 354.2d0) Trellakis is used)
         (The calculation for ABS(x) from 353.6 up to 354.2 fails with Goano on SUN.)
       - Trellakis

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

 

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.
     Ec' = Ec - e phi
     Ev' = Ev - 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 %a as 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.0d0  ! [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 edge Ec(x) = Ec,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 to 0 [V]
                = no   ! (default is 'no')

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 with flow-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.

 

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 Ec(x) has a minimum and where the valence band edge Ev(x) has a maximum.
Note that the position x might not be identical for MIN(Ec) and MAX(Ev).
So we can tell the eigenvalue solver to search for energies around the middle of the band gap, i.e. Eseparate = (Ec - Ev) / 2.
Then all eigenvalues that are larger than Eseparate are electron eigenvalues and those that are smaller in energy are hole eigenvalues.

The actual value of Eseparate 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 that E_separate is (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-electrons is used) are sought above separation-energy-shift-eV [eV]
-      hole eigenvalues (if $quantum-model-holes is 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 (dx2,dy2,dz2). 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(kx,ky,kz) in a semiconductor:

  • Luttinger parameters:                              gamma1, gamma2, gamma3 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: gamma1, gamma2, gamma3, Deltaso, (kappa)
  • DKK parameters:       L, M, N,                               Deltaso, (kappa)

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

  • Luttinger parameters: gamma1', gamma2', gamma3', Deltaso, (kappa'), Egap, Ep, S ( = 1 + 2 F), B
  • DKK parameters:       L', M'=M, N',                          Deltaso, (kappa'), Egap, Ep, 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 gamma1, gamma2, gamma3 (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 Egap.
    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 Egap is automatically taken into account.
    This is important, as L', N', kappa' (or gamma1', gamma2', gamma3', 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 Ep 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  ! 
    EP in units of [eV]           (default)
                                   P    ! 
    P   in units of [eV Angstrom]
    This choice affects the value of Ep 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 gamma1, gamma2 and gamma3 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 = gamma1   gamma2   gamma3   ---.
                      = 6x6kp-kappa      Use database entries Luttinger-parameters = gamma1   gamma2   gamma3  kappa.
                      = 6x6kp-kappa-only
Use database entries Luttinger-parameters = ---      ---      ---    kappa, and L, M, N (or L', M', N') parameters.
                      = 8x8kp           
Use database entries Luttinger-parameters = gamma1'  gamma2'  gamma3' ---. Note: Here, modified (i.e. 8-band k.p) Luttinger parameters have to be specified.
                      = 8x8kp-kappa     
Use database entries Luttinger-parameters = gamma1'  gamma2'  gamma3' kappa'. Note: Here, modified (i.e. 8-band k.p) Luttinger parameters have to be specified.
                      = 8x8kp-kappa-only
Use database entries Luttinger-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 (gamma1,  gamma2, gamma3, i.e. Luttinger-parameters, and optionally also kappa) instead of the DKK notation (L,M,N, i.e. 6x6kp-parameters and 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 Egap, Deltaso, EP and me.
                                  = 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 gamma1',  gamma2', gamma3' directly if Luttinger-parameters = 8x8kp or = 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 Egap, Deltaso, EP and me.
                                = LMN            ! Calculate L',M'=M,N'. Note: S is taken from database entry.
                                = S              ! S is calculated from Egap, Deltaso, EP and me.
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 take 8x8kp-parameters from database or input file (apart from B, EP).
Note: S is calculated from Egap, Deltaso, EP and me.

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, EP, S).
Note: S is taken from database and is not calculated.

S: S is calculated from Egap, Deltaso, EP and me.
Don't take 8x8kp-parameters from database or input file (apart from L',M'=M,N',B, EP).

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,Egap,Deltaso,EP,me).

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,EP,L',N',N+',N-' that are calculated anew.
ZERO: Rescale 8-band k.p parameters so that S=0. This affects S,EP,L',N',N+',N-' that are calculated anew.

This flag will overwrite the 8-band k.p parameters (L',N',N+',N-',EP,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 Deltax where kmax = 2pi /Deltax.
                     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  ! EP in units of [eV]           (default)
                                P    ! 
P   in units of [eV Angstrom]

Refers to EP parameter specified in 8x8kp-parameters in input file or database.
All parameters in the database are EP. 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 Hcv part of the k.p Hamiltonian.
vv
is related to the Hvv 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 the z direction.
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 a default direction 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), see SUBROUTINE 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?)

 

 

Calculate band edge energies

There are two algorithms implemented to calculate the valence band edge energies that are used as input for the 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 = (kx,ky,kz) = 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.

calculate-bandedge-energies = vb-method-1  ! (default)  (method implemented by M. Povolotskyi for zincblende, generalized by S. Birner for wurtzite) - This is the default since 2013-09-17.
                            = vb-method-2  !        (method implemented by M. Povolotskyi for zincblende and the method of M. Sabathil generalized by S. Birner for wurtzite) - implemented since 2013-09-17.
                            = vb-method-3  !        (method implemented by M. Povolotskyi for zincblende and the method of M. Sabathil for wurtzite (which is actually a zincblende version) - This was the default until 2013-09-17.
                            = vb-method-4  !        (method implemented by M. Sabathil for zincblende and the method of M. Sabathil for wurtzite (which is actually a zincblende version) - This was the default a long time ago.

 

 

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-electrons or $quantum-model-holes has to be specified.
           = full-band-density-electrons      ! 
FB-EFA: The k.p Hamiltonian is solved only once. $quantum-model-electrons has 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-electrons and $quantum-model-holes should be specified.
           = full-band-density-holes-only     ! 
FB-EFA: The k.p Hamiltonian is solved only once. Both $quantum-model-electrons and $quantum-model-holes should 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 !
similar to 2nd-order but uses Tse-Pal model and Tse-Pal material parameters
                    = 4th-order-Tse-Pal !
Tse-Pal model and Tse-Pal material parameters (4th order for zinc blende, 2nd order for wurtzite)

See e.g. - First- and second-order piezoelectricity in III-V semiconductors
         A. Beya-Wakata, P.-Y. Prodhomme, G. Bester
         Phys. Rev. B 84, 195207 (2011)
       - G. Bester et al., PRL 96, 187602 (2006) (zinc blende)
       - L. Pedesseau et al., APL 100, 031903 (2012) (wurtzite)

Here, one can switch on the second order effect of piezoelectricity. Obviously several piezoelectric parameters are needed.
  - zinc blende: e14              (1st order effect); B114  B124  B156 (2nd order effect)
  - wurtzite:      e33  e31  e15 (1st order effect); B311  B312  B313  B333  B115  B125  B135  B344 (2nd 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 2nd order model of A. Beya-Wakata et al., Phys. Rev. B 84, 195207 (2011) is used.
      4 parameters are needed.
      e14  B114  B124  B156   (default)
      The equations for the second order terms read:
      Px = 2 * B114 * eps_xx * eps_yz + 2 * B124 * eps_yz * (eps_yy + eps_zz) + 4 * B156 * eps_xz * eps_xy
      Px = 2 * B114 * eps_yy * eps_xz + 2 * B124 * eps_xz * (eps_zz + eps_xx) + 4 * B156 * eps_yz * eps_xy
      Pz = 2 * B114 * eps_zz * eps_xy + 2 * B124 * eps_xy * (eps_xx + eps_yy) + 4 * B156 * eps_yz * eps_xz
    - For wurtzite materials, the 2nd order model of L. Pedesseau et al., Appl. Phys. Lett. 100, 031903 (2012) is used.
      11 parameters are needed.
      e33  e31  e15  B311  B312  B313  B333  B115  B125  B135  B344   (default)
      The equations for the second order terms read
      Px = ?
      Py = ?
      Pz = ( B311 + B312 ) * eps_xx * eps_yy + 2 * B313 * SQRT( eps_xx * eps_yy ) * eps_zz + 0.5 * B333 * eps_zz2
  • 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 3rd and 4th order terms.
      4 parameters are needed.
      e14  B114  B124  B156   (default)
      Note that the Tse-Pal model assumes B114 = e114, B124 = e124 and B156 = e156 = 0.
    - For wurtzite materials, the 2nd 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).
      e33 e31 e15  e311 e313 e333
      The wurtzite material parameters for this model can only be specified in the input file but not in the database file database_nn3.in.
     
  • 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, 2nd, 3rd and 4th order terms for zinc blende.
      13 parameters are needed.
      e14  B114  B124  B156   e1114 e1124 e1224 e1234  e11114 e11124 e12224 e12234 e11234
      Note that the Tse-Pal model assumes B114 = e114, B124 = e124 and B156 = e156 = 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 mo 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 * 1018 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 KSL values are calculated.

 superlattice-save-wavefunctions = yes                     !
(default)
                                 = no                      !
 no
: Program will not save wave functions for KSL 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)

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 exction 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    ! 1D/3D electron ground state
                                        = 2    ! 1D/3D
2nd electron state
                                        = ...
 exciton-hole-state-number              = 1    ! 1D/3D
hole ground state
                                        = 2    ! 1D/3D
2nd 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 2nd electron state with the hole ground state.
 - Calculate exciton energy of the electron ground state with the 2nd hole state.
 - Calculate exciton energy of the 2nd electron state with the 2nd hole 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 1st k.p hole state corresponds to the 1st single-band holes state.
        The 3rd k.p hole state corresponds to the 2nd single-band holes state.
        The 5rd k.p hole state corresponds to the 3nd single-band holes state.



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

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.


2D:
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 1st k.p hole state corresponds to the 1st single-band holes state.
        The 3rd k.p hole state corresponds to the 2nd single-band holes state.
        The 5rd k.p hole state corresponds to the 3nd single-band holes state.

 

In 1D we currently have a routine that calculates the exciton correction in type-I quantum wells (single-band only, suited for electron and hole ground states 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.0d0   ! [nm] 1D only (optional)
 exciton-lambda-step-length = 0.05d0  ! [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) |2
                    
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) |2
                    
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
                     = 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.

tight-binding-parameters = -3.53284d0 ! Esa         (GaAs)
                            0.27772d0 ! Epa
                           -8.11499d0 ! Esc
                            4.57341d0 ! Epc
                           12.33930d0 ! Es_a
                            4.31241d0 ! Es_c
                           -6.87653d0 ! Vss
                            1.33572d0 ! Vxx
                            5.07596d0 ! Vxy
                            0d0       ! Vs_s_
                            2.85929d0 ! Vsa_pc
                           11.09774d0 ! Vsc_pa
                            6.31619d0 ! Vs_a_pc
                            5.02335d0 ! Vs_c_pa
                            0.32703d0 ! Delta_so_a
                            0.12000d0 ! Delta_so_c

Note: Special case for graphene (3 parameters only, nearest-neighbor tight-binding approximation):

 tight-binding-parameters =  0d0     ! [eV] E2p: site energy of the 2pz atomic orbital (orbital energy)
                            -3.013d0 ! [eV] gamma0: C-C transfer energy (nearest-neighbor, nn)
                             0.129d0 ! []   s0: 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.28d0  ! [eV] E2p: site energy of the 2pz atomic orbital (orbital energy)
                            -2.97d0  ! [eV] gamma0: C-C transfer energy (nearest-neighbor, nn)
                             0.073d0 ! []   s0: denotes the overlap of the electronic wave function on adjacent sites (nn)
                            -0.073d0 ! [eV] gamma1: (2nd-nn)
                             0.018d0 ! []   s1:          (2nd-nn)
                            -0.33d0  ! [eV] gamma2: (3rd-nn)
                             0.026d0 ! []   s2:          (3rd-nn)

Special option for third-nearest-neighbor tight-binding approximation in graphene:

tight-binding-calculate-parameter = no     ! (default) use 7 parameters (E2p,gamma0,gamma1,gamma2,s0,s1,s2)
                                  = E2p    !
            use 6 parameters (       gamma0,gamma1,gamma2,s0,s1,s2)
                                  = gamma1 !
            use 6 parameters (E2p,gamma0,             gamma2,s0,s1,s2)
E2p or gamma1 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 = 0d0                   0.0001d0               4d0         ! [eV] [eV] [] (default)
                               = 0d0                   0.0001d0                           ! [eV] [eV]
                               = 0d0                                                      ! [eV]
(default, i.e. "ideal" graphene layer)
                               = 0.100d0                                                  ! [eV]
                               = 0.100d0               0.001d0                3d0         ! [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 1st 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 2nd value of these array entries (optional) corresponds to the energy grid resolution of the Gaussian integration in units of [eV].
The 3rd 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 2nd and 3rd value has implications on CPU time. 3d0 is a very reasonable value for sigma_factor and also an energy grid resolution of 1 meV (0.001d0) 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 - El,m ) / 2 t - 1 )

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