TABLE OF CONTENTS
ABINIT/pawdfptenergy [ Functions ]
NAME
pawdfptenergy
FUNCTION
This routine compute the Hartree+XC PAW on-site contributions to a 1st-order or 2nd-order energy. These contributions are equal to: E_onsite= Int{ VHxc[n1_a^(1);nc^(1)].n1_b } -Int{ VHxc[tild_n1_a^(j1)+hat_n1_a^(j1);tild_n_c^(j1)].(tild_n1_b+n1_b)^(j2) } Some typical uses: A-Contribution to non-stationary expression of the 2nd-order total energy: In that case, n1_a^(1)[r]=n1^(j1)[r] and n1_b[r]=delta_n1^(j2)[r] where j1 and j2 are two given perturbations, and delta_n1^(j)[r] is the 1s-order density only due to change of WF overlap. See PRB 78, 035105 (2008), Eq.(80) E_onsite= Int{ VHxc[n1^(j1);nc^(j1)].delta_n1^(j2) } -Int{ VHxc[tild_n1^(j1)+hat_n1^(j1);tild_n_c^(j1)].delta_(tild_n1+hat_n1)^(j2) } B-Contribution to first-order Fermi energy: In that case, n1_a^(1)[r]=n1^(j1)[r] and n1_b[r]=n1[r,EFermi] where j1 is the current perturbation, and n1[r,EFermi] is the density at Fermi level. E_onsite= Int{ VHxc[n1^(j1);nc^(j1)].n1[r,EFermi] } -Int{ VHxc[tild_n1^(j1)+hat_n1^(j1);tild_n_c^(j1)].(tild_n1+hat_n1)[r,EFermi] }
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
ipert1,ipert2=indexes of perturbations (j1) and (j2)
if ipert2<=0, we compute a first-order energy
if ipert2> 0, we compute a second-order energy
ixc= choice of exchange-correlation scheme
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
ntypat=number of types of atoms in unit cell.
nzlmopt_a= For the n1_a density:
if -1, compute all LM-moments of the density and use non-zero LM-moments
if 0, compute all LM-moments of the density and use all LM-moments
if +1, compute only non-zero LM-moments of the density (stored before)
nzlmopt_b= For the n1_b density:
if -1, compute all LM-moments of the density and use non-zero LM-moments
if 0, compute all LM-moments of the density and use all LM-moments
if +1, compute only non-zero LM-moments of the density (stored before)
paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
paw_an1(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh
This corresponds to (j1) perturbation
paw_ij1(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
This corresponds to (j1) perturbation
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_a(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies for the (j1) perturbation
pawrhoij_b(natom) <type(pawrhoij_type)>=
if ipert2> 0: paw rhoij 1st-order occupancies for the (j2) perturbation corrsponding to n1_b^(j2)[r]
if ipert2<=0: paw rhoij occupancies corresponding to n1_b[r]
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)
xclevel= XC functional level
OUTPUT
delta_energy(2)= real and imaginary parts of contributions to non-stationary expression for the
second derivative of the total energy
SIDE EFFECTS
==== if paw_an1(:)%has_vxc<2, compute 1st-order XC potentials
paw_an1(natom)%vxc1(cplex_a*mesh_size,:,nspden) =AE 1st-order XC potential Vxc^(j1)
paw_an1(natom)%vxct1(cplex_a*mesh_size,:,nspden)=PS 1st-order XC potential tVxc^(j1)
==== if paw_ij1(:)%has_dijhartree<2, compute 1st-order Dij_hartree
paw_ij1(natom)%dijhartree(cplex_a*lmn2_size)=Hartree contribution to Dij^(j1)
PARENTS
dfpt_nstpaw,newfermie1
CHILDREN
free_my_atmtab,get_my_atmtab,pawdensities,pawdijhartree,pawxc_dfpt pawxcm_dfpt,timab,xmpi_sum
SOURCE
89 #if defined HAVE_CONFIG_H 90 #include "config.h" 91 #endif 92 93 #include "abi_common.h" 94 95 subroutine pawdfptenergy(delta_energy,ipert1,ipert2,ixc,my_natom,natom,ntypat,nzlmopt_a,nzlmopt_b,& 96 & paw_an0,paw_an1,paw_ij1,pawang,pawprtvol,pawrad,pawrhoij_a,pawrhoij_b,& 97 & pawtab,pawxcdev,xclevel, & 98 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 99 100 101 use defs_basis 102 use m_profiling_abi 103 use m_errors 104 use m_xmpi, only : xmpi_comm_self,xmpi_sum 105 106 use m_pawang, only : pawang_type 107 use m_pawrad, only : pawrad_type 108 use m_pawtab, only : pawtab_type 109 use m_paw_an, only : paw_an_type 110 use m_paw_ij, only : paw_ij_type 111 use m_pawrhoij, only : pawrhoij_type 112 use m_pawdij, only : pawdijhartree 113 use m_pawxc, only : pawxc_dfpt, pawxcm_dfpt 114 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 115 116 !This section has been created automatically by the script Abilint (TD). 117 !Do not modify the following lines by hand. 118 #undef ABI_FUNC 119 #define ABI_FUNC 'pawdfptenergy' 120 use interfaces_18_timing 121 use interfaces_65_paw, except_this_one => pawdfptenergy 122 !End of the abilint section 123 124 implicit none 125 126 !Arguments --------------------------------------------- 127 !scalars 128 integer,intent(in) :: ipert1,ipert2,ixc,my_natom,natom,ntypat,nzlmopt_a,nzlmopt_b 129 integer,intent(in) :: pawprtvol,pawxcdev,xclevel 130 integer,optional,intent(in) :: comm_atom 131 type(pawang_type),intent(in) :: pawang 132 !arrays 133 integer,optional,target,intent(in) :: mpi_atmtab(:) 134 real(dp),intent(out) :: delta_energy(2) 135 type(paw_an_type),intent(in) :: paw_an0(my_natom) 136 type(paw_an_type),intent(inout) :: paw_an1(my_natom) 137 type(paw_ij_type),intent(inout) :: paw_ij1(my_natom) 138 type(pawrad_type),intent(in) :: pawrad(ntypat) 139 type(pawrhoij_type),intent(in) :: pawrhoij_a(my_natom),pawrhoij_b(my_natom) 140 type(pawtab_type),intent(in) :: pawtab(ntypat) 141 142 !Local variables --------------------------------------- 143 !scalars 144 integer :: cplex_a,cplex_b,cplex_dijh1,iatom,iatom_tot,ierr,irhoij,ispden,itypat,jrhoij 145 integer :: klmn,lm_size_a,lm_size_b,mesh_size,my_comm_atom,nspden,nspdiag,opt_compch,optexc,optvxc 146 integer :: usecore,usetcore,usexcnhat 147 logical :: my_atmtab_allocated,paral_atom 148 real(dp) :: compch,eexc,eexc_im 149 character(len=500) :: msg 150 !arrays 151 integer,pointer :: my_atmtab(:) 152 logical,allocatable :: lmselect_a(:),lmselect_b(:),lmselect_tmp(:) 153 real(dp) :: dij(2),delta_energy_h(2),delta_energy_xc(2),ro(2),tsec(2) 154 real(dp),allocatable :: kxc_dum(:,:,:),nhat1(:,:,:),rho1(:,:,:),trho1(:,:,:) 155 156 ! ************************************************************************* 157 158 DBG_ENTER("COLL") 159 160 call timab(567,1,tsec) 161 162 if (.not.(ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11 & 163 & .or.ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11)) then 164 if((abs(nzlmopt_a)/=1.and.nzlmopt_a/=0).or.(abs(nzlmopt_b)/=1.and.nzlmopt_b/=0)) then 165 msg='invalid value for nzlmopt !' 166 MSG_BUG(msg) 167 end if 168 if (my_natom>0) then 169 if(paw_ij1(1)%has_dijhartree==0) then 170 msg='dijhartree must be allocated !' 171 MSG_BUG(msg) 172 end if 173 if(paw_an1(1)%has_vxc==0) then 174 msg='vxc1 and vxct1 must be allocated !' 175 MSG_BUG(msg) 176 end if 177 if(paw_an0(1)%has_kxc==0) then 178 msg='kxc1 must be allocated !' 179 MSG_BUG(msg) 180 end if 181 if ((ipert1<=natom.or.ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11).and.paw_an0(1)%has_kxc/=2) then 182 msg='XC kernels for ground state must be in memory !' 183 MSG_BUG(msg) 184 end if 185 if (paw_ij1(1)%cplex/=paw_an1(1)%cplex) then 186 msg='paw_ij1()%cplex and paw_an1()%cplex must be equal !' 187 MSG_BUG(msg) 188 end if 189 if (pawrhoij_a(1)%cplex<paw_an1(1)%cplex.or.pawrhoij_b(1)%cplex<paw_an1(1)%cplex) then 190 msg='pawrhoij()%cplex must be >=paw_an1()%cplex !' 191 MSG_BUG(msg) 192 end if 193 if (pawrhoij_a(1)%nspden/=pawrhoij_b(1)%nspden) then 194 msg='pawrhoij_a()%nspden must =pawrhoij_b()%nspden !' 195 MSG_BUG(msg) 196 end if 197 end if 198 end if 199 200 !Set up parallelism over atoms 201 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 202 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 203 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 204 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 205 206 !Init contribution to 1st-order (or 2nd-order) energy 207 delta_energy(1:2)=zero 208 209 !For some perturbations, nothing else to do 210 if (ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11 .or. & 211 & ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11) return 212 213 !Various inits 214 opt_compch=0;optvxc=1;optexc=3 215 usecore=0;usetcore=0 ! This is true for phonons and Efield pert. 216 usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat) 217 delta_energy_xc(1:2)=zero;delta_energy_h(1:2)=zero 218 dij(1:2)=zero;ro(1:2)=zero 219 220 221 !================ Loop on atomic sites ======================= 222 do iatom=1,my_natom 223 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 224 225 itypat=pawrhoij_a(iatom)%itypat 226 mesh_size=pawtab(itypat)%mesh_size 227 nspden=pawrhoij_a(iatom)%nspden 228 cplex_a=pawrhoij_a(iatom)%cplex 229 cplex_b=pawrhoij_b(iatom)%cplex 230 cplex_dijh1=paw_ij1(iatom)%cplex 231 lm_size_a=paw_an1(iatom)%lm_size 232 if (ipert2<=0) lm_size_b=paw_an0(iatom)%lm_size 233 if (ipert2> 0) lm_size_b=paw_an1(iatom)%lm_size 234 235 ! If Vxc potentials are not in memory, compute them 236 if (paw_an1(iatom)%has_vxc/=2) then 237 ABI_ALLOCATE(rho1 ,(cplex_a*mesh_size,lm_size_a,nspden)) 238 ABI_ALLOCATE(trho1,(cplex_a*mesh_size,lm_size_a,nspden)) 239 ABI_ALLOCATE(nhat1,(cplex_a*mesh_size,lm_size_a,nspden*usexcnhat)) 240 ABI_ALLOCATE(lmselect_a,(lm_size_a)) 241 lmselect_a(:)=paw_an1(iatom)%lmselect(:) 242 ABI_ALLOCATE(lmselect_tmp,(lm_size_a)) 243 lmselect_tmp(:)=.true. 244 if (nzlmopt_a==1) lmselect_tmp(:)=lmselect_a(:) 245 ! Compute on-site 1st-order densities 246 call pawdensities(compch,cplex_a,iatom_tot,lmselect_tmp,lmselect_a,& 247 & lm_size_a,nhat1,nspden,nzlmopt_a,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,& 248 & pawrad(itypat),pawrhoij_a(iatom),pawtab(itypat),rho1,trho1) 249 ABI_DEALLOCATE(lmselect_tmp) 250 ! Compute on-site 1st-order xc potentials 251 if (pawxcdev/=0) then 252 call pawxcm_dfpt(pawtab(itypat)%coredens,cplex_a,cplex_a,eexc,ixc,paw_an0(iatom)%kxc1,& 253 & lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,optvxc,& 254 & pawang,pawrad(itypat),rho1,usecore,0,& 255 & paw_an1(iatom)%vxc1,xclevel) 256 call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),& 257 & cplex_a,cplex_a,eexc,ixc,paw_an0(iatom)%kxct1,& 258 & lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,optvxc,& 259 & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,& 260 & paw_an1(iatom)%vxct1,xclevel) 261 else 262 call pawxc_dfpt(pawtab(itypat)%coredens,cplex_a,cplex_a,eexc,ixc,paw_an0(iatom)%kxc1,& 263 & lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,optvxc,& 264 & pawang,pawrad(itypat),rho1,usecore,0,& 265 & paw_an0(iatom)%vxc1,paw_an1(iatom)%vxc1,xclevel) 266 call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),& 267 & cplex_a,cplex_a,eexc,ixc,paw_an0(iatom)%kxct1,& 268 & lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,optvxc,& 269 & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,& 270 & paw_an0(iatom)%vxct1,paw_an1(iatom)%vxct1,xclevel) 271 end if 272 273 paw_an1(iatom)%has_vxc=2 274 ABI_DEALLOCATE(lmselect_a) 275 ABI_DEALLOCATE(rho1) 276 ABI_DEALLOCATE(trho1) 277 ABI_DEALLOCATE(nhat1) 278 end if ! has_vxc 279 280 ! If Dij_hartree are not in memory, compute them 281 if (paw_ij1(iatom)%has_dijhartree/=2) then 282 call pawdijhartree(cplex_dijh1,paw_ij1(iatom)%dijhartree,paw_ij1(iatom)%nspden,& 283 & pawrhoij_a(iatom),pawtab(itypat)) 284 paw_ij1(iatom)%has_dijhartree=2 285 end if 286 287 ! Compute contribution to 1st-order (or 2nd-order) energy from 1st-order XC potential 288 ABI_ALLOCATE(rho1 ,(cplex_b*mesh_size,lm_size_b,nspden)) 289 ABI_ALLOCATE(trho1,(cplex_b*mesh_size,lm_size_b,nspden)) 290 ABI_ALLOCATE(nhat1,(cplex_b*mesh_size,lm_size_b,nspden*usexcnhat)) 291 ABI_ALLOCATE(lmselect_b,(lm_size_b)) 292 if (ipert2<=0) lmselect_b(:)=paw_an0(iatom)%lmselect(:) 293 if (ipert2> 0) lmselect_b(:)=paw_an1(iatom)%lmselect(:) 294 ABI_ALLOCATE(lmselect_tmp,(lm_size_b)) 295 lmselect_tmp(:)=.true. 296 if (nzlmopt_b==1) lmselect_tmp(:)=lmselect_b(:) 297 ! Compute on-site 1st-order densities 298 call pawdensities(compch,cplex_b,iatom_tot,lmselect_tmp,lmselect_b,& 299 & lm_size_b,nhat1,nspden,nzlmopt_b,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,& 300 & pawrad(itypat),pawrhoij_b(iatom),pawtab(itypat),rho1,trho1) 301 ABI_DEALLOCATE(lmselect_tmp) 302 ! Compute contributions to 1st-order (or 2nd-order) energy 303 if (pawxcdev/=0) then 304 ABI_ALLOCATE(kxc_dum,(mesh_size,pawang%angl_size,0)) 305 call pawxcm_dfpt(pawtab(itypat)%coredens,cplex_b,cplex_a,eexc,ixc,kxc_dum,& 306 & lm_size_b,lmselect_b,nhat1,0,mesh_size,nspden,optexc,pawang,pawrad(itypat),& 307 & rho1,usecore,0,paw_an1(iatom)%vxc1,xclevel,d2enxc_im=eexc_im) 308 309 delta_energy_xc(1)=delta_energy_xc(1)+eexc 310 delta_energy_xc(2)=delta_energy_xc(2)+eexc_im 311 call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),& 312 & cplex_b,cplex_a,eexc,ixc,kxc_dum,& 313 & lm_size_b,lmselect_b,nhat1,0,mesh_size,nspden,optexc,pawang,pawrad(itypat),& 314 & trho1,usetcore,2*usexcnhat,paw_an1(iatom)%vxct1,xclevel,& 315 & d2enxc_im=eexc_im) 316 ABI_DEALLOCATE(kxc_dum) 317 delta_energy_xc(1)=delta_energy_xc(1)-eexc 318 delta_energy_xc(2)=delta_energy_xc(2)-eexc_im 319 else 320 ABI_ALLOCATE(kxc_dum,(mesh_size,lm_size_b,0)) 321 call pawxc_dfpt(pawtab(itypat)%coredens,cplex_b,cplex_a,eexc,ixc,kxc_dum,& 322 & lm_size_b,lmselect_b,nhat1,0,mesh_size,nspden,optexc,pawang,pawrad(itypat),& 323 & rho1,usecore,0,paw_an0(iatom)%vxc1,paw_an1(iatom)%vxc1,xclevel,d2enxc_im=eexc_im) 324 delta_energy_xc(1)=delta_energy_xc(1)+eexc 325 delta_energy_xc(2)=delta_energy_xc(2)+eexc_im 326 call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),& 327 & cplex_b,cplex_a,eexc,ixc,kxc_dum,& 328 & lm_size_b,lmselect_b,nhat1,0,mesh_size,nspden,optexc,pawang,pawrad(itypat),& 329 & trho1,usetcore,2*usexcnhat,paw_an0(iatom)%vxct1,paw_an1(iatom)%vxct1,xclevel,& 330 & d2enxc_im=eexc_im) 331 ABI_DEALLOCATE(kxc_dum) 332 delta_energy_xc(1)=delta_energy_xc(1)-eexc 333 delta_energy_xc(2)=delta_energy_xc(2)-eexc_im 334 end if 335 ABI_DEALLOCATE(lmselect_b) 336 ABI_DEALLOCATE(rho1) 337 ABI_DEALLOCATE(trho1) 338 ABI_DEALLOCATE(nhat1) 339 340 ! Compute contribution to 1st-order(or 2nd-order) energy from 1st-order Hartree potential 341 nspdiag=1;if (nspden==2) nspdiag=2 342 do ispden=1,nspdiag 343 if (cplex_dijh1==1) then 344 jrhoij=1 345 do irhoij=1,pawrhoij_b(iatom)%nrhoijsel 346 klmn=pawrhoij_b(iatom)%rhoijselect(irhoij) 347 dij(1)=paw_ij1(iatom)%dijhartree(klmn) 348 ro(1)=pawrhoij_b(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn) 349 delta_energy_h(1)=delta_energy_h(1)+ro(1)*dij(1) 350 if (cplex_b==2) then 351 ro(2)=pawrhoij_b(iatom)%rhoijp(jrhoij+1,ispden)*pawtab(itypat)%dltij(klmn) 352 delta_energy_h(2)=delta_energy_h(2)+ro(2)*dij(1) 353 end if 354 jrhoij=jrhoij+cplex_b 355 end do 356 else ! cplex_dijh1==2 357 jrhoij=1 358 do irhoij=1,pawrhoij_b(iatom)%nrhoijsel 359 klmn=pawrhoij_b(iatom)%rhoijselect(irhoij) 360 dij(1:2)=paw_ij1(iatom)%dijhartree(2*klmn-1:2*klmn) 361 ro(1)=pawrhoij_b(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn) 362 delta_energy_h(1)=delta_energy_h(1)+ro(1)*dij(1) 363 delta_energy_h(2)=delta_energy_h(2)-ro(1)*dij(2) 364 if (cplex_b==2) then 365 ro(2)=pawrhoij_b(iatom)%rhoijp(jrhoij+1,ispden)*pawtab(itypat)%dltij(klmn) 366 delta_energy_h(1)=delta_energy_h(1)+ro(2)*dij(2) 367 delta_energy_h(2)=delta_energy_h(2)+ro(2)*dij(1) 368 end if 369 jrhoij=jrhoij+cplex_b 370 end do 371 end if 372 end do 373 374 ! ================ End loop oon atomic sites ======================= 375 end do 376 377 !Final building of 1st-order (or 2nd-order) energy 378 delta_energy(1:2)=delta_energy_xc(1:2)+delta_energy_h(1:2) 379 380 !Reduction in case of parallelism 381 if (paral_atom) then 382 call xmpi_sum(delta_energy,my_comm_atom,ierr) 383 end if 384 385 !Destroy atom table used for parallelism 386 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 387 388 call timab(567,2,tsec) 389 390 DBG_EXIT("COLL") 391 392 end subroutine pawdfptenergy