TABLE OF CONTENTS
- ABINIT/cross_mkcore
- ABINIT/cross_mkcore_alt
- ABINIT/indpos_mkcore_alt
- ABINIT/mkcore
- ABINIT/mkcore_alt
ABINIT/cross_mkcore [ Functions ]
NAME
cross_mkcore
FUNCTION
Define magnitude of cross product of two vectors
PARENTS
CHILDREN
CHILDREN
SOURCE
530 function cross_mkcore(xx,yy,zz,aa,bb,cc) 531 532 533 !This section has been created automatically by the script Abilint (TD). 534 !Do not modify the following lines by hand. 535 #undef ABI_FUNC 536 #define ABI_FUNC 'cross_mkcore' 537 !End of the abilint section 538 539 real(dp) :: cross_mkcore 540 real(dp),intent(in) :: xx,yy,zz,aa,bb,cc 541 ! ************************************************************************* 542 cross_mkcore=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2) 543 end function cross_mkcore 544 545 end subroutine mkcore
ABINIT/cross_mkcore_alt [ Functions ]
NAME
cross_mkcore_alt
FUNCTION
Define magnitude of cross product of two vectors
PARENTS
CHILDREN
SOURCE
1091 function cross_mkcore_alt(xx,yy,zz,aa,bb,cc) 1092 1093 1094 !This section has been created automatically by the script Abilint (TD). 1095 !Do not modify the following lines by hand. 1096 #undef ABI_FUNC 1097 #define ABI_FUNC 'cross_mkcore_alt' 1098 !End of the abilint section 1099 1100 real(dp) :: cross_mkcore_alt 1101 real(dp),intent(in) :: xx,yy,zz,aa,bb,cc 1102 ! ************************************************************************* 1103 cross_mkcore_alt=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2) 1104 end function cross_mkcore_alt
ABINIT/indpos_mkcore_alt [ Functions ]
NAME
indpos_mkcore_alt
FUNCTION
Find the grid index of a given position in the cell according to the BC Determine also whether the index is inside or outside the box for free BC
PARENTS
mkcore
CHILDREN
SOURCE
1124 subroutine indpos_mkcore_alt(periodic,ii,nn,jj,inside) 1125 ! Find the grid index of a given position in the cell according to the BC 1126 ! Determine also whether the index is inside or outside the box for free BC 1127 1128 !This section has been created automatically by the script Abilint (TD). 1129 !Do not modify the following lines by hand. 1130 #undef ABI_FUNC 1131 #define ABI_FUNC 'indpos_mkcore_alt' 1132 !End of the abilint section 1133 1134 integer, intent(in) :: ii,nn 1135 integer, intent(out) :: jj 1136 logical, intent(in) :: periodic 1137 logical, intent(out) :: inside 1138 ! ************************************************************************* 1139 if (periodic) then 1140 inside=.true. ; jj=modulo(ii-1,nn)+1 1141 else 1142 jj=ii ; inside=(ii>=1.and.ii<=nn) 1143 end if 1144 end subroutine indpos_mkcore_alt 1145 1146 end subroutine mkcore_alt
ABINIT/mkcore [ Functions ]
NAME
mkcore
FUNCTION
Optionally compute: (1) pseudo core electron density throughout unit cell (2) pseudo-core contribution to forces (3) pseudo-core contribution to stress tensor (4) pseudo-core contrib. to frozen-wf part the dynamical matrix (part 2)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (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
natom=number of atoms in cell. nfft=(effective) number of FFT grid points (for this processor) nspden=number of spin-density components 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 option: 1 for computing xccc3d (core charge density), 2 for computing core charge contribution to $d(E_{xc})/d(tau)$, 3 for computing core charge contribution to stress tensor corstr, 4 for contribution to frozen-wavefunction part of dynamical matrix rprimd(3,3)=dimensional primitive translation vectors (bohr) typat(natom)=integer type for each atom in cell ucvol=unit cell volume (bohr**3). vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space--only used when option=2,3, or 4, else ignored 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
corstr(6)=core charge contribution to stress tensor, only if option=3 dyfrx2(3,3,natom)=non-linear xc core correction part of the frozen-wavefunction part of the dynamical matrix, only for option=4 grxc(3,natom)=d(Exc)/d(xred), hartree (only computed when option=2, else ignored)
SIDE EFFECTS
xccc3d(n1*n2*n3)=3D core electron density for XC core correction, bohr^-3 (computed and returned when option=1, needed as input when option=3)
NOTES
Note that this routine is tightly connected to the dfpt_mkcore.f routine
PARENTS
dfpt_dyfro,forces,nonlinear,prcref,prcref_PMA,respfn,setvtr,stress xchybrid_ncpp_cc
CHILDREN
SOURCE
62 #if defined HAVE_CONFIG_H 63 #include "config.h" 64 #endif 65 66 #include "abi_common.h" 67 68 subroutine mkcore(corstr,dyfrx2,grxc,mpi_enreg,natom,nfft,nspden,ntypat,n1,n1xccc,& 69 & n2,n3,option,rprimd,typat,ucvol,vxc,xcccrc,xccc1d,xccc3d,xred) 70 71 use defs_basis 72 use defs_abitypes 73 use m_profiling_abi 74 use m_xmpi 75 use m_errors 76 use m_linalg_interfaces 77 78 use m_mpinfo, only : ptabs_fourdp 79 80 !This section has been created automatically by the script Abilint (TD). 81 !Do not modify the following lines by hand. 82 #undef ABI_FUNC 83 #define ABI_FUNC 'mkcore' 84 use interfaces_18_timing 85 use interfaces_41_geometry 86 !End of the abilint section 87 88 implicit none 89 90 !Arguments ------------------------------------ 91 !scalars 92 integer,intent(in) :: n1,n1xccc,n2,n3,natom,nfft,nspden,ntypat,option 93 real(dp),intent(in) :: ucvol 94 type(mpi_type),intent(in) :: mpi_enreg 95 !arrays 96 integer,intent(in) :: typat(natom) 97 real(dp),intent(in) :: rprimd(3,3),vxc(nfft,nspden),xccc1d(n1xccc,6,ntypat) 98 real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom) 99 real(dp),intent(inout) :: xccc3d(nfft) 100 real(dp),intent(out) :: corstr(6),dyfrx2(3,3,natom) 101 real(dp),intent(inout) :: grxc(3,natom) 102 103 !Local variables------------------------------- 104 !scalars 105 integer :: i1,i2,i3,iatom,ier,ifft,ishift,ishift1,ishift2 106 integer :: ishift3,itypat,ixp,jj,me_fft,mrange,mu,nfftot,nu 107 real(dp) :: dd,delta,delta2div6,deltam1,diff,difmag,difmag2 108 real(dp) :: difmag2_fact,difmag2_part,fact,func,grxc1,grxc2,grxc3,range,range2 109 real(dp) :: rangem1,rdiff1,rdiff2,rdiff3,strdia,t1,t2,t3,term,term1,term2 110 character(len=500) :: message 111 !arrays 112 integer :: igrid(3),irange(3),ngfft(3) 113 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 114 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 115 integer,allocatable :: ii(:,:) 116 real(dp) :: yy,aa,bb,cc 117 real(dp) :: corfra(3,3),lencp(3),rmet(3,3),scale(3),tau(3),tsec(2),tt(3) 118 real(dp),allocatable :: rrdiff(:,:),work(:,:,:) 119 120 !************************************************************************ 121 122 call timab(12,1,tsec) 123 124 !Make sure option is acceptable 125 if (option<0.or.option>4) then 126 write(message, '(a,i12,a,a,a)' )& 127 & 'option=',option,' is not allowed.',ch10,& 128 & 'Must be 1, 2, 3 or 4.' 129 MSG_BUG(message) 130 end if 131 132 !Zero out only the appropriate array according to option: 133 !others are dummies with no storage 134 135 if (option==1) then 136 ! Zero out array to permit accumulation over atom types below: 137 xccc3d(:)=zero 138 else if (option==2) then 139 ! Zero out gradient of Exc array 140 grxc(:,:)=zero 141 else if (option==3) then 142 ! Zero out locally defined stress array 143 corfra(:,:)=zero 144 strdia=zero 145 else if (option==4) then 146 ! Zero out fr-wf part of the dynamical matrix 147 dyfrx2(:,:,:)=zero 148 else 149 MSG_BUG(" Can't be here! (bad option)") 150 end if 151 152 !Compute lengths of cross products for pairs of primitive 153 !translation vectors (used in setting index search range below) 154 lencp(1)=cross_mkcore(rprimd(1,2),rprimd(2,2),rprimd(3,2),& 155 & rprimd(1,3),rprimd(2,3),rprimd(3,3)) 156 lencp(2)=cross_mkcore(rprimd(1,3),rprimd(2,3),rprimd(3,3),& 157 & rprimd(1,1),rprimd(2,1),rprimd(3,1)) 158 lencp(3)=cross_mkcore(rprimd(1,1),rprimd(2,1),rprimd(3,1),& 159 & rprimd(1,2),rprimd(2,2),rprimd(3,2)) 160 161 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3 162 !(recall ucvol=R1.(R2xR3)) 163 scale(:)=ucvol/lencp(:) 164 165 !Compute metric tensor in real space rmet 166 do nu=1,3 167 rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu) 168 end do 169 170 ngfft(1)=n1 171 ngfft(2)=n2 172 ngfft(3)=n3 173 nfftot=n1*n2*n3 174 me_fft = mpi_enreg%me_fft 175 176 ! Get the distrib associated with this fft_grid 177 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 178 179 delta=one/(n1xccc-1) 180 deltam1=n1xccc-1 181 delta2div6=delta**2/6.0d0 182 183 if (option>=2) then 184 ABI_ALLOCATE(work,(n1,n2,n3)) 185 ! For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down)) 186 ! For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22}) 187 if (nspden>=2) then 188 ifft=1 189 do i3=1,n3 190 if(me_fft==fftn3_distrib(i3)) then 191 do i2=1,n2 192 do i1=1,n1 193 work(i1,i2,i3)=half*(vxc(ifft,1)+vxc(ifft,2)) 194 ifft=ifft+1 195 end do 196 end do 197 end if 198 end do 199 else 200 ifft=1 201 do i3=1,n3 202 if(me_fft==fftn3_distrib(i3)) then 203 do i2=1,n2 204 do i1=1,n1 205 work(i1,i2,i3)=vxc(ifft,1) 206 ifft=ifft+1 207 end do 208 end do 209 end if 210 end do 211 ! call DCOPY(nfft,vxc,1,work,1) 212 end if 213 end if 214 215 !Loop over atoms in unit cell 216 do iatom=1,natom 217 218 if(option==2)then 219 grxc1=zero 220 grxc2=zero 221 grxc3=zero 222 end if 223 224 ! Set search range (density cuts off perfectly beyond range) 225 itypat=typat(iatom) 226 range=xcccrc(itypat) 227 228 ! Skip loop if this atom has no core charge 229 if (abs(range)<1.d-16) cycle 230 231 range2=range**2 232 rangem1=one/range 233 234 ! Consider each component in turn : compute range 235 do mu=1,3 236 237 ! Convert reduced coord of given atom to [0,1) 238 tau(mu)=mod(xred(mu,iatom)+one-aint(xred(mu,iatom)),one) 239 240 ! Use tau to find nearest grid point along R(mu) 241 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 242 igrid(mu)=nint(tau(mu)*dble(ngfft(mu))) 243 244 ! Use range to compute an index range along R(mu) 245 ! (add 1 to make sure it covers full range) 246 irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu))) 247 248 end do 249 250 ! Allocate arrays that depends on the range 251 mrange=maxval(irange(1:3)) 252 ABI_ALLOCATE(ii,(2*mrange+1,3)) 253 ABI_ALLOCATE(rrdiff,(2*mrange+1,3)) 254 255 ! Set up counters that explore the relevant range 256 ! of points around the atom 257 do mu=1,3 258 ishift=0 259 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 260 ishift=ishift+1 261 ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu)) 262 rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu) 263 end do 264 end do 265 266 ! Conduct triple loop over restricted range of grid points for iatom 267 do ishift3=1,1+2*irange(3) 268 ! map back to [1,ngfft(3)] for usual fortran index in unit cell 269 i3=ii(ishift3,3) 270 if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle 271 ! find vector from atom location to grid point (reduced) 272 rdiff3=rrdiff(ishift3,3) 273 274 do ishift2=1,1+2*irange(2) 275 i2=ii(ishift2,2) 276 rdiff2=rrdiff(ishift2,2) 277 ! Prepare the computation of difmag2 278 difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2& 279 & +2.0d0*rmet(3,2)*rdiff3*rdiff2 280 difmag2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2) 281 282 do ishift1=1,1+2*irange(1) 283 rdiff1=rrdiff(ishift1,1) 284 285 ! Compute (rgrid-tau-Rprim)**2 286 difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1) 287 288 ! Only accept contribution inside defined range 289 if (difmag2<range2-tol12) then 290 291 ! Prepare computation of core charge function and derivative, 292 ! using splines 293 i1=ii(ishift1,1) 294 difmag=sqrt(difmag2) 295 yy=difmag*rangem1 296 297 ! Compute index of yy over 1 to n1xccc scale 298 jj=1+int(yy*(n1xccc-1)) 299 diff=yy-(jj-1)*delta 300 301 ! Will evaluate spline fit (p. 86 Numerical Recipes, Press et al; 302 ! NOTE error in book for sign of "aa" term in derivative; 303 ! also see splfit routine). 304 bb = diff*deltam1 305 aa = one-bb 306 cc = aa*(aa**2-one)*delta2div6 307 dd = bb*(bb**2-one)*delta2div6 308 309 310 ! Test first for option 2, the most frequently used 311 if (option==2) then 312 313 ! Accumulate contributions to Exc gradients 314 315 if (difmag>1.0d-10) then 316 317 ! Evaluate spline fit of 1st der of core charge density 318 ! from xccc1d(:,2,:) and (:,4,:) 319 func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 320 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 321 term=work(i1,i2,i3)*func/difmag 322 grxc1=grxc1+rdiff1*term 323 grxc2=grxc2+rdiff2*term 324 grxc3=grxc3+rdiff3*term 325 end if 326 327 else if (option==1) then 328 329 ! Evaluate spline fit of core charge density 330 ! from xccc1d(:,1,:) and (:,3,:) 331 func=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +& 332 & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat) 333 334 ! Accumulate contributions to core electron density 335 ! throughout unit cell 336 ifft=i1+n1*(i2-1+n2*(ffti3_local(i3)-1)) 337 xccc3d(ifft)=xccc3d(ifft)+func 338 339 else if (option==3) then 340 341 ! Accumulate contributions to stress tensor 342 ! in reduced coordinates 343 344 if (difmag>1.0d-10) then 345 346 ! Evaluate spline fit of 1st der of core charge density 347 ! from xccc1d(:,2,:) and (:,4,:) 348 func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 349 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 350 term=work(i1,i2,i3)*func*rangem1/difmag/dble(n1*n2*n3) 351 ! Write out the 6 symmetric components 352 corfra(1,1)=corfra(1,1)+term*rdiff1**2 353 corfra(2,2)=corfra(2,2)+term*rdiff2**2 354 corfra(3,3)=corfra(3,3)+term*rdiff3**2 355 corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2 356 corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1 357 corfra(2,1)=corfra(2,1)+term*rdiff2*rdiff1 358 ! (the above still needs to be transformed to cartesian coords) 359 360 end if 361 362 ! Also compute a diagonal term 363 ! Evaluate spline fit of core charge density 364 ! from xccc1d(:,1,:) and (:,3,:) 365 func=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +& 366 & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat) 367 strdia=strdia+work(i1,i2,i3)*func 368 369 else if (option==4) then 370 371 ! Compute frozen-wf contribution to Dynamical matrix 372 373 tt(1)=rmet(1,1)*rdiff1+rmet(1,2)*rdiff2+rmet(1,3)*rdiff3 374 tt(2)=rmet(2,1)*rdiff1+rmet(2,2)*rdiff2+rmet(2,3)*rdiff3 375 tt(3)=rmet(3,1)*rdiff1+rmet(3,2)*rdiff2+rmet(3,3)*rdiff3 376 377 if (difmag>1.d-10) then 378 379 ! Accumulate contributions to dynamical matrix 380 term=(ucvol/dble(nfftot))*work(i1,i2,i3)*rangem1/difmag 381 ! Evaluate spline fit of 1st der of core charge density 382 ! from xccc1d(:,2,:) and (:,4,:) 383 func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 384 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 385 term1=term*func 386 ! Evaluate spline fit of 2nd der of core charge density 387 ! from xccc1d(:,3,:) and (:,5,:) 388 func=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +& 389 & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat) 390 term2=term*func*rangem1/difmag 391 do mu=1,3 392 do nu=1,3 393 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)& 394 & +(term2-term1/difmag**2)*tt(mu)*tt(nu)& 395 & +term1*rmet(mu,nu) 396 end do 397 end do 398 399 else 400 401 ! There is a contribution from difmag=zero ! 402 ! Evaluate spline fit of 2nd der of core charge density 403 ! from xccc1d(:,3,:) and (:,5,:) 404 func=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +& 405 & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat) 406 term=(ucvol/dble(nfftot))*work(i1,i2,i3)*func*rangem1**2 407 do mu=1,3 408 do nu=1,3 409 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu) 410 end do 411 end do 412 413 ! End of condition not to be precisely on the point (difmag=zero) 414 end if 415 416 ! If option is not 1, 2, 3, or 4. 417 else 418 MSG_BUG("Can't be here in mkcore") 419 ! End of choice of option 420 end if 421 422 ! End of condition on the range 423 end if 424 425 ! End loop on ishift1 426 end do 427 428 ! End loop on ishift2 429 end do 430 431 ! End loop on ishift3 432 end do 433 434 ABI_DEALLOCATE(ii) 435 ABI_DEALLOCATE(rrdiff) 436 437 if(option==2)then 438 fact=-(ucvol/dble(nfftot))/range 439 grxc(1,iatom)=grxc1*fact 440 grxc(2,iatom)=grxc2*fact 441 grxc(3,iatom)=grxc3*fact 442 end if 443 444 ! End big loop on atoms 445 end do 446 447 if (option==2) then 448 449 ! Apply rmet as needed to get reduced coordinate gradients 450 do iatom=1,natom 451 t1=grxc(1,iatom) 452 t2=grxc(2,iatom) 453 t3=grxc(3,iatom) 454 grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3 455 456 end do 457 end if 458 459 if (option==3) then 460 461 ! Transform stress tensor from full storage mode to symmetric storage mode 462 corstr(1)=corfra(1,1) 463 corstr(2)=corfra(2,2) 464 corstr(3)=corfra(3,3) 465 corstr(4)=corfra(3,2) 466 corstr(5)=corfra(3,1) 467 corstr(6)=corfra(2,1) 468 469 ! Transform stress tensor from reduced coordinates to cartesian coordinates 470 call strconv(corstr,rprimd,corstr) 471 472 ! Compute diagonal contribution to stress tensor (need input xccc3d) 473 ! strdia = (1/N) Sum(r) [mu_xc_avg(r) * rho_core(r)] 474 ifft=0 ; strdia=zero 475 do i3=1,n3 476 if(me_fft==fftn3_distrib(i3)) then 477 do i2=1,n2 478 do i1=1,n1 479 ifft=ifft+1 480 strdia=strdia+work(i1,i2,i3)*xccc3d(ifft) 481 end do 482 end do 483 end if 484 end do 485 strdia=strdia/dble(nfftot) 486 ! strdia=DDOT(nfft,work,1,xccc3d,1)/dble(nfftot) 487 488 ! Add diagonal term to stress tensor 489 corstr(1)=corstr(1)+strdia 490 corstr(2)=corstr(2)+strdia 491 corstr(3)=corstr(3)+strdia 492 end if 493 494 if(option>=2) then 495 ABI_DEALLOCATE(work) 496 end if 497 498 if(mpi_enreg%nproc_fft > 1) then 499 call timab(539,1,tsec) 500 if(option==2) then 501 call xmpi_sum(grxc,mpi_enreg%comm_fft,ier) 502 end if 503 if(option==3) then 504 call xmpi_sum(corstr,mpi_enreg%comm_fft,ier) 505 end if 506 if(option==4) then 507 call xmpi_sum(dyfrx2,mpi_enreg%comm_fft,ier) 508 end if 509 call timab(539,2,tsec) 510 end if 511 512 call timab(12,2,tsec) 513 514 contains
ABINIT/mkcore_alt [ Functions ]
NAME
mkcore_alt
FUNCTION
Optionally compute: (1) pseudo core electron density throughout unit cell (2) pseudo-core contribution to forces (3) pseudo-core contribution to stress tensor (4) pseudo-core contrib. to frozen-wf part the dynamical matrix (part 2) This routine is an alternative to mkcore, to be used for PAW and/or WVL.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (TRangel,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 .
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx
icoulomb= periodic treatment of Hartree potential: 0=periodic, 1=free BC, 2=surface BC
mpi_enreg=informations about MPI parallelization
natom=number of atoms in cell.
nfft=(effective) number of FFT grid points (for this processor)
nspden=number of spin-density components
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
option: 1 for computing core charge density
2 for computing core charge contribution to forces
3 for computing core charge contribution to stress tensor
4 for computing contribution to frozen-wf part of dynamical matrix
pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
rprimd(3,3)=dimensional primitive translation vectors (bohr)
ucvol=unit cell volume (bohr**3)
usepaw=flag for PAW method
vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space
xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
xccc1d(n1xccc,6,ntypat)=1D core charge function and 5 derivatives for each atom type
xred(3,natom)=reduced coordinates for atoms in unit cell
OUTPUT
=== if option==1 === xccc3d(n1*n2*n3)=3D core electron density for XC core correction (bohr^-3) === if option==2 === grxc(3,natom)=core charge contribution to forces === if option==3 === corstr(6)=core charge contribution to stress tensor === if option==4 === dyfrx2(3,3,natom)=non-linear xc core correction part of the frozen-wavefunction part of the dynamical matrix
SIDE EFFECTS
xccc3d(n1*n2*n3)=3D core electron density for XC core correction (bohr^-3) (computed and returned when option=1, needed as input when option=3)
NOTES
Based on mkcore.F90
PARENTS
forces,setvtr,stress
CHILDREN
SOURCE
617 #if defined HAVE_CONFIG_H 618 #include "config.h" 619 #endif 620 621 #include "abi_common.h" 622 623 subroutine mkcore_alt(atindx1,corstr,dyfrx2,grxc,icoulomb,mpi_enreg,natom,nfft,nspden,& 624 & nattyp,ntypat,n1,n1xccc,n2,n3,option,rprimd,ucvol,vxc,xcccrc,xccc1d,& 625 & xccc3d,xred,pawrad,pawtab,usepaw) 626 627 use defs_basis 628 use defs_abitypes 629 use m_profiling_abi 630 use m_xmpi 631 use m_errors 632 use m_linalg_interfaces 633 634 use m_mpinfo, only : ptabs_fourdp 635 use m_sort, only : sort_dp 636 use m_pawrad, only : pawrad_type,pawrad_init,pawrad_free 637 use m_pawtab, only : pawtab_type 638 use m_paw_numeric, only : paw_splint 639 640 !This section has been created automatically by the script Abilint (TD). 641 !Do not modify the following lines by hand. 642 #undef ABI_FUNC 643 #define ABI_FUNC 'mkcore_alt' 644 use interfaces_18_timing 645 use interfaces_41_geometry 646 !End of the abilint section 647 648 implicit none 649 650 !Arguments ------------------------------------ 651 !scalars 652 integer,intent(in) :: icoulomb,n1,n1xccc,n2,n3,natom,nfft,nspden,ntypat,option,usepaw 653 real(dp),intent(in) :: ucvol 654 type(mpi_type),intent(in) :: mpi_enreg 655 type(pawrad_type),intent(in) :: pawrad(:) 656 type(pawtab_type),intent(in) :: pawtab(:) 657 !arrays 658 integer,intent(in) :: atindx1(natom),nattyp(ntypat) 659 real(dp),intent(in) :: rprimd(3,3),xccc1d(n1xccc,6,ntypat) 660 real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom) 661 real(dp),intent(in),target :: vxc(nfft,nspden) 662 real(dp),intent(out) :: corstr(6),grxc(3,natom),dyfrx2(3,3,natom) 663 real(dp),intent(inout) :: xccc3d(nfft) 664 665 !Local variables------------------------------- 666 !scalars 667 integer :: i1,i2,i3,iat,iatm,iatom,ier,ipts 668 integer :: ishift,ishift1,ishift2,ishift3 669 integer :: itypat,ixp,jj,jpts,me_fft,mrange,mu 670 integer :: nfftot,npts,npts12,nu 671 logical :: letsgo 672 real(dp) :: aa,bb,cc,dd,delta,delta2div6,deltam1 673 real(dp) :: diff,difmag,fact,range,range2 674 real(dp) :: rangem1,rdiff1,rdiff2,rdiff3 675 real(dp) :: rnorm2,rnorm2_fact,rnorm2_part 676 real(dp) :: strdia,t1,t2,t3,term,term1,term2,yy 677 character(len=1) :: geocode 678 character(len=500) :: message 679 type(pawrad_type) :: core_mesh 680 !arrays 681 integer :: igrid(3),irange(3),ishiftmax(3),ngfft(3) 682 integer,allocatable :: ii(:,:),iindex(:),indx1(:),indx2(:) 683 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 684 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 685 logical :: per(3) 686 real(dp) :: corfra(3,3),corgr(3),lencp(3),rmet(3,3) 687 real(dp) :: scale(3),tau(3),tsec(2),tt(3) 688 real(dp),allocatable :: dtcore(:),d2tcore(:),rnorm(:) 689 real(dp),allocatable :: rrdiff(:,:),tcore(:) 690 real(dp), ABI_CONTIGUOUS pointer :: vxc_eff(:) 691 692 !************************************************************************ 693 694 call timab(12,1,tsec) 695 696 !Make sure option is acceptable 697 if (option<0.or.option>4) then 698 write(message, '(a,i12,a,a,a)' )& 699 & 'option=',option,' is not allowed.',ch10,& 700 & 'Must be 1, 2, 3 or 4.' 701 MSG_BUG(message) 702 end if 703 704 !Zero out only the appropriate array according to option: 705 if (option==1) then 706 xccc3d(:)=zero 707 else if (option==2) then 708 grxc(:,:)=zero 709 else if (option==3) then 710 corfra(:,:)=zero 711 strdia=zero 712 else if (option==4) then 713 dyfrx2(:,:,:)=zero 714 end if 715 716 !Conditions for periodicity in the three directions 717 geocode='P' 718 if (icoulomb==1) geocode='F' 719 if (icoulomb==2) geocode='S' 720 per(1)=(geocode /= 'F') 721 per(2)=(geocode == 'P') 722 per(3)=(geocode /= 'F') 723 724 !Compute lengths of cross products for pairs of primitive 725 !translation vectors (used in setting index search range below) 726 lencp(1)=cross_mkcore_alt(rprimd(1,2),rprimd(2,2),rprimd(3,2),& 727 & rprimd(1,3),rprimd(2,3),rprimd(3,3)) 728 lencp(2)=cross_mkcore_alt(rprimd(1,3),rprimd(2,3),rprimd(3,3),& 729 & rprimd(1,1),rprimd(2,1),rprimd(3,1)) 730 lencp(3)=cross_mkcore_alt(rprimd(1,1),rprimd(2,1),rprimd(3,1),& 731 & rprimd(1,2),rprimd(2,2),rprimd(3,2)) 732 733 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3 734 !(recall ucvol=R1.(R2xR3)) 735 scale(:)=ucvol/lencp(:) 736 737 !Compute metric tensor in real space rmet 738 do nu=1,3 739 rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu) 740 end do 741 742 !Get the distrib associated with this fft_grid 743 ngfft(1)=n1;ngfft(2)=n2;ngfft(3)=n3 744 nfftot=n1*n2*n3 ; me_fft=mpi_enreg%me_fft 745 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 746 747 delta=one/(n1xccc-1) 748 deltam1=n1xccc-1 749 delta2div6=delta**2/6.0_dp 750 751 if (option>=2) then 752 ! For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down)) 753 ! For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22}) 754 if (nspden>=2) then 755 ABI_ALLOCATE(vxc_eff,(nfft)) 756 do jj=1,nfft 757 vxc_eff(jj)=half*(vxc(jj,1)+vxc(jj,2)) 758 end do 759 else 760 vxc_eff => vxc(1:nfft,1) 761 end if 762 end if 763 764 !Loop over atom types 765 iatm=0 766 do itypat=1,ntypat 767 768 ! Set search range (density cuts off perfectly beyond range) 769 range=xcccrc(itypat);if (usepaw==1) range=pawtab(itypat)%rcore 770 range2=range**2 ; rangem1=one/range 771 772 ! Skip loop if this type has no core charge 773 if (abs(range)<1.d-16) cycle 774 775 ! PAW: create mesh for core density 776 if (usepaw==1) then 777 call pawrad_init(core_mesh,mesh_size=pawtab(itypat)%core_mesh_size,& 778 & mesh_type=pawrad(itypat)%mesh_type,& 779 & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep) 780 end if 781 782 ! Loop over atoms of the type 783 do iat=1,nattyp(itypat) 784 iatm=iatm+1;iatom=atindx1(iatm) 785 786 if(option==2) corgr(:)=zero 787 788 ! Consider each component in turn : compute range 789 do mu=1,3 790 ! Convert reduced coord of given atom to [0,1) 791 tau(mu)=mod(xred(mu,iatom)+one-aint(xred(mu,iatom)),one) 792 ! Use tau to find nearest grid point along R(mu) 793 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 794 igrid(mu)=nint(tau(mu)*real(ngfft(mu),dp)) 795 ! Use range to compute an index range along R(mu) 796 ! (add 1 to make sure it covers full range) 797 irange(mu)=1+nint((range/scale(mu))*real(ngfft(mu),dp)) 798 end do 799 800 ! Allocate arrays that depends on the range 801 mrange=maxval(irange(1:3)) 802 ABI_ALLOCATE(ii,(2*mrange+1,3)) 803 ABI_ALLOCATE(rrdiff,(2*mrange+1,3)) 804 805 ! Set up counters that explore the relevant range of points around the atom 806 if (geocode=='P') then 807 ! Fully periodic version 808 do mu=1,3 809 ishift=0 810 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 811 ishift=ishift+1 812 ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu)) 813 rrdiff(ishift,mu)=real(ixp,dp)/real(ngfft(mu),dp)-tau(mu) 814 end do 815 ishiftmax(mu)=ishift 816 end do 817 else 818 ! Free or surface conditions 819 do mu=1,3 820 ishift=0 821 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 822 call indpos_mkcore_alt(per(mu),ixp,ngfft(mu),jj,letsgo) 823 if (letsgo) then 824 ishift=ishift+1;ii(ishift,mu)=1+jj 825 rrdiff(ishift,mu)=real(ixp,dp)/real(ngfft(mu),dp)-tau(mu) 826 end if 827 end do 828 ishiftmax(mu)=ishift 829 end do 830 end if 831 npts12=ishiftmax(1)*ishiftmax(2) 832 ABI_ALLOCATE(indx1,(npts12)) 833 ABI_ALLOCATE(indx2,(npts12)) 834 ABI_ALLOCATE(iindex,(npts12)) 835 ABI_ALLOCATE(rnorm,(npts12)) 836 if (option==1.or.option==3) then 837 ABI_ALLOCATE(tcore,(npts12)) 838 end if 839 if (option>=2) then 840 ABI_ALLOCATE(dtcore,(npts12)) 841 end if 842 if (option==4) then 843 ABI_ALLOCATE(d2tcore,(npts12)) 844 end if 845 846 ! Conduct loop over restricted range of grid points for iatom 847 do ishift3=1,ishiftmax(3) 848 i3=ii(ishift3,3) 849 rdiff3=rrdiff(ishift3,3) 850 851 if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle 852 853 ! Select the vectors located around the current atom 854 ! TR: all of the following could be done inside or 855 ! outside the loops (i2,i1,i3). 856 ! Outside: the memory consumption increases. 857 ! Inside: the time of calculation increases. 858 ! Here, I choose to do it here, somewhere in the middle. 859 npts=0 860 do ishift2=1,ishiftmax(2) 861 i2=ii(ishift2,2) ; rdiff2=rrdiff(ishift2,2) 862 rnorm2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2 & 863 & +2.0d0*rmet(3,2)*rdiff3*rdiff2 864 rnorm2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2) 865 do ishift1=1,ishiftmax(1) 866 i1=ii(ishift1,1) ; rdiff1=rrdiff(ishift1,1) 867 rnorm2=rnorm2_part+rdiff1*(rnorm2_fact+rmet(1,1)*rdiff1) 868 ! Only accept contributions inside defined range 869 if (rnorm2<range2-tol12) then 870 npts=npts+1 ; iindex(npts)=npts 871 indx1(npts)=ishift1;indx2(npts)=ishift2 872 rnorm(npts)=sqrt(rnorm2) 873 end if 874 end do 875 end do 876 if (npts==0) cycle 877 if (npts>npts12) then 878 message='npts>npts12!' 879 MSG_BUG(message) 880 end if 881 882 ! Evaluate core density (and derivatives) on the set of selected points 883 if (usepaw==1) then 884 ! PAW: use splint routine 885 call sort_dp(npts,rnorm(1:npts),iindex(1:npts),tol16) 886 if (option==1.or.option==3) then 887 ! Evaluate fit of core density 888 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 889 & pawtab(itypat)%tcoredens(:,1), & 890 & pawtab(itypat)%tcoredens(:,3),& 891 & npts,rnorm(1:npts),tcore(1:npts)) 892 end if 893 if (option>=2) then 894 ! Evaluate fit of 1-der of core density 895 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 896 & pawtab(itypat)%tcoredens(:,2), & 897 & pawtab(itypat)%tcoredens(:,4),& 898 & npts,rnorm(1:npts),dtcore(1:npts)) 899 end if 900 if (option==4) then 901 ! Evaluate fit of 2nd-der of core density 902 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 903 & pawtab(itypat)%tcoredens(:,3), & 904 & pawtab(itypat)%tcoredens(:,5),& 905 & npts,rnorm(1:npts),d2tcore(1:npts)) 906 end if 907 else 908 ! Norm-conserving PP: 909 ! Evaluate spline fit with method from Numerical Recipes 910 ! (p. 86 Numerical Recipes, Press et al; 911 ! NOTE error in book for sign of "aa" term in derivative) 912 do ipts=1,npts 913 yy=rnorm(ipts)*rangem1 914 jj=1+int(yy*(n1xccc-1)) 915 diff=yy-(jj-1)*delta 916 bb = diff*deltam1 ; aa = one-bb 917 cc = aa*(aa**2-one)*delta2div6 918 dd = bb*(bb**2-one)*delta2div6 919 if (option==1.or.option==3) then 920 tcore(ipts)=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +& 921 & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat) 922 end if 923 if (option>=2) then 924 dtcore(ipts)=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 925 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 926 end if 927 if (option==4) then 928 d2tcore(ipts)=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +& 929 & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat) 930 end if 931 end do 932 end if 933 934 ! Now, perform the loop over selected grid points 935 do ipts=1,npts 936 ishift1=indx1(iindex(ipts)) 937 ishift2=indx2(iindex(ipts)) 938 difmag=rnorm(ipts) 939 940 rdiff1=rrdiff(ishift1,1);rdiff2=rrdiff(ishift2,2) 941 jpts=ii(ishift1,1)+n1*(ii(ishift2,2)-1+n2*(ffti3_local(i3)-1)) 942 943 ! === Evaluate charge density 944 if (option==1) then 945 xccc3d(jpts)=xccc3d(jpts)+tcore(ipts) 946 947 ! === Accumulate contributions to forces 948 else if (option==2) then 949 if (difmag>tol10) then 950 term=vxc_eff(jpts)*dtcore(ipts)/difmag 951 corgr(1)=corgr(1)+rdiff1*term 952 corgr(2)=corgr(2)+rdiff2*term 953 corgr(3)=corgr(3)+rdiff3*term 954 end if 955 956 ! === Accumulate contributions to stress tensor (in red. coordinates) 957 else if (option==3) then 958 if (difmag>tol10) then 959 term=vxc_eff(jpts)*dtcore(ipts)*rangem1/difmag/real(nfftot,dp) 960 ! Write out the 6 symmetric components 961 corfra(1,1)=corfra(1,1)+term*rdiff1*rdiff1 962 corfra(2,2)=corfra(2,2)+term*rdiff2*rdiff2 963 corfra(3,3)=corfra(3,3)+term*rdiff3*rdiff3 964 corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2 965 corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1 966 corfra(2,1)=corfra(2,1)+term*rdiff2*rdiff1 967 ! (the above still needs to be transformed to cartesian coords) 968 end if 969 ! Also compute a diagonal term 970 strdia=strdia+vxc_eff(jpts)*tcore(ipts) 971 972 ! === Compute frozen-wf contribution to Dynamical matrix 973 else if (option==4) then 974 tt(1)=rmet(1,1)*rdiff1+rmet(1,2)*rdiff2+rmet(1,3)*rdiff3 975 tt(2)=rmet(2,1)*rdiff1+rmet(2,2)*rdiff2+rmet(2,3)*rdiff3 976 tt(3)=rmet(3,1)*rdiff1+rmet(3,2)*rdiff2+rmet(3,3)*rdiff3 977 if (difmag>tol10) then 978 term=(ucvol/real(nfftot,dp))*vxc_eff(jpts)*rangem1/difmag 979 term1=term*tcore(ipts) 980 term2=term*d2tcore(ipts)*rangem1/difmag 981 do mu=1,3 982 do nu=1,3 983 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)& 984 & +(term2-term1/difmag**2)*tt(mu)*tt(nu)& 985 & +term1*rmet(mu,nu) 986 end do 987 end do 988 else 989 ! There is a contribution from difmag=zero ! 990 term=(ucvol/real(nfftot,dp))*vxc_eff(jpts)*d2tcore(ipts)*rangem1**2 991 do mu=1,3 992 do nu=1,3 993 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu) 994 end do 995 end do 996 end if 997 end if ! Choice of option 998 999 end do ! Loop on ipts (ishift1, ishift2) 1000 1001 end do ! Loop on ishift3 1002 1003 ABI_DEALLOCATE(ii) 1004 ABI_DEALLOCATE(rrdiff) 1005 ABI_DEALLOCATE(indx1) 1006 ABI_DEALLOCATE(indx2) 1007 ABI_DEALLOCATE(iindex) 1008 ABI_DEALLOCATE(rnorm) 1009 if (allocated(tcore)) then 1010 ABI_DEALLOCATE(tcore) 1011 end if 1012 if (allocated(dtcore)) then 1013 ABI_DEALLOCATE(dtcore) 1014 end if 1015 if (allocated(d2tcore)) then 1016 ABI_DEALLOCATE(d2tcore) 1017 end if 1018 1019 if (option==2) then 1020 fact=-(ucvol/real(nfftot,dp))/range 1021 grxc(:,iatom)=corgr(:)*fact 1022 end if 1023 1024 ! End loop on atoms 1025 end do 1026 1027 if (usepaw==1) then 1028 call pawrad_free(core_mesh) 1029 end if 1030 1031 !End loop over atom types 1032 end do 1033 1034 if(option>=2.and.nspden>=2) then 1035 ABI_DEALLOCATE(vxc_eff) 1036 end if 1037 1038 !Forces: translate into reduced coordinates 1039 if (option==2) then 1040 do iatom=1,natom 1041 t1=grxc(1,iatom);t2=grxc(2,iatom);t3=grxc(3,iatom) 1042 grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3 1043 end do 1044 end if 1045 1046 !Stress tensor: symmetrize, translate into cartesian coord., add diagonal part 1047 if (option==3) then 1048 corstr(1)=corfra(1,1) ; corstr(2)=corfra(2,2) 1049 corstr(3)=corfra(3,3) ; corstr(4)=corfra(3,2) 1050 corstr(5)=corfra(3,1) ; corstr(6)=corfra(2,1) 1051 call strconv(corstr,rprimd,corstr) 1052 corstr(1)=corstr(1)+strdia/real(nfftot,dp) 1053 corstr(2)=corstr(2)+strdia/real(nfftot,dp) 1054 corstr(3)=corstr(3)+strdia/real(nfftot,dp) 1055 end if 1056 1057 !If needed sum over MPI processes 1058 if(mpi_enreg%nproc_fft>1) then 1059 call timab(539,1,tsec) 1060 if (option==2) then 1061 call xmpi_sum(grxc,mpi_enreg%comm_fft,ier) 1062 end if 1063 if (option==3) then 1064 call xmpi_sum(corstr,mpi_enreg%comm_fft,ier) 1065 end if 1066 if (option==4) then 1067 call xmpi_sum(dyfrx2,mpi_enreg%comm_fft,ier) 1068 end if 1069 call timab(539,2,tsec) 1070 end if 1071 1072 call timab(12,2,tsec) 1073 1074 contains