TABLE OF CONTENTS
ABINIT/dfpt_mkcore [ Functions ]
NAME
dfpt_mkcore
FUNCTION
Compute the derivative of the core electron density with respect to one specific atom displacement In case of derivative with respect to k or electric (magnetic) field perturbation, the 1st-order core electron density vanishes.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, DRH) 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
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX idir=direction of atomic displacement (=1,2 or 3 : displacement of atom ipert along the 1st, 2nd or 3rd axis) or cartesian coordinate pair for strain perturbation ipert=number of the atom being displaced or natom+3,4 for strain perturbation natom=number of atoms in cell. 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 qphon(3)=wavevector of the phonon rprimd(3,3)=dimensional primitive translation vectors (bohr) typat(natom)=integer type for each atom in cell ucvol=unit cell volume (bohr**3). 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
xccc3d1(cplex*n1*n2*n3)=3D core electron density for XC core correction,
bohr^-3
NOTES
Note that this routine is tightly connected to the mkcore.f routine
PARENTS
dfpt_dyxc1,dfpt_eltfrxc,dfpt_looppert,dfpt_nselt,dfpt_nstdy,dfpt_nstpaw dfptnl_loop
CHILDREN
SOURCE
55 #if defined HAVE_CONFIG_H 56 #include "config.h" 57 #endif 58 59 #include "abi_common.h" 60 61 62 subroutine dfpt_mkcore(cplex,idir,ipert,natom,ntypat,n1,n1xccc,& 63 & n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xccc3d1,xred) 64 65 use defs_basis 66 use m_profiling_abi 67 use m_errors 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'dfpt_mkcore' 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: cplex,idir,ipert,n1,n1xccc,n2,n3,natom,ntypat 80 real(dp),intent(in) :: ucvol 81 !arrays 82 integer,intent(in) :: typat(natom) 83 real(dp),intent(in) :: qphon(3),rprimd(3,3),xccc1d(n1xccc,6,ntypat) 84 real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom) 85 real(dp),intent(out) :: xccc3d1(cplex*n1*n2*n3) 86 87 !Local variables------------------------------- 88 !scalars 89 integer,parameter :: mshift=401 90 integer :: i1,i2,i3,iatom,ifft,ishift,ishift1,ishift2,ishift3,istr 91 integer :: ixp,jj,ka,kb,mrange,mu,nu 92 real(dp) :: aa,bb,cc,dd,delta,delta2div6,deltam1,diff,difmag 93 real(dp) :: difmag2,difmag2_fact,difmag2_part,func,phase,phi,phr,prod 94 real(dp) :: range,range2,rangem1,rdiff1,rdiff2,rdiff3,term 95 real(dp) :: yy 96 character(len=500) :: message 97 !arrays 98 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 99 integer :: igrid(3),irange(3),ngfft(3) 100 integer,allocatable :: ii(:,:) 101 real(dp) :: drmetds(3,3),lencp(3),rmet(3,3),scale(3),tau(3) 102 real(dp),allocatable :: rrdiff(:,:) 103 104 ! ************************************************************************* 105 106 ! if( ipert<1 .or. ipert> natom+7) then 107 ! write(message,'(a,i0,a,a,a,i0,a)')& 108 !& ' The argument ipert must be between 1 and natom+7=',natom+7,',',ch10,& 109 !& ' while it is ipert=',ipert,'.' 110 ! MSG_BUG(message) 111 ! end if 112 113 if( (ipert==natom+3 .or. ipert==natom+4) .and. cplex/=1) then 114 write(message,'(3a,i4,a)')& 115 & 'The argument cplex must be 1 for strain perturbationh',ch10,& 116 & 'while it is cplex=',cplex,'.' 117 MSG_BUG(message) 118 end if 119 120 !Zero out array 121 xccc3d1(:)=0.0_dp 122 123 !For a non-linear XC core correction, the perturbation must be phonon-type or strain type 124 if(ipert<=natom .or. ipert==natom+3 .or. ipert==natom+4) then 125 126 if( idir<1 .or. idir> 3) then 127 write(message,'(a,a,a,i4,a)')& 128 & 'The argument idir must be between 1 and 3,',ch10,& 129 & 'while it is idir=',idir,'.' 130 MSG_BUG(message) 131 end if 132 133 ! Compute lengths of cross products for pairs of primitive 134 ! translation vectors (used in setting index search range below) 135 lencp(1)=cross_mk(rprimd(1,2),rprimd(2,2),rprimd(3,2),& 136 & rprimd(1,3),rprimd(2,3),rprimd(3,3)) 137 lencp(2)=cross_mk(rprimd(1,3),rprimd(2,3),rprimd(3,3),& 138 & rprimd(1,1),rprimd(2,1),rprimd(3,1)) 139 lencp(3)=cross_mk(rprimd(1,1),rprimd(2,1),rprimd(3,1),& 140 & rprimd(1,2),rprimd(2,2),rprimd(3,2)) 141 142 ! Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3 143 ! (recall ucvol=R1.(R2xR3)) 144 scale(:)=ucvol/lencp(:) 145 146 ! Compute metric tensor in real space rmet 147 do nu=1,3 148 rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+& 149 & rprimd(3,:)*rprimd(3,nu) 150 end do 151 152 ! Section to be executed only for strain perturbation 153 ! Compute derivative of metric tensor wrt strain component istr 154 if(ipert==natom+3 .or. ipert==natom+4) then 155 istr=idir + 3*(ipert-natom-3) 156 157 ka=idx(2*istr-1);kb=idx(2*istr) 158 do jj = 1,3 159 drmetds(:,jj)=(rprimd(ka,:)*rprimd(kb,jj)+rprimd(kb,:)*rprimd(ka,jj)) 160 end do 161 ! For historical reasons: 162 drmetds(:,:)=0.5_dp*drmetds(:,:) 163 164 ! end of strain perturbation section 165 end if 166 167 ngfft(1)=n1 168 ngfft(2)=n2 169 ngfft(3)=n3 170 171 delta=1.0_dp/(n1xccc-1) 172 deltam1=n1xccc-1 173 delta2div6=delta**2/6.0_dp 174 175 ! Loop over atoms in unit cell 176 ! Note that we cycle immediately for all except the displaced atom 177 ! for such a perturbation. The loop is executed over all the 178 ! atoms for a strain peturbation. 179 do iatom=1,natom 180 if(ipert<=natom .and. iatom/=ipert) cycle 181 ! Set search range (density cuts off perfectly beyond range) 182 ! Cycle if no range. 183 range=0.0_dp 184 range=xcccrc(typat(iatom)) 185 if(range<1.d-16) cycle 186 187 range2=range**2 188 rangem1=1.0_dp/range 189 190 ! compute mrange for ii(:,3), rrdiff(:,3), inserted by MM (2005/12/06) 191 ! Consider each component in turn : compute range 192 do mu=1,3 193 194 ! Convert reduced coord of given atom to [0,1) 195 tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp) 196 197 ! Use tau to find nearest grid point along R(mu) 198 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 199 igrid(mu)=nint(tau(mu)*dble(ngfft(mu))) 200 201 ! Use range to compute an index range along R(mu) 202 ! (add 1 to make sure it covers full range) 203 irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu))) 204 205 end do 206 207 ! Allocate arrays that depends on the range 208 mrange=maxval(irange(1:3)) 209 ABI_ALLOCATE(ii,(2*mrange+1,3)) 210 ABI_ALLOCATE(rrdiff,(2*mrange+1,3)) 211 212 ! Consider each component in turn 213 do mu=1,3 214 215 ! temporarily suppressed by MM (2005/12/02) 216 ! Convert reduced coord of given atom to [0,1) 217 ! tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp) 218 219 ! Use tau to find nearest grid point along R(mu) 220 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 221 ! igrid(mu)=nint(tau(mu)*dble(ngfft(mu))) 222 223 ! Use range to compute an index range along R(mu) 224 ! (add 1 to make sure it covers full range) 225 ! irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu))) 226 227 ! Check that the largest range is smallest than the maximum 228 ! allowed one 229 ! if(2*irange(mu)+1 > mshift)then 230 ! write(message, '(a,a,a,a,i6,a)' ) ch10,& 231 ! & ' dfpt_mkcore : BUG -',ch10,& 232 ! & ' The range around atom',iatom,' is too large.' 233 ! MSG_BUG(message) 234 ! end if 235 236 ! Set up a counter that explore the relevant range 237 ! of points around the atom 238 ishift=0 239 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 240 ishift=ishift+1 241 ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu)) 242 rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu) 243 end do 244 245 ! End loop on mu 246 end do 247 248 ! Conduct triple loop over restricted range of grid points for iatom 249 250 do ishift3=1,1+2*irange(3) 251 ! map back to [1,ngfft(3)] for usual fortran index in unit cell 252 i3=ii(ishift3,3) 253 ! find vector from atom location to grid point (reduced) 254 rdiff3=rrdiff(ishift3,3) 255 256 do ishift2=1,1+2*irange(2) 257 i2=ii(ishift2,2) 258 rdiff2=rrdiff(ishift2,2) 259 ! Prepare the computation of difmag2 260 difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2& 261 & +2.0_dp*rmet(3,2)*rdiff3*rdiff2 262 difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2) 263 264 do ishift1=1,1+2*irange(1) 265 rdiff1=rrdiff(ishift1,1) 266 267 ! Compute (rgrid-tau-Rprim)**2 268 difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1) 269 270 ! Only accept contribution inside defined range 271 if (difmag2<range2) then 272 273 ! Prepare computation of core charge function and derivative, 274 ! using splines 275 difmag=sqrt(difmag2) 276 if (difmag>=1.0d-10) then 277 i1=ii(ishift1,1) 278 yy=difmag*rangem1 279 280 ! Compute index of yy over 1 to n1xccc scale 281 jj=1+int(yy*(n1xccc-1)) 282 diff=yy-(jj-1)*delta 283 284 ! Will evaluate spline fit (p. 86 Numerical Recipes, Press et al; 285 ! NOTE error in book for sign of "aa" term in derivative; 286 ! also see splfit routine). 287 bb = diff*deltam1 288 aa = 1.0_dp-bb 289 cc = aa*(aa**2-1.0_dp)*delta2div6 290 dd = bb*(bb**2-1.0_dp)*delta2div6 291 292 ! Evaluate spline fit of 1st der of core charge density 293 ! from xccc1d(:,2,:) and (:,4,:) 294 func=aa*xccc1d(jj,2,typat(iatom))+bb*xccc1d(jj+1,2,typat(iatom)) +& 295 & cc*xccc1d(jj,4,typat(iatom))+dd* xccc1d(jj+1,4,typat(iatom)) 296 297 if(ipert<=natom) then 298 phase=2*pi*(qphon(1)*(rdiff1+xred(1,iatom)) & 299 & +qphon(2)*(rdiff2+xred(2,iatom)) & 300 & +qphon(3)*(rdiff3+xred(3,iatom))) 301 prod=rmet(idir,1)*rdiff1+rmet(idir,2)*rdiff2+rmet(idir,3)*rdiff3 302 303 term=-func*rangem1/difmag*prod 304 ifft=i1+n1*(i2-1+n2*(i3-1)) 305 phr=cos(phase) 306 if(cplex==1)then 307 xccc3d1(ifft)=xccc3d1(ifft)+term*phr 308 else 309 phi=sin(phase) 310 xccc3d1(2*ifft-1)=xccc3d1(2*ifft-1)+term*phr 311 xccc3d1(2*ifft )=xccc3d1(2*ifft )-term*phi 312 end if 313 else 314 prod=& 315 & (rdiff1*(drmetds(1,1)*rdiff1+drmetds(1,2)*rdiff2+drmetds(1,3)*rdiff3)& 316 & +rdiff2*(drmetds(2,1)*rdiff1+drmetds(2,2)*rdiff2+drmetds(2,3)*rdiff3)& 317 & +rdiff3*(drmetds(3,1)*rdiff1+drmetds(3,2)*rdiff2+drmetds(3,3)*rdiff3)) 318 term=prod*func*rangem1/difmag 319 320 ifft=i1+n1*(i2-1+n2*(i3-1)) 321 xccc3d1(ifft)=xccc3d1(ifft)+term 322 323 end if 324 325 ! End of the condition for the distance not to vanish 326 end if 327 328 ! End of condition to be inside the range 329 end if 330 331 ! End loop on ishift1 332 end do 333 334 ! End loop on ishift2 335 end do 336 337 ! End loop on ishift3 338 end do 339 340 ABI_DEALLOCATE(ii) 341 ABI_DEALLOCATE(rrdiff) 342 ! End loop on atoms 343 end do 344 345 ! End of the condition ipert corresponds to a phonon type perturbation 346 ! or strain type perturbation 347 end if 348 349 contains 350 351 function cross_mk(xx,yy,zz,aa,bb,cc) 352 353 354 !This section has been created automatically by the script Abilint (TD). 355 !Do not modify the following lines by hand. 356 #undef ABI_FUNC 357 #define ABI_FUNC 'cross_mk' 358 !End of the abilint section 359 360 real(dp) :: cross_mk 361 real(dp),intent(in) :: xx,yy,zz,aa,bb,cc 362 cross_mk=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2) 363 end function cross_mk 364 365 end subroutine dfpt_mkcore