! !---------------------------------------------------------------- program gwhs !---------------------------------------------------------------- ! ! pilot code the GW-HS method ! empirical pseudopotential method implementation ! for nonpolar tetrahedral semiconductors ! ! 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). ! ! IMPORTANT: every wfs has the same G-vects, I am not changing ! the cutoff for every k-point: |G|^2 is used instead of |k+G|^2 ! so mind when using k beyond the first BZ ! !---------------------------------------------------------------- ! use gspace use parameters use constants use kspace implicit none ! ! variables ! integer :: ig, ik, ig1, ig2, iq, nk, ik0, i, j, k, ios integer :: iw, iwp, iwmwp, count, ipol, ir integer :: recl, unf_recl, irp, igp integer :: nwcoul, nwgreen, nwalloc, nwsigma real(dbl) :: ui, uj, uk, wgreenmin, wgreenmax real(dbl) :: gcutm, gcutms, arg, wp real(dbl) :: k0mq(3), kplusg(3), xk0(3,nk0), xxq(3) real(dbl), allocatable :: ss(:), vs(:) real(DP), parameter :: eps8 = 1.0D-8 real(DP) :: et(nbnd_occ, nq) real(DP), allocatable :: g2kin (:) real(DP), allocatable :: wtmp(:), wcoul(:), wgreen(:), wsigma(:) complex(DP), allocatable :: psi(:,:) complex(dbl), allocatable :: vr(:), aux(:) complex(dbl), allocatable :: scrcoul (:,:,:), greenf (:,:,:), sigma(:,:,:) logical :: allowed ! ! !---------------------------------------------------------------- ! 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) ! !---------------------------------------------------------------- ! FIND THE G-VECTORS FOR THE SMALL SIGMA CUTOFF !---------------------------------------------------------------- ! ! 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 ! gcutms = four * ecuts / tpiba2 ! gcutms = ecuts / tpiba2 ! ! set the fft grid ! ! estimate nr1 and check if it is an allowed value for FFT ! nr1s = 1 + int (2 * sqrt (gcutms) * sqrt( at(1,1)**2 + at(2,1)**2 + at(3,1)**2 ) ) nr2s = 1 + int (2 * sqrt (gcutms) * sqrt( at(1,2)**2 + at(2,2)**2 + at(3,2)**2 ) ) nr3s = 1 + int (2 * sqrt (gcutms) * sqrt( at(1,3)**2 + at(2,3)**2 + at(3,3)**2 ) ) ! do while (.not.allowed(nr1s)) nr1s = nr1s + 1 enddo do while (.not.allowed(nr2s)) nr2s = nr2s + 1 enddo do while (.not.allowed(nr3s)) nr3s = nr3s + 1 enddo ! call ggens ( gcutms ) ! !---------------------------------------------------------------- ! ! read from file the k-points for the self-energy \Sigma(k) ! do ik = 1, nk0 read (5,*) xk0(1,ik), xk0(2,ik), xk0(3,ik) enddo ! ! construct the empirical pseudopotential ! allocate ( ss(ngm), vs(ngm) ) 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 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 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) ! ! allocate ( g2kin (ngm) ) allocate ( aux (nrs) ) allocate ( xq(3,nq), wq(nq), eval_occ(nbnd_occ,nq) ) allocate ( gmap(ngm,27) ) ! ! generate the uniform {q} grid for the Coulomb interaction ! no symmetry-reduction for now - uniform and Gamma-centered ! (I was going insane with the folding of the MP mesh, I am not sure ! it's self-contained) ! count = 0 do i = 1, nq1 ui = (i - 1.d0) / float (nq1) do j = 1, nq2 uj = (j - 1.d0) / float (nq2) do k = 1, nq3 uk = (k - 1.d0) / float (nq3) count = count + 1 xq (:, count) = ui * bg(:,1) + uj * bg(:,2) + uk * bg(:,3) enddo enddo enddo wq = one / float ( count ) if (count.ne.nq) call error ('gwhs','q-point count',count) ! the {k} grid is taken to coincide with the {q} grid ! nks = 2 * nq allocate ( xk (3,nks), wk(nks) ) ! write(6,'(/4x,a)') repeat('-',67) write(6,'(4x,"Uniform q-point grid for the screened Coulomb interaction"/)') do iq = 1, nq write ( 6, '(4x,"q(",i3," ) = (",3f12.7," ), wq =",f12.7)') & iq, (xq (ipol, iq) , ipol = 1, 3) , wq (iq) enddo write(6,'(4x,a/)') repeat('-',67) ! ! generate the occupied eigenstates on the uniform grid ! this will be needed for the screened Coulomb below ! recl = 2 * nbnd_occ * ngm ! 2 stands for complex unf_recl = DIRECT_IO_FACTOR * recl open ( iunwfc, file = "./silicon.wfc", iostat = ios, form = 'unformatted', & status = 'unknown', access = 'direct', recl = unf_recl) ! write(6,'(/4x,a)') repeat('-',67) write(6,'(4x,"Occupied eigenvalues (eV)")') write(6,'(4x,a/)') repeat('-',67) ! allocate ( psi (ngm, nbnd_occ) ) do iq = 1, nq ! ! the k-dependent kinetic energy in Ry ! [Eq. (14) of Ihm,Zunger,Cohen J Phys C 12, 4409 (1979)] ! do ig = 1, ngm kplusg = xq(:, iq) + g(:,ig) g2kin ( ig ) = tpiba2 * dot_product ( kplusg, kplusg ) enddo ! call eigenstates2 ( xq(:, iq), vr, g2kin, psi, eval_occ(:,iq) ) ! ! direct write to file - take into account the k/k+q alternation ! write ( iunwfc, rec = 2 * iq - 1, iostat = ios) psi ! write ( 6, '(4x,"q(",i3," )",10(3x,f7.3))') iq, eval_occ(:,iq)*ryd2ev ! enddo deallocate( psi ) write(6,'(4x,a/)') repeat('-',67) ! ! here we generate the G-map for the folding into the first BZ ! call refold ( ) ! ! now prepare the unit to write the Haydock's coefficients ! recl = 2 * 4 * ( nstep+2 ) + 2 * 4 + 4 ! this correspond to a(nstep+2,4), b(nstep+2,4), a_term(4), b_term(4) and norm(4) ! (all of them are real) unf_recl = DIRECT_IO_FACTOR * recl open ( iuncoeff, file = "./silicon.coeff", iostat = ios, form = 'unformatted', & status = 'unknown', access = 'direct', recl = unf_recl) ! ! ------------------------------------------------ ! MAIN: CALCULATE G AND W AND PERFORM CONVOLUTIONS ! ------------------------------------------------ ! ! ---------------------------------------------------------------- ! generate frequency bins ! ---------------------------------------------------------------- ! ! Here I assume Sigma is needed for w0 between wsigmamin and wsigmamax ! The convolution requires W for positive frequencies w up to wcoulmax ! (even function - cf Shishkin and Kress) and the GF spanning w0+-w. ! Therefore the freq. range of GF is ! from (wsigmamin-wcoulmax) to (wsigmamax+wcoulmax) ! the freq. dependence of the GF is inexpensive, so we use the same spacing ! ! NB: I assume wcoulmax>0, wsigmamin=<0, wsigmamax>0 and zero of energy at the Fermi level ! wgreenmin = wsigmamin-wcoulmax wgreenmax = wsigmamax+wcoulmax ! nwalloc = 1 + ceiling( (wgreenmax-wgreenmin)/deltaw ) allocate ( wtmp(nwalloc), wcoul(nwalloc), wgreen(nwalloc), wsigma(nwalloc) ) wcoul = zero wgreen = zero wsigma = zero ! do iw = 1, nwalloc wtmp(iw) = wgreenmin + (wgreenmax-wgreenmin)/float(nwalloc-1)*float(iw-1) enddo ! align the bins with the zero of energy wtmp = wtmp - minval ( abs ( wgreen) ) ! nwgreen = 0 nwcoul = 0 nwsigma = 0 ! do iw = 1, nwalloc if ( ( wtmp(iw) .ge. wgreenmin ) .and. ( wtmp(iw) .le. wgreenmax) ) then nwgreen = nwgreen + 1 wgreen(nwgreen) = wtmp(iw) endif if ( ( wtmp(iw) .ge. zero ) .and. ( wtmp(iw) .le. wcoulmax) ) then nwcoul = nwcoul + 1 wcoul(nwcoul) = wtmp(iw) endif if ( ( wtmp(iw) .ge. wsigmamin ) .and. ( wtmp(iw) .le. wsigmamax) ) then nwsigma = nwsigma + 1 wsigma(nwsigma) = wtmp(iw) endif enddo ! allocate ( scrcoul (nrs, nrs, nwcoul) ) allocate ( greenf (nrs, nrs, nwgreen) ) allocate ( sigma (nrs, nrs, nwsigma) ) ! ! prepare the unit to write the Coulomb potential ! each q-point is associated with one record ! recl = 2 * nrs * nrs * nwcoul unf_recl = DIRECT_IO_FACTOR * recl open ( iuncoul, file = "./silicon.coul", iostat = ios, form = 'unformatted', & status = 'unknown', access = 'direct', recl = unf_recl) ! ! ! loop over {q} for the screened Coulomb interaction ! do iq = 1, 3 !@ nq ! write(6,'(4x,3x,"iq = ",i3)') iq ! ! the grids {k} and {k+q} for the dVscf will be obtained ! by shuffling the {q} grid ! call coulomb ( xq(:,iq), nwcoul, wcoul, scrcoul ) ! write ( iuncoul, rec = iq, iostat = ios) scrcoul ! write (6,'(4x,"Written scrcoul for iq = ",i3)') iq ! enddo ! ! loop over the {k0} set for the Self-Energy ! do ik0 = 1, 1 !@ nk0 ! write(6,'(4x,"ik0 = ",i3)') ik0 ! ! loop over the {k0-q} grid for the Green's function ! do iq = 1, 3 !@ nq ! write(6,'(4x,3x,"iq = ",i3)') iq ! ! k0mq = k0 - q ! k0mq = xk0(:,ik0) - xq(:,iq) ! ! the k-dependent kinetic energy in Ry ! [Eq. (14) of Ihm,Zunger,Cohen J Phys C 12, 4409 (1979)] ! do ig = 1, ngm kplusg = k0mq + g(:,ig) g2kin ( ig ) = tpiba2 * dot_product ( kplusg, kplusg ) enddo ! call green_coeff ( iq, g2kin, vr, nwgreen, wgreen) ! frequency passed for test purposes call green_fraction ( iq, nwgreen, wgreen ) !@ need to bring greenf() out ! ! now greenf (nrs,nrs,nwgreen) contains the Green's function ! for this k0mq point, all frequencies, and in G-space: !@@ !@@ write on file here !@@ ! ! end loop on {k0-q} and {q} enddo ! ! end loop on {k0} enddo ! ! G TIMES W PRODUCT ! ! *** NOT TESTED *** ! ! sigma = czero do ik0 = 1, 1 !@ nk0 ! ! now sum over {q} the products G(k0-q)W(q) ! do iq = 1, 3 !@ nq ! ! go to R-space to perform the direct product with W ! both indeces are in SIZE order of G-vectors !@ !@ here need to read Green function !@ do iw = 1, 2 !@ nwgreen do ig = 1, ngms aux = czero do igp = 1, ngms aux(nls(igp)) = greenf(igp,ig,iw) enddo call cfft3s ( aux, nr1s, nr2s, nr3s, 1) greenf(:,ig,iw) = aux enddo do ir = 1, nrs aux = czero do igp = 1, ngms aux(nls(igp)) = greenf(ir,igp,iw) enddo call cfft3s ( aux, nr1s, nr2s, nr3s, 1) greenf(ir,:,iw) = aux enddo enddo write(6,'(4x,"FFTW of Green''s function passed"/)') ! read ( iuncoul, rec = iq, iostat = ios) scrcoul ! ! now scrcoul(nrs,nrs,nwcoul) contains the screened Coulomb ! interaction for this q point, all frequencies, and in G-space: ! go to R-space to perform the direct product with G ! the 1st index is in SIZE order of G-vectors ! the 2nd index is in FFT order of G-vectors (look at coulomb.f90) ! do iw = 1, 2 !@ nwcoul do ig = 1, ngms aux = czero do igp = 1, ngms aux(nls(igp)) = scrcoul(igp,ig,iw) enddo call cfft3s ( aux, nr1s, nr2s, nr3s, 1) scrcoul(:,ig,iw) = aux enddo do ir = 1, nrs aux = czero do igp = 1, ngms aux(igp) = scrcoul(ir,igp,iw) enddo call cfft3s ( aux, nr1s, nr2s, nr3s, 1) scrcoul(ir,:,iw) = aux enddo enddo write(6,'(4x,"FFTW of Coulomb passed"/)') ! ! combine Green's function and screened Coulomb ( sum_q wq = 1 ) ! do iw = 1, nwsigma ! do ir = 1, nrs do irp = 1, nrs ! do iwp = 1, nwsigma ! ! iwmwp = ... iwmwp = 4 !@ temporary ! ! wp = ... wp = 2.d0 !@ temporary ! ! ! not sure about the sign of the deltas inside G and W... ! sigma(ir,irp,iw) = sigma(ir,irp,iw) & + wq (iq) * ci / twopi * exp ( -ci * delta * wp ) & * greenf(ir,irp,iwmwp) * scrcoul (ir,irp,iwp) ! enddo ! enddo enddo ! enddo ! ! end loop on {k0-q} and {q} enddo ! ! end loop on {k0} enddo ! ! Now we have summed over q in G(k0-q)W(q) and we can go back ! to G-space before calculating the sandwitches with the wavefunctions ! note: we go to SIZE order of G-vectors ! do iw = 1, 2 !@ nwsigma do ir = 1, nrs aux = czero do irp = 1, nrs aux(irp) = sigma(ir,irp,iw) enddo call cfft3s ( aux, nr1s, nr2s, nr3s, -1) do ig = 1, ngms sigma (ir,ig,iw) = aux(nls(ig)) enddo enddo do ig = 1, ngms aux = czero do irp = 1, nrs aux(irp) = sigma(irp,ig,iw) enddo call cfft3s ( aux, nr1s, nr2s, nr3s, -1) do igp = 1, ngms sigma (igp,ig,iw) = aux(nls(igp)) enddo enddo enddo ! ! everything beyond ngms is garbage ! do ig = ngms+1, nrs do igp = ngms+1, nrs do iw = 1, nwsigma sigma (ig,igp,iw) = czero enddo enddo enddo ! deallocate( vr, g2kin, greenf, scrcoul, sigma, gmap) close (iunwfc, status = 'delete') close (iuncoul, status = 'keep') close (iuncoeff, status = 'keep') 100 format (4x,"Green's function: q G G' status") 101 format (4x,15x,3i6," done") ! end program gwhs !---------------------------------------------------------------- !