1! 2! Copyright (C) 2015 GWW 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 SUBROUTINE adduspos_gamma_r(iw,jw,r_ij,ik,becp_iw,becp_jw) 10 !---------------------------------------------------------------------- 11 ! 12 ! This routine adds the US term < Psi_iw|r><r|Psi_jw> 13 ! to the array r_ij 14 ! this is a GAMMA only routine (i.e. r_ij is real) 15 ! 16 USE kinds, ONLY : DP 17 USE realus, ONLY : tabp 18 USE ions_base, ONLY : nat, ntyp => nsp, ityp 19 USE fft_base, ONLY : dfftp 20 USE gvect, ONLY : ngm, gg, g 21 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk 22 USE uspp, ONLY : okvan, becsum, nkb, ijtoh, indv_ijkb0 23 USE uspp_param, ONLY : upf, lmaxq, nh 24 USE wvfct, ONLY : wg 25 USE control_flags, ONLY : gamma_only 26 USE wavefunctions, ONLY : psic 27 USE io_global, ONLY : stdout 28 USE cell_base, ONLY : omega 29 ! 30 USE mp_bands, ONLY : intra_bgrp_comm 31 USE mp, ONLY : mp_sum 32 ! 33 IMPLICIT NONE 34 ! 35 ! 36 INTEGER, INTENT(in) :: iw,jw!the states indices 37 REAL(kind=DP), INTENT(inout) :: r_ij(dfftp%nnr)!where to add the us term 38 INTEGER, INTENT(in) :: ik!spin index for spin polarized calculations NOT IMPLEMENTED YET 39 REAL(kind=DP), INTENT(in) :: becp_iw( nkb)!overlap of wfcs with us projectors 40 REAL(kind=DP), INTENT(in) :: becp_jw( nkb)!overlap of wfcs with us projectors 41 42 ! here the local variables 43 ! 44 45 INTEGER :: na, nt, nhnt, ir, ih, jh, is , ia, mbia, irb, iqs, sizeqsave 46 INTEGER :: ikb, jkb, np 47 ! counters 48 49 ! work space for rho(G,nspin) 50 ! Fourier transform of q 51 52 IF (.not.okvan) RETURN 53 IF( .not.gamma_only) THEN 54 WRITE(stdout,*) ' adduspos_gamma_r is a gamma ONLY routine' 55 STOP 56 ENDIF 57 58 DO is=1,nspin 59 ! 60 DO np = 1, ntyp 61 ! 62 iqs = 0 63 ! 64 IF ( upf(np)%tvanp ) THEN 65 ! 66 DO ia = 1, nat 67 ! 68 mbia = tabp(ia)%maxbox 69 nt = ityp(ia) 70 nhnt = nh(nt) 71 ! 72 IF ( ityp(ia) /= np ) iqs=iqs+(nhnt+1)*nhnt*mbia/2 73 IF ( ityp(ia) /= np ) CYCLE 74 ! 75 DO ih = 1, nhnt 76 ! 77 ikb = indv_ijkb0(ia) + ih 78 ! 79 DO jh = ih, nhnt 80 ! 81 jkb = indv_ijkb0(ia) + jh 82 ! 83 DO ir = 1, mbia 84 ! 85 irb = tabp(ia)%box(ir) 86 iqs = iqs + 1 87 ! 88 r_ij(irb) = r_ij(irb) + tabp(ia)%qr(ir,ijtoh(ih,jh,np))& 89 *becp_iw(ikb)*becp_jw(jkb)*omega 90 ! 91 IF ( ih /= jh ) THEN 92 r_ij(irb) = r_ij(irb) + tabp(ia)%qr(ir,ijtoh(ih,jh,np))& 93 *becp_iw(jkb)*becp_jw(ikb)*omega 94 ENDIF 95 ENDDO 96 ENDDO 97 ENDDO 98 ! 99 ENDDO 100 ! 101 ENDIF 102 ENDDO 103 ENDDO 104 ! 105 RETURN 106 ! 107 END SUBROUTINE adduspos_gamma_r 108 ! 109 SUBROUTINE adduspos_r(r_ij,becp_iw,becp_jw) 110 !---------------------------------------------------------------------- 111 ! 112 ! This routine adds the US term < Psi_iw|r><r|Psi_jw> 113 ! to the array r_ij 114 USE kinds, ONLY : DP 115 USE realus, ONLY : tabp 116 USE ions_base, ONLY : nat, ntyp => nsp, ityp 117 USE fft_base, ONLY : dfftp 118 USE gvect, ONLY : ngm, gg, g 119 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk 120 USE uspp, ONLY : okvan, becsum, nkb, ijtoh, indv_ijkb0 121 USE uspp_param, ONLY : upf, lmaxq, nh 122 USE wvfct, ONLY : wg 123 USE control_flags, ONLY : gamma_only 124 USE wavefunctions, ONLY : psic 125 USE cell_base, ONLY : omega 126 ! 127 IMPLICIT NONE 128 ! 129 COMPLEX(kind=DP), INTENT(inout) :: r_ij(dfftp%nnr)!where to add the us term 130 COMPLEX(kind=DP), INTENT(in) :: becp_iw( nkb)!overlap of wfcs with us projectors 131 COMPLEX(kind=DP), INTENT(in) :: becp_jw( nkb)!overlap of wfcs with us projectors 132 ! here the local variables 133 ! 134 INTEGER :: na, ia, nt, nhnt, ir, ih, jh, is, mbia, irb, iqs 135 INTEGER :: ikb, jkb, np 136 ! counters 137 138 ! work space for rho(G,nspin) 139 ! Fourier transform of q 140 141 IF (.not.okvan) RETURN 142 143 DO is=1,nspin 144 ! 145 DO np = 1, ntyp 146 ! 147 iqs = 0 148 ! 149 IF ( upf(np)%tvanp ) THEN 150 ! 151 DO ia = 1, nat 152 ! 153 mbia = tabp(ia)%maxbox 154 nt = ityp(ia) 155 nhnt = nh(nt) 156 ! 157 IF ( ityp(ia) /= np ) iqs=iqs+(nhnt+1)*nhnt*mbia/2 158 IF ( ityp(ia) /= np ) CYCLE 159 ! 160 DO ih = 1, nhnt 161 ! 162 ikb = indv_ijkb0(ia) + ih 163 DO jh = ih, nhnt 164 ! 165 jkb = indv_ijkb0(ia) + jh 166 ! 167 DO ir = 1, mbia 168 ! 169 irb = tabp(ia)%box(ir) 170 iqs = iqs + 1 171 ! 172 r_ij(irb) = r_ij(irb) + tabp(ia)%qr(ir,ijtoh(ih,jh,np))& 173 *conjg(becp_iw(ikb))*becp_jw(jkb)*omega 174 ! 175 IF ( ih /= jh ) THEN 176 r_ij(irb) = r_ij(irb) + tabp(ia)%qr(ir,ijtoh(ih,jh,np))& 177 *conjg(becp_iw(jkb))*becp_jw(ikb)*omega 178 ENDIF 179 ENDDO 180 ENDDO 181 ENDDO 182 ! 183 ENDDO 184 ! 185 ENDIF 186 ENDDO 187 ENDDO 188 ! 189 RETURN 190 END SUBROUTINE adduspos_r 191 ! 192 SUBROUTINE adduspos_real(sca,qq_op,becp_iw,becp_jw) 193 !---------------------------------------------------------------------- 194 ! 195 ! This routine adds the US term < Psi_iw|r><r|Psi_jw> 196 ! to the array r_ij 197 USE kinds, ONLY : DP 198 USE realus, ONLY : tabp 199 USE ions_base, ONLY : nat, ntyp => nsp, ityp 200 USE gvect, ONLY : ngm, gg, g 201 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk 202 USE uspp, ONLY : okvan, becsum, nkb, indv_ijkb0 203 USE uspp_param, ONLY : upf, lmaxq, nh, nhm 204 USE wvfct, ONLY : wg 205 USE control_flags, ONLY : gamma_only 206 USE wavefunctions, ONLY : psic 207 USE cell_base, ONLY : omega 208 ! 209 USE mp_bands, ONLY : intra_bgrp_comm 210 USE mp, ONLY : mp_sum 211 ! 212 IMPLICIT NONE 213 ! 214 REAL(kind=DP), INTENT(inout) :: sca!where to add the us term 215 REAL(kind=DP), INTENT(in) :: becp_iw( nkb)!overlap of wfcs with us projectors 216 REAL(kind=DP), INTENT(in) :: becp_jw( nkb)!overlap of wfcs with us projectors 217 REAL(kind=DP), INTENT(in) :: qq_op(nhm, nhm,nat)!US augmentation charges 218 219 ! here the local variables 220 ! 221 222 INTEGER :: na, ia, nhnt, nt, ih, jh, is, mbia 223 INTEGER :: ikb, jkb, np 224 ! counters 225 226 ! work space for rho(G,nspin) 227 ! Fourier transform of q 228 229 IF (.not.okvan) RETURN 230 231 DO is=1,nspin 232 ! 233 DO np = 1, ntyp 234 ! 235 IF ( upf(np)%tvanp ) THEN 236 ! 237 DO ia = 1, nat 238 ! 239 IF ( ityp(ia) /= np ) CYCLE 240 ! 241 mbia = tabp(ia)%maxbox 242 nt = ityp(ia) 243 nhnt = nh(nt) 244 ! 245 DO ih = 1, nhnt 246 ! 247 ikb = indv_ijkb0(ia) + ih 248 DO jh = ih, nhnt 249 ! 250 jkb = indv_ijkb0(ia) + jh 251 ! 252 sca = sca + qq_op(ih,jh,ia) * becp_iw(ikb)*becp_jw(jkb) 253 ! 254 IF ( ih /= jh ) THEN 255 sca = sca + qq_op(jh,ih,ia) * becp_iw(ikb)*becp_jw(jkb) 256 ENDIF 257 ! 258 ENDDO 259 ENDDO 260 ! 261 ENDDO 262 ! 263 ENDIF 264 ENDDO 265 ENDDO 266 ! 267 RETURN 268 ! 269 END SUBROUTINE adduspos_real 270 ! 271