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