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