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