1!
2! Copyright (C) 2003 PWSCF group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9!-----------------------------------------------------------------------
10SUBROUTINE A_h(npw,e,h,ah)
11  !-----------------------------------------------------------------------
12  USE kinds, ONLY: DP
13  USE cell_base,ONLY : alat, omega, tpiba2
14  USE uspp,     ONLY : vkb, nkb
15  USE lsda_mod, ONLY : current_spin, nspin
16  USE wvfct, ONLY: nbnd, npwx, g2kin
17  USE wavefunctions,  ONLY: evc, psic
18  USE scf,      ONLY : vrs, rho
19  USE fft_base, ONLY : dffts, dfftp
20  USE fft_interfaces, ONLY : fwfft, invfft
21  USE gvect,    ONLY : gstart, g, gg
22  USE constants,  ONLY: degspin, e2, fpi
23  USE becmod, ONLY: bec_type, becp, calbec
24  USE gc_lr, ONLY:  grho, dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s
25  USE cgcom, ONLY: auxr, aux2, aux3, dmuxc
26  USE funct, ONLY: dft_is_gradient
27  !
28  IMPLICIT NONE
29  INTEGER :: npw, j, jkb, ibnd, na,nt,ih
30  real(DP) :: e(nbnd)
31  COMPLEX(DP) :: h(npwx,nbnd), ah(npwx,nbnd)
32  !
33  COMPLEX(DP) :: fp, fm
34  COMPLEX(DP), POINTER :: dpsic(:), drhoc(:)
35  REAL(dp), allocatable :: dv(:)
36  real(DP), POINTER  :: drho(:)
37  !
38  CALL start_clock('a_h')
39  !
40  drho  => auxr
41  dpsic => aux2
42  drhoc => aux3
43  !
44  drho(:) = 0.d0
45  !
46  ! [(k+G)^2 - e ]psi
47  DO ibnd = 1,nbnd
48     ! set to zero the imaginary part of h at G=0
49     ! needed for numerical stability
50     IF (gstart==2) h(1,ibnd) = cmplx( dble(h(1,ibnd)),0.d0,kind=DP)
51     DO j = 1,npw
52        ah(j,ibnd) = (g2kin(j)-e(ibnd)) * h(j,ibnd)
53     ENDDO
54  ENDDO
55  !     V_Loc psi
56  DO ibnd = 1,nbnd, 2
57     dpsic(:)= (0.d0, 0.d0)
58     psic(:) = (0.d0, 0.d0)
59     IF (ibnd<nbnd) THEN
60        ! two ffts at the same time
61        DO j = 1,npw
62           psic (dfftp%nl (j)) = evc(j,ibnd) + (0.0d0,1.d0)* evc(j,ibnd+1)
63           dpsic(dfftp%nl (j)) =   h(j,ibnd) + (0.0d0,1.d0)*   h(j,ibnd+1)
64           psic (dfftp%nlm(j))= conjg(evc(j,ibnd)-(0.0d0,1.d0)* evc(j,ibnd+1))
65           dpsic(dfftp%nlm(j))= conjg(  h(j,ibnd)-(0.0d0,1.d0)*   h(j,ibnd+1))
66        ENDDO
67     ELSE
68        DO j = 1,npw
69           psic (dfftp%nl (j)) = evc(j,ibnd)
70           dpsic(dfftp%nl (j)) =   h(j,ibnd)
71           psic (dfftp%nlm(j)) = conjg( evc(j,ibnd))
72           dpsic(dfftp%nlm(j)) = conjg(   h(j,ibnd))
73        ENDDO
74     ENDIF
75     CALL invfft ('Wave', psic, dffts)
76     CALL invfft ('Wave',dpsic, dffts)
77     DO j = 1,dfftp%nnr
78        drho(j) = drho(j) - 2.0d0*degspin/omega *   &
79              dble(psic(j)*conjg(dpsic(j)))
80        dpsic(j) = dpsic(j) * vrs(j,current_spin)
81     ENDDO
82     CALL fwfft ('Wave',dpsic, dffts)
83     IF (ibnd<nbnd) THEN
84        ! two ffts at the same time
85        DO j = 1,npw
86           fp = (dpsic (dfftp%nl(j)) + dpsic (dfftp%nlm(j)))*0.5d0
87           fm = (dpsic (dfftp%nl(j)) - dpsic (dfftp%nlm(j)))*0.5d0
88           ah(j,ibnd  ) = ah(j,ibnd)  +cmplx( dble(fp), aimag(fm),kind=DP)
89           ah(j,ibnd+1) = ah(j,ibnd+1)+cmplx(aimag(fp),- dble(fm),kind=DP)
90        ENDDO
91     ELSE
92        DO j = 1,npw
93           ah(j,ibnd) = ah(j,ibnd)  + dpsic (dfftp%nl(j))
94        ENDDO
95     ENDIF
96  ENDDO
97  !
98  ! V_NL psi
99  CALL calbec ( npw, vkb, h, becp)
100  IF (nkb > 0) CALL add_vuspsi (npwx, npw, nbnd, ah)
101  !
102  DO j = 1,dfftp%nnr
103     drhoc(j) = cmplx(drho(j),0.d0,kind=DP)
104  ENDDO
105  CALL fwfft ('Rho', drhoc, dfftp)
106  DO j = 1,dfftp%ngm
107     dpsic(j) = drhoc(dfftp%nl(j))
108  ENDDO
109  !
110  ! drho is deltarho(r)
111  ! drhoc is deltarho(g) on the FFT grid
112  ! dpsic is deltarho(g) on the G-vector grid
113  !
114  !  mu'(n(r)) psi(r) delta psi(r)
115  !
116  ALLOCATE (dv(dfftp%nnr))
117  DO j = 1,dfftp%nnr
118     dv(j) = drho(j)*dmuxc(j)
119  ENDDO
120  !
121  !  add gradient correction contribution (if any)
122  !
123  CALL start_clock('dgradcorr')
124  IF (dft_is_gradient() ) THEN
125     !
126     CALL dgradcor1  &
127         (dfftp, rho%of_r, grho, dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s,  &
128          drho, dpsic, nspin, g, dv)
129     !
130  ENDIF
131  CALL stop_clock('dgradcorr')
132  NULLIFY(dpsic)
133  NULLIFY (drho)
134  !
135  !  1/|r-r'| * psi(r') delta psi(r')
136  !
137  ! gstart is the first nonzero G vector (needed for parallel execution)
138  !
139  IF (gstart==2) drhoc(dfftp%nl(1)) = 0.d0
140  !
141  DO j = gstart,dfftp%ngm
142     drhoc(dfftp%nl (j)) = e2*fpi*drhoc(dfftp%nl(j))/ (tpiba2*gg(j))
143     drhoc(dfftp%nlm(j)) = conjg(drhoc(dfftp%nl (j)))
144  ENDDO
145  CALL invfft ('Rho', drhoc, dfftp)
146  !
147  ! drhoc now contains deltaV_hartree
148  !
149  DO j = 1,dfftp%nnr
150     dv(j) = - dv(j) - dble(drhoc(j))
151  ENDDO
152  !
153  CALL vloc_psi_gamma(npwx, npw, nbnd, evc, dv, ah)
154  !
155  NULLIFY(drhoc)
156  DEALLOCATE (dv)
157  !
158  ! set to zero the imaginary part of ah at G=0
159  ! needed for numerical stability
160  IF (gstart==2) THEN
161     DO ibnd = 1, nbnd
162        ah(1,ibnd) = cmplx( dble(ah(1,ibnd)),0.d0,kind=DP)
163     ENDDO
164  ENDIF
165  !
166  CALL stop_clock('a_h')
167  !
168  RETURN
169END SUBROUTINE A_h
170