Brillouin zone integration
ln MODULE
quantum_solutions variables for both Schrödinger (1-band) and k.p (multi-band) solutions
are provided.
TYPE sg_solution:
num_K_z -->
Actual number of kz values
max_num_K_Z --> Maximum
number of kz values for superlattices
otherwise = 1 (for kz=0)
Note: Superlattice not valid any more :
num_K_z = max_num_K_Z = 1
number_of_ev --> Number of
eigenvalues calculated
max_number_of_ev --> Maximum number of eigenvalues
num_deg -->
Number of different sg_equations for that energy level.
Must be consistent with def_quantum_modelsM.
eigenvaluesV (1..num_deg,1..max_num_K_Z,1..number_of_ev)
--> eigenvalues [eV]
eigenfunctionsM(1..num_deg,1..max_num_K_Z,1..number_of_ev,1..dim_qr)
--> eigenfunctions in [1/sqrt(AA)]
K_zV(1..num_K_z) --> Special kz for
superlattices (1/AA), otherwise not allocated.
weight_kzV(1..num_K_z)
--> Weights for kz for superlattices, otherwise not allocated.
energy_edge_input -> Separation quantum mechanical
<--> classical electrons
relative to lowest point of conduction band
in [eV]
energy_edge --> See above,
here the value that is actually used. [eV]
SeparationModelC ...
.. 'eigenvalue' --> All eigenvalues up to
max_number_ev are considered.
Above topmost eigenvalue electrons are considered to be classical.
.. 'max_energy' --> energy_edge is valid.
It is read from input and does not change during calculation.
.. 'edge_model' --> energy_edge is valid.
It is permanently updated according to the separation model.
See also subroutine
determine_edge_qm for further explanation of separation models
and the glossary and
link therein.
TYPE :: sg_diff_energy
num_sg_diff_en --> = num_sg_el
TYPE (sg_solution),DIMENSION(:),POINTER :: data_deg_sgM
NOTE: num_sg_el = num_sg_diff_en contains the number of
Schrödinger equations which have different energy edge. In data_deg_sgM
one must additionally distinguish between Schrödinger equations which have the
same energy but different effective masses (=num_deg). See
documentation of M. Sabathil for detailed information.
Units
Energies in eV
EigenfunctionsM in 1/sqrt(Angstrom)
K_zV in 1/Angstrom
TYPE kp_solution:
input_num_kp_par --> number of k|| for whole
2D-BZ (via input)
max_num_kp_par --> maximum number of k||
max_num_K_Z --> maximum
number of kz values for superlattices
otherwise =1 (for kz=0)
max_number_of_ev --> maximum number of eigenvalues
num_kp_par --> actual
number of k||
num_K_z -->
actual number of kz values
number_of_ev --> number of
eigenvalues calculated
dim_boundary --> =1, if Dirichlet or
Neumann or periodic
=2, if mixed Dirichlet/Neumann
eigenvaluesV(1..dim_boundary,1..max_num_K_Z,1..max_num_kp_par,1..number_of_ev)
--> eigenvalues [eV]
spinorM
(1..dim_boundary,1..max_num_K_Z,1..max_num_kp_par,1..number_of_ev,
1..dim_qr,1:6/8)
--> spinor [1/sqrt(AA)]
kx_kpV,ky_kpV(1..num_kp_par) --> special k|| points
(1/AA)
kz_kpV (1..num_K_z)
--> special kz for superlattices (1/AA)
otherwise not allocated
weight_kpV (1..num_kp_par) --> weights of
k-points
weight_kzV (1..num_K_z) -->
weights for kz for superlattices,
otherwise not allocated
dim_intBZ --> '1D'
or '2D' : for integration routine over BZ
'1D' --> one-dimensional integration
'2D' --> two-dimensional integration
'DO' --> first calculate density of states, then integrate over energy
ispecify_case --> See
subroutine get_case_integrate_BZ in
MODULE mod_get_dim_2D_BZ for detailed information.
isymmetry, num_wedges, wedgesV(1..12), nkgamx, nkgamy
...
only for dim_intBZ='DO':
These are the data provided by subroutine
get_k_par1D
See description there for detailed information.
k_par_max --> for
integration with gendos [ 1/AA ]
energy_edge_input -> Separation
quantummechanical <-> classical electrons
relative to lowest point of conduction band
in [eV].
energy_edge --> See above,
here the value that is actually used [eV].
SeparationModelC
.. 'eigenvalue' --> All eigenvalues up to
max_number_ev are considered.
Above topmost eigenvalue, electrons are considered to be classical.
.. 'max_energy' --> energy_edge is valid.
It is read from input and does not change during calculation.
energy_edge is relative to extreme bandedge.
All eigenstates up to this value are considered in quantum mechanical
calculation.
Above this value carriers are considered to be classical.
.. 'edge_model' --> energy_edge is valid.
It is permanently updated according to the separation model.
See also subroutine
determine_edge_qm for further explanation of separation models
and the glossary and
link therein.
boundaryC:
= 'per' -->
Periodic boundary conditions (superlattice)
= 'neu' -->
Neumann boundary conditions
= 'dir' -->
Dirichlet boundary conditions
= 'mix' -->
Mixed Dirichlet/Neumann conditions are taken.
readjust_kp_par --> In order to
determine the range for k|| in 2D BZ, the bulk Hamiltonian is
diagonalized in the two directions perpendicular to the quantization
direction:
After the discretization the program checks, if this range is adequate
and sets
readjust_kp_par for the next iteration:
Shift the range about readjust_kp_par% (percent)
in the next loop.
Not used thus far.
mult_num_ev_arnoldi --> see
numeric_control: arnoldi_kp
NOTE: Up to now only Dirichlet values are implemented.
kind_kpC = '8x8' , '6x6'
kpChargeV(1)%qcV,kpChargeV(2)%qcV(1..num_qr) of
TYPE(kp_solution)
Units
Energies in eV
spinorM in 1/sqrt(AA)
kx_kpV, ky_kpV in 1/AA
kparallel = k||
k||² = kx² + ky²
For k|| integration the whole 2D Brillouin zone (BZ) is
discretized on a k grid of max_num_kp_par points.
input_num_kp_par is the number of k|| specified via
input in
$quantum-model-electrons/$quantum-model-holes
(num-kp-parallel = ...).
--> It always refers to the total number of k||
points in the whole 2D Brillouin zone.
Subroutine get_dim_2D_BZ(num_qr,input_num_kp_par,dim_intBZ) provides for num_qr(quantum region) and input_num_kp_par
and dim_intBZ(via input)
1) dim_intBZ
2) max_num_kp_par (<= input_num_kp_par)
Example:
On input it is checked if the problem is isotropic (the case for wurtzite
growth direction [0001] (= quantization direction)).
If so: dim_intBZ --> '1D' and
max_num_kp_par --> SQRT(input_num_kp_par)/2
If for a 2D integration the BZ is reduced to an irreducible part, the
number of k|| points is renormalized accordingly:
So if only 1/8th of the BZ is used (square, irreducible wedge)
--> max_num_kp_par = input_num_kp_par/8
max_num_kp_par is always the dimension of the arrays
kx_kpV, ...
-
allocate_quantum_states calls subroutine
get_dim_2D_BZ
(when
first time called: flag calculate_kL = .FALSE.)
CALL get_dim_2D_BZ
(num_qr,'el',input_num_kp_par,max_num_kp_par,int_method,ispecify_case,.FALSE.)
num_qr -->
Number of quantum cluster
input_num_kp_par -->
Number of k|| for whole 2D BZ
(unchanged, from input)
chargeC
= 'el' or 'hl'
Optional variables are not present.
(num_kp_par,isymmetry,num_wedges,wedgeV,
nkgamx,nkgamy,kxV,kyV,weightV)
with input input_num_kp_par which is the number of k||
points over the whole Brillouin zone that should be integrated. Then this
routines determines from this information how many k|| points really
have to be calculated (i.e. it checks if the problem can be simplified
since special orientations need less k|| points to be considered)
and returns the actual number (which can be the same, or e.g. 1/8th) that
should be calculated.
dim_intBZ ... Input: The option from Input Parser,
may be changed here. The value on output is written to MODULE
quantum_solutions in SUBROUTINE
allocate_quantum_states.
max_num_kp_par
ispecify_case ... Output: Stores the valid value.
- When this subroutine is called the second time (flag
calculate_kL =
.TRUE.), then the optional variables ARE present.
dim_intBZ ... only input:
It is the value that has been stored in MODULE
quantum_solutions.
max_num_kp_par ... is input: Contains the
maximum number of k|| in order to allocate the arrays.
ispecify_case ... only: input
This SUBROUTINE discretizes a BZ which extends from k|| = [0d0 ...
1d0] [1/AA].
NOTE: The k points and the weights must be scaled accordingly.
num_kp_par --> Number of k|| that are actually
used for discretization.
If one has wurtzite structure and growth direction [0001]=quantization
direction:
--> This is THE ONLY CASE when one-dimensional integration is
possible.
So if dim_intBZ = '1D'
--> k points along line ky=0, kx in{0,...,1.0d0} are
taken.
(isymmetry = 0: not defined)
The eigenstates are calculated for density calculation.
--> Only a small part of the BZ is needed.
--> So we discretize this part on a SQUARE of length
[1/Angstrom] for wurtzite and zincblende structures.
if dim_intBZ = '2D'
(isymmetry = 0: not defined)
k points may have any symmetry.
--> dim_intBZ = 'DO'
We discretize this part on a square/rectangular/hexagonal
BZ of length [1/Angstrom] for wurtzite and
zincblende structures.
Then, in isymmetry additional
information is provided on the structure of the discretized region:
isymmetry
1 = square
2 = rectangular
3 = hexagonal BZ (according to
gendos)
nkgamx, nkgamy --> for
gen2Ddos
num_wedges
= 1 , wedgeV(1..num_wedges)
\ 3 |
2 /
\ |
/
\ |
/
4 \ | /
1
---------+---------
5 / | \
8
/ |
\
/ |
\
/ 6 | 7
\
Example:
wedge 1 as special k points -->
num_wedges=1
wedgeV(1)=1, wedgeV(2:8)
undefined
wedge 1 6 7 8
--> num_wedges=4
wedgeV(1 .. 4) =
(1,6,7,8)
Note: For each point a weight
factor is provided.
If the whole BZ is discretized -->
weight=1
Note: The wedges must be in
ascending order. If wedge i has been written, the next wedge only
contains k points not contained in wedge
i. The last (specified) wedge only contains k points not
contained in the first (specified) wedge.
isymmetry = 2 <-->
rectangular --> one has only 4 wedges
= 3 <--> hexagonal --> one has 12 wedges
The points must be provided in the same
order as is required for gendos2D.
Note: In this version the whole Brillouin zone is discretized on a
rectangle.
--> num_wedges = 4
isymmetry =
2
i_kpoints per wedge =
max_num_kp_par/num_wedges
For dim_intBC='2D'
and ='DO' the same k points are provided.
Note: isymmetry, num_wedges, wedgeV,
nkgamx, nkgamy only play a role for dim_intBZ='DO',
weightV only for dim_intBZ='2D'.
- We
CALL get_dim_2D_BZ for the first time (flag calculate_kL = .FALSE.):
Now we check the case to determine which integration to use, i.e. if
it is wurtzite or zincblende in a certain growth/quantization direction (namely
[0001] or [001]). This will be done by calling the subroutine
get_case_integrate_BZ.
Input: num_qr, dim_intBZ,
Output: returns for 2D integration ispecify_case = ...
dim_intBZ='1D' :
crystC='wz' and
(rotM = I or rotM=-I --> [0001]) : ispecify_case = 1
crystC='zb'
and (rotM = I or rotM=-I --> [001] ) : ispecify_case =
2
all other cases
: ispecify_case = 3
dim_intBZ='2D' :
crystC='wz' and
(rotM = I or rotM=-I --> [0001]) : ispecify_case = 4
crystC='zb'
and (rotM = I or rotM=-I --> [001] ) : ispecify_case =
5
all other cases
: ispecify_case = 6
dim_intBZ='DO' :
crystC='wz' and
(rotM = I or rotM=-I --> [0001]) : ispecify_case = 7
crystC='zb'
and (rotM = I or rotM=-I --> [001] ) : ispecify_case =
8
all other cases
: ispecify_case = 9
SUBROUTINE get_case_integrate_BZ calls SUBROUTINE
get_crystal_and_rotation
in order to obtain the crystal structure (wurtzite or zincblende) and the
rotation matrix that transforms the crystal coordinate system into the simulation system. With the
help of this matrix, one can check for these cases:
- case wurtzite: Is discretization
(i.e. quantization) direction along [0001] or [000-1]?
- case zincblende: Is discretization (i.e. quantization) direction along [001] or [00-1]?
(Of course, the other two directions (parallel space axes) must be checked
as well
in order to avoid cases like this one:
x = [ 0 0 1] - quantization direction
y = [ 2 3 0]
z = [-3 2 0]
- all other cases
Input for SUBROUTINE
get_crystal_and_rotation
is num_qr = number of quantum cluster.
Output is crystC = 'wz' or 'zb'
and rotM = (3x3) rotation matrix.
For setup: Now we redefine dim_intBZ and determine
max_num_kp_par according to the returned values of ispecify_case.
Case 1:
Note: One dimensional integration (wurtzite -> [0001])
Integration over Brillouin zone: One-dimensional
--> max_num_kp_par = SQRT(input_num_kp_par)/2
dim_intBZ = '1D'
Case 2, 5:
Note: 2D BZ is discretized only on irreducible wedge.
--> max_num_kp_par = input_num_kp_par/8
dim_intBZ = '2D'
Case 3, 4, 6: (Case 4 should get an extra treatment in the near future ->
gendos "1/9th")
Note: Whole BZ is discretized.
--> max_num_kp_par = input_num_kp_par
dim_intBZ = '2D'
Case 8:
Note: 2D BZ is discretized only on irreducible wedge.
--> max_num_kp_par = input_num_kp_par/8
dim_intBZ = 'DO'
Case 7, 9: (Case 7 should get an extra treatment in the near future ->
gendos "1/9th")
Note: Whole BZ is discretized.
--> max_num_kp_par = input_num_kp_par
dim_intBZ = 'DO'
- Now we
CALL get_dim_2D_BZ for the second time
(flag calculate_kL = .TRUE.):
For calculation of k||: dim_intBZ as input.
Note: Range of k|| [0,...,1.0] [1/AA]
Depending on ispecify we are calling one of the following
subroutines:
-
ispecify = 1
--> CALL discretize_1D - 1D
integration
which itself calls subroutine
spk_1D(max_num_kp_par,ispn,xkar)
(see description further below)
only possible for wurtzite materials
in [0001] growth direction
-
ispecify = 3, 4, 6, 7, 9
--> CALL discretize_whole_2D - 2D integration whole Brillouin zone
4 and 7 should get a special
subroutine in the future ("1/9th" gendos)
for rectangular lattice, 4 wedges
-
ispecify = 2, 5, 8
--> CALL discretize_2D_irr8 - 2D integration BZ, only
irreducible part (1/8th),
for square lattice, 1 wedge
Note: Here the program discretizes the irreducible
wedge of a cubic BZ (on a square).
spk_1D(max_num_kp_par,ispn,xkar)
Calculates special k points for hexagonal lattice in [100] direction. Here
adapted for GaN.
Input: N
Number of special k points in kx direction.
max_dim is defined for consistency reasons but of no importance.
max_dim points are generated.
Output: ispn ... number of k points (=N)
xkar ... (2,N) .. (kx,ky)
Note: Range of k|| from 0 to 1.0
MODULE mod_get_dim_2D_BZ
contains
SUBROUTINE get_dim_2D_BZ
SUBROUTINE discretize_1D
SUBROUTINE discretize_whole_2D
SUBROUTINE discretize_2D_irr8
SUBROUTINE spk_1D(N,ispn,xkar)
SUBROUTINE
get_crystal_and_rotation
SUBROUTINE get_case_integrate_BZ
END MODULE
SUBROUTINE gen2ddos
| 
