TABLE OF CONTENTS


ABINIT/cross_mkcore [ Functions ]

[ Top ] [ Functions ]

NAME

  cross_mkcore

FUNCTION

  Define magnitude of cross product of two vectors

PARENTS

CHILDREN

SOURCE

563    function cross_mkcore(xx,yy,zz,aa,bb,cc)
564 
565 
566 !This section has been created automatically by the script Abilint (TD).
567 !Do not modify the following lines by hand.
568 #undef ABI_FUNC
569 #define ABI_FUNC 'cross_mkcore'
570 !End of the abilint section
571 
572    real(dp) :: cross_mkcore
573    real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
574 ! *************************************************************************
575    cross_mkcore=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
576  end function cross_mkcore
577 
578 end subroutine mkcore

ABINIT/cross_mkcore_alt [ Functions ]

[ Top ] [ Functions ]

NAME

  cross_mkcore_alt

FUNCTION

  Define magnitude of cross product of two vectors

PARENTS

CHILDREN

SOURCE

1099    function cross_mkcore_alt(xx,yy,zz,aa,bb,cc)
1100 
1101 
1102 !This section has been created automatically by the script Abilint (TD).
1103 !Do not modify the following lines by hand.
1104 #undef ABI_FUNC
1105 #define ABI_FUNC 'cross_mkcore_alt'
1106 !End of the abilint section
1107 
1108     real(dp) :: cross_mkcore_alt
1109     real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
1110 ! *************************************************************************
1111    cross_mkcore_alt=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
1112  end function cross_mkcore_alt

ABINIT/dfpt_mkcore [ Functions ]

[ Top ] [ 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

PARENTS

      dfpt_dyxc1,dfpt_eltfrxc,dfpt_looppert,dfpt_nselt,dfpt_nstdy,dfpt_nstpaw
      dfptnl_loop

CHILDREN

SOURCE

1204 subroutine dfpt_mkcore(cplex,idir,ipert,natom,ntypat,n1,n1xccc,&
1205 & n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xccc3d1,xred)
1206 
1207 
1208 !This section has been created automatically by the script Abilint (TD).
1209 !Do not modify the following lines by hand.
1210 #undef ABI_FUNC
1211 #define ABI_FUNC 'dfpt_mkcore'
1212 !End of the abilint section
1213 
1214  implicit none
1215 
1216 !Arguments ------------------------------------
1217 !scalars
1218  integer,intent(in) :: cplex,idir,ipert,n1,n1xccc,n2,n3,natom,ntypat
1219  real(dp),intent(in) :: ucvol
1220 !arrays
1221  integer,intent(in) :: typat(natom)
1222  real(dp),intent(in) :: qphon(3),rprimd(3,3),xccc1d(n1xccc,6,ntypat)
1223  real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom)
1224  real(dp),intent(out) :: xccc3d1(cplex*n1*n2*n3)
1225 
1226 !Local variables-------------------------------
1227 !scalars
1228  integer,parameter :: mshift=401
1229  integer :: i1,i2,i3,iatom,ifft,ishift,ishift1,ishift2,ishift3,istr
1230  integer :: ixp,jj,ka,kb,mrange,mu,nu
1231  real(dp) :: aa,bb,cc,dd,delta,delta2div6,deltam1,diff,difmag
1232  real(dp) :: difmag2,difmag2_fact,difmag2_part,func,phase,phi,phr,prod
1233  real(dp) :: range,range2,rangem1,rdiff1,rdiff2,rdiff3,term
1234  real(dp) :: yy
1235  character(len=500) :: message
1236 !arrays
1237  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1238  integer :: igrid(3),irange(3),ngfft(3)
1239  integer,allocatable :: ii(:,:)
1240  real(dp) :: drmetds(3,3),lencp(3),rmet(3,3),scale(3),tau(3)
1241  real(dp),allocatable :: rrdiff(:,:)
1242 
1243 ! *************************************************************************
1244 
1245 ! if( ipert<1 .or. ipert> natom+7) then
1246 !   write(message,'(a,i0,a,a,a,i0,a)')&
1247 !&   ' The argument ipert must be between 1 and natom+7=',natom+7,',',ch10,&
1248 !&   ' while it is ipert=',ipert,'.'
1249 !   MSG_BUG(message)
1250 ! end if
1251 
1252  if( (ipert==natom+3 .or. ipert==natom+4) .and. cplex/=1) then
1253    write(message,'(3a,i4,a)')&
1254 &   'The argument cplex must be 1 for strain perturbationh',ch10,&
1255 &   'while it is cplex=',cplex,'.'
1256    MSG_BUG(message)
1257  end if
1258 
1259 !Zero out array
1260  xccc3d1(:)=0.0_dp
1261 
1262 !For a non-linear XC core correction, the perturbation must be phonon-type or strain type
1263  if(ipert<=natom .or. ipert==natom+3 .or. ipert==natom+4) then
1264 
1265    if( idir<1 .or. idir> 3) then
1266      write(message,'(a,a,a,i4,a)')&
1267 &     'The argument idir must be between 1 and 3,',ch10,&
1268 &     'while it is idir=',idir,'.'
1269      MSG_BUG(message)
1270    end if
1271 
1272 !  Compute lengths of cross products for pairs of primitive
1273 !  translation vectors (used in setting index search range below)
1274    lencp(1)=cross_mk(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
1275 &   rprimd(1,3),rprimd(2,3),rprimd(3,3))
1276    lencp(2)=cross_mk(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
1277 &   rprimd(1,1),rprimd(2,1),rprimd(3,1))
1278    lencp(3)=cross_mk(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
1279 &   rprimd(1,2),rprimd(2,2),rprimd(3,2))
1280 
1281 !  Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
1282 !  (recall ucvol=R1.(R2xR3))
1283    scale(:)=ucvol/lencp(:)
1284 
1285 !  Compute metric tensor in real space rmet
1286    do nu=1,3
1287      rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+&
1288 &     rprimd(3,:)*rprimd(3,nu)
1289    end do
1290 
1291 !  Section to be executed only for strain perturbation
1292 !  Compute derivative of metric tensor wrt strain component istr
1293    if(ipert==natom+3 .or. ipert==natom+4) then
1294      istr=idir + 3*(ipert-natom-3)
1295 
1296      ka=idx(2*istr-1);kb=idx(2*istr)
1297      do jj = 1,3
1298        drmetds(:,jj)=(rprimd(ka,:)*rprimd(kb,jj)+rprimd(kb,:)*rprimd(ka,jj))
1299      end do
1300 !    For historical reasons:
1301      drmetds(:,:)=0.5_dp*drmetds(:,:)
1302 
1303 !    end of strain perturbation section
1304    end if
1305 
1306    ngfft(1)=n1
1307    ngfft(2)=n2
1308    ngfft(3)=n3
1309 
1310    delta=1.0_dp/(n1xccc-1)
1311    deltam1=n1xccc-1
1312    delta2div6=delta**2/6.0_dp
1313 
1314 !  Loop over atoms in unit cell
1315 !  Note that we cycle immediately for all except the displaced atom
1316 !  for such a perturbation.  The loop is executed over all the
1317 !  atoms for a strain peturbation.
1318    do iatom=1,natom
1319      if(ipert<=natom .and. iatom/=ipert) cycle
1320 !    Set search range (density cuts off perfectly beyond range)
1321 !    Cycle if no range.
1322      range=0.0_dp
1323      range=xcccrc(typat(iatom))
1324      if(range<1.d-16) cycle
1325 
1326      range2=range**2
1327      rangem1=1.0_dp/range
1328 
1329 !    compute mrange for ii(:,3), rrdiff(:,3), inserted by MM (2005/12/06)
1330 !    Consider each component in turn : compute range
1331      do mu=1,3
1332 
1333 !      Convert reduced coord of given atom to [0,1)
1334        tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp)
1335 
1336 !      Use tau to find nearest grid point along R(mu)
1337 !      (igrid=0 is the origin; shift by 1 to agree with usual index)
1338        igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
1339 
1340 !      Use range to compute an index range along R(mu)
1341 !      (add 1 to make sure it covers full range)
1342        irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
1343 
1344      end do
1345 
1346 !    Allocate arrays that depends on the range
1347      mrange=maxval(irange(1:3))
1348      ABI_ALLOCATE(ii,(2*mrange+1,3))
1349      ABI_ALLOCATE(rrdiff,(2*mrange+1,3))
1350 
1351 !    Consider each component in turn
1352      do mu=1,3
1353 
1354 !      temporarily suppressed by MM (2005/12/02)
1355 !      Convert reduced coord of given atom to [0,1)
1356 !      tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp)
1357 
1358 !      Use tau to find nearest grid point along R(mu)
1359 !      (igrid=0 is the origin; shift by 1 to agree with usual index)
1360 !      igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
1361 
1362 !      Use range to compute an index range along R(mu)
1363 !      (add 1 to make sure it covers full range)
1364 !      irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
1365 
1366 !      Check that the largest range is smallest than the maximum
1367 !      allowed one
1368 !      if(2*irange(mu)+1 > mshift)then
1369 !      write(message, '(a,a,a,a,i6,a)' ) ch10,&
1370 !      &    ' dfpt_mkcore : BUG -',ch10,&
1371 !      &    '  The range around atom',iatom,' is too large.'
1372 !      MSG_BUG(message)
1373 !      end if
1374 
1375 !      Set up a counter that explore the relevant range
1376 !      of points around the atom
1377        ishift=0
1378        do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
1379          ishift=ishift+1
1380          ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
1381          rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu)
1382        end do
1383 
1384 !      End loop on mu
1385      end do
1386 
1387 !    Conduct triple loop over restricted range of grid points for iatom
1388 
1389      do ishift3=1,1+2*irange(3)
1390 !      map back to [1,ngfft(3)] for usual fortran index in unit cell
1391        i3=ii(ishift3,3)
1392 !      find vector from atom location to grid point (reduced)
1393        rdiff3=rrdiff(ishift3,3)
1394 
1395        do ishift2=1,1+2*irange(2)
1396          i2=ii(ishift2,2)
1397          rdiff2=rrdiff(ishift2,2)
1398 !        Prepare the computation of difmag2
1399          difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
1400 &         +2.0_dp*rmet(3,2)*rdiff3*rdiff2
1401          difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
1402 
1403          do ishift1=1,1+2*irange(1)
1404            rdiff1=rrdiff(ishift1,1)
1405 
1406 !          Compute (rgrid-tau-Rprim)**2
1407            difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
1408 
1409 !          Only accept contribution inside defined range
1410            if (difmag2<range2) then
1411 
1412 !            Prepare computation of core charge function and derivative,
1413 !            using splines
1414              difmag=sqrt(difmag2)
1415              if (difmag>=1.0d-10) then
1416                i1=ii(ishift1,1)
1417                yy=difmag*rangem1
1418 
1419 !              Compute index of yy over 1 to n1xccc scale
1420                jj=1+int(yy*(n1xccc-1))
1421                diff=yy-(jj-1)*delta
1422 
1423 !              Will evaluate spline fit (p. 86 Numerical Recipes, Press et al;
1424 !              NOTE error in book for sign of "aa" term in derivative;
1425 !              also see splfit routine).
1426                bb = diff*deltam1
1427                aa = 1.0_dp-bb
1428                cc = aa*(aa**2-1.0_dp)*delta2div6
1429                dd = bb*(bb**2-1.0_dp)*delta2div6
1430 
1431 !              Evaluate spline fit of 1st der of core charge density
1432 !              from xccc1d(:,2,:) and (:,4,:)
1433                func=aa*xccc1d(jj,2,typat(iatom))+bb*xccc1d(jj+1,2,typat(iatom)) +&
1434 &               cc*xccc1d(jj,4,typat(iatom))+dd* xccc1d(jj+1,4,typat(iatom))
1435 
1436                if(ipert<=natom) then
1437                  phase=2*pi*(qphon(1)*(rdiff1+xred(1,iatom))  &
1438 &                 +qphon(2)*(rdiff2+xred(2,iatom))  &
1439 &                 +qphon(3)*(rdiff3+xred(3,iatom)))
1440                  prod=rmet(idir,1)*rdiff1+rmet(idir,2)*rdiff2+rmet(idir,3)*rdiff3
1441 
1442                  term=-func*rangem1/difmag*prod
1443                  ifft=i1+n1*(i2-1+n2*(i3-1))
1444                  phr=cos(phase)
1445                  if(cplex==1)then
1446                    xccc3d1(ifft)=xccc3d1(ifft)+term*phr
1447                  else
1448                    phi=sin(phase)
1449                    xccc3d1(2*ifft-1)=xccc3d1(2*ifft-1)+term*phr
1450                    xccc3d1(2*ifft  )=xccc3d1(2*ifft  )-term*phi
1451                  end if
1452                else
1453                  prod=&
1454 &                 (rdiff1*(drmetds(1,1)*rdiff1+drmetds(1,2)*rdiff2+drmetds(1,3)*rdiff3)&
1455 &                 +rdiff2*(drmetds(2,1)*rdiff1+drmetds(2,2)*rdiff2+drmetds(2,3)*rdiff3)&
1456 &                 +rdiff3*(drmetds(3,1)*rdiff1+drmetds(3,2)*rdiff2+drmetds(3,3)*rdiff3))
1457                  term=prod*func*rangem1/difmag
1458 
1459                  ifft=i1+n1*(i2-1+n2*(i3-1))
1460                  xccc3d1(ifft)=xccc3d1(ifft)+term
1461 
1462                end if
1463 
1464 !              End of the condition for the distance not to vanish
1465              end if
1466 
1467 !            End of condition to be inside the range
1468            end if
1469 
1470 !          End loop on ishift1
1471          end do
1472 
1473 !        End loop on ishift2
1474        end do
1475 
1476 !      End loop on ishift3
1477      end do
1478 
1479      ABI_DEALLOCATE(ii)
1480      ABI_DEALLOCATE(rrdiff)
1481 !    End loop on atoms
1482    end do
1483 
1484 !  End of the condition ipert corresponds to a phonon type perturbation
1485 !  or strain type perturbation
1486  end if
1487 
1488  contains
1489 
1490    function cross_mk(xx,yy,zz,aa,bb,cc)
1491 
1492 
1493 !This section has been created automatically by the script Abilint (TD).
1494 !Do not modify the following lines by hand.
1495 #undef ABI_FUNC
1496 #define ABI_FUNC 'cross_mk'
1497 !End of the abilint section
1498 
1499    real(dp) :: cross_mk
1500    real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
1501    cross_mk=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
1502  end function cross_mk
1503 
1504 end subroutine dfpt_mkcore

ABINIT/indpos_mkcore_alt [ Functions ]

[ Top ] [ 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

1132    subroutine indpos_mkcore_alt(periodic,ii,nn,jj,inside)
1133 !    Find the grid index of a given position in the cell according to the BC
1134 !    Determine also whether the index is inside or outside the box for free BC
1135 
1136 !This section has been created automatically by the script Abilint (TD).
1137 !Do not modify the following lines by hand.
1138 #undef ABI_FUNC
1139 #define ABI_FUNC 'indpos_mkcore_alt'
1140 !End of the abilint section
1141 
1142     integer, intent(in) :: ii,nn
1143     integer, intent(out) :: jj
1144     logical, intent(in) :: periodic
1145     logical, intent(out) :: inside
1146 ! *************************************************************************
1147    if (periodic) then
1148      inside=.true. ; jj=modulo(ii-1,nn)+1
1149    else
1150      jj=ii ; inside=(ii>=1.and.ii<=nn)
1151    end if
1152  end subroutine indpos_mkcore_alt
1153 
1154 end subroutine mkcore_alt

ABINIT/m_mkcore [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mkcore

FUNCTION

  Routines related to non-linear core correction.

COPYRIGHT

  Copyright (C) 1998-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_mkcore
28 
29  use defs_basis
30  use defs_abitypes
31  use m_abicore
32  use m_xmpi
33  use m_errors
34  use m_linalg_interfaces
35 
36  use m_geometry,   only : strconv
37  use m_time,       only : timab
38  use m_mpinfo,     only : ptabs_fourdp
39  use m_sort,        only : sort_dp
40  use m_pawrad,      only : pawrad_type, pawrad_init, pawrad_free
41  use m_pawtab,      only : pawtab_type
42  use m_paw_numeric, only : paw_splint
43 
44  implicit none
45 
46  private

ABINIT/mkcore [ Functions ]

[ Top ] [ 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

PARENTS

      dfpt_dyfro,forces,nonlinear,prcref,prcref_PMA,respfn,setvtr,stress
      xchybrid_ncpp_cc

CHILDREN

SOURCE

111 subroutine mkcore(corstr,dyfrx2,grxc,mpi_enreg,natom,nfft,nspden,ntypat,n1,n1xccc,&
112 & n2,n3,option,rprimd,typat,ucvol,vxc,xcccrc,xccc1d,xccc3d,xred)
113 
114 
115 !This section has been created automatically by the script Abilint (TD).
116 !Do not modify the following lines by hand.
117 #undef ABI_FUNC
118 #define ABI_FUNC 'mkcore'
119 !End of the abilint section
120 
121  implicit none
122 
123 !Arguments ------------------------------------
124 !scalars
125  integer,intent(in) :: n1,n1xccc,n2,n3,natom,nfft,nspden,ntypat,option
126  real(dp),intent(in) :: ucvol
127  type(mpi_type),intent(in) :: mpi_enreg
128 !arrays
129  integer,intent(in) :: typat(natom)
130  real(dp),intent(in) :: rprimd(3,3),vxc(nfft,nspden),xccc1d(n1xccc,6,ntypat)
131  real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom)
132  real(dp),intent(inout) :: xccc3d(nfft)
133  real(dp),intent(out) :: corstr(6),dyfrx2(3,3,natom)
134  real(dp),intent(inout) :: grxc(3,natom)
135 
136 !Local variables-------------------------------
137 !scalars
138  integer :: i1,i2,i3,iatom,ier,ifft,ishift,ishift1,ishift2
139  integer :: ishift3,itypat,ixp,jj,me_fft,mrange,mu,nfftot,nu
140  real(dp) :: dd,delta,delta2div6,deltam1,diff,difmag,difmag2
141  real(dp) :: difmag2_fact,difmag2_part,fact,func,grxc1,grxc2,grxc3,range,range2
142  real(dp) :: rangem1,rdiff1,rdiff2,rdiff3,strdia,t1,t2,t3,term,term1,term2
143  character(len=500) :: message
144 !arrays
145  integer :: igrid(3),irange(3),ngfft(3)
146  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
147  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
148  integer,allocatable :: ii(:,:)
149  real(dp) :: yy,aa,bb,cc
150  real(dp) :: corfra(3,3),lencp(3),rmet(3,3),scale(3),tau(3),tsec(2),tt(3)
151  real(dp),allocatable :: rrdiff(:,:),work(:,:,:)
152 
153 !************************************************************************
154 
155  call timab(12,1,tsec)
156 
157 !Make sure option is acceptable
158  if (option<0.or.option>4) then
159    write(message, '(a,i12,a,a,a)' )&
160 &   'option=',option,' is not allowed.',ch10,&
161 &   'Must be 1, 2, 3 or 4.'
162    MSG_BUG(message)
163  end if
164 
165 !Zero out only the appropriate array according to option:
166 !others are dummies with no storage
167 
168  if (option==1) then
169 !  Zero out array to permit accumulation over atom types below:
170    xccc3d(:)=zero
171  else if (option==2) then
172 !  Zero out gradient of Exc array
173    grxc(:,:)=zero
174  else if (option==3) then
175 !  Zero out locally defined stress array
176    corfra(:,:)=zero
177    strdia=zero
178  else if (option==4) then
179 !  Zero out fr-wf part of the dynamical matrix
180    dyfrx2(:,:,:)=zero
181  else
182    MSG_BUG(" Can't be here! (bad option)")
183  end if
184 
185 !Compute lengths of cross products for pairs of primitive
186 !translation vectors (used in setting index search range below)
187  lencp(1)=cross_mkcore(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
188 & rprimd(1,3),rprimd(2,3),rprimd(3,3))
189  lencp(2)=cross_mkcore(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
190 & rprimd(1,1),rprimd(2,1),rprimd(3,1))
191  lencp(3)=cross_mkcore(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
192 & rprimd(1,2),rprimd(2,2),rprimd(3,2))
193 
194 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
195 !(recall ucvol=R1.(R2xR3))
196  scale(:)=ucvol/lencp(:)
197 
198 !Compute metric tensor in real space rmet
199  do nu=1,3
200    rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu)
201  end do
202 
203  ngfft(1)=n1
204  ngfft(2)=n2
205  ngfft(3)=n3
206  nfftot=n1*n2*n3
207  me_fft = mpi_enreg%me_fft
208 
209  ! Get the distrib associated with this fft_grid
210  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
211 
212  delta=one/(n1xccc-1)
213  deltam1=n1xccc-1
214  delta2div6=delta**2/6.0d0
215 
216  if (option>=2) then
217    ABI_ALLOCATE(work,(n1,n2,n3))
218 !  For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down))
219 !  For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22})
220    if (nspden>=2) then
221      ifft=1
222      do i3=1,n3
223        if(me_fft==fftn3_distrib(i3)) then
224          do i2=1,n2
225            do i1=1,n1
226              work(i1,i2,i3)=half*(vxc(ifft,1)+vxc(ifft,2))
227              ifft=ifft+1
228            end do
229          end do
230        end if
231      end do
232    else
233      ifft=1
234      do i3=1,n3
235        if(me_fft==fftn3_distrib(i3)) then
236          do i2=1,n2
237            do i1=1,n1
238              work(i1,i2,i3)=vxc(ifft,1)
239              ifft=ifft+1
240            end do
241          end do
242        end if
243      end do
244 !    call DCOPY(nfft,vxc,1,work,1)
245    end if
246  end if
247 
248 !Loop over atoms in unit cell
249  do iatom=1,natom
250 
251    if(option==2)then
252      grxc1=zero
253      grxc2=zero
254      grxc3=zero
255    end if
256 
257 !  Set search range (density cuts off perfectly beyond range)
258    itypat=typat(iatom)
259    range=xcccrc(itypat)
260 
261 !  Skip loop if this atom has no core charge
262    if (abs(range)<1.d-16) cycle
263 
264    range2=range**2
265    rangem1=one/range
266 
267 !  Consider each component in turn : compute range
268    do mu=1,3
269 
270 !    Convert reduced coord of given atom to [0,1)
271      tau(mu)=mod(xred(mu,iatom)+one-aint(xred(mu,iatom)),one)
272 
273 !    Use tau to find nearest grid point along R(mu)
274 !    (igrid=0 is the origin; shift by 1 to agree with usual index)
275      igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
276 
277 !    Use range to compute an index range along R(mu)
278 !    (add 1 to make sure it covers full range)
279      irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
280 
281    end do
282 
283 !  Allocate arrays that depends on the range
284    mrange=maxval(irange(1:3))
285    ABI_ALLOCATE(ii,(2*mrange+1,3))
286    ABI_ALLOCATE(rrdiff,(2*mrange+1,3))
287 
288 !  Set up counters that explore the relevant range
289 !  of points around the atom
290    do mu=1,3
291      ishift=0
292      do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
293        ishift=ishift+1
294        ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
295        rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu)
296      end do
297    end do
298 
299 !  Conduct triple loop over restricted range of grid points for iatom
300    do ishift3=1,1+2*irange(3)
301 !    map back to [1,ngfft(3)] for usual fortran index in unit cell
302      i3=ii(ishift3,3)
303      if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle
304 !    find vector from atom location to grid point (reduced)
305      rdiff3=rrdiff(ishift3,3)
306 
307      do ishift2=1,1+2*irange(2)
308        i2=ii(ishift2,2)
309        rdiff2=rrdiff(ishift2,2)
310 !      Prepare the computation of difmag2
311        difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
312 &       +2.0d0*rmet(3,2)*rdiff3*rdiff2
313        difmag2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
314 
315        do ishift1=1,1+2*irange(1)
316          rdiff1=rrdiff(ishift1,1)
317 
318 !        Compute (rgrid-tau-Rprim)**2
319          difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
320 
321 !        Only accept contribution inside defined range
322          if (difmag2<range2-tol12) then
323 
324 !          Prepare computation of core charge function and derivative,
325 !          using splines
326            i1=ii(ishift1,1)
327            difmag=sqrt(difmag2)
328            yy=difmag*rangem1
329 
330 !          Compute index of yy over 1 to n1xccc scale
331            jj=1+int(yy*(n1xccc-1))
332            diff=yy-(jj-1)*delta
333 
334 !          Will evaluate spline fit (p. 86 Numerical Recipes, Press et al;
335 !          NOTE error in book for sign of "aa" term in derivative;
336 !          also see splfit routine).
337            bb = diff*deltam1
338            aa = one-bb
339            cc = aa*(aa**2-one)*delta2div6
340            dd = bb*(bb**2-one)*delta2div6
341 
342 
343 !          Test first for option 2, the most frequently used
344            if (option==2) then
345 
346 !            Accumulate contributions to Exc gradients
347 
348              if (difmag>1.0d-10) then
349 
350 !              Evaluate spline fit of 1st der of core charge density
351 !              from xccc1d(:,2,:) and (:,4,:)
352                func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +&
353 &               cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
354                term=work(i1,i2,i3)*func/difmag
355                grxc1=grxc1+rdiff1*term
356                grxc2=grxc2+rdiff2*term
357                grxc3=grxc3+rdiff3*term
358              end if
359 
360            else if (option==1) then
361 
362 !            Evaluate spline fit of core charge density
363 !            from xccc1d(:,1,:) and (:,3,:)
364              func=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +&
365 &             cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat)
366 
367 !            Accumulate contributions to core electron density
368 !            throughout unit cell
369              ifft=i1+n1*(i2-1+n2*(ffti3_local(i3)-1))
370              xccc3d(ifft)=xccc3d(ifft)+func
371 
372            else if (option==3) then
373 
374 !            Accumulate contributions to stress tensor
375 !            in reduced coordinates
376 
377              if (difmag>1.0d-10) then
378 
379 !              Evaluate spline fit of 1st der of core charge density
380 !              from xccc1d(:,2,:) and (:,4,:)
381                func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +&
382 &               cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
383                term=work(i1,i2,i3)*func*rangem1/difmag/dble(n1*n2*n3)
384 !              Write out the 6 symmetric components
385                corfra(1,1)=corfra(1,1)+term*rdiff1**2
386                corfra(2,2)=corfra(2,2)+term*rdiff2**2
387                corfra(3,3)=corfra(3,3)+term*rdiff3**2
388                corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2
389                corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1
390                corfra(2,1)=corfra(2,1)+term*rdiff2*rdiff1
391 !              (the above still needs to be transformed to cartesian coords)
392 
393              end if
394 
395 !            Also compute a diagonal term
396 !            Evaluate spline fit of core charge density
397 !            from xccc1d(:,1,:) and (:,3,:)
398              func=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +&
399 &             cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat)
400              strdia=strdia+work(i1,i2,i3)*func
401 
402            else if (option==4) then
403 
404 !            Compute frozen-wf contribution to Dynamical matrix
405 
406              tt(1)=rmet(1,1)*rdiff1+rmet(1,2)*rdiff2+rmet(1,3)*rdiff3
407              tt(2)=rmet(2,1)*rdiff1+rmet(2,2)*rdiff2+rmet(2,3)*rdiff3
408              tt(3)=rmet(3,1)*rdiff1+rmet(3,2)*rdiff2+rmet(3,3)*rdiff3
409 
410              if (difmag>1.d-10) then
411 
412 !              Accumulate contributions to dynamical matrix
413                term=(ucvol/dble(nfftot))*work(i1,i2,i3)*rangem1/difmag
414 !              Evaluate spline fit of 1st der of core charge density
415 !              from xccc1d(:,2,:) and (:,4,:)
416                func=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +&
417 &               cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
418                term1=term*func
419 !              Evaluate spline fit of 2nd der of core charge density
420 !              from xccc1d(:,3,:) and (:,5,:)
421                func=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +&
422 &               cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat)
423                term2=term*func*rangem1/difmag
424                do mu=1,3
425                  do nu=1,3
426                    dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)&
427 &                   +(term2-term1/difmag**2)*tt(mu)*tt(nu)&
428 &                   +term1*rmet(mu,nu)
429                  end do
430                end do
431 
432              else
433 
434 !              There is a contribution from difmag=zero !
435 !              Evaluate spline fit of 2nd der of core charge density
436 !              from xccc1d(:,3,:) and (:,5,:)
437                func=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +&
438 &               cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat)
439                term=(ucvol/dble(nfftot))*work(i1,i2,i3)*func*rangem1**2
440                do mu=1,3
441                  do nu=1,3
442                    dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu)
443                  end do
444                end do
445 
446 !              End of condition not to be precisely on the point (difmag=zero)
447              end if
448 
449 !            If option is not 1, 2, 3, or 4.
450            else
451              MSG_BUG("Can't be here in mkcore")
452 !            End of choice of option
453            end if
454 
455 !          End of condition on the range
456          end if
457 
458 !        End loop on ishift1
459        end do
460 
461 !      End loop on ishift2
462      end do
463 
464 !    End loop on ishift3
465    end do
466 
467    ABI_DEALLOCATE(ii)
468    ABI_DEALLOCATE(rrdiff)
469 
470    if(option==2)then
471      fact=-(ucvol/dble(nfftot))/range
472      grxc(1,iatom)=grxc1*fact
473      grxc(2,iatom)=grxc2*fact
474      grxc(3,iatom)=grxc3*fact
475    end if
476 
477 !  End big loop on atoms
478  end do
479 
480  if (option==2) then
481 
482 !  Apply rmet as needed to get reduced coordinate gradients
483    do iatom=1,natom
484      t1=grxc(1,iatom)
485      t2=grxc(2,iatom)
486      t3=grxc(3,iatom)
487      grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3
488 
489    end do
490  end if
491 
492  if (option==3) then
493 
494 !  Transform stress tensor from full storage mode to symmetric storage mode
495    corstr(1)=corfra(1,1)
496    corstr(2)=corfra(2,2)
497    corstr(3)=corfra(3,3)
498    corstr(4)=corfra(3,2)
499    corstr(5)=corfra(3,1)
500    corstr(6)=corfra(2,1)
501 
502 !  Transform stress tensor from reduced coordinates to cartesian coordinates
503    call strconv(corstr,rprimd,corstr)
504 
505 !  Compute diagonal contribution to stress tensor (need input xccc3d)
506 !  strdia = (1/N) Sum(r) [mu_xc_avg(r) * rho_core(r)]
507    ifft=0 ; strdia=zero
508    do i3=1,n3
509      if(me_fft==fftn3_distrib(i3)) then
510        do i2=1,n2
511          do i1=1,n1
512            ifft=ifft+1
513            strdia=strdia+work(i1,i2,i3)*xccc3d(ifft)
514          end do
515        end do
516      end if
517    end do
518    strdia=strdia/dble(nfftot)
519 !  strdia=DDOT(nfft,work,1,xccc3d,1)/dble(nfftot)
520 
521 !  Add diagonal term to stress tensor
522    corstr(1)=corstr(1)+strdia
523    corstr(2)=corstr(2)+strdia
524    corstr(3)=corstr(3)+strdia
525  end if
526 
527  if(option>=2)  then
528    ABI_DEALLOCATE(work)
529  end if
530 
531  if(mpi_enreg%nproc_fft > 1) then
532    call timab(539,1,tsec)
533    if(option==2) then
534      call xmpi_sum(grxc,mpi_enreg%comm_fft,ier)
535    end if
536    if(option==3) then
537      call xmpi_sum(corstr,mpi_enreg%comm_fft,ier)
538    end if
539    if(option==4) then
540      call xmpi_sum(dyfrx2,mpi_enreg%comm_fft,ier)
541    end if
542    call timab(539,2,tsec)
543  end if
544 
545  call timab(12,2,tsec)
546 
547  contains

ABINIT/mkcore_alt [ Functions ]

[ Top ] [ 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=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

 644 subroutine mkcore_alt(atindx1,corstr,dyfrx2,grxc,icoulomb,mpi_enreg,natom,nfft,nspden,&
 645 &          nattyp,ntypat,n1,n1xccc,n2,n3,option,rprimd,ucvol,vxc,xcccrc,xccc1d,&
 646 &          xccc3d,xred,pawrad,pawtab,usepaw)
 647 
 648 
 649 !This section has been created automatically by the script Abilint (TD).
 650 !Do not modify the following lines by hand.
 651 #undef ABI_FUNC
 652 #define ABI_FUNC 'mkcore_alt'
 653 !End of the abilint section
 654 
 655  implicit none
 656 
 657 !Arguments ------------------------------------
 658 !scalars
 659  integer,intent(in) :: icoulomb,n1,n1xccc,n2,n3,natom,nfft,nspden,ntypat,option,usepaw
 660  real(dp),intent(in) :: ucvol
 661  type(mpi_type),intent(in) :: mpi_enreg
 662  type(pawrad_type),intent(in) :: pawrad(:)
 663  type(pawtab_type),intent(in) :: pawtab(:)
 664 !arrays
 665  integer,intent(in) :: atindx1(natom),nattyp(ntypat)
 666  real(dp),intent(in) :: rprimd(3,3),xccc1d(n1xccc,6,ntypat)
 667  real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom)
 668  real(dp),intent(in),target :: vxc(nfft,nspden)
 669  real(dp),intent(out) :: corstr(6),grxc(3,natom),dyfrx2(3,3,natom)
 670  real(dp),intent(inout) :: xccc3d(nfft)
 671 
 672 !Local variables-------------------------------
 673 !scalars
 674  integer :: i1,i2,i3,iat,iatm,iatom,ier,ipts
 675  integer :: ishift,ishift1,ishift2,ishift3
 676  integer :: itypat,ixp,jj,jpts,me_fft,mrange,mu
 677  integer :: nfftot,npts,npts12,nu
 678  logical :: letsgo
 679  real(dp) :: aa,bb,cc,dd,delta,delta2div6,deltam1
 680  real(dp) :: diff,difmag,fact,range,range2
 681  real(dp) :: rangem1,rdiff1,rdiff2,rdiff3
 682  real(dp) :: rnorm2,rnorm2_fact,rnorm2_part
 683  real(dp) :: strdia,t1,t2,t3,term,term1,term2,yy
 684  character(len=1) :: geocode
 685  character(len=500) :: message
 686  type(pawrad_type) :: core_mesh
 687 !arrays
 688  integer :: igrid(3),irange(3),ishiftmax(3),ngfft(3)
 689  integer,allocatable :: ii(:,:),iindex(:),indx1(:),indx2(:)
 690  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
 691  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
 692  logical :: per(3)
 693  real(dp) :: corfra(3,3),corgr(3),lencp(3),rmet(3,3)
 694  real(dp) :: scale(3),tau(3),tsec(2),tt(3)
 695  real(dp),allocatable :: dtcore(:),d2tcore(:),rnorm(:)
 696  real(dp),allocatable :: rrdiff(:,:),tcore(:)
 697  real(dp), ABI_CONTIGUOUS pointer :: vxc_eff(:)
 698 
 699 !************************************************************************
 700 
 701  call timab(12,1,tsec)
 702 
 703 !Make sure option is acceptable
 704  if (option<0.or.option>4) then
 705    write(message, '(a,i12,a,a,a)' )&
 706 &   'option=',option,' is not allowed.',ch10,&
 707 &   'Must be 1, 2, 3 or 4.'
 708    MSG_BUG(message)
 709  end if
 710 
 711 !Zero out only the appropriate array according to option:
 712  if (option==1) then
 713    xccc3d(:)=zero
 714  else if (option==2) then
 715    grxc(:,:)=zero
 716  else if (option==3) then
 717    corfra(:,:)=zero
 718    strdia=zero
 719  else if (option==4) then
 720    dyfrx2(:,:,:)=zero
 721  end if
 722 
 723 !Conditions for periodicity in the three directions
 724  geocode='P'
 725  if (icoulomb==1) geocode='F'
 726  if (icoulomb==2) geocode='S'
 727  per(1)=(geocode /= 'F')
 728  per(2)=(geocode == 'P')
 729  per(3)=(geocode /= 'F')
 730 
 731 !Compute lengths of cross products for pairs of primitive
 732 !translation vectors (used in setting index search range below)
 733  lencp(1)=cross_mkcore_alt(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
 734 & rprimd(1,3),rprimd(2,3),rprimd(3,3))
 735  lencp(2)=cross_mkcore_alt(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
 736 & rprimd(1,1),rprimd(2,1),rprimd(3,1))
 737  lencp(3)=cross_mkcore_alt(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
 738 & rprimd(1,2),rprimd(2,2),rprimd(3,2))
 739 
 740 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
 741 !(recall ucvol=R1.(R2xR3))
 742  scale(:)=ucvol/lencp(:)
 743 
 744 !Compute metric tensor in real space rmet
 745  do nu=1,3
 746    rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+rprimd(3,:)*rprimd(3,nu)
 747  end do
 748 
 749 !Get the distrib associated with this fft_grid
 750  ngfft(1)=n1;ngfft(2)=n2;ngfft(3)=n3
 751  nfftot=n1*n2*n3 ; me_fft=mpi_enreg%me_fft
 752  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
 753 
 754  delta=one/(n1xccc-1)
 755  deltam1=n1xccc-1
 756  delta2div6=delta**2/6.0_dp
 757 
 758  if (option>=2) then
 759 !  For spin-polarization, replace vxc by (1/2)*(vxc(up)+vxc(down))
 760 !  For non-collinear magnetism, replace vxc by (1/2)*(vxc^{11}+vxc^{22})
 761    if (nspden>=2) then
 762      ABI_ALLOCATE(vxc_eff,(nfft))
 763      do jj=1,nfft
 764        vxc_eff(jj)=half*(vxc(jj,1)+vxc(jj,2))
 765      end do
 766    else
 767      vxc_eff => vxc(1:nfft,1)
 768    end if
 769  end if
 770 
 771 !Loop over atom types
 772  iatm=0
 773  do itypat=1,ntypat
 774 
 775 !  Set search range (density cuts off perfectly beyond range)
 776    range=xcccrc(itypat);if (usepaw==1) range=pawtab(itypat)%rcore
 777    range2=range**2 ; rangem1=one/range
 778 
 779 !  Skip loop if this type has no core charge
 780    if (abs(range)<1.d-16) cycle
 781 
 782 !  PAW: create mesh for core density
 783    if (usepaw==1) then
 784      call pawrad_init(core_mesh,mesh_size=pawtab(itypat)%core_mesh_size,&
 785 &     mesh_type=pawrad(itypat)%mesh_type,&
 786 &     rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep)
 787    end if
 788 
 789 !  Loop over atoms of the type
 790    do iat=1,nattyp(itypat)
 791      iatm=iatm+1;iatom=atindx1(iatm)
 792 
 793      if(option==2) corgr(:)=zero
 794 
 795 !    Consider each component in turn : compute range
 796      do mu=1,3
 797 !      Convert reduced coord of given atom to [0,1)
 798        tau(mu)=mod(xred(mu,iatom)+one-aint(xred(mu,iatom)),one)
 799 !      Use tau to find nearest grid point along R(mu)
 800 !      (igrid=0 is the origin; shift by 1 to agree with usual index)
 801        igrid(mu)=nint(tau(mu)*real(ngfft(mu),dp))
 802 !      Use range to compute an index range along R(mu)
 803 !      (add 1 to make sure it covers full range)
 804        irange(mu)=1+nint((range/scale(mu))*real(ngfft(mu),dp))
 805      end do
 806 
 807 !    Allocate arrays that depends on the range
 808      mrange=maxval(irange(1:3))
 809      ABI_ALLOCATE(ii,(2*mrange+1,3))
 810      ABI_ALLOCATE(rrdiff,(2*mrange+1,3))
 811 
 812 !    Set up counters that explore the relevant range of points around the atom
 813      if (geocode=='P') then
 814 !      Fully periodic version
 815        do mu=1,3
 816          ishift=0
 817          do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
 818            ishift=ishift+1
 819            ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
 820            rrdiff(ishift,mu)=real(ixp,dp)/real(ngfft(mu),dp)-tau(mu)
 821          end do
 822          ishiftmax(mu)=ishift
 823        end do
 824      else
 825 !      Free or surface conditions
 826        do mu=1,3
 827          ishift=0
 828          do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
 829            call indpos_mkcore_alt(per(mu),ixp,ngfft(mu),jj,letsgo)
 830            if (letsgo) then
 831              ishift=ishift+1;ii(ishift,mu)=1+jj
 832              rrdiff(ishift,mu)=real(ixp,dp)/real(ngfft(mu),dp)-tau(mu)
 833            end if
 834          end do
 835          ishiftmax(mu)=ishift
 836        end do
 837      end if
 838      npts12=ishiftmax(1)*ishiftmax(2)
 839      ABI_ALLOCATE(indx1,(npts12))
 840      ABI_ALLOCATE(indx2,(npts12))
 841      ABI_ALLOCATE(iindex,(npts12))
 842      ABI_ALLOCATE(rnorm,(npts12))
 843      if (option==1.or.option==3) then
 844        ABI_ALLOCATE(tcore,(npts12))
 845      end if
 846      if (option>=2) then
 847        ABI_ALLOCATE(dtcore,(npts12))
 848      end if
 849      if (option==4) then
 850        ABI_ALLOCATE(d2tcore,(npts12))
 851      end if
 852 
 853 !    Conduct loop over restricted range of grid points for iatom
 854      do ishift3=1,ishiftmax(3)
 855        i3=ii(ishift3,3)
 856        rdiff3=rrdiff(ishift3,3)
 857 
 858        if(fftn3_distrib(i3)/=mpi_enreg%me_fft) cycle
 859 
 860 !      Select the vectors located around the current atom
 861 !        TR: all of the following  could be done inside or
 862 !        outside the loops (i2,i1,i3).
 863 !        Outside: the memory consumption increases.
 864 !        Inside: the time of calculation increases.
 865 !        Here, I choose to do it here, somewhere in the middle.
 866        npts=0
 867        do ishift2=1,ishiftmax(2)
 868          i2=ii(ishift2,2) ; rdiff2=rrdiff(ishift2,2)
 869          rnorm2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2 &
 870 &         +2.0d0*rmet(3,2)*rdiff3*rdiff2
 871          rnorm2_fact=2.0d0*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
 872          do ishift1=1,ishiftmax(1)
 873            i1=ii(ishift1,1) ; rdiff1=rrdiff(ishift1,1)
 874            rnorm2=rnorm2_part+rdiff1*(rnorm2_fact+rmet(1,1)*rdiff1)
 875 !          Only accept contributions inside defined range
 876            if (rnorm2<range2-tol12) then
 877              npts=npts+1 ; iindex(npts)=npts
 878              indx1(npts)=ishift1;indx2(npts)=ishift2
 879              rnorm(npts)=sqrt(rnorm2)
 880            end if
 881          end do
 882        end do
 883        if (npts==0) cycle
 884        if (npts>npts12) then
 885          message='npts>npts12!'
 886          MSG_BUG(message)
 887        end if
 888 
 889 !      Evaluate core density (and derivatives) on the set of selected points
 890        if (usepaw==1) then
 891 !        PAW: use splint routine
 892          call sort_dp(npts,rnorm(1:npts),iindex(1:npts),tol16)
 893          if (option==1.or.option==3) then
 894 !          Evaluate fit of core density
 895            call paw_splint(core_mesh%mesh_size,core_mesh%rad, &
 896 &           pawtab(itypat)%tcoredens(:,1), &
 897 &           pawtab(itypat)%tcoredens(:,3),&
 898 &           npts,rnorm(1:npts),tcore(1:npts))
 899          end if
 900          if (option>=2) then
 901 !          Evaluate fit of 1-der of core density
 902            call paw_splint(core_mesh%mesh_size,core_mesh%rad, &
 903 &           pawtab(itypat)%tcoredens(:,2), &
 904 &           pawtab(itypat)%tcoredens(:,4),&
 905 &           npts,rnorm(1:npts),dtcore(1:npts))
 906          end if
 907          if (option==4) then
 908 !          Evaluate fit of 2nd-der of core density
 909            call paw_splint(core_mesh%mesh_size,core_mesh%rad, &
 910 &           pawtab(itypat)%tcoredens(:,3), &
 911 &           pawtab(itypat)%tcoredens(:,5),&
 912 &           npts,rnorm(1:npts),d2tcore(1:npts))
 913          end if
 914        else
 915 !        Norm-conserving PP:
 916 !          Evaluate spline fit with method from Numerical Recipes
 917 !          (p. 86 Numerical Recipes, Press et al;
 918 !          NOTE error in book for sign of "aa" term in derivative)
 919          do ipts=1,npts
 920            yy=rnorm(ipts)*rangem1
 921            jj=1+int(yy*(n1xccc-1))
 922            diff=yy-(jj-1)*delta
 923            bb = diff*deltam1 ; aa = one-bb
 924            cc = aa*(aa**2-one)*delta2div6
 925            dd = bb*(bb**2-one)*delta2div6
 926            if (option==1.or.option==3) then
 927              tcore(ipts)=aa*xccc1d(jj,1,itypat)+bb*xccc1d(jj+1,1,itypat) +&
 928 &             cc*xccc1d(jj,3,itypat)+dd*xccc1d(jj+1,3,itypat)
 929            end if
 930            if (option>=2) then
 931              dtcore(ipts)=aa*xccc1d(jj,2,itypat)+bb*xccc1d(jj+1,2,itypat) +&
 932 &             cc*xccc1d(jj,4,itypat)+dd*xccc1d(jj+1,4,itypat)
 933            end if
 934            if (option==4) then
 935              d2tcore(ipts)=aa*xccc1d(jj,3,itypat)+bb*xccc1d(jj+1,3,itypat) +&
 936 &             cc*xccc1d(jj,5,itypat)+dd*xccc1d(jj+1,5,itypat)
 937            end if
 938          end do
 939        end if
 940 
 941 !      Now, perform the loop over selected grid points
 942        do ipts=1,npts
 943          ishift1=indx1(iindex(ipts))
 944          ishift2=indx2(iindex(ipts))
 945          difmag=rnorm(ipts)
 946 
 947          rdiff1=rrdiff(ishift1,1);rdiff2=rrdiff(ishift2,2)
 948          jpts=ii(ishift1,1)+n1*(ii(ishift2,2)-1+n2*(ffti3_local(i3)-1))
 949 
 950 !        === Evaluate charge density
 951          if (option==1) then
 952            xccc3d(jpts)=xccc3d(jpts)+tcore(ipts)
 953 
 954 !        === Accumulate contributions to forces
 955          else if (option==2) then
 956            if (difmag>tol10) then
 957              term=vxc_eff(jpts)*dtcore(ipts)/difmag
 958              corgr(1)=corgr(1)+rdiff1*term
 959              corgr(2)=corgr(2)+rdiff2*term
 960              corgr(3)=corgr(3)+rdiff3*term
 961            end if
 962 
 963 !        === Accumulate contributions to stress tensor (in red. coordinates)
 964          else if (option==3) then
 965            if (difmag>tol10) then
 966              term=vxc_eff(jpts)*dtcore(ipts)*rangem1/difmag/real(nfftot,dp)
 967 !            Write out the 6 symmetric components
 968              corfra(1,1)=corfra(1,1)+term*rdiff1*rdiff1
 969              corfra(2,2)=corfra(2,2)+term*rdiff2*rdiff2
 970              corfra(3,3)=corfra(3,3)+term*rdiff3*rdiff3
 971              corfra(3,2)=corfra(3,2)+term*rdiff3*rdiff2
 972              corfra(3,1)=corfra(3,1)+term*rdiff3*rdiff1
 973              corfra(2,1)=corfra(2,1)+term*rdiff2*rdiff1
 974 !            (the above still needs to be transformed to cartesian coords)
 975            end if
 976 !          Also compute a diagonal term
 977            strdia=strdia+vxc_eff(jpts)*tcore(ipts)
 978 
 979 !        === Compute frozen-wf contribution to Dynamical matrix
 980          else if (option==4) then
 981            tt(1)=rmet(1,1)*rdiff1+rmet(1,2)*rdiff2+rmet(1,3)*rdiff3
 982            tt(2)=rmet(2,1)*rdiff1+rmet(2,2)*rdiff2+rmet(2,3)*rdiff3
 983            tt(3)=rmet(3,1)*rdiff1+rmet(3,2)*rdiff2+rmet(3,3)*rdiff3
 984            if (difmag>tol10) then
 985              term=(ucvol/real(nfftot,dp))*vxc_eff(jpts)*rangem1/difmag
 986              term1=term*tcore(ipts)
 987              term2=term*d2tcore(ipts)*rangem1/difmag
 988              do mu=1,3
 989                do nu=1,3
 990                  dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)&
 991 &                 +(term2-term1/difmag**2)*tt(mu)*tt(nu)&
 992 &                 +term1*rmet(mu,nu)
 993                end do
 994              end do
 995            else
 996 !            There is a contribution from difmag=zero !
 997              term=(ucvol/real(nfftot,dp))*vxc_eff(jpts)*d2tcore(ipts)*rangem1**2
 998              do mu=1,3
 999                do nu=1,3
1000                  dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu)
1001                end do
1002              end do
1003            end if
1004          end if ! Choice of option
1005 
1006        end do ! Loop on ipts (ishift1, ishift2)
1007 
1008      end do ! Loop on ishift3
1009 
1010      ABI_DEALLOCATE(ii)
1011      ABI_DEALLOCATE(rrdiff)
1012      ABI_DEALLOCATE(indx1)
1013      ABI_DEALLOCATE(indx2)
1014      ABI_DEALLOCATE(iindex)
1015      ABI_DEALLOCATE(rnorm)
1016      if (allocated(tcore)) then
1017        ABI_DEALLOCATE(tcore)
1018      end if
1019      if (allocated(dtcore)) then
1020        ABI_DEALLOCATE(dtcore)
1021      end if
1022      if (allocated(d2tcore)) then
1023        ABI_DEALLOCATE(d2tcore)
1024      end if
1025 
1026      if (option==2) then
1027        fact=-(ucvol/real(nfftot,dp))/range
1028        grxc(:,iatom)=corgr(:)*fact
1029      end if
1030 
1031 !  End loop on atoms
1032    end do
1033 
1034    if (usepaw==1) then
1035      call pawrad_free(core_mesh)
1036    end if
1037 
1038 !End loop over atom types
1039  end do
1040 
1041  if(option>=2.and.nspden>=2)  then
1042    ABI_DEALLOCATE(vxc_eff)
1043  end if
1044 
1045 !Forces: translate into reduced coordinates
1046  if (option==2) then
1047    do iatom=1,natom
1048      t1=grxc(1,iatom);t2=grxc(2,iatom);t3=grxc(3,iatom)
1049      grxc(:,iatom)=rmet(:,1)*t1+rmet(:,2)*t2+rmet(:,3)*t3
1050    end do
1051  end if
1052 
1053 !Stress tensor: symmetrize, translate into cartesian coord., add diagonal part
1054  if (option==3) then
1055    corstr(1)=corfra(1,1) ; corstr(2)=corfra(2,2)
1056    corstr(3)=corfra(3,3) ; corstr(4)=corfra(3,2)
1057    corstr(5)=corfra(3,1) ; corstr(6)=corfra(2,1)
1058    call strconv(corstr,rprimd,corstr)
1059    corstr(1)=corstr(1)+strdia/real(nfftot,dp)
1060    corstr(2)=corstr(2)+strdia/real(nfftot,dp)
1061    corstr(3)=corstr(3)+strdia/real(nfftot,dp)
1062  end if
1063 
1064 !If needed sum over MPI processes
1065  if(mpi_enreg%nproc_fft>1) then
1066    call timab(539,1,tsec)
1067    if (option==2) then
1068      call xmpi_sum(grxc,mpi_enreg%comm_fft,ier)
1069    end if
1070    if (option==3) then
1071      call xmpi_sum(corstr,mpi_enreg%comm_fft,ier)
1072    end if
1073    if (option==4) then
1074      call xmpi_sum(dyfrx2,mpi_enreg%comm_fft,ier)
1075    end if
1076    call timab(539,2,tsec)
1077  end if
1078 
1079  call timab(12,2,tsec)
1080 
1081  contains