k.p densities
1) 1D k.p
The densities in k.p are calculated as follows:
n(z) = sum_{kz} 1/(2 pi)^2 integral dkx
dky f(E_n) sum_{i=1..8} |Chi_{i,n}|^2 (see
script)
There are several options to carry out this
integration:
A) one-dimensional integration
When the problem is isotropic
(i.e. depends only on absolute value of k_parallel),
the integration over the 2D-BZ can be reduced to a
one-dimensional integration:
1/(2 pi)^2 integal dkx dky
f(E_n) sum_{i=1..8} |Chi_{i,n}|^2 =
= 1/(2 pi) integal dk_par k_par f(E_n) sum_{i=1..8}
|Chi_{i,n}|^2
This integration is carried
out in kp_hole_densities1D, when dim_intBC='1D' is
specified:
!------------------------------------------------------------------------------------------------------------------
integal dk_par -->
loop i=1...num_k_points , dens=dens+ (value)*space_kV(i)
k_absV = (/
(dsqrt(kpChargeV(2)%qcV(num_qr)%kx_kpV(i)**2+
&
kpChargeV(2)%qcV(num_qr)%ky_kpV(i)**2),i=1,num_k) /)
space_kV(1)=0.5d0*(k_absV(2))
space_kV(2:num_k-1)=
&
(/ (0.5d0*(k_absV(i+1)-k_absV(i-1)),i=2,num_k-1) /)
space_kV(num_k)=0.5d0*(k_absV(num_k)-k_absV(num_k-1))
! space_kV(i) = space in k-space which is occupied by k-point i
spinor_dens8V=kpChargeV(2)%qcV(num_qr)%spinorM(dim_bnd,1,i,j,iqr,1:8)
spinor_dens=sum(spinor_dens8V*dconjg(spinor_dens8V))
! calculate sum_{i=1..8} |Chi_{i,n}|^2
ev =
(kpChargeV(2)%qcV(num_qr)%eigenvaluesV(dim_bnd,1,i,j)-dphi)*abs(electron_charge)
if (ev>=E_edge) then
argV=-beta*(ev-fermi)
ergV=exponentV(argV,'fermi')
! ergV(1) = f(E_n)
select case(dimC)
case('1D')
dens=dens+spinor_dens*space_kV(i)*k_absV(i)*ergV(1)*div2pi
! div2pi=1/(2 pi)
....
! k_absV(i) = k_par
end select
end if
!-----------------------------------------------------------------------------------------------------------------
B) two-dimensionl integration
The integration over
the 2D-BZ is carried out as follows:
1/(2 pi)^2 integral
dkx dky f(E_n) sum_{i=1..8} |Chi_{i,n}|^2 =
= sum_{all different (kx,ky)} f(E_n) sum_{i=1..8} |Chi_{i,n}|^2
* weight(kx,ky)
weight(kx,ky) =
= area_kxky{which this point occupies in
k-space} * 1/(2 pi)^2 * multiplicity of that point
!----------------------------------------------------------------------------------------------------------------
spinor_dens8V=kpChargeV(2)%qcV(num_qr)%spinorM(dim_bnd,1,i,j,iqr,1:8)
spinor_dens=sum(spinor_dens8V*dconjg(spinor_dens8V))
ev =
(kpChargeV(2)%qcV(num_qr)%eigenvaluesV(dim_bnd,1,i,j)-dphi)*abs(electron_charge)
if (ev>=E_edge) then
argV=-beta*(ev-fermi)
ergV=exponentV(argV,'fermi')
select case(dimC)
.....
case('2D')
dens=dens+spinor_dens*kpChargeV(2)%qcV(num_qr)%weight_kpV(i)*ergV(1)
end select
!----------------------------------------------------------------------------------------------------------------
C) Integration via density of states
(gendos2D)
The k-points (where the Hamiltonian
is diagonalized) are provided by subroutine get_k_par1D:
subroutine
get_k_par1D(num_qr,max_point,maxval_k,num_k,kxV,kyV,weightV,dim_intBZ, &
isymmetry,num_wedges,wedgeV,nkgamx,nkgamy)
!-----------------------------------------------------------------------------------
! input:
! num_qr
... number of quantum_region
!
max_point ... maximum number of k points
!
maxval_k ... maximum absolute value of k-vector,
!
k-vectors with greater value are not taken into account
! num_k
... number of points to be returned
!
! output:
! num_k
... number of actually generated k-points
!
num_k <= max_point
!
kxV(1..max_point) .. kx - component of k-vector, only defined for
!
1..num_k [1/AA]
!
kyV(1..max_point) .. ky - component of k-vector, only defined for
!
1..num_k [1/AA]
!
weightV(1..max_point) weight factors for corresponding k-points
!
! it must be: num_k >= 1
!
kxV(1) = 0.0d0
!
kyV(1) = 0.0d0
!
!----------------------------------------------------------------------------------
!
! 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,...,k_par_max}
!
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 k_par_max
! for WZ and ZB-structure
! if dim_intBZ='2D'
! (isymmetry=0 : not
defined)
! k_points may have any
symmetry
!
! --> dim_intBZ='DO'
! we discretize this part
on a SQUARE/RECTANGLE/HEXAGONAL BZ of length k_par_max
! for WZ and ZB-structure
! 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 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 BZ is discretized
on a rectangle
! -->
num_wedges = 4
!
isymmetry= 2
!
i_kpoints per wedge = max_point/num_wedges
!
for dim_intBC='2D' and ='DO' the same k-points are provided
!-----------------------------------------------------------------------------------
USE math_constants,only:pi
implicit none
integer,intent(in) :: max_point,num_qr
character(len=2),intent(in out) :: dim_intBZ
real(8),intent(in) :: maxval_k
integer,intent(in out) :: num_k
real(8),dimension(1:max_point),intent(out) :: kxV,kyV,weightV
integer,intent(out) ::
isymmetry,num_wedges,nkgamx,nkgamy
integer,dimension(12),intent(out) :: wedgeV
integer :: max_dim,N,istatus,count,i,j
integer :: ispn
real(8),dimension(:,:),allocatable ::
xkar
real(8) :: weight,delta_k
!----------------------------------------------------------------------------------
select case(dim_intBZ)
case('1D')
! in '1D' --> weights do not play any role
! isymmetry,.... of no importance
allocate(xkar(2,max_point),STAT=istatus)
if (istatus/=0) stop 'error get_k_points:
alloc. failed'
isymmetry=0
wedgeV=0
num_wedges=0
nkgamx=0
nkgamy=0
call spk_1D(max_point,maxval_k,ispn,xkar)
do i=1,ispn
kabs=dsqrt(xkar(1,i)**2+xkar(2,i)**2)
kxV(i)=xkar(1,i)
kyV(i)=xkar(2,i)
end do
num_k=ispn
weightV=0.0d0
!
!
weights are defined for consistency reasons but of no
!
importance; are set to 1d0
deallocate(xkar,w)
case('2D','DO')
!---------------------------------------------------------------------------
! note: here the program discretizes the
whole BZ
! on a
rectangle
!---------------------------------------------------------------------------
isymmetry=2
num_wedges=4
wedgeV(1)=1
wedgeV(2)=2
wedgeV(3)=3
wedgeV(4)=4
wedgeV(5:12)=0
N=sqrt(real(int(max_point/num_wedges)))
nkgamx=N
nkgamy=N
if (N<=1) stop 'error get_k_points:
number of k-points not sufficint'
delta_k=maxval_k/(N-1)
count=1
weight=delta_k**2/(2d0*pi)**2
! NOTE: weight = area the point occupies in k-space * 1/(2 pi)^2 *
muliplicity
!----------------------------------
! first wedge
!----------------------------------
do i=1,N
do j=1,N
kxV(count)=(i-1)*delta_k
kyV(count)=(j-1)*delta_k
if
(((i==1).and.(j>1)).or.((i>1).and.(j==1)) then
weightV(count)=weight*2d0
else if
((i==1).and.(j==1)) then
weightV(count)=weight*4d0
else
weightV(count)=weight
end if
count=count+1
end do
end do
!----------------------------------
! second wedge
!----------------------------------
do i=2,N
do j=1,N
kxV(count)=-(i-1)*delta_k
kyV(count)=(j-1)*delta_k
if (j==1) then
weightV(count)=weight*2d0
else
weightV(count)=weight
end if
count=count+1
end do
end do
!----------------------------------
! third wedge
!----------------------------------
do i=1,N
do j=2,N
kxV(count)=(i-1)*delta_k
kyV(count)=-(j-1)*delta_k
if (i==1) then
weightV(count)=weight*2d0
else
weightV(count)=weight
end if
count=count+1
end do
end do
!----------------------------------
! fourth wedge
!----------------------------------
do i=2,N
do j=2,N
kxV(count)=-(i-1)*delta_k
kyV(count)=-(j-1)*delta_k
weightV(count)=weight
count=count+1
end do
end do
num_k=count-1
!------------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------
subroutine spk_1D(N,k_par,ispn,xkar)
!------------------------------------------------------------------------------
! calcultes special k-points for hexagonal lattice
in [100] direction
! here adapted for GaN
!
! input
! N
... number of special k-points in k_x 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)
!------------------------------------------------------------------------------
implicit none
integer,intent(in) :: N
integer,intent(out) :: ispn
real(8),intent(in) :: k_par
real(8),dimension(2,max_dim),intent(out)
:: xkar
!--------------------------------------------------
real(8) :: k_anf,k_end,delta_k
integer :: i
!--------------------------------------------------
if (N<=1) stop 'error : N for spk_1D must be
>1'
k_end = k_par
k_anf = 0.0d0
if (N>1) then
delta_k = (k_end-k_anf)/dble(N-1)
xkar(1,1) = 0.0d0
xkar(2,1) = 0.0d0
do i=2,N
xkar(1,i) = delta_k*dble(i-1)
xkar(2,i) = 0.0d0
end do
else
xkar(1,1) = 0.0d0
xkar(2,1) = 0.0d0
end if
!--------------------
ispn = N
end subroutine spk_1D
!------------------------------------------------------------------------------
end subroutine get_k_par1D
!-----------------------------------------------------------------------------------
Superlattice
In a superlattice one has additionally a
superlattice vector.
The miniband dispersion is also discretized and the
k_z-points are obtained via
!-----------------------------------------------------------------------------------
subroutine get_K_z_points1D(num_k,kzV,weightV)
!-----------------------------------------------------------------------------------
! num_k is input --> number of k-points
!
which are equidistant from 0 to pi/L
!
! kzV --> contains k in AA
! weightV --> such that sum (weightV) =1
!----------------------------------------------------------------------------------
USE gitter,only:n_1,h_1V
USE math_constants,only: pi
implicit none
integer,intent(in) :: num_k
real(8),dimension(num_k),intent(out) ::
kzV,weightV
!------------------------------------------
real(8) :: L,period,w
integer :: n
!---------------------------------------------------------------------------------
if (num_k<2) stop 'error get_K_z_points1D: for
superlattice num_K_z must be >1'
weightV(1)=1.0d0
weightV(num_k)=1.0d0
if (num_k>2) weightV(2:num_k-1) = 2.0d0
L = sum(h_1V(1:n_1-1)) *1d10 ! --
AA --
period = pi/L
kzV(1) = 0.0d0
do n=2,num_k
kzV(n) = period/dble(num_k-1)*dble(n-1)
end do
w=sum(weightV)
weightV=weightV/w
end subroutine get_K_z_points1D
!-------------------------------------------------------------------------------------
In density calculation this means an additional
integration over K_z:
do n=1,num_K_z
spinor_dens8V=kpChargeV(2)%qcV(num_qr)%spinorM(dim_bnd,n,i,j,iqr,1:8)
spinor_dens=sum(spinor_dens8V*dconjg(spinor_dens8V))
ev =
(kpChargeV(2)%qcV(num_qr)%eigenvaluesV(dim_bnd,n,i,j)-dphi)*abs(electron_charge)
if (ev>=E_edge) then
argV=-beta*(ev-fermi)
ergV=exponentV(argV,'fermi')
weight =
kpChargeV(2)%qcV(num_qr)%weight_kzV(n) ! weight
for superlattice,
select case(dimC)
case('1D')
dens=dens+spinor_dens*space_kV(i)*k_absV(i)*ergV(1)*div2pi*weight
case('2D')
dens=dens+spinor_dens*kpChargeV(2)%qcV(num_qr)%weight_kpV(i)*ergV(1)*weight
end select
end if
end do
| 
