1! 2! Copyright (C) 2006 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! 9!---------------------------------------------------------------------- 10SUBROUTINE compute_ppsi (ppsi, ppsi_us, ik, ipol, nbnd_occ, current_spin) 11 !---------------------------------------------------------------------- 12 ! 13 ! On output: ppsi contains P_c^+ p | psi_ik > for the ipol cartesian 14 ! coordinate 15 ! ppsi_us contains the additional term required for US PP. 16 ! See J. Chem. Phys. 120, 9935 (2004) Eq. 10. 17 ! 18 ! (important: vkb and evc must have been initialized for this k-point) 19 ! 20 USE kinds, ONLY : DP 21 USE ions_base, ONLY : nat, ityp, ntyp => nsp 22 USE cell_base, ONLY : tpiba 23 USE io_global, ONLY : stdout 24 USE wavefunctions, ONLY : evc 25 USE wvfct, ONLY : et, nbnd, npwx 26 USE uspp, ONLY : nkb, vkb, deeq, qq_nt, qq_so, deeq_nc, okvan 27 USE spin_orb, ONLY : lspinorb 28 USE lsda_mod, ONLY : nspin 29 USE gvect, ONLY : g 30 USE klist, ONLY : xk, nks, ngk, igk_k 31 USE noncollin_module, ONLY : noncolin, npol 32 USE becmod, ONLY : bec_type, becp, calbec 33 USE uspp_param, ONLY : nh, nhm 34 IMPLICIT NONE 35 ! 36 INTEGER, INTENT(in) :: ipol, ik, nbnd_occ, current_spin 37 ! 38 COMPLEX(DP) :: ppsi(npwx,npol,nbnd_occ), ppsi_us(npwx,npol,nbnd_occ) 39 ! Local variables 40 ! 41 INTEGER :: npw, ig, na, ibnd, ikb, jkb, nt, ih, jh, ip, ijkb0 42 ! counters 43 44 REAL(DP), ALLOCATABLE :: gk (:,:) 45 ! the derivative of |k+G| 46 REAL(DP) :: vers(3), gk2 47 48 COMPLEX(DP), ALLOCATABLE :: ps2(:,:,:), dvkb (:,:), dvkb1 (:,:), & 49 work (:,:), becp2(:,:), becp2_nc(:,:,:), psc(:,:,:,:), ps(:), & 50 ps_nc(:,:), dpqq_so(:,:,:,:,:) 51 52 REAL(DP), ALLOCATABLE :: dpqq(:,:,:,:) 53 54 ALLOCATE (work ( npwx, max(nkb,1))) 55 ALLOCATE (gk ( 3, npwx)) 56 IF (nkb > 0) THEN 57 IF (noncolin) THEN 58 ALLOCATE (becp2_nc (nkb, npol, nbnd)) 59 ELSE 60 ALLOCATE (becp2 (nkb, nbnd)) 61 ENDIF 62 63 ALLOCATE (dvkb (npwx, nkb)) 64 ALLOCATE (dvkb1(npwx, nkb)) 65 dvkb (:,:) = (0.d0, 0.d0) 66 dvkb1(:,:) = (0.d0, 0.d0) 67 ENDIF 68 npw = ngk(ik) 69 DO ig = 1, npw 70 gk (1, ig) = (xk (1, ik) + g (1, igk_k(ig,ik) ) ) * tpiba 71 gk (2, ig) = (xk (2, ik) + g (2, igk_k(ig,ik) ) ) * tpiba 72 gk (3, ig) = (xk (3, ik) + g (3, igk_k(ig,ik) ) ) * tpiba 73 ENDDO 74 ! 75 ! this is the kinetic contribution to p : (k+G)_ipol * psi 76 ! 77 DO ip=1,npol 78 DO ibnd = 1, nbnd_occ 79 DO ig = 1, npw 80 ppsi(ig,ip,ibnd)=gk(ipol,ig)*evc(ig+npwx*(ip-1),ibnd) 81 ENDDO 82 ENDDO 83 ENDDO 84 ! 85 ! from now on we need (k+G)_ipol / |k+G| 86 ! 87 DO ig = 1, npw 88 gk2 = gk (1, ig) **2 + gk (2, ig) **2 + gk (3, ig) **2 89 IF (gk2 < 1.0d-10) THEN 90 gk (:, ig) = 0.d0 91 ELSE 92 gk (:, ig) = gk (:, ig) / sqrt (gk2 ) 93 ENDIF 94 ENDDO 95 96 ! 97 ! and this is the contribution from nonlocal pseudopotentials 98 ! 99 CALL gen_us_dj (ik, dvkb) 100 vers=0.d0 101 vers(ipol)=1.d0 102 CALL gen_us_dy (ik, vers, dvkb1) 103 104 jkb = 0 105 DO nt = 1, ntyp 106 DO na = 1, nat 107 IF (nt == ityp (na)) THEN 108 DO ikb = 1, nh (nt) 109 jkb = jkb + 1 110 DO ig = 1, npw 111 work (ig,jkb)=dvkb1(ig,jkb)+dvkb(ig,jkb)*gk(ipol,ig) 112 ENDDO 113 ENDDO 114 ENDIF 115 ENDDO 116 ENDDO 117 DEALLOCATE (gk) 118 119 IF (noncolin) THEN 120 CALL calbec ( npw, work, evc, becp2_nc ) 121 ELSE 122 CALL calbec ( npw, work, evc, becp2 ) 123 ENDIF 124 125 ijkb0 = 0 126 IF (noncolin) THEN 127 ALLOCATE (psc( nkb, 2, nbnd_occ, 2)) 128 psc=(0.d0,0.d0) 129 ELSE 130 ALLOCATE (ps2( nkb, nbnd_occ, 2)) 131 ps2=(0.d0,0.d0) 132 ENDIF 133 DO nt = 1, ntyp 134 DO na = 1, nat 135 IF (nt == ityp (na)) THEN 136 DO ih = 1, nh (nt) 137 ikb = ijkb0 + ih 138 DO jh = 1, nh (nt) 139 jkb = ijkb0 + jh 140 DO ibnd = 1, nbnd_occ 141 IF (noncolin) THEN 142 IF (lspinorb) THEN 143 psc(ikb,1,ibnd,1)=psc(ikb,1,ibnd,1)+(0.d0,-1.d0)* & 144 (becp2_nc(jkb,1,ibnd)*(deeq_nc(ih,jh,na,1) & 145 -et(ibnd,ik)*qq_so(ih,jh,1,nt) )+ & 146 becp2_nc(jkb,2,ibnd)*(deeq_nc(ih,jh,na,2)- & 147 et(ibnd,ik)* qq_so(ih,jh,2,nt) ) ) 148 psc(ikb,2,ibnd,1)=psc(ikb,2,ibnd,1)+(0.d0,-1.d0)* & 149 (becp2_nc(jkb,1,ibnd)*(deeq_nc(ih,jh,na,3) & 150 -et(ibnd,ik)*qq_so(ih,jh,3,nt) )+ & 151 becp2_nc(jkb,2,ibnd)*(deeq_nc(ih,jh,na,4)- & 152 et(ibnd,ik)* qq_so(ih,jh,4,nt) ) ) 153 psc(ikb,1,ibnd,2)=psc(ikb,1,ibnd,2)+(0.d0,-1.d0)* & 154 (becp%nc(jkb,1,ibnd)*(deeq_nc(ih,jh,na,1) & 155 -et(ibnd,ik)*qq_so(ih,jh,1,nt) )+ & 156 becp%nc(jkb,2,ibnd)*(deeq_nc(ih,jh,na,2)- & 157 et(ibnd,ik)* qq_so(ih,jh,2,nt) ) ) 158 psc(ikb,2,ibnd,2)=psc(ikb,2,ibnd,2)+(0.d0,-1.d0)* & 159 (becp%nc(jkb,1,ibnd)*(deeq_nc(ih,jh,na,3) & 160 -et(ibnd,ik)*qq_so(ih,jh,3,nt) )+ & 161 becp%nc(jkb,2,ibnd)*(deeq_nc(ih,jh,na,4)- & 162 et(ibnd,ik)* qq_so(ih,jh,4,nt) ) ) 163 ELSE 164 psc(ikb,1,ibnd,1)=psc(ikb,1,ibnd,1)+ (0.d0,-1.d0)* & 165 ( becp2_nc(jkb,1,ibnd)*(deeq_nc(ih,jh,na,1) & 166 -et(ibnd,ik)*qq_nt(ih,jh,nt)) + & 167 becp2_nc(jkb,2,ibnd)*deeq_nc(ih,jh,na,2) ) 168 psc(ikb,2,ibnd,1)=psc(ikb,2,ibnd,1)+ (0.d0,-1.d0)* & 169 ( becp2_nc(jkb,2,ibnd)*(deeq_nc(ih,jh,na,4) & 170 -et(ibnd,ik)*qq_nt(ih,jh,nt))+ & 171 becp2_nc(jkb,1,ibnd)*deeq_nc(ih,jh,na,3) ) 172 psc(ikb,1,ibnd,2)=psc(ikb,1,ibnd,2)+ (0.d0,-1.d0)* & 173 ( becp%nc(jkb,1,ibnd)*(deeq_nc(ih,jh,na,1) & 174 -et(ibnd,ik)*qq_nt(ih,jh,nt))+ & 175 becp%nc(jkb,2,ibnd)*deeq_nc(ih,jh,na,2) ) 176 psc(ikb,2,ibnd,2)=psc(ikb,2,ibnd,2)+ (0.d0,-1.d0)* & 177 ( becp%nc(jkb,2,ibnd)*(deeq_nc(ih,jh,na,4) & 178 -et(ibnd,ik)*qq_nt(ih,jh,nt))+ & 179 becp%nc(jkb,1,ibnd)*deeq_nc(ih,jh,na,3) ) 180 ENDIF 181 ELSE 182 ps2(ikb,ibnd,1) = ps2(ikb,ibnd,1)+ becp2(jkb,ibnd)* & 183 (0.d0,-1.d0)*(deeq(ih,jh,na,current_spin) & 184 -et(ibnd,ik)*qq_nt(ih,jh,nt)) 185 ps2(ikb,ibnd,2) = ps2(ikb,ibnd,2) +becp%k(jkb,ibnd) * & 186 (0.d0,-1.d0)*(deeq(ih,jh,na,current_spin)& 187 -et(ibnd,ik)*qq_nt(ih,jh,nt)) 188 ENDIF 189 ENDDO 190 ENDDO 191 ENDDO 192 ijkb0=ijkb0+nh(nt) 193 ENDIF 194 ENDDO 195 ENDDO 196 IF (ikb /= nkb .or. jkb /= nkb) CALL errore ('compute_ppsi', & 197 'unexpected error',1) 198 199 IF (nkb>0) THEN 200 IF (noncolin) THEN 201 CALL zgemm( 'N', 'N', npwx, nbnd_occ*npol, nkb, & 202 (0.d0,0.5d0), vkb, npwx, psc(1,1,1,1), nkb, (1.d0,0.d0), & 203 ppsi, npwx ) 204 CALL zgemm( 'N', 'N', npwx, nbnd_occ*npol, nkb, & 205 (0.d0,0.5d0), work, npwx, psc(1,1,1,2), nkb, (1.d0,0.d0), & 206 ppsi, npwx ) 207 ELSE 208 CALL zgemm( 'N', 'N', npw, nbnd_occ, nkb, & 209 (0.d0,0.5d0), vkb(1,1), npwx, ps2(1,1,1), nkb, (1.d0,0.0d0), & 210 ppsi, npwx ) 211 CALL zgemm( 'N', 'N', npw, nbnd_occ, nkb, & 212 (0.d0,0.5d0), work(1,1), npwx, ps2(1,1,2), nkb, (1.d0,0.0d0), & 213 ppsi, npwx ) 214 ENDIF 215 ENDIF 216 IF (noncolin) THEN 217 DEALLOCATE (psc) 218 ELSE 219 DEALLOCATE (ps2) 220 ENDIF 221! 222! ppsi contains p - i/2 [x, V_{nl}-eS] psi_v for the ipol polarization 223! 224! In the US case there is another term in the matrix element. 225! This term has to be multiplied by the difference of the eigenvalues, 226! so it is calculated separately here and multiplied in the calling 227! routine. 228 229 IF (okvan) THEN 230 ppsi_us=(0.d0,0.d0) 231 ALLOCATE (dpqq( nhm, nhm, 3, ntyp)) 232 CALL compute_qdipol(dpqq,ipol) 233 IF (noncolin) THEN 234 ALLOCATE (ps_nc(nbnd_occ,npol)) 235 IF (lspinorb) THEN 236 ALLOCATE (dpqq_so( nhm, nhm, nspin, 3, ntyp)) 237 CALL compute_qdipol_so(dpqq, dpqq_so,ipol) 238 ENDIF 239 ELSE 240 ALLOCATE (ps(nbnd_occ)) 241 ENDIF 242 ijkb0 = 0 243 DO nt = 1, ntyp 244 DO na = 1, nat 245 IF (ityp(na)==nt) THEN 246 DO ih = 1, nh (nt) 247 ikb = ijkb0 + ih 248 IF (noncolin) THEN 249 ps_nc = (0.d0,0.d0) 250 ELSE 251 ps = (0.d0,0.d0) 252 ENDIF 253 DO jh = 1, nh (nt) 254 jkb = ijkb0 + jh 255 DO ibnd=1, nbnd_occ 256 IF (noncolin) THEN 257 DO ip=1,npol 258 IF (lspinorb) THEN 259 ps_nc(ibnd,ip)=ps_nc(ibnd,ip) + & 260 (0.d0,1.d0)*(becp2_nc(jkb,1,ibnd)* & 261 qq_so(ih,jh,1+(ip-1)*2,nt) + & 262 becp2_nc(jkb,2,ibnd) * & 263 qq_so(ih,jh,2+(ip-1)*2,nt) ) & 264 + becp%nc(jkb,1,ibnd)* & 265 dpqq_so(ih,jh,1+(ip-1)*2,ipol,nt) & 266 + becp%nc(jkb,2,ibnd)* & 267 dpqq_so(ih,jh,2+(ip-1)*2,ipol,nt) 268 ELSE 269 ps_nc(ibnd,ip)=ps_nc(ibnd,ip)+ & 270 becp2_nc(jkb,ip,ibnd)*(0.d0,1.d0)* & 271 qq_nt(ih,jh,nt)+becp%nc(jkb,ip,ibnd) & 272 *dpqq(ih,jh,ipol,nt) 273 ENDIF 274 ENDDO 275 ELSE 276 ps(ibnd) = ps(ibnd) + becp2(jkb,ibnd) * & 277 (0.d0,1.d0) * qq_nt(ih,jh,nt) + & 278 becp%k(jkb,ibnd) * dpqq(ih,jh,ipol,nt) 279 ENDIF 280 ENDDO 281 ENDDO 282 DO ibnd = 1, nbnd_occ 283 IF (noncolin) THEN 284 DO ip=1,npol 285 CALL zaxpy(npw,ps_nc(ibnd,ip),vkb(1,ikb),1,& 286 ppsi_us(1,ip,ibnd),1) 287 ENDDO 288 ELSE 289 CALL zaxpy(npw,ps(ibnd),vkb(1,ikb),1,ppsi_us(1,1,ibnd),1) 290 ENDIF 291 ENDDO 292 ENDDO 293 ijkb0=ijkb0+nh(nt) 294 ENDIF 295 ENDDO 296 ENDDO 297 IF (jkb/=nkb) CALL errore ('compute_ppsi', 'unexpected error', 1) 298 IF (noncolin) THEN 299 DEALLOCATE(ps_nc) 300 ELSE 301 DEALLOCATE(ps) 302 ENDIF 303 ENDIF 304 305 306 IF (nkb > 0) THEN 307 DEALLOCATE (dvkb1, dvkb) 308 IF (noncolin) THEN 309 DEALLOCATE(becp2_nc) 310 ELSE 311 DEALLOCATE(becp2) 312 ENDIF 313 ENDIF 314 DEALLOCATE (work) 315 316 RETURN 317END SUBROUTINE compute_ppsi 318