1!=========================================================================== 2! 3! Routines() 4! 5! (1) shiftenergy() Originally by ? Last Edited: Apr/2016 (FHJ) 6! 7! Computes and symmetrizes the quasiparticle spectrum, and 8! perform linear interpolation/extrapolation to solve Dyson`s equation. 9! Note that, since we may use more than two frequency points to evaluate 10! Sigma(omega), our "eqp1" in general goes beyond the linear correction 11! scheme in Eq. (37) of Hybertsen & Louie PRB. 12! 13!=========================================================================== 14 15#include "f_defs.h" 16 17subroutine shiftenergy(sig, wfnk, alda, asx, ach, achcor, ach_n1, achcor_n1, & 18 ax, efsto, asig, enew, zrenorm, nfreqgpp, ncore_excl) 19 20 use global_m 21 implicit none 22 23 type (siginfo), intent(in) :: sig 24 type (wfnkstates), intent(in) :: wfnk 25 integer, intent(in) :: nfreqgpp 26 SCALAR, intent(inout) :: & 27 alda(sig%ndiag+sig%noffdiag,sig%nspin), & !< vxc 28 asx(nfreqgpp,sig%ndiag+sig%noffdiag,sig%nspin), & !< sx 29 ach(nfreqgpp,sig%ndiag+sig%noffdiag,sig%nspin) !< ch 30 complex(DPC), intent(inout) :: achcor(sig%ndiag+sig%noffdiag,sig%nspin) !< static remainder 31 SCALAR, intent(inout) :: ach_n1(sig%ntband,sig%ndiag+sig%noffdiag,sig%nspin) 32 SCALAR, intent(inout) :: achcor_n1(sig%ntband,sig%ndiag+sig%noffdiag,sig%nspin) 33 SCALAR, intent(inout) :: ax(sig%ndiag+sig%noffdiag,sig%nspin) !< x 34 real(DP), intent(out) :: efsto(sig%ndiag,sig%nspin) !< eqp0 35 SCALAR, intent(inout) :: asig(sig%ndiag+sig%noffdiag,sig%nspin) !< sig 36 real(DP), intent(out) :: enew(sig%ndiag,sig%nspin) !< eqp1 37 real(DP), intent(out) :: zrenorm(sig%ndiag,sig%nspin) !< Znk 38 integer, intent(in) :: ncore_excl !< number of core states excluded 39 40 integer :: ii,jj,istart,istop,nsubs,ispin,iwlda0,iwlda1,iwlda2,iw 41 integer :: ndeg(sig%ndiag) 42 real(DP) :: fact,dek,dele,diffmin,delta_omega,diff,e_lk,freq0 43 SCALAR :: aldai,axi 44 SCALAR, dimension(nfreqgpp) :: asigi, asxi, achi 45 SCALAR, dimension(sig%ntband) :: ach_n1i, achcor_n1i 46 complex(DPC) :: achcori 47 48 PUSH_SUB(shiftenergy) 49 50 do ispin=1,sig%nspin 51 52 nsubs = 1 53 ndeg(nsubs) = 1 54 ! This loop is used to find degeneracy 55 do ii=2,sig%ndiag 56! DVF : ncore_excl has to be substracted here because wfnk%elda is defined in the 57! read_wavefunction subroutine in input.f90 to be referenced to the case with 58! no core states. Same with wfnk%el/elda below. 59 dek = wfnk%elda(sig%diag(ii)-ncore_excl,ispin) - wfnk%elda(sig%diag(ii-1)-ncore_excl,ispin) 60 if (abs(dek)<sig%tol .and. sig%symmetrize) then 61 ! Band ii is degenerate to ii-1 62 ndeg(nsubs) = ndeg(nsubs) + 1 63 else 64 ! Start a new degenerate subspace 65 nsubs = nsubs + 1 66 ndeg(nsubs) = 1 67 endif 68 enddo 69 70 istop = 0 71 ! loop degeneracy 72 do ii=1,nsubs 73 istart = istop + 1 74 istop = istart + ndeg(ii) - 1 75 aldai = ZERO 76 axi = ZERO 77 asxi = ZERO 78 achi = ZERO 79 achcori = (0.0d0, 0.0d0) 80 ach_n1i(:) = ZERO 81 achcor_n1i(:) = ZERO 82 ! start real loop over bands 83 do jj=istart,istop 84 ! average VXC in the same degen group: combine first 85 aldai = aldai + alda(jj,ispin) 86 ! "i" is intermediate steps in Ry 87 axi = axi + ax(jj,ispin) 88 asxi(:) = asxi(:) + asx(:,jj,ispin) 89 achi(:) = achi(:) + ach(:,jj,ispin) 90 achcori = achcori + achcor(jj,ispin) 91 ach_n1i(:) = ach_n1i(:) + ach_n1(:,jj,ispin) 92 achcor_n1i(:) = achcor_n1i(:) + achcor_n1(:,jj,ispin) 93 enddo 94 95 fact = ryd / dble(ndeg(ii)) 96 do jj=istart,istop 97 ! then AVERAGE 98 ! common quantities in eV 99 alda(jj,ispin) = aldai * fact 100 ax(jj,ispin) = axi * fact 101 asx(:,jj,ispin) = asxi(:) * fact 102 ach(:,jj,ispin) = achi(:) * fact 103 achcor(jj,ispin) = achcori * fact 104 ach_n1(:,jj,ispin) = ach_n1i(:) * fact 105 achcor_n1(:,jj,ispin) = achcor_n1i(:) * fact 106 ! intermediate the quantity, sigma in Ry, unnecessary i 107 asigi(:) = ax(jj,ispin) + asx(:,jj,ispin) + ach(:,jj,ispin) 108 109 110 ! FHJ: Determine the expansion points for the linear interpolation: iwlda0 111 ! is where we evaluate eqp0, and iwlda1, iwlda2 are the expansion points 112 ! we used in the finite difference linearization (ie, used to compute the 113 ! derivatives) 114 if (sig%fdf==-3) then 115 ! FHJ: Find iw closest to e_lk 116 iwlda1 = 1 117 iwlda2 = 2 118 diffmin = INF 119 e_lk = wfnk%ek(sig%diag(jj)-ncore_excl,ispin) 120 ! FHJ: Figure out starting frequency for freq. grid 121 if (sig%freq_grid_shift<2) then 122 freq0 = sig%freqevalmin 123 else 124 freq0 = e_lk - sig%freqevalstep*(sig%nfreqeval-1)/2 125 endif 126 do iw=1,sig%nfreqeval 127 diff = abs(freq0 + (iw-1)*sig%freqevalstep - e_lk) 128 if (diff .lt. diffmin) then 129 diffmin=diff 130 iwlda2=iwlda1 131 iwlda1=iw 132 endif 133 enddo 134 iwlda0 = iwlda1 135 delta_omega = (iwlda2 - iwlda1)*sig%freqevalstep 136 else 137 !iwlda0 is always 2 here 138 iwlda0 = 2 139 iwlda1 = 2 140 iwlda2 = 2 141 select case (sig%fdf) 142 case (-1) 143 iwlda1 = 1 144 iwlda2 = 2 145 case (0) 146 iwlda1 = 1 147 iwlda2 = 3 148 case (1,2) 149 iwlda1 = 2 150 iwlda2 = 3 151 end select 152 delta_omega = (iwlda2 - iwlda1)*sig%dw 153 endif 154 asig(jj,ispin) = asigi(iwlda0) 155 156!FHJ: Perform linear interpolation/extrapolation for all finite difference 157! schemes at once. enew is the interpolated off-shell answer = "eqp1", 158! and efsto is the on-shell answer = "eqp0". 159 160 efsto(jj,ispin) = wfnk%elda(sig%diag(jj)-ncore_excl,ispin) - & 161 alda(jj,ispin) + asig(jj,ispin) 162 dele = efsto(jj,ispin) - wfnk%ek(sig%diag(jj)-ncore_excl,ispin) 163 164 if (iwlda1/=iwlda2) then 165 enew(jj,ispin) = efsto(jj,ispin) + & 166 (asigi(iwlda2)-asigi(iwlda1)) / (delta_omega - asigi(iwlda2) + asigi(iwlda1)) * dele 167 zrenorm(jj, ispin) = 1d0 / (1d0 - (asigi(iwlda2) - asigi(iwlda1))/delta_omega) 168 else 169 enew(jj,ispin) = efsto(jj,ispin) 170 zrenorm(jj, ispin) = 1d0 171 endif 172 173 enddo ! jj 174 enddo ! ii 175 enddo ! ispin 176 177 POP_SUB(shiftenergy) 178 179 return 180end subroutine shiftenergy 181