1! 2! Copyright (C) 2001-2008 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 dvqpsi_us_only (ik, uact, becp1, alphap) 11 !---------------------------------------------------------------------- 12 ! 13 ! This routine calculates dV_bare/dtau * psi for one perturbation 14 ! with a given q. The displacements are described by a vector uact. 15 ! The result is stored in dvpsi. The routine is called for each k point 16 ! and for each pattern u. It computes simultaneously all the bands. 17 ! This routine implements Eq. B29 of PRB 64, 235118 (2001). 18 ! Only the contribution of the nonlocal potential is calculated here. 19 ! 20 ! 21 USE kinds, only : DP 22 USE cell_base, ONLY : tpiba 23 USE gvect, ONLY : g 24 USE klist, ONLY : xk, ngk, igk_k 25 USE ions_base, ONLY : nat, ityp, ntyp => nsp 26 USE lsda_mod, ONLY : lsda, current_spin, isk, nspin 27 USE spin_orb, ONLY : lspinorb 28 USE wvfct, ONLY : nbnd, npwx, et 29 USE noncollin_module, ONLY : noncolin, npol 30 USE uspp, ONLY: okvan, nkb, vkb 31 USE uspp_param, ONLY: nh, nhm 32 USE phus, ONLY : int1, int1_nc, int2, int2_so 33 34 USE qpoint, ONLY : nksq, ikks, ikqs 35 USE becmod, ONLY : bec_type 36 USE eqv, ONLY : dvpsi 37 USE control_lr, ONLY : lgamma 38 39 implicit none 40 ! 41 ! The dummy variables 42 ! 43 44 integer :: ik 45 ! input: the k point 46 complex(DP) :: uact (3 * nat) 47 ! input: the pattern of displacements 48 TYPE(bec_type) :: becp1(nksq), alphap(3,nksq) 49 ! 50 ! And the local variables 51 ! 52 53 integer :: na, nb, mu, nu, ikk, ikq, ig, igg, nt, ibnd, ijkb0, & 54 ikb, jkb, ih, jh, ipol, is, js, ijs, npwq 55 ! counter on atoms 56 ! counter on modes 57 ! the point k 58 ! the point k+q 59 ! counter on G vectors 60 ! auxiliary counter on G vectors 61 ! counter on atomic types 62 ! counter on bands 63 ! auxiliary variable for counting 64 ! counter on becp functions 65 ! counter on becp functions 66 ! counter on n index 67 ! counter on m index 68 ! counter on polarizations 69 70 real(DP), parameter :: eps = 1.d-12 71 72 complex(DP), allocatable :: ps1 (:,:), ps2 (:,:,:), aux (:), deff_nc(:,:,:,:) 73 real(DP), allocatable :: deff(:,:,:) 74 complex(DP), allocatable :: ps1_nc (:,:,:), ps2_nc (:,:,:,:) 75 ! work space 76 77 logical :: ok 78 79 call start_clock ('dvqpsi_us_on') 80 if (noncolin) then 81 allocate (ps1_nc(nkb , npol, nbnd)) 82 allocate (ps2_nc(nkb , npol, nbnd , 3)) 83 allocate (deff_nc(nhm, nhm, nat, nspin)) 84 else 85 allocate (ps1 ( nkb , nbnd)) 86 allocate (ps2 ( nkb , nbnd , 3)) 87 allocate (deff(nhm, nhm, nat)) 88 end if 89 allocate (aux ( npwx)) 90 ikk = ikks(ik) 91 ikq = ikqs(ik) 92 if (lsda) current_spin = isk (ikk) 93 ! 94 ! we first compute the coefficients of the vectors 95 ! 96 if (noncolin) then 97 ps1_nc(:,:,:) = (0.d0, 0.d0) 98 ps2_nc(:,:,:,:) = (0.d0, 0.d0) 99 else 100 ps1(:,:) = (0.d0, 0.d0) 101 ps2(:,:,:) = (0.d0, 0.d0) 102 end if 103 do ibnd = 1, nbnd 104 IF (noncolin) THEN 105 CALL compute_deff_nc(deff_nc,et(ibnd,ikk)) 106 ELSE 107 CALL compute_deff(deff,et(ibnd,ikk)) 108 ENDIF 109 ijkb0 = 0 110 do nt = 1, ntyp 111 do na = 1, nat 112 if (ityp (na) .eq.nt) then 113 mu = 3 * (na - 1) 114 do ih = 1, nh (nt) 115 ikb = ijkb0 + ih 116 do jh = 1, nh (nt) 117 jkb = ijkb0 + jh 118 do ipol = 1, 3 119 if ( abs (uact (mu + 1) ) + & 120 abs (uact (mu + 2) ) + & 121 abs (uact (mu + 3) ) > eps) then 122 IF (noncolin) THEN 123 ijs=0 124 DO is=1,npol 125 DO js=1,npol 126 ijs=ijs+1 127 ps1_nc(ikb,is,ibnd)=ps1_nc(ikb,is,ibnd) + & 128 deff_nc(ih,jh,na,ijs) * & 129 alphap(ipol, ik)%nc(jkb,js,ibnd)* & 130 uact(mu + ipol) 131 ps2_nc(ikb,is,ibnd,ipol)= & 132 ps2_nc(ikb,is,ibnd,ipol)+ & 133 deff_nc(ih,jh,na,ijs) * & 134 becp1(ik)%nc(jkb,js,ibnd) * & 135 (0.d0,-1.d0) * uact(mu+ipol) * tpiba 136 END DO 137 END DO 138 ELSE 139 ps1 (ikb, ibnd) = ps1 (ikb, ibnd) + & 140 deff(ih, jh, na) * & 141 alphap(ipol, ik)%k(jkb, ibnd) * uact (mu + ipol) 142 ps2 (ikb, ibnd, ipol) = ps2 (ikb, ibnd, ipol) +& 143 deff(ih,jh,na)*becp1(ik)%k (jkb, ibnd) * & 144 (0.0_DP,-1.0_DP) * uact (mu + ipol) * tpiba 145 ENDIF 146 IF (okvan) THEN 147 IF (noncolin) THEN 148 ijs=0 149 DO is=1,npol 150 DO js=1,npol 151 ijs=ijs+1 152 ps1_nc(ikb,is,ibnd)=ps1_nc(ikb,is,ibnd)+ & 153 int1_nc(ih,jh,ipol,na,ijs) * & 154 becp1(ik)%nc(jkb,js,ibnd)*uact(mu+ipol) 155 END DO 156 END DO 157 ELSE 158 ps1 (ikb, ibnd) = ps1 (ikb, ibnd) + & 159 (int1 (ih, jh, ipol,na, current_spin) * & 160 becp1(ik)%k (jkb, ibnd) ) * uact (mu +ipol) 161 END IF 162 END IF 163 END IF ! uact>0 164 if (okvan) then 165 do nb = 1, nat 166 nu = 3 * (nb - 1) 167 IF (noncolin) THEN 168 IF (lspinorb) THEN 169 ijs=0 170 DO is=1,npol 171 DO js=1,npol 172 ijs=ijs+1 173 ps1_nc(ikb,is,ibnd)= & 174 ps1_nc(ikb,is,ibnd)+ & 175 int2_so(ih,jh,ipol,nb,na,ijs)* & 176 becp1(ik)%nc(jkb,js,ibnd)*uact(nu+ipol) 177 END DO 178 END DO 179 ELSE 180 DO is=1,npol 181 ps1_nc(ikb,is,ibnd)=ps1_nc(ikb,is,ibnd)+ & 182 int2(ih,jh,ipol,nb,na) * & 183 becp1(ik)%nc(jkb,is,ibnd)*uact(nu+ipol) 184 END DO 185 END IF 186 ELSE 187 ps1 (ikb, ibnd) = ps1 (ikb, ibnd) + & 188 (int2 (ih, jh, ipol, nb, na) * & 189 becp1(ik)%k (jkb, ibnd) ) * uact (nu + ipol) 190 END IF 191 enddo 192 endif ! okvan 193 enddo ! ipol 194 enddo ! jh 195 enddo ! ih 196 ijkb0 = ijkb0 + nh (nt) 197 endif 198 enddo ! na 199 enddo ! nt 200 enddo ! nbnd 201 ! 202 ! This term is proportional to beta(k+q+G) 203 ! 204 npwq = ngk(ikq) 205 if (nkb.gt.0) then 206 if (noncolin) then 207 call zgemm ('N', 'N', npwq, nbnd*npol, nkb, & 208 (1.d0, 0.d0), vkb, npwx, ps1_nc, nkb, (1.d0, 0.d0) , dvpsi, npwx) 209 else 210 call zgemm ('N', 'N', npwq, nbnd, nkb, & 211 (1.d0, 0.d0) , vkb, npwx, ps1, nkb, (1.d0, 0.d0) , dvpsi, npwx) 212 end if 213 end if 214 ! 215 ! This term is proportional to (k+q+G)_\alpha*beta(k+q+G) 216 ! 217 do ikb = 1, nkb 218 do ipol = 1, 3 219 ok = .false. 220 IF (noncolin) THEN 221 do ibnd = 1, nbnd 222 ok = ok.or.(abs (ps2_nc (ikb, 1, ibnd, ipol) ).gt.eps).or. & 223 (abs (ps2_nc (ikb, 2, ibnd, ipol) ).gt.eps) 224 end do 225 ELSE 226 do ibnd = 1, nbnd 227 ok = ok.or. (abs (ps2 (ikb, ibnd, ipol) ) .gt.eps) 228 enddo 229 ENDIF 230 if (ok) then 231 do ig = 1, npwq 232 igg = igk_k (ig,ikq) 233 aux (ig) = vkb(ig, ikb) * (xk(ipol, ikq) + g(ipol, igg) ) 234 enddo 235 do ibnd = 1, nbnd 236 IF (noncolin) THEN 237 call zaxpy(npwq,ps2_nc(ikb,1,ibnd,ipol),aux,1,dvpsi(1,ibnd),1) 238 call zaxpy(npwq,ps2_nc(ikb,2,ibnd,ipol),aux,1, & 239 dvpsi(1+npwx,ibnd),1) 240 ELSE 241 call zaxpy (npwq, ps2(ikb,ibnd,ipol), aux, 1, dvpsi(1,ibnd), 1) 242 END IF 243 enddo 244 endif 245 enddo 246 247 enddo 248 deallocate (aux) 249 IF (noncolin) THEN 250 deallocate (ps2_nc) 251 deallocate (ps1_nc) 252 deallocate (deff_nc) 253 ELSE 254 deallocate (ps2) 255 deallocate (ps1) 256 deallocate (deff) 257 END IF 258 259 call stop_clock ('dvqpsi_us_on') 260 return 261end subroutine dvqpsi_us_only 262