|
| |
Electronic Structure
This chapter comprises:
- Storage of data and setup of arrays on electronic structure
- Information on how to extract information on electronic structure
- Rotation of Hamilton
3.1 Storage of rotation
3.2 Setup of rotation
- Information on how to assemble the bulk Hamiltonian
- Setup of k.p Hamiltonian
- Calculate new k.p states
-
Setup of single-band Schrödinger equation including magnetic field
- Separation models:
Eigenvalues
<--> Thomas-Fermi
Separation between quantum and continuum states
- Calculate eigenvalues
of k.p Schrödinger equation
- Calculate
eigenvalues of single-band Schrödinger equation
- Environment where these routines are called
In the main file the arrays for the electronic structure are allocated and
defined in the sequence:
CALL
setup_quantum_models
CALL setup_kp_parameters
CALL allocate_quantum_states
This subroutine allocates the pointer def_quantum_modelsM in MODULE
quantum_models.
Here the full information about degeneracies and effective masses is
stored. See the comments in the source code for more information.
The structures sgChargeV(1)%qcV, sgChargeV(2)%qcV in
MODULE quantum_solutions are allocated and set to default. The variables
sgChargeV(1)%qcV, sgChargeV(2)%qcV contain
the solution of the single-band Schrödinger equation.
More details here.
(1) = el
(2) = hl
The corresponding Fermi levels for the single-band states are defined
in
MODULE
fermi_level:
fermi_sg_elV
fermi_sg_hlV
They are allocated in SUBROUTINE
setup_quantum_models.
Note: In order to simulate population inversion in lasers, we use different
Fermi levels for the quantum states.
So each energy level could have a different Fermi level, in principle.
However, this is not implemented so far. So all electron quantum states have
the same electron Fermi level and all hole quantum states have the same hole
Fermi level.
To be more flexible each wavefunction has its own position-dependent Fermi
level, stored in fermi_sg_elV, fermi_sg_hlV. These are taken for the density calculation.
Up to now,
we assume one Fermi level for electrons and one for holes.
--> SUBROUTINE
SetNewFermiLevel in
MODULE ModifyFermiLevel is used after solving the drift-diffusion equation
in order to update the Fermi levels
accordingly.
The k.p parameters in MODULE
kp_parameters
are defined. These parameters are provided by routines in MODULE
input_kp_parameters.
Look up documentation there. See script for detailed information on
k.p parameters. The storage format is described in more detail in
"Input of k.p parameters".
Units: [eV], [Angstrom]
Allocates k.p quantum states and defines them. The
k.p solutions are stored in MODULE
quantum_solutions: kpChargeV.
More details here.
Allocates their Fermi levels and sets them to default: MODULE
fermi_level:
fermi_kp_elV,
fermi_kp_hlV.
For the Fermi levels for each spinor of the k.p equation it holds.
-> Each quantum cluster has its own Fermi level.
-> Each eigenvalue has its own Fermi level.
-> Each boundary condition has its own Fermi level.
So far, all of these different Fermi levels are identical, i.e. we
assume one Fermi level for electrons and one for holes.
The necessary information is
provided by the subroutine get_kp_data in MODULE input_quantummodels.
The
information about number of k|| and the density integration is obtained
via subroutine
get_dim_BZ_2D
with option calculate_kL = .FALSE.
In MODULE extract_kp_parameters there are various (auxiliary) functions that provide information about k.p:
MODULE extract_kp_parameters
SUBROUTINE get_kp_guenther_parametersSUBROUTINE check_kpSUBROUTINE Get_kp_MaterialParametersSUBROUTINE get_which_kp
!--------------------------------------------------------------------
SUBROUTINE check_kp(i,i_qr_guent,num_qr,kind_kpC,kpL)
!--------------------------------------------------------------------
! Input:
! i
!
!
! kpL = .TRUE. i
is contained within a quantum
region.
! num_qr
! kind_kpC = '8x8_zb','8x8_wz','6x6_zb','6x6_wz'
! i_qr_guent =
!
!
! map_qr_guentherV(num_qr)%number_of_points [MODULE gitter1D]
! num_qr.
! map_qr_guentherV(num_qr)%guenther_grid_points
!
! --> map_qr_guentherV(num_qr)%guenther_grid_points(i_qr_guent)=i
!
!----------------------------------------------------------------------
This subroutine returns for any point on
material grid the quantum region, if kp is done, if so what kind of kp, and in
i_qr_guent the position within the point numbering of that specific
quantum region [see above].
!--------------------------------------------------------------------
SUBROUTINE get_kp_guenther_parameters(i,parV,kpL)
!---------------------------------------------------------------------
! Returns for point i
on material grid k.p parameters
which are used for matrix assembly.
!
!
! L1 = L2 = L
! M1 = M2 = M3 = M
! N1 = N2 = N
! B1 = B2 = B3 = B
! P1 = P2 = P
! s1 = s2 = s
! d1 = 0
! d2 = d3 = d/3
!
! If this point does not lie in a k.p
region --> kpL = .FALSE.
!
! Input:
! i ...
!
!
!
! kpL .TRUE.
! parV
!
! parV(1) = L1
! parV(2) = L2
! parV(3) = M1
! parV(4) = M2
! parV(5) = M3
! parV(6) = N1
! parV(7) = N2
! parV(8) = B1
! parV(9) = B2
! parV(10) = B3
! parV(11) = P1
! parV(12) = P2
! parV(13) = s1
! parV(14) = s2
! parV(15) = E_c
! parV(16) = E_v
! parV(17) = d1
! parV(18) = d2
! parV(19) = d3
!
! [eV]
and [Angstrom]
!
!---------------------------------------------------------------------
Please note: This subroutine is valid for
wurtzite and for zincblende (see script for meaning of parameters). Same
routine for 2D and 3D.
par_minV -->
material parameters between i-1
and
i
par_andV -->
material parameters between i
and
i+1
parV -->
material
parameters on i
!-------------------------------------------------------------------------
SUBROUTINE Get_kp_MaterialParameters(i_phys_min,i_phys,i_phys_and , &
j_phys_min,j_phys,j_phys_and , &
k_phys_min,k_phys,k_phys_and , &
par_i_minV,par_j_minV,par_k_minV, &
par_i_andV,par_j_andV,par_k_andV,
&
parV,par_octM)
!-------------------------------------------------------------------------
!
Input:
!
! 1D: i_phys_min,i_phys,i_phys_and
! 2D: i_phys_min,i_phys,i_phys_and,
! j_phys_min,j_phys,j_phys_and
! 3D: i_phys_min,i_phys,i_phys_and,
! j_phys_min,j_phys,j_phys_and,
! k_phys_min,k_phys,k_phys_and
!
(i_phys_min = i_phys = i_phys_and) just parV
is valid.
!
!
!
! par_i_minV -> i_phys_min,i_phys
! par_i_andV -> between i_phys ,i_phys_and
! par_j_minV -> between j_phys_min,j_phys
! par_j_andV -> between j_phys ,j_phys_and
! par_k_minV -> between k_phys_min,k_phys
! par_k_andV -> between k_phys ,k_phys_and
! parV ->
on (i_phys,j_phys)
!
! par_andV(1) = L1
! par_andV(2) = L2
! par_andV(3) = M1
! par_andV(4) = M2
! par_andV(5) = M3
! par_andV(6) = N1
! par_andV(7) = N2
! par_andV(8) = B1
! par_andV(9) = B2
! par_andV(10) = B3
! par_andV(11) = P1
! par_andV(12) = P2
! par_andV(13) = s1
! par_andV(14) = s2
!
! only for parV
!
! parV(15) = E_c (including electrostatic potential)
! parV(16) = E_v (including electrostatic potential)
! parV(17) = d1
! parV(18) = d2
! parV(19) = d3
!
! [eV]
and [Angstrom]
!
! Note: The electrostatic potential
phiV is contained in MODULE
potentials.
!
!----------------------------------------------------------------------
!---------------------------------------------------------------------
SUBROUTINE get_which_kp(num_qr,kpL,kind_kpC,chargeC)
!---------------------------------------------------------------------
! Input:
! num_qr ...
!
!
! kpL,kind_kpC,chargeC
!
!---------------------------------------------------
!
! num_qr.
! --> kpL = .TRUE. (
.FALSE. otherwise)
! kind_kpC = '8x8','6x6'
! ('---' if kpL=.FALSE.)
! chargeC = 'el','hl','bo'
! 'bo' --> k.p done for both carriers
! '--' --> if no k.p done (kpL=.false.)
!
!---------------------------------------------------------------------
3.1 Storage of rotation / Rotation of k.p
Hamiltonian
The rotation is carried out from the crystal to
the calculation system as described in 'Rotations'. The rotated
k.p-Hamiltonian has the form as
specified in MODULE
rotate_kp.
/ r11 r12 r13 \
rotation matrix : R = | r21 r22 r23 |
\ r31 r32 r33 /
For definiton of k.p wurtzite matrix, follow
this link.
REAL(8) :: r11,r12,r13, &
r21,r22,r23,
&
r31,r32,r33
Valence band
1. Rotation of Hvv
(3x3)-matrix with L,M,N in basis |x>, |y>, |z>.
The rotation of Hvv is carried out as follows:
H'(k') = RH(RTk')RT
k': k vector in new coordinate system which is the simulation
system.
H': Hamiltonian in new coordinate system which is the simulation system.
H: Hamiltonian in crystal system.
RTk' = k: k vector in crystal system.
Since the basis states |x>, |y>, |z> are functions in the crystal sytem, they
must be transferred into the simulation system.
H'vv is a
(3x3) matrix which has 9
entries which are functions of 6 variables k'x², k'y², k'z²,
k'xk'y=k'yk'x, k'xk'z=k'zk'x,
k'yk'z=k'zk'y.
9 variables k'x², k'y², k'z²,
k'xk'y, k'xk'z,
k'yk'z,k'yk'x, k'zk'x, k'zk'y
We are writing the rotated (3x3) H'vv matrix as the sum of
several (3x3) matrices ordered by k'x², k'y², k'z², k'xk'y,
k'xk'z, k'yk'z:
ordered by k'x², k'y², k'z², k'xk'y,
k'xk'z, k'yk'z, k'yk'x,
k'zk'x, k'zk'y:
H'vv = ( 3x3 ) k'x² + ( 3x3 ) k'y² + (
3x3 ) k'z² + ( 3x3 ) k'xk'y + ( 3x3 ) k'xk'z
+ ( 3x3 ) k'yk'z
( 3x3 ) k'x² + ( 3x3 ) k'y² + (
3x3 ) k'z² + ( 3x3 ) k'xk'y + ( 3x3 ) k'xk'z
+ ( 3x3 ) k'yk'z + ( 3x3 ) k'yk'x + ( 3x3 ) k'zk'x
+ ( 3x3 ) k'zk'y
or
H'vv = RM1.k'x² + RM2.k'y²
+ RM3.k'z² + RM4.k'xk'y
+ RM5.k'xk'z + RM6.k'yk'z
RM1.k'x² + RM2.k'y²
+ RM3.k'z² + RM4.k'xk'y
+ RM5.k'xk'z + RM6.k'yk'z
+ RM7.k'yk'x + RM8.k'zk'x + RM9.k'zk'y
RMi
(i=1...6)
(i=1...9)
i=1...6
i=1...9
i=1 <--> k'x²
i=2 <--> k'y²
i=3 <--> k'z²
i=4 <--> k'xk'y
i=5 <--> k'xk'z
i=6 <--> k'yk'z
i=7 <--> k'yk'x
i=8 <--> k'zk'x
i=9 <--> k'zk'y
RMi is a
(3x3) matrix which has 9
entries which are functions of 7 variables L1, L2, M1,
M2, M3, N1, N2.
Again, we are writing the rotated (3x3) RMi matrix as the sum of
several (3x3) matrices ordered by L1, L2, M1, M2,
M3, N1, N2:
RMi = ( 3x3 ) L1 + ( 3x3 ) L2 + (
3x3 ) M1 + ( 3x3 ) M2 + ( 3x3 ) M3 + ( 3x3 ) N1
+ ( 3x3 ) N2
or
RMi = RMiL1.L1 + RMiL2.L2
+ RMiM1.M1 + RMiM2.M2
+ RMiM3.M3 + RMiN1.N1
+ RMiN2.N2
(L,M,N Dresselhaus parameters)
RMij
(i=1...6, j=1...7)
(i=1...9, j=1...7)
j=1...7
j=1 <--> L1
j=2 <--> L2
j=3 <--> M1
j=4 <--> M2
j=5 <--> M3
j=6 <--> N1
j=7 <--> N2
rotationCoeff_wurM(1:6,1:7)%matrixM(1:3,1:3) = rotationCoeffM(i,j)%matrixM(3x3) = RMij
rotationCoeff_wurM(1:9,1:7)%matrixM(1:3,1:3) = rotationCoeffM(i,j)%matrixM(3x3) = RMij
REAL(8),DIMENSION(6,7,3,3) :: rotationCoeff_wurM
Representation for wurtzite is also used for zincblende:
L=L1=L2
M=M1=M2=M3
N=N1=N2
The diagonal elements Ev+hbar²/(2m0) k² do
not change under rotation.
Conduction band:
2. Rotation of Hcc
For definiton of Hcc wurtzite matrix, follow
this link.
Hcc = (Hcc)ij
H'11 = H11 + Ec + hbar²/(2m0).s2.(k'x²+k'y²+k'z²)
+
+ hbar²/(2m0).(s1-s2).(r13²k'x²
+ r23²k'y² + r33²k'z² + 2 r13r23²k'xk'y
+ 2 r13r33²k'xk'z + 2 r23r33k'yk'z
)
3. Rotation of Hcv
Follow this link
for the structure of the matrix Hcv.
H'cv is a
(2x6) matrix which can
be split into four (1x3) matrices. Two of these are zero, the other two are
identical and have these three components Hcvsx, Hcvsy,
Hcvsz which are functions of 6 variables k'x²,
k'y², k'z², k'xk'y=k'yk'x, k'xk'z=k'zk'x,
k'yk'z=k'zk'y.
6 variables k'xk'y=k'yk'x, k'xk'z=k'zk'x,
k'yk'z=k'zk'y,
k'x,
k'y,
k'z. (???)
Hcvsx = B ky kz + i P kx
Hcvsy = B kx kz + i P ky
Hcvsz = B kx ky + i P kz
The B terms are functions in second order of k, the P terms are functions in
first order of k, thus we divide the matrix components into 2 parts, one for B
and one for P.
H'cvsx = H'cvsx
(only B part) + H'cvsx (only P part)
B Part
x' = R . x
This is how the rotated H'cvsn should look
like:
H'cvsx (only B part) = < s | H | x' > = r11
B1 ky kz + r12
B2 kx
kz + r13 B3 kx ky
H'cvsy (only B part) = < s | H | y' > = r21
B1 ky kz + r22
B2 kx
kz + r23 B3 kx ky
H'cvsz (only B part) = < s | H | z' > = r31
B1 ky kz + r32
B2 kx
kz + r33 B3 kx ky
However, we are writing the rotated Hcvsn (n=x,y,z) component (=scalar) as the sum of
several scalar components ordered by k'x², k'y², k'z², k'xk'y,
k'xk'z, k'yk'z:
H'cvsx (only B part) = < s | H | x' > = RMx1.k'x² + RMx2.k'y²
+ RMx3.k'z² + RMx4.k'xk'y
+ RMx5.k'xk'z + RMx6.k'yk'z
H'cvsy (only B part) = < s | H | y' > = RMy1.k'x² + RMy2.k'y²
+ RMy3.k'z² + RMy4.k'xk'y
+ RMy5.k'xk'z + RMy6.k'yk'z
H'cvsz (only B part) = < s | H | z' > = RMz1.k'x² + RMz2.k'y²
+ RMz3.k'z² + RMz4.k'xk'y
+ RMz5.k'xk'z + RMz6.k'yk'z
H'cvsn = RMn1.k'x² + RMn2.k'y²
+ RMn3.k'z² + RMn4.k'xk'y
+ RMn5.k'xk'z + RMn6.k'yk'z
RMni is a scalar! (i=1...6)
(n=1...3)
RMx1 = RMx1B1
B1 + RMx1B2 B2 + RMx1B3 B3;
RMx2 = ...; RMx3 = ... ; RMx4 =
...; RMx5 = ... ; RMx6 = ...;
RMy1 = RMy1B1
B1 + RMy1B2 B2 + RMy1B3 B3;
RMy2 = ...; RMy3 = ... ; RMy4 =
...; RMy5 = ... ; RMy6 = ...;
RMz1 = RMz1B1
B1 + RMz1B2 B2 + RMz1B3 B3;
RMz2 = ...; RMz3 = ... ; RMz4 =
...; RMz5 = ... ; RMz6 = ...;
RMni = RMniB1 B1 + RMniB2 B2
+ RMniB3 B2 (Bm=1...3) (wurtzite)
RMniBm is a scalar! (Bm=1...3)
(i=1...6) (n=1...3)
rotationCondM --> conduction band (for B)
A vector of dimension 3 stores this information:
[ H'cvsy (only B part), H'cvsz (only B part), H'cvsz (only B part)
]
--> vector(n=1..3) = (Hsx Hsy Hsz)
REAL(8),DIMENSION(6,3,3) :: rotationCond_wurM
Representation for wurtzite is also used for zincblende:
B = B1 = B2 =B3
P Part
x' = R . x
This is how the rotated H'cvsn should look like:
H'cvsx (only P part) = < s | H | x' > = r11
P2 kx + r12
P2 ky + r13
P1 kz
H'cvsy (only P part) = < s | H | y' > = r21
P2 kx + r22
P2 ky + r23
P1 kz
H'cvsz (only P part) = < s | H | z' > = r31
P2 kx + r32
P2 ky + r33
P1 kz
However, we are writing the rotated Hcvsn (n=x,y,z) component (=scalar) as the sum of
several scalar components ordered by kx,
ky, kz:
H'cvsx (only P part) = < s | H | x' > = RMx1.k'x + RMx2.k'y
+ RMx3.k'z
H'cvsy (only P part) = < s | H | y' > = RMy1.k'x + RMy2.k'y
+ RMy3.k'z
H'cvsz (only P part) = < s | H | z' > = RMz1.k'x + RMz2.k'y
+ RMz3.k'z
H'cvsn
= RMn1.k'x + RMn2.k'y
+ RMn3.k'z
RMni is a scalar! (i=1...3)
(n=1...3)
RMx1 = RMx1P1 P1 + RMx1P2 P2;
RMx2 = ...; RMx3 = ... ; RMx4 =
...; RMx5 = ... ; RMx6 = ...;
RMy1 = RMy1P1 P1 + RMy1P2 P2;
RMy2 = ...; RMy3 = ... ; RMy4 =
...; RMy5 = ... ; RMy6 = ...;
RMz1 = RMz1P1 P1 + RMz1P2 P2;
RMz2 = ...; RMz3 = ... ; RMz4 =
...; RMz5 = ... ; RMz6 = ...;
RMni = RMniP1 P1 + RMniP2 P2 (Pm=1...2) (wurtzite)
RMniPm is a scalar! (Pm=1...2)
(i=1...3) (n=1...3)
rotation_P_M --> conduction band (for P)
A vector of dimension 3 stores this information:
[ H'cvsx (only P part), H'cvsy (only
P part), H'cvsz (only P part) ]
--> vector(n=1..3) = (Hsx Hsy Hsz)
COMPLEX(8),DIMENSION(3,2,3) :: rotation_P_wurM
complex!!!
Representation for wurtzite is also used for zincblende:
P = P1 = P2
Crystal field splitting
This
link for details.
Here, we don't need to split the rotated matrix D' into the k'
components. We just store the rotated matrix after we have done the basis
transformation:
D' = R.D.RT
Crystal field splitting: crysSplit_wurM (d1)
wurtzite: d1
zincblende: d1=0
REAL(8),DIMENSION(3,3) :: crysSplit_wurM
Spin-orbit interaction
This
link for details.
Here, we don't need to split the rotated matrix HSO' into
the k' components. We just store the rotated matrix after we have done the basis
transformation:
HSO' = R.HSO.RT
rotationSpin_wurM rotation of spin matrix
for wurtzitze: d2, d3
for zincblende: d2=d/3 , d3=d/3
COMPLEX(8),DIMENSION(6,6) :: rotationSpin_wurM
Strain
This link for
details.
In subroutine SetupStrainHamiltonianVB (MODULE
mod_setup_rot_kp)
the strain matrix Hstrain will be set up in the crystal
coordinate system. Then it will be rotated into the simulation coordinate system
leading to H'strain. Here we don't need to split H'strain
into the k' components.
H'strain = R.Hstrain.RT
The strain Hamiltonian Hstrain
is completely set up in
the crystal system and then rotated directly to the calculation system. This is done in
subroutine set_strain_mat_physgrid_point / subroutine
set_strain_mat_guent_point in MODULE
set_strain_mat.
Strain for holes: zincblende and wurtzite in same matrix.
Holes: StrainHamiltonianVB_3x3M(3,3)
COMPLEX(8),DIMENSION(3,3) :: StrainHamiltonianVB_3x3M
NOTE: The data in this module only depend on the rotation matrix.
--> Only
calculated anew when rotation has changed.
--> Stores
for which direction last
rotation has been calculated in.
last_rot_matM(3,3)
If
rotation not defined: last_rot_matM = 0.0d0 ( = starting value)
REAL(8),DIMENSION(3,3) :: last_rot_matM=0.0d0
3.2 Setup of rotation
The rotation matrix is defined in
update_rotation /
update_rotation_guent in MODULE
mod_update_rotation. These subroutines update the rotation
matrix in MODULE
rotate_kp
and the arrays rotationCoeff_wurM,
rotationCond_wurM, rotation_P_wurM
in MODULE
rotate_kp for a given point on
physical/material grid.
The arrays in MODULE
rotate_kp are
defined with the help of the routines in MODULE
mod_setup_rot_kp:
!--------------------------------------------------!
!
! 1. subroutine setup_rotationCoeffM
!
! 2. subroutine setup_rotationSpin_wurM
!
! 3. subroutine setup_rotationCondM
!
! 4. subroutine setup_rotation_P_M
!
! 5. subroutine setup_crysSplit_wurM
!
! 6. subroutine SetupStrainHamiltonianVB
!
!--------------------------------------------------!
NOTE: Before these subroutines are
called, the rotation matrix in MODULE
rotate_kp
(i.e. r11,r12,r13,
r21,r22,r23,
r31,r32,r33) must have been defined.
In subroutine
setup_kpbulk
[MODULE kp_bulk] the rotated bulk k.p Hamiltonian is assembled. One can study in this subroutine, how the assembly
works:
!-----------------------------------------------------
subroutine
setup_kpbulk(i_guent,k1,k2,k3,hamiltonM)
!-----------------------------------------------------
For a given point i_guent on
material grid we are calculating th k.p matrix matM.
After that we are adding the strain matrix.
Note: k1,k2,k3 in
Angstrom (AA).
It sets up the strain matrix using directly strain from
material grid.
!---------------------------------------------------------------------------
USE derived_constants ,ONLY:h2b2m_evAA2
USE
rotate_kp ,ONLY:rotationCond_wurM,rotationCoeff_wurM,rotation_P_wurM, &
rotationSpin_wurM,crysSplit_wurM,StrainHamiltonianVB_3x3M,
&
r11,r12,r13,r22,r23,r33,r21,r31,r32
USE
mod_update_rotation ,ONLY:update_rotation_guent
USE mod_setup_rot_kp ,ONLY:setup_rotationCoeffM,setup_rotationSpin_wurM,
&
setup_rotationCondM,setup_rotation_P_M,setup_crysSplit_wurM
USE gitter
,ONLY:addressM,PointInfoM
USE set_strain_mat ,ONLY:set_strain_mat_guent_point
USE structure_functions ,ONLY:convert_guenther_to_phys2D
REAL(8),intent(in) :: k1,k2,k3
LOGICAL :: kpL
!-----------------------------------------------------------
REAL(8) :: ind1,ind
INTEGER :: i,j,adr1,adr2
INTEGER ,DIMENSION(4,3) :: coordM
! 2D
Here we are retrieving the k.p parameters for the points i_guent
from MODULE extract_kp_parameters.
USE
extract_kp_parameters,ONLY: get_kp_guenther_parameters
integer,intent(in) :: i_guent
real(8),dimension(19) :: parV
real(8) :: E_c,E_v,s1,s2,d1,d2,d3,l1,l2,
m1,m2,m3,n1,n2,b1,b2,b3,p1,p2,pot
CALL
get_kp_guenther_parameters(i_guent,parV,kpL)
Note: The definition of parV is the same for
wurtzite and zincblende.
l1 = parV(1) ! Dresselhaus
parameters
l2 = parV(2)
m1 = parV(3)
m2 = parV(4)
m3 = parV(5)
n1 = parV(6)
n2 = parV(7)
b1 = parV(8) !
b2 = parV(9)
b3 = parV(10)
p1 = parV(11) ! s
p2 = parV(12)
s1 = parV(13) ! Inverse effective masses of electrons
s2 = parV(14)
E_c = parV(15) !
E_v = parV(16) !
d1 = parV(17) !
d2 = parV(18) !
d3 = parV(19) !
The electrostatic potential is taken from
phiV.
USE define_geometry,ONLY : n_material_points
USE gitter ,ONLY :
addressM
IF (domain_dimension==1) THEN
if (i_guent==1) then
adr1=addressM(1,1,1)
adr2=adr1
else if (i_guent==n_material_points) then
adr1=addressM(n_material_points+1,1,1)
adr2=adr1
else
if
(PointInfoM(i_guent,1,1)==0) then
adr1=addressM(i_guent,1,1)
else
adr1=addressM(i_guent,1,1)+1
end if
adr2=addressM(i_guent+1,1,1)
end if
pot = 0.5d0 * (phiV(adr1) +
phiV(adr2))
ELSE IF (domain_dimension==2) THEN
SetupStrainHamiltonianVB CALL convert_guenther_to_phys2D(i_guent,coordM)
pot = 0.25d0 * (phiV(coordM(1,3)) + phiV(coordM(2,3)) + &
phiV(coordM(3,3)) + phiV(coordM(4,3)) )
END IF
E_c=E_c-pot
E_v=E_v-pot
Now E_c and E_v
include the potential.
Conduction band energy (including potential and
including strain)
Here we define the rotation matrix and the rotated
k.p Hamiltonian.
CALL
update_rotation_guent(i_guent)
(Here the three routines
call
setup_rotationCoeffM
call
setup_rotationCondM
call
setup_rotation_P_M
are called.)
CALL
setup_rotationSpin_wurM(d2,d3)
CALL
setup_crysSplit_wurM(d1)
CALL
set_strain_mat_guent_point(i_guent)
Here we are retrieving the
deformation potentials.
call
get_deformation_pot_vb(i_guent,l1,l2,m1,m2,m3,n1,n2,a_abs)
a_abs
is ac but it isn't used as it is already included in E_c
which is defined in the k.p parameters
described above.
Here we are extracting the strain tensor for grid
point i_guent (specified in
crystal coordinate system).
call
get_strain_voigt_notation_cxyz(i_guent,strainM)
call
SetupStrainHamiltonianVB(epsM,l1,l2,n1,n2,m1,m2,m3)
Sets up
strain matrix H'strain
for H'vv
part [ (6*6) = two 3*3 blocks ] and
carries out rotation numerically.
For
rotationCoeff_wurM(1:6,1:7,1:3,1:3) = rotationCoeffM(i,,j,3x3) = RMij
| RMi
(i=1...6) |
RMij
(i=1...6, j=1...7) |
i=1...6
i=1 <--> k'x²
i=2 <--> k'y²
i=3 <--> k'z²
i=4 <--> k'xk'y
i=5 <--> k'xk'z
i=6 <--> k'yk'z |
j=1...7
j=1 <--> L1
j=2 <--> L2
j=3 <--> M1
j=4 <--> M2
j=5 <--> M3
j=6 <--> N1
j=7 <--> N2 |
(for details, see above!)
we are looping over i and j.
do i=1,6
select case(i)
case(1)
ind=k1**2
case(2)
ind=k2**2
case(3)
ind=k3**2
case(4)
ind=k1*k2
case(5)
ind=k1*k3
case(6)
ind=k2*k3
end select
do j=1,7
select case(j)
case(1)
ind1=l1
case(2)
ind1=l2
case(3)
ind1=m1
case(4)
ind1=m2
case(5)
ind1=m3
case(6)
ind1=n1
case(7)
ind1=n2
end select
The Hamiltonian matrix is always 8x8, even in the case for 6x6 kp. In the
latter case, the first two rows and first two columns are not used. Here we
calculate H'vv and H'vv which are
identical.
complex(8),dimension(8,8),intent(out) :: hamiltonM
hamiltonM = cmplx(0d0)
hamiltonM(3:5,3:5) =
hamiltonM(3:5,3:5) +
dcmplx(ind*ind1*rotationCoeff_wurM(i,j,:,:))
! H'vv
hamiltonM(6:8,6:8) = hamiltonM(6:8,6:8) +
dcmplx(ind*ind1*rotationCoeff_wurM(i,j,:,:))
!
end do
Now we consider the inversion symmetry parameters
B leading to H'cv and
H'cvT
do j=1,3
select case(j)
case(1)
ind1=b1
case(2)
ind1=b2
case(3)
ind1=b3
end select
H'cv:
hamiltonM(1,3)=hamiltonM(1,3)+dcmplx(ind*ind1*rotationCond_wurM(i,j,1))
hamiltonM(1,4)=hamiltonM(1,4)+dcmplx(ind*ind1*rotationCond_wurM(i,j,2))
hamiltonM(1,5)=hamiltonM(1,5)+dcmplx(ind*ind1*rotationCond_wurM(i,j,3))
hamiltonM(2,6)=hamiltonM(2,6)+dcmplx(ind*ind1*rotationCond_wurM(i,j,1))
hamiltonM(2,7)=hamiltonM(2,7)+dcmplx(ind*ind1*rotationCond_wurM(i,j,2))
hamiltonM(2,8)=hamiltonM(2,8)+dcmplx(ind*ind1*rotationCond_wurM(i,j,3))
H'cvT:
hamiltonM(3,1)=hamiltonM(3,1)+ &
dconjg(dcmplx(ind*ind1*rotationCond_wurM(i,j,1)))
hamiltonM(4,1)=hamiltonM(4,1)+ &
dconjg(dcmplx(ind*ind1*rotationCond_wurM(i,j,2)))
hamiltonM(5,1)=hamiltonM(5,1)+ &
dconjg(dcmplx(ind*ind1*rotationCond_wurM(i,j,3)))
hamiltonM(6,2)=hamiltonM(6,2)+ &
dconjg(dcmplx(ind*ind1*rotationCond_wurM(i,j,1)))
hamiltonM(7,2)=hamiltonM(7,2)+ &
dconjg(dcmplx(ind*ind1*rotationCond_wurM(i,j,2)))
hamiltonM(8,2)=hamiltonM(8,2)+ &
dconjg(dcmplx(ind*ind1*rotationCond_wurM(i,j,3)))
end do
end do
Now we consider the P matrix elements that are
mixing the conduction and valence bands.
do i=1,3
select case(i)
case(1)
ind=k1
case(2)
ind=k2
case(3)
ind=k3
end select
do j=1,2
select case(j)
case(1)
ind1=p1
case(2)
ind1=p2
end select
H'cv:
hamiltonM(1,3)=hamiltonM(1,3)+ind*ind1*rotation_P_wurM(i,j,1)
hamiltonM(1,4)=hamiltonM(1,4)+ind*ind1*rotation_P_wurM(i,j,2)
hamiltonM(1,5)=hamiltonM(1,5)+ind*ind1*rotation_P_wurM(i,j,3)
hamiltonM(2,6)=hamiltonM(2,6)+ind*ind1*rotation_P_wurM(i,j,1)
hamiltonM(2,7)=hamiltonM(2,7)+ind*ind1*rotation_P_wurM(i,j,2)
hamiltonM(2,8)=hamiltonM(2,8)+ind*ind1*rotation_P_wurM(i,j,3)
H'cvT:
hamiltonM(3,1)=hamiltonM(3,1)+dconjg(ind*ind1*rotation_P_wurM(i,j,1))
hamiltonM(4,1)=hamiltonM(4,1)+dconjg(ind*ind1*rotation_P_wurM(i,j,2))
hamiltonM(5,1)=hamiltonM(5,1)+dconjg(ind*ind1*rotation_P_wurM(i,j,3))
hamiltonM(6,2)=hamiltonM(6,2)+dconjg(ind*ind1*rotation_P_wurM(i,j,1))
hamiltonM(7,2)=hamiltonM(7,2)+dconjg(ind*ind1*rotation_P_wurM(i,j,2))
hamiltonM(8,2)=hamiltonM(8,2)+dconjg(ind*ind1*rotation_P_wurM(i,j,3))
end do
end do
Now we consider the spin-orbit interaction and the crystal field splitting (crystal
field splitting = 0 for zincblende).
H'vv: hamiltonM(3:8,3:8) = hamiltonM(3:8,3:8) + rotationSpin_wurM
H'vv: hamiltonM(3:5,3:5) = hamiltonM(3:5,3:5) + crysSplit_wurM
H'vv: hamiltonM(6:8,6:8) = hamiltonM(6:8,6:8) + crysSplit_wurM
Inverse effective masses of electrons (S
parameters).
H'cc: hamiltonM(1,1)=hamiltonM(1,1)+ &
dcmplx(E_c+h2b2m_evAA2*s2*(k1**2+k2**2+k3**2))+
dcmplx(h2b2m_evAA2*(s1-s2)*(r13**2*k1**2+r23**2*k2**2+
r33**2*k3**2+2*r13*r23*k1*k2+2*r13*r33*k1*k3+2*r23*r33*k2*k3))
H'cc: hamiltonM(2,2)=hamiltonM(2,2)+ &
dcmplx(E_c+h2b2m_evAA2*s2*(k1**2+k2**2+k3**2))+
dcmplx(h2b2m_evAA2*(s1-s2)*(r13**2*k1**2+r23**2*k2**2+
r33**2*k3**2+2*r13*r23*k1*k2+2*r13*r33*k1*k3+2*r23*r33*k2*k3))
H'vv: hamiltonM(3,3)=hamiltonM(3,3)+dcmplx(E_v+h2b2m_evAA2*(k1**2+k2**2+k3**2))
H'vv: hamiltonM(4,4)=hamiltonM(4,4)+dcmplx(E_v+h2b2m_evAA2*(k1**2+k2**2+k3**2))
H'vv: hamiltonM(5,5)=hamiltonM(5,5)+dcmplx(E_v+h2b2m_evAA2*(k1**2+k2**2+k3**2))
H'vv: hamiltonM(6,6)=hamiltonM(6,6)+dcmplx(E_v+h2b2m_evAA2*(k1**2+k2**2+k3**2))
H'vv: hamiltonM(7,7)=hamiltonM(7,7)+dcmplx(E_v+h2b2m_evAA2*(k1**2+k2**2+k3**2))
H'vv: hamiltonM(8,8)=hamiltonM(8,8)+dcmplx(E_v+h2b2m_evAA2*(k1**2+k2**2+k3**2))
Now strain:
H'vv: hamiltonM(3:5,3:5)
= hamiltonM(3:5,3:5) + StrainHamiltonianVB_3x3M
H'vv: hamiltonM(6:8,6:8) = hamiltonM(6:8,6:8) +
StrainHamiltonianVB_3x3M
Final matrix must be hermitian.
end subroutine setup_kpbulk
The setup of the k.p Hamiltonian is done in subroutine
setup_kp_mat (setup_kp_mat1D).
!---------------------------------------------
subroutine
setup_kp_mat(num_qr,k_z,bndC)
!---------------------------------------------
! Creates k.p matrix for given region num_qr.
! MODULE sparse_matrix must be deallocated.
!
! num_qr = number of quantum region
! bndC = 'per','dir','neu'
!
!--------------------------------------------
Note: If Arnoldi is used the symmetry is
disturbed by adding a number disturb_arnoldi [MODULE:
numeric_control].
The discretized Hamiltonian can be written in the form D.H, where D
is a diagonal matrix and H is a hermitian one. It can be shown (see e.g.
in the
glossary) that the eigenvalues of the matrix D.H
are the same as
those of the matrix D1/2.H.D1/2. The matrix H is
stored as a sparse matrix in MODULE
sparse_matrix (or in 1D alternatively in MODULE
lapack_mat). The diagonal matrix D1/2 is stored as vector diag_kpV in MODULE
normalization_factors_sg. The
dimension is 6xdim_qr or 8xdim_qr for 6x6 or 8x8 k.p theory, respectively (dim_qr
is the dimension of the quantum region, i.e. the number of grid points on the
physical grid contained in the quantum region).
3D
The program loops through all grid points of the quantum region (1..dim_qr)
and gets the discretization of the differential operator for the corresponding point via subroutine
get_differentialkp3D. The
output (for grid point (i,j,k) on the physical grid) consists of
the following
components:
differentialM(n,8,8) , n=1..19 is specified as follows:
1 <--> (i , j , k-1)
2 <--> (i , j-1, k )
3 <--> (i-1, j , k )
4 <--> (i , j , k )
5 <--> (i+1, j , k )
6 <--> (i , j+1, k )
7 <--> (i , j , k+1)
8 <--> (i-1, j-1, k )
9 <--> (i+1, j-1, k )
10 <--> (i-1, j+1, k )
11 <--> (i+1, j+1, k )
12 <--> (i-1, j , k-1)
13 <--> (i+1, j , k-1)
14 <--> (i-1, j , k+1)
15 <--> (i+1, j , k+1)
16 <--> (i , j-1, k-1)
17 <--> (i , j+1, k-1)
18 <--> (i , j-1, k+1)
19 <--> (i , j+1, k+1)
So for each of the 19 adjoining grid points on the Physical mesh an 8x8 block is
provided which contains the discretized Hamiltonian on this point. For
6x6 k.p theory only the (3:8,3:8) submatrix is valid.- In order to write these elements to the matrix, one must know the position
of the points
n=1..19 within the numbering in the quantum region (1..dim_qr).
This information is provided in coord_in_qrM(1:19).
- The component in the vector
diag_kpV for the regarded grid point is
provided by diag_kp.
In principle the setup is done twice. The first time only the nonvanishing
matrix elements are counted in order to be able to allocate the storage space
for the sparse matrix.
Then the real setup is done.
See the source code for further information (setup_kp_mat).
In 3D, the subroutine
get_differentialkp3D first gets the k.p parameters and then discretizes the
differentials.
!---------------------------------------------------------------------
subroutine
get_differentialkp3D(iqr,num_qr,bndC,differentialM,coord_in_qrM,diag_kp)
!---------------------------------------------------------------------
!
! Input:
! iqr --> num_qr
!
! bndC = 'neu' or 'dir'
!
! Output:
!
! differentialM(1..19,8,8)
!
!
! coord_in_qrM(1:19) contains for each point in differentialM
! the position within the quantum region and the matrix.
!
! coord_in_qrM = 0 if the corresponding point does not lie within quantum region.
!
! Note: For 6x6 k.p theory only the (3:8,3:8)
submatrix is valid.
! diag_kp diag_kpV for the point
iqr.
!
!----------------------------------------------------------------------
First, the routine gets the k.p parameters: call
Get_kp_MaterialParameters.
(More
details here.)
Then the dicretization routine is called:
call
discretize_kp3D
Note:
The hermitian matrix D^{1/2} H D^{1/2} is assembled. D is
returned in diag_kp.- The differential operator is set up in
!----------------------------------------------------------------------
subroutine
discretize_kp3D(i_phys,j_phys,k_phys,adr,hV,parV,
&
par_betweenM,neuVL,differentialM)
!----------------------------------------------------------------------
!
! par_betweenM(1:6,1:14):
! -----------------------
! par_betweenM(n_par,1:14) <--> n_par=1 (i ,j ,k-1)
! n_par=2 (i ,j-1,k )
! n_par=3 (i-1,j ,k )
! n_par=4 (i+1,j ,k )
! n_par=5 (i ,j+1,k )
! n_par=6 (i ,j ,k+1)
! parV(1:19)
!
! parV, par_minV and
par_andV see: Information on how to extract information on electronic structure
!--------------------------------------------
!
!
! --------
! differentialM(1..19,8,8)
!
!
!--------------------------------------------
!
! neuVL (1) <--> i-
! (2) <--> i+
! (3) <--> j-
! (4) <--> j+
! (5) <--> k-
! (6) <--> k+
!
! = .TRUE.
!
!-----------------------------------------------------------------
2D
2D is analogous to 3D. See also source code for further information (setup_kp_mat).
Note: Discretization is done in x- and y-plane. kz must be provided via input
(1/AA).
In 2D, the subroutine
get_differentialkp2D first gets the k.p parameters and then discretizes the
differentials.
The program loops through all grid points of the quantum region (1..dim_qr)
and gets the discretization for the corresponding point via subroutine
get_differentialkp2D. The
output (for grid point (i,j) on the physical grid) consists of
the following
components (same scheme as in 3D but k-component discarded):
differentialM(n,8,8) , n=1..19 is specified as follows:
1 <--> ---------
2 <--> (i , j-1)
3 <--> (i-1, j )
4 <--> (i , j )
5 <--> (i+1, j )
6 <--> (i , j+1)
7 <--> ---------
8 <--> (i-1, j-1)
9 <--> (i+1, j-1)
10 <--> (i-1, j+1)
11 <--> (i+1, j+1)
12 <--> ---------
13 <--> ---------
14 <--> ---------
15 <--> ---------
16 <--> ---------
17 <--> ---------
18 <--> ---------
19 <--> ---------
So for each of the 9 (3D: 19) adjoining grid points on the physical mesh an 8x8 block is
provided which contains the discretized Hamiltonian on this point. For
6x6 k.p theory only the (3:8,3:8) submatrix is valid.
!---------------------------------------------------------------------
subroutine
get_differentialkp2D(iqr,num_qr,k_z,bndC,differentialM,
coord_in_qrM,diag_kp)
!---------------------------------------------------------------------
!
! discretizes in (i,j) plane,
! k-component discarded
!
!--------------------
! k_z [1/AA]
!---------------------------------------------------------------------!---------------------------------------------------------------------
subroutine
discretize_kp2D(i_phys,j_phys,adr,k_z,hV,parV,
&
par_betweenM,neuVL,differentialM)
!---------------------------------------------------------------------
!!
! par_betweenM(1:6,1:14) :
!------------------------
! par_betweenM(n_par,1:14) <--> n_par=1 ---------
! n_par=2 (i ,j-1)
! n_par=3 (i-1,j
)
! n_par=4 (i+1,j
)
! n_par=5 (i ,j+1)
! n_par=6 ---------
!
5
!
3 4
!
2
!
! parV(1:19)
! parV, par_minV and
par_andV see: Information on how to extract information on electronic structure
!------------------------------------------------------
!
!
!--------------
! differentialM(1..19,8,8)
!
!--------------------------------------------
!
! neuVL (1) <--> i-
! (2) <--> i+
! (3) <--> j-
! (4) <--> j+
! (5) <--> ------------
! (6) <--> ------------
!
! = .TRUE.
!
!-----------------------------------------------------------------
1D
In 1D there is the additional option to assemble a full matrix (instead of a
sparse matrix) which can be diagonalized by LAPACK.
The k.p matrix is assembled in the following steps:
- Loop through all points in quantum region (
1..dim_qr).
- Get the k.p parameters for this point via subroutine
Get_kp_MaterialParameters.
- Get the discretization of the Hamiltonian (differential star): subroutine
discretize_kp1D.
In 1D, only the two neighboring points are involved in discretization --> so
the discretization on point iqr(=1..dim_qr) can be expressed with three
8x8 matrices leftM, onM, rightM. Since the matrix is hermitian, only onM
and
rightM are provided.
- Counts off-diagonal elements (upper triangle matrix) different from 0
- Allocates matrices: either for LAPACK (i.e.
full matrix) or for iterative method (i.e.
sparse matrix)
- Does the same loop as under i., ii. and iii. and defines k.p matrix
under iv.
- Write the nonvanishing matrix elements of
onM and rightM to the
k.p matrix.
See also source code for further information (setup_kp_mat1D).
-
!-------------------------------------------------------------
subroutine
setup_kp_mat1D(num_qr,bndC,kx,ky,K_z,matrix_formC)
!-------------------------------------------------------------
!
!
for given region num_reg.
!
MODULE sparse_matrix
must be
deallocated.
!
!
ija, sa
in MODULE sparse_matrix
if matrix_formC = 'iterat'
! kp_lapackM
in MODULE lapack_mat
if matrix_formC = 'lapack'.
!
!
!
!
!
! num_qr =
! bndC = 'per','dir','neu'
!
! kx,ky,K_z
!
! matrix_formC = 'lapack' -->
MODULE lapack_mat
!
'iterat' -->
sparse matrix
!
!
! kx,ky
! K_z
!
! norm_kpV (MODULE normalization_factors_sg)
!
!-------------------------------------------------------------
Note: If Arnoldi is used the symmetry is
disturbed by adding a number disturb_arnoldi [MODULE
numeric_control].
The differential operator in discretized form
for a given point on the physical grid is provided by subroutine
discretize_kp1D in
MODULE discretize_kp_1D.
Note:
The hermitian matrix D^{1/2} H D^{1/2} is assembled.
First the rotation is updated (see
Rotation of k.p Hamiltonian). Then the
rotation matrices for crystal-splitting and spin-orbit-coupling are defined.
Finally the strain k.p matrix for the considered point is set up:
call
update_rotation(adr)
d1=parV(17)
d2=parV(18)
d3=parV(19)
call
setup_crysSplit_wurM(d1)
call
setup_rotationSpin_wurM(d2,d3)
call
set_strain_mat_physgrid_point(i_phys,1,1)
!------------------------------------------------------------------------
subroutine
discretize_kp1D(K_z,i_phys,adr,h_imin,h_i,kx,ky,par_minV,parV, &
par_andV,onM,rightM)
!------------------------------------------------------------------------
!
!
i=i_phys.
!
In any case, the 8x8-Hamiltonian is discretized.
!
If only 6x6-kp is done, only the part (3:8,3:8)
is relevant.
!
In this version, the superlattice vector K_z
is
also taken into account.
!
If no superlattice --> K_z must be 0.
!
The discretization provides a diagonal (8x8)
block [called onM here]
!
and two off-diagonal (8x8)-blocks.
!
Since total matrix is hermitian, only the
off-diagonal block to the right
!
[called rightM here] is calculated:
!
!
i-1 --- -- hi-1 ----- i ----- hi
----- i+1
!
i-1 --- -- h_imin ----- i ----- h_i ----- i+1
!
par_minV parV par_andV
!
!
onM rightM
!
(8x8) (8x8)
!
!
!
!
H = RM1 kx2 + RM2 ky2
+ RM3 kz2 + RM4 kxky
+ RM5 kxkz + RM6 kykz
!
!
!
!
1) subroutine
put_together(num_diff,coeff_lM,coeff_rM)
!
! num_diff = 1 --> kx2
!
= 2 --> ky2
!
= 3 --> kz2
!
= 4 --> kxky
!
= 5 --> kxkz
!
= 6 --> kykz
!
= 7 --> constant
!
= 8 --> kx
!
= 9 --> ky
!
= 10 --> kz
!
!
Creates corresponding coefficient matrix out of
par_andV
and par_minV
for the specified differential.
!
coeff_lM par_minV
is valid, for coeff_rM par_andV.
!
!
2) subroutine put_together_on(num_diff,coeff_M)
!
! num_diff = as above
!
!
parV for the
specified differential.
!
!
put_together
and put_together_on
are in principal identical.
!
1..10.
!
!
grid point is used, sometimes the value between two points.
!
So the values on and between the points are set up
!
separately and delivered to the subroutines that setup the
differentials.
!
!
It holds the relation (with coeff_iM
from
put_together_on(i,coeff_iM)
!
and num_diff_i
as defined above):
!
! H = sum_i
[ coeff_iM num_diff_i ] (i
1
to 10)
!
!------------------------------------------------------------------------
!
!
For kx,ky
the values from input are directly used.
!
! kz
kz^2
are discretized:
!
!
1) subroutine
set_kzkz(h_imin,h_i,mat_onM,mat_lM,mat_rM,on_1M,right_1M)
!
!
kzkz
!
on input:
mat_lM -> i-1/2 , mat_rM -> i+1/2 , mat_onM -> i
! kz -->
-i d/dz
!
!
----------------------------------------------------
! (kz)^2 = -d^2/dz^2
+ 2K_z kz + K_z^2 (K_z
superlattice vector,
!
=0
!
----------------------------------------------------
!
mat_onM,mat_lM,mat_rM -->
!
on_1M,right_1M -->
!
!
!
!2) subroutine
set_kz(h_imin,h_i,mat_onM,mat_lM,mat_rM,on_2M,right_2M)
!
!
kz
!
on input:
mat_lM -> i-1/2 , mat_rM -> i+1/2 , mat_onM -> i
! kz -->
-i d/dz
!
!
-----------------------------------------------------
! kz = -i d/dz
+ K_z
!
-----------------------------------------------------
!
!
mat_onM,mat_lM,mat_rM -->
input
!
on_2M,right_2M -->
!
!------------------------------------------------------------------------------
!
!
onM -->
!
rightM -->
!
!
i_phys
(1..n_1)
!
adr
is point number on grid (1..lengvecV).
!
-->
.
!
h_imin
and h_i
distance to neighboring point [in
AA] !!!!!!
!
kx,ky,K_z
1/AA !!!!
!
! par_minV -->
i-1
and i
! par_andV -->
material parameters between i
and i+1
! parV -->
material parameters on i
!
! par_andV( 1) = L1
! par_andV( 2) = L2
! par_andV( 3) = M1
! par_andV( 4) = M2
! par_andV( 5) = M3
! par_andV( 6) = N1
! par_andV( 7) = N2
! par_andV( 8) = B1
! par_andV( 9) = B2
! par_andV(10) = B3
! par_andV(11) = P1
! par_andV(12) = P2
! par_andV(13) = s1
! par_andV(14) = s2
!
!
only for parV
!
! parV(15) = E_c
! parV(16) = E_v
! parV(17) = d1
! parV(18) = d2
! parV(19) = d3
!
!
For definiton of parV, par_minV and par_andV
see: Information on how to extract information on electronic structure
!
!
eV , AA
!
! neumann : flag --> = 0
(normal: inner point or periodic or Dirichlet boundary)
!
= 1
!
= 2
!
!------------------------------------------------------------------------------
The subroutine
update_kp
calculates the new k.p states for the potential phiV, which is provided as
input. Additional input is whether to calculate it for electrons (el)
or holes (hl) and the number of the quantum region (num_qr).
The following steps are carried out:
Store potential for k.p parameters
Since the k.p parameters Ec and Ev must include
the electrostatic potential, the array
phiV in MODULE
potentials is taken into account (see
subroutine Get_kp_MaterialParameters
in How to extract information on electronic structure).
Determine energy where to search for electron/hole states
The internal
subroutine
prepare_det_edge determines the values E_cond and E_val. E_cond
is the
lowest conduction band value (of Gamma band which is assumed to have band number
1) throughout the quantum region and E_val is the highest valence band value (of band numbers 1 to 3 = hh, lh, split-off hole).
From E_cond and E_val the energy level E_separate
is calculated.
NOTE: This model treats electrons and holes as strictly separated. This allows
to combine different models for electrons and holes. This model fails,
however, if E_cond<E_val. Then the program stops. In order
to handle such a situation within 8x8 k.p theory a model must be implemented that does not
distinguish between electron and hole states.
If n eigenvalues are requested for the conduction band, the
subroutine
calculate_kp_eigenvalues looks for the n eigenvalues
above E_separate (below E_separate for the valence
band respectively).
Solve eigenvalue problem for k|| = 0
This is done in two steps:
- Setup k.p matrix in subroutine
setup_kp_mat1D (see Setup of k.p Hamiltonian).
- Calculate the specified eigenvalues of that matrix in subroutine
calculate_kp_eigenvalues (see
Calculate k.p eigenvalues).
Determine separation energy (SeparationEnergyEdgeV) quantum-mechanical <--> classical
density anew
Now the separation energy is determined anew in subroutine
determine_edge_qm.
This subroutine uses pot_condM / pot_valM (el/hl), which are generated by the
internal subroutine
prepare_det_edge,
as input (see Separation models:
Eigenvalues <--> Thomas-Fermi).
Problem of defining k||
The range of k|| is determined in subroutine
determine_range_for_kpar. Input is the energy of the ground
eigenstate for k||=0. Furthermore, the maximum energy which separates
classical and quantum mechanical treatment (which is known from subroutine
determine_edge_qm (which was described in the paragraph above) is provided.
This subroutine is structured as follows:
- The Fermi level in the quantum region is retrieved for electrons or holes,
respectively.
- The energy is calculated up to which the eigenstates contribute to the
density. This energy depends on the Fermi level.
- The minimum/maximum value (el/hl) of this energy and the separation energy
(input) is finally the energy which must be reached by k|| states.
- In the following the program employs the bulk dispersion (in directions
100, 010 and 110. Note: Quantization in the rotated system is always in
z direction) in order to determine k||. In this algorithm electrons/holes
are supposed to have an upward/downward dispersion. If this is not the case
due to the k.p parameters, the algorithm fails.
Discretize Brillouin zone
The k points are provided by subroutine
get_k_par1D. Here the k points are retrieved by calling subroutine get_dim_2D_BZ.
Details on it can be found on the documentation of the
Brillouin zone integration.
Note: Since subroutine
get_dim_2D_BZ discretizes a Brillouin zone with kx and ky varying from 0 to
1.0/AA, the k vectors must be scaled to k_parallel_max accordingly.
Redefine number of eigenvalues to be calculated
It may turn out that fewer eigenvalues are necessary for the density
calculation than the specified number. Then the number of eigenvalues to be
calculated is reduced.
Calculate eigenstates for nonvanishing k vectors
This is done in the internal subroutine
calculate_kp_par.
The subroutine
setup_sg defines the Schrödinger matrix H' in MODULE sparse_matrix_real
(i.e. as a
sparse matrix) as well as the diagonal matrix D in diag_sgV (which
is stored as a vector) in MODULE
normalization_factors_sg. Note that the whole
eigenvalue problem is the solution of H: H=D.H'
The subroutine
setup_sg (num_qr,num_sg,num_deg,chargeC,phiV,max_edge)
sets up
the Schrödinger equation, specified by:
Input:
num_qr ... Number of quantum cluster
num_sg ... In a quantum cluster there are several
Schrödinger equations to be solved.
First this set of equations can be distinguished with
respect to the energy edge
which is valid for the corresponding equation.
(Gamma, X, L, heavy hole, light hole, split-off hole; including splitting of
them due to strain.)
-> num_sg ... specifies the band edge.
num_deg ... For the same energy edge, there are several Schrödinger
equations to be solved
with different masses
(-> different mass tensor, e.g. due to anisotropic mass tensor) or boundary conditions.
-> num_deg ... specifies one of them.
(See MODULE
quantum_models for further explanation.)
chargeC = 'el' , 'hl' (electrons or holes)
phiV ... Electrostatic potential which is valid for the setup.
As a result, the matrix in MODULE sparse_matrix_real is allocated and contains the
Schrödinger matrix.
MODULE
normalization_factors_sg: diag_sgV are allocated and defined.
Output:
max_edge : Minimum (maximum) band edge energy that enters the Schrödinger equation for electrons (holes).
--> Eigenvalues lie above (below) max_edge for electrons (holes).
The routine basically consists of the following steps:
- Loop through all points of the quantum cluster (
1..dim_qr).
- For each point extract the parameters needed for the setup of the differential
at that point (effective masses and band edge):
CALL get_parameters_for_setup
- With these parameters define computational molecule (finite difference
scheme):
CALL set_differentials_sg
Here:
valV(1..4) Matrix elements for the following points (of computational molecule):
With (i,j,k) the coordinates on the physical grid of the point under
consideration, it holds
for the components of valV:
1 --> diagonal (i,j,k)
2 --> off-diagonal for i+1
3 --> off-diagonal for j+1
4 --> off-diagonal for k+1
coordM(1..4,2) --> Coordinates (row, column) in the matrix for the
corresponding matrix elements (1..4) defined above.
- Write nonvanishing elements to the matrix.
Internal details:
The routine uses several internal auxiliary subroutines:
SUBROUTINE neighbor_points_sg3D:
Routine that provides grid information on
neighboring points.SUBROUTINE write_sg_masses:
Effective masses are stored on material grid,
for setup of Hamiltonian matrix they are used very often.
--> Store them on an auxiliary variable:
mass_sgM(1..dim_qr,1..8,1..3) contains for each point within
the quantum cluster (on the physical grid) [1..dim_qr]
for all 8 adjoining octants (which can be ascribed to
one of material grid points) [1..8] the effective masses
in x-,y- and z-direction.SUBROUTINE get_mass_qr:
Provides mass for Schrödinger equation
between the points num1 and num2.SUBROUTINE get_parameter_sg:
Auxiliary to extract parameters (energy
edge, effective mass).FUNCTION octant_in_qr:
Auxiliary subroutine to decide if adjoining points
are within/outside quantum cluster. It is used to handle boundary conditions.SUBROUTINE set_differentials_sg:
Only the part of the differential for the upper offdiagonal matrix is
assembled since the Hamiltonian is symmetric.
1D
In 1D, the Schrödinger equation is solved with LAPACK. Consequently, a full
matrix is defined and it is solved immediately. This is done in subroutine
SG_1D. The parameters
effective mass, energy
edge and grid spacing must be provided as input. Eigenvalues and eigenfunctions
are provided as output.
Details:
Setup of Schrödinger equation
General hints on Schrödinger equation
The data which enter the Schrödinger equation
are just from grid points which lie within quantum clusters.
Quantum clusters are primarily defined on material
grid.
So on boundary points one has to decide which
octants are valid and which are not.
This is done in the internal subroutine
octants_in_qr.
!==============================================================
The Schrödinger matrix is defined in
SUBROUTINE
setup_sg.
[In 1D it is defined and solved in
SUBROUTINE
SG_1D].
The matrix is written to MODULE
sparse_matrix_rout_real.
1) First the masses for the Schrödinger equation
are extracted and written to mass_sgM.
mass_sgM(1..dim_qr,1..8,1..3) contains for each point of the quantum
cluster (1..dim_qr) and
each octant(1..8) the masses in x-,y- and
z-direction (1..3).
This is done in
internal
SUBROUTINE write_sg_masses:
dim_qr=dim_qrV(num_qr)
! dim_qr Dimension of quantum region
allocate(mass_sgM(dim_qr,8,3),STAT=istatus)
if (istatus/=0) stop 'error write_sg_masses: alloc.
failed'
mass_sgM=0.0d0
num_guent=map_qr_guentherV(num_qr)%number_of_points
! num_guent number of material grid points which defined this quantum region
! note: quantum regions originally defined on material grid
do ncount=1,num_guent
! noe loop through all the points
i_j_k=map_qr_guentherV(num_qr)%guenther_grid_points(ncount)
! i_j_k position on material grid
call
convert_guenther_to_phys3D(i_j_k,coordM)
! coordM(8,3) contains in (n,1:3) the coordinates (i,j,k) of the 8 (n=1..8)
adjoining points
! on the Physical grid.
i=coordM(1,1)
j=coordM(1,2)
! now treat first of these 8 points: (i,j,k)
k=coordM(1,3)
iqr=quantumreg_numM(i,j,k)
! iqr is position of (i,j,k) in quantum region (1..dim_qr)
nqr=quantum_gridM(i,j,k)
! nqr is number of quantum region, in which point (i,j,k) lies.
if (nqr==num_qr) then
! if this point lies within this quantum region, the corresponding octant (1..8)
OF THIS POINT
select case(chargeC)
! is ascribed the corresponding mass
case('el')
mass_sgM(iqr,8,:)=def_quantum_modelsM(num_qr)%
&
data_sg_elM(num_sg)%mV(num_deg,ncount,:)
case('hl')
mass_sgM(iqr,8,:)=def_quantum_modelsM(num_qr)%
&
data_sg_hlM(num_sg)%mV(num_deg,ncount,:)
end select
end if
.........
the same for the remaining 7 points
end do
!-----------------------------------------------------------------
do i=1,dim_qr
if (all(mass_sgM(i,:,:)==0.0d0))
stop 'error write_sg_masses: &
! check: if all zero --> internal error
& grid point without mass within quantum region'
if (any(mass_sgM(i,:,:)==0.0d0))
then
! note: on boundaries there may be octants which do not lie within the quantum
region any more
massV=0.0d0
! then they are ascribed the
average masses of the other octants
su=0
do n=1,8
if (any(mass_sgM(i,n,:)/=0.0d0))
then
su=su+1
massV=massV+mass_sgM(i,n,:)
end if
end do
massV=massV/dble(su)
do n=1,8
if (any(mass_sgM(i,n,:)==0.0d0))
mass_sgM(i,n,:)=massV
end do
end if
end do
===============================================================
2) There are several auxiliary internal
subroutines:
2.1) internal subroutine
neighbor_points_sg3D:
!------------------------------------------------------------------------------
subroutine
neighbor_points_sg3D(num_qr,num_sg,num_deg,chargeC,num,num_max, &
i1,i2,j1,j2,k1,k2,hi1,hi2,hj1,hj2,hk1,hk2)
!-----------------------------------------------------------------------------
input:
num is the position of the mesh
point in the whole vector
of the
schroedinger equation,i.e. in 1..dim_qrV
num_max is the number of mesh
points=dim_qrV(num_qr)
num_sg ... number of
schroedinger equation with different energy levels
num_deg ... number of
schroedinger equation which belongs to that
energy level
chargeC = 'el' or 'hl'
output:
positions of the neighboring
points in the whole vector (1..dim_qr)
- one step back in x-direction
: i1 , distance = hi1
- one step forward in
x-direction : i2 , distance = hi2
- one step back in y-direction
: j1 , distance = hj1
- one step forward in
y-direction : j2 , distance = hj2
- one step back in z-direction
: k1 , distance = hk1
- one step forward in
z-direction : k2 , distance = hk2
for boundary points
- periodic: number of point to
which it is connected
- Neumann <==>
value=0 , then distance =0
- Dirichlet <==> value=-1 ,
then distance =0
it must be:
num,i1,i2,j1,j2,k1,k2 =<
num_max (=dim_qr)
num >= 1
problem must be symmetric, i.e.
if num=a --> i1=b then
also num=b --> i2=a and so on
!-----------------------------------------------------------------
units: lengths in AA
!------------------------------------------------------------------------
So the subroutine provides the
positions of the adjoining points in the Schrödinger matrix.
2.2) internal subroutine
get_mass_qr(num_qr,num1,num2,dirC,mass)
!------------------------------------------------------------------
subroutine get_mass_qr(num_qr,num1,num2,dirC,mass)
!------------------------------------------------------------------
! provides mass for schroedinger equation between
the points
! num1 and num2
!
! input:
!
! num_qr ... number of quantum
region
! dirC ... 'x','y','z'
!
! output:
!
! mass ... effective mass
in specified direction
!
for electrons positive,
!
for holes positive
!
!-------------------------------------------------------------------
i1=inv_quantumreg_numM(num_qr,num1,1)
j1=inv_quantumreg_numM(num_qr,num1,2)
k1=inv_quantumreg_numM(num_qr,num1,3)
i2=inv_quantumreg_numM(num_qr,num2,1)
j2=inv_quantumreg_numM(num_qr,num2,2)
! (i1,j1,k1),(i2,j2,k2) coordinates of corresponding points
k2=inv_quantumreg_numM(num_qr,num2,3)
! on Physical grid
indV(1)=i2-i1
indV(2)=j2-j1
indV(3)=k2-k1
if (.not.((any(abs(indV)==1)).and.(sum(abs(indV))==1)))
&
stop 'error in subroutine
get_mass_qr: in this version just adjoining points are allowed'
locV = maxloc(abs(indV))
! locV(1)= 1 --> num1 --- x-direction --- num2
select case (locV(1))
case (1)
if (dirC/='x')
stop 'error get_mass_qr: internal inconsistency'
l_octV(1)=2
! get corresponding octants
l_octV(2)=4
l_octV(3)=6
l_octV(4)=8
r_octV(1)=1
r_octV(2)=3
r_octV(3)=5
r_octV(4)=7
if
(indV(1)==-1) then
! if i2<i1 swap
z_octV=r_octV
r_octV=l_octV
l_octV=z_octV
end if
l_massV = (/
(mass_sgM(num1,l_octV(l),1) , l=1,4 ) /)
! get left and right mass
r_massV = (/
(mass_sgM(num2,r_octV(l),1) , l=1,4 ) /)
case (2)
! locV(1)= 2 --> num1 --- y-direction --- num2
! locV(1)= 3 --> num1 --- z-direction --- num2
... end select
!--------------------------------------------
mass = (sum(l_massV)+sum(r_massV))/8d0
! take the average
end subroutine get_mass_qr
!------------------------------------------------------------------
2.3) internal function octant_in_qr
!------------------------------------------------------------------
function octant_in_qr(i,j,k,oct,num_qr)
result(inL)
!------------------------------------------------------------------
! --> true if all adjoining points on physical
grid are within
! quantum region
num_qr
!------------------------------------------------------------------
Note: for point (i,j,k) and oct=1 (i.e. octant i-,j-,k-)
the adjoining points are:
(i-1,j-1,k-1),(i,j-1,k-1),(i-1,j,k-1),(i,j,k-1),(i-1,j-1,k),(i,j-1,k),(i-1,j,k),(i,j,k)
2.4) internal subroutine
get_parameter_sg:
!-------------------------------------------------------------------
subroutine
get_parameter_sg(iwhich,whereC,num1,num2,val,num_qr, &
num_sg,chargeC)
!-------------------------------------------------------------------
! input:
! iwhich = 1
==> E_edge
!
= 2 ==> m_eff_x
!
= 3 ==> m_eff_y
!
= 4 ==> m_eff_z
!
! whereC = 'on' ==>
val = value on point with number num1
!
in this case num2 is of no importance
!
!
= 'be' ==> val = value between points
!
with number num1 and num2
!
! output:
! val value
returned
!
! note:
! it always must be : 1 =< num1 , num2
=< n_sg3D
!
! num_qr ... number of quantum
region
! chargeC ... = 'el' or 'hl'
! mass is negative for holes !!!
!
! units : eV and electron_mass
!--------------------------------------------------------------------
select case(whereC)
case('on')
i=inv_quantumreg_numM(num_qr,num1,1)
j=inv_quantumreg_numM(num_qr,num1,2)
k=inv_quantumreg_numM(num_qr,num1,3)
! (i,j,k) of the point on Physical grid
select case(chargeC)
case('el')
select case(iwhich)
case(1)
inVL=.false.
call
get_addresses_from_phys3D(i,j,k,oct1V)
! positions of the eight octants on Stefan's grid
do n=1,8
mat=mat_numberV(oct1V(n))
! material number in numbering on whole device
num_mat=map_qr_to_materialV(mat,num_qr)
! material number in numbering in quantum region
num_band=def_quantum_modelsM(num_qr)%data_sg_elM(num_sg)%bandnumV(num_mat)
! get corresponding band parameter
vec1V(n)=E_cM(num_band,oct1V(n))/abs(electron_charge)-phiV(oct1V(n))
inVL(n)=octant_in_qr(i,j,k,n,num_qr)
end do
val=0.0d0
su=0
do n=1,8
! take the average of all bandedges which are within quantum region
if (inVL(n)) then
su=su+1
val=val+(vec1V(n))
end if
end do
if (su==0) stop 'error
get_parameter_sg: internal error'
val=val/dble(su)
case(2,3,4)
stop 'error in
get_parameter_sg: masses just between points implemented'
case default
stop 'error in
get_parameter_sg: iwhichC wrong'
end select
.................
case('be')
select case(chargeC)
case('el')
select case(iwhich)
case(1)
stop 'error in
get_parameter_sg: energies just on points implemented'
case(2)
call
get_mass_qr(num_qr,num1,num2,'x',val)
case(3)
call
get_mass_qr(num_qr,num1,num2,'y',val)
case(4)
call
get_mass_qr(num_qr,num1,num2,'z',val)
case default
stop 'error in
get_parameter_sg: iwhichC wrong'
end select
...........
note : mass negative for holes
2.5) internal subroutine
get_parameters_for_setup
!-----------------------------------------------------------------
subroutine
get_parameters_for_setup(num_qr,num_sg,num_deg,dim_qr,&
iqr,chargeC,neighborpointsV,dist_neighborV,kind_of_pointV,&
E_edge,m_effV)
!------------------------------------------------------------------
!
! input:
! num_qr,num_sg,num_deg ... see above
! dim_qr ... dimension of quantum
region
! iqr ... number of point in
quantum region
! chargeC='el' , 'hl'
!
! output:
! neighborpointsV(1:6) =
(i1,i2,j1,j2,k1,k2)
! dist_neighborV(1:6) =
(hi1,hi2,hj1,hj2,hk1,hk2) in [AA]
! kind_of_pointV(1:3) = 'i' --> inner
point
!
'n' --> neumann point
!
'd' --> dirichlet point
!
note: for periodic boundary conditions
!
inner point must be specified
! note: 1 --> x-direction
(i+1)
!
2 --> y-direction (j+1)
!
3 --> z-direction (k+1)
!
! E_edge ... energy in [eV] on (i,j,k)
! m_effV(1..6) effective masses
in [electron_mass]
!
(mi1,mi2,mj1,mj2,mk1,mk2)
!
!---------------------------------------------------------
!
!
i1
!
!
|
!
hi1 , mi1
!
|
!
!
j1 -hj1- (i,j,k) -hj2- j2
!
mj1 mj2
'i'
!
|
!
hi2 , mi2
!
|
!
!
i2 'i'
!
!------------------------------------------------------------------
================================================================
3) internal subroutine
set_differentials_sg:
provides differential (difference star)
for Schrödinger equation:
Note: Matrix is symmetric, so only
off-diagonal elements for upper half are calculated
!------------------------------------------------------------------
subroutine
set_differentials_sg(iqr,neighborpointsV,dist_neighborV,&
kind_of_pointV,E_edge,m_effV,valV,coordM)
!------------------------------------------------------------------
!
! input:
! iqr .. = row in matrix
! neighborpointsV(1:6) =
(i1,i2,j1,j2,k1,k2)
! dist_neighborV(1:6) =
(hi1,hi2,hj1,hj2,hk1,hk2) in [AA]
! kind_of_pointV(1:3) = 'i'
--> inner point
!
'n' --> neumann point
!
'd' --> dirichlet point
! note: for periodic
boundary conditions inner point must be specified
! note: 1 --> x-direction
(i+1)
!
2 --> y-direction (j+1)
!
3 --> z-direction (k+1)
!
! E_edge ... energy in [eV] on (i,j,k)
! m_effV(1..6) effective masses
in [electron_mass]
!
(mi1,mi2,mj1,mj2,mk1,mk2)
!
! output:
! valV(1..4)
!
1 --> diagonal
!
2 --> off-diagonal for i+1
!
3 --> off-diagonal for j+1
!
4 --> off-diagonal for k+1
!
! coordM(1..4,2) --> (row,column)
of corresponding element
!------------------------------------------------------------------
4) main subroutine:
!-------- define matrix
--------------------------------
call setup_sparse_matrix_real(dim_qr,num_points)
do iqr=1,dim_qr
! loops through all points
call get_parameters_for_setup(num_qr,num_sg,num_deg,dim_qr,&
! gets parameters
iqr,chargeC,neighborpointsV,dist_neighborV,kind_of_pointV,&
E_edge,m_effV)
call
set_differentials_sg(iqr,neighborpointsV,dist_neighborV,&
! calculates differentials
kind_of_pointV,E_edge,m_effV,valV,coordM)
do n=1,4
! write matrix elements
if (valV(n)/=0.0d0) then
row=coordM(n,1)
column=coordM(n,2)
value=valV(n)
call
add_element_sparse_real(row,column,value,.TRUE.)
end if
end do
end do
!--------- allocate normalization and
trafovector ---------------
allocates and defines diag_sgV in
MODULE
normalization_factors_sg.
Separation between quantum and continuum states
The separation energy is the energy level that separates quantized and
continuum states. See the
glossary for more details.
The information about these models is internally stored in MODULE
quantum_solutions.
The separation energy is determined in SUBROUTINE
determine_edge_qm.
SUBROUTINE determine_edge_qm(dim_qr,potM,num_ev,eigenvalV,chargeC,SeparationEnergyEdgeV,SeparationModelC)
For given potential
- potM(1..dim_qr,1:2) (1D)
- potM(1..dim_qr,1:4) (2D)
- potM(1..dim_qr,1:8) (3D)
eigenvalV(1..num_ev) [where eigenvalV(1..num_ev) are used for
the density calculation)], SeparationEnergyEdgeV (i.e. the energy which separates classical and quantum mechanical
density) is determined.
chargeC = 'el' --> eigenvaluesV must be given in ascending order
= 'hl' --> eigenvaluesV must be given in descending order
potM contains for each point 1..dim_qr the potental value for each octant [potM(.,1:8)]
(i.e. the conduction (or valence) band edge minus the electrostatic potential).
1 2 12 13 14 25
26 ... dim_qr
* * ... *
* * ... *
* ... [grid points]
----------------| |----------- energy = 10.0
|
|
|
|
|
|
|----------------------|
energy = 0.0
potM(1,1) potM(12,1) potM(25,1)
= 10.0 = 10.0 = 0.0
potM(1,2) potM(12,2)
potM(25,2)
= 10.0 = 0.0
= 10.0
potM(13,1)
= 0.0
potM(13,2)
= 0.0
Note: Ordering of eigenvalues
electrons: ev(1) = lowest eigenvalue
holes: ev(1) = highest eigenvalue
SeparationModelC
... 'eigenvalue' -->
SeparationEnergyEdgeV = topmost eigenvalue (num_ev)
... 'max_energy' -->
SeparationEnergyEdgeV =
lowest energy point i this region
+ SeparationEnergyEdgeV
[Here SeparationEnergyEdgeV as input and output.]
(for electrons)
... 'edge_model' -->
SeparationEnergyEdgeV is the highest eigenvalue that is below minimum of
conduction band edges
on left and right boundary (1D) (for
electrons)
For holes accordingly.
More details on these models can be found in the
Glossary.
All potentials must have the same units (they are just compared relative to each other).
After the matrix has been defined, the k.p eigenvalues are calculated in
subroutine calculate_kp_eigenvalues (MODULE
mod_calculate_kp_eigenvalues).
Up to now three methods are implemented:
- Conjugate gradient method of J. Majewski
--> Can only be used for extreme
eigenvalues.
Iterative solver using conjugate gradient method
(up to now only for holes in 6x6 k.p).
- Arnoldi (ARPACK package)
- Jacobi-Davidson (default method)
In 1D, there is additionally LAPACK (which is the default method). The method,
which is used, is specified in SchroedingerkpEvSolv1D in MODULE
control_numeric.
The eigenstates from the previous calculation are taken as
starting values. After calculation the eigenstates are backtransformed and
normalized. Finally the sparse matrices are deallocated. In 2D and 3D the
situation is analogous but there is no LAPACK any more.
See the source code and the documentation on the eigenvalues for detailed
information.
!------------------------------------------------------------------------------------------------------------------------
subroutine
calculate_kp_eigenvalues(dim_qr,num_qr,num_states,E_separate, &
chargeC,kind_kpC,eigenvalueV,spinorM)
!-----------------------------------------------------------------------------------------------------------------------
!
! calculates num_states kp-eigenstates
! dim_qr : dimension of quantum region
! num_qr ... number of quantum region
! kind_kpC = '8x8' or '6x6'
! chargeC = 'el' or 'hl'
!
! kind_kpC = '8x8' and chargeC = 'el' : the
num_states from E_separate
!
upwards are calculated
!
! kind_kpC = '8x8' and chargeC = 'hl' : the
num_states from E_separate
!
downwards are calculated
!
! kind_kpC = '6x6' and chargeC = 'hl' : the
num_states maximum states
!
are calculated
!
!--------------------------------------------------------------------
!
! uses method specified in
control_numeric: SchroedingerkpEvSolv1D:
!
! - 'lapack' diagonalizes
matrix with lapack
! - 'it_jam' uses
diagonlization routine of Jacek:
!
!
init=init_itjam
!
itere_max=itere_stagesV(stage)
!
eps_itere=epsilon_it_jam
!
itsub_max=itsub_stagesV(stage)
!
eps_itsub=epsilon_it_jam
!
iwrite=iwtite_it_jam
!
!
stage and epsilon_it_jam
!
must be defined in the upper level:
!
--> level of precision
!
!
(thus far only for 6x6 - holes)
!
! - 'arnoldi' precison =
SchroedingerkpResidual1D
!
!
! SchroedingerkpIterations1D --> number of
iterations
!
set to a default value, does not change during
!
simulation
! SchroedingerkpResidual1D --> precision for
arnoldi method
!
is defined according to precision of outer loop
!
! note: for extremous eigenvalues -->
arnoldi with option 'LR' or 'SR'
!
! for
inner eigenvalues --> arnoldi with option 'SM'
! NOTE:
all eigenvalues around E_separate are calculated
!
due to H --> H-E_separate[eV]
!
This shift is carried out here explicitly.
!
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 kp-region
!
If the calculated eigenvalues are not sufficient,
!
--> mult_num_ev_arnoldi (in the kp-region) is
!
increased by add_num_ev_arnoldi, otherwise
!
decreased by dec_num_ev_arnoldi
!
If the calculated eigenvalues are not sufficient,
!
the calculation is started anew
!
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 degeneracies
!
!------------------------------------------------------------------------------------------------------------------------
1D
In 1D they are calculated immediately after the setup of the Schrödinger matrix in
subroutine
SG_1D for the given potential phiV (INPUT).
2D/3D
The eigenvalues are calculated in subroutine
calculate_sg for the given potential
phiV (INPUT).
First the sparse matrix is defined: CALL
setup_sg).
Then the eigenvalues are calculated.
Up to now three methods are implemented:
- Conjugate gradient method of J. Majewski
--> Can only be used for extreme
eigenvalues.
- Arnoldi (ARPACK package)
- Jacobi-Davidson (default method)
At the end of this subroutine, the sparse matrices are deallocated.
subroutine calculate_sg (num_qr,num_sg,num_deg,chargeC,phiV,num_ev, &
dim_qr,eigenvaluesV,eigenfunctionsM)
!--------------------------------------------------------------------------------------------------------
! solves schroedinger equation in specified region and band
!
! input:
! num_qr ... number of quantum region
! num_sg ... number of schroedinger equation
!
(refers to those with different bandedges)
! num_deg ... for the specified class of schroedinger equations
!
(num_deg refers to different masses and boundary conditions
! chargeC = 'el' or 'hl'
! phiV (1..lengvecV) ... potential
! num_ev ... number of eigenvalues
! dim_qr ... dimension of the quantum region
!
! output:
! eigenvaluesV(1..num_ev) ... eigenvalues [in eV]
! eigenfunctionsM(1..num_ev,1..dim_qr) ... eigenfunctions
! --> in [1/sqrt(AA)]
! scaling factor is taken into account.i.e. wavefunction is scaled with this
factor
!
!------------------------------------------------------------------------
!
! note:
!
! - for energy edge on boundaries of quantum region, it holds:
! if multiple point, only the octant which
is within quantum region
! is valid
! to decide this the program only accesses
quantum_grid3DV
! but no information of guenther
! (see setup_sg for detailed information)
!----------------------------------------------------------------------------------------------------------
1)
setup_sg:
Setup of Schrödinger matrix.
Note: The Schroedinger matrix H can be
transformed to a symmetric matrix H via
H = DH ,where D is a positive definite diagonal matrix.
In setup_sg the matrix H is defined and written as a real,
symmetric sparse matrix
in MODULE
sparse_matrix_rout_real. The diagonal matrix D is stored in
diag_sgV in
MODULE
normalization_factors_sg.
(see
also
glossary)
2) The eigenvalue problem D^{1/2} H D^{1/2} is solved : eigenfunctions
F
Up to now, only ARNOLDI method is implemented:
eigensolver_real_arpack.
The eigensolver needs a subroutine, which carries out
the matrix vector product:
subroutine
matvec_sg_arpack
The precision and number of iterations are obtained
from epsilon_it_arnoldi and
max_iter_arnoldi
in MODULE
numeric_control.
3) The eigenfunctions are backtransformed via Y =
D^{1/2} F and are finally normalized.
All eigenvalues are calculated anew for given potential (phiV as input) in
subroutine
calculate_eigenstates.
The subroutine carries out following steps:
- Loop through all quantum regions.
- For each quantum region: Check if k.p is done:
get_which_kp.
If so, calculate k.p states: CALL
update_kp.
- Now treat single-band Schrödinger equation:
In each quantum region, the different Schrödinger equations can be split into
different sets, for which the same band edge is valid, e.g. the degenerate
L-valleys.
For each set, there are several Schrödinger equations, which share the same
band edge, but differ in effective masses and/or boundary conditions.
--> There are two inner loops:
do
num_sg=1,num_sg_el
do n_deg=1,num_deg
Now the specified Schrödinger equation is solved: call
calculate_sg (calculate_sg1D).
Finally, the solution is stored in the pointer sgChargeV(1)%qcV(...)%... in MODULE
quantum_solutions.
- The energy edge, as described in "8. Separation models: eigenvalues
<--> Thomas-Fermi",
is determined for single-band states. For k.p states this is done in
update_kp.
|
