TABLE OF CONTENTS
ABINIT/eltxccore [ Functions ]
NAME
eltxccore
FUNCTION
Compute the core charge contributions to the 2nd derivatives of the exchange-correlation energy with respect to all pairs of strain or strain and atomic displacement for the frozen wavefunction contribution to the elastic tensor. 1st-order potentials representing the perturbation by one strain are supplied, and the routine loops over the second strain and over all atomic displacements.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DRH, DCA, XG, GMR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
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. nfft=number of fft grid points ntypat=number of types of atoms in cell. n1,n2,n3=fft grid dimensions. n1xccc=dimension of xccc1d ; 0 if no XC core correction is used rprimd(3,3)=dimensional primitive translation vectors (bohr) typat(natom)=integer type for each atom in cell ucvol=unit cell volume (bohr**3). vxc_core(nfft)=spin-averaged xc potential vxc10_core(nfft)=spin-averaged 1st-order xc potential for elastic tensor vxc1is_core(nfft)=spin-averaged 1st-order xc potential for internal strain xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives, for each type of atom, from psp xred(3,natom)=reduced coordinates for atoms in unit cell
OUTPUT
eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the elastic tensor
SIDE EFFECTS
eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the elastic and internal-strain tensor. One column is incremented by the core contribution.
NOTES
Note that this routine is related to the mkcore.f routine
PARENTS
dfpt_eltfrxc
CHILDREN
free_my_atmtab,get_my_atmtab,timab,xmpi_sum
SOURCE
60 #if defined HAVE_CONFIG_H 61 #include "config.h" 62 #endif 63 64 #include "abi_common.h" 65 66 subroutine eltxccore(eltfrxc,is2_in,my_natom,natom,nfft,ntypat,& 67 & n1,n1xccc,n2,n3,rprimd,typat,ucvol,vxc_core,vxc10_core,vxc1is_core,& 68 & xcccrc,xccc1d,xred, & 69 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 70 71 use defs_basis 72 use m_errors 73 use m_profiling_abi 74 75 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 76 use m_xmpi, only : xmpi_comm_self,xmpi_sum 77 78 !This section has been created automatically by the script Abilint (TD). 79 !Do not modify the following lines by hand. 80 #undef ABI_FUNC 81 #define ABI_FUNC 'eltxccore' 82 use interfaces_18_timing 83 !End of the abilint section 84 85 implicit none 86 87 !Arguments ------------------------------------ 88 !scalars 89 integer,intent(in) :: is2_in,n1,n1xccc,n2,n3,my_natom,natom,nfft,ntypat 90 integer,optional,intent(in) :: comm_atom 91 real(dp),intent(in) :: ucvol 92 !arrays 93 integer,intent(in) :: typat(natom) 94 integer,optional,target,intent(in) :: mpi_atmtab(:) 95 real(dp),intent(in) :: rprimd(3,3),vxc10_core(nfft),vxc1is_core(nfft) 96 real(dp),intent(in) :: vxc_core(nfft),xccc1d(n1xccc,6,ntypat) 97 real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom) 98 real(dp),intent(inout) :: eltfrxc(6+3*natom,6) 99 100 !Local variables------------------------------- 101 !scalars 102 integer,parameter :: mshift=401 103 integer :: i1,i2,i3,iat,iatom,ierr,ifft,is1,is2,ishift,ishift1,ishift2 104 integer :: ishift3,ixp,jj,js,ka,kb,kd,kg,mu,my_comm_atom,nu 105 logical :: my_atmtab_allocated,paral_atom 106 real(dp) :: aa,bb,cc,d2rss,dd,delta,delta2div6,deltam1,diff 107 real(dp) :: difmag,difmag2,difmag2_fact,difmag2_part,drss1,drss2,func1 108 real(dp) :: func2,range,range2,rangem1,rdiff1,rdiff2,rdiff3 109 real(dp) :: term1,term2,yy 110 character(len=500) :: message 111 !arrays 112 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 113 integer :: igrid(3),ii(mshift,3),irange(3),ngfft(3) 114 integer,pointer :: my_atmtab(:) 115 real(dp) :: drm(3,3,6),eltfrxc_core(6+3*natom,6),lencp(3),rmet(3,3),rrdiff(mshift,3) 116 real(dp) :: scale(3),tau(3),ts2(3),tsec(2),tt(3) 117 real(dp),allocatable :: d2rm(:,:,:,:) 118 119 ! ************************************************************************* 120 121 !Compute lengths of cross products for pairs of primitive 122 !translation vectors (used in setting index search range below) 123 lencp(1)=cross_elt(rprimd(1,2),rprimd(2,2),rprimd(3,2),& 124 & rprimd(1,3),rprimd(2,3),rprimd(3,3)) 125 lencp(2)=cross_elt(rprimd(1,3),rprimd(2,3),rprimd(3,3),& 126 & rprimd(1,1),rprimd(2,1),rprimd(3,1)) 127 lencp(3)=cross_elt(rprimd(1,1),rprimd(2,1),rprimd(3,1),& 128 & rprimd(1,2),rprimd(2,2),rprimd(3,2)) 129 130 !Set up parallelism over atoms 131 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 132 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 133 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 134 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 135 136 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3 137 !(recall ucvol=R1.(R2xR3)) 138 scale(:)=ucvol/lencp(:) 139 140 !Compute metric tensor in real space rmet 141 do nu=1,3 142 rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+& 143 & rprimd(3,:)*rprimd(3,nu) 144 end do 145 146 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components 147 !and store for use in inner loop below. 148 149 ABI_ALLOCATE(d2rm,(3,3,6,6)) 150 151 !Loop over 2nd strain index 152 do is2=1,6 153 kg=idx(2*is2-1);kd=idx(2*is2) 154 do jj = 1,3 155 drm(:,jj,is2)=rprimd(kg,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kg,jj) 156 end do 157 158 ! Loop over 1st strain index 159 do is1=1,6 160 161 ka=idx(2*is1-1);kb=idx(2*is1) 162 d2rm(:,:,is1,is2)=0._dp 163 do jj = 1,3 164 if(ka==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)& 165 & +rprimd(kb,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kb,jj) 166 if(ka==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)& 167 & +rprimd(kb,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(kb,jj) 168 if(kb==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)& 169 & +rprimd(ka,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(ka,jj) 170 if(kb==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)& 171 & +rprimd(ka,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(ka,jj) 172 end do 173 end do !is1 174 end do !is2 175 176 ngfft(1)=n1 177 ngfft(2)=n2 178 ngfft(3)=n3 179 delta=1.0_dp/(n1xccc-1) 180 deltam1=n1xccc-1 181 delta2div6=delta**2/6.0_dp 182 183 !Loop over atoms in unit cell 184 eltfrxc_core(:,:)=zero 185 186 do iat=1,my_natom 187 iatom=iat;if (paral_atom) iatom=my_atmtab(iat) 188 js=7+3*(iatom-1) 189 ! Set search range (density cuts off perfectly beyond range) 190 ! Cycle if no range. 191 range=0.0_dp 192 range=xcccrc(typat(iatom)) 193 if(range<1.d-16) cycle 194 195 range2=range**2 196 rangem1=1.0_dp/range 197 198 ! Consider each component in turn 199 do mu=1,3 200 tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp) 201 202 ! Use tau to find nearest grid point along R(mu) 203 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 204 igrid(mu)=nint(tau(mu)*dble(ngfft(mu))) 205 206 ! Use range to compute an index range along R(mu) 207 ! (add 1 to make sure it covers full range) 208 irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu))) 209 210 ! Check that the largest range is smallest than the maximum 211 ! allowed one 212 if(2*irange(mu)+1 > mshift)then 213 write(message, '(a,i0,a)' )' The range around atom',iatom,' is too large.' 214 MSG_BUG(message) 215 end if 216 217 ! Set up a counter that explore the relevant range 218 ! of points around the atom 219 ishift=0 220 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 221 ishift=ishift+1 222 ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu)) 223 rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu) 224 end do 225 226 ! End loop on mu 227 end do 228 229 ! Conduct triple loop over restricted range of grid points for iatom 230 231 do ishift3=1,1+2*irange(3) 232 ! map back to [1,ngfft(3)] for usual fortran index in unit cell 233 i3=ii(ishift3,3) 234 ! find vector from atom location to grid point (reduced) 235 rdiff3=rrdiff(ishift3,3) 236 237 do ishift2=1,1+2*irange(2) 238 i2=ii(ishift2,2) 239 rdiff2=rrdiff(ishift2,2) 240 ! Prepare the computation of difmag2 241 difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2& 242 & +2.0_dp*rmet(3,2)*rdiff3*rdiff2 243 difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2) 244 245 do ishift1=1,1+2*irange(1) 246 rdiff1=rrdiff(ishift1,1) 247 248 ! Compute (rgrid-tau-Rprim)**2 249 difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1) 250 251 ! Only accept contribution inside defined range 252 if (difmag2<range2) then 253 254 ! Prepare computation of core charge function and derivatives, 255 ! using splines 256 difmag=sqrt(difmag2) 257 if (difmag>=1.0d-10) then 258 i1=ii(ishift1,1) 259 yy=difmag*rangem1 260 261 ! Compute index of yy over 1 to n1xccc scale 262 jj=1+int(yy*(n1xccc-1)) 263 diff=yy-(jj-1)*delta 264 265 ! Will evaluate spline fit (p. 86 Numerical Recipes, Press et al; 266 ! NOTE error in book for sign of "aa" term in derivative; 267 ! also see splfit routine). 268 bb = diff*deltam1 269 aa = 1.0_dp-bb 270 cc = aa*(aa**2-1.0_dp)*delta2div6 271 dd = bb*(bb**2-1.0_dp)*delta2div6 272 273 ! Evaluate spline fit of 1st der of core charge density 274 ! from xccc1d(:,2,:) and (:,4,:) 275 func1=aa*xccc1d(jj,2,typat(iatom))+bb*xccc1d(jj+1,2,typat(iatom)) +& 276 & cc*xccc1d(jj,4,typat(iatom))+dd*xccc1d(jj+1,4,typat(iatom)) 277 term1=func1*rangem1 278 ! Evaluate spline fit of 2nd der of core charge density 279 ! from xccc1d(:,3,:) and (:,5,:) 280 func2=aa*xccc1d(jj,3,typat(iatom))+bb*xccc1d(jj+1,3,typat(iatom)) +& 281 & cc*xccc1d(jj,5,typat(iatom))+dd*xccc1d(jj+1,5,typat(iatom)) 282 term2=func2*rangem1**2 283 284 ifft=i1+n1*(i2-1+n2*(i3-1)) 285 tt(:)=rmet(:,1)*rdiff1+rmet(:,2)*rdiff2+rmet(:,3)*rdiff3 286 287 ! Add contributions to 2nd derivative tensor 288 drss2=& 289 & (rdiff1*(drm(1,1,is2_in)*rdiff1+drm(1,2,is2_in)*rdiff2& 290 & +drm(1,3,is2_in)*rdiff3)& 291 & +rdiff2*(drm(2,1,is2_in)*rdiff1+drm(2,2,is2_in)*rdiff2& 292 & +drm(2,3,is2_in)*rdiff3)& 293 & +rdiff3*(drm(3,1,is2_in)*rdiff1+drm(3,2,is2_in)*rdiff2& 294 & +drm(3,3,is2_in)*rdiff3)) 295 296 ! Loop over 1st strain index 297 do is1=1,6 298 299 drss1=& 300 & (rdiff1*(drm(1,1,is1)*rdiff1+drm(1,2,is1)*rdiff2& 301 & +drm(1,3,is1)*rdiff3)& 302 & +rdiff2*(drm(2,1,is1)*rdiff1+drm(2,2,is1)*rdiff2& 303 & +drm(2,3,is1)*rdiff3)& 304 & +rdiff3*(drm(3,1,is1)*rdiff1+drm(3,2,is1)*rdiff2& 305 & +drm(3,3,is1)*rdiff3)) 306 307 d2rss=& 308 & (rdiff1*(d2rm(1,1,is1,is2_in)*rdiff1+d2rm(1,2,is1,is2_in)*rdiff2& 309 & +d2rm(1,3,is1,is2_in)*rdiff3)& 310 & +rdiff2*(d2rm(2,1,is1,is2_in)*rdiff1+d2rm(2,2,is1,is2_in)*rdiff2& 311 & +d2rm(2,3,is1,is2_in)*rdiff3)& 312 & +rdiff3*(d2rm(3,1,is1,is2_in)*rdiff1+d2rm(3,2,is1,is2_in)*rdiff2& 313 & +d2rm(3,3,is1,is2_in)*rdiff3)) 314 315 ! Vall(0) X Rhocore(2) term 316 eltfrxc_core(is1,is2_in)=eltfrxc_core(is1,is2_in)+0.25_dp*& 317 & (vxc_core(ifft)*(term1*(d2rss/difmag& 318 & -drss1*drss2/difmag**3)& 319 & +term2*drss1*drss2/difmag**2)) 320 321 ! Vall(1) X Rhocore(1) term 322 eltfrxc_core(is1,is2_in)=eltfrxc_core(is1,is2_in)+0.25_dp*& 323 & vxc10_core(ifft)*drss1*term1/difmag 324 eltfrxc_core(is2_in,is1)=eltfrxc_core(is2_in,is1)+0.25_dp*& 325 & vxc10_core(ifft)*drss1*term1/difmag 326 327 ! End loop in is1 328 end do 329 ! Internal strain contributions 330 ts2(:)=drm(:,1,is2_in)*rdiff1+drm(:,2,is2_in)*rdiff2& 331 & +drm(:,3,is2_in)*rdiff3 332 333 eltfrxc_core(js:js+2,is2_in)=eltfrxc_core(js:js+2,is2_in)& 334 & -(vxc1is_core(ifft)*term1/difmag& 335 & +0.5_dp*vxc_core(ifft)*(term2-term1/difmag)*drss2/difmag**2)*tt(:)& 336 & -(vxc_core(ifft)*term1/difmag)*ts2(:) 337 338 ! End of the condition for the distance not to vanish 339 end if 340 341 ! End of condition to be inside the range 342 end if 343 344 ! End loop on ishift1 345 end do 346 347 ! End loop on ishift2 348 end do 349 350 ! End loop on ishift3 351 end do 352 353 ! End loop on atoms 354 end do 355 356 !In case of parallelism over atoms: communicate 357 if (paral_atom) then 358 call timab(48,1,tsec) 359 call xmpi_sum(eltfrxc_core,my_comm_atom,ierr) 360 call timab(48,2,tsec) 361 end if 362 363 !Add core contribution to XC elastic tensor 364 eltfrxc(:,:)=eltfrxc(:,:)+eltfrxc_core(:,:) 365 366 !Destroy atom table used for parallelism 367 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 368 369 ABI_DEALLOCATE(d2rm) 370 371 contains 372 373 function cross_elt(xx,yy,zz,aa,bb,cc) 374 !Define magnitude of cross product of two vectors 375 376 !This section has been created automatically by the script Abilint (TD). 377 !Do not modify the following lines by hand. 378 #undef ABI_FUNC 379 #define ABI_FUNC 'cross_elt' 380 !End of the abilint section 381 382 real(dp) :: cross_elt 383 real(dp),intent(in) :: xx,yy,zz,aa,bb,cc 384 cross_elt=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2) 385 end function cross_elt 386 387 end subroutine eltxccore