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 drhodvnl (ik, ikk, nper, nu_i0, wdyn, becp1, alphap, & 11 dbecq, dalpq) 12 !----------------------------------------------------------------------- 13 ! 14 ! This subroutine computes the electronic term 2 <dpsi|dv-e ds|psi> of 15 ! the dynamical matrix. It can be used both for KB and for US 16 ! pseudopotentials. All the nonlocal (and overlap matrix) terms 17 ! are computed here. The contribution of the local potential is not 18 ! computed here. This routine must be called for each k point and 19 ! accumulates in wdyn the contribution of each k point. 20 ! 21 USE kinds, ONLY : DP 22 USE ions_base, ONLY : nat, ntyp => nsp, ityp 23 USE noncollin_module, ONLY : noncolin, npol 24 USE uspp, ONLY : okvan, nkb 25 USE uspp_param,ONLY : nh, nhm 26 USE becmod, ONLY : bec_type 27 USE wvfct, ONLY : nbnd, et 28 USE klist, ONLY : wk 29 USE lsda_mod, ONLY : current_spin, nspin 30 USE spin_orb, ONLY : lspinorb 31 USE phus, ONLY : int1, int1_nc, int2, int2_so 32 USE qpoint, ONLY : nksq 33 USE mp_bands, ONLY: intra_bgrp_comm 34 USE mp, ONLY: mp_sum 35 36 implicit none 37 integer :: ik, ikk, nper, nu_i0 38 ! input: the current k point 39 ! input: the number of perturbations 40 ! input: the initial mode 41 42 TYPE(bec_type) :: dbecq(nper), dalpq(3,nper), becp1(nksq), alphap(3,nksq) 43 ! input: the becp with psi_{k+q} 44 ! input: the alphap with psi_{k} 45 complex(DP) :: wdyn (3 * nat, 3 * nat) 46 ! output: the term of the dynamical matryx 47 48 complex(DP) :: ps, ps_nc(npol), dynwrk (3 * nat, 3 * nat) 49 ! dynamical matrix 50 complex(DP) , allocatable :: ps1 (:,:), ps2 (:,:,:) 51 complex(DP) , allocatable :: ps1_nc (:,:,:), ps2_nc (:,:,:,:), & 52 deff_nc(:,:,:,:) 53 real(DP), allocatable :: deff(:,:,:) 54 55 integer :: ibnd, ijkb0, ijkb0b, ih, jh, ikb, jkb, ipol, & 56 startb, lastb, iper, na, nb, nt, ntb, mu, nu, is, js, ijs 57 ! counters 58 59 IF (noncolin) THEN 60 allocate (ps1_nc ( nkb, npol, nbnd)) 61 allocate (ps2_nc ( nkb, npol, nbnd, 3)) 62 allocate (deff_nc ( nhm, nhm, nat, nspin )) 63 ps1_nc = (0.d0, 0.d0) 64 ps2_nc = (0.d0, 0.d0) 65 ELSE 66 allocate (ps1 ( nkb , nbnd)) 67 allocate (ps2 ( nkb , nbnd , 3)) 68 allocate (deff ( nhm, nhm, nat )) 69 ps1 = (0.d0, 0.d0) 70 ps2 = (0.d0, 0.d0) 71 END IF 72 73 dynwrk (:, :) = (0.d0, 0.d0) 74 75 call divide (intra_bgrp_comm, nbnd, startb, lastb) 76 ! 77 ! Here we prepare the two terms 78 ! 79 do ibnd = startb, lastb 80 IF (noncolin) THEN 81 CALL compute_deff_nc(deff_nc,et(ibnd,ikk)) 82 ELSE 83 CALL compute_deff(deff,et(ibnd,ikk)) 84 ENDIF 85 ijkb0 = 0 86 do nt = 1, ntyp 87 do na = 1, nat 88 if (ityp (na) == nt) then 89 do ih = 1, nh (nt) 90 ikb = ijkb0 + ih 91 do jh = 1, nh (nt) 92 jkb = ijkb0 + jh 93 IF (noncolin) THEN 94 ijs=0 95 DO is=1, npol 96 DO js=1, npol 97 ijs=ijs+1 98 ps1_nc(ikb,is,ibnd)=ps1_nc(ikb,is,ibnd) + & 99 deff_nc(ih,jh,na,ijs) * becp1(ik)%nc(jkb,js,ibnd) 100 END DO 101 END DO 102 ELSE 103 ps1 (ikb, ibnd) = ps1 (ikb, ibnd) + & 104 deff(ih,jh,na)*becp1(ik)%k(jkb,ibnd) 105 END IF 106 do ipol = 1, 3 107 IF (noncolin) THEN 108 ijs=0 109 DO is=1, npol 110 DO js=1, npol 111 ijs=ijs+1 112 ps2_nc(ikb,is,ibnd,ipol) = & 113 ps2_nc(ikb,is,ibnd,ipol)+ & 114 deff_nc(ih,jh,na,ijs) * & 115 alphap(ipol,ik)%nc(jkb,js,ibnd) 116 END DO 117 END DO 118 ELSE 119 ps2 (ikb, ibnd, ipol) = ps2 (ikb, ibnd, ipol) + & 120 deff(ih,jh,na) * alphap(ipol,ik)%k(jkb,ibnd) 121 END IF 122 123 IF (okvan) THEN 124 IF (noncolin) THEN 125 ijs=0 126 DO is=1, npol 127 DO js=1, npol 128 ijs=ijs+1 129 ps2_nc (ikb, is, ibnd, ipol) = & 130 ps2_nc (ikb, is, ibnd, ipol) + & 131 int1_nc(ih, jh, ipol, na, ijs) * & 132 becp1(ik)%nc (jkb, js, ibnd) 133 END DO 134 END DO 135 ELSE 136 ps2 (ikb, ibnd, ipol) = & 137 ps2 (ikb, ibnd, ipol) + & 138 int1 (ih, jh, ipol, na, current_spin) * & 139 becp1(ik)%k (jkb, ibnd) 140 END IF 141 END IF 142 enddo ! ipol 143 enddo ! jh 144 enddo ! ih 145 ijkb0 = ijkb0 + nh (nt) 146 endif 147 enddo ! na 148 enddo ! nt 149 enddo ! nbnd 150 ! 151 ! Here starts the loop on the atoms (rows) 152 ! 153 ijkb0 = 0 154 do nt = 1, ntyp 155 do na = 1, nat 156 if (ityp (na) == nt) then 157 do ipol = 1, 3 158 mu = 3 * (na - 1) + ipol 159 do ibnd = startb, lastb 160 do ih = 1, nh (nt) 161 ikb = ijkb0 + ih 162 do iper = 1, nper 163 nu = nu_i0 + iper 164 IF (noncolin) THEN 165 DO is=1, npol 166 dynwrk (nu, mu) = dynwrk (nu, mu) +2.d0*wk(ikk)* & 167 (ps2_nc(ikb,is,ibnd,ipol)* & 168 CONJG(dbecq(iper)%nc(ikb,is,ibnd))+ & 169 ps1_nc(ikb,is,ibnd)*CONJG( & 170 dalpq(ipol,iper)%nc(ikb,is,ibnd)) ) 171 END DO 172 ELSE 173 dynwrk (nu, mu) = dynwrk (nu, mu) + & 174 2.d0 * wk (ikk) * (ps2 (ikb, ibnd, ipol) * & 175 CONJG(dbecq(iper)%k(ikb, ibnd) ) + & 176 ps1(ikb,ibnd) * CONJG(dalpq(ipol,iper)%k(ikb,ibnd))) 177 END IF 178 enddo 179 enddo 180 if (okvan) then 181 ijkb0b = 0 182 do ntb = 1, ntyp 183 do nb = 1, nat 184 if (ityp (nb) == ntb) then 185 do ih = 1, nh (ntb) 186 ikb = ijkb0b + ih 187 IF (noncolin) THEN 188 ps_nc = (0.d0, 0.d0) 189 ELSE 190 ps = (0.d0, 0.d0) 191 END IF 192 do jh = 1, nh (ntb) 193 jkb = ijkb0b + jh 194 IF (noncolin) THEN 195 IF (lspinorb) THEN 196 ijs=0 197 DO is=1, npol 198 DO js=1, npol 199 ijs=ijs+1 200 ps_nc(is)=ps_nc(is)+ & 201 int2_so(ih,jh,ipol,na,nb,ijs)*& 202 becp1(ik)%nc(jkb, js, ibnd) 203 END DO 204 END DO 205 ELSE 206 DO is=1, npol 207 ps_nc(is)=ps_nc(is)+ & 208 int2(ih,jh,ipol,na,nb)*& 209 becp1(ik)%nc(jkb, is, ibnd) 210 END DO 211 END IF 212 ELSE 213 ps = ps + int2 (ih, jh, ipol, na, nb) * & 214 becp1(ik)%k (jkb, ibnd) 215 ENDIF 216 enddo 217 do iper = 1, nper 218 nu = nu_i0 + iper 219 IF (noncolin) THEN 220 DO is=1, npol 221 dynwrk (nu, mu) = dynwrk (nu, mu) + & 222 2.d0 * wk (ikk) * ps_nc(is) * & 223 CONJG(dbecq(iper)%nc(ikb, is, ibnd)) 224 END DO 225 ELSE 226 dynwrk (nu, mu) = dynwrk (nu, mu) + & 227 2.d0 * wk (ikk) * ps * & 228 CONJG(dbecq(iper)%k(ikb,ibnd) ) 229 END IF 230 enddo 231 enddo 232 ijkb0b = ijkb0b + nh (ntb) 233 endif 234 enddo 235 enddo 236 endif 237 enddo 238 enddo 239 ijkb0 = ijkb0 + nh (nt) 240 endif 241 enddo 242 enddo 243 call mp_sum ( dynwrk, intra_bgrp_comm ) 244 wdyn (:,:) = wdyn (:,:) + dynwrk (:,:) 245 246 IF (noncolin) THEN 247 deallocate (ps2_nc) 248 deallocate (ps1_nc) 249 deallocate (deff_nc) 250 ELSE 251 deallocate (ps2) 252 deallocate (deff) 253 END IF 254 return 255end subroutine drhodvnl 256