TABLE OF CONTENTS
- ABINIT/cross_mkcore
- ABINIT/cross_mkcore_alt
- ABINIT/dfpt_mkcore
- ABINIT/indpos_mkcore_alt
- ABINIT/m_mkcore
- ABINIT/mkcore
- ABINIT/mkcore_alt
ABINIT/cross_mkcore [ Functions ]
NAME
cross_mkcore
FUNCTION
Define magnitude of cross product of two vectors
SOURCE
531 function cross_mkcore(xx,yy,zz,aa,bb,cc) 532 533 real(dp) :: cross_mkcore 534 real(dp),intent(in) :: xx,yy,zz,aa,bb,cc 535 ! ************************************************************************* 536 cross_mkcore=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2) 537 end function cross_mkcore 538 539 end subroutine mkcore
ABINIT/cross_mkcore_alt [ Functions ]
NAME
cross_mkcore_alt
FUNCTION
Define magnitude of cross product of two vectors
SOURCE
1061 function cross_mkcore_alt(xx,yy,zz,aa,bb,cc) 1062 1063 real(dp) :: cross_mkcore_alt 1064 real(dp),intent(in) :: xx,yy,zz,aa,bb,cc 1065 ! ************************************************************************* 1066 cross_mkcore_alt=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2) 1067 end function cross_mkcore_alt
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.
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
SOURCE
1141 subroutine dfpt_mkcore(cplex,idir,ipert,natom,ntypat,n1,n1xccc,& 1142 & n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xccc3d1,xred) 1143 1144 !Arguments ------------------------------------ 1145 !scalars 1146 integer,intent(in) :: cplex,idir,ipert,n1,n1xccc,n2,n3,natom,ntypat 1147 real(dp),intent(in) :: ucvol 1148 !arrays 1149 integer,intent(in) :: typat(natom) 1150 real(dp),intent(in) :: qphon(3),rprimd(3,3),xccc1d(n1xccc,6,ntypat) 1151 real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom) 1152 real(dp),intent(out) :: xccc3d1(cplex*n1*n2*n3) 1153 1154 !Local variables------------------------------- 1155 !scalars 1156 integer,parameter :: mshift=401 1157 integer :: i1,i2,i3,iatom,ifft,ishift,ishift1,ishift2,ishift3,istr 1158 integer :: ixp,jj,ka,kb,mrange,mu,nu 1159 real(dp) :: aa,bb,cc,dd,delta,delta2div6,deltam1,diff,difmag 1160 real(dp) :: difmag2,difmag2_fact,difmag2_part,func,phase,phi,phr,prod 1161 real(dp) :: range,range2,rangem1,rdiff1,rdiff2,rdiff3,term 1162 real(dp) :: yy 1163 character(len=500) :: message 1164 !arrays 1165 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 1166 integer :: igrid(3),irange(3),ngfft(3) 1167 integer,allocatable :: ii(:,:) 1168 real(dp) :: drmetds(3,3),lencp(3),rmet(3,3),scale(3),tau(3) 1169 real(dp),allocatable :: rrdiff(:,:) 1170 1171 ! ************************************************************************* 1172 1173 ! if( ipert<1 .or. ipert> natom+7) then 1174 ! write(message,'(a,i0,a,a,a,i0,a)')& 1175 !& ' The argument ipert must be between 1 and natom+7=',natom+7,',',ch10,& 1176 !& ' while it is ipert=',ipert,'.' 1177 ! ABI_BUG(message) 1178 ! end if 1179 1180 if( (ipert==natom+3 .or. ipert==natom+4) .and. cplex/=1) then 1181 write(message,'(3a,i4,a)')& 1182 & 'The argument cplex must be 1 for strain perturbationh',ch10,& 1183 & 'while it is cplex=',cplex,'.' 1184 ABI_BUG(message) 1185 end if 1186 1187 !Zero out array 1188 xccc3d1(:)=0.0_dp 1189 1190 !For a non-linear XC core correction, the perturbation must be phonon-type or strain type 1191 if(ipert<=natom .or. ipert==natom+3 .or. ipert==natom+4) then 1192 1193 if( idir<1 .or. idir> 3) then 1194 write(message,'(a,a,a,i4,a)')& 1195 & 'The argument idir must be between 1 and 3,',ch10,& 1196 & 'while it is idir=',idir,'.' 1197 ABI_BUG(message) 1198 end if 1199 1200 ! Compute lengths of cross products for pairs of primitive 1201 ! translation vectors (used in setting index search range below) 1202 lencp(1)=cross_mk(rprimd(1,2),rprimd(2,2),rprimd(3,2),& 1203 & rprimd(1,3),rprimd(2,3),rprimd(3,3)) 1204 lencp(2)=cross_mk(rprimd(1,3),rprimd(2,3),rprimd(3,3),& 1205 & rprimd(1,1),rprimd(2,1),rprimd(3,1)) 1206 lencp(3)=cross_mk(rprimd(1,1),rprimd(2,1),rprimd(3,1),& 1207 & rprimd(1,2),rprimd(2,2),rprimd(3,2)) 1208 1209 ! Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3 1210 ! (recall ucvol=R1.(R2xR3)) 1211 scale(:)=ucvol/lencp(:) 1212 1213 ! Compute metric tensor in real space rmet 1214 do nu=1,3 1215 rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+& 1216 & rprimd(3,:)*rprimd(3,nu) 1217 end do 1218 1219 ! Section to be executed only for strain perturbation 1220 ! Compute derivative of metric tensor wrt strain component istr 1221 if(ipert==natom+3 .or. ipert==natom+4) then 1222 istr=idir + 3*(ipert-natom-3) 1223 1224 ka=idx(2*istr-1);kb=idx(2*istr) 1225 do jj = 1,3 1226 drmetds(:,jj)=(rprimd(ka,:)*rprimd(kb,jj)+rprimd(kb,:)*rprimd(ka,jj)) 1227 end do 1228 ! For historical reasons: 1229 drmetds(:,:)=0.5_dp*drmetds(:,:) 1230 1231 ! end of strain perturbation section 1232 end if 1233 1234 ngfft(1)=n1 1235 ngfft(2)=n2 1236 ngfft(3)=n3 1237 1238 delta=1.0_dp/(n1xccc-1) 1239 deltam1=n1xccc-1 1240 delta2div6=delta**2/6.0_dp 1241 1242 ! Loop over atoms in unit cell 1243 ! Note that we cycle immediately for all except the displaced atom 1244 ! for such a perturbation. The loop is executed over all the 1245 ! atoms for a strain peturbation. 1246 do iatom=1,natom 1247 if(ipert<=natom .and. iatom/=ipert) cycle 1248 ! Set search range (density cuts off perfectly beyond range) 1249 ! Cycle if no range. 1250 range=0.0_dp 1251 range=xcccrc(typat(iatom)) 1252 if(range<1.d-16) cycle 1253 1254 range2=range**2 1255 rangem1=1.0_dp/range 1256 1257 ! compute mrange for ii(:,3), rrdiff(:,3), inserted by MM (2005/12/06) 1258 ! Consider each component in turn : compute range 1259 do mu=1,3 1260 1261 ! Convert reduced coord of given atom to [0,1) 1262 tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp) 1263 1264 ! Use tau to find nearest grid point along R(mu) 1265 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 1266 igrid(mu)=nint(tau(mu)*dble(ngfft(mu))) 1267 1268 ! Use range to compute an index range along R(mu) 1269 ! (add 1 to make sure it covers full range) 1270 irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu))) 1271 1272 end do 1273 1274 ! Allocate arrays that depends on the range 1275 mrange=maxval(irange(1:3)) 1276 ABI_MALLOC(ii,(2*mrange+1,3)) 1277 ABI_MALLOC(rrdiff,(2*mrange+1,3)) 1278 1279 ! Consider each component in turn 1280 do mu=1,3 1281 1282 ! temporarily suppressed by MM (2005/12/02) 1283 ! Convert reduced coord of given atom to [0,1) 1284 ! tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp) 1285 1286 ! Use tau to find nearest grid point along R(mu) 1287 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 1288 ! igrid(mu)=nint(tau(mu)*dble(ngfft(mu))) 1289 1290 ! Use range to compute an index range along R(mu) 1291 ! (add 1 to make sure it covers full range) 1292 ! irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu))) 1293 1294 ! Check that the largest range is smallest than the maximum 1295 ! allowed one 1296 ! if(2*irange(mu)+1 > mshift)then 1297 ! write(message, '(a,a,a,a,i6,a)' ) ch10,& 1298 ! & ' dfpt_mkcore : BUG -',ch10,& 1299 ! & ' The range around atom',iatom,' is too large.' 1300 ! ABI_BUG(message) 1301 ! end if 1302 1303 ! Set up a counter that explore the relevant range 1304 ! of points around the atom 1305 ishift=0 1306 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 1307 ishift=ishift+1 1308 ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu)) 1309 rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu) 1310 end do 1311 1312 ! End loop on mu 1313 end do 1314 1315 ! Conduct triple loop over restricted range of grid points for iatom 1316 1317 do ishift3=1,1+2*irange(3) 1318 ! map back to [1,ngfft(3)] for usual fortran index in unit cell 1319 i3=ii(ishift3,3) 1320 ! find vector from atom location to grid point (reduced) 1321 rdiff3=rrdiff(ishift3,3) 1322 1323 do ishift2=1,1+2*irange(2) 1324 i2=ii(ishift2,2) 1325 rdiff2=rrdiff(ishift2,2) 1326 ! Prepare the computation of difmag2 1327 difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2& 1328 & +2.0_dp*rmet(3,2)*rdiff3*rdiff2 1329 difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2) 1330 1331 do ishift1=1,1+2*irange(1) 1332 rdiff1=rrdiff(ishift1,1) 1333 1334 ! Compute (rgrid-tau-Rprim)**2 1335 difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1) 1336 1337 ! Only accept contribution inside defined range 1338 if (difmag2<range2) then 1339 1340 ! Prepare computation of core charge function and derivative, 1341 ! using splines 1342 difmag=sqrt(difmag2) 1343 if (difmag>=1.0d-10) then 1344 i1=ii(ishift1,1) 1345 yy=difmag*rangem1 1346 1347 ! Compute index of yy over 1 to n1xccc scale 1348 jj=1+int(yy*(n1xccc-1)) 1349 diff=yy-(jj-1)*delta 1350 1351 ! Will evaluate spline fit (p. 86 Numerical Recipes, Press et al; 1352 ! NOTE error in book for sign of "aa" term in derivative; 1353 ! also see splfit routine). 1354 bb = diff*deltam1 1355 aa = 1.0_dp-bb 1356 cc = aa*(aa**2-1.0_dp)*delta2div6 1357 dd = bb*(bb**2-1.0_dp)*delta2div6 1358 1359 ! Evaluate spline fit of 1st der of core charge density 1360 ! from xccc1d(:,2,:) and (:,4,:) 1361 func=aa*xccc1d(jj,2,typat(iatom))+bb*xccc1d(jj+1,2,typat(iatom)) +& 1362 & cc*xccc1d(jj,4,typat(iatom))+dd* xccc1d(jj+1,4,typat(iatom)) 1363 1364 if(ipert<=natom) then 1365 phase=2*pi*(qphon(1)*(rdiff1+xred(1,iatom)) & 1366 & +qphon(2)*(rdiff2+xred(2,iatom)) & 1367 & +qphon(3)*(rdiff3+xred(3,iatom))) 1368 prod=rmet(idir,1)*rdiff1+rmet(idir,2)*rdiff2+rmet(idir,3)*rdiff3 1369 1370 term=-func*rangem1/difmag*prod 1371 ifft=i1+n1*(i2-1+n2*(i3-1)) 1372 phr=cos(phase) 1373 if(cplex==1)then 1374 xccc3d1(ifft)=xccc3d1(ifft)+term*phr 1375 else 1376 phi=sin(phase) 1377 xccc3d1(2*ifft-1)=xccc3d1(2*ifft-1)+term*phr 1378 xccc3d1(2*ifft )=xccc3d1(2*ifft )-term*phi 1379 end if 1380 else 1381 prod=& 1382 & (rdiff1*(drmetds(1,1)*rdiff1+drmetds(1,2)*rdiff2+drmetds(1,3)*rdiff3)& 1383 & +rdiff2*(drmetds(2,1)*rdiff1+drmetds(2,2)*rdiff2+drmetds(2,3)*rdiff3)& 1384 & +rdiff3*(drmetds(3,1)*rdiff1+drmetds(3,2)*rdiff2+drmetds(3,3)*rdiff3)) 1385 term=prod*func*rangem1/difmag 1386 1387 ifft=i1+n1*(i2-1+n2*(i3-1)) 1388 xccc3d1(ifft)=xccc3d1(ifft)+term 1389 1390 end if 1391 1392 ! End of the condition for the distance not to vanish 1393 end if 1394 1395 ! End of condition to be inside the range 1396 end if 1397 1398 ! End loop on ishift1 1399 end do 1400 1401 ! End loop on ishift2 1402 end do 1403 1404 ! End loop on ishift3 1405 end do 1406 1407 ABI_FREE(ii) 1408 ABI_FREE(rrdiff) 1409 ! End loop on atoms 1410 end do 1411 1412 ! End of the condition ipert corresponds to a phonon type perturbation 1413 ! or strain type perturbation 1414 end if 1415 1416 contains 1417 1418 function cross_mk(xx,yy,zz,aa,bb,cc) 1419 1420 real(dp) :: cross_mk 1421 real(dp),intent(in) :: xx,yy,zz,aa,bb,cc 1422 cross_mk=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2) 1423 end function cross_mk 1424 1425 end subroutine dfpt_mkcore
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
SOURCE
1082 subroutine indpos_mkcore_alt(periodic,ii,nn,jj,inside) 1083 ! Find the grid index of a given position in the cell according to the BC 1084 ! Determine also whether the index is inside or outside the box for free BC 1085 integer, intent(in) :: ii,nn 1086 integer, intent(out) :: jj 1087 logical, intent(in) :: periodic 1088 logical, intent(out) :: inside 1089 ! ************************************************************************* 1090 if (periodic) then 1091 inside=.true. ; jj=modulo(ii-1,nn)+1 1092 else 1093 jj=ii ; inside=(ii>=1.and.ii<=nn) 1094 end if 1095 end subroutine indpos_mkcore_alt 1096 1097 end subroutine mkcore_alt
ABINIT/m_mkcore [ Modules ]
NAME
m_mkcore
FUNCTION
Routines related to non-linear core correction.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_mkcore 23 24 use defs_basis 25 use m_abicore 26 use m_xmpi 27 use m_errors 28 use m_linalg_interfaces 29 30 use defs_abitypes, only : mpi_type 31 use m_geometry, only : strconv 32 use m_time, only : timab 33 use m_mpinfo, only : ptabs_fourdp 34 use m_sort, only : sort_dp 35 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 36 use m_pawtab, only : pawtab_type 37 use m_paw_numeric, only : paw_splint 38 39 implicit none 40 41 private
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)
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
SOURCE
100 subroutine mkcore(corstr,dyfrx2,grxc,mpi_enreg,natom,nfft,nspden,ntypat,n1,n1xccc,& 101 & n2,n3,option,rprimd,typat,ucvol,vxc,xcccrc,xccc1d,xccc3d,xred) 102 103 !Arguments ------------------------------------ 104 !scalars 105 integer,intent(in) :: n1,n1xccc,n2,n3,natom,nfft,nspden,ntypat,option 106 real(dp),intent(in) :: ucvol 107 type(mpi_type),intent(in) :: mpi_enreg 108 !arrays 109 integer,intent(in) :: typat(natom) 110 real(dp),intent(in) :: rprimd(3,3),vxc(nfft,nspden),xccc1d(n1xccc,6,ntypat) 111 real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom) 112 real(dp),intent(inout) :: xccc3d(nfft) 113 real(dp),intent(out) :: corstr(6),dyfrx2(3,3,natom) 114 real(dp),intent(inout) :: grxc(3,natom) 115 116 !Local variables------------------------------- 117 !scalars 118 integer :: i1,i2,i3,iatom,ier,ifft,ishift,ishift1,ishift2 119 integer :: ishift3,itypat,ixp,jj,me_fft,mrange,mu,nfftot,nu 120 real(dp) :: dd,delta,delta2div6,deltam1,diff,difmag,difmag2 121 real(dp) :: difmag2_fact,difmag2_part,fact,func,grxc1,grxc2,grxc3,range,range2 122 real(dp) :: rangem1,rdiff1,rdiff2,rdiff3,strdia,t1,t2,t3,term,term1,term2 123 character(len=500) :: message 124 !arrays 125 integer :: igrid(3),irange(3),ngfft(3) 126 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 127 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 128 integer,allocatable :: ii(:,:) 129 real(dp) :: yy,aa,bb,cc 130 real(dp) :: corfra(3,3),lencp(3),rmet(3,3),scale(3),tau(3),tsec(2),tt(3) 131 real(dp),allocatable :: rrdiff(:,:),work(:,:,:) 132 133 !************************************************************************ 134 135 call timab(12,1,tsec) 136 137 !Make sure option is acceptable 138 if (option<0.or.option>4) then 139 write(message, '(a,i12,a,a,a)' )& 140 'option=',option,' is not allowed.',ch10,& 141 'Must be 1, 2, 3 or 4.' 142 ABI_BUG(message) 143 end if 144 145 !Zero out only the appropriate array according to option: 146 !others are dummies with no storage 147 148 if (option==1) then 149 ! Zero out array to permit accumulation over atom types below: 150 xccc3d(:)=zero 151 else if (option==2) then 152 ! Zero out gradient of Exc array 153 grxc(:,:)=zero 154 else if (option==3) then 155 ! Zero out locally defined stress array 156 corfra(:,:)=zero 157 strdia=zero 158 else if (option==4) then 159 ! Zero out fr-wf part of the dynamical matrix 160 dyfrx2(:,:,:)=zero 161 else 162 ABI_BUG(" Can't be here! (bad option)") 163 end if 164 165 !Compute lengths of cross products for pairs of primitive 166 !translation vectors (used in setting index search range below) 167 lencp(1)=cross_mkcore(rprimd(1,2),rprimd(2,2),rprimd(3,2),& 168 & rprimd(1,3),rprimd(2,3),rprimd(3,3)) 169 lencp(2)=cross_mkcore(rprimd(1,3),rprimd(2,3),rprimd(3,3),& 170 & rprimd(1,1),rprimd(2,1),rprimd(3,1)) 171 lencp(3)=cross_mkcore(rprimd(1,1),rprimd(2,1),rprimd(3,1),& 172 & rprimd(1,2),rprimd(2,2),rprimd(3,2)) 173 174 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3 175 !(recall ucvol=R1.(R2xR3)) 176 scale(:)=ucvol/lencp(:) 177 178 !Compute metric tensor in real space rmet 179 do nu=1,3 180 rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu) 181 end do 182 183 ngfft(1)=n1 184 ngfft(2)=n2 185 ngfft(3)=n3 186 nfftot=n1*n2*n3 187 me_fft = mpi_enreg%me_fft 188 189 ! Get the distrib associated with this fft_grid 190 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 191 192 delta=one/(n1xccc-1) 193 deltam1=n1xccc-1 194 delta2div6=delta**2/6.0d0 195 196 if (option>=2) then 197 ABI_MALLOC(work,(n1,n2,n3)) 198 ! For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down)) 199 ! For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22}) 200 if (nspden>=2) then 201 ifft=1 202 do i3=1,n3 203 if(me_fft==fftn3_distrib(i3)) then 204 do i2=1,n2 205 do i1=1,n1 206 work(i1,i2,i3)=half*(vxc(ifft,1)+vxc(ifft,2)) 207 ifft=ifft+1 208 end do 209 end do 210 end if 211 end do 212 else 213 ifft=1 214 do i3=1,n3 215 if(me_fft==fftn3_distrib(i3)) then 216 do i2=1,n2 217 do i1=1,n1 218 work(i1,i2,i3)=vxc(ifft,1) 219 ifft=ifft+1 220 end do 221 end do 222 end if 223 end do 224 ! call DCOPY(nfft,vxc,1,work,1) 225 end if 226 end if 227 228 !Loop over atoms in unit cell 229 do iatom=1,natom 230 231 if(option==2)then 232 grxc1=zero 233 grxc2=zero 234 grxc3=zero 235 end if 236 237 ! Set search range (density cuts off perfectly beyond range) 238 itypat=typat(iatom) 239 range=xcccrc(itypat) 240 241 ! Skip loop if this atom has no core charge 242 if (abs(range)<1.d-16) cycle 243 244 range2=range**2 245 rangem1=one/range 246 247 ! Consider each component in turn : compute range 248 do mu=1,3 249 250 ! Convert reduced coord of given atom to [0,1) 251 tau(mu)=mod(xred(mu,iatom)+one-aint(xred(mu,iatom)),one) 252 253 ! Use tau to find nearest grid point along R(mu) 254 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 255 igrid(mu)=nint(tau(mu)*dble(ngfft(mu))) 256 257 ! Use range to compute an index range along R(mu) 258 ! (add 1 to make sure it covers full range) 259 irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu))) 260 261 end do 262 263 ! Allocate arrays that depends on the range 264 mrange=maxval(irange(1:3)) 265 ABI_MALLOC(ii,(2*mrange+1,3)) 266 ABI_MALLOC(rrdiff,(2*mrange+1,3)) 267 268 ! Set up counters that explore the relevant range 269 ! of points around the atom 270 do mu=1,3 271 ishift=0 272 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 273 ishift=ishift+1 274 ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu)) 275 rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu) 276 end do 277 end do 278 279 ! Conduct triple loop over restricted range of grid points for iatom 280 do ishift3=1,1+2*irange(3) 281 ! map back to [1,ngfft(3)] for usual fortran index in unit cell 282 i3=ii(ishift3,3) 283 if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle 284 ! find vector from atom location to grid point (reduced) 285 rdiff3=rrdiff(ishift3,3) 286 287 do ishift2=1,1+2*irange(2) 288 i2=ii(ishift2,2) 289 rdiff2=rrdiff(ishift2,2) 290 ! Prepare the computation of difmag2 291 difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2& 292 & +2.0d0*rmet(3,2)*rdiff3*rdiff2 293 difmag2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2) 294 295 do ishift1=1,1+2*irange(1) 296 rdiff1=rrdiff(ishift1,1) 297 298 ! Compute (rgrid-tau-Rprim)**2 299 difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1) 300 301 ! Only accept contribution inside defined range 302 if (difmag2<range2-tol12) then 303 304 ! Prepare computation of core charge function and derivative, 305 ! using splines 306 i1=ii(ishift1,1) 307 difmag=sqrt(difmag2) 308 yy=difmag*rangem1 309 310 ! Compute index of yy over 1 to n1xccc scale 311 jj=1+int(yy*(n1xccc-1)) 312 diff=yy-(jj-1)*delta 313 314 ! Will evaluate spline fit (p. 86 Numerical Recipes, Press et al; 315 ! NOTE error in book for sign of "aa" term in derivative; 316 ! also see splfit routine). 317 bb = diff*deltam1 318 aa = one-bb 319 cc = aa*(aa**2-one)*delta2div6 320 dd = bb*(bb**2-one)*delta2div6 321 322 323 ! Test first for option 2, the most frequently used 324 if (option==2) then 325 326 ! Accumulate contributions to Exc gradients 327 328 if (difmag>1.0d-10) then 329 330 ! Evaluate spline fit of 1st der of core charge density 331 ! from xccc1d(:,2,:) and (:,4,:) 332 func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 333 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 334 term=work(i1,i2,i3)*func/difmag 335 grxc1=grxc1+rdiff1*term 336 grxc2=grxc2+rdiff2*term 337 grxc3=grxc3+rdiff3*term 338 end if 339 340 else if (option==1) then 341 342 ! Evaluate spline fit of core charge density 343 ! from xccc1d(:,1,:) and (:,3,:) 344 func=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +& 345 & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat) 346 347 ! Accumulate contributions to core electron density 348 ! throughout unit cell 349 ifft=i1+n1*(i2-1+n2*(ffti3_local(i3)-1)) 350 xccc3d(ifft)=xccc3d(ifft)+func 351 352 else if (option==3) then 353 354 ! Accumulate contributions to stress tensor 355 ! in reduced coordinates 356 357 if (difmag>1.0d-10) then 358 359 ! Evaluate spline fit of 1st der of core charge density 360 ! from xccc1d(:,2,:) and (:,4,:) 361 func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 362 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 363 term=work(i1,i2,i3)*func*rangem1/difmag/dble(n1*n2*n3) 364 ! Write out the 6 symmetric components 365 corfra(1,1)=corfra(1,1)+term*rdiff1**2 366 corfra(2,2)=corfra(2,2)+term*rdiff2**2 367 corfra(3,3)=corfra(3,3)+term*rdiff3**2 368 corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2 369 corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1 370 corfra(2,1)=corfra(2,1)+term*rdiff2*rdiff1 371 ! (the above still needs to be transformed to cartesian coords) 372 373 end if 374 375 ! Also compute a diagonal term 376 ! Evaluate spline fit of core charge density 377 ! from xccc1d(:,1,:) and (:,3,:) 378 func=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +& 379 & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat) 380 strdia=strdia+work(i1,i2,i3)*func 381 382 else if (option==4) then 383 384 ! Compute frozen-wf contribution to Dynamical matrix 385 386 tt(1)=rmet(1,1)*rdiff1+rmet(1,2)*rdiff2+rmet(1,3)*rdiff3 387 tt(2)=rmet(2,1)*rdiff1+rmet(2,2)*rdiff2+rmet(2,3)*rdiff3 388 tt(3)=rmet(3,1)*rdiff1+rmet(3,2)*rdiff2+rmet(3,3)*rdiff3 389 390 if (difmag>1.d-10) then 391 392 ! Accumulate contributions to dynamical matrix 393 term=(ucvol/dble(nfftot))*work(i1,i2,i3)*rangem1/difmag 394 ! Evaluate spline fit of 1st der of core charge density 395 ! from xccc1d(:,2,:) and (:,4,:) 396 func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 397 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 398 term1=term*func 399 ! Evaluate spline fit of 2nd der of core charge density 400 ! from xccc1d(:,3,:) and (:,5,:) 401 func=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +& 402 & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat) 403 term2=term*func*rangem1/difmag 404 do mu=1,3 405 do nu=1,3 406 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)& 407 & +(term2-term1/difmag**2)*tt(mu)*tt(nu)& 408 & +term1*rmet(mu,nu) 409 end do 410 end do 411 412 else 413 414 ! There is a contribution from difmag=zero ! 415 ! Evaluate spline fit of 2nd der of core charge density 416 ! from xccc1d(:,3,:) and (:,5,:) 417 func=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +& 418 & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat) 419 term=(ucvol/dble(nfftot))*work(i1,i2,i3)*func*rangem1**2 420 do mu=1,3 421 do nu=1,3 422 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu) 423 end do 424 end do 425 426 ! End of condition not to be precisely on the point (difmag=zero) 427 end if 428 429 ! If option is not 1, 2, 3, or 4. 430 else 431 ABI_BUG("Can't be here in mkcore") 432 ! End of choice of option 433 end if 434 435 end if ! End of condition on the range 436 end do ! End loop on ishift1 437 end do ! End loop on ishift2 438 end do ! End loop on ishift3 439 440 ABI_FREE(ii) 441 ABI_FREE(rrdiff) 442 443 if(option==2)then 444 fact=-(ucvol/dble(nfftot))/range 445 grxc(1,iatom)=grxc1*fact 446 grxc(2,iatom)=grxc2*fact 447 grxc(3,iatom)=grxc3*fact 448 end if 449 450 end do ! End big loop on atoms 451 452 if (option==2) then 453 454 ! Apply rmet as needed to get reduced coordinate gradients 455 do iatom=1,natom 456 t1=grxc(1,iatom) 457 t2=grxc(2,iatom) 458 t3=grxc(3,iatom) 459 grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3 460 461 end do 462 end if 463 464 if (option==3) then 465 466 ! Transform stress tensor from full storage mode to symmetric storage mode 467 corstr(1)=corfra(1,1) 468 corstr(2)=corfra(2,2) 469 corstr(3)=corfra(3,3) 470 corstr(4)=corfra(3,2) 471 corstr(5)=corfra(3,1) 472 corstr(6)=corfra(2,1) 473 474 ! Transform stress tensor from reduced coordinates to cartesian coordinates 475 call strconv(corstr,rprimd,corstr) 476 477 ! Compute diagonal contribution to stress tensor (need input xccc3d) 478 ! strdia = (1/N) Sum(r) [mu_xc_avg(r) * rho_core(r)] 479 ifft=0 ; strdia=zero 480 do i3=1,n3 481 if(me_fft==fftn3_distrib(i3)) then 482 do i2=1,n2 483 do i1=1,n1 484 ifft=ifft+1 485 strdia=strdia+work(i1,i2,i3)*xccc3d(ifft) 486 end do 487 end do 488 end if 489 end do 490 strdia=strdia/dble(nfftot) 491 ! strdia=DDOT(nfft,work,1,xccc3d,1)/dble(nfftot) 492 493 ! Add diagonal term to stress tensor 494 corstr(1)=corstr(1)+strdia 495 corstr(2)=corstr(2)+strdia 496 corstr(3)=corstr(3)+strdia 497 end if 498 499 if(option>=2) then 500 ABI_FREE(work) 501 end if 502 503 if(mpi_enreg%nproc_fft > 1) then 504 call timab(539,1,tsec) 505 if(option==2) then 506 call xmpi_sum(grxc,mpi_enreg%comm_fft,ier) 507 end if 508 if(option==3) then 509 call xmpi_sum(corstr,mpi_enreg%comm_fft,ier) 510 end if 511 if(option==4) then 512 call xmpi_sum(dyfrx2,mpi_enreg%comm_fft,ier) 513 end if 514 call timab(539,2,tsec) 515 end if 516 517 call timab(12,2,tsec) 518 519 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.
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=information 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
[usekden]= --optional-- if TRUE, output the kinetic enrgy density instead of the density
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
SOURCE
601 subroutine mkcore_alt(atindx1,corstr,dyfrx2,grxc,icoulomb,mpi_enreg,natom,nfft,nspden,& 602 & nattyp,ntypat,n1,n1xccc,n2,n3,option,rprimd,ucvol,vxc,xcccrc,xccc1d,& 603 & xccc3d,xred,pawrad,pawtab,usepaw,& 604 & usekden) ! optional argument 605 606 !Arguments ------------------------------------ 607 !scalars 608 integer,intent(in) :: icoulomb,n1,n1xccc,n2,n3,natom,nfft,nspden,ntypat,option,usepaw 609 logical,intent(in),optional :: usekden 610 real(dp),intent(in) :: ucvol 611 type(mpi_type),intent(in) :: mpi_enreg 612 type(pawrad_type),intent(in) :: pawrad(:) 613 type(pawtab_type),target,intent(in) :: pawtab(:) 614 !arrays 615 integer,intent(in) :: atindx1(natom),nattyp(ntypat) 616 real(dp),intent(in) :: rprimd(3,3),xccc1d(n1xccc,6,ntypat) 617 real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom) 618 real(dp),intent(in),target :: vxc(nfft,nspden) 619 real(dp),intent(out) :: corstr(6),grxc(3,natom),dyfrx2(3,3,natom) 620 real(dp),intent(inout) :: xccc3d(nfft) 621 622 !Local variables------------------------------- 623 !scalars 624 integer :: i1,i2,i3,iat,iatm,iatom,ier,ipts 625 integer :: ishift,ishift1,ishift2,ishift3 626 integer :: itypat,ixp,jj,jpts,me_fft,mrange,msz,mu 627 integer :: nfftot,npts,npts12,nu 628 logical :: letsgo,usekden_ 629 real(dp) :: aa,bb,cc,dd,delta,delta2div6,deltam1 630 real(dp) :: diff,difmag,fact,range,range2 631 real(dp) :: rangem1,rdiff1,rdiff2,rdiff3 632 real(dp) :: rnorm2,rnorm2_fact,rnorm2_part 633 real(dp) :: strdia,t1,t2,t3,term,term1,term2,yy 634 character(len=1) :: geocode 635 character(len=500) :: message 636 type(pawrad_type) :: core_mesh 637 !arrays 638 integer :: igrid(3),irange(3),ishiftmax(3),ngfft(3) 639 integer,allocatable :: ii(:,:),iindex(:),indx1(:),indx2(:) 640 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 641 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 642 logical :: per(3) 643 real(dp) :: corfra(3,3),corgr(3),lencp(3),rmet(3,3) 644 real(dp) :: scale(3),tau(3),tsec(2),tt(3) 645 real(dp),allocatable :: dtcore(:),d2tcore(:),rnorm(:) 646 real(dp),allocatable :: rrdiff(:,:),tcore(:) 647 real(dp),allocatable,target :: tcoretau(:,:) 648 real(dp), ABI_CONTIGUOUS pointer :: corespl(:,:),vxc_eff(:) 649 650 !************************************************************************ 651 652 call timab(12,1,tsec) 653 654 !Make sure options are acceptable 655 usekden_=.false.;if (present(usekden)) usekden_=usekden 656 if (option<0.or.option>4) then 657 write(message, '(a,i12,a,a,a)' )& 658 'option=',option,' is not allowed.',ch10,& 659 'Must be 1, 2, 3 or 4.' 660 ABI_BUG(message) 661 end if 662 if (usekden_) then 663 message='usekden=1 not yet allowed!' 664 ABI_BUG(message) 665 end if 666 667 668 !Zero out only the appropriate array according to option: 669 if (option==1) then 670 xccc3d(:)=zero 671 else if (option==2) then 672 grxc(:,:)=zero 673 else if (option==3) then 674 corfra(:,:)=zero 675 strdia=zero 676 else if (option==4) then 677 dyfrx2(:,:,:)=zero 678 end if 679 680 !Conditions for periodicity in the three directions 681 geocode='P' 682 if (icoulomb==1) geocode='F' 683 if (icoulomb==2) geocode='S' 684 per(1)=(geocode /= 'F') 685 per(2)=(geocode == 'P') 686 per(3)=(geocode /= 'F') 687 688 !Compute lengths of cross products for pairs of primitive 689 !translation vectors (used in setting index search range below) 690 lencp(1)=cross_mkcore_alt(rprimd(1,2),rprimd(2,2),rprimd(3,2),& 691 & rprimd(1,3),rprimd(2,3),rprimd(3,3)) 692 lencp(2)=cross_mkcore_alt(rprimd(1,3),rprimd(2,3),rprimd(3,3),& 693 & rprimd(1,1),rprimd(2,1),rprimd(3,1)) 694 lencp(3)=cross_mkcore_alt(rprimd(1,1),rprimd(2,1),rprimd(3,1),& 695 & rprimd(1,2),rprimd(2,2),rprimd(3,2)) 696 697 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3 698 !(recall ucvol=R1.(R2xR3)) 699 scale(:)=ucvol/lencp(:) 700 701 !Compute metric tensor in real space rmet 702 do nu=1,3 703 rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu) 704 end do 705 706 !Get the distrib associated with this fft_grid 707 ngfft(1)=n1;ngfft(2)=n2;ngfft(3)=n3 708 nfftot=n1*n2*n3 ; me_fft=mpi_enreg%me_fft 709 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 710 711 delta=one/(n1xccc-1) 712 deltam1=n1xccc-1 713 delta2div6=delta**2/6.0_dp 714 715 if (option>=2) then 716 ! For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down)) 717 ! For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22}) 718 if (nspden>=2) then 719 ABI_MALLOC(vxc_eff,(nfft)) 720 do jj=1,nfft 721 vxc_eff(jj)=half*(vxc(jj,1)+vxc(jj,2)) 722 end do 723 else 724 vxc_eff => vxc(1:nfft,1) 725 end if 726 end if 727 728 !Loop over atom types 729 iatm=0 730 do itypat=1,ntypat 731 732 ! Set search range (density cuts off perfectly beyond range) 733 range=xcccrc(itypat);if (usepaw==1) range=pawtab(itypat)%rcore 734 range2=range**2 ; rangem1=one/range 735 736 ! Skip loop if this type has no core charge 737 if (abs(range)<1.d-16) cycle 738 739 ! PAW: select core density type and create mesh 740 if (usepaw==1) then 741 if (usekden_) then 742 msz=pawtab(itypat)%coretau_mesh_size 743 ABI_MALLOC(tcoretau,(msz,1)) 744 tcoretau(:,1)=pawtab(itypat)%coretau(:) 745 corespl => tcoretau 746 else 747 msz=pawtab(itypat)%core_mesh_size 748 corespl => pawtab(itypat)%tcoredens 749 end if 750 call pawrad_init(core_mesh,mesh_size=msz,& 751 & mesh_type=pawrad(itypat)%mesh_type,& 752 & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep) 753 end if 754 755 ! Loop over atoms of the type 756 do iat=1,nattyp(itypat) 757 iatm=iatm+1;iatom=atindx1(iatm) 758 759 if(option==2) corgr(:)=zero 760 761 ! Consider each component in turn : compute range 762 do mu=1,3 763 ! Convert reduced coord of given atom to [0,1) 764 tau(mu)=mod(xred(mu,iatom)+one-aint(xred(mu,iatom)),one) 765 ! Use tau to find nearest grid point along R(mu) 766 ! (igrid=0 is the origin; shift by 1 to agree with usual index) 767 igrid(mu)=nint(tau(mu)*real(ngfft(mu),dp)) 768 ! Use range to compute an index range along R(mu) 769 ! (add 1 to make sure it covers full range) 770 irange(mu)=1+nint((range/scale(mu))*real(ngfft(mu),dp)) 771 end do 772 773 ! Allocate arrays that depends on the range 774 mrange=maxval(irange(1:3)) 775 ABI_MALLOC(ii,(2*mrange+1,3)) 776 ABI_MALLOC(rrdiff,(2*mrange+1,3)) 777 778 ! Set up counters that explore the relevant range of points around the atom 779 if (geocode=='P') then 780 ! Fully periodic version 781 do mu=1,3 782 ishift=0 783 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 784 ishift=ishift+1 785 ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu)) 786 rrdiff(ishift,mu)=real(ixp,dp)/real(ngfft(mu),dp)-tau(mu) 787 end do 788 ishiftmax(mu)=ishift 789 end do 790 else 791 ! Free or surface conditions 792 do mu=1,3 793 ishift=0 794 do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu) 795 call indpos_mkcore_alt(per(mu),ixp,ngfft(mu),jj,letsgo) 796 if (letsgo) then 797 ishift=ishift+1;ii(ishift,mu)=1+jj 798 rrdiff(ishift,mu)=real(ixp,dp)/real(ngfft(mu),dp)-tau(mu) 799 end if 800 end do 801 ishiftmax(mu)=ishift 802 end do 803 end if 804 npts12=ishiftmax(1)*ishiftmax(2) 805 ABI_MALLOC(indx1,(npts12)) 806 ABI_MALLOC(indx2,(npts12)) 807 ABI_MALLOC(iindex,(npts12)) 808 ABI_MALLOC(rnorm,(npts12)) 809 if (option==1.or.option==3) then 810 ABI_MALLOC(tcore,(npts12)) 811 end if 812 if (option>=2) then 813 ABI_MALLOC(dtcore,(npts12)) 814 end if 815 if (option==4) then 816 ABI_MALLOC(d2tcore,(npts12)) 817 end if 818 819 ! Conduct loop over restricted range of grid points for iatom 820 do ishift3=1,ishiftmax(3) 821 i3=ii(ishift3,3) 822 rdiff3=rrdiff(ishift3,3) 823 824 if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle 825 826 ! Select the vectors located around the current atom 827 ! TR: all of the following could be done inside or 828 ! outside the loops (i2,i1,i3). 829 ! Outside: the memory consumption increases. 830 ! Inside: the time of calculation increases. 831 ! Here, I choose to do it here, somewhere in the middle. 832 npts=0 833 do ishift2=1,ishiftmax(2) 834 i2=ii(ishift2,2) ; rdiff2=rrdiff(ishift2,2) 835 rnorm2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2 & 836 & +2.0d0*rmet(3,2)*rdiff3*rdiff2 837 rnorm2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2) 838 do ishift1=1,ishiftmax(1) 839 i1=ii(ishift1,1) ; rdiff1=rrdiff(ishift1,1) 840 rnorm2=rnorm2_part+rdiff1*(rnorm2_fact+rmet(1,1)*rdiff1) 841 ! Only accept contributions inside defined range 842 if (rnorm2<range2-tol12) then 843 npts=npts+1 ; iindex(npts)=npts 844 indx1(npts)=ishift1;indx2(npts)=ishift2 845 rnorm(npts)=sqrt(rnorm2) 846 end if 847 end do 848 end do 849 if (npts==0) cycle 850 if (npts>npts12) then 851 message='npts>npts12!' 852 ABI_BUG(message) 853 end if 854 855 ! Evaluate core density (and derivatives) on the set of selected points 856 if (usepaw==1) then 857 ! PAW: use splint routine 858 call sort_dp(npts,rnorm(1:npts),iindex(1:npts),tol16) 859 if (option==1.or.option==3) then 860 ! Evaluate fit of core density 861 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 862 & corespl(:,1),corespl(:,3),& 863 & npts,rnorm(1:npts),tcore(1:npts)) 864 end if 865 if (option>=2) then 866 ! Evaluate fit of 1-der of core density 867 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 868 & corespl(:,2),corespl(:,4),& 869 & npts,rnorm(1:npts),dtcore(1:npts)) 870 end if 871 if (option==4) then 872 ! Evaluate fit of 2nd-der of core density 873 call paw_splint(core_mesh%mesh_size,core_mesh%rad, & 874 & corespl(:,3),corespl(:,5),& 875 & npts,rnorm(1:npts),d2tcore(1:npts)) 876 end if 877 else 878 ! Norm-conserving PP: 879 ! Evaluate spline fit with method from Numerical Recipes 880 ! (p. 86 Numerical Recipes, Press et al; 881 ! NOTE error in book for sign of "aa" term in derivative) 882 do ipts=1,npts 883 yy=rnorm(ipts)*rangem1 884 jj=1+int(yy*(n1xccc-1)) 885 diff=yy-(jj-1)*delta 886 bb = diff*deltam1 ; aa = one-bb 887 cc = aa*(aa**2-one)*delta2div6 888 dd = bb*(bb**2-one)*delta2div6 889 if (option==1.or.option==3) then 890 tcore(ipts)=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +& 891 & cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat) 892 end if 893 if (option>=2) then 894 dtcore(ipts)=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +& 895 & cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat) 896 end if 897 if (option==4) then 898 d2tcore(ipts)=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +& 899 & cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat) 900 end if 901 end do 902 end if 903 904 ! Now, perform the loop over selected grid points 905 do ipts=1,npts 906 ishift1=indx1(iindex(ipts)) 907 ishift2=indx2(iindex(ipts)) 908 difmag=rnorm(ipts) 909 910 rdiff1=rrdiff(ishift1,1);rdiff2=rrdiff(ishift2,2) 911 jpts=ii(ishift1,1)+n1*(ii(ishift2,2)-1+n2*(ffti3_local(i3)-1)) 912 913 ! === Evaluate charge density 914 if (option==1) then 915 xccc3d(jpts)=xccc3d(jpts)+tcore(ipts) 916 917 ! === Accumulate contributions to forces 918 else if (option==2) then 919 if (difmag>tol10) then 920 term=vxc_eff(jpts)*dtcore(ipts)/difmag 921 corgr(1)=corgr(1)+rdiff1*term 922 corgr(2)=corgr(2)+rdiff2*term 923 corgr(3)=corgr(3)+rdiff3*term 924 end if 925 926 ! === Accumulate contributions to stress tensor (in red. coordinates) 927 else if (option==3) then 928 if (difmag>tol10) then 929 term=vxc_eff(jpts)*dtcore(ipts)*rangem1/difmag/real(nfftot,dp) 930 ! Write out the 6 symmetric components 931 corfra(1,1)=corfra(1,1)+term*rdiff1*rdiff1 932 corfra(2,2)=corfra(2,2)+term*rdiff2*rdiff2 933 corfra(3,3)=corfra(3,3)+term*rdiff3*rdiff3 934 corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2 935 corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1 936 corfra(2,1)=corfra(2,1)+term*rdiff2*rdiff1 937 ! (the above still needs to be transformed to cartesian coords) 938 end if 939 ! Also compute a diagonal term 940 strdia=strdia+vxc_eff(jpts)*tcore(ipts) 941 942 ! === Compute frozen-wf contribution to Dynamical matrix 943 else if (option==4) then 944 tt(1)=rmet(1,1)*rdiff1+rmet(1,2)*rdiff2+rmet(1,3)*rdiff3 945 tt(2)=rmet(2,1)*rdiff1+rmet(2,2)*rdiff2+rmet(2,3)*rdiff3 946 tt(3)=rmet(3,1)*rdiff1+rmet(3,2)*rdiff2+rmet(3,3)*rdiff3 947 if (difmag>tol10) then 948 term=(ucvol/real(nfftot,dp))*vxc_eff(jpts)*rangem1/difmag 949 term1=term*tcore(ipts) 950 term2=term*d2tcore(ipts)*rangem1/difmag 951 do mu=1,3 952 do nu=1,3 953 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)& 954 & +(term2-term1/difmag**2)*tt(mu)*tt(nu)& 955 & +term1*rmet(mu,nu) 956 end do 957 end do 958 else 959 ! There is a contribution from difmag=zero ! 960 term=(ucvol/real(nfftot,dp))*vxc_eff(jpts)*d2tcore(ipts)*rangem1**2 961 do mu=1,3 962 do nu=1,3 963 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu) 964 end do 965 end do 966 end if 967 end if ! Choice of option 968 969 end do ! Loop on ipts (ishift1, ishift2) 970 971 end do ! Loop on ishift3 972 973 ABI_FREE(ii) 974 ABI_FREE(rrdiff) 975 ABI_FREE(indx1) 976 ABI_FREE(indx2) 977 ABI_FREE(iindex) 978 ABI_FREE(rnorm) 979 if (allocated(tcore)) then 980 ABI_FREE(tcore) 981 end if 982 if (allocated(tcoretau)) then 983 ABI_FREE(tcoretau) 984 end if 985 if (allocated(dtcore)) then 986 ABI_FREE(dtcore) 987 end if 988 if (allocated(d2tcore)) then 989 ABI_FREE(d2tcore) 990 end if 991 992 if (option==2) then 993 fact=-(ucvol/real(nfftot,dp))/range 994 grxc(:,iatom)=corgr(:)*fact 995 end if 996 997 ! End loop on atoms 998 end do 999 1000 if (usepaw==1) then 1001 call pawrad_free(core_mesh) 1002 end if 1003 1004 !End loop over atom types 1005 end do 1006 1007 if(option>=2.and.nspden>=2) then 1008 ABI_FREE(vxc_eff) 1009 end if 1010 1011 !Forces: translate into reduced coordinates 1012 if (option==2) then 1013 do iatom=1,natom 1014 t1=grxc(1,iatom);t2=grxc(2,iatom);t3=grxc(3,iatom) 1015 grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3 1016 end do 1017 end if 1018 1019 !Stress tensor: symmetrize, translate into cartesian coord., add diagonal part 1020 if (option==3) then 1021 corstr(1)=corfra(1,1) ; corstr(2)=corfra(2,2) 1022 corstr(3)=corfra(3,3) ; corstr(4)=corfra(3,2) 1023 corstr(5)=corfra(3,1) ; corstr(6)=corfra(2,1) 1024 call strconv(corstr,rprimd,corstr) 1025 corstr(1)=corstr(1)+strdia/real(nfftot,dp) 1026 corstr(2)=corstr(2)+strdia/real(nfftot,dp) 1027 corstr(3)=corstr(3)+strdia/real(nfftot,dp) 1028 end if 1029 1030 !If needed sum over MPI processes 1031 if(mpi_enreg%nproc_fft>1) then 1032 call timab(539,1,tsec) 1033 if (option==2) then 1034 call xmpi_sum(grxc,mpi_enreg%comm_fft,ier) 1035 end if 1036 if (option==3) then 1037 call xmpi_sum(corstr,mpi_enreg%comm_fft,ier) 1038 end if 1039 if (option==4) then 1040 call xmpi_sum(dyfrx2,mpi_enreg%comm_fft,ier) 1041 end if 1042 call timab(539,2,tsec) 1043 end if 1044 1045 call timab(12,2,tsec) 1046 1047 contains