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