TABLE OF CONTENTS
ABINIT/pawdenpot [ Functions ]
NAME
pawdenpot
FUNCTION
Compute different (PAW) energies, densities and potentials (or potential-like quantities) inside PAW spheres Can also compute first-order densities potentials and second-order energies (RF calculations).
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (FJ, 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
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument) [hyb_mixing, hyb_mixing_sr]= -- optional-- mixing factors for the global (resp. screened) XC hybrid functional ipert=index of perturbation (used only for RF calculation ; set ipert<=0 for GS calculations. ixc= choice of exchange-correlation scheme (see above, and below) 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=total number of atoms in cell nspden=number of spin-density components ntypat=number of types of atoms in unit cell. nucdipmom(3,my_natom) nuclear dipole moments nzlmopt= if -1, compute all LM-moments of densities initialize "lmselect" (index of non-zero LM-moments of densities) if 0, compute all LM-moments of densities force "lmselect" to .true. (index of non-zero LM-moments of densities) if 1, compute only non-zero LM-moments of densities (stored before) option=0: compute both energies and potentials 1: compute only potentials 2: compute only energies paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh paw_an0(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh for Ground-State paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawprtvol=control print volume and debugging output for PAW pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawspnorb=flag: 1 if spin-orbit coupling is activated pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) ucvol=unit cell volume (bohr^3) xclevel= XC functional level xc_denpos= lowest allowed density (usually for the computation of the XC functionals) znucl(ntypat)=gives the nuclear charge for all types of atoms
OUTPUT
paw_ij(my_natom)%dijhartree(cplex*lmn2_size)=Hartree contribution to dij; Enters into calculation of hartree energy ==== if option=0 or 2 epaw=contribution to total energy from the PAW "on-site" part epawdc=contribution to total double-counting energy from the PAW "on-site" part ==== if option=0 or 2 and ipert<=0 compch_sph=compensation charge integral inside spheres computed over spherical meshes ==== if (option=0 or 1) and paw_an(:)%has_vxc=1 paw_an(my_natom)%vxc1[m](cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" density paw_an(my_natom)%vxct1[m](cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" pseudo density ==== if paw_an(iatom_tot)%has_vxcval==1 compute also XC potentials neglecting core charge paw_an(my_natom)%vxc1_val[m](cplex*mesh_size,:nspden)=XC potential calculated from spherical valence density paw_an(my_natom)%vxct1_val[m](cplex*mesh_size,:nspden)=XC potential calculated from spherical valence pseudo density ==== if nzlmopt==-1, paw_an(iatom_tot)%lnmselect(lm_size,nspden)=select the non-zero LM-moments of rho1 and trho1 ==== if paw_an(:)%has_vhartree=1 paw_an(my_natom)%vh1(cplex*mesh_size,1,1)=Hartree total potential calculated from "on-site" density ==== if pawspnorb>0 paw_ij(my_natom)%dijso(cplex_dij*lmn2_size,nspden)=spin-orbit contribution to dij
NOTES
Response function calculations:
In order to compute first- or second-order qunatities, paw_an (resp. paw_ij) datastructures
must contain first-order quantities, namely paw_an1 (resp. paw_ij1).
PARENTS
bethe_salpeter,dfpt_scfcv,odamix,paw_qpscgw,respfn,scfcv,screening sigma
CHILDREN
free_my_atmtab,get_my_atmtab,pawdensities,pawdijfock,pawdijhartree pawdijnd,pawdijso,pawrad_deducer0,pawuenergy,pawxc,pawxc_dfpt,pawxcm pawxcm_dfpt,pawxcmpositron,pawxcpositron,pawxenergy,pawxpot,poisson setnoccmmp,timab,wrtout,xmpi_sum
SOURCE
92 #if defined HAVE_CONFIG_H 93 #include "config.h" 94 #endif 95 96 #include "abi_common.h" 97 98 subroutine pawdenpot(compch_sph,epaw,epawdc,ipert,ixc,& 99 & my_natom,natom,nspden,ntypat,nucdipmom,nzlmopt,option,paw_an,paw_an0,& 100 & paw_ij,pawang,pawprtvol,pawrad,pawrhoij,pawspnorb,pawtab,pawxcdev,spnorbscl,xclevel,xc_denpos,ucvol,znucl,& 101 & electronpositron,mpi_atmtab,comm_atom,vpotzero,hyb_mixing,hyb_mixing_sr) ! optional arguments 102 103 use defs_basis 104 use m_profiling_abi 105 use m_errors 106 use m_xmpi 107 108 use m_pawang, only: pawang_type 109 use m_pawrad, only: pawrad_type, pawrad_deducer0, poisson, simp_gen 110 use m_pawtab, only: pawtab_type 111 use m_paw_an, only : paw_an_type 112 use m_paw_ij, only : paw_ij_type 113 use m_pawrhoij,only: pawrhoij_type 114 use m_pawdij, only: pawdijhartree, pawdijnd, pawdijso, pawxpot ,pawdijfock 115 use m_pawxc, only: pawxc, pawxc_dfpt, pawxcm, pawxcm_dfpt, pawxcpositron, pawxcmpositron 116 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 117 use m_electronpositron, only: electronpositron_type,electronpositron_calctype 118 119 !This section has been created automatically by the script Abilint (TD). 120 !Do not modify the following lines by hand. 121 #undef ABI_FUNC 122 #define ABI_FUNC 'pawdenpot' 123 use interfaces_14_hidewrite 124 use interfaces_18_timing 125 use interfaces_65_paw, except_this_one => pawdenpot 126 !End of the abilint section 127 128 implicit none 129 130 !Arguments --------------------------------------------- 131 !scalars 132 integer,intent(in) :: ipert,ixc,my_natom,natom,nspden,ntypat,nzlmopt,option,pawprtvol 133 integer,intent(in) :: pawspnorb,pawxcdev,xclevel 134 integer,optional,intent(in) :: comm_atom 135 real(dp), intent(in) :: spnorbscl,xc_denpos,ucvol 136 real(dp),intent(in),optional :: hyb_mixing,hyb_mixing_sr 137 real(dp),intent(out) :: compch_sph,epaw,epawdc 138 type(electronpositron_type),pointer,optional :: electronpositron 139 type(pawang_type),intent(in) :: pawang 140 !arrays 141 integer,optional,target,intent(in) :: mpi_atmtab(:) 142 real(dp),intent(in) :: nucdipmom(3,my_natom),znucl(ntypat) 143 real(dp),intent(out),optional :: vpotzero(2) 144 type(paw_an_type),intent(inout) :: paw_an(my_natom) 145 type(paw_an_type), intent(in) :: paw_an0(my_natom) 146 type(paw_ij_type),intent(inout) :: paw_ij(my_natom) 147 type(pawrad_type),intent(in) :: pawrad(ntypat) 148 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 149 type(pawtab_type),intent(in) :: pawtab(ntypat) 150 151 !Local variables --------------------------------------- 152 !scalars 153 integer :: cplex,cplex_dij,has_kxc,iatom,iatom_tot,idum,ierr,ipositron,irhoij,ispden,itypat,itypat0 154 integer :: jrhoij,kklmn,klmn,lm_size,lmn2_size,mesh_size,my_comm_atom,ndij,nkxc1,nspdiag,nsppol,opt_compch 155 integer :: usecore,usepawu,usetcore,usexcnhat,usenhat,usefock 156 logical :: keep_vhartree,my_atmtab_allocated,need_kxc,non_magnetic_xc,paral_atom,temp_vxc 157 real(dp) :: e1t10,e1xc,e1xcdc,efock,efockdc,eexc,eexcdc,eexdctemp 158 real(dp) :: eexc_val,eexcdc_val,eexex,eexexdc,eextemp,eh2 159 real(dp) :: eldaumdc,eldaumdcdc,enucdip,espnorb,etild1xc,etild1xcdc 160 real(dp) :: exccore,exchmix,hyb_mixing_,hyb_mixing_sr_,rdum 161 character(len=3) :: pertstrg 162 character(len=500) :: msg 163 !arrays 164 integer :: idum1(0),idum3(0,0,0) 165 integer,pointer :: my_atmtab(:) 166 logical,allocatable :: lmselect_cur(:),lmselect_cur_ep(:),lmselect_ep(:),lmselect_tmp(:) 167 real(dp) :: ro(2),mpiarr(7),tsec(2) 168 real(dp),allocatable :: dij_ep(:),dijfock_vv(:,:),dijfock_cv(:,:),one_over_rad2(:),kxc_tmp(:,:,:),nhat1(:,:,:),nhat1_ep(:,:,:) 169 real(dp) :: rdum2(0,0),rdum3(0,0,0),rdum3a(0,0,0),rdum4(0,0,0,0) 170 real(dp),allocatable :: rho(:),rho1(:,:,:),rho1_ep(:,:,:),rho1xx(:,:,:) 171 real(dp),allocatable :: trho1(:,:,:),trho1_ep(:,:,:),vxc_tmp(:,:,:) 172 173 ! ************************************************************************* 174 175 DBG_ENTER("COLL") 176 177 call timab(560,1,tsec) 178 179 if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then 180 msg='invalid value for variable "nzlmopt".' 181 MSG_BUG(msg) 182 end if 183 184 if (my_natom>0) then 185 if(paw_ij(1)%has_dijhartree==0.and. .not.(ipert==natom+1.or.ipert==natom+10)) then 186 msg='dijhartree must be allocated !' 187 MSG_BUG(msg) 188 end if 189 if (paw_ij(1)%cplex/=paw_an(1)%cplex) then 190 msg='paw_ij()%cplex and paw_an()%cplex must be equal !' 191 MSG_BUG(msg) 192 end if 193 if (pawrhoij(1)%cplex<paw_an(1)%cplex) then 194 msg='pawrhoij()%cplex must be >=paw_an()%cplex !' 195 MSG_BUG(msg) 196 end if 197 if (ipert<=0.and.paw_ij(1)%cplex/=1) then 198 msg='cplex must be 1 for GS calculations !' 199 MSG_BUG(msg) 200 end if 201 if (ipert>0.and.(ipert<=natom.or.ipert==natom+2).and.paw_an0(1)%has_kxc/=2) then 202 msg='XC kernels for ground state must be in memory !' 203 MSG_BUG(msg) 204 end if 205 if(paw_an(1)%has_vxc==0.and.(option==0.or.option==1).and. & 206 & .not.(ipert==natom+1.or.ipert==natom+10)) then 207 msg='vxc1 and vxct1 must be allocated !' 208 MSG_BUG(msg) 209 end if 210 if (ipert>0.and.paw_an(1)%has_vhartree==1) then 211 msg='computation of vhartree not compatible with RF (ipert>0) !' 212 MSG_BUG(msg) 213 end if 214 if (ipert>0.and.paw_an(1)%has_vxcval==1.and.(option==0.or.option==1)) then 215 msg='computation of vxc_val not compatible with RF (ipert>0) !' 216 MSG_BUG(msg) 217 end if 218 end if 219 220 ipositron=0 221 if (present(electronpositron)) then 222 ipositron=electronpositron_calctype(electronpositron) 223 if (ipositron==1.and.pawtab(1)%has_kij/=2) then 224 msg='kij must be in memory for electronpositron%calctype=1 !' 225 MSG_BUG(msg) 226 end if 227 if (ipert>0) then 228 msg='electron-positron calculation not available for ipert>0 !' 229 MSG_ERROR(msg) 230 end if 231 end if 232 233 !Set up parallelism over atoms 234 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 235 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 236 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 237 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 238 & my_natom_ref=my_natom) 239 240 !For some perturbations, nothing to do 241 if (ipert==natom+1.or.ipert==natom+10) then 242 if (option/=1) then 243 epaw=zero;epawdc=zero 244 end if 245 return 246 end if 247 248 !Various inits 249 hyb_mixing_ =zero ; if(present(hyb_mixing)) hyb_mixing_ =hyb_mixing 250 hyb_mixing_sr_=zero ; if(present(hyb_mixing_sr)) hyb_mixing_sr_=hyb_mixing_sr 251 usefock=0;if (abs(hyb_mixing_)>tol8.or.abs(hyb_mixing_sr_)>tol8) usefock=1 252 usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat) 253 usenhat = usexcnhat 254 keep_vhartree=(maxval(paw_an(:)%has_vhartree)>0) 255 if(keep_vhartree) usenhat = 1 256 usepawu=maxval(pawtab(1:ntypat)%usepawu) 257 non_magnetic_xc=(usepawu==4).or.(usepawu==14) 258 compch_sph=-1.d5 259 opt_compch=0;if (option/=1.and.ipert<=0) opt_compch=1 260 if (opt_compch==1) compch_sph=zero 261 nspdiag=1;if (nspden==2) nspdiag=2 262 nsppol=1;if (my_natom>0) nsppol=pawrhoij(1)%nsppol 263 pertstrg=" ";if (ipert>0) pertstrg="(1)" 264 ro(:)=zero 265 266 !Init energies 267 if (option/=1) then 268 e1xc=zero ; e1xcdc=zero 269 etild1xc=zero ; etild1xcdc=zero 270 exccore=zero ; eh2=zero ; e1t10=zero 271 eldaumdc=zero ; eldaumdcdc=zero 272 eexex=zero ; eexexdc=zero 273 eextemp=zero ; eexdctemp=zero 274 espnorb=zero ; enucdip=zero 275 efock=zero ; efockdc=zero 276 if (ipositron/=0) then 277 electronpositron%e_paw =zero 278 electronpositron%e_pawdc=zero 279 end if 280 end if 281 !vpotzero is needed for both the energy and the potential 282 if (present(vpotzero)) vpotzero(:)=zero 283 284 !if PAW+U, compute noccmmp^{\sigma}_{m,m'} occupation matrix 285 if (usepawu>0.and.ipert<=0.and.ipositron/=1) then 286 if (paral_atom) then 287 call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,nsppol,0,ntypat,& 288 & paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu,& 289 & comm_atom=my_comm_atom,mpi_atmtab=mpi_atmtab) 290 else 291 call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,nsppol,0,ntypat,& 292 & paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu) 293 end if 294 end if 295 296 !Print some titles 297 if (abs(pawprtvol)>=2) then 298 if (nzlmopt<1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,& 299 ' ====== Moments of (n1-tn1)',trim(pertstrg),' =========' 300 if (nzlmopt==1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,& 301 ' ==== Non-zero Moments of (n1-tn1)',trim(pertstrg),' ====' 302 call wrtout(std_out,msg,'COLL') 303 if (usexcnhat/=0) then 304 write(msg, '(6a)')' The moments of (n1-tn1-nhat1)',trim(pertstrg),' must be very small...' 305 call wrtout(std_out,msg,'COLL') 306 end if 307 end if 308 309 !================ Big loop on atoms ======================= 310 !========================================================== 311 312 do iatom=1,my_natom 313 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 314 itypat=pawrhoij(iatom)%itypat 315 exchmix=pawtab(itypat)%exchmix 316 lmn2_size=paw_ij(iatom)%lmn2_size 317 lm_size=paw_an(iatom)%lm_size 318 mesh_size=pawtab(itypat)%mesh_size 319 320 usecore=1;usetcore =pawtab(itypat)%usetcore 321 if (ipert/=0) usecore=0 ! This is true for phonons and Efield pert. 322 if (ipert/=0) usetcore=0 ! This is true for phonons and Efield pert. 323 has_kxc=paw_an(iatom)%has_kxc;need_kxc=(has_kxc==1) 324 cplex=1;if (ipert>0) cplex=pawrhoij(iatom)%cplex 325 cplex_dij=paw_ij(iatom)%cplex_dij 326 ndij=paw_ij(iatom)%ndij 327 328 ! Allocations of "on-site" densities 329 ABI_ALLOCATE(rho1 ,(cplex*mesh_size,lm_size,nspden)) 330 ABI_ALLOCATE(trho1,(cplex*mesh_size,lm_size,nspden)) 331 ABI_ALLOCATE(nhat1,(cplex*mesh_size,lm_size,nspden*usenhat)) 332 rho1(:,:,:)=zero;trho1(:,:,:)=zero;nhat1(:,:,:)=zero 333 if (ipositron/=0) then ! Additional allocation for the electron-positron case 334 ABI_ALLOCATE(rho1_ep ,(cplex*mesh_size,lm_size,nspden)) 335 ABI_ALLOCATE(trho1_ep,(cplex*mesh_size,lm_size,nspden)) 336 ABI_ALLOCATE(nhat1_ep,(cplex*mesh_size,lm_size,nspden*usenhat)) 337 end if 338 ABI_ALLOCATE(lmselect_cur,(lm_size)) 339 lmselect_cur(:)=.true. 340 if (nzlmopt==1) lmselect_cur(:)=paw_an(iatom)%lmselect(:) 341 342 ! Store some usefull quantities 343 itypat0=0;if (iatom>1) itypat0=pawrhoij(iatom-1)%itypat 344 if (itypat/=itypat0) then 345 ABI_ALLOCATE(one_over_rad2,(mesh_size)) 346 one_over_rad2(1)=zero 347 one_over_rad2(2:mesh_size)=one/pawrad(itypat)%rad(2:mesh_size)**2 348 end if 349 350 ! Need to allocate vxc1 in particular cases 351 if (pawspnorb>0.and.ipert==0.and.option==2.and.ipositron/=1.and. & 352 & pawrhoij(iatom)%cplex==2.and.paw_an(iatom)%has_vxc==0) then 353 ! these should already be allocated in paw_an_init! 354 if (pawxcdev==0)then 355 if (allocated(paw_an(iatom)%vxc1)) then 356 ABI_DEALLOCATE(paw_an(iatom)%vxc1) 357 end if 358 ABI_ALLOCATE(paw_an(iatom)%vxc1,(cplex*mesh_size,paw_an(iatom)%angl_size,nspden)) 359 end if 360 if (pawxcdev/=0 ) then 361 if (allocated(paw_an(iatom)%vxc1)) then 362 ABI_DEALLOCATE(paw_an(iatom)%vxc1) 363 end if 364 ABI_ALLOCATE(paw_an(iatom)%vxc1,(cplex*mesh_size,lm_size,nspden)) 365 end if 366 paw_an(iatom)%has_vxc=1 367 temp_vxc=.true. 368 else 369 temp_vxc=.false. 370 end if 371 372 ! ===== Compute "on-site" densities (n1, ntild1, nhat1) ===== 373 ! ========================================================== 374 375 call pawdensities(compch_sph,cplex,iatom_tot,lmselect_cur,paw_an(iatom)%lmselect,lm_size,& 376 & nhat1,nspden,nzlmopt,opt_compch,1-usenhat,-1,1,pawang,pawprtvol,pawrad(itypat),& 377 & pawrhoij(iatom),pawtab(itypat),rho1,trho1,one_over_rad2=one_over_rad2) 378 379 if (ipositron/=0) then 380 ! Electron-positron calculation: need additional on-site densities: 381 ! if ipositron==1, need electronic on-site densities 382 ! if ipositron==2, need positronic on-site densities 383 ABI_ALLOCATE(lmselect_ep,(lm_size)) 384 ABI_ALLOCATE(lmselect_cur_ep,(lm_size)) 385 lmselect_cur_ep(:)=.true. 386 if (nzlmopt==1) lmselect_cur_ep(:)=electronpositron%lmselect_ep(1:lm_size,iatom) 387 388 call pawdensities(rdum,cplex,iatom_tot,lmselect_cur_ep,lmselect_ep,& 389 & lm_size,nhat1_ep,nspden,nzlmopt,0,1-usenhat,-1,0,pawang,0,pawrad(itypat),& 390 & electronpositron%pawrhoij_ep(iatom),pawtab(itypat),& 391 & rho1_ep,trho1_ep,one_over_rad2=one_over_rad2) 392 393 if (nzlmopt<1) electronpositron%lmselect_ep(1:lm_size,iatom)=lmselect_ep(1:lm_size) 394 ABI_DEALLOCATE(lmselect_ep) 395 ABI_DEALLOCATE(lmselect_cur_ep) 396 end if 397 398 ! =========== Compute XC potentials and energies =========== 399 ! ========================================================== 400 401 ! Temporary storage 402 nkxc1=0;if (paw_an(iatom)%has_kxc/=0) nkxc1=paw_an(iatom)%nkxc1 403 if (pawxcdev/=0) then 404 ABI_ALLOCATE(vxc_tmp,(cplex*mesh_size,lm_size,nspden)) 405 if (need_kxc) then 406 ABI_ALLOCATE(kxc_tmp,(mesh_size,lm_size,nkxc1)) 407 end if 408 end if 409 if (pawxcdev==0) then 410 ABI_ALLOCATE(vxc_tmp,(cplex*mesh_size,pawang%angl_size,nspden)) 411 vxc_tmp(:,:,:)=zero 412 if (need_kxc) then 413 ABI_ALLOCATE(kxc_tmp,(mesh_size,pawang%angl_size,nkxc1)) 414 end if 415 end if 416 idum=0 417 if (.not.allocated(kxc_tmp)) then 418 ABI_ALLOCATE(kxc_tmp,(0,0,0)) 419 end if 420 421 ! ===== Vxc1 term ===== 422 if (ipositron/=1) then 423 if (pawxcdev/=0) then 424 if (ipert==0) then 425 call pawxcm(pawtab(itypat)%coredens,eexc,eexcdc,idum,ixc,kxc_tmp,lm_size,& 426 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 427 & pawang,pawrad(itypat),pawxcdev,rho1,usecore,0,vxc_tmp,xclevel,xc_denpos) 428 else 429 call pawxcm_dfpt(pawtab(itypat)%coredens,cplex,cplex,eexc,ixc,paw_an0(iatom)%kxc1,lm_size,& 430 & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,option,& 431 & pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel) 432 eexcdc=zero 433 end if 434 else 435 if (ipert==0) then 436 call pawxc(pawtab(itypat)%coredens,eexc,eexcdc,ixc,kxc_tmp,lm_size,& 437 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 438 & pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel,xc_denpos) 439 else 440 call pawxc_dfpt(pawtab(itypat)%coredens,cplex,cplex,eexc,ixc,paw_an0(iatom)%kxc1,lm_size,& 441 & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,option,& 442 & pawang,pawrad(itypat),rho1,usecore,0,paw_an0(iatom)%vxc1,vxc_tmp,xclevel) 443 eexcdc=zero 444 end if 445 end if 446 447 if (option/=1) then 448 e1xc=e1xc+eexc 449 e1xcdc=e1xcdc+eexcdc 450 end if 451 if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=vxc_tmp(:,:,:) 452 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=kxc_tmp(:,:,:) 453 else ! ipositron==1 454 if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=zero 455 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=zero 456 end if 457 458 459 ! Additional electron-positron XC term (if ipositron/=0) 460 if (ipositron/=0) then 461 if (pawxcdev/=0) then 462 call pawxcmpositron(ipositron,pawtab(itypat)%coredens,eexc,eexcdc,electronpositron%ixcpositron,& 463 & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),& 464 & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,& 465 & electronpositron%posdensity0_limit,rho1,rho1_ep,usecore,0,vxc_tmp,xc_denpos) 466 else 467 call pawxcpositron(ipositron,pawtab(itypat)%coredens,eexc,eexcdc,electronpositron%ixcpositron,& 468 & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),& 469 & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),& 470 & electronpositron%posdensity0_limit,rho1,rho1_ep,usecore,0,vxc_tmp,xc_denpos) 471 end if 472 if (option/=1) then 473 electronpositron%e_paw =electronpositron%e_paw +eexc 474 electronpositron%e_pawdc=electronpositron%e_pawdc+eexcdc 475 end if 476 if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=paw_an(iatom)%vxc1(:,:,:)+vxc_tmp(:,:,:) 477 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=paw_an(iatom)%kxc1(:,:,:)+kxc_tmp(:,:,:) 478 end if 479 480 ! ===== tVxc1 term ===== 481 if (ipositron/=1) then 482 if (pawxcdev/=0) then 483 if (ipert==0) then 484 call pawxcm(pawtab(itypat)%tcoredens(:,1),& 485 & eexc,eexcdc,idum,ixc,kxc_tmp,lm_size,& 486 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 487 & pawang,pawrad(itypat),pawxcdev,trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel,xc_denpos) 488 else 489 call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),& 490 & cplex,cplex,eexc,ixc,paw_an0(iatom)%kxct1,lm_size,& 491 & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,option,& 492 & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel) 493 eexcdc=zero 494 end if 495 else 496 if (ipert==0) then 497 call pawxc(pawtab(itypat)%tcoredens(:,1),& 498 & eexc,eexcdc,ixc,kxc_tmp,lm_size,& 499 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 500 & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel,xc_denpos) 501 else 502 call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),& 503 & cplex,cplex,eexc,ixc,paw_an0(iatom)%kxct1,lm_size,& 504 & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,option,& 505 & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,paw_an0(iatom)%vxct1,vxc_tmp,xclevel) 506 eexcdc=zero 507 end if 508 end if 509 if (option/=1) then 510 etild1xc=etild1xc+eexc 511 etild1xcdc=etild1xcdc+eexcdc 512 end if 513 if (option<2) paw_an(iatom)%vxct1(:,:,:)=vxc_tmp(:,:,:) 514 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=kxc_tmp(:,:,:) 515 else ! ipositron==1 516 if (option<2) paw_an(iatom)%vxct1(:,:,:)=zero 517 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=zero 518 end if 519 520 ! Additional electron-positron XC term (if ipositron/=0) 521 if (ipositron/=0) then 522 if (pawxcdev/=0) then 523 call pawxcmpositron(ipositron,pawtab(itypat)%tcoredens(:,1),& 524 & eexc,eexcdc,electronpositron%ixcpositron,& 525 & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),& 526 & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,& 527 & electronpositron%posdensity0_limit,trho1,trho1_ep,usetcore,2*usexcnhat,vxc_tmp,xc_denpos) 528 else 529 call pawxcpositron(ipositron,pawtab(itypat)%tcoredens(:,1),& 530 & eexc,eexcdc,electronpositron%ixcpositron,& 531 & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),& 532 & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),& 533 & electronpositron%posdensity0_limit,trho1,trho1_ep,usetcore,2*usexcnhat,vxc_tmp,xc_denpos) 534 end if 535 if (option/=1) then 536 electronpositron%e_paw =electronpositron%e_paw -eexc 537 electronpositron%e_pawdc=electronpositron%e_pawdc-eexcdc 538 end if 539 if (option<2) paw_an(iatom)%vxct1(:,:,:)=paw_an(iatom)%vxct1(:,:,:)+vxc_tmp(:,:,:) 540 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=paw_an(iatom)%kxct1(:,:,:)+kxc_tmp(:,:,:) 541 end if 542 543 ! Update flags defining the state of vxc and kxc 544 if (option<2) paw_an(iatom)%has_vxc=2 545 if (need_kxc.and.nkxc1>0) paw_an(iatom)%has_kxc=2 546 547 ! Update core XC conjtribution to energy 548 if (option/=1.and.ipositron/=1) exccore=exccore+pawtab(itypat)%exccore 549 if (ipositron/=0) then 550 ABI_DEALLOCATE(rho1_ep) 551 ABI_DEALLOCATE(trho1_ep) 552 ABI_DEALLOCATE(nhat1_ep) 553 end if 554 555 ! =========== Compute valence-only XC potentials =========== 556 ! ========================================================== 557 if (ipert==0.and.paw_an(iatom)%has_vxcval==1.and.(option==0.or.option==1)) then 558 if (.not.allocated(paw_an(iatom)%vxc1_val).or..not.allocated(paw_an(iatom)%vxct1_val)) then 559 msg=' vxc1_val and vxct1_val must be associated' 560 MSG_BUG(msg) 561 end if 562 ! ===== Vxc1_val term, vxc[n1] ===== 563 if (pawxcdev/=0) then 564 write(msg,'(4a,es16.6)')ch10,& 565 & ' pawdenpot : Computing valence-only v_xc[n1] using moments ',ch10,& 566 & ' Min density rho1 = ',MINVAL(rho1) 567 call wrtout(std_out,msg,'COLL') 568 call pawxcm(pawtab(itypat)%coredens,eexc_val,eexcdc_val,idum,ixc,kxc_tmp,lm_size,& 569 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 570 & pawang,pawrad(itypat),pawxcdev,rho1,0,0,vxc_tmp,xclevel,xc_denpos) 571 else 572 write(msg,'(2a)')ch10,' pawdenpot : Computing valence-only v_xc[n1] using angular mesh ' 573 call wrtout(std_out,msg,'COLL') 574 575 call pawxc(pawtab(itypat)%coredens,eexc_val,eexcdc_val,ixc,kxc_tmp,lm_size,& 576 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 577 & pawang,pawrad(itypat),rho1,0,0,vxc_tmp,xclevel,xc_denpos) 578 end if 579 if (option<2) paw_an(iatom)%vxc1_val(:,:,:)=vxc_tmp(:,:,:) 580 581 ! ===== tVxc1_val term ===== 582 if (pawxcdev/=0) then 583 if (usexcnhat/=0) then 584 write(msg,'(4a,e16.6,2a,es16.6)')ch10,& 585 & ' pawdenpot : Computing valence-only v_xc[tn1+nhat] using moments ',ch10,& 586 & ' Min density trho1 = ',MINVAL(trho1),ch10,& 587 & ' Min density trho1 + nhat = ',MINVAL(trho1+nhat1) 588 else 589 write(msg,'(4a,e16.6)')ch10,& 590 & ' pawdenpot : Computing valence-only v_xc[tn1] using moments ',ch10,& 591 & ' Min density trho1 = ',MINVAL(trho1) 592 end if 593 call wrtout(std_out,msg,'COLL') 594 call pawxcm(pawtab(itypat)%tcoredens(:,1),& 595 & eexc_val,eexcdc_val,idum,ixc,kxc_tmp,lm_size,& 596 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 597 & pawang,pawrad(itypat),pawxcdev,trho1,0,2*usexcnhat,vxc_tmp,xclevel,xc_denpos) 598 else 599 write(msg,'(2a)')ch10,' pawdenpot : Computing valence-only v_xc[tn1+nhat] using angular mesh' 600 call wrtout(std_out,msg,'COLL') 601 call pawxc(pawtab(itypat)%tcoredens(:,1),& 602 & eexc_val,eexcdc_val,ixc,kxc_tmp,lm_size,& 603 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 604 & pawang,pawrad(itypat),trho1,0,2*usexcnhat,vxc_tmp,xclevel,xc_denpos) 605 end if 606 if (option<2) then 607 paw_an(iatom)%vxct1_val(:,:,:)=vxc_tmp(:,:,:) 608 paw_an(iatom)%has_vxcval=2 609 end if 610 end if ! valence-only XC potentials 611 612 ABI_DEALLOCATE(vxc_tmp) 613 ABI_DEALLOCATE(kxc_tmp) 614 615 ! ===== Compute first part of local exact-exchange energy term ===== 616 ! ===== Also compute corresponding potential ===== 617 ! ================================================================== 618 619 if (pawtab(itypat)%useexexch>0.and.ipert==0.and.ipositron/=1) then 620 621 ! ===== Re-compute a partial "on-site" density n1 (only l=lexexch contrib.) 622 ABI_ALLOCATE(rho1xx,(mesh_size,lm_size,nspden)) 623 ABI_ALLOCATE(lmselect_tmp,(lm_size)) 624 lmselect_tmp(:)=lmselect_cur(:) 625 call pawdensities(rdum,cplex,iatom_tot,lmselect_cur,lmselect_tmp,lm_size,rdum3,nspden,& 626 & 1,0,2,pawtab(itypat)%lexexch,0,pawang,pawprtvol,pawrad(itypat),& 627 & pawrhoij(iatom),pawtab(itypat),rho1xx,rdum3a,one_over_rad2=one_over_rad2) 628 ABI_DEALLOCATE(lmselect_tmp) 629 ! ===== Re-compute Exc1 and Vxc1; for local exact-exchange, this is done in GGA only 630 ABI_ALLOCATE(vxc_tmp,(mesh_size,lm_size,nspden)) 631 ABI_ALLOCATE(kxc_tmp,(mesh_size,lm_size,nkxc1)) 632 call pawxcm(pawtab(itypat)%coredens,eextemp,eexdctemp,pawtab(itypat)%useexexch,ixc,kxc_tmp,lm_size,& 633 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,& 634 & rho1xx,0,0,vxc_tmp,xclevel,xc_denpos) 635 if (option/=1) then 636 e1xc=e1xc-eextemp*exchmix 637 e1xcdc=e1xcdc-eexdctemp*exchmix 638 end if 639 if (option<2) then 640 paw_an(iatom)%vxc_ex(:,:,:)=vxc_tmp(:,:,:) 641 paw_an(iatom)%has_vxc_ex=2 642 end if 643 ABI_DEALLOCATE(rho1xx) 644 ABI_DEALLOCATE(vxc_tmp) 645 ABI_DEALLOCATE(kxc_tmp) 646 647 end if ! useexexch 648 649 itypat0=0;if (iatom<my_natom) itypat0=pawrhoij(iatom+1)%itypat 650 if (itypat/=itypat0) then 651 ABI_DEALLOCATE(one_over_rad2) 652 end if 653 654 ABI_DEALLOCATE(lmselect_cur) 655 656 ! ==== Compute Hartree potential terms and some energy terms ==== 657 ! =============================================================== 658 659 ! Electron-positron calculation: compute Dij due to fixed particles (elec. or pos. depending on calctype) 660 if (ipositron/=0) then 661 ABI_ALLOCATE(dij_ep,(cplex*lmn2_size)) 662 call pawdijhartree(paw_ij(iatom)%cplex,dij_ep,nspden,electronpositron%pawrhoij_ep(iatom),pawtab(itypat)) 663 if (option/=1) then 664 do ispden=1,nspdiag 665 jrhoij=1 666 do irhoij=1,pawrhoij(iatom)%nrhoijsel 667 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 668 ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn) 669 electronpositron%e_paw =electronpositron%e_paw -ro(1)*dij_ep(klmn) 670 electronpositron%e_pawdc=electronpositron%e_pawdc-ro(1)*dij_ep(klmn) 671 if (ipositron==1) e1t10=e1t10+ro(1)*two*(pawtab(itypat)%kij(klmn)-pawtab(itypat)%dij0(klmn)) 672 jrhoij=jrhoij+pawrhoij(iatom)%cplex 673 end do 674 end do 675 end if 676 end if 677 678 ! Hartree Dij computation 679 if (ipositron/=1) then 680 call pawdijhartree(paw_ij(iatom)%cplex,paw_ij(iatom)%dijhartree,nspden,pawrhoij(iatom),pawtab(itypat)) 681 else 682 paw_ij(iatom)%dijhartree(:)=zero 683 end if 684 paw_ij(iatom)%has_dijhartree=2 685 686 ! Hartree energy computation 687 if (option/=1) then 688 if (cplex==1.or.ipert==0) then 689 do ispden=1,nspdiag 690 jrhoij=1 691 do irhoij=1,pawrhoij(iatom)%nrhoijsel 692 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 693 ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn) 694 eh2=eh2 +ro(1)*paw_ij(iatom)%dijhartree(klmn) 695 e1t10=e1t10+ro(1)*pawtab(itypat)%dij0(klmn) 696 jrhoij=jrhoij+pawrhoij(iatom)%cplex 697 end do 698 end do 699 else 700 do ispden=1,nspdiag 701 jrhoij=1 702 do irhoij=1,pawrhoij(iatom)%nrhoijsel 703 klmn=pawrhoij(iatom)%rhoijselect(irhoij);kklmn=2*klmn-1 704 ro(1:2)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+1,ispden)*pawtab(itypat)%dltij(klmn) 705 eh2=eh2+ro(1)*paw_ij(iatom)%dijhartree(kklmn)+ro(2)*paw_ij(iatom)%dijhartree(kklmn+1) 706 ! Imaginary part (not used) 707 ! eh2=eh2+ro(2)*paw_ij(iatom)%dijhartree(kklmn)-ro(1)*paw_ij(iatom)%dijhartree(kklmn+1) 708 jrhoij=jrhoij+pawrhoij(iatom)%cplex 709 end do 710 end do 711 end if 712 end if 713 714 ! Electron-positron calculation: add electron and positron 715 if (ipositron/=0) then 716 paw_ij(iatom)%dijhartree(:)=paw_ij(iatom)%dijhartree(:)-dij_ep(:) 717 ABI_DEALLOCATE(dij_ep) 718 end if 719 720 ! Compute 1st moment of total Hartree potential VH(n_Z+n_core+n1) 721 ! equation 10 (density) and up to 43 (Hartree potential of density) 722 ! of Kresse and Joubert PRB 59 1758 (1999) 723 keep_vhartree=(paw_an(iatom)%has_vhartree>0) 724 if ((pawspnorb>0.and.ipert==0.and.ipositron/=1).or.keep_vhartree) then 725 ! in the first clause case, would it not be simpler just to turn on has_vhartree? 726 if (.not. allocated(paw_an(iatom)%vh1)) then 727 ABI_ALLOCATE(paw_an(iatom)%vh1,(mesh_size,1,1)) 728 end if 729 if (.not. allocated(paw_an(iatom)%vht1)) then 730 ABI_ALLOCATE(paw_an(iatom)%vht1,(mesh_size,1,1)) 731 end if 732 733 ! construct vh1 734 ! the sqrt(4pi) factor comes from the fact we are calculating the spherical moments, 735 ! and for the 00 channel the prefactor of Y_00 = 2 sqrt(pi) 736 ABI_ALLOCATE(rho,(mesh_size)) 737 rho(1:mesh_size)=(rho1(1:mesh_size,1,1) + sqrt(four_pi)*pawtab(itypat)%coredens(1:mesh_size)) & 738 & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2 739 ! rho(1:mesh_size)=rho1(1:mesh_size,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2 740 741 call poisson(rho,0,pawrad(itypat),paw_an(iatom)%vh1(:,1,1)) 742 743 paw_an(iatom)%vh1(2:mesh_size,1,1)=(paw_an(iatom)%vh1(2:mesh_size,1,1) & 744 & -sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size) 745 ! TODO: check this is equivalent to the previous version (commented) which explicitly recalculated VH(coredens) 746 ! DONE: numerically there are residual differences on abiref (7th digit). 747 ! paw_an(iatom)%vh1(2:mesh_size,1,1)=paw_an(iatom)%vh1(2:mesh_size,1,1)/pawrad(itypat)%rad(2:mesh_size) & 748 !& + sqrt(four_pi) * pawtab(itypat)%VHnZC(2:mesh_size) 749 750 call pawrad_deducer0(paw_an(iatom)%vh1(:,1,1),mesh_size,pawrad(itypat)) 751 752 ! same for vht1 753 rho = zero 754 if (usenhat /= 0) then 755 rho(1:mesh_size)=nhat1(1:mesh_size,1,1) 756 end if 757 rho(1:mesh_size)=(rho(1:mesh_size) + trho1(1:mesh_size,1,1) + sqrt(four_pi)*pawtab(itypat)%tcoredens(1:mesh_size,1)) & 758 & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2 759 ! rho(1:mesh_size)=(rho(1:mesh_size) + trho1(1:mesh_size,1,1))*four_pi*pawrad(itypat)%rad(1:mesh_size)**2 760 761 call poisson(rho,0,pawrad(itypat),paw_an(iatom)%vht1(:,1,1)) 762 763 paw_an(iatom)%vht1(2:mesh_size,1,1)=(paw_an(iatom)%vht1(2:mesh_size,1,1) & 764 & -sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size) 765 ! paw_an(iatom)%vht1(2:mesh_size,1,1)=paw_an(iatom)%vht1(2:mesh_size,1,1)/pawrad(itypat)%rad(2:mesh_size) & 766 !& + sqrt(four_pi)*pawtab(itypat)%vhtnzc(2:mesh_size) 767 call pawrad_deducer0(paw_an(iatom)%vht1(:,1,1),mesh_size,pawrad(itypat)) 768 769 paw_an(iatom)%has_vhartree=2 770 ABI_DEALLOCATE(rho) 771 end if 772 773 ! ========= Compute PAW+U and energy contribution ========= 774 ! ========================================================== 775 776 if (pawtab(itypat)%usepawu>0.and.ipert==0.and.ipositron/=1.and.option/=1.and.pawtab(itypat)%usepawu<10) then 777 call pawuenergy(iatom_tot,eldaumdc,eldaumdcdc,pawprtvol,pawtab(itypat),paw_ij(iatom)) 778 end if 779 780 ! ========= Compute nuclear dipole moment energy contribution ======== 781 ! ========================================================== 782 783 ! Compute nuclear dipole contribution to Dij if necessary 784 if (any(abs(nucdipmom(:,iatom))>tol8).and.ipert==0.and.ipositron/=1.and.option/=1) then 785 if (paw_ij(iatom)%has_dijnd/=2) then 786 call pawdijnd(cplex_dij,paw_ij(iatom)%dijnd,ndij,nucdipmom(:,iatom),pawrad(itypat),pawtab(itypat)) 787 paw_ij(iatom)%has_dijnd=2 788 end if 789 end if ! end compute dijnd if necessary 790 791 ! Compute contribution to on-site energy 792 if (any(abs(nucdipmom(:,iatom))>tol8).and.ipert==0.and.ipositron/=1.and.option/=1) then 793 if(pawrhoij(iatom)%cplex/=2) then 794 msg=' pawrhoij must be complex for nuclear dipole moments !' 795 MSG_BUG(msg) 796 end if 797 jrhoij=2 !Select imaginary part of rhoij 798 do irhoij=1,pawrhoij(iatom)%nrhoijsel 799 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 800 ! select imaginary part of dijnd (only nonzero part) 801 kklmn=cplex_dij*(klmn-1)+2 802 ! because dijnd involves no spin flips, following formula is correct for all values of nspden 803 enucdip=enucdip-pawrhoij(iatom)%rhoijp(jrhoij,1)*paw_ij(iatom)%dijnd(kklmn,1) & 804 & *pawtab(itypat)%dltij(klmn) 805 jrhoij=jrhoij+pawrhoij(iatom)%cplex 806 end do 807 end if 808 809 ! ========= Compute spin-orbit energy contribution ======== 810 ! ========================================================== 811 812 ! Compute spin-orbit contribution to Dij 813 if (pawspnorb>0.and.ipert==0.and.ipositron/=1.and.& 814 & (option/=2.or.pawrhoij(iatom)%cplex==2)) then 815 call pawdijso(cplex_dij,paw_ij(iatom)%dijso,ndij,nspden,pawang,pawrad(itypat),pawtab(itypat), & 816 & pawxcdev,spnorbscl,paw_an(iatom)%vh1,paw_an(iatom)%vxc1) 817 paw_ij(iatom)%has_dijso=2 818 if (.not.keep_vhartree) then 819 paw_an(iatom)%has_vhartree=0 820 ABI_DEALLOCATE(paw_an(iatom)%vh1) 821 end if 822 if (temp_vxc) then 823 paw_an(iatom)%has_vxc=0 824 ABI_DEALLOCATE(paw_an(iatom)%vxc1) 825 end if 826 end if 827 828 ! Compute contribution to on-site energy 829 if (option/=1.and.pawspnorb>0.and.ipert==0.and.ipositron/=1.and.pawrhoij(iatom)%cplex==2) then 830 if(pawrhoij(iatom)%nspden/=4) then 831 msg=' pawrhoij must have 4 components !' 832 MSG_BUG(msg) 833 end if 834 jrhoij=2 !Select imaginary part of rhoij 835 do irhoij=1,pawrhoij(iatom)%nrhoijsel 836 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 837 kklmn=cplex_dij*(klmn-1)+1 838 espnorb=espnorb-pawrhoij(iatom)%rhoijp(jrhoij,3)*paw_ij(iatom)%dijso(kklmn,3) & 839 & *pawtab(itypat)%dltij(klmn) 840 if (cplex_dij==2) then 841 kklmn=kklmn+1 842 espnorb=espnorb-(pawrhoij(iatom)%rhoijp(jrhoij,2)*paw_ij(iatom)%dijso(kklmn,3) & 843 & +half*pawrhoij(iatom)%rhoijp(jrhoij,4)*(paw_ij(iatom)%dijso(kklmn,1) & 844 & -paw_ij(iatom)%dijso(kklmn,2))) & 845 & *pawtab(itypat)%dltij(klmn) 846 end if 847 jrhoij=jrhoij+pawrhoij(iatom)%cplex 848 end do 849 end if 850 851 ! === Compute 2nd part of local exact-exchange energy and potential === 852 ! ====================================================================== 853 854 if (pawtab(itypat)%useexexch>0.and.ipert==0.and.ipositron/=1) then 855 856 if(paw_ij(iatom)%nspden==4) then 857 msg=' Local exact-exch. not implemented for nspden=4 !' 858 MSG_ERROR(msg) 859 end if 860 861 if (option<2) then 862 call pawxpot(ndij,pawprtvol,pawrhoij(iatom),pawtab(itypat),paw_ij(iatom)%vpawx) 863 paw_ij(iatom)%has_exexch_pot=2 864 end if 865 if (option/=1) then 866 if (abs(pawprtvol)>=2) then 867 write(msg, '(2a)' )ch10,'======= PAW local exact exchange terms (in Hartree) ====' 868 call wrtout(std_out, msg,'COLL') 869 write(msg, '(2a,i4)' )ch10,' For Atom',iatom_tot 870 call wrtout(std_out, msg,'COLL') 871 end if 872 call pawxenergy(eexex,pawprtvol,pawrhoij(iatom),pawtab(itypat)) 873 end if 874 875 end if ! useexexch 876 877 ! ==== Compute Fock Dij term and Fock energy terms ==== 878 ! ===================================================== 879 880 ! Computation of Fock contribution to Dij 881 if (usefock==1) then 882 883 if (ipositron/=1) then 884 885 ! Fock contribution to Dij 886 ABI_ALLOCATE(dijfock_vv,(cplex_dij*lmn2_size,ndij)) 887 ABI_ALLOCATE(dijfock_cv,(cplex_dij*lmn2_size,ndij)) 888 call pawdijfock(cplex,cplex_dij,dijfock_vv,dijfock_cv,hyb_mixing_,hyb_mixing_sr_, & 889 & ndij,nspden,nsppol,pawrhoij(iatom),pawtab(itypat)) 890 paw_ij(iatom)%dijfock(:,:)=dijfock_vv(:,:)+dijfock_cv(:,:) 891 892 ! Fock contribution to energy 893 894 if (option/=1) then 895 if ((cplex==1).or.(ipert==0)) then 896 do ispden=1,pawrhoij(iatom)%nspden 897 jrhoij=1 898 do irhoij=1,pawrhoij(iatom)%nrhoijsel 899 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 900 ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn) 901 efockdc=efockdc+ro(1)*half*dijfock_vv(klmn,ispden) 902 efock=efock+ro(1)*(half*dijfock_vv(klmn,ispden)+dijfock_cv(klmn,ispden)) 903 jrhoij=jrhoij+pawrhoij(iatom)%cplex 904 end do 905 end do 906 else 907 do ispden=1,pawrhoij(iatom)%nspden 908 jrhoij=1 909 do irhoij=1,pawrhoij(iatom)%nrhoijsel 910 klmn=pawrhoij(iatom)%rhoijselect(irhoij);kklmn=2*klmn-1 911 ro(1:2)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+1,ispden)*pawtab(itypat)%dltij(klmn) 912 efockdc=efockdc+half*(ro(1)*dijfock_vv(kklmn,ispden)+ro(2)*dijfock_vv(kklmn+1,ispden)) 913 efock=efock+ro(1)*(half*dijfock_vv(kklmn, ispden)+dijfock_cv(kklmn, ispden)) & 914 & +ro(2)*(half*dijfock_vv(kklmn+1,ispden)+dijfock_cv(kklmn+1,ispden)) 915 jrhoij=jrhoij+pawrhoij(iatom)%cplex 916 end do 917 end do 918 end if 919 end if 920 921 ABI_DEALLOCATE(dijfock_vv) 922 ABI_DEALLOCATE(dijfock_cv) 923 924 else 925 ! Not Fock contribution for positron 926 paw_ij(iatom)%dijfock(:,:)=zero 927 end if 928 paw_ij(iatom)%has_dijfock=2 929 end if 930 931 ! ========================================================== 932 ! No more need of densities 933 ABI_DEALLOCATE(rho1) 934 ABI_DEALLOCATE(trho1) 935 ABI_DEALLOCATE(nhat1) 936 937 ! === Compute the zero of the potentials if requested ================== 938 ! ====================================================================== 939 940 if (pawtab(itypat)%usepotzero==1 .AND. present(vpotzero)) then 941 942 ! term 1 : beta 943 vpotzero(1)=vpotzero(1)-pawtab(itypat)%beta/ucvol 944 945 ! term 2 : \sum_ij rho_ij gamma_ij 946 do ispden=1,nspdiag 947 jrhoij=1 948 do irhoij=1,pawrhoij(iatom)%nrhoijsel 949 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 950 vpotzero(2) = vpotzero(2) - & 951 & pawrhoij(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn)*pawtab(itypat)%gammaij(klmn)/ucvol 952 jrhoij=jrhoij+pawrhoij(iatom)%cplex 953 end do 954 end do 955 956 end if 957 958 ! =========== End loop on atoms ============================ 959 ! ========================================================== 960 961 end do 962 963 !========== Assemble "on-site" energy terms =============== 964 !========================================================== 965 966 if (option/=1) then 967 if (ipert==0) then 968 epaw=e1xc+half*eh2+e1t10-exccore-etild1xc+eldaumdc+eexex+espnorb+efock+enucdip 969 epawdc=e1xc-e1xcdc-half*eh2-exccore-etild1xc+etild1xcdc+eldaumdcdc-eexex-efockdc 970 else 971 epaw=e1xc-etild1xc+eh2 972 epawdc=zero 973 end if 974 end if 975 976 !========== Reduction in case of parallelism ============== 977 !========================================================== 978 979 if (paral_atom) then 980 if (option/=1) then 981 call timab(48,1,tsec) 982 mpiarr=zero 983 mpiarr(1)=compch_sph;mpiarr(2)=epaw;mpiarr(3)=epawdc 984 if (ipositron/=0) then 985 mpiarr(4)=electronpositron%e_paw 986 mpiarr(5)=electronpositron%e_pawdc 987 end if 988 if (present(vpotzero)) then 989 mpiarr(6)=vpotzero(1) 990 mpiarr(7)=vpotzero(2) 991 end if 992 call xmpi_sum(mpiarr,my_comm_atom,ierr) 993 compch_sph=mpiarr(1);epaw=mpiarr(2);epawdc=mpiarr(3) 994 if (ipositron/=0) then 995 electronpositron%e_paw=mpiarr(4) 996 electronpositron%e_pawdc=mpiarr(5) 997 end if 998 if (present(vpotzero)) then 999 vpotzero(1)=mpiarr(6) 1000 vpotzero(2)=mpiarr(7) 1001 end if 1002 call timab(48,2,tsec) 1003 end if 1004 end if 1005 1006 !Destroy atom table used for parallelism 1007 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1008 1009 call timab(560,2,tsec) 1010 1011 DBG_EXIT("COLL") 1012 1013 end subroutine pawdenpot