TABLE OF CONTENTS
ABINIT/pawaccrhoij [ Functions ]
NAME
pawaccrhoij
FUNCTION
Accumulate the PAW quantities rhoij (augmentation occupancies) or their 1-st order change or their gradient vs r Add the contribution of a given k-point and band Remember: for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
atindx(natom)=index table for atoms (sorted-->random), inverse of atindx.
cplex: if 1, WFs (or 1st-order WFs) are REAL, if 2, COMPLEX
cwaveprj(natom,nspinor)= LEFT wave function at given n,k
projected with non-local projectors: cwaveprj=<p_i|Cnk>
cwaveprj1(natom,nspinor)= RIGHT wave function at n,k,q
projected with non-local projectors: cwaveprj1=<p_i|C1nk,q>
* USED for RF : C1nk is the first-order wave function
* USED for DMFT: C1nk is the RIGHT wave function
* NOT USED in usual GS case; can be set to cwaveprj in that case
ipert=index of perturbation (RF only, i.e. option=2)
isppol=index of current spin component
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
comm_atom=--optional-- MPI communicator over atoms
my_natom=number of atoms treated by current processor
natom=number of atoms in cell
nspinor=number of spinorial components (on current proc)
occ_k=occupation number for current band n,k
option: choice of calculation:
1: update rhoij (Ground-State)
2: update 1-st order rhoij (Response Function) according to ipert
3: update gradients of rhoij with respect to r,strain of both
usetimerev=.TRUE. if time-reversal symmetry is used (WF(-k)=Conjg[WF(k)])
wtk_k=weight assigned to current k-point
SIDE EFFECTS
pawrhoij(natom) <type(pawrhoij_type)>= GS: paw rhoij occupancies and related data RF: 1-st order paw rhoij occupancies and related data On output, has been updated with the contribution of current n,k === option=1: pawrhoij(:)%rhoij_(lmn2_size,nspden)= (non symetrized) Sum_{n,k} {occ(n,k)*conjugate[cprj_nk(ii)].cprj_nk(jj)} === option=2: pawrhoij(:)%rhoij_(lmn2_size,nspden)= (non symetrized) Sum_{n,k} {occ(n,k)*(conjugate[cprj_nk(ii)].cprj1_nk,q(jj) conjugate[cprj_nk(jj)].cprj1_nk,q(ii)} + Sum_{n,k} {occ(n,k)*(conjugate[dcprj_nk(ii)/dlambda].cprj_nk(jj) +conjugate[cprj_nk(ii)].dcprj_nk(jj)/dlambda)} === option=3: pawrhoij(:)%grhoij(lmn2_size,mu,nspden)= (non symetrized) Sum_{n,k} {occ(n,k)*(conjugate[dcprj_nk(ii)/dr_mu].cprj_nk(jj) +conjugate[cprj_nk(ii)].dcprj_nk(jj)/dr_mu)}
PARENTS
d2frnl,dfpt_accrho,energy,pawmkrhoij,posdoppler,wfd_pawrhoij
CHILDREN
free_my_atmtab,get_my_atmtab
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 subroutine pawaccrhoij(atindx,cplex,cwaveprj,cwaveprj1,ipert,isppol,my_natom,natom,& 78 & nspinor,occ_k,option,pawrhoij,usetimerev,wtk_k,occ_k_2, & 79 & comm_atom,mpi_atmtab ) ! optional (parallelism) 80 81 82 use defs_basis 83 use m_profiling_abi 84 use m_errors 85 use m_xmpi, only : xmpi_comm_self 86 87 use m_pawrhoij, only : pawrhoij_type 88 use m_pawcprj, only : pawcprj_type 89 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 90 91 !This section has been created automatically by the script Abilint (TD). 92 !Do not modify the following lines by hand. 93 #undef ABI_FUNC 94 #define ABI_FUNC 'pawaccrhoij' 95 !End of the abilint section 96 97 implicit none 98 99 !Arguments --------------------------------------------- 100 !scalars 101 integer,intent(in) :: cplex,ipert,isppol,my_natom,natom,nspinor,option 102 integer,optional,intent(in) :: comm_atom 103 logical,intent(in) :: usetimerev 104 real(dp),intent(in) :: occ_k,wtk_k 105 real(dp),optional,intent(in) :: occ_k_2 106 !arrays 107 integer,intent(in) :: atindx(natom) 108 integer,optional,target,intent(in) :: mpi_atmtab(:) 109 type(pawcprj_type),intent(in) :: cwaveprj(natom,nspinor),cwaveprj1(natom,nspinor) 110 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 111 112 !Local variables --------------------------------------- 113 !scalars 114 integer :: cplex_rhoij,iatm,iatom,iatom1,ilmn,iplex,j0lmn,jlmn,klmn,klmn_im,klmn_re 115 integer :: mu,my_comm_atom,ncpgr,nspden_rhoij 116 logical :: compute_impart,compute_impart_cplex,substract_diagonal 117 logical :: my_atmtab_allocated,paral_atom 118 real(dp) :: ro11_im,ro11_re,ro12_im,ro12_re,ro21_im,ro21_re,ro22_im,ro22_re,weight,weight_2 119 character(len=500) :: message 120 !arrays 121 integer,pointer :: my_atmtab(:) 122 real(dp) :: cpi0(2,nspinor),cpi1(2,nspinor),cpj0(2,nspinor),cpj1(2,nspinor) 123 real(dp) :: dcpi0(2,nspinor,9),dcpj0(2,nspinor,9) 124 125 ! *********************************************************************** 126 127 DBG_ENTER("COLL") 128 129 if (my_natom==0) return 130 131 ncpgr=0 132 if (option==2.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4)) ncpgr=1 133 if (option==3) ncpgr=cwaveprj(1,1)%ncpgr 134 !Tests 135 if(option==2.and.(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11)) then 136 message = ' not relevant for ipert=natom+1 or ipert=natom+10 or ipert=natom+11 !' 137 MSG_BUG(message) 138 end if 139 if(option==2.and.cwaveprj(1,1)%ncpgr<ncpgr) then 140 message = ' Error on cwaveprj1 factors derivatives !' 141 MSG_BUG(message) 142 end if 143 if(option==3.and.cwaveprj(1,1)%ncpgr/=ncpgr) then 144 message = ' Error on cwaveprj factors derivatives !' 145 MSG_BUG(message) 146 end if 147 148 !Set up parallelism over atoms 149 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 150 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 151 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 152 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 153 & my_natom_ref=my_natom) 154 155 weight=wtk_k*occ_k 156 weight_2=zero 157 if(present(occ_k_2).and.nspinor==2) weight_2=wtk_k*occ_k_2 158 if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1) weight=half*weight 159 if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1.and.present(occ_k_2)) weight_2=half*weight_2 160 if (option==1) then 161 162 ! ================================================================== 163 ! === OPTION 1: Accumulate (n,k) contribution to rhoij ============= 164 ! ================================================================== 165 compute_impart=((.not.usetimerev).and.(pawrhoij(1)%cplex==2)) 166 compute_impart_cplex=((compute_impart).and.(cplex==2)) 167 if (nspinor==1) then 168 cplex_rhoij=pawrhoij(1)%cplex 169 if (cplex_rhoij==1) then 170 do iatom=1,my_natom 171 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 172 iatm=atindx(iatom1) 173 do jlmn=1,pawrhoij(iatom)%lmn_size 174 j0lmn=jlmn*(jlmn-1)/2 175 cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn) 176 do ilmn=1,jlmn 177 klmn=j0lmn+ilmn 178 cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn) 179 ro11_re=zero 180 do iplex=1,cplex 181 ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1) 182 end do 183 pawrhoij(iatom)%rhoij_(klmn,isppol)=pawrhoij(iatom)%rhoij_(klmn,isppol)+weight*ro11_re 184 end do 185 end do 186 end do 187 else 188 do iatom=1,my_natom 189 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 190 iatm=atindx(iatom1) 191 do jlmn=1,pawrhoij(iatom)%lmn_size 192 j0lmn=jlmn*(jlmn-1)/2 193 cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn) 194 do ilmn=1,jlmn 195 klmn=j0lmn+ilmn 196 klmn_re=cplex_rhoij*(klmn-1)+1 197 cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn) 198 ro11_re=zero 199 do iplex=1,cplex 200 ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1) 201 end do 202 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 203 if (compute_impart_cplex) then 204 klmn_im=klmn_re+1 205 ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1) 206 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 207 end if 208 end do 209 end do 210 end do 211 end if 212 else ! nspinor=2 213 do iatom=1,my_natom 214 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 215 iatm=atindx(iatom1) 216 cplex_rhoij=pawrhoij(iatom)%cplex 217 nspden_rhoij=pawrhoij(iatom)%nspden 218 do jlmn=1,pawrhoij(iatom)%lmn_size 219 j0lmn=jlmn*(jlmn-1)/2 220 cpj0(1:cplex,1)=cwaveprj(iatm,1)%cp(1:cplex,jlmn) 221 cpj0(1:cplex,2)=cwaveprj(iatm,2)%cp(1:cplex,jlmn) 222 do ilmn=1,jlmn 223 klmn=j0lmn+ilmn 224 klmn_re=cplex_rhoij*(klmn-1)+1 225 cpi0(1:cplex,1)=cwaveprj1(iatm,1)%cp(1:cplex,ilmn) 226 cpi0(1:cplex,2)=cwaveprj1(iatm,2)%cp(1:cplex,ilmn) 227 ro11_re=zero;ro22_re=zero 228 ro12_re=zero;ro21_re=zero 229 ro12_im=zero;ro21_im=zero 230 do iplex=1,cplex 231 ro11_re=ro11_re+cpi0(iplex,1)*cpj0(iplex,1) 232 ro22_re=ro22_re+cpi0(iplex,2)*cpj0(iplex,2) 233 end do 234 pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re) 235 if (nspden_rhoij>1) then 236 do iplex=1,cplex 237 ro12_re=ro12_re+cpi0(iplex,2)*cpj0(iplex,1) 238 ro21_re=ro21_re+cpi0(iplex,1)*cpj0(iplex,2) 239 end do 240 pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re) 241 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re) 242 if (cplex==2) then 243 ro12_im=cpi0(1,2)*cpj0(2,1)-cpi0(2,2)*cpj0(1,1) 244 ro21_im=cpi0(1,1)*cpj0(2,2)-cpi0(2,1)*cpj0(1,2) 245 pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im) 246 end if 247 end if 248 if (compute_impart) then 249 klmn_im=klmn_re+1 250 if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re) 251 if (cplex==2) then 252 ro11_im=cpi0(1,1)*cpj0(2,1)-cpi0(2,1)*cpj0(1,1) 253 ro22_im=cpi0(1,2)*cpj0(2,2)-cpi0(2,2)*cpj0(1,2) 254 pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im) 255 if (nspden_rhoij>1) then 256 pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im) 257 pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im) 258 end if 259 if (present(occ_k_2).and.nspinor==2) then 260 pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight_2*(-ro11_im-ro22_im) 261 pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight_2*( ro11_re+ro22_re) 262 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight_2*(-ro21_im-ro12_im) 263 pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight_2*( ro21_re+ro12_re) 264 pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight_2*(-ro12_re+ro21_re) 265 pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight_2*(-ro12_im+ro21_im) 266 pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight_2*(-ro11_im+ro22_im) 267 pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight_2*( ro11_re-ro22_re) 268 end if 269 end if 270 end if 271 end do 272 end do 273 end do 274 end if 275 276 else if (option==2) then 277 278 ! ================================================================== 279 ! === OPTION 2: Accumulate (n,k) contribution to 1st-order rhoij === 280 ! ================================================================== 281 282 compute_impart=(pawrhoij(1)%cplex==2) 283 compute_impart_cplex=((pawrhoij(1)%cplex==2).and.(cplex==2)) 284 substract_diagonal=(ipert==natom+3) 285 286 ! Accumulate (n,k) contribution to rhoij1 287 ! due to derivative of wave-function 288 if (nspinor==1) then 289 do iatom=1,my_natom 290 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 291 iatm=atindx(iatom1) 292 cplex_rhoij=pawrhoij(iatom)%cplex 293 do jlmn=1,pawrhoij(iatom)%lmn_size 294 j0lmn=jlmn*(jlmn-1)/2 295 cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn) 296 cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn) 297 do ilmn=1,jlmn 298 klmn=j0lmn+ilmn 299 klmn_re=cplex_rhoij*(klmn-1)+1 300 cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn) 301 cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn) 302 ro11_re=zero 303 do iplex=1,cplex 304 ro11_re=ro11_re+cpi0(iplex,1)*cpj1(iplex,1)+cpj0(iplex,1)*cpi1(iplex,1) 305 end do 306 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 307 if (compute_impart_cplex) then 308 klmn_im=klmn_re+1 309 ro11_im=cpi0(1,1)*cpj1(2,1)-cpi0(2,1)*cpj1(1,1)+cpj0(1,1)*cpi1(2,1)-cpj0(2,1)*cpi1(1,1) 310 pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 311 end if 312 end do 313 end do 314 end do 315 else ! nspinor=2 316 do iatom=1,my_natom 317 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 318 iatm=atindx(iatom1) 319 cplex_rhoij=pawrhoij(iatom)%cplex 320 nspden_rhoij=pawrhoij(iatom)%nspden 321 do jlmn=1,pawrhoij(iatom)%lmn_size 322 j0lmn=jlmn*(jlmn-1)/2 323 cpj0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,jlmn) 324 cpj0(1:2,2)=cwaveprj (iatm,2)%cp(1:2,jlmn) 325 cpj1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,jlmn) 326 cpj1(1:2,2)=cwaveprj1(iatm,2)%cp(1:2,jlmn) 327 do ilmn=1,jlmn 328 klmn=j0lmn+ilmn 329 klmn_re=cplex_rhoij*(klmn-1)+1 330 cpi0(1:2,1)=cwaveprj (iatm,1)%cp(1:2,ilmn) 331 cpi0(1:2,2)=cwaveprj (iatm,2)%cp(1:2,ilmn) 332 cpi1(1:2,1)=cwaveprj1(iatm,1)%cp(1:2,ilmn) 333 cpi1(1:2,2)=cwaveprj1(iatm,2)%cp(1:2,ilmn) 334 ro11_re=zero;ro22_re=zero 335 ro12_re=zero;ro21_re=zero 336 ro12_im=zero;ro21_im=zero 337 do iplex=1,cplex 338 ro11_re=cpj0(iplex,1)*cpi1(iplex,1)+cpi0(iplex,1)*cpj1(iplex,1) 339 ro22_re=cpj0(iplex,2)*cpi1(iplex,2)+cpi0(iplex,2)*cpj1(iplex,2) 340 end do 341 pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re) 342 if (nspden_rhoij>1) then 343 do iplex=1,cplex 344 ro12_re=cpj0(iplex,1)*cpi1(iplex,2)+cpi0(iplex,2)*cpj1(iplex,1) 345 ro21_re=cpj0(iplex,2)*cpi1(iplex,1)+cpi0(iplex,1)*cpj1(iplex,2) 346 end do 347 pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re) 348 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re) 349 if (cplex==2) then 350 ro12_im=cpj0(1,1)*cpi1(2,2)-cpi1(1,1)*cpj0(2,2)+cpi0(1,2)*cpj1(2,1)-cpj1(1,2)*cpi0(2,1) 351 ro21_im=cpj0(1,2)*cpi1(2,1)-cpi1(1,2)*cpj0(2,1)+cpi0(1,1)*cpj1(2,2)-cpj1(1,1)*cpi0(2,2) 352 pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im) 353 end if 354 end if 355 if (compute_impart) then 356 klmn_im=klmn_re+1 357 if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro12_re-ro21_re) 358 if (cplex==2) then 359 ro11_im=cpj0(1,1)*cpi1(2,1)-cpi1(1,1)*cpj0(2,1)+cpi0(1,1)*cpj1(2,1)-cpj1(1,1)*cpi0(2,1) 360 ro22_im=cpj0(1,2)*cpi1(2,2)-cpi1(1,2)*cpj0(2,2)+cpi0(1,2)*cpj1(2,2)-cpj1(1,2)*cpi0(2,2) 361 pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im) 362 if (nspden_rhoij>1) then 363 pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im) 364 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_im+ro21_im) 365 end if 366 end if 367 end if 368 end do 369 end do 370 end do 371 end if 372 373 ! Accumulate (n,k) contribution to rhoij1 374 ! due to derivative of projectors 375 if (ipert/=natom+2) then 376 if (nspinor==1) then 377 do iatom=1,my_natom 378 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 379 iatm=atindx(iatom1) 380 if (ipert<=natom.and.iatom/=ipert) cycle 381 cplex_rhoij=pawrhoij(iatom)%cplex 382 do jlmn=1,pawrhoij(iatom)%lmn_size 383 j0lmn=jlmn*(jlmn-1)/2 384 cpj0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,jlmn) 385 dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn) 386 do ilmn=1,jlmn 387 klmn=j0lmn+ilmn 388 klmn_re=cplex_rhoij*(klmn-1)+1 389 cpi0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,ilmn) 390 dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn) 391 ro11_re=zero 392 do iplex=1,cplex 393 ro11_re=ro11_re+dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1) 394 end do 395 if (substract_diagonal) then 396 do iplex=1,cplex 397 ro11_re=ro11_re-cpi0(iplex,1)*cpj0(iplex,1) 398 end do 399 end if 400 pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re 401 ! This imaginary part does not have to be computed 402 ! It is cancelled because rho_ij+rho_ji is stored in rho_ij 403 ! if (compute_impart_cplex) then 404 ! klmn_im=klmn_re+1 405 ! ro11_im=dcpi0(1,1,1)*cpj0(2,1)-dcpi0(2,1,1)*cpj0(1,1)+cpi0(1,1)*dcpj0(2,1,1)-cpi0(2,1)*dcpj0(1,1,1) 406 ! if (substract_diagonal) then 407 ! ro11_im=ro11_im-cpi0(1,1)*cpj0(2,1)+cpi0(2,1)*cpj0(1,1) 408 ! end if 409 ! pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im 410 ! end if 411 end do 412 end do 413 end do 414 else ! nspinor=2 415 do iatom=1,my_natom 416 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 417 iatm=atindx(iatom1) 418 if (ipert<=natom.and.iatom/=ipert) cycle 419 cplex_rhoij=pawrhoij(iatom)%cplex 420 nspden_rhoij=pawrhoij(iatom)%nspden 421 do jlmn=1,pawrhoij(iatom)%lmn_size 422 j0lmn=jlmn*(jlmn-1)/2 423 cpj0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,jlmn) 424 dcpj0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,jlmn) 425 cpj0 (1:2,2) =cwaveprj(iatm,2)%cp (1:2 ,jlmn) 426 dcpj0(1:2,2,1)=cwaveprj(iatm,2)%dcp(1:2,1,jlmn) 427 do ilmn=1,jlmn 428 klmn=j0lmn+ilmn 429 klmn_re=cplex_rhoij*(klmn-1)+1 430 cpi0 (1:2,1) =cwaveprj(iatm,1)%cp (1:2 ,ilmn) 431 dcpi0(1:2,1,1)=cwaveprj(iatm,1)%dcp(1:2,1,ilmn) 432 cpi0 (1:2,2) =cwaveprj(iatm,2)%cp (1:2 ,ilmn) 433 dcpi0(1:2,2,1)=cwaveprj(iatm,2)%dcp(1:2,1,ilmn) 434 ro11_re=zero;ro22_re=zero 435 ro12_re=zero;ro21_re=zero 436 ro12_im=zero;ro21_im=zero 437 do iplex=1,cplex 438 ro11_re=dcpi0(iplex,1,1)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,1) 439 ro22_re=dcpi0(iplex,2,1)*cpj0(iplex,2)+cpi0(iplex,2)*dcpj0(iplex,2,1) 440 end do 441 if (substract_diagonal) then 442 do iplex=1,cplex 443 ro11_re=ro11_re-cpi0(iplex,1)*cpj0(iplex,1) 444 ro22_re=ro22_re-cpi0(iplex,2)*cpj0(iplex,2) 445 end do 446 end if 447 pawrhoij(iatom)%rhoij_(klmn_re,1)=pawrhoij(iatom)%rhoij_(klmn_re,1)+weight*(ro11_re+ro22_re) 448 if (nspden_rhoij>1) then 449 do iplex=1,cplex 450 ro12_re=dcpi0(iplex,2,1)*cpj0(iplex,1)+cpi0(iplex,2)*dcpj0(iplex,1,1) 451 ro21_re=dcpi0(iplex,1,1)*cpj0(iplex,2)+cpi0(iplex,1)*dcpj0(iplex,2,1) 452 end do 453 if (substract_diagonal) then 454 do iplex=1,cplex 455 ro12_re=ro12_re-cpi0(iplex,2)*cpj0(iplex,1) 456 ro21_re=ro21_re-cpi0(iplex,1)*cpj0(iplex,2) 457 end do 458 end if 459 pawrhoij(iatom)%rhoij_(klmn_re,4)=pawrhoij(iatom)%rhoij_(klmn_re,4)+weight*(ro11_re-ro22_re) 460 pawrhoij(iatom)%rhoij_(klmn_re,2)=pawrhoij(iatom)%rhoij_(klmn_re,2)+weight*(ro12_re+ro21_re) 461 if (cplex==2) then 462 ro12_im=dcpi0(1,2,1)*cpj0(2,1)-dcpi0(2,2,1)*cpj0(1,1)+cpi0(1,2)*dcpj0(2,1,1)-cpi0(2,2)*dcpj0(1,1,1) 463 ro21_im=dcpi0(1,1,1)*cpj0(2,2)-dcpi0(2,1,1)*cpj0(1,2)+cpi0(1,1)*dcpj0(2,2,1)-cpi0(2,1)*dcpj0(1,2,1) 464 if (substract_diagonal) then 465 ro12_im=ro12_im-cpi0(1,2)*cpj0(2,1)+cpi0(2,2)*cpj0(1,1) 466 ro21_im=ro21_im-cpi0(1,1)*cpj0(2,2)+cpi0(2,1)*cpj0(1,2) 467 end if 468 pawrhoij(iatom)%rhoij_(klmn_re,3)=pawrhoij(iatom)%rhoij_(klmn_re,3)+weight*(ro21_im-ro12_im) 469 end if 470 end if 471 if (compute_impart) then 472 klmn_im=klmn_re+1 473 if (nspden_rhoij>1) pawrhoij(iatom)%rhoij_(klmn_im,3)=pawrhoij(iatom)%rhoij_(klmn_im,3)+weight*(ro12_re-ro21_re) 474 if (cplex==2) then 475 ro11_im=dcpi0(1,1,1)*cpj0(2,1)-dcpi0(2,1,1)*cpj0(1,1)+cpi0(1,1)*dcpj0(2,1,1)-cpi0(2,1)*dcpj0(1,1,1) 476 ro22_im=dcpi0(1,2,1)*cpj0(2,2)-dcpi0(2,2,1)*cpj0(1,2)+cpi0(1,2)*dcpj0(2,2,1)-cpi0(2,2)*dcpj0(1,2,1) 477 if (substract_diagonal) then 478 ro11_im=ro11_im-cpi0(1,1)*cpj0(2,1)+cpi0(2,1)*cpj0(1,1) 479 ro22_im=ro22_im-cpi0(1,2)*cpj0(2,2)+cpi0(2,2)*cpj0(1,2) 480 end if 481 pawrhoij(iatom)%rhoij_(klmn_im,1)=pawrhoij(iatom)%rhoij_(klmn_im,1)+weight*(ro11_im+ro22_im) 482 if (nspden_rhoij>1) then 483 pawrhoij(iatom)%rhoij_(klmn_im,4)=pawrhoij(iatom)%rhoij_(klmn_im,4)+weight*(ro11_im-ro22_im) 484 pawrhoij(iatom)%rhoij_(klmn_im,2)=pawrhoij(iatom)%rhoij_(klmn_im,2)+weight*(ro12_im+ro21_im) 485 end if 486 end if 487 end if 488 end do 489 end do 490 end do 491 end if 492 end if 493 494 else if (option==3) then 495 496 ! ================================================================== 497 ! === OPTION 3: Accumulate (n,k) contribution to drhoij/dr ========= 498 ! ================================================================== 499 500 compute_impart=((.not.usetimerev).and.(pawrhoij(1)%cplex==2)) 501 compute_impart_cplex=((compute_impart).and.(cplex==2)) 502 if (nspinor==1) then 503 do iatom=1,my_natom 504 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 505 iatm=atindx(iatom1) 506 cplex_rhoij=pawrhoij(iatom)%cplex 507 do jlmn=1,pawrhoij(iatom)%lmn_size 508 j0lmn=jlmn*(jlmn-1)/2 509 cpj0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,jlmn) 510 dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn) 511 do ilmn=1,jlmn 512 klmn=j0lmn+ilmn 513 klmn_re=cplex_rhoij*(klmn-1)+1 514 cpi0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,ilmn) 515 dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn) 516 do mu=1,ncpgr 517 ro11_re=zero 518 do iplex=1,cplex 519 ro11_re=ro11_re+dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu) 520 end do 521 pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_re,isppol)+weight*ro11_re 522 end do 523 if (compute_impart_cplex) then 524 klmn_im=klmn_re+1 525 do mu=1,ncpgr 526 ro11_im=dcpi0(1,1,mu)*cpj0(2,1)+cpi0(1,1)*dcpj0(2,1,mu)-dcpi0(2,1,mu)*cpj0(1,1)-cpi0(2,1)*dcpj0(1,1,mu) 527 pawrhoij(iatom)%grhoij(mu,klmn_im,isppol)=pawrhoij(iatom)%grhoij(mu,klmn_im,isppol)+weight*ro11_im 528 end do 529 end if 530 end do 531 end do 532 end do 533 else ! nspinor=2 534 do iatom=1,my_natom 535 iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom) 536 iatm=atindx(iatom1) 537 cplex_rhoij=pawrhoij(iatom)%cplex 538 nspden_rhoij=pawrhoij(iatom)%nspden 539 do jlmn=1,pawrhoij(iatom)%lmn_size 540 j0lmn=jlmn*(jlmn-1)/2 541 cpj0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,jlmn) 542 cpj0(1:cplex,2) =cwaveprj(iatm,2)%cp (1:cplex,jlmn) 543 dcpj0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,jlmn) 544 dcpj0(1:cplex,2,1:ncpgr)=cwaveprj(iatm,2)%dcp(1:cplex,1:ncpgr,jlmn) 545 do ilmn=1,jlmn 546 klmn=j0lmn+ilmn 547 klmn_re=cplex_rhoij*(klmn-1)+1 548 klmn_im=klmn_re+1 549 cpi0(1:cplex,1) =cwaveprj(iatm,1)%cp (1:cplex,ilmn) 550 cpi0(1:cplex,2) =cwaveprj(iatm,2)%cp (1:cplex,ilmn) 551 dcpi0(1:cplex,1,1:ncpgr)=cwaveprj(iatm,1)%dcp(1:cplex,1:ncpgr,ilmn) 552 dcpi0(1:cplex,2,1:ncpgr)=cwaveprj(iatm,2)%dcp(1:cplex,1:ncpgr,ilmn) 553 do mu=1,ncpgr 554 ro11_re=zero;ro22_re=zero 555 ro12_re=zero;ro21_re=zero 556 ro12_im=zero;ro21_im=zero 557 do iplex=1,cplex 558 ro11_re=dcpi0(iplex,1,mu)*cpj0(iplex,1)+cpi0(iplex,1)*dcpj0(iplex,1,mu) 559 ro22_re=dcpi0(iplex,2,mu)*cpj0(iplex,2)+cpi0(iplex,2)*dcpj0(iplex,2,mu) 560 end do 561 pawrhoij(iatom)%grhoij(mu,klmn_re,1)=pawrhoij(iatom)%grhoij(mu,klmn_re,1)+weight*(ro11_re+ro22_re) 562 if (nspden_rhoij>1) then 563 do iplex=1,cplex 564 ro12_re=ro12_re+dcpi0(iplex,2,mu)*cpj0(iplex,1)+cpi0(iplex,2)*dcpj0(iplex,1,mu) 565 ro21_re=ro21_re+dcpi0(iplex,1,mu)*cpj0(iplex,2)+cpi0(iplex,1)*dcpj0(iplex,2,mu) 566 end do 567 pawrhoij(iatom)%grhoij(mu,klmn_re,4)=pawrhoij(iatom)%grhoij(mu,klmn_re,4)+weight*(ro11_re-ro22_re) 568 pawrhoij(iatom)%grhoij(mu,klmn_re,2)=pawrhoij(iatom)%grhoij(mu,klmn_re,2)+weight*(ro12_re+ro21_re) 569 if (cplex==2) then 570 ro12_im=dcpi0(1,2,mu)*cpj0(2,1)+cpi0(1,2)*dcpj0(2,1,mu)-dcpi0(2,2,mu)*cpj0(1,1)-cpi0(2,2)*dcpj0(1,1,mu) 571 ro21_im=dcpi0(1,1,mu)*cpj0(2,2)+cpi0(1,1)*dcpj0(2,2,mu)-dcpi0(2,1,mu)*cpj0(1,2)-cpi0(2,1)*dcpj0(1,2,mu) 572 pawrhoij(iatom)%grhoij(mu,klmn_re,3)=pawrhoij(iatom)%grhoij(mu,klmn_re,3)+weight*(ro21_im-ro12_im) 573 end if 574 end if 575 if (compute_impart) then 576 if (nspden_rhoij>1) then 577 pawrhoij(iatom)%grhoij(mu,klmn_im,3)=pawrhoij(iatom)%grhoij(mu,klmn_im,3)+weight*(ro12_re-ro21_re) 578 end if 579 if (cplex==2) then 580 ro11_im=dcpi0(1,1,mu)*cpj0(2,1)+cpi0(1,1)*dcpj0(2,1,mu)-dcpi0(2,1,mu)*cpj0(1,1)-cpi0(2,1)*dcpj0(1,1,mu) 581 ro22_im=dcpi0(1,2,mu)*cpj0(2,2)+cpi0(1,2)*dcpj0(2,2,mu)-dcpi0(2,2,mu)*cpj0(1,2)-cpi0(2,2)*dcpj0(1,2,mu) 582 pawrhoij(iatom)%grhoij(mu,klmn_im,1)=pawrhoij(iatom)%grhoij(mu,klmn_im,1)+weight*(ro11_im+ro22_im) 583 if (nspden_rhoij>1) then 584 pawrhoij(iatom)%grhoij(mu,klmn_im,4)=pawrhoij(iatom)%grhoij(mu,klmn_im,4)+weight*(ro11_im-ro22_im) 585 pawrhoij(iatom)%grhoij(mu,klmn_im,2)=pawrhoij(iatom)%grhoij(mu,klmn_im,2)+weight*(ro12_im+ro21_im) 586 end if 587 end if 588 end if 589 end do 590 end do 591 end do 592 end do 593 end if 594 595 ! End 596 end if ! option 597 598 !Destroy atom table used for parallelism 599 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 600 601 DBG_EXIT("COLL") 602 603 end subroutine pawaccrhoij