forcefields::ZhuPhilpott< P > Class Template Reference

Forcefield for water and platinum interactions. More...

#include <zhu_philpott.hpp>

Inheritance diagram for forcefields::ZhuPhilpott< P >:

Public Member Functions
	ZhuPhilpott ()
	ZhuPhilpott (double cutoff, double switchingWidth)
	ZhuPhilpott (ZhuPhilpott const &)
void	operator= (ZhuPhilpott const &)
	~ZhuPhilpott ()
void	computeHH_O_Pt_ (const int nWater, const int nPt, const double r[], double f[], double &energy, double const b[], bool const fixed[])
	Compute water-platinum forcefield.
void	computeHH_O_ (const int nWater, const double r[], double f[], double &energy, double const b[], bool const fixed[])
	Compute water-platinum interactions, call with water's positions only.
int	nPlatinum () const
	Number of platinum atoms.
void	setPlatinum (int nPlatinum, double const positions[])
	Initialises the positions of the atoms of platinum.
template<int H, int O, int H3, int O3>
void	computeTemplate (const int nWater, const double(*const rh1)[H3], const double(*const rh2)[H3], const double(*const ro)[O3], double(*const fh1)[H3], double(*const fh2)[H3], double(*const fo)[O3], const int nPt, const double rPt[][3], double fPt[][3], double &energy, double const b[], bool const (*const xh1)[H], bool const (*const xh2)[H], bool const (*const xo)[O], bool const xPt[])
Public Member Functions inherited from forcefields::SpceCcl
	SpceCcl ()
	SpceCcl (double cutoff, double switchingWidth)
void	computeHH_O_ (const int nAtoms, const double R[], double F[], double &U, const double b[])
	Compute the forces and the energy.
void	computeHH_O_ (const int nAtoms, const double R[], double F[], double &U, const double b[], const bool fixed[])
	Compute the forces and the energy (with fixed atom optimisation).
char const *	getName () const
	Name of the potential.
Public Member Functions inherited from forcefields::Ccl
	Ccl ()
void	computeHH_O_ (const int nAtoms, const double R[], double F[], double &U, const double b[])
	Compute the forces and the energy.
void	computeHH_O_ (const int nAtoms, const double R[], double F[], double &U, const double b[], const bool fixed[])
	Compute the forces and the energy (with fixed atom optimisation).
char const *	getName () const
	Name of the potential.
Public Member Functions inherited from forcefields::PotentialBase
	PotentialBase ()
	Non bond interaction cutoff.
	PotentialBase (double cutoff, double switchingWidth)
virtual	~PotentialBase ()
double	getCutoff () const
void	setCutoff (double cutoff)
double	getSwitchingWidth () const
	Width of the switching zone.
void	setSwitchingWidth (double width)

Static Public Member Functions
static char const *	getName ()
	Name of the potential.
static double	applyPeriodicity0 (double r, double const period)
	Minimum image representation.
static void	applyPeriodicity0 (double r[], double const periods[])
	Minimum image representation.

Private Member Functions
template<int H, int O, int H3, int O3>
void	computeTemplate (const int nWater, const double(*const rh1)[H3], const double(*const rh2)[H3], const double(*const ro)[O3], double(*const fh1)[H3], double(*const fh2)[H3], double(*const fo)[O3], const int nPt, const double rPt[][3], double fPt[][3], double &energy, double const b[], bool const (*const xh1)[H]=0, bool const (*const xh2)[H]=0, bool const (*const xo)[O]=0, bool const *xPt=0)
void	interactWithCorePt (Water &water, int const nPt, double const rPt[][3], double fPt[][3], double &energy)
	Interaction of one molecule of water with the whole platinum.
void	interactionPtO (double const R1[], double const R2[], double F1[], double F2[], double &energy)
void	interactionPtH (double const R1[], double const R2[], double F1[], double F2[], double &energy)
void	anisotropic (const double distance[], double force[], double &energy, double const epsilon, double const sigma, double const alpha)
	Anisotropic interaction between water and platinum.
void	isotropic10 (double const distance, double &force, double &energy, double const epsilon, double const sigma, double const C10)
	Isotropic interaction between water and platinum.
void	interactWithImage (Water &w1, Water &w2, double &U)
	Interactions of a molecule with an image.
void	coulombWithCutoff (Water &w1, Water &w2, double &u, double const relativePermittivity)
	Coulomb interaction between two molecules within the swithcing zone.
void	coulombFull (Water &w1, Water &w2, double &U, double const relativePermittivity)
	Coulomb interaction between two molecules with cutoff.

Private Attributes
int	nPlatinum_ {0}
std::vector< double >	positions_
std::vector< double >	forces_

Additional Inherited Members
Protected Member Functions inherited from forcefields::SpceCcl
void	intramolecular (Water &water, double &U)
	Interactions within a molecules.
void	lennardJonesWithCutoff (Water &w1, Water &w2, double &U)
	Interactions between two molecules.
void	coulombWithCutoff (Water &w1, Water &w2, double &U)
	Interactions between two molecules.
void	coulombFull (Water &w1, Water &w2, double &U)
	Interactions between two molecules.
Protected Member Functions inherited from forcefields::Ccl
	Ccl (double cutoff, double switchingWidth)
	Constructor with cutoff for derived classes.
void	intramolecular (double const rh1[], double const rh2[], double const ro[], double fh1[], double fh2[], double fo[], double &energy)
	Interactions inside one molecules.
void	initialiseRho (Vector3 const &v, Rho &r)
	Initialise Rho.
void	initialiseDtheta (Vector3 const &v1, Vector3 const &v2, Dtheta &dth, double const thetaEquilibrium)
	Initialise Dtheta.
void	intramolecular (Rho const &ro1, Rho const &ro2, Dtheta const &dth, double &energy, double fh1[], double fh2[], double fo[])
	Compute intramolecular energy and forces.
Protected Member Functions inherited from forcefields::PotentialBase
void	restrainAngle (double const r1[], double const r2[], double const r3[], double f1[], double f2[], double f3[], double &u, const double k, const double aeq)
	Angular quadratic restraints.
void	calculateCentre (double const r1[], double const r2[], double rc[])
	Calculate centre of two points.
void	calculateCentre (double const r1[], double const r2[], double const r3[], double rc[])
	Calculate centre of three points.
void	calculateWeightedCentre (double const w1, double const w2, double const w3, double const r1[], double const r2[], double const r3[], double rc[])
	Calculate barycentre.
void	coulombWithCutoff (const double r1[], const double r2[], double f1[], double f2[], double &u, double const qq)
	Coulomb interaction with single charge based cutoff.
void	coulomb (const double r1[], const double r2[], double f1[], double f2[], double &u, double const qq)
	Compute Coulomb interaction between two charges.
void	coulomb (double distance, double &force, double &energy, double const qq)
	Compute Coulomb interaction between two charges.
void	spreadForce (double f1[], double f2[], double const fc[])
	Spread force on centre to other points.
void	spreadForce (double f1[], double f2[], double f3[], double const fc[])
	Spread force on centre to other points.
void	spreadWeightedForce (double const w1, double const w2, double const w3, double f1[], double f2[], double f3[], double const fc[])
	Spread force on barycentre to atoms.
template<int N, class R, class F>
void	addForces (ChargeGroup< N, R, F > const &g1, ChargeGroup< N, R, F > &g2)
	Increment forces.
template<int N, class R, class F>
void	switching (ChargeGroup< N, R, F > &g1, ChargeGroup< N, R, F > &g2, double &energy, double cutoff, double switchingWidth)
void	computePt (int const nAtoms, double positions[], double forces[], double &energy, double const periods[], bool const fixed[])
	Potential for Platinum.
void	restrainLength (const double R1[], const double R2[], double F1[], double F2[], double &u, const double k, const double r0)
	Compute quadratic restraints between two atoms.
void	lennardJones (double const distance, double &force, double &energy, double const epsilon, double const sigma)
	Lennard-Jones 12-6 between two atoms.
void	lennardJones (const double R1[], const double R2[], double F1[], double F2[], double &E, double const epsilon, double const sigma)
	Lennard-Jones 12-6 between two atoms.
void	lennardJonesWithCutoff (double const r1[], double const r2[], double f1[], double f2[], double &energy, double const epsilon, double const sigma)
	Lennard Jones 12-6 Potential with cutoff.
void	switching (double const distance, double &force, double &energy)
	Smooth cutoff switch off.
void	switching (double const r1[], double const r2[], double f1[], double f2[], double &energy)
	Smooth cutoff switch off.
double	applyPeriodicity1 (double r, int const axis)
	Minimum image representation.
void	applyPeriodicity1 (double r[])
	Minimum image representation.
void	distance (const double x[], const double y[], double z[], double &z1, double &z2)
	Distance vector, norm and norm square.
void	distance (const double x[], const double y[], double z[], double &z1)
	Distance.
void	distance (const double x[], const double y[], double z[])
	Distance vector, norm.
void	distance (const double x[], const double y[], double &z1)
	Distance.
void	distance (const double x[], const double y[], Vector3 &z)
	Distance vector, norm and norm square.
void	setPeriodicity (const double periods[])
	Set periodicity.
double	unBreak1 (double const r, double const ref, int const axis)
	Undo the separation of two atoms created by the periodic boundaries.
void	unBreak1 (double r[], double const ref[])
	Undo the separation of two atoms created by the periodic boundaries.
Static Protected Member Functions inherited from forcefields::PotentialBase
static void	divide (double v[], double const divisor)
	Divide vector.
static double *	crossProduct (const double v[], const double w[], double e[])
	Cross product.
static double	dotProduct (double const v[], double const w[])
	Dot product.
static void	multiply (double v[], double const factor)
	Multiply vector.
static double	norm (double const v[])
	Norm of 3D vector.
static void	normalise (double v[])
	Normalise 3D vector \( \frac{\mathbf v}{ |\mathbf v|} \).
static double	isoscelesBase (double length, double angle)
	Calculate the base of an isosceles triangle.
static double	epsilon (double const A, double const B)
	Conversion for Lennard-Jones.
static double	sigma (double const A, double const B)
	Conversion for Lennard-Jones.
static double	smithKongEpsilon (double sigma1, double epsilon1, double sigma2, double epsilon2)
	Smith and Kong combination rules.
static double	smithKongSigma (double sigma1, double epsilon1, double sigma2, double epsilon2)
	Smith and Kong combination rules.
static double	unBreak0 (double const r, double const ref, double const period)
	Undo the separation of two atoms created by the periodic boundaries.
static void	unBreak0 (double r[], double const ref[], double const periods[])
	Undo the separation of two atoms created by the periodic boundaries.
Protected Attributes inherited from forcefields::PotentialBase
double	cutoff_
double	periods_ [3]
double	switchingWidth_
Static Protected Attributes inherited from forcefields::SpceCcl
static const double	roh_ = 1.0
	Distance OH.
static const double	theta_ = 1.91063
	Angle HOH.
static const double	rhh_ = 1.63299
	Distance HH.
static const double	charge_ = 0.4238
	Charge on one hydrogen.
static const double	charge2_ = 0.179606
	Square of # charge_.
static const double	A_ = 27291.6
	Lennard-Jones.
static const double	B_ = 27.1223
	Lennard-Jones.
static const double	sigma_ = 3.16556
	Lennard-Jones. See PotentialBase::lennardJones() for definition.
static const double	epsilon_ = 0.00673853
	Lennard-Jones. See PotentialBase::lennardJones() for definition.
static double const	polarisationEnergy_ = 0.0541015
	Polarisation correction.
Static Protected Attributes inherited from forcefields::Ccl
static double const	re_ = ::re_
	Distance OH at equilibrium.
static double const	thetae_ = ::thetae_
	Distance OH at equilibrium.
static const double	ONE_OVER_4_PI_EPSILON0
static const double	EPSILON_PT = 0.68165797577788501
	Platinum Lennard-Jones.
static const double	SIGMA_PT = 2.54
Static Private Attributes inherited from forcefields::zhu_philpott_parameters::Standard
static double const	sigmaO_ = 2.86
static double const	epsilonO_ = 0.0023734176137013181
static double const	sigmaH_ = 2.56
static double const	epsilonH_ = 0.00087578073518673099
static double const	C10_O_ = 1.28
static double const	C10_H_ = 1.2
static double const	alpha_ = 0.8
static double const	sigmaHPt_ = 2.730249677569295
static double const	sigmaOPt_ = 2.7735458150747108
static double const	epsilonHPt_ = 0.016217645873043762
static double const	epsilonOPt_ = 0.03387330127291549

Detailed Description

template<class P = zhu_philpott_parameters::Standard>
class forcefields::ZhuPhilpott< P >

Forcefield for water and platinum interactions.

This forcefield is the A2 water-platinum potential invented by Zhuand Philpott. It includes the SPC/E+CCL potential for interaction between the molecules of water. The SPC/E is a potential for constrained water. This implementation includes restraints so it may be used without constraint (see class SpceCcl).
The potential uses Kong's rules to combine the Lennard-Jones parameters of platinum with oxygen and hydrogen. A review of combination rules including Kong's rule can be find on the web.
The system of unit used by the class is (eV, Angstrom, fs, e).

References

Interaction of water with metal surfaces, S.-B. Zhu and M.R. Philpott J. Chem. Phys. (1994) vol. 100, No 9, p. 6961.
A Molecular Dynamics Simulation of Water Droplet in Contact with a Platinum Surface, T. Kimura and S. Maruyama, University of Tokyo.
Unlike Lennard-Jones Parameters for Vapor-Liquid Equilibria, Thorsten Schnabel, Jadran Vrabec , Hans Hasse, Institut fur Technische Thermodynamik und Thermische Verfahrenstechnik, Universitat Stuttgart, D-70550 Stuttgart, Germany, http://www.itt.uni-stuttgart.de/~schnabel/CR.pdf.

Definition at line 25 of file zhu_philpott.hpp.

Constructor & Destructor Documentation

ZhuPhilpott() [1/3]

template<class P>

ZhuPhilpott::ZhuPhilpott

(

)

Definition at line 66 of file zhu_philpott.cpp.

    : SpceCcl() {
  nPlatinum_ = 0;
  positions_ = 0;
  forces_ = 0;
}

ZhuPhilpott() [2/3]

template<class P>

ZhuPhilpott::ZhuPhilpott	(	double	cutoff,
		double	switchingWidth )

Definition at line 74 of file zhu_philpott.cpp.

    : SpceCcl(cutoff, switchingWidth) {
  nPlatinum_ = 0;
  positions_ = 0;
  forces_ = 0;
}

ZhuPhilpott() [3/3]

template<class P>

ZhuPhilpott::ZhuPhilpott

(

ZhuPhilpott< P > const &

zhuPhilpott

)

Definition at line 81 of file zhu_philpott.cpp.

                                                                             {
  operator=(zhuPhilpott);
}

~ZhuPhilpott()

template<class P>

ZhuPhilpott::~ZhuPhilpott ( )

default

Member Function Documentation

anisotropic()

template<class P>

void ZhuPhilpott::anisotropic ( const double distance[], double force[], double & energy, double const epsilon, double const sigma, double const alpha )

private

Anisotropic interaction between water and platinum.

The isotropic iteraction between a atom Pt and an atom O or H of a molecule of water has the function form:

\[ E=4\epsilon\left[ \left(\frac{\sigma^2}{\alpha^2\rho^2+z^2}\right)^6-\left(\frac{\sigma^2}{\rho^2/\alpha^2+z^2}\right)^3 \right] \]

Parameters

[in]	distance	Vector distance between to atom.
[in,out]	force	Vector force created by the interaction between the two atoms.
[in,out]	energy	Energy created by the interaction.
[in]	epsilon	See equation.
[in]	sigma	See equation.
[in]	alpha	See equation.

Definition at line 388 of file zhu_philpott.cpp.

                                                                         {
  double z2, z, a, b, A, B, dE_da, dE_db, rho, rho2, alpha2;
  // WARNING: F1 and F2 are incremented.
  alpha2 = alpha * alpha;
  double f;
  rho2 = distance[0] * distance[0] + distance[1] * distance[1];
  rho = std::sqrt(rho2);
  z = distance[2];
  z2 = z * z;
  double sigma2 = sigma * sigma;
  a = sigma2 / (rho2 * alpha2 + z2);
  b = sigma2 / (rho2 / alpha2 + z2);
  double a2 = a * a, a3 = a2 * a;
  double b2 = b * b;
  A = a3 * a3; // a^6
  B = b * b2;  // b^3
  energy += 4.0 * epsilon * (A - B);
  dE_da = 24.0 * epsilon * a2 * a3; // a^5
  dE_db = -12.0 * epsilon * b2;     // b^2
  double invSigma2 = 1.0 / sigma2;
  f = 2.0 * (dE_da * a2 * alpha2 + dE_db * b2 / alpha2) * invSigma2;
  force[0] += f * distance[0];
  force[1] += f * distance[1];
  f = 2.0 * (dE_da * a2 + dE_db * b2) * invSigma2;
  force[2] += f * distance[2];
}

computeHH_O_()

template<class P>

void ZhuPhilpott::computeHH_O_	(	const int	nWater,
		const double	r[],
		double	f[],
		double &	energy,
		double const	b[],
		bool const	fixed[] )

Compute water-platinum interactions, call with water's positions only.

This function is called with the coordinates of the atoms of water only. The positions of the platinum atoms must have been previously provided through function setPlatinum().

Parameters

[in]	nWater	Number of molecules of water.
[in]	r	Positions of atoms in water.
[out]	f	Forces on the atoms in water.
[out]	energy	Potential energy.
[in]	b	Periodic boundaries. Length: 3
[in]	fixed	True for fixed atoms, false otherwise. The length of arrays r and f is \( 9\times nWater \). The length of array fixed is \( 3 \times nWater \).

Note

The function does not return the forces on platinum atoms as they as assumed to be fixed.

See also

setPlatinum() and computeHH_O_Pt_().

Definition at line 161 of file zhu_philpott.cpp.

                                                      {
  const double (*const rh1)[6] = reinterpret_cast<const double (*)[6]>(r);
  const double (*const rh2)[6] = reinterpret_cast<const double (*)[6]>(&r[3]);
  const double (*const ro)[3] =
      reinterpret_cast<const double (*)[3]>(&r[nWater * 6]);
  const double (*const rPt)[3] =
      reinterpret_cast<const double (*)[3]>(positions_.data());
 
  double (*const fh1)[6] = reinterpret_cast<double (*)[6]>(f);
  double (*const fh2)[6] = reinterpret_cast<double (*)[6]>(&f[3]);
  double (*const fo)[3] = reinterpret_cast<double (*)[3]>(&f[nWater * 6]);
  double (*const fPt)[3] = reinterpret_cast<double (*)[3]>(forces_.data());
 
  if (fixed) {
    bool const(*const xh1)[2] = reinterpret_cast<bool const(*)[2]>(fixed);
    bool const(*const xh2)[2] = reinterpret_cast<bool const(*)[2]>(&fixed[1]);
    bool const(*const xo)[1] =
        reinterpret_cast<bool const(*)[1]>(&fixed[nWater * 2]);
    computeTemplate(nWater, rh1, rh2, ro, fh1, fh2, fo, nPlatinum_, rPt, fPt,
                    energy, b, xh1, xh2, xo, 0);
  } else {
    bool const(*const xh1)[2] = 0;
    bool const(*const xh2)[2] = 0;
    bool const(*const xo)[1] = 0;
    computeTemplate(nWater, rh1, rh2, ro, fh1, fh2, fo, nPlatinum_, rPt, fPt,
                    energy, b, xh1, xh2, xo, 0);
  };
}

computeHH_O_Pt_()

template<class P>

void ZhuPhilpott::computeHH_O_Pt_	(	const int	nWater,
		const int	nPt,
		const double	r[],
		double	f[],
		double &	energy,
		double const	b[],
		bool const	fixed[] )

Compute water-platinum forcefield.

Parameters

[in]	nWater	Number of molecules of water.
[in]	nPt	Number of platinum atoms.
[in]	r	Positions of the atoms.
[out]	f	Forces on the atoms.
[out]	energy	Potential energy.
[in]	b	Periodic boundaries. Length: 3
[in]	fixed	True for fixed atoms, false otherwise. The order of the atoms is all hydrogens (hydrogens belonging to the same molecule are next to each others), all oxygen, all platinum. The length of arrays r and f is \( 3\times(3 \times nWater+nPt) \). The length of array fixed is \( 3 \times nWater+nPt \). When fixed is provided the interaction between two fixed atoms may be skipped to speed up the calculation. The parameter is optional.

Warning

The conductivity is simulated by interaction a mirror images of charges. The mirror surface is set at z=0, which should correspond to the location of the highest layer of platinum atoms.

Definition at line 111 of file zhu_philpott.cpp.

                                                         {
  const double (*const rh1)[6] = reinterpret_cast<const double (*)[6]>(r);
  const double (*const rh2)[6] = reinterpret_cast<const double (*)[6]>(&r[3]);
  const double (*const ro)[3] =
      reinterpret_cast<const double (*)[3]>(&r[nWater * 6]);
  const double (*const rPt)[3] =
      reinterpret_cast<const double (*)[3]>(&r[nWater * 9]);
 
  double (*const fh1)[6] = reinterpret_cast<double (*)[6]>(f);
  double (*const fh2)[6] = reinterpret_cast<double (*)[6]>(&f[3]);
  double (*const fo)[3] = reinterpret_cast<double (*)[3]>(&f[nWater * 6]);
  double (*const fPt)[3] = reinterpret_cast<double (*)[3]>(&f[nWater * 9]);
 
  if (fixed) {
    bool const(*const xh1)[2] = reinterpret_cast<bool const(*)[2]>(fixed);
    bool const(*const xh2)[2] = reinterpret_cast<bool const(*)[2]>(&fixed[1]);
    bool const(*const xo)[1] =
        reinterpret_cast<bool const(*)[1]>(&fixed[nWater * 2]);
    bool const *const xPt = &fixed[nWater * 3];
    computeTemplate(nWater, rh1, rh2, ro, fh1, fh2, fo, nPt, rPt, fPt, energy,
                    b, xh1, xh2, xo, xPt);
  } else {
    bool const(*const xh1)[2] = 0;
    bool const(*const xh2)[2] = 0;
    bool const(*const xo)[1] = 0;
    computeTemplate(nWater, rh1, rh2, ro, fh1, fh2, fo, nPt, rPt, fPt, energy,
                    b, xh1, xh2, xo, 0);
  };
}

computeTemplate() [1/2]

template<class P = zhu_philpott_parameters::Standard>

template<int H, int O, int H3, int O3>

void forcefields::ZhuPhilpott< P >::computeTemplate	(	const int	nWater,
		const double(*)	rh1[H3],
		const double(*)	rh2[H3],
		const double(*)	ro[O3],
		double(*)	fh1[H3],
		double(*)	fh2[H3],
		double(*)	fo[O3],
		const int	nPt,
		const double	rPt[][3],
		double	fPt[][3],
		double &	energy,
		double const	b[],
		bool const (*)	xh1[H],
		bool const (*)	xh2[H],
		bool const (*)	xo[O],
		bool const	xPt[] )

Definition at line 215 of file zhu_philpott.cpp.

                                                 {
  for (int i = 0; i < nWater; ++i) {
    for (int a = 0; a < 3; ++a) {
      fh1[i][a] = 0.0;
      fh2[i][a] = 0.0;
      fo[i][a] = 0.0;
    };
  };
  for (int i = 0; i < nPt; ++i) {
    for (int a = 0; a < 3; ++a)
      fPt[i][a] = 0.0;
  };
  energy = 0.0;
  setPeriodicity(b);
  for (int i = nWater - 1; i >= 0; --i) {
    double centre1[3] = {0};
    calculateCentre(rh1[i], rh2[i], ro[i], centre1);
    Water w1(rh1[i], rh2[i], ro[i], centre1, fh1[i], fh2[i], fo[i]);
    intramolecular(w1, energy);
    assert(not std::isnan(energy));
    DEBUG_LEVEL_RETURN(0)
    // Two next lines for interaction with platinum
    interactWithCorePt(w1, nPt, rPt, fPt,
                       energy); // called Vw-core in Zhu and Philpott
    assert(not std::isnan(energy));
    DEBUG_LEVEL_RETURN(1)
    interactWithImage(w1, w1, energy); // called Vw-cond in Zhu and Philpott
    DEBUG_LEVEL_RETURN(2)
    assert(not std::isnan(energy));
    for (int j = i - 1; j >= 0; --j) {
      bool areFixed = false;
      if (xh1 and xh2 and xo) {
        areFixed = xh1[i][0] and xh2[i][0] and xo[i][0];
        // check if all the atoms of molecule j are fixed.
        areFixed &= xh1[j][0] and xh2[j][0] and xo[j][0];
      };
      // if both molecules are fixed skip force calculation
      if (not areFixed) {
        double centre2[3] = {0};
        calculateCentre(rh1[j], rh2[j], ro[j], centre2);
        Water w2(rh1[j], rh2[j], ro[j], centre2, fh1[j], fh2[j], fo[j]);
        lennardJonesWithCutoff(w1, w2, energy);
        DEBUG_LEVEL_RETURN(3)
        coulombWithCutoff(w1, w2, energy, 1.0);
        DEBUG_LEVEL_RETURN(4)
        // Two next lines for interaction with platinum
        // called Vw-cond in Zhu and Philpott
        interactWithImage(w1, w2, energy); // w1 interacts with image of w2
        interactWithImage(w2, w1, energy); // w2 interacts with image of w1
        DEBUG_LEVEL_RETURN(5)
      };
    };
#ifndef NDEBUG
    for (int a = 0; a < 3; a++) {
      assert(not std::isnan(fh1[i][a]));
      assert(not std::isnan(fh2[i][a]));
      assert(not std::isnan(fo[i][a]));
    };
#endif
  };
  assert(not std::isnan(energy) and not std::isinf(energy));
}

computeTemplate() [2/2]

template<class P = zhu_philpott_parameters::Standard>

template<int H, int O, int H3, int O3>

void forcefields::ZhuPhilpott< P >::computeTemplate ( const int nWater, const double(*) rh1[H3], const double(*) rh2[H3], const double(*) ro[O3], double(*) fh1[H3], double(*) fh2[H3], double(*) fo[O3], const int nPt, const double rPt[][3], double fPt[][3], double & energy, double const b[], bool const (*) xh1[H] = 0, bool const (*) xh2[H] = 0, bool const (*) xo[O] = 0, bool const * xPt = 0 )

private

coulombFull()

template<class P>

void ZhuPhilpott::coulombFull ( Water & w1, Water & w2, double & U, double const relativePermittivity )

private

Coulomb interaction between two molecules with cutoff.

Coulomb interaction between two molecules of water.

Parameters

[in,out]	w1, w2	Water molecules.
[in,out]	U	Add the potential energy to U.
[in]	relativePermittivity

See also

getCutOff() and getSwitchingWidth().

Definition at line 525 of file zhu_philpott.cpp.

                                                                    {
  // q^2/relativePermittivity;
  const double qq2overEr = charge2_ / relativePermittivity;
  // Coulomb interactions between hydrogens
  coulomb(w1.rh1_, w2.rh1_, w1.fh1_, w2.fh1_, U, qq2overEr);
  coulomb(w1.rh1_, w2.rh2_, w1.fh1_, w2.fh2_, U, qq2overEr);
  coulomb(w1.rh2_, w2.rh1_, w1.fh2_, w2.fh1_, U, qq2overEr);
  coulomb(w1.rh2_, w2.rh2_, w1.fh2_, w2.fh2_, U, qq2overEr);
  // interactions between H and O.
  coulomb(w1.ro_, w2.rh1_, w1.fo_, w2.fh1_, U, -2.0 * qq2overEr);
  coulomb(w1.ro_, w2.rh2_, w1.fo_, w2.fh2_, U, -2.0 * qq2overEr);
  coulomb(w1.rh1_, w2.ro_, w1.fh1_, w2.fo_, U, -2.0 * qq2overEr);
  coulomb(w1.rh2_, w2.ro_, w1.fh2_, w2.fo_, U, -2.0 * qq2overEr);
  // interactions between O1, O2
  coulomb(w1.ro_, w2.ro_, w1.fo_, w2.fo_, U, 4.0 * qq2overEr);
}

coulombWithCutoff()

template<class P>

void ZhuPhilpott::coulombWithCutoff ( Water & w1, Water & w2, double & U, double const relativePermittivity )

private

Coulomb interaction between two molecules within the swithcing zone.

The precondition for using the function is that the distance between w1 and w2 defined as \( |\mathbf r_{n1}-\mathbf r_{n2}| \) verifies the equation:

\[ C_{utoff}-S_{witchingWidth} < |\mathbf r_{n1}-\mathbf r_{n2}| < C_{utoff} \]

Parameters

[in,out]	w1	Molecule 1.
[in,out]	w2	Molecule 2.
[in,out]	U	Add the potential energy to U.
[in]	relativePermittivity	Relative permittivity.

See also

getCutOff() and getSwitchingWidth().

Definition at line 496 of file zhu_philpott.cpp.

                                                                          {
  double z[3], z1, z2;
  distance(w1.rc_, w2.rc_, z, z1, z2);
  if (z1 <= cutoff_ - switchingWidth_) {
    coulombFull(w1, w2, U, relativePermittivity);
  } else if (z1 < cutoff_) {
    double f1[3][3] = {{0}}, f2[3][3] = {{0}};
    Water v1(w1, f1[0], f1[1], f1[2]);
    Water v2(w2, f2[0], f2[1], f2[2]);
    double energy = 0.0;
    coulombFull(v1, v2, energy, relativePermittivity);
    ChargeGroup<3> g1 = {v1.rc_, 0, f1};
    ChargeGroup<3> g2 = {v2.rc_, 0, f2};
    switching(g1, g2, energy, cutoff_, switchingWidth_);
    U += energy;
    w1.addForces(v1);
    w2.addForces(v2);
  };
}

getName()

template<class P>

char const * ZhuPhilpott::getName ( )

static

Name of the potential.

Definition at line 209 of file zhu_philpott.cpp.

                                                       {
  return "ZhuPhilpott";
}

interactionPtH()

template<class P>

void ZhuPhilpott::interactionPtH ( double const R1[], double const R2[], double F1[], double F2[], double & energy )

private

Definition at line 359 of file zhu_philpott.cpp.

                                                                              {
  double R12[3], F12[3] = {0}, r1, f1 = 0;
  distance(R1, R2, R12, r1);
  anisotropic(R12, F12, energy, P::epsilonHPt_, P::sigmaHPt_, P::alpha_);
  isotropic10(r1, f1, energy, P::epsilonHPt_, P::sigmaHPt_, P::C10_H_);
  for (int j = 0; j < 3; j++) {
    F1[j] += F12[j] + f1 * R12[j] / r1;
    F2[j] -= F12[j] + f1 * R12[j] / r1;
  };
}

interactionPtO()

template<class P>

void ZhuPhilpott::interactionPtO ( double const R1[], double const R2[], double F1[], double F2[], double & energy )

private

Definition at line 346 of file zhu_philpott.cpp.

                                                                              {
  double R12[3], F12[3] = {0}, r1, f1 = 0;
  distance(R1, R2, R12, r1);
  anisotropic(R12, F12, energy, P::epsilonOPt_, P::sigmaOPt_, P::alpha_);
  isotropic10(r1, f1, energy, P::epsilonOPt_, P::sigmaOPt_, P::C10_O_);
  for (int j = 0; j < 3; j++) {
    F1[j] += F12[j] + f1 * R12[j] / r1;
    F2[j] -= F12[j] + f1 * R12[j] / r1;
  };
}

interactWithCorePt()

template<class P>

void ZhuPhilpott::interactWithCorePt ( Water & water, int const nPt, double const rPt[][3], double fPt[][3], double & energy )

private

Interaction of one molecule of water with the whole platinum.

Parameters

[in,out]	water	Molecule of water. The forces are incremented.
[in]	nPt	Number of Platinum atom.
[in]	rPt	Positions of the platinum.
[in,out]	fPt	Forces on the platinum. The forces are incremented
[in,out]	energy	Add the energy.

Definition at line 292 of file zhu_philpott.cpp.

                                                        {
  for (int i = 0; i < nPt; ++i) {
    double r[3], r1; // for distance between two atoms
                     // Pt - O
    distance(water.ro_, rPt[i], r, r1);
    if (r1 <= cutoff_ - switchingWidth_) {
      interactionPtO(water.ro_, rPt[i], water.fo_, fPt[i], energy);
    } else if (r1 < cutoff_) {
      double en = 0;
      double fo[3] = {0}, fPtTmp[3] = {0};
      interactionPtO(water.ro_, rPt[i], fo, fPtTmp, en);
      switching(water.ro_, rPt[i], fo, fPtTmp, en);
      energy += en;
      for (int j = 0; j < 3; j++) {
        water.fo_[j] += fo[j];
        fPt[i][j] += fPtTmp[j];
      };
    };
    // Pt - H1
    distance(water.rh1_, rPt[i], r, r1);
    if (r1 <= cutoff_ - switchingWidth_) {
      interactionPtH(water.rh1_, rPt[i], water.fh1_, fPt[i], energy);
    } else if (r1 < cutoff_) {
      double en = 0;
      double fh[3] = {0}, fPtTmp[3] = {0};
      interactionPtH(water.rh1_, rPt[i], fh, fPtTmp, en);
      switching(water.rh1_, rPt[i], fh, fPtTmp, en);
      energy += en;
      for (int j = 0; j < 3; j++) {
        water.fh1_[j] += fh[j];
        fPt[i][j] += fPtTmp[j];
      };
    };
    // Pt - H2
    distance(water.rh2_, rPt[i], r, r1);
    if (r1 <= cutoff_ - switchingWidth_) {
      interactionPtH(water.rh2_, rPt[i], water.fh2_, fPt[i], energy);
    } else if (r1 < cutoff_) {
      double en = 0;
      double fh[3] = {0}, fPtTmp[3] = {0};
      interactionPtH(water.rh2_, rPt[i], fh, fPtTmp, en);
      switching(water.rh2_, rPt[i], fh, fPtTmp, en);
      energy += en;
      for (int j = 0; j < 3; j++) {
        water.fh2_[j] += fh[j];
        fPt[i][j] += fPtTmp[j];
      };
    };
  };
}

interactWithImage()

template<class P>

void ZhuPhilpott::interactWithImage ( Water & w1, Water & w2, double & U )

private

Interactions of a molecule with an image.

Coulomb interaction between w1 and the image of w2. Includes cutoff.

Parameters

[in,out]	w1	Molecule 1.
[in,out]	w2	Molecule 2.
[in,out]	U	Add the potential energy to U.

Note

Does not include the interaction of w2 with the image of w1. For this you must explicitely call interactWithImage(w2, w1, U).

Definition at line 449 of file zhu_philpott.cpp.

                                                                      {
  // f is used to store the force on the atoms' images. These forces have no use
  // and are at the end discarded. we use the same vector to store of the forces
  // even when they are on different atoms' images.
  double r[4][3], f[3][3] = {{0}};
  for (int k = 0; k < 2; k++) {
    r[0][k] = w2.rh1_[k];
    r[1][k] = w2.rh2_[k];
    r[2][k] = w2.ro_[k];
    r[3][k] = w2.rc_[k];
  };
  // mirror about plane at z=0.
  r[0][2] = -w2.rh1_[2];
  r[1][2] = -w2.rh2_[2];
  r[2][2] = -w2.ro_[2];
  r[3][2] = -w2.rc_[2];
  Water w2m(r[0], r[1], r[2], r[3], f[0], f[1], f[2]);
  // -2 is given as the relative permittivity. This value is not the
  // permittivity of the metal which is infinite in the case of a perfect metal.
  // It is the value (1+epsilon_r)/(1-epsilon_r) which is equal to -1 for a
  // perfect metal. This value is the 'effective' permittivity for the
  // interaction of a charge with a mirror image.
  coulombWithCutoff(w1, w2m, U, -2.0);
  for (int k = 0; k < 2; k++) {
    w2.fh1_[k] += f[0][k];
    w2.fh2_[k] += f[1][k];
    w2.fo_[k] += f[2][k];
  };
  // mirror about plane
  w2.fh1_[2] -= f[0][2];
  w2.fh2_[2] -= f[1][2];
  w2.fo_[2] -= f[2][2];
}

isotropic10()

template<class P>

void ZhuPhilpott::isotropic10 ( double const distance, double & force, double & energy, double const epsilon, double const sigma, double const C10 )

private

Isotropic interaction between water and platinum.

The isotropic iteraction between a atom Pt and an atom O or H of a molecule of water has the functional form:

\[ E=-4\epsilon\ C_{10}\left(\frac{\sigma}{r}\right)^{10} \]

Parameters

[in]	distance	Distance between two atoms.
[in,out]	force	Force create by the interaction between the two atoms.
[in,out]	energy	Energy created by the interaction.
[in]	epsilon	See equation.
[in]	sigma	See equation.
[in]	C10	See equation.

Definition at line 430 of file zhu_philpott.cpp.

                                                                       {
  double const s = sigma / distance;
  double s10 = s * s * s * s * s;
  s10 *= s10;
  energy += -4.0 * epsilon * C10 * s10;
  force -= 40.0 * epsilon * C10 * s10 / distance;
}

nPlatinum()

template<class P>

int ZhuPhilpott::nPlatinum

(

)

const

Number of platinum atoms.

See also

setPlatinum().

Definition at line 194 of file zhu_philpott.cpp.

{ return nPlatinum_; }

operator=()

template<class P>

void ZhuPhilpott::operator=

(

ZhuPhilpott< P > const &

zhuPhilpott

)

Definition at line 86 of file zhu_philpott.cpp.

                                                           {
  setPlatinum(zhuPhilpott.nPlatinum_, zhuPhilpott.positions_.data());
}

setPlatinum()

template<class P>

void ZhuPhilpott::setPlatinum	(	int	nPlatinum,
		double const	positions[] )

Initialises the positions of the atoms of platinum.

Use before calling computeHH_O_().

Definition at line 200 of file zhu_philpott.cpp.

                                                                        {
  assert(nPlatinum >= 0);
  nPlatinum_ = nPlatinum;
  int const n = nPlatinum * 3;
  positions_.assign(positions, positions + n);
  forces_.assign(n, 0.0);
}

Member Data Documentation

forces_

template<class P = zhu_philpott_parameters::Standard>

std::vector<double> forcefields::ZhuPhilpott< P >::forces_

private

Definition at line 70 of file zhu_philpott.hpp.

nPlatinum_

template<class P = zhu_philpott_parameters::Standard>

int forcefields::ZhuPhilpott< P >::nPlatinum_ {0}

private

Definition at line 68 of file zhu_philpott.hpp.

{0};

positions_

template<class P = zhu_philpott_parameters::Standard>

std::vector<double> forcefields::ZhuPhilpott< P >::positions_

private

Definition at line 69 of file zhu_philpott.hpp.

---

The documentation for this class was generated from the following files:

• /home/runner/work/eOn/eOn/include/eon/potentials/Water_Pt/zhu_philpott.hpp
• /home/runner/work/eOn/eOn/client/potentials/Water_Pt/zhu_philpott.cpp