qmcplusplus::ewaldref Namespace Reference

Classes
class	KspaceEwaldTerm
	Functor for term within the k-space sum in Drummond 2008 formula 6. More...
class	KspaceMadelungTerm
	Functor for term within the k-space sum in Drummond 2008 formula 7. More...
class	RspaceEwaldTerm
	Functor for term within the real-space sum in Drummond 2008 formula 6. More...
class	RspaceMadelungTerm
	Functor for term within the real-space sum in Drummond 2008 formula 7. More...

Typedefs
using	int_t = int
	Type for integers. More...
using	real_t = double
	Type for floating point numbers. More...
using	IntVec = TinyVector< int_t, DIM >
	Type for integer vectors of length DIM. More...
using	RealVec = TinyVector< real_t, DIM >
	Type for floating point vectors of length DIM. More...
using	RealMat = Tensor< real_t, DIM >
	Type for floating point matrices of shape DIM,DIM. More...
using	PosArray = std::vector< RealVec >
	Type for lists of particle positions. More...
using	ChargeArray = std::vector< real_t >
	Type for lists of particle charges. More...

Enumerations
enum	{ DIM = 3 }
	Reference Ewald implemented for 3D only. More...

Functions
template<typename T >
real_t	gridSum (T &function, bool zero=true, real_t tol=1e-11)
	Perform a sum over successively larger cubic integer grids in DIM dimensional space for arbitrary functors. More...
real_t	getKappaMadelung (real_t volume)
	Find the optimal kappa for Madelung sums. More...
real_t	madelungSum (const RealMat &a, real_t tol=1e-10)
	Compute the Madelung constant to a given tolerance. More...
real_t	getKappaEwald (real_t volume)
	Find the optimal kappa for Ewald pair sums. More...
real_t	ewaldSum (const RealVec &r, const RealMat &a, real_t tol=1e-10)
	Compute the Ewald interaction of a particle pair to a given tolerance. More...
real_t	ewaldEnergy (const RealMat &a, const PosArray &R, const ChargeArray &Q, real_t tol)
	Compute the total Ewald potential energy for a collection of charges. More...

Typedef Documentation

ChargeArray

using ChargeArray = std::vector<real_t>

Type for lists of particle charges.

Definition at line 54 of file EwaldRef.h.

int_t

using int_t = int

Type for integers.

Definition at line 42 of file EwaldRef.h.

IntVec

using IntVec = TinyVector<int_t, DIM>

Type for integer vectors of length DIM.

Definition at line 46 of file EwaldRef.h.

PosArray

using PosArray = std::vector<RealVec>

Type for lists of particle positions.

Definition at line 52 of file EwaldRef.h.

real_t

using real_t = double

Type for floating point numbers.

Definition at line 44 of file EwaldRef.h.

RealMat

using RealMat = Tensor<real_t, DIM>

Type for floating point matrices of shape DIM,DIM.

Definition at line 50 of file EwaldRef.h.

RealVec

using RealVec = TinyVector<real_t, DIM>

Type for floating point vectors of length DIM.

Definition at line 48 of file EwaldRef.h.

Enumeration Type Documentation

anonymous enum

anonymous enum

Reference Ewald implemented for 3D only.

Enumerator
DIM

Definition at line 36 of file EwaldRef.h.

{
  DIM = 3
};

Function Documentation

ewaldEnergy()

real_t ewaldEnergy	(	const RealMat &	a,
		const PosArray &	R,
		const ChargeArray &	Q,
		real_t	tol
	)

Compute the total Ewald potential energy for a collection of charges.

Corresponds to the entirety of Drummond 2008 formula 5, but for arbitrary charges.

Parameters

a	Real-space cell axes.
R	List of particle coordinates.
R	List of particle charges.
tol	Tolerance for the total potential energy in Ha.

Definition at line 337 of file EwaldRef.cpp.

References qmcplusplus::abs(), qmcplusplus::createGlobalTimer(), DIM, qmcplusplus::dot(), ewaldSum(), qmcplusplus::floor(), qmcplusplus::inverse(), madelungSum(), and qmcplusplus::Units::force::N.

Referenced by CoulombPBCAA::CoulombPBCAA(), and qmcplusplus::TEST_CASE().

{
  // Timer for EwaldRef
  ScopedTimer totalEwaldTimer(createGlobalTimer("EwaldRef"));

  // Number of particles
  const size_t N = R.size();

  // Total Ewald potential energy
  real_t ve = 0.0;

  {
    // Sum Madelung contributions
    ScopedTimer totalMadelungTimer(createGlobalTimer("MadelungSum"));
    // Maximum self-interaction charge product
    real_t qqmax = 0.0;
    for (size_t i = 0; i < N; ++i)
      qqmax = std::max(std::abs(Q[i] * Q[i]), qqmax);

    // Compute the Madelung term (Drummond 2008 formula 7)
    real_t vm = madelungSum(a, tol * 2. / qqmax);

    // Sum the Madelung self interaction for each particle
    for (size_t i = 0; i < N; ++i)
      ve += Q[i] * Q[i] * vm / 2;
  }

  {
    // Sum the interaction terms for all particle pairs
    ScopedTimer EwaldSumTimer(createGlobalTimer("EwaldSum"));

#pragma omp parallel for reduction(+ : ve)
    for (size_t i = 1; i < N / 2 + 1; ++i)
    {
      for (size_t j = 0; j < i; ++j)
      {
        real_t qq       = Q[i] * Q[j];
        RealVec reduced = dot(R[i] - R[j], inverse(a));
        for (size_t dim = 0; dim < DIM; dim++)
          reduced[dim] -= std::floor(reduced[dim]);
        RealVec rr = dot(reduced, a);
        ve += qq * ewaldSum(rr, a, tol / qq);
      }

      const size_t i_reverse = N - i;
      if (i == i_reverse)
        continue;

      for (size_t j = 0; j < i_reverse; ++j)
      {
        real_t qq       = Q[i_reverse] * Q[j];
        RealVec reduced = dot(R[i_reverse] - R[j], inverse(a));
        for (size_t dim = 0; dim < DIM; dim++)
          reduced[dim] -= std::floor(reduced[dim]);
        RealVec rr = dot(reduced, a);
        ve += qq * ewaldSum(rr, a, tol / qq);
      }
    }
  }

  return ve;
}

ewaldSum()

real_t qmcplusplus::ewaldref::ewaldSum	(	const RealVec &	r,
		const RealMat &	a,
		real_t	tol = 1e-10
	)

Compute the Ewald interaction of a particle pair to a given tolerance.

Corresponds to the entirety of Drummond 2008 formula 6.

Parameters

r	Inter-particle separation vector.
a	Real-space cell axes.
tol	Tolerance for the Ewald pair interaction in Ha.

Definition at line 292 of file EwaldRef.cpp.

References qmcplusplus::abs(), qmcplusplus::det(), getKappaEwald(), gridSum(), qmcplusplus::inverse(), qmcplusplus::pow(), qmcplusplus::sqrt(), and qmcplusplus::transpose().

Referenced by ewaldEnergy().

{
  // Real-space cell volume
  real_t volume = std::abs(det(a));
  // k-space cell axes
  RealMat b = 2 * M_PI * transpose(inverse(a));
  // k-space cutoff (kappa)
  real_t kconv = getKappaEwald(volume);

  // Set constants for real-/k-space Ewald pair functors
  real_t rconst  = 1. / (std::sqrt(2.) * kconv);
  real_t kconst  = -std::pow(kconv, 2) / 2;
  real_t kfactor = 4 * M_PI / volume;

  // Create real-/k-space fuctors for terms within the sums in formula 6
  RspaceEwaldTerm rfunc(r, a, rconst);
  KspaceEwaldTerm kfunc(r, b, kconst, kfactor);

  // Compute the constant term
  real_t cval = -2 * M_PI * std::pow(kconv, 2) / volume;
  // Compute the real-space sum (includes zero)
  real_t rval = gridSum(rfunc, true, tol);
  // Compute the k-space sum (excludes zero)
  real_t kval = gridSum(kfunc, false, tol);

  // Sum all contributions to get the Ewald pair interaction
  real_t es = cval + rval + kval;

  return es;
}

getKappaEwald()

real_t qmcplusplus::ewaldref::getKappaEwald

(

real_t

volume

)

Find the optimal kappa for Ewald pair sums.

The optimal kappa balances the number of points within a given isosurface of the Gaussians (or limiting Gaussians from erfc) in the real-space and k-space Ewald pair terms. The balancing condition is made under isotropic assumptions, as reflected by the use of a sphere equal in volume to the simulation cell to determine the radius.

Parameters

volume

Volume of the real space cell.

Definition at line 275 of file EwaldRef.cpp.

References qmcplusplus::pow(), and qmcplusplus::sqrt().

Referenced by ewaldSum().

{
  real_t radius = std::pow(3. * volume / (4 * M_PI), 1. / 3);
  return radius / std::sqrt(2 * M_PI);
}

getKappaMadelung()

real_t qmcplusplus::ewaldref::getKappaMadelung

(

real_t

volume

)

Find the optimal kappa for Madelung sums.

The optimal kappa balances the number of points within a given isosurface of the Gaussians (or limiting Gaussians from erfc) in the real-space and k-space Madelung terms. The balancing condition is made under isotropic assumptions, as reflected by the use of a sphere equal in volume to the simulation cell to determine the radius.

Parameters

volume

Volume of the real space cell.

Definition at line 217 of file EwaldRef.cpp.

References qmcplusplus::pow(), and qmcplusplus::sqrt().

Referenced by madelungSum().

{
  real_t radius = std::pow(3. * volume / (4 * M_PI), 1. / 3);
  return std::sqrt(M_PI) / radius;
}

gridSum()

real_t qmcplusplus::ewaldref::gridSum	(	T &	function,
		bool	zero = true,
		real_t	tol = 1e-11
	)

Perform a sum over successively larger cubic integer grids in DIM dimensional space for arbitrary functors.

Parameters

function	A functor accepting a point in the grid and returning the real-valued contribution to the sum from that point.
zero	Include the origin in the sum (or not).
tol	Tolerance for the sum. Summation ceases when the contribution to the sum from the outermost cubic shell of points is less than tol.

Definition at line 140 of file EwaldRef.cpp.

References qmcplusplus::abs().

Referenced by ewaldSum(), and madelungSum().

{
  real_t dv  = std::numeric_limits<real_t>::max();
  real_t dva = std::numeric_limits<real_t>::max();
  real_t v   = 0.0;
  int_t im   = 0;
  int_t jm   = 0;
  int_t km   = 0;
  IntVec iv;
  // Add the value at the origin, if requested.
  if (zero)
  {
    iv = 0;
    v += function(iv);
  }
  // Sum over cubic surface shells until the tolerance is reached.
  while (std::abs(dva) > tol)
  {
    dva = 0.0;
    dv  = 0.0; // Surface shell contribution.
    // Sum over new surface planes perpendicular to the x direction.
    im += 1;
    for (int_t i : {-im, im})
      for (int_t j = -jm; j < jm + 1; ++j)
        for (int_t k = -km; k < km + 1; ++k)
        {
          iv[0]    = i;
          iv[1]    = j;
          iv[2]    = k;
          real_t f = function(iv);
          dv += f;
          dva += std::abs(f);
        }
    // Sum over new surface planes perpendicular to the y direction.
    jm += 1;
    for (int_t j : {-jm, jm})
      for (int_t k = -km; k < km + 1; ++k)
        for (int_t i = -im; i < im + 1; ++i)
        {
          iv[0]    = i;
          iv[1]    = j;
          iv[2]    = k;
          real_t f = function(iv);
          dv += f;
          dva += std::abs(f);
        }
    // Sum over new surface planes perpendicular to the z direction.
    km += 1;
    for (int_t k : {-km, km})
      for (int_t i = -im; i < im + 1; ++i)
        for (int_t j = -jm; j < jm + 1; ++j)
        {
          iv[0]    = i;
          iv[1]    = j;
          iv[2]    = k;
          real_t f = function(iv);
          dv += f;
          dva += std::abs(f);
        }
    v += dv;
  }

  return v;
}

madelungSum()

real_t qmcplusplus::ewaldref::madelungSum	(	const RealMat &	a,
		real_t	tol = 1e-10
	)

Compute the Madelung constant to a given tolerance.

Corresponds to the entirety of Drummond 2008 formula 7.

Parameters

a	Real-space cell axes.
tol	Tolerance for the Madelung constant in Ha.

Definition at line 232 of file EwaldRef.cpp.

References qmcplusplus::abs(), qmcplusplus::det(), getKappaMadelung(), gridSum(), qmcplusplus::inverse(), qmcplusplus::Units::time::ms, qmcplusplus::pow(), qmcplusplus::sqrt(), and qmcplusplus::transpose().

Referenced by ewaldEnergy().

{
  // Real-space cell volume
  real_t volume = std::abs(det(a));
  // k-space cell axes
  RealMat b = 2 * M_PI * transpose(inverse(a));
  // k-space cutoff (kappa)
  real_t kconv = getKappaMadelung(volume);

  // Set constants for real-/k-space Madelung functors
  real_t rconst  = kconv;
  real_t kconst  = -1. / (4 * std::pow(kconv, 2));
  real_t kfactor = 4 * M_PI / volume;

  // Create real-/k-space fuctors for terms within the sums in formula 7
  RspaceMadelungTerm rfunc(a, rconst);
  KspaceMadelungTerm kfunc(b, kconst, kfactor);

  // Compute the constant term
  real_t cval = -M_PI / (std::pow(kconv, 2) * volume) - 2 * kconv / std::sqrt(M_PI);
  // Compute the real-space sum (excludes zero)
  real_t rval = gridSum(rfunc, false, tol);
  // Compute the k-space sum (excludes zero)
  real_t kval = gridSum(kfunc, false, tol);

  // Sum all contributions to get the Madelung constant
  real_t ms = cval + rval + kval;

  return ms;
}