1! 2! Copyright (C) 2001-2018 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! 8!-------------------------------------------------------------- 9subroutine zstar_eu_us 10 !-------------------------------------------------------------- 11 ! 12 ! Calculates the additional part of the Born effective charges 13 ! in the case of USPP. 14 ! 15 USE kinds, ONLY : DP 16 USE mp, ONLY : mp_sum 17 USE mp_pools, ONLY : inter_pool_comm 18 USE mp_bands, ONLY : intra_bgrp_comm 19 USE cell_base, ONLY : omega 20 USE ions_base, ONLY : nat, ntyp => nsp, ityp 21 USE buffers, ONLY : get_buffer 22 USE klist, ONLY : xk, wk, ngk, igk_k 23 USE gvecs, ONLY : doublegrid 24 USE fft_base, ONLY : dfftp, dffts 25 USE fft_interfaces, ONLY : fft_interpolate 26 USE lsda_mod, ONLY : nspin, current_spin, isk, lsda 27 USE uspp, ONLY : okvan, nkb, vkb, nlcc_any 28 USE wvfct, ONLY : nbnd, npwx 29 USE paw_variables, ONLY : okpaw 30 USE wavefunctions, ONLY : evc 31 USE uspp_param, ONLY : upf, nhm, nh 32 USE noncollin_module, ONLY : noncolin, npol, nspin_mag 33 USE efield_mod, ONLY : zstareu0 34 USE phus, ONLY : becsumort 35 USE modes, ONLY : u, npert, nirr 36 USE units_ph, ONLY : lrdwf, iucom, lrcom, lrebar, iuebar, & 37 lrdrhous, iudrhous, iudwf 38 USE units_lr, ONLY : iuwfc, lrwfc 39 USE mp_pools, ONLY : nproc_pool, npool 40 USE control_lr, ONLY : nbnd_occ 41 USE lrus, ONLY : int3, int3_paw 42 USE eqv, ONLY : dvpsi, dpsi 43 USE qpoint, ONLY : nksq, ikks 44 USE ldaU, ONLY : lda_plus_u 45 USE dv_of_drho_lr 46 ! 47 implicit none 48 integer :: npw, ibnd, jbnd, ipol, jpol, imode0, irr, imode, nrec, mode 49 integer :: ik, ikk, ig, ir, is, i, j, mu, ipert 50 integer :: ih, jh, ijh 51 integer :: iuhxc, lrhxc 52 ! 53 real(DP) :: weight, fact 54 ! 55 complex(DP), allocatable :: dbecsum(:,:,:,:), aux1 (:) 56 COMPLEX(DP), ALLOCATABLE :: dbecsum_nc(:, :, :, :, :) 57 ! the becsum with dpsi 58 ! auxillary work space for fft 59 complex(DP) , pointer :: & 60 dvscf(:,:,:) 61 complex(DP), allocatable :: pdsp(:,:) 62 complex(DP), allocatable :: drhoscfh (:,:) 63 complex(DP), allocatable :: dvkb (:,:,:) 64 integer :: npe, irr1, imode1, na, nt 65 66#ifdef TIMINIG_ZSTAR_US 67 call start_clock('zstar_eu_us') 68 call start_clock('zstar_us_1') 69#endif 70 71 ! auxiliary space for <psi|ds/du|psi> 72 allocate (dvscf( dfftp%nnr , nspin_mag, 3)) 73 allocate (dbecsum( nhm*(nhm+1)/2, nat, nspin_mag, 3)) 74 if (noncolin) allocate (dbecsum_nc( nhm, nhm, nat, nspin, 3)) 75 allocate (aux1( dffts%nnr)) 76 allocate (pdsp(nbnd,nbnd)) 77 78 ! 79 ! Set the initial values to zero 80 ! 81 pdsp = (0.d0,0.d0) 82 dvscf = (0.d0,0.d0) 83 dbecsum = (0.d0,0.d0) 84 if (noncolin) dbecsum_nc=(0.d0,0.d0) 85 ! 86 ! first we calculate the perturbed charge density and the perturbed 87 ! Hartree and exchange and correlation potential , which we need later 88 ! for the calculation of the Hartree and xc part 89 ! 90 do ik = 1, nksq 91 ikk = ikks(ik) 92 npw = ngk(ikk) 93 if (nksq.gt.1) call get_buffer (evc, lrwfc, iuwfc, ikk) 94 if (lsda) current_spin = isk (ikk) 95 call init_us_2 (npw, igk_k(1,ikk), xk(1,ikk), vkb) 96 weight = wk (ikk) 97 do jpol = 1, 3 98 nrec = (jpol - 1) * nksq + ik 99 call get_buffer(dpsi, lrdwf, iudwf, nrec) 100 if (noncolin) then 101 call incdrhoscf_nc (dvscf(1,1,jpol),weight,ik, & 102 dbecsum_nc(1,1,1,1,jpol), dpsi, 1.0d0) 103 else 104 call incdrhoscf (dvscf(1,current_spin,jpol),weight,ik, & 105 dbecsum(1,1,current_spin,jpol), dpsi) 106 endif 107 end do 108 end do 109 110 IF (noncolin) THEN 111 call mp_sum ( dbecsum_nc, intra_bgrp_comm ) 112 ELSE 113 call mp_sum ( dbecsum, intra_bgrp_comm ) 114 END IF 115#ifdef TIMINIG_ZSTAR_US 116 call stop_clock('zstar_us_1') 117 call start_clock('zstar_us_2') 118#endif 119 120 if (doublegrid) then 121 do is = 1, nspin_mag 122 do ipol = 1, 3 123 call fft_interpolate(dffts, dvscf(:,is,ipol), dfftp, dvscf(:,is,ipol)) 124 end do 125 end do 126 end if 127 128 IF (noncolin.and.okvan) CALL set_dbecsum_nc(dbecsum_nc, dbecsum, 3) 129 130 call addusddense (dvscf, dbecsum) 131 132 call mp_sum ( dvscf, inter_pool_comm ) 133 134#ifdef TIMINIG_ZSTAR_US 135 call stop_clock('zstar_us_2') 136 call start_clock('zstar_us_3') 137#endif 138 139 if (nlcc_any) call addnlcc_zstar_eu_us (dvscf) 140 141 do ipol = 1, 3 142 ! 143 ! Instead of recalculating the perturbed charge density, 144 ! it can also be read from file 145 ! NB: Then the variable fildrho must be set 146 ! 147 ! call davcio_drho(dvscf(1,1,ipol),lrdrho,iudrho,ipol,-1) 148 ! 149 call dv_of_drho (dvscf (:, :, ipol), .false.) 150 enddo 151 call psyme (dvscf) 152 153#ifdef TIMINIG_ZSTAR_US 154 call stop_clock('zstar_us_3') 155 call start_clock('zstar_us_4') 156#endif 157! 158! Calculate the parts with the perturbed Hartree and exchange and correlation 159! potenial 160! 161 imode0 = 0 162 allocate(drhoscfh(dfftp%nnr,nspin_mag)) 163 do irr = 1, nirr 164 npe = npert(irr) 165 do imode = 1, npe 166 mode = imode0 + imode 167 call get_buffer(drhoscfh, lrdrhous, iudrhous, mode) 168 do jpol = 1, 3 169 do is=1,nspin_mag 170 zstareu0(jpol,mode) = zstareu0(jpol,mode) - & 171 dot_product(dvscf(1:dfftp%nnr,is,jpol),drhoscfh(1:dfftp%nnr,is)) & 172 * omega / DBLE(dfftp%nr1*dfftp%nr2*dfftp%nr3) 173 end do 174 end do 175 end do 176 imode0 = imode0 + npe 177 end do 178 deallocate (drhoscfh) 179#ifdef TIMINIG_ZSTAR_US 180 call stop_clock('zstar_us_4') 181 call start_clock('zstar_us_5') 182#endif 183 ! 184 ! Calculate the part with the position operator 185 ! 186 allocate (dvkb(npwx,nkb,3)) 187 do ik = 1, nksq 188 ikk = ikks(ik) 189 npw = ngk(ikk) 190 weight = wk (ikk) 191 if (nksq.gt.1) call get_buffer (evc, lrwfc, iuwfc, ikk) 192 call init_us_2 (npw, igk_k(1,ikk), xk (1, ikk), vkb) 193 call dvkb3(ik, dvkb) 194 imode0 = 0 195 do irr = 1, nirr 196 do imode = 1, npert (irr) 197 mode = imode+imode0 198 do jpol = 1, 3 199 dvpsi = (0.d0,0.d0) 200 ! 201 ! read the Commutator+add. terms 202 ! 203 nrec = (jpol - 1) * nksq + ik 204 call get_buffer(dvpsi, lrebar, iuebar, nrec) 205 ! 206 pdsp = (0.d0,0.d0) 207 call psidspsi (ik, u (1, mode), pdsp ) 208#if defined(__MPI) 209 call mp_sum( pdsp, intra_bgrp_comm ) 210#endif 211 ! 212 ! DFPT+U: add to dvpsi the scf term 213 ! dV_Hub(is)/dE_jpol|psi(ik)> contributing to 214 ! the derivative of the eigenvalue w.r.t. E_jpol 215 ! 216 if (lda_plus_u) call adddvhubscf (jpol, ik) 217 ! 218 ! add the term of the double summation 219 ! 220 do ibnd = 1, nbnd_occ(ikk) 221 do jbnd = 1, nbnd_occ(ikk) 222 zstareu0(jpol,mode)=zstareu0(jpol, mode) + & 223 weight * & 224 dot_product(evc(1:npwx*npol,ibnd), & 225 dvpsi(1:npwx*npol,jbnd))*pdsp(jbnd,ibnd) 226 enddo 227 enddo 228 dvpsi = (0.d0,0.d0) 229 dpsi = (0.d0,0.d0) 230 ! 231 ! For the last part, we read the commutator from disc, 232 ! but this time we calculate 233 ! dS/du P_c [H-eS]|psi> + (dK(r)/du - dS/du)r|psi> 234 ! 235 ! first we read P_c [H-eS]|psi> and store it in dpsi 236 ! 237 nrec = (jpol - 1) * nksq + ik 238 call get_buffer(dpsi, lrcom, iucom, nrec) 239 ! 240 ! Apply the matrix dS/du 241 ! 242 call add_for_charges(ik, u(1,mode)) 243 ! 244 ! Add (dK(r)/du - dS/du) r | psi> 245 ! 246 call add_dkmds(ik, u(1,mode), jpol, dvkb) 247 ! 248 ! And calculate finally the scalar product 249 ! 250 do ibnd = 1, nbnd_occ(ikk) 251 zstareu0(jpol,mode)=zstareu0(jpol, mode) - weight * & 252 dot_product(evc(1:npwx*npol,ibnd),dvpsi(1:npwx*npol,ibnd)) 253 enddo 254 enddo 255 enddo 256 imode0 = imode0 + npert (irr) 257 enddo 258 enddo 259 deallocate (dvkb) 260 261 deallocate (pdsp) 262 deallocate (dbecsum) 263 if (noncolin) deallocate(dbecsum_nc) 264 deallocate (dvscf) 265 deallocate (aux1) 266 267 fact=1.0_DP 268#if defined(__MPI) 269 fact=1.0_DP/nproc_pool/npool 270#endif 271 IF (okpaw) THEN 272 imode0 = 0 273 do irr1 = 1, nirr 274 do ipert = 1, npert (irr1) 275 mode = imode0 + ipert 276 do nt=1,ntyp 277 if (upf(nt)%tpawp) then 278 ijh=0 279 do ih=1,nh(nt) 280 do jh=ih,nh(nt) 281 ijh=ijh+1 282 do na=1,nat 283 if (ityp(na)==nt) then 284 do jpol = 1, 3 285 do is=1,nspin_mag 286 zstareu0(jpol,mode)=zstareu0(jpol,mode) & 287 -fact*int3_paw(ih,jh,na,is,jpol)* & 288 becsumort(ijh,na,is,mode) 289 enddo 290 enddo 291 endif 292 enddo 293 enddo 294 enddo 295 endif 296 enddo 297 enddo 298 imode0 = imode0 + npert (irr1) 299 enddo 300 endif 301 302#ifdef TIMINIG_ZSTAR_US 303 call stop_clock('zstar_us_5') 304 call stop_clock('zstar_eu_us') 305#endif 306 307 return 308end subroutine zstar_eu_us 309