1!
2! Copyright (C) 2001-2016 Quantum ESPRESSO 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!
8MODULE dv_of_drho_lr
9
10CONTAINS
11
12!-----------------------------------------------------------------------
13subroutine dv_of_drho (dvscf, add_nlcc, drhoc)
14  !-----------------------------------------------------------------------
15  !
16  !  This routine computes the change of the self consistent potential
17  !  (Hartree and XC) due to the perturbation.
18  !  Note: gamma_only is disregarded for PHonon calculations,
19  !  TDDFPT purposes only.
20  !
21  USE kinds,             ONLY : DP
22  USE constants,         ONLY : e2, fpi
23  USE fft_base,          ONLY : dfftp
24  USE fft_interfaces,    ONLY : fwfft, invfft
25  USE gvect,             ONLY : ngm, g, gstart
26  USE cell_base,         ONLY : tpiba2, omega
27  USE noncollin_module,  ONLY : nspin_lsda, nspin_mag, nspin_gga
28  USE funct,             ONLY : dft_is_gradient, dft_is_nonlocc
29  USE scf,               ONLY : rho, rho_core
30  USE uspp,              ONLY : nlcc_any
31  USE control_flags,     ONLY : gamma_only
32  USE martyna_tuckerman, ONLY : wg_corr_h, do_comp_mt
33  USE Coul_cut_2D,       ONLY : do_cutoff_2D
34  USE Coul_cut_2D_ph,    ONLY : cutoff_dv_of_drho
35  USE qpoint,            ONLY : xq
36  USE gc_lr,             ONLY : grho, dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s
37  USE control_lr,        ONLY : lrpa
38  USE eqv,               ONLY : dmuxc
39
40  IMPLICIT NONE
41  COMPLEX(DP), INTENT(INOUT) :: dvscf(dfftp%nnr, nspin_mag)
42  ! input:  response charge density
43  ! output: response Hartree-and-XC potential
44  LOGICAL, INTENT(IN) :: add_nlcc
45  ! input: if true add core charge density
46  COMPLEX(DP), INTENT(IN), OPTIONAL :: drhoc(dfftp%nnr)
47  ! input: response core charge density
48  ! (needed only for PHonon when add_nlcc=.true.)
49
50  INTEGER :: ir, is, is1, ig
51  ! counter on r vectors
52  ! counter on spin polarizations
53  ! counter on g vectors
54  REAL(DP) :: qg2, fac, eh_corr
55  ! qg2: the modulus of (q+G)^2
56  ! fac: the structure factor
57  ! eh_corr: the correction to response Hartree energy due
58  ! to Martyna-Tuckerman correction (calculated, but not used).
59  COMPLEX(DP), ALLOCATABLE :: dvaux(:,:), dvhart(:,:), &
60                              dvaux_mt(:), rgtot(:)
61  ! dvaux: response XC potential
62  ! dvhart: response Hartree potential
63  ! dvaux_mt: auxiliary array for Martyna-Tuckerman correction
64  ! rgtot: total response density
65
66  CALL start_clock ('dv_of_drho')
67  !
68  allocate (dvaux( dfftp%nnr, nspin_mag))
69  dvaux (:,:) = (0.d0, 0.d0)
70  !
71  if (add_nlcc .and. .not.present(drhoc)) &
72     & CALL errore( 'dv_of_drho', 'drhoc is not present in the input of the routine', 1 )
73  !
74  ! 1) The exchange-correlation contribution is computed in real space
75  !
76  if (lrpa) goto 111
77  !
78  fac = 1.d0 / DBLE (nspin_lsda)
79  !
80  if (nlcc_any.and.add_nlcc) then
81     rho%of_r(:, 1) = rho%of_r(:, 1) + rho_core (:)
82     do is = 1, nspin_lsda
83        dvscf(:, is) = dvscf(:, is) + fac * drhoc (:)
84     enddo
85  endif
86  !
87  do is = 1, nspin_mag
88     do is1 = 1, nspin_mag
89        do ir = 1, dfftp%nnr
90           dvaux(ir,is) = dvaux(ir,is) + dmuxc(ir,is,is1) * dvscf(ir,is1)
91        enddo
92     enddo
93  enddo
94  !
95  ! Add gradient correction to the response XC potential.
96  ! NB: If nlcc=.true. we need to add here its contribution.
97  ! grho contains already the core charge
98  !
99  if ( dft_is_gradient() ) call dgradcorr(dfftp, rho%of_r, grho, dvxc_rr, &
100                                dvxc_sr, dvxc_ss, dvxc_s, xq, dvscf, &
101                                nspin_mag, nspin_gga, g, dvaux)
102  !
103  if ( dft_is_nonlocc() )  call dnonloccorr(rho%of_r, dvscf, xq, dvaux)
104  !
105  if (nlcc_any.and.add_nlcc) then
106     rho%of_r(:, 1) = rho%of_r(:, 1) - rho_core (:)
107     do is = 1, nspin_lsda
108        dvscf(:, is) = dvscf(:, is) - fac * drhoc (:)
109     enddo
110  endif
111  !
112111 continue
113  !
114  ! Copy the total (up+down) delta rho in dvscf(*,1) and go to G-space
115  !
116  if (nspin_mag == 2) then
117     dvscf(:,1) = dvscf(:,1) + dvscf(:,2)
118  end if
119  !
120  CALL fwfft ('Rho', dvscf(:,1), dfftp)
121  !
122  ! 2) Hartree contribution is computed in reciprocal space
123  !
124  IF (do_comp_mt) THEN
125      !
126      ! Response Hartree potential with the Martyna-Tuckerman correction
127      !
128      allocate(dvhart(dfftp%nnr,nspin_mag))
129      dvhart(:,:) = (0.d0,0.d0)
130      !
131      do is = 1, nspin_lsda
132        do ig = gstart, ngm
133          qg2 = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2 + (g(3,ig)+xq(3))**2
134          dvhart(dfftp%nl(ig),is) = e2 * fpi * dvscf(dfftp%nl(ig),1) / (tpiba2 * qg2)
135        enddo
136      enddo
137      !
138      ! Add Martyna-Tuckerman correction to response Hartree potential
139      !
140      allocate( dvaux_mt( ngm ), rgtot(ngm) )
141      !
142      ! Total response density
143      !
144      do ig = 1, ngm
145         rgtot(ig) = dvscf(dfftp%nl(ig),1)
146      enddo
147      !
148      CALL wg_corr_h (omega, ngm, rgtot, dvaux_mt, eh_corr)
149      !
150      do is = 1, nspin_lsda
151        !
152        do ig = 1, ngm
153           dvhart(dfftp%nl(ig),is)  = dvhart(dfftp%nl(ig),is)  + dvaux_mt(ig)
154        enddo
155        if (gamma_only) then
156           do ig = 1, ngm
157              dvhart(dfftp%nlm(ig),is) = conjg(dvhart(dfftp%nl(ig),is))
158           enddo
159        endif
160        !
161        ! Transform response Hartree potential to real space
162        !
163        CALL invfft ('Rho', dvhart (:,is), dfftp)
164        !
165      enddo
166      !
167      ! At the end the two contributions (XC+Hartree) are added
168      !
169      dvscf = dvaux + dvhart
170      !
171      deallocate( dvaux_mt, rgtot )
172      deallocate(dvhart)
173      !
174  ELSE
175   !
176   ! Response Hartree potential (without Martyna-Tuckerman correction)
177   !
178   If (gamma_only) then
179      !
180      ! Gamma_only case
181      !
182      allocate(dvhart(dfftp%nnr,nspin_mag))
183      dvhart(:,:) = (0.d0,0.d0)
184      !
185      do is = 1, nspin_lsda
186        do ig = 1, ngm
187           qg2 = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2 + (g(3,ig)+xq(3))**2
188           if (qg2 > 1.d-8) then
189              dvhart(dfftp%nl(ig),is) = e2 * fpi * dvscf(dfftp%nl(ig),1) / (tpiba2 * qg2)
190              dvhart(dfftp%nlm(ig),is)=conjg(dvhart(dfftp%nl(ig),is))
191           endif
192        enddo
193        !
194        ! Transformed back to real space
195        !
196        CALL invfft ('Rho', dvhart (:, is), dfftp)
197        !
198      enddo
199      !
200      ! At the end the two contributes are added
201      !
202      dvscf = dvaux + dvhart
203      !
204      deallocate(dvhart)
205      !
206   else
207      !
208      ! General k points implementation
209      !
210      do is = 1, nspin_lsda
211         CALL fwfft ('Rho', dvaux (:, is), dfftp)
212         IF (do_cutoff_2D) THEN
213            call cutoff_dv_of_drho(dvaux, is, dvscf)
214         ELSE
215            do ig = 1, ngm
216               qg2 = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2 + (g(3,ig)+xq(3))**2
217               if (qg2 > 1.d-8) then
218                  dvaux(dfftp%nl(ig),is) = dvaux(dfftp%nl(ig),is) + &
219                                 & e2 * fpi * dvscf(dfftp%nl(ig),1) / (tpiba2 * qg2)
220               endif
221            enddo
222         ENDIF
223         !
224         ! Transformed back to real space
225         !
226         CALL invfft ('Rho', dvaux (:, is), dfftp)
227         !
228      enddo
229      !
230      ! At the end the two contributes are added
231      !
232      dvscf (:,:) = dvaux (:,:)
233      !
234   endif
235   !
236  ENDIF
237  !
238  deallocate (dvaux)
239  !
240  CALL stop_clock ('dv_of_drho')
241  !
242  RETURN
243  !
244end subroutine dv_of_drho
245
246END MODULE dv_of_drho_lr
247