!
!----------------------------------------------------------------
program epm
!----------------------------------------------------------------
!
! test-case for the GW-DFPT proposal
! empirical pseudopotential method implementation
! for nonpolar tetrahedral semiconductors
!
! Fri May 30, 4.26 PM:
! great! it gives exactly the same Si bandstructure as in
! Cohen&Bergstresser!
!
! NOTE: I use only one G-grid. The reason is that, while V(G-G')
! has in principle components 2G_max, in the EPS scheme the
! largest value happens for G2=11, therefore we can simply use
! the smooth grid...
! To perform a true scaling test we need to include the double
! grid (and count the operations on the smooth and on the dense
! grids).
!
! NOTE2: the k-dependent Hamiltonian is real and symmetric
! in reciprocal space
! (but not in real space because of the momentum component)
! therefore we can find eigenstates with c(G) all real
! (the eigenvectors of a rs matrix form a orthogonal matrix)
! of course psi(r) does not need to be real
!
! Note that the Hamiltonian in G-space is real for the
! diamond structure
!
!----------------------------------------------------------------
!
use gspace
use parameters
use constants
implicit none
!
! variables
!
real(dbl) :: gcutm, arg
real(dbl) :: xk (3,nk), deltag(3), kplusg(3), xq(3,nq)
integer :: ig, ik, ig1, ig2, ideltag, notfound, ierr, ibnd, iq
logical :: allowed, equiv
real(dbl), allocatable :: t(:,:), ss(:), vs(:), h(:,:,:)
complex(dbl), allocatable :: vr(:)
real(kind=DP), parameter :: eps8 = 1.0D-8
real(kind=DP), allocatable :: eval(:), u(:,:), fv1(:), fv2(:), hk(:,:)
real(kind=DP), allocatable :: psi(:,:), e(:)
!
! temporary variables - debug
integer :: i, j, k, ir, ir1, ir2
real(kind=DP), allocatable :: psincr(:), hpsincr(:)
complex(kind=DP), allocatable :: psincc(:)
!
!
!----------------------------------------------------------------
! DEFINE THE CRYSTAL LATTICE
!----------------------------------------------------------------
!
! direct lattice vectors (cart. coord. in units of a_0)
! [Yu and Cardona, pag 23]
!
at( :, 1) = (/ 0.0, 0.5, 0.5 /) ! a1
at( :, 2) = (/ 0.5, 0.0, 0.5 /) ! a2
at( :, 3) = (/ 0.5, 0.5, 0.0 /) ! a3
!
! reciprocal lattice vectors (cart. coord. in units 2 pi/a_0)
!
bg( :, 1) = (/ -1.0, 1.0, 1.0 /) ! b1
bg( :, 2) = (/ 1.0, -1.0, 1.0 /) ! b2
bg( :, 3) = (/ 1.0, 1.0, -1.0 /) ! b3
!
! atomic coordinates (cart. coord. in units of a_0)
!
tau( :, 1) = (/ 0.125, 0.125, 0.125 /)
tau( :, 2) = (/ -0.125, -0.125, -0.125 /)
!
!----------------------------------------------------------------
! GENERATE THE G-VECTORS
!----------------------------------------------------------------
!
! the G^2 cutoff in units of 2pi/a_0
! Note that in Ry units the kinetic energy is G^2, not G^2/2
! (note for the Hamiltonian we need to double the size, 2Gmax, hence the factor 4)
!
gcutm = four * ecutwfc / tpiba2
! gcutm = ecutwfc / tpiba2
!
! set the fft grid
!
! estimate nr1 and check if it is an allowed value for FFT
!
nr1 = 1 + int (2 * sqrt (gcutm) * sqrt( at(1,1)**2 + at(2,1)**2 + at(3,1)**2 ) )
nr2 = 1 + int (2 * sqrt (gcutm) * sqrt( at(1,2)**2 + at(2,2)**2 + at(3,2)**2 ) )
nr3 = 1 + int (2 * sqrt (gcutm) * sqrt( at(1,3)**2 + at(2,3)**2 + at(3,3)**2 ) )
!
do while (.not.allowed(nr1))
nr1 = nr1 + 1
enddo
do while (.not.allowed(nr2))
nr2 = nr2 + 1
enddo
do while (.not.allowed(nr3))
nr3 = nr3 + 1
enddo
!
call ggen( gcutm)
!@@
goto 10
!@@
!
!----------------------------------------------------------------
! GENERATE THE k-POINTS
!----------------------------------------------------------------
!
! brutally defined by hand - here I should read from file
!
do ik = 1, nk
read (5,*) xk(1,ik), xk(2,ik), xk(3,ik)
enddo
!
!----------------------------------------------------------------
! CONSTRUCT THE HAMILTONIAN
!----------------------------------------------------------------
!
allocate ( t(ngm, nk), ss(ngm), vs(ngm) )
!
! the structure factor
!
do ig = 1, ngm
arg = twopi * ( g(1,ig) * tau( 1, 1) + g(2,ig) * tau( 2, 1) + g(3,ig) * tau( 3, 1) )
ss (ig) = cos ( arg )
enddo
!
! the empirical pseudopotential
!
vs = zero
! integer comparison - careful with other structures
do ig = 1, ngm
if ( int ( gl(igtongl(ig)) ) .eq. 3 ) then
vs (ig) = v3
elseif ( int ( gl(igtongl(ig)) ) .eq. 8 ) then
vs (ig) = v8
elseif ( int ( gl(igtongl(ig)) ) .eq. 11 ) then
vs (ig) = v11
endif
enddo
!
! the k-dependent kinetic energy in Ry
! [Eq. (14) of Ihm,Zunger,Cohen J Phys C 12, 4409 (1979)]
!
do ik = 1, nk
do ig = 1, ngm
kplusg = xk(:,ik) + g(:,ig)
t ( ig, ik) = tpiba2 * dot_product ( kplusg, kplusg )
enddo
enddo
!
allocate ( h (ngm, ngm, nk) )
do ig = 1, ngm
do ik = 1, nk
h ( ig, ig, ik) = t ( ig, ik)
enddo
enddo
!
do ig1 = 1, ngm
do ig2 = ig1+1, ngm
!
! define ideltag
!
deltag = g(:,ig1) - g(:,ig2)
ideltag = 1
do while ( .not. equiv ( deltag, g(:,ideltag) ) .and. ideltag.lt.ngm )
ideltag = ideltag + 1
enddo
!
! the free-electron dispersions look ok (vs=0)
! when compared to Yu/Cardona Fig. 2.8
!
do ik = 1, nk
if (ideltag.ne.ngm) &
h ( ig1, ig2, ik) = ss (ideltag) * vs (ideltag)
h ( ig2, ig1, ik) = h ( ig1, ig2, ik)
enddo
!
enddo
enddo
!
! the empirical pseudopotential in real space
! for further use in h_psi
!
allocate ( vr(nr) )
vr = czero
do ig = 1, ngm
vr ( nl ( ig ) ) = dcmplx ( ss (ig) * vs (ig), zero )
enddo
call cfft3 ( vr, nr1, nr2, nr3, 1)
!
deallocate ( ss, vs)
!
! !----------------------------------------------------------------
! ! BRUTE-FORCE DIAGONALIZATION
! !----------------------------------------------------------------
! !
! allocate ( eval(ngm), hk(ngm,ngm), u(ngm,ngm), fv1(ngm), fv2(ngm) )
! !
! do ik = 1, nk
! !
! hk = h(:,:,ik)
! !
! ! 0 = no eigenvectors
! call rs ( ngm, ngm, hk, eval, 0, u, fv1, fv2, ierr)
! !
! write(10,'(i3,10(3x,f13.5))') ik, eval(1:nbnd)*ryd2ev
! !
! enddo
! !
! deallocate ( eval, hk, u, fv1, fv2 )
!
!----------------------------------------------------------------
! CONJUGATE GRADIENTS AND GREEN'S FUNCTION FROM RECURSION METHOD
!----------------------------------------------------------------
!
! I checked that the band-structure from CG is identical to the
! one obtained by a full diagonalization -> FFT and H|psi> are ok.
!
allocate ( e (nbnd), psi (ngm, nbnd) )
!
write(6,'(/4x,a)') repeat('-',67)
write(6,'(4x,"Eigenvalues (eV) from preconditioned conjugate gradient")')
write(6,'(4x,a/)') repeat('-',67)
!
allocate ( eval(ngm0), hk(ngm0,ngm0), u(ngm0,ngm0), fv1(ngm0), fv2(ngm0) )
!
do ik = 1, nk
!
!----------------------------------------------------------------
! CONJUGATE GRADIENTS DIAGONALIZATION
!----------------------------------------------------------------
!
! trial wavefunctions for CG algorithm, cutoff ecut0 (modules.f90)
!
hk = h(1:ngm0,1:ngm0,ik)
call rs ( ngm0, ngm0, hk, eval, 1, u, fv1, fv2, ierr)
psi = zero
do ibnd = 1, nbnd
do ig = 1, ngm0
psi(ig,ibnd) = u(ig,ibnd)
enddo
enddo
!
call cgdiag (nbnd, psi, e, t(1,ik), vr )
!
write(6,'(4x,i3,8(2x,f6.3))') ik, e(1:8)*ryd2ev
write(6,'(4x,a/)') repeat('-',67)
!
!----------------------------------------------------------------
! GREEN'S FUNCTION G(r,r',k,w) IN THE PSINC BASIS
!----------------------------------------------------------------
!
! I put this in the CG loop since I want to compare the GF
! calculated by expanding on empty states with the one obtained
! by the recursion method
!
!@ do ir1 = 1, nr
!@ do ir2 = 1, nr
do ir1 = 134, 134
do ir2 = 1000, 1000
!
!! call green ( ir1, ir2, t(1,ik), vr )
call green ( ir1, ir2, t(1,ik), vr, psi, e ) ! @@@ debug
!
enddo
enddo
!
enddo
!@@
10 continue
!@@
!
!----------------------------------------------------------------
! SCREENED COULOMB INTERACTION
!----------------------------------------------------------------
!
! q point for test purposes - really need a loop here
!
xq(:, 1) = (/ 0.001, 0.0, 0.0/)
xq(:, 2) = (/ 0.010, 0.0, 0.0/)
xq(:, 3) = (/ 0.050, 0.0, 0.0/)
xq(:, 4) = (/ 0.100, 0.0, 0.0/)
xq(:, 5) = (/ 0.150, 0.0, 0.0/)
xq(:, 6) = (/ 0.200, 0.0, 0.0/)
xq(:, 7) = (/ 0.300, 0.0, 0.0/)
xq(:, 8) = (/ 0.500, 0.0, 0.0/)
xq(:, 9) = (/ 0.750, 0.0, 0.0/)
xq(:,10) = (/ 1.000, 0.0, 0.0/)
xq(:,11) = (/ 1.500, 0.0, 0.0/)
xq(:,12) = (/ 2.000, 0.0, 0.0/)
xq(:,13) = (/ 6.000, 0.0, 0.0/)
xq(:,14) = (/ 10.000, 0.0, 0.0/)
xq(:,15) = (/ 2.0, 0.0, 0.0/)
xq(:,16) = (/ 1.0, 1.0, 1.0/)
!
! do iq = 1, nq
! do iq = 1, 14
do iq = 2, 2
call coulomb ( xq(:,iq) )
! call dielec_mat ( xq(:,iq) )
enddo
!
!----------------------------------------------------------------
! SELF-ENERGY
!----------------------------------------------------------------
!
end
!