! !---------------------------------------------------------------- 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, nw_sigma, count, ipol integer :: recl, unf_recl real(dbl) :: ui, uj, uk real(dbl) :: gcutm, 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 (:) complex(DP), allocatable :: psi(:,:) complex(dbl), allocatable :: vr(:) 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) ! !---------------------------------------------------------------- ! ! 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) ! ! !@ nw_sigma = 10 !@ allocate ( psi (ngm, nbnd_occ) ) allocate ( g2kin (ngm) ) allocate ( scrcoul (ngm_sigma, ngm_sigma, fbins) ) allocate ( greenf (ngm_sigma, ngm_sigma, nw) ) allocate ( xq(3,nq), wq(nq), eval(nbnd_occ,nq) ) allocate ( sigma (ngm_sigma, ngm_sigma, nw_sigma) ) allocate ( gmap(ngm,27) ) sigma = czero ! ! 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) ! 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(:,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(:,iq)*ryd2ev ! enddo write(6,'(4x,a/)') repeat('-',67) ! ! here we generate the G-map for the folding into the first BZ ! call refold ( ) ! do ik0 = 1, 2 !@ nk0 ! write(6,'(4x,"ik0 = ",i3)') ik0 ! ! the grids {k} and {k+q} for the dVscf will be obtained ! by shuffling the {q} grid ! ! now loop over {k0-q} grid for the Green's function ! and on {q} for the screened Coulomb interaction ! do iq = 1, 3 !@ nq ! xxq = xq(:,iq) ! write(6,'(4x,3x,"iq = ",i3)') iq ! ! k0mq = k0 - q ! k0mq = xk0(:,ik0) - xxq ! write ( 6, 100) do ig1 = 1, 2 !@ ngm_sigma do ig2 = 1, 2 !@ ngm_sigma ! ! 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 ( ig1, ig2, g2kin, vr, greenf (ig1, ig2, :) ) write ( 6, 101) iq, ig1, ig2 ! enddo enddo write(6,'(4x,a/)') ! call coulomb ( xxq, scrcoul ) ! ! combine Green's function and screened Coulomb ( sum_q wq = 1 ) ! do iw = 1, nw_sigma ! do ig1 = 1, ngm_sigma do ig2 = 1, ngm_sigma ! do iwp = 1, nw_sigma ! ! iwmwp = ... iwmwp = 4 !@ temporary ! ! wp = ... wp = 2.d0 !@ temporary ! ! ! not sure about the sign of the deltas inside G and W... ! ! cannot multiply in G-space... need fftw ! sigma(ig1,ig2,iw) = sigma(ig1,ig2,iw) & + wq (iq) * ci / twopi * exp ( -ci * delta * wp ) & * greenf(ig1,ig2,iwmwp) * scrcoul (ig1,ig2,iwp) ! enddo ! enddo enddo ! enddo ! ! end loop on {k0-q} and {q} enddo ! ! end loop on {k0} enddo ! deallocate( vr, psi, g2kin, greenf, scrcoul, sigma, gmap) close (iunwfc, status = 'delete') 100 format (4x,"Green's function: q G G' status") 101 format (4x,15x,3i6," done") ! end program gwhs !---------------------------------------------------------------- !