|||
module global
implicit none
! all unit used in the program is atomic unit
! energy : Hartree
! length : Bohr
real(kind=8),parameter::Ry=0.5d0,pi=3.14159265358979,Ang=1.8903592
integer,parameter::Num_elem=14,Near_neighbor=10,N_bands=12,N_point=40,NMax_pw=200
integer,parameter::vect_norm(Near_neighbor)=(/ 0, 3, 4, 8, 11, 12, 16, 19, 20, 24/)
character(len=4),parameter::semi(Num_elem)=(/ 'Si','Ge','Sn','GaP','GaAs','AlSb','InP',&
&'GaSb','InAs','InSb','ZnS','ZnSe','ZnTe','CdTe' /)
! a0 lattice constant of fcc structure semiconductor
real(kind=8),parameter::a0_semi(Num_elem)=(/5.43d0,5.66d0,6.49d0,5.44d0,5.64d0,6.13d0,5.86d0,&
& 6.12d0,6.04d0,6.48d0,5.41d0,5.65d0,6.07d0,6.41d0/),&
str_fac(6,Num_elem)=(/-0.21d0, 0.04d0, 0.08d0, 0.00d0, 0.00d0, 0.00d0,&
-0.23d0, 0.01d0, 0.06d0, 0.00d0, 0.00d0, 0.00d0,&
-0.20d0, 0.00d0, 0.04d0, 0.00d0, 0.00d0, 0.00d0,&
-0.22d0, 0.03d0, 0.07d0, 0.12d0, 0.07d0, 0.02d0,&
-0.23d0, 0.01d0, 0.06d0, 0.07d0, 0.05d0, 0.01d0,&
-0.21d0, 0.02d0, 0.06d0, 0.06d0, 0.04d0, 0.02d0,&
-0.23d0, 0.01d0, 0.06d0, 0.07d0, 0.05d0, 0.01d0,&
-0.22d0, 0.00d0, 0.05d0, 0.06d0, 0.05d0, 0.01d0,&
-0.22d0, 0.00d0, 0.05d0, 0.08d0, 0.05d0, 0.03d0,&
-0.20d0, 0.00d0, 0.04d0, 0.06d0, 0.05d0, 0.01d0,&
-0.22d0, 0.03d0, 0.07d0, 0.24d0, 0.14d0, 0.04d0,&
-0.23d0, 0.01d0, 0.06d0, 0.18d0, 0.12d0, 0.03d0,&
-0.22d0, 0.00d0, 0.05d0, 0.13d0, 0.10d0, 0.01d0,&
-0.20d0, 0.00d0, 0.04d0, 0.15d0, 0.09d0, 0.04d0 /)
end module
program main
! This program use the empirical method to sovle the Schrodinger Equation ,and use the result to plot the
! band structure of some semicondutor ,inlcude 'Si','Ge','Sn','GaP','GaAs','AlSb','InP', 'GaSb','InAs',
! 'InSb','ZnS','ZnSe','ZnTe','CdTe' . In orde to compile this program ,you have to install MKL on your systerm .
use global
implicit none
integer::i,choose,s,icount,ia,ic,ig2,igx,igx2,jgy,&
jgy2,kgz,kgz2,igxy2,igxyz2,ik_curr,ik,aux,&
G_1D
real(kind=8)::a,sumx,X0,Y0,Z0,dkx,dky,dkz,sumx0,sumx1,sumx2
real(kind=8),intrinsic::sqrt,floor
integer,allocatable::G_vect(:,:),G_vect1(:,:)
real(kind=8)::k_vec(200)
write(*,*)"This program use Empirical Potential Method to calculate the bandstructure"
write(*,*)" of semicontucors ,please choose one of the following semiconductors :"
do i=1,Num_elem
write(*,"(i2,' ==> ',a4)")i,semi(i)
end do
write(*,*)"choose a number ===>"
read(*,*)choose
! potentials in terms of symmetric (VS(G)=V1+V2)) and antisymmetric (VA(G)=V1-V2)) dom factors
! V0(G) is given by V0(G)=cos(G*t)VS(G)+i*sin(G*t)VA(G)
allocate(G_vect(NMax_pw,3))
icount=0
do i=1,Near_neighbor
ig2=vect_norm(i)
a=sqrt(ig2*1.00d0)
ia=floor(a)
ic=0
do igx=-ia,ia
igx2=igx*igx
do jgy=-ia,ia
jgy2=jgy*jgy
igxy2=igx2+jgy2
do kgz=-ia,ia
kgz2=kgz*kgz
igxyz2=igxy2+kgz2
if (igxyz2==ig2) then
ic=ic+1
G_vect(ic+icount,1)=igx
G_vect(ic+icount,2)=jgy
G_vect(ic+icount,3)=kgz
end if
end do
end do
end do
!call gen_gvecs(i,ia,ig2,ic)
icount=icount+ic
write(*,"(i21x,i21x,i21x)"),ig2,ic,icount
end do
allocate(G_vect1(icount,3))
Write(*,*)"All G vectors have been used "
do i=1,icount
G_vect1(i,1)=G_vect(i,1)
G_vect1(i,2)=G_vect(i,2)
G_vect1(i,3)=G_vect(i,3)
write(*,"((3i3))")G_vect1(i,1),G_vect1(i,2),G_vect1(i,3)
end do
deallocate(G_vect)
G_1D=icount
ik_curr=0
sumx=0.0d0
X0=0.5d0
Y0=0.5d0
Z0=0.5d0
! print*
print*,'===================================================================='
print*,'L -> Gamma'
! print*
do ik=1,N_point+1
aux=ik+ik_curr
print*,'ik=',ik
dkx=0.50d0*(1.0d0-(ik-1.0d0)/real(N_point,kind=8))
dky=0.50d0*(1.0d0-(ik-1.0d0)/real(N_point,kind=8))
dkz=0.50d0*(1.0d0-(ik-1.0d0)/real(N_point,kind=8))
! print*,dkx,dky,dkz
! print*,"entry" used for debug
! pause
call Hamilt(G_vect1,choose,G_1D,aux,dkx,dky,dkz)
sumx=sqrt(((X0-dkx)*(X0-dkx)+(Y0-dky)*(Y0-dky)+(Z0-dkz)*(Z0-dkz))*1.0d0)
k_vec(aux)=sumx
end do
ik_curr=ik
! ! Gamma -> X
sumx0=sumx
X0=0.0d0
Y0=0.0d0
Z0=0.0d0
! print*
print*,'===================================================================='
print*,'Gamma -> X'
! print*
do ik=1,N_point
!print*,ik
aux=ik+ik_curr
dkx=real(ik)/real(N_point,kind=8)
dky=0.0d0
dkz=0.0d0
call Hamilt(G_vect1,choose,G_1D,aux,dkx,dky,dkz)
sumx=sqrt(((X0-dkx)*(X0-dkx)+(Y0-dky)*(Y0-dky)+(Z0-dkz)*(Z0-dkz))*1.0d0)
k_vec(aux)=sumx+sumx0
end do
ik_curr=ik+ik_curr
! ! X -> U
sumx1=sumx
X0=1.0d0
Y0=0.0d0
Z0=0.0d0
! print*
print*,'===================================================================='
print*,'X -> U'
! print*
do ik=1,N_point
!print*,ik
aux=ik+ik_curr
dkx=1.0d0
dky=(ik)/real(N_point,kind=8)*1.0d0/4.0d0
dkz=(ik)/real(N_point,kind=8)*1.0d0/4.0d0
call Hamilt(G_vect1,choose,G_1D,aux,dkx,dky,dkz)
sumx=sqrt(((X0-dkx)*(X0-dkx)+(Y0-dky)*(Y0-dky)+(Z0-dkz)*(Z0-dkz))*1.0d0)
k_vec(aux)=sumx+sumx1+sumx0
end do
ik_curr=ik+ik_curr
! ! K -> Gamma
sumx2=sumx
X0=0.75d0
Y0=0.75d0
Z0=0.0d0
! print*
print*,'===================================================================='
print*,'K -> Gamma'
! print*
do ik=1,N_point
!print*,ik
aux=ik+ik_curr
dkx=0.75d0*(1.0d0-(ik)/real(N_point,kind=8))
dky=0.75d0*(1.0d0-(ik)/real(N_point,kind=8))
dkz=0.0d0
call Hamilt(G_vect1,choose,G_1D,aux,dkx,dky,dkz)
sumx=sqrt((X0-dkx)*(X0-dkx)+(Y0-dky)*(Y0-dky)+(Z0-dkz)*(Z0-dkz))
k_vec(aux)=sumx+sumx1+sumx0
end do
end program
subroutine Hamilt(G_vect,choose,G_1D,aux,dkx,dky,dkz)
use global
implicit none
integer::i,j,G_1D,idgx,idgy,idgz,id2,choose,aux
real(kind=8)::dkx,dky,dkz,V_s,V_a,gx_i,gy_i,gz_i,gx_j,gy_j,gz_j,&
gkx2,gky2,gkz2,Hreal,Himag,tau_fac,ct,st,&
Vs3,Vs8,Vs11,Va3,Va8,Va11,dgx,dgy,dgz,a0
integer::G_vect(G_1D,*)
real(kind=8)::E(N_bands),E_sort(G_1D)
complex(kind=8)::H_EPM(G_1D,G_1D)
!mkl variable
!======================================================
character(len=1),parameter::jobvl='N',jobvr='N'
integer::LDA,LDVL,LDVR,INFO,LWORK
integer(kind=8),parameter::LWMAX=100000
complex(kind=8),allocatable::VR(:,:),W(:),RWORK(:),VL(:,:)
complex(kind=8)::WORK(LWMAX)
!=====================================================
intrinsic::cmplx,floor,anint,sin,cos,int,min
external::ZGEEV
Vs3=str_fac(1,choose)*Ry
Vs8=str_fac(2,choose)*Ry
Vs11=str_fac(3,choose)*Ry
Va3=str_fac(4,choose)*Ry
Va8=str_fac(5,choose)*Ry
Va11=str_fac(6,choose)*Ry
a0=a0_semi(choose)*Ang
H_EPM(G_1D,G_1D)=cmplx(0.0d0,0.0d0)
do i=1,G_1D
gx_i=real(G_vect(i,1),kind=8)
gy_i=real(G_vect(i,2),kind=8)
gz_i=real(G_vect(i,3),kind=8)
do j=1,G_1D
gx_j=real(G_vect(j,1),kind=8)
gy_j=real(G_vect(j,2),kind=8)
gz_j=real(G_vect(j,3),kind=8)
if (i==j) then
gkx2=(gx_i+dkx)*(gx_i+dkx)
gky2=(gy_i+dky)*(gy_i+dky)
gkz2=(gz_i+dkz)*(gz_i+dkz)
Hreal= (gkx2+gky2+gkz2)/2.0d0*(pi/a0)**2.0d0
H_EPM(i,j)=cmplx(Hreal,0.0d0)
else
dgx=gx_i-gx_j
dgy=gy_i-gy_j
dgz=gz_i-gz_j
idgx=anint(dgx)
idgy=anint(dgy)
idgz=anint(dgz)
id2=idgx*idgx+idgy*idgy+idgz*idgz
tau_fac=2.0d0*pi/8.0d0
if ((id2.eq.3)) then
V_s=Vs3
V_a=Va3
elseif ((id2.eq.8)) then
V_s=Vs8
V_a=Va8
elseif ((id2.eq.11)) then
V_s=Vs11
V_a=Va11
else
V_s=0.0d0
V_a=0.0d0
endif
ct=cos((dgx+dgy+dgz)*tau_fac)
st=sin((dgx+dgy+dgz)*tau_fac)
Hreal= V_s*ct
Himag= V_a*st
H_EPM(i,j)=cmplx(Hreal,Himag)
end if
end do
end do
! call the MKL subroutine
LDA=G_1D
LDVL=G_1D
LDVR=G_1D
allocate(RWORK(2*G_1D))
allocate(VL(LDVL,G_1D))
allocate(VR(LDVR,G_1D))
allocate(W(G_1D))
!query the optimal worksapce
LWORK=-1
call ZGEEV('Vectors','Vectors',G_1D,H_EPM,LDA,W,VL,LDVL,VR,LDVR,WORK,LWORK,RWORK,INFO)
LWORK=min(LWMAX,int(WORK(1)))
!solve eigenproblem
call ZGEEV('Vectors','Vectors',G_1D,H_EPM,LDA,W,VL,LDVL,VR,LDVR,WORK,LWORK,RWORK,INFO)
!check for convergence
If(INFO.GT.0)then
write(*,*)'The algorithm failed to compute eigenvalues.'
stop
end if
E_sort=1.0d0
call sort(W,E_sort,G_1D)
E=E_sort(1:N_bands)
write(*,"('print the K-point coordinate: (',f7.4,',',f7.4,',',f7.4,')' )")dkx,dky,dkz
do i=1,N_bands
!ek(aux,1:N_bands)=E(N_bands)
write(*,"(d15.8)")E(i)
end do
!print*,"the next K-point!"
return
end subroutine
!Auxiliary routine: printing a matrix.
subroutine PRINT_MATRIX( DESC, M, N, A, LDA )
character*(*) DESC
integer M, N, LDA
complex A( LDA, * )
integer I, J
write(*,*)
write(*,*) DESC
do I = 1, M
write(*,200) ( A( I, J ), J = 1, N )
end do
200 format( 5(:,1X,'(',d15.8,',',d15.8,')') )
return
end subroutine
subroutine sort(W,E_sort,G_1D)
implicit none
complex(kind=8)::W(*)
real(kind=8):: temp,E_sort(*)
integer :: i,j,G_1D
E_sort(1:G_1D)=real(W(1:G_1D))
do i=1,G_1D-1
do j=i+1,G_1D
if (E_sort(i) .gt. E_sort(j)) then
temp = E_sort(i)
E_sort(i) =E_sort(j)
E_sort(j) = temp
endif
enddo
enddo
return
end subroutine
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-6-4 06:56
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社