forcefields::Tip4p Class Reference

TIP4P forcefield with fexible molecules and group based cutoff. More...

#include <tip4p_ccl.hpp>

Inheritance diagram for forcefields::Tip4p:

Classes
struct	Water
	Pointers for one molecule of water. More...

Public Member Functions
	Tip4p ()
	Tip4p (double cutoff, double switchingWidth)
	~Tip4p ()
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).
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>
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.
void	coulombWithCutoff (Water &w1, Water &w2, double &U)
	Interactions between two molecules.
void	coulombFull (Water &w1, Water &w2, double &U)
	Interactions between two molecules.
void	lennardJonesWithCutoff (Water &w1, Water &w2, double &U)
	Interactions between two molecules.

Additional Inherited Members
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::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

Detailed Description

TIP4P forcefield with fexible molecules and group based cutoff.

The function computing the forces and energy is compute(). Parameters for the forcefield were taken from Abascal et al.. The parameters are those of TIP4P (original version). TIP4P is a potential for constrained molecules of water. In this implementation quadratic restraints were added to bonds OH and HH inside molecules of water. So the potential can be used without constraints. The potential has a cutoff for long range interaction. The cutoff can be changed by setCutoff(). There is also and switching zone at the edge of the cutoff to cut off the interactions smoothly. The width of this switching zone is controlled by setSwitchingWidth().
The unit system for this potential is eV (electron volt), Angstrom, e (e charge).

References

A general purpose model for the condensed phases of water: TIP4P/2005, J.L.F. Abascal and C. Vega, J. Chem. Phys. (2005) vol. 123, p. 234505.

Definition at line 23 of file tip4p_ccl.hpp.

Constructor & Destructor Documentation

Tip4p() [1/2]

Tip4p::Tip4p

(

)

Definition at line 66 of file tip4p_ccl.cpp.

    : Ccl() {}

Tip4p() [2/2]

Tip4p::Tip4p	(	double	cutoff,
		double	switchingWidth )

Definition at line 69 of file tip4p_ccl.cpp.

    : Ccl(cutoff, switchingWidth) {}

~Tip4p()

forcefields::Tip4p::~Tip4p ( )

inline

Definition at line 27 of file tip4p_ccl.hpp.

{}

Member Function Documentation

computeHH_O_() [1/2]

void Tip4p::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 84 of file tip4p_ccl.cpp.

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

computeHH_O_() [2/2]

void Tip4p::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 102 of file tip4p_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 Tip4p::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 156 of file tip4p_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], rn1[3], fn1[3] = {0}; // rn position charge N, rc centre of
                                         // charge for groupbased cutoff.
    calculateWeightedCentre(wo_, wh_, wh_, ro[i], rh1[i], rh2[i], rn1);
    calculateCentre(rn1, rh1[i], rh2[i], rc1);
    Water w1 = {rh1[i], rh2[i], ro[i], rn1, rc1, fh1[i], fh2[i], fo[i], fn1};
    intramolecular(rh1[i], rh2[i], ro[i], fh1[i], fh2[i], fo[i], 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 rc2[3], rn2[3], fn2[3] = {0};
        calculateWeightedCentre(wo_, wh_, wh_, ro[j], rh1[j], rh2[j], rn2);
        calculateCentre(rn2, rh1[j], rh2[j], rc2);
        Water w2 = {rh1[j], rh2[j], ro[j], rn2, rc2,
                    fh1[j], fh2[j], fo[j], fn2};
        coulombWithCutoff(w1, w2, energy);
        spreadWeightedForce(wo_, wh_, wh_, w2.fo_, w2.fh1_, w2.fh2_, w2.fn_);
        lennardJonesWithCutoff(w1, w2, energy);
      };
    };
    spreadWeightedForce(wo_, wh_, wh_, w1.fo_, w1.fh1_, w1.fh2_, w1.fn_);
  };
  assert(not std::isnan(energy) and not std::isinf(energy));
}

coulombFull()

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

private

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().

Note

The function add the force on charge N to Water::fn_ but does not spread it to Water::fh1_, etc ...

See also

spreadN().

Definition at line 251 of file tip4p_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 N.
  coulomb(w1.rn_, w2.rh1_, w1.fn_, w2.fh1_, U, -2.0 * charge2_);
  coulomb(w1.rn_, w2.rh2_, w1.fn_, w2.fh2_, U, -2.0 * charge2_);
  coulomb(w1.rh1_, w2.rn_, w1.fh1_, w2.fn_, U, -2.0 * charge2_);
  coulomb(w1.rh2_, w2.rn_, w1.fh2_, w2.fn_, U, -2.0 * charge2_);
  // interactions between N1, N2
  coulomb(w1.rn_, w2.rn_, w1.fn_, w2.fn_, U, 4.0 * charge2_);
}

coulombWithCutoff()

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

private

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.

Definition at line 209 of file tip4p_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_) {
    double f1[3][3] = {{0}}, f2[3][3] = {{0}};
    // f1[0], f1[1], 0, f1[2] are respectively, H1, H2, O, N. There is no charge
    // on O.
    Water v1 = {w1.rh1_, w1.rh2_, w1.ro_, w1.rn_, w1.rc_,
                f1[0],   f1[1],   0,      f1[2]};
    Water v2 = {w2.rh1_, w2.rh2_, w2.ro_, w2.rn_, w2.rc_,
                f2[0],   f2[1],   0,      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.fn_[i] += v1.fn_[i];
      w2.fh1_[i] += v2.fh1_[i];
      w2.fh2_[i] += v2.fh2_[i];
      w2.fn_[i] += v2.fn_[i];
    };
  }
}

getName()

char const * Tip4p::getName ( )

static

Name of the potential.

Definition at line 120 of file tip4p_ccl.cpp.

{ return "Tip4p"; }

lennardJonesWithCutoff()

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

private

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 273 of file tip4p_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];
    };
  };
}

---

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

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