forcefields::SpceCcl Class Reference

SPC/E+CCl water potential. More...

#include <spce_ccl.hpp>

Inheritance diagram for forcefields::SpceCcl:

Classes
struct	Water
	Pointers to molecule of water. More...

Public Member Functions
	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)

Protected Member Functions
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 Attributes
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

Private Member Functions
void	initialiseRho (Vector3 const &v, Rho &r)
	Initialise Rho.
template<int H, int O>
void	computeTemplate (const int nMolecules, const double(*const rh1)[H *3], const double(*const rh2)[H *3], const double(*const ro)[O *3], double(*const fh1)[H *3], double(*const fh2)[H *3], double(*const fo)[O *3], double &energy, double const b[], bool const (*const xh1)[H]=0, bool const (*const xh2)[H]=0, bool const (*const xo)[O]=0)
	Compute forces and energy.

Additional Inherited Members
static double	applyPeriodicity0 (double r, double const period)
	Minimum image representation.
static void	applyPeriodicity0 (double r[], double const periods[])
	Minimum image representation.
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_

Detailed Description

SPC/E+CCl water potential.

This potential combines potential SPC/E (for intermolecular interactions) with potential CCL (for intramolecular interactions). The functions and parameters for the SPC/E are from Berendsen, et al.. References for potential CCL can be found in in class Ccl 's documentation. The rules used to merge these two potentials were taken from Amira, et al.

References

The Missing Term in Effective Pair Potentials. HJC Berendsen, JR Grigera et al. J. Phys. Chem. 1987, vol. 91, p. 6269-6271
Derivation and evaluation of a flexible SPC model for liquid water. S Amira, D Spangberg, K Hermansson, Chemical Physics, 2004, vol. 303, p. 327-334

Definition at line 36 of file spce_ccl.hpp.

Constructor & Destructor Documentation

SpceCcl() [1/2]

SpceCcl::SpceCcl

(

)

Definition at line 56 of file spce_ccl.cpp.

    : Ccl() {}

SpceCcl() [2/2]

SpceCcl::SpceCcl	(	double	cutoff,
		double	switchingWidth )

Definition at line 59 of file spce_ccl.cpp.

    : Ccl(cutoff, switchingWidth) {}

Member Function Documentation

computeHH_O_() [1/2]

void SpceCcl::computeHH_O_	(	const int	nAtoms,
		const double	R[],
		double	F[],
		double &	U,
		const double	b[] )

Compute the forces and the energy.

The order of the atoms is very important. For this function the order is H1, H1, H2, H2, etc ... , O1, O2, etc ... The numbers are for the molecules. The suffix HH_O_ is to reminds the ordering of the atoms. It indicates the position of the atoms of the first molecule, while the underscore marks the locations of other atoms.

Parameters

[in]	nAtoms	Number of Atoms.
[in]	R	coordinates.
[out]	F	Forces.
[out]	U	Potential energy.
[in]	b	Periodic boundaries.

Warning

Be careful with the order of the atoms.

Definition at line 62 of file spce_ccl.cpp.

                                                        {
  computeHH_O_(nAtoms, R, F, U, b, 0);
}

computeHH_O_() [2/2]

void SpceCcl::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).

In a simulation, some atoms may be fixed (not allowed to move). Computing the interaction between two fixed atoms is useless. When the potential knows which atoms are fixed, it can avoid these useless computations.

Parameters

[in]	nAtoms	Number of Atoms.
[in]	R	coordinates.
[out]	F	Forces.
[out]	U	Potential energy.
[in]	b	Periodic boundaries.
[in]	fixed	The length of the array must be equal to nAtoms in compute(). True when atom is fixed, false otherwise.

See also

Check compute() for the order of the atoms.

Definition at line 67 of file spce_ccl.cpp.

                                                                            {
  int const nMolecules = nAtoms / 3;
  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[nMolecules * 6]);
  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[nMolecules * 6]);
  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[nMolecules * 2]);
  computeTemplate(nMolecules, rh1, rh2, ro, fh1, fh2, fo, U, b, xh1, xh2, xo);
}

computeTemplate()

template<int H, int O>

void SpceCcl::computeTemplate ( const int nMolecules, const double(*) rh1[H *3], const double(*) rh2[H *3], const double(*) ro[O *3], double(*) fh1[H *3], double(*) fh2[H *3], double(*) fo[O *3], double & energy, double const b[], bool const (*) xh1[H] = 0, bool const (*) xh2[H] = 0, bool const (*) xo[O] = 0 )

private

Compute forces and energy.

The function calculates and returns the energy and forces applied on each atom. The function has preconditions with which the user must comply before calling this function. The periodic boundaries must set with (setPeriodicity()). The content of and arrays fh1 , fh2 , fo must be zero for all elements.
Arrays ending in h1 are for the hydrogen one of the molecules and those ending in h2 for hydrogen 2.

Parameters

[in]	nMolecules	Number of molecules.
[in]	rh1, rh2, ro	Positions of hydrogens and oxygens.
[in,out]	fh1, fh2, fo	Forces on hydrogens and oxygens.
[in,out]	energy	Potential energy.
[in]	b	Periodic boundaries.
[in]	xh1, xh2, xo	Tell which atoms are fixed and which are movable. This parameter are optional. When provided, interaction between fixed atoms are skipped. The three arrays xh1, xh2, xo must be provided for the optimisation to work.

Warning

Remember the preconditions.

Definition at line 175 of file spce_ccl.cpp.

                               {
  for (int i = 0; i < nMolecules; ++i) {
    for (int a = 0; a < 3; a++) {
      fh1[i][a] = 0.0;
      fh2[i][a] = 0.0;
      fo[i][a] = 0.0;
    };
  };
  energy = 0.0;
  setPeriodicity(b);
 
  for (int i = nMolecules - 1; i >= 0; --i) {
    double rc1[3] = {0};
    calculateCentre(rh1[i], rh2[i], ro[i], rc1);
    Water w1(rh1[i], rh2[i], ro[i], rc1, fh1[i], fh2[i], fo[i]);
    intramolecular(w1, energy);
    energy += polarisationEnergy_;
    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 rc2[3] = {0};
        calculateCentre(rh1[j], rh2[j], ro[j], rc2);
        Water w2(rh1[j], rh2[j], ro[j], rc2, fh1[j], fh2[j], fo[j]);
        lennardJonesWithCutoff(w1, w2, energy);
        coulombWithCutoff(w1, w2, energy);
      };
    };
  };
  assert(not std::isnan(energy) and not std::isinf(energy));
}

coulombFull()

void SpceCcl::coulombFull ( Water & w1, Water & w2, double & U )

protected

Interactions between two molecules.

Coulomb interaction between two molecules. Full interaction (i.e. no cutoff).

Parameters

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

See also

intermolecularFull() and intermolecularSwitching().

Definition at line 133 of file spce_ccl.cpp.

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

coulombWithCutoff()

void SpceCcl::coulombWithCutoff ( Water & w1, Water & w2, double & U )

protected

Interactions between two molecules.

Coulomb interaction between two molecules of water with molecules-based cutoff

Parameters

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

See also

intermolecularFull() and intermolecularSwitching().

Definition at line 104 of file spce_ccl.cpp.

                                                               {
  double z[3], z1, z2;
  distance(w1.rc_, w2.rc_, z, z1, z2);
  if (z1 <= cutoff_ - switchingWidth_) {
    coulombFull(w1, w2, U);
  } else if (z1 < cutoff_) {
    // store forces and energy of the full interaction in temporaries
    double f1[3][3] = {{0}}, f2[3][3] = {{0}};
    Water v1(w1.rh1_, w1.rh2_, w1.ro_, w1.rc_, f1[0], f1[1], f1[2]);
    Water v2(w2.rh1_, w2.rh2_, w2.ro_, w2.rc_, f2[0], f2[1], f2[2]);
    double energy = 0.0;
    coulombFull(v1, v2, energy);
    ChargeGroup<3> g1 = {v1.rc_, 0, f1};
    ChargeGroup<3> g2 = {v2.rc_, 0, f2};
    // Calculate the weakened forces and energy.
    switching(g1, g2, energy, cutoff_, switchingWidth_);
    // add weakened force and energy to those of other interactions.
    U += energy;
    for (int i = 0; i < 3; ++i) {
      w1.fh1_[i] += v1.fh1_[i];
      w1.fh2_[i] += v1.fh2_[i];
      w1.fo_[i] += v1.fo_[i];
      w2.fh1_[i] += v2.fh1_[i];
      w2.fh2_[i] += v2.fh2_[i];
      w2.fo_[i] += v2.fo_[i];
    };
  };
}

getName()

char const * SpceCcl::getName

(

)

const

Name of the potential.

Definition at line 84 of file spce_ccl.cpp.

{ return "SpceCcl"; }

initialiseRho()

void SpceCcl::initialiseRho ( Vector3 const & v, Rho & r )

private

Initialise Rho.

Compute Rho coordinates according to Amira's grafting rules:

\[ \rho^{SPC+CCL}=\frac{r-r^{SPC}_e}{r-r^{SPC}_e+r^{CCL}_e} \]

\[ \frac{d\rho}{d\mathbf r}=\frac{-r_e^{CCL}}{\left(r-r^{SPC}_e+r^{CCL}_e\right)^2}\frac{\mathbf r}{r} \]

Definition at line 148 of file spce_ccl.cpp.

                                                     {
  // re_ = r_CCL and roh_ = r_SPC
  ro._1 = (v._1 - roh_) / (v._1 - roh_ + re_);
  ro._2 = ro._1 * ro._1;
  ro._3 = ro._2 * ro._1;
  double const a = v._1 - roh_ + re_;
  double const d = re_ / a / a / v._1;
  for (int i = 0; i < 3; ++i)
    ro.n[i] = v.v[i] * d;
}

intramolecular()

void SpceCcl::intramolecular ( Water & water, double & U )

protected

Interactions within a molecules.

CCL interaction with Amira's grafting rules.

Parameters

[in,out]	water	Add forces.
[in,out]	U	Add energy to U.

Definition at line 159 of file spce_ccl.cpp.

                                                     {
  Vector3 v1, v2;
  distance(w.rh1_, w.ro_, v1);
  distance(w.rh2_, w.ro_, v2);
  // ------ prepare --------------------
  Rho ro1, ro2;
  initialiseRho(v1, ro1);
  initialiseRho(v2, ro2);
 
  // ------- prepare Delta theta
  Dtheta dth;
  initialiseDtheta(v1, v2, dth, theta_); // theta_ SPC equilibrium angle.
  Ccl::intramolecular(ro1, ro2, dth, energy, w.fh1_, w.fh2_, w.fo_);
}

lennardJonesWithCutoff()

void SpceCcl::lennardJonesWithCutoff ( Water & w1, Water & w2, double & U )

protected

Interactions between two molecules.

Lennard Jones between oxygen only with cutoff

Parameters

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

See also

intermolecularFull() and intermolecularSwitching().

Definition at line 86 of file spce_ccl.cpp.

                                                                    {
  double z[3], z1;
  distance(w1.ro_, w2.ro_, z, z1);
  if (z1 <= cutoff_ - switchingWidth_) {
    lennardJones(w1.ro_, w2.ro_, w1.fo_, w2.fo_, U, epsilon_, sigma_);
  } else if (z1 < cutoff_) {
    double f1[3] = {0}, f2[3] = {0};
    double energy = 0.0;
    lennardJones(w1.ro_, w2.ro_, f1, f2, energy, epsilon_, sigma_);
    switching(w1.ro_, w2.ro_, f1, f2, energy);
    U += energy;
    for (int i = 0; i < 3; i++) {
      w1.fo_[i] += f1[i];
      w2.fo_[i] += f2[i];
    };
  };
}

Member Data Documentation

A_

const double SpceCcl::A_ = 27291.6

staticprotected

Lennard-Jones.

See PotentialBase::sigma() and PotentialBase::epsilon() for definition

Definition at line 160 of file spce_ccl.hpp.

B_

const double SpceCcl::B_ = 27.1223

staticprotected

Lennard-Jones.

See PotentialBase::sigma() and PotentialBase::epsilon() for definition

Definition at line 163 of file spce_ccl.hpp.

charge2_

const double SpceCcl::charge2_ = 0.179606

staticprotected

Square of # charge_.

Definition at line 157 of file spce_ccl.hpp.

charge_

const double SpceCcl::charge_ = 0.4238

staticprotected

Charge on one hydrogen.

Definition at line 155 of file spce_ccl.hpp.

epsilon_

const double SpceCcl::epsilon_ = 0.00673853

staticprotected

Lennard-Jones. See PotentialBase::lennardJones() for definition.

Definition at line 167 of file spce_ccl.hpp.

polarisationEnergy_

const double SpceCcl::polarisationEnergy_ = 0.0541015

staticprotected

Polarisation correction.

Energy added per molecule (see original publication for definition).

Definition at line 170 of file spce_ccl.hpp.

rhh_

const double SpceCcl::rhh_ = 1.63299

staticprotected

Distance HH.

Definition at line 153 of file spce_ccl.hpp.

roh_

const double SpceCcl::roh_ = 1.0

staticprotected

Distance OH.

Definition at line 149 of file spce_ccl.hpp.

sigma_

const double SpceCcl::sigma_ = 3.16556

staticprotected

Lennard-Jones. See PotentialBase::lennardJones() for definition.

Definition at line 165 of file spce_ccl.hpp.

theta_

const double SpceCcl::theta_ = 1.91063

staticprotected

Angle HOH.

Definition at line 151 of file spce_ccl.hpp.

---

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

• /home/runner/work/eOn/eOn/include/eon/potentials/Water/spce_ccl.hpp
• /home/runner/work/eOn/eOn/client/potentials/Water/spce_ccl.cpp