1! 2! Copyright (C) 2001-2018 Quantum ESPRESSO 3! This file is distributed under the terms 4! GNU General Public License. See the file 5! in the root directory of the present dis 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9!----------------------------------------------------------------------- 10SUBROUTINE dvqhub_barepsi_us (ik, uact) 11 !----------------------------------------------------------------------- 12 ! 13 ! DFPT+U 14 ! This routines calculates the BARE derivative of the Hubbard potential times psi. 15 ! is = current_spin 16 ! isi = opposite of the current_spin 17 ! 18 ! |Delta V_BARE_(k+q,is) psi(ibnd,k,is)> = 19 ! + \sum_(I,m1,m2) Hubbard_U(I) * [0.5\delta_(m1,m2) - ns(m1,m2,is,I)] * 20 ! { |dqsphi(imode,I,k+q,m1)><S\phi(I,k,m2)|psi(ibnd,k,is)> + 21 ! |S\phi(I,k+q,m1)><dmqsphi(imode,I,k,m2)|psi(ibnd,k,is)> } 22 ! - \sum_(I,m1,m2) Hubbard_U(I) * dnsbare(m1,m2,is,I,imode) * 23 ! |S\phi(I,k+q,m1)><S\phi(I,k,m2)|psi(ibnd,k,is)> 24 ! 25 ! Addition of the J0 terms: 26 ! 27 ! + \sum_(I,m1,m2) Hubbard_J0(I) * ns(m1,m2,isi,I) * 28 ! { |dqsphi(imode,I,k+q,m1)><S\phi(I,k,m2)|psi(ibnd,k,is)> + 29 ! |S\phi(I,k+q,m1)><dmqsphi(imode,I,k,m2)|psi(ibnd,k,is)> } 30 ! + \sum_(I,m1,m2) Hubbard_J0(I) * dnsbare(m1,m2,isi,I,imode) * 31 ! |S\phi(I,k+q,m1)><S\phi(I,k,m2)|psi(ibnd,k,is)> 32 ! 33 ! Important: in this routine vkb is a beta function at k+q, and vkb_ is beta at k. 34 ! This is done so because vkb is calculated at k+q in solve_linter (i.e. before calling 35 ! this routine), so we keep the same attribution here. 36 ! 37 ! Written by A. Floris 38 ! Modified by I. Timrov (01.10.2018) 39 ! 40 USE kinds, ONLY : DP 41 USE io_global, ONLY : stdout, ionode 42 USE io_files, ONLY : nwordwfcU 43 USE ions_base, ONLY : nat, ityp, ntyp => nsp 44 USE klist, ONLY : xk, ngk, igk_k 45 USE ldaU, ONLY : U_projection, Hubbard_l, is_hubbard, Hubbard_J0, offsetU, nwfcU 46 USE ldaU_ph, ONLY : wfcatomk, wfcatomkpq, swfcatomk, swfcatomkpq, dwfcatomkpq, & 47 sdwfcatomk, sdwfcatomkpq, dvkb, vkbkpq, dvkbkpq, & 48 proj1, proj2, dnsbare, effU 49 USE wvfct, ONLY : npwx, nbnd 50 USE uspp, ONLY : vkb,nkb 51 USE qpoint, ONLY : nksq, ikks, ikqs 52 USE control_lr, ONLY : lgamma, ofsbeta 53 USE units_lr, ONLY : iuatwfc, iuatswfc 54 USE uspp_param, ONLY : nh 55 USE lsda_mod, ONLY : lsda, current_spin, isk 56 USE wavefunctions, ONLY : evc 57 USE eqv, ONLY : dvpsi 58 USE scf, ONLY : rho 59 USE mp_bands, ONLY : intra_bgrp_comm 60 USE mp, ONLY : mp_sum 61 USE buffers, ONLY : get_buffer 62 ! 63 IMPLICIT NONE 64 ! 65 INTEGER, INTENT(IN) :: ik 66 ! the k point under consideration 67 COMPLEX(DP), INTENT(IN) :: uact(3*nat) 68 ! the pattern of displacements 69 ! 70 ! Local variables 71 ! 72 INTEGER :: i, j, k, icart, counter, na, nt, l, ih, n, mu, ig, & 73 ihubst, ihubst1, ihubst2, nah, m, m1, m2, ibnd, op_spin, & 74 ikk, ikq, npw, npwq, ibeta 75 COMPLEX(DP), ALLOCATABLE :: aux1(:), aux2(:), aux3(:), aux4(:), aux5(:), & 76 dqsphi(:,:), dmqsphi(:,:), dvqi(:,:), dvqhbar(:,:,:,:), & 77 vkb_(:,:), dwfcatom_(:) 78 COMPLEX(DP), EXTERNAL :: ZDOTC 79 ! 80 ALLOCATE (proj1(nbnd,nwfcU)) 81 ALLOCATE (proj2(nbnd,nwfcU)) 82 ALLOCATE (aux1(npwx)) 83 ALLOCATE (aux2(npwx)) 84 ALLOCATE (aux3(npwx)) 85 ALLOCATE (aux4(npwx)) 86 ALLOCATE (aux5(npwx)) 87 ALLOCATE (dqsphi(npwx,nwfcU)) 88 ALLOCATE (dmqsphi(npwx,nwfcU)) 89 ALLOCATE (dvqi(npwx,nbnd)) 90 ALLOCATE (dvqhbar(npwx,nbnd,3,nat)) 91 ALLOCATE (vkb_(npwx,nkb)) 92 ALLOCATE (dwfcatom_(npwx)) 93 ! 94 proj1 = (0.d0, 0.d0) 95 proj2 = (0.d0, 0.d0) 96 ! 97 ikk = ikks(ik) 98 ikq = ikqs(ik) 99 npw = ngk(ikk) 100 npwq= ngk(ikq) 101 ! 102 IF (lsda) THEN 103 current_spin = isk(ikk) 104 IF (current_spin==1) THEN 105 op_spin = 2 106 ELSE 107 op_spin = 1 108 ENDIF 109 ELSE 110 op_spin = 1 111 ENDIF 112 ! 113 ! Compute the beta function at k and put the result in vkb_ 114 ! 115 IF (.NOT.lgamma) THEN 116 CALL init_us_2 (npw, igk_k(1,ikk), xk(:,ikk), vkb_) 117 ELSE 118 vkb_ = vkb 119 ENDIF 120 ! 121 ! The beta function at k+q. Let us put it in the proper array vkbkpq, 122 ! because in the following of this routine the array vkb will be 123 ! used as a workspace in the routine swfc. 124 ! 125 vkbkpq = vkb 126 ! 127 ! Calculate the derivatives of beta functions 128 ! d^{icart}beta at k and k+q for all the bands and for 129 ! the 3 cartesian directions 130 ! 131 DO icart = 1, 3 132 DO na = 1, nat 133 nt = ityp(na) 134 DO ih = 1, nh(nt) 135 ! 136 ibeta = ofsbeta(na) + ih 137 ! 138 CALL dwfc (npw, igk_k(1,ikk), ikk, icart, & 139 vkb_(:,ibeta), dvkb(:,ibeta,icart)) 140 IF (.NOT.lgamma) & 141 CALL dwfc (npwq, igk_k(1,ikq), ikq, icart, & 142 vkbkpq(:,ibeta), dvkbkpq(:,ibeta,icart)) 143 ! 144 ENDDO 145 ENDDO 146 ENDDO 147 ! 148 ! Read \phi at k and k+q from file (unit iuatwfc) 149 ! 150 CALL get_buffer (wfcatomk, nwordwfcU, iuatwfc, ikk) 151 IF (.NOT.lgamma) CALL get_buffer (wfcatomkpq, nwordwfcU, iuatwfc, ikq) 152 ! 153 ! Read S*\phi at k and k+q from file (unit iuatswfc) 154 ! 155 CALL get_buffer (swfcatomk, nwordwfcU, iuatswfc, ikk) 156 IF (.NOT.lgamma) CALL get_buffer (swfcatomkpq, nwordwfcU, iuatswfc, ikq) 157 ! 158 dvqhbar = (0.d0, 0.d0) 159 ! 160 DO na = 1, nat 161 ! 162 DO icart = 1, 3 163 ! 164 dqsphi = (0.d0, 0.d0) 165 dmqsphi = (0.d0, 0.d0) 166 ! 167 DO nah = 1, nat 168 ! 169 nt = ityp(nah) 170 ! 171 ! For Hubbard_U - Hubbard_J0 172 ! 173 IF (is_hubbard(nt)) THEN 174 ! 175 DO m = 1, 2*Hubbard_l(nt)+1 176 ! 177 ihubst = offsetU(nah) + m ! I m index 178 ! 179 IF (nah==na) THEN 180 ! 181 ! Calculate | d_icart\phi_(k,I,m)) > 182 ! 183 CALL dwfc (npw, igk_k(1,ikk), ikk, icart, & 184 wfcatomk(:,ihubst), dwfcatom_) 185 ! 186 ! Apply the S matrix: | S d_^(I,icart)\phi_(k,I,m) > 187 ! 188 CALL swfc (npw, 1, vkb_, dwfcatom_, sdwfcatomk(:,ihubst)) 189 ! 190 IF (.NOT.lgamma) THEN 191 ! 192 ! Calculate |d_icart\phi_(k+q,I,m)) > 193 ! 194 CALL dwfc (npwq, igk_k(1,ikq), ikq, icart, & 195 wfcatomkpq(:,ihubst), dwfcatom_) 196 ! 197 ! Calculate | S d_^(I,icart)\phi_(k+q,I,m) > 198 ! 199 CALL swfc (npwq, 1, vkbkpq, dwfcatom_, sdwfcatomkpq(:,ihubst)) 200 ! 201 ENDIF 202 ! 203 ENDIF 204 ! 205 ! Calculate |\Delta_q(S_k \phi_(k,I,m)) > 206 ! and |\Delta_{-q}(S_{k+q} \phi_(k+q,I,m)) > 207 ! 208 CALL delta_sphi (ikk, ikq, na, icart, nah, ihubst, wfcatomk, wfcatomkpq, & 209 sdwfcatomk, sdwfcatomkpq, vkb_, vkbkpq, dvkb(:,:,icart), & 210 dvkbkpq(:,:,icart), dqsphi, dmqsphi, 1) 211 ! 212 ! Calculate: 213 ! proj1 (ihubst,ibnd) = < S_{k}\phi_(k,I,m)| psi(ibnd,k) > 214 ! proj2 (ihubst,ibnd) = < \Delta_{-q}(S_{k+q} \phi_(k+q,I,m)) | psi(ibnd,k) > 215 ! 216 DO ibnd = 1, nbnd 217 proj1(ibnd,ihubst) = ZDOTC (npw, swfcatomk(:,ihubst), 1, evc(:,ibnd), 1) 218 proj2(ibnd,ihubst) = ZDOTC (npw, dmqsphi(:,ihubst), 1, evc(:,ibnd), 1) 219 ENDDO 220 ! 221 ENDDO ! m 222 ! 223 ENDIF 224 ! 225 ENDDO ! nah 226 ! 227 CALL mp_sum (proj1, intra_bgrp_comm) 228 CALL mp_sum (proj2, intra_bgrp_comm) 229 ! 230 DO nah = 1, nat 231 ! 232 nt = ityp(nah) 233 ! 234 dvqi = (0.d0, 0.d0) 235 ! 236 IF (is_hubbard(nt)) THEN 237 ! 238 DO ibnd = 1, nbnd 239 ! 240 DO m1 = 1, 2*Hubbard_l(nt)+1 241 ! 242 ihubst1 = offsetU(nah) + m1 243 ! 244 DO ig = 1, npwq 245 ! 246 aux1(ig) = dqsphi(ig,ihubst1) * proj1(ibnd,ihubst1) 247 ! 248 aux3(ig) = swfcatomkpq(ig,ihubst1) * proj2(ibnd,ihubst1) 249 ! 250 dvqi(ig,ibnd) = dvqi(ig,ibnd) + 0.5d0 * (aux1(ig)+aux3(ig)) 251 ! 252 ENDDO 253 ! 254 DO m2 = 1, 2*Hubbard_l(nt)+1 255 ! 256 ihubst2 = offsetU(nah) + m2 257 ! 258 DO ig = 1, npwq 259 ! 260 aux2(ig) = dqsphi(ig,ihubst1) * rho%ns(m1,m2,current_spin,nah) & 261 * proj1(ibnd, ihubst2) 262 aux4(ig) = swfcatomkpq(ig,ihubst1) * rho%ns(m1,m2,current_spin,nah) & 263 * proj2(ibnd, ihubst2) 264 aux5(ig) = swfcatomkpq(ig,ihubst1) & 265 * dnsbare(m1,m2,current_spin,nah,icart,na) & 266 * proj1(ibnd, ihubst2) 267 ! 268 dvqi(ig, ibnd) = dvqi(ig, ibnd) - aux2(ig) - aux4(ig) - aux5(ig) 269 ! 270 ENDDO 271 ! 272 ENDDO ! m2 273 ! 274 ENDDO ! m1 275 ! 276 ENDDO ! ibnd 277 ! 278 ! effU = Hubbard_U - Hubbard_J0 279 ! 280 dvqi = dvqi * effU(nt) 281 ! 282 DO ig = 1, npwq 283 dvqhbar(ig,:,icart,na) = dvqhbar(ig,:,icart,na) + dvqi(ig,:) 284 ENDDO 285 ! 286 ENDIF 287 ! 288 ! Hubbard_J0 \= 0 case 289 ! 290 dvqi = (0.d0, 0.d0) 291 ! 292 IF (Hubbard_J0(nt).NE.0.d0) THEN 293 ! 294 DO ibnd = 1, nbnd 295 ! 296 DO m1 = 1, 2*Hubbard_l(nt)+1 297 ! 298 ihubst1 = offsetU(nah) + m1 299 ! 300 ! No diagonal term for J0 301 ! 302 DO m2 = 1, 2*Hubbard_l(nt)+1 303 ! 304 ihubst2 = offsetU(nah) + m2 305 ! 306 DO ig = 1, npwq 307 aux2(ig) = dqsphi(ig, ihubst1) * rho%ns(m1,m2,op_spin,nah) & 308 * proj1(ibnd, ihubst2) 309 aux4(ig) = swfcatomkpq(ig,ihubst1) * rho%ns(m1,m2,op_spin,nah) & 310 * proj2(ibnd, ihubst2) 311 aux5(ig) = swfcatomkpq(ig,ihubst1) & 312 * dnsbare (m1,m2,op_spin,nah,icart,na) & 313 * proj1 (ibnd, ihubst2) 314 ! 315 ! Note the sign change w.r.t. the case above 316 ! 317 dvqi(ig, ibnd) = dvqi(ig, ibnd) + aux2(ig) + aux4(ig) + aux5(ig) 318 ! 319 ENDDO 320 ! 321 ENDDO ! m2 322 ! 323 ENDDO ! m1 324 ! 325 ENDDO ! ibnd 326 ! 327 dvqi = dvqi * Hubbard_J0(nt) 328 ! 329 DO ig = 1, npwq 330 dvqhbar(ig,:,icart,na) = dvqhbar(ig,:,icart,na) + dvqi(ig,:) 331 ENDDO 332 ! 333 ENDIF 334 ! 335 ENDDO ! nah 336 ! 337 ENDDO ! icart 338 ! 339 ENDDO ! na 340 ! 341 ! Compute the displacement along the pattern uact (for each band). 342 ! The result is stored in aux1. 343 ! 344 DO ibnd = 1, nbnd 345 ! 346 aux1 = (0.d0, 0.d0) 347 ! 348 DO na = 1, nat 349 mu = 3 * (na - 1) 350 ! Here is the basis transformation from cartesian to pattern uact 351 DO icart = 1, 3 352 DO ig = 1, npwq 353 aux1(ig) = aux1(ig) + dvqhbar(ig,ibnd,icart,na) * uact(mu+icart) 354 ENDDO 355 ENDDO 356 ENDDO 357 ! 358 ! Add the result to dvpsi 359 ! 360 DO ig = 1, npwq 361 dvpsi(ig,ibnd) = dvpsi(ig,ibnd) + aux1(ig) 362 ENDDO 363 ! 364 ENDDO 365 ! 366 DEALLOCATE (proj1) 367 DEALLOCATE (proj2) 368 DEALLOCATE (aux1) 369 DEALLOCATE (aux2) 370 DEALLOCATE (aux3) 371 DEALLOCATE (aux4) 372 DEALLOCATE (aux5) 373 DEALLOCATE (dqsphi) 374 DEALLOCATE (dmqsphi) 375 DEALLOCATE (dvqi) 376 DEALLOCATE (dvqhbar) 377 DEALLOCATE (vkb_) 378 DEALLOCATE (dwfcatom_) 379 ! 380 RETURN 381 ! 382END SUBROUTINE dvqhub_barepsi_us 383