forcefields::PotentialBase Class Reference

Functions and tools commonly used by potentials. More...

#include <potential_base.hpp>

Inheritance diagram for forcefields::PotentialBase:

Classes
struct	ChargeGroup
struct	Vector3

Public Member Functions
	PotentialBase ()
	Non bond interaction cutoff.
	PotentialBase (double cutoff, double switchingWidth)
virtual	~PotentialBase ()
double	getCutoff () const
void	setCutoff (double cutoff)

Protected Member Functions
Angle
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.

Static Protected Member Functions
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|} \).

Protected Attributes
double	cutoff_
double	periods_ [3]
double	switchingWidth_

Coulomb
static const double	ONE_OVER_4_PI_EPSILON0
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)
static const double	EPSILON_PT = 0.68165797577788501
	Platinum Lennard-Jones.
static const double	SIGMA_PT = 2.54
double	getSwitchingWidth () const
	Width of the switching zone.
void	setSwitchingWidth (double width)
static double	applyPeriodicity0 (double r, double const period)
	Minimum image representation.
static void	applyPeriodicity0 (double r[], double const periods[])
	Minimum image representation.
void	computePt (int const nAtoms, double positions[], double forces[], double &energy, double const periods[], bool const fixed[])
	Potential for Platinum.

Bond Interactions
static double	isoscelesBase (double length, double angle)
	Calculate the base of an isosceles triangle.
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.

Lennard-Jones
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.
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.

Periodic Boundaries, Distance...
3-D vector.
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.
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.

Detailed Description

Functions and tools commonly used by potentials.

Contains functions for quadratic retrainst, Lennard-Jones, Coulomb, to manage the periodic boundaries and interface for cutoff. The system of units used is (eV, Angstrom, fs, e). Rules to combine different Lennard-Jones parameters were taken from:
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 24 of file potential_base.hpp.

Constructor & Destructor Documentation

PotentialBase() [1/2]

PotentialBase::PotentialBase

(

)

Non bond interaction cutoff.

When the distance between two molecules is over getCutoff(), van der Waals and Coulomb interactions between the two molecules are ignored.

See also

getSwitchingWidth().

Definition at line 40 of file potential_base.cpp.

    : cutoff_(6.5),
      switchingWidth_(2.0) {
  periods_[0] = 0.0;
  periods_[1] = 0.0;
  periods_[2] = 0.0;
}

PotentialBase() [2/2]

PotentialBase::PotentialBase	(	double	cutoff,
		double	switchingWidth )

Definition at line 48 of file potential_base.cpp.

    : cutoff_(cutoff),
      switchingWidth_(switchingWidth) {
  if (switchingWidth_ > cutoff_) {
    std::cerr << "Error: getSwitchingWidth() > getCutoff()" << std::endl;
    std::exit(EXIT_FAILURE);
  };
  periods_[0] = 0.0;
  periods_[1] = 0.0;
  periods_[2] = 0.0;
}

~PotentialBase()

virtual forcefields::PotentialBase::~PotentialBase ( )

inlinevirtual

Definition at line 28 of file potential_base.hpp.

{};

Member Function Documentation

addForces()

template<int N, class R, class F>

void forcefields::PotentialBase::addForces ( ChargeGroup< N, R, F > const & g1, ChargeGroup< N, R, F > & g2 )

protected

Increment forces.

Increment forces in g2 by forces in g1.

Definition at line 157 of file potential_base.hpp.

                                                                     {
  increment(g2.force_, g1.force_);
}

applyPeriodicity0() [1/2]

double PotentialBase::applyPeriodicity0 ( double r, double const period )

static

Minimum image representation.

Returns

Minimum image representation of r. The value returned is \(\frac{-period}{2} \le r \le \frac{period}{2} \)

Parameters

[in]	r	Cyclic coordinate.
[in]	period	Period of the coordinates.

Definition at line 85 of file potential_base.cpp.

                                                                     {
  if (not std::isinf(period)) {
    double n = r / period + 0.5;
    // This is slightly more efficient than using function floor().
    int m = static_cast<int>(n);
    if (n < 0)
      --m;
    r -= m * period;
  };
  return r;
}

applyPeriodicity0() [2/2]

void PotentialBase::applyPeriodicity0 ( double r[], double const periods[] )

static

Minimum image representation.

Returns

Minimum image representation of r. The value returned is \(\frac{-period}{2} \le r \le \frac{period}{2} \).

Parameters

[in]	r	Three-dimension vector.
[in]	periods	Three-dimension vector containing periods along each of the axes.

Definition at line 104 of file potential_base.cpp.

                                                                        {
  for (int i = 0; i < 3; ++i)
    r[i] = applyPeriodicity0(r[i], periods[i]);
}

applyPeriodicity1() [1/2]

double PotentialBase::applyPeriodicity1 ( double r, int const axis )

protected

Minimum image representation.

Returns

Minimum image representation of r. The value returned is \(\frac{-period}{2} \le r \le \frac{period}{2} \). The periods should be set with setPeriodicity();

Parameters

[in]	r	Cyclic coordinate.
[in]	axis	Either 0, 1 or 2.

See also

setPeriodicity().

Definition at line 545 of file potential_base.cpp.

                                                                {
  return applyPeriodicity0(r, periods_[axis]);
}

applyPeriodicity1() [2/2]

void PotentialBase::applyPeriodicity1 ( double r[] )

protected

Minimum image representation.

Returns

Minimum image representation of r. The value returned is \(\frac{-period}{2} \le r \le \frac{period}{2} \). The periods should be set with setPeriodicity();

Parameters

[in]

r

Three-dimension vector.

See also

setPeriodicity().

Definition at line 534 of file potential_base.cpp.

                                                {
  applyPeriodicity0(r, periods_);
}

calculateCentre() [1/2]

void PotentialBase::calculateCentre ( double const r1[], double const r2[], double const r3[], double rc[] )

protected

Calculate centre of three points.

See also

calculateCentre().

Definition at line 196 of file potential_base.cpp.

                                                                    {
  for (int i = 0; i < 3; ++i) {
    rc[i] =
        (r1[i] + unBreak1(r2[i], r1[i], i) + unBreak1(r3[i], r1[i], i)) / 3.0;
  };
}

calculateCentre() [2/2]

void PotentialBase::calculateCentre ( double const r1[], double const r2[], double rc[] )

protected

Calculate centre of two points.

Calculate the centre of r1 and r2 and store the result in rc.

Note

The centre is calculated for the image or r2 the closed to r1.

Definition at line 188 of file potential_base.cpp.

                                                 {
  for (int i = 0; i < 3; ++i) {
    rc[i] = (r1[i] + unBreak1(r2[i], r1[i], i)) * 0.5;
  };
}

calculateWeightedCentre()

void PotentialBase::calculateWeightedCentre ( double const w1, double const w2, double const w3, double const r1[], double const r2[], double const r3[], double rc[] )

protected

Calculate barycentre.

The inverse of this functions the forces is spreadWeightedForce()

Precondition

Sum of weigths w1, w2, w3 must be 1.

Definition at line 208 of file potential_base.cpp.

                                                                            {
  assert(std::fabs(w1 + w2 + w3 - 1.0) < 1e-9);
  for (int i = 0; i < 3; ++i) {
    rc[i] = w1 * r1[i] + w2 * unBreak1(r2[i], r1[i], i) +
            w3 * unBreak1(r3[i], r1[i], i);
  };
}

computePt()

void PotentialBase::computePt ( int const nAtoms, double positions[], double forces[], double & energy, double const periods[], bool const fixed[] )

protected

Potential for Platinum.

Simple 12-6 Lennard-Jones potential for Platinum.
Parameters from: C. Kittel, Introduction to Solid State Physics, (Wiley, New York, 1986).

Parameters

[in]	nAtoms	Number of atoms.
[in]	positions	Positions of atoms.
[out]	forces	Forces on atoms.
[out]	energy	Potential energy.
[in]	periods	Periods for periodic boundary conditions.
[in]	fixed	Array of boolean of length equal to nAtoms. The array shall contain true is the atom is fixed, false is movable.

Definition at line 658 of file potential_base.cpp.

                                                                          {
  int const nCoord = 3 * nAtoms;
  for (int i = 0; i < nCoord; ++i)
    forces[i] = 0.0;
  energy = 0.0;
  double const(*r)[3] = reinterpret_cast<double const(*)[3]>(positions);
  double (*f)[3] = reinterpret_cast<double (*)[3]>(forces);
  setPeriodicity(periods); // Set Periodic boundaries. Essential in order for
                           // some functions to work.
  for (int i = nAtoms - 1; i > 0; --i) {
    for (int j = i - 1; j >= 0; --j) {
      if (not fixed[i] and not fixed[j]) {
        lennardJonesWithCutoff(r[i], r[j], f[i], f[j], energy, EPSILON_PT,
                               SIGMA_PT);
      };
    };
  };
}

coulomb() [1/2]

void PotentialBase::coulomb ( const double r1[], const double r2[], double f1[], double f2[], double & energy, double const qq )

protected

Compute Coulomb interaction between two charges.

\( U = \frac{q_1q_2}{|\mathbf r_1- \mathbf r_2|} \)

Parameters

[in]	r1	Positions of atom 1.
[in]	r2	Positions of atom 2.
[in,out]	f1	Incremented by force on atom 1.
[in,out]	f2	Incremented by force on atom 2.
[in,out]	energy	Incremented by energy.
[in]	qq	Product of the charges.

Warning

f1, f2 and u are incremented.

Definition at line 256 of file potential_base.cpp.

                                                                    {
  double d1, d[3], f, en;
  // d: vector distance (r1-r2), d1 distance and d2 squared distance
  distance(r1, r2, d, d1);
  coulomb(d1, f, en, qq);
  for (int l = 0; l < 3; ++l) {
    f1[l] += f * d[l] / d1;
    f2[l] -= f * d[l] / d1;
  };
  energy += en;
}

coulomb() [2/2]

void PotentialBase::coulomb ( double distance, double & force, double & energy, double const qq )

protected

Compute Coulomb interaction between two charges.

\( U = \frac{q_1q_2}{|\mathbf r_1- \mathbf r_2|} \)

Parameters

[in]	distance	between two charges
[out]	force	Force Resulting from the coulom interaction.
[out]	energy	Potential energy.
[in]	qq	Product of the charges.

Definition at line 275 of file potential_base.cpp.

                                             {
  energy = ONE_OVER_4_PI_EPSILON0 * qq / distance;
  force = energy / distance;
}

coulombWithCutoff()

void PotentialBase::coulombWithCutoff ( const double r1[], const double r2[], double f1[], double f2[], double & energy, double const qq )

protected

Coulomb interaction with single charge based cutoff.

\( U = \frac{q_1q_2}{|\mathbf r_1- \mathbf r_2|} \)

Parameters

[in]	r1	Positions of atom 1.
[in]	r2	Positions of atom 2.
[in,out]	f1	Incremented by force on atom 1.
[in,out]	f2	Incremented by force on atom 2.
[in,out]	energy	Incremented by energy.
[in]	qq	Product of the charges.

Warning

f1, f2 and u are incremented.

See also

getCutoff() and getSwitchingWidth().

Definition at line 229 of file potential_base.cpp.

                                                 {
  double d1, d[3];
  // d: vector distance (r1-r2), d1 distance and d2 squared distance
  distance(r1, r2, d, d1);
  if (d1 < cutoff_) {
    double f, en;
    coulomb(d1, f, en, qq);
    switching(d1, f, en);
    for (int l = 0; l < 3; ++l) {
      f1[l] += f * d[l] / d1;
      f2[l] -= f * d[l] / d1;
    };
    energy += en;
  }
}

distance() [1/5]

void PotentialBase::distance ( const double x[], const double y[], double & z1 )

protected

Distance.

See also

distance

Definition at line 586 of file potential_base.cpp.

                                                                           {
  double z[3];
  distance(x, y, z, z1);
}

distance() [2/5]

void PotentialBase::distance ( const double x[], const double y[], double z[] )

protected

Distance vector, norm.

See also

distance(const double x[], const double y[], double z[], double & z1, double & z2) and setPeriodicity().

Definition at line 581 of file potential_base.cpp.

                                                                           {
  double z1;
  distance(x, y, z, z1);
}

distance() [3/5]

void PotentialBase::distance ( const double x[], const double y[], double z[], double & z1 )

protected

Distance.

Parameters

[in]	x, y	Two three dimension vector.
[out]	z[]
[out]	z1	\(| \mathbf x -\mathbf y | \) with pbc applied.

See also

setPeriodicity().

Definition at line 572 of file potential_base.cpp.

                                         {
  double z2;
  distance(x, y, z, z1, z2);
}

distance() [4/5]

void PotentialBase::distance ( const double x[], const double y[], double z[], double & z1, double & z2 )

protected

Distance vector, norm and norm square.

Parameters

[in]	x, y	Two three dimension vector.
[out]	z	\( \mathbf x -\mathbf y \) with pbc applied.
[out]	z1	\( |\mathbf z | \)
[out]	z2	\( |\mathbf z |^2 \)

See also

setPeriodicity().

Definition at line 555 of file potential_base.cpp.

                                                     {
  assert(x);
  assert(y);
  assert(z);
  for (int i = 0; i < 3; ++i)
    z[i] = x[i] - y[i];
  applyPeriodicity1(z);
  z2 = z[0] * z[0] + z[1] * z[1] + z[2] * z[2];
  z1 = std::sqrt(z2);
}

distance() [5/5]

void PotentialBase::distance ( const double x[], const double y[], Vector3 & z )

protected

Distance vector, norm and norm square.

Parameters

[in]	x, y	Two three dimension vector.
[out]	z	Return x-y, |x-y| and |x-y|^2.

See also

setPeriodicity().

Definition at line 595 of file potential_base.cpp.

                                                                           {
  distance(x, y, z.v, z._1, z._2);
}

epsilon()

double PotentialBase::epsilon ( double const A, double const B )

staticprotected

Conversion for Lennard-Jones.

The Lennard-Jones is sometimes expressed as:

\[ \frac{A}{r^{12}}-\frac{B}{r^6} \]

whereas the functions lennardJones() use the expression:

\[ 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right] \]

The function calculates the parameter epsilon.

Definition at line 333 of file potential_base.cpp.

                                                            {
  return B * B / A / 4.0;
}

getCutoff()

double PotentialBase::getCutoff

(

)

const

Definition at line 60 of file potential_base.cpp.

{ return cutoff_; }

getSwitchingWidth()

double PotentialBase::getSwitchingWidth

(

)

const

Width of the switching zone.

In the switching zone the potential is replaced by \( V=U\times S(x) \) where U is the real potential V the switching potential and S the switching function. Variable x=0 at the beginning of the switching zone and x=1 at the end the function S is \( S=2x^3-3x^2+1 \) which has the following properties: S(0)=1, S(1)=0, S'(0)=0 and S'(1)=0. The switching when the distance between the two molecules is getCutoff()-getSwitchingWidth() and ends when the distance is getCutoff();

See also

getCutoff().

Definition at line 74 of file potential_base.cpp.

{ return switchingWidth_; }

isoscelesBase()

double PotentialBase::isoscelesBase ( double length, double angle )

staticprotected

Calculate the base of an isosceles triangle.

An isosceles triangle is a triangle with two sides of equal length. If you know the length of the two equal sides and the angle between the two, the function can calculate the length of the third side.

Parameters

[in]	length	Length of the two equal sides of the triangle.
[in]	angle	Angle between the two equal side in radiant.

Definition at line 116 of file potential_base.cpp.

                                                               {
  return 2.0 * length * std::sin(angle / 2.0);
}

lennardJones() [1/2]

void PotentialBase::lennardJones ( const double R1[], const double R2[], double F1[], double F2[], double & E, double const epsilon, double const sigma )

protected

Lennard-Jones 12-6 between two atoms.

Parameters

[in]	R1, R2	Positions of atom 1 and 2.
[in,out]	F1, F2	Force on atom 1 and 2.
[in,out]	E	Energy.
[in]	epsilon, sigma	Parameters of the Lennard-Jones potential.

Warning

F1, F2 and u are incremented.

See also

lennardJones(distance, ...).

Definition at line 363 of file potential_base.cpp.

                                                                           {
  // WARNING: F1 and F2 are incremented.
  double F, Fx, en;
  double R12[3];
  double L, L2;
  distance(R2, R1, R12, L, L2);
  lennardJones(L, F, en, epsilon, sigma);
  E += en;
  Fx = F * R12[0] / L;
  F1[0] -= Fx;
  F2[0] += Fx;
  Fx = F * R12[1] / L;
  F1[1] -= Fx;
  F2[1] += Fx;
  Fx = F * R12[2] / L;
  F1[2] -= Fx;
  F2[2] += Fx;
}

lennardJones() [2/2]

void PotentialBase::lennardJones ( double const distance, double & force, double & energy, double const epsilon, double const sigma )

protected

Lennard-Jones 12-6 between two atoms.

Return the energy and scalar force resulting from a Lennard-Jones potential. \( E=4*\epsilon*(x^{12}-x^6) \) with \( x=\frac{\sigma}{r} \) and \(F=-\frac{dE}{dr} \) .

Parameters

[in]	distance	Distance between two atoms.
[out]	force	Force.
[out]	energy	Potential energy.
[in]	epsilon, sigma	Lennard-Jones parameters.

See also

sigma() and epsilon()

Definition at line 347 of file potential_base.cpp.

                                                                           {
  double x = sigma * sigma / (distance * distance);
  x *= x * x;
  energy = 4 * epsilon * (x - 1) * x;
  force = 24 * epsilon * (2 * x - 1) * x / distance;
}

lennardJonesWithCutoff()

void PotentialBase::lennardJonesWithCutoff ( double const r1[], double const r2[], double f1[], double f2[], double & energy, double const epsilon, double const sigma )

protected

Lennard Jones 12-6 Potential with cutoff.

Interaction between two atoms. Cutoff and switching width should be set with setCutoff() and setSwitchingWidth().

Parameters

[in]	r1, r2	Positions of the two atoms.
[in,out]	f1, f2	Forces on the two atoms.
[in,out]	energy	Potential energy.
[in]	epsilon, sigma	See lennardJones() for definition of sigma and esplion.

See also

lennardJones().

Definition at line 397 of file potential_base.cpp.

                                                               {
  double r12[3] = {0}, d = 0.0, f, en; // distance
  distance(r2, r1, r12, d);
  if (d < cutoff_) {
    lennardJones(d, f, en, epsilon, sigma);
    switching(d, f, en);
    for (int i = 0; i < 3; ++i) {
      f1[i] -= f * r12[i] / d;
      f2[i] += f * r12[i] / d;
    };
    energy += en;
  };
}

restrainAngle()

void PotentialBase::restrainAngle ( double const r1[], double const r2[], double const r3[], double f1[], double f2[], double f3[], double & u, const double k, const double aeq )

protected

Angular quadratic restraints.

\( U=k(a-a_0)^2 \)

Parameters

[in]	r1[], r2[], r3[]	Positions of atoms.
[in,out]	f1[], f2[], f3[]	Force on atoms. The forces are incremeneted by the function.
[in,out]	u	Incremented by energy.
[in]	k	Potential curvature
[in]	aeq	Angle made by the three atoms 1,2,3 at equilibrium.

Definition at line 160 of file potential_base.cpp.

                                                    {
  double r12[3], r12_1, r12_2, r23[3], r23_1, r23_2, cosa, a, m, d12, d23;
  distance(r2, r1, r12, r12_1, r12_2);
  distance(r3, r2, r23, r23_1, r23_2);
  cosa = -dotProduct(r12, r23) / r12_1 / r23_1;
  a = std::acos(cosa);
  u += k * (a - aeq) * (a - aeq);
  m = -2 * k * (a - aeq);               // moment
  m /= -1 * std::sqrt(1 - cosa * cosa); //*d(a)/d(cos a)
  for (int i = 0; i < 3; ++i) {
    d12 = -cosa * r12[i] / r12_2 - r23[i] / r12_1 / r23_1; // d(cos a)/d(r12)
    d23 = -cosa * r23[i] / r23_2 - r12[i] / r12_1 / r23_1; // d(cos a)/d(r23)
    f1[i] -= m * d12;
    f2[i] += m * (d12 - d23);
    f3[i] += m * d23;
  };
}

restrainLength()

void PotentialBase::restrainLength ( const double R1[], const double R2[], double F1[], double F2[], double & u, const double k, const double req )

protected

Compute quadratic restraints between two atoms.

\( U=k(r-r_0)^2 \)

Parameters

[in]	R1[]	Positions of atom 1.
[in]	R2[]	Positions of atom 2.
[in,out]	F1[]	Incremented by force on atom 1.
[in,out]	F2[]	Incremented by force on atom 2.
[in,out]	u	Incremented by energy.
[in]	k	Potential curvature
[in]	req	Distance of the two atoms at equilibrium.

Warning

F1, F2 and u are incremented.

Definition at line 130 of file potential_base.cpp.

                                                                     {
  double f, fx;
  double r[3], r1, r2;
  distance(R2, R1, r, r1, r2);
  double d = r1 - req;
  f = -2.0 * k * d;
  u += k * d * d;
  fx = f * r[0] / r1;
  F1[0] -= fx;
  F2[0] += fx;
  fx = f * r[1] / r1;
  F1[1] -= fx;
  F2[1] += fx;
  fx = f * r[2] / r1;
  F1[2] -= fx;
  F2[2] += fx;
}

setCutoff()

void PotentialBase::setCutoff

(

double

cutoff

)

Definition at line 62 of file potential_base.cpp.

{ cutoff_ = cutoff; }

setPeriodicity()

void PotentialBase::setPeriodicity ( const double periods[] )

protected

Set periodicity.

Must be used before calling certain functions: distance(), applyPeriodicity(), etc...

Definition at line 603 of file potential_base.cpp.

                                                         {
  periods_[0] = periods[0];
  periods_[1] = periods[1];
  periods_[2] = periods[2];
}

setSwitchingWidth()

void PotentialBase::setSwitchingWidth

(

double

width

)

Definition at line 76 of file potential_base.cpp.

{ switchingWidth_ = width; }

sigma()

double PotentialBase::sigma ( double const A, double const B )

staticprotected

Conversion for Lennard-Jones.

The Lennard-Jones is sometimes expressed as:

\[ \frac{A}{r^{12}}-\frac{B}{r^6} \]

whereas the functions lennardJones() use the expression:

\[ 4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right] \]

The function calculates the parameter sigma.

Definition at line 423 of file potential_base.cpp.

                                                          {
  return std::pow(A / B, 1.0 / 6.0);
}

smithKongEpsilon()

double PotentialBase::smithKongEpsilon ( double sigma1, double epsilon1, double sigma2, double epsilon2 )

staticprotected

Smith and Kong combination rules.

\[ \epsilon^{SK}_{12}= \frac{2^{13}\epsilon_1\sigma_1^6\epsilon_2\sigma_2^6}{\left[(\epsilon_1\sigma_1^{12})^\frac{1}{13}+(\epsilon_1\sigma_1^{12})^\frac{1}{13}\right]^{13}} \]

See also

smithKongSigma()

Definition at line 433 of file potential_base.cpp.

                                                                       {
  double n, d;
  double s1_2 = sigma1 * sigma1, sa6 = s1_2 * s1_2 * s1_2;
  double s2_2 = sigma2 * sigma2, sb6 = s2_2 * s2_2 * s2_2;
  n = 8192.0 * epsilon1 * sa6 * epsilon2 * sb6; // 2^13 = 8192
  d = std::pow(epsilon1 * sa6 * sa6, 1.0 / 13.0) +
      std::pow(epsilon2 * sb6 * sb6, 1.0 / 13.0);
  d = std::pow(d, 13);
  return n / d;
}

smithKongSigma()

double PotentialBase::smithKongSigma ( double sigma1, double epsilon1, double sigma2, double epsilon2 )

staticprotected

Smith and Kong combination rules.

Combination of unlike Lennard-Jones parameters

\[ \sigma^{SK}= \left[ \frac{(\epsilon_1\sigma_1^{12})^{\frac{1}{13}}+(\epsilon_2\sigma_2^{12})^{\frac{1}{13}}} {2^{13}\sqrt{\epsilon_1\sigma_1^6\epsilon_2\sigma_2^6}} \right]^\frac{1}{6} \]

Definition at line 455 of file potential_base.cpp.

                                                                     {
  double s1_2 = sigma1 * sigma1;
  double const sa6 = s1_2 * s1_2 * s1_2;
  double s2_2 = sigma2 * sigma2;
  double const sb6 = s2_2 * s2_2 * s2_2;
  double n, d;
  n = std::pow(epsilon1 * sa6 * sa6, 1.0 / 13.0) +
      std::pow(epsilon2 * sb6 * sb6, 1.0 / 13.0);
  n = std::pow(n, 13.0);
  d = 8192.0 * std::sqrt(epsilon1 * sa6 * epsilon2 * sb6); // 2^13 = 8192
  return std::pow(n / d, 1.0 / 6.0);
}

spreadForce() [1/2]

void PotentialBase::spreadForce ( double f1[], double f2[], double const fc[] )

protected

Spread force on centre to other points.

Spread force fc to f1, f2. Each receives a half of the force.

Definition at line 284 of file potential_base.cpp.

                                                                           {
  for (int i = 0; i < 3; ++i) {
    double const f = fc[i] / 2;
    f1[i] += f;
    f2[i] += f;
  };
}

spreadForce() [2/2]

void PotentialBase::spreadForce ( double f1[], double f2[], double f3[], double const fc[] )

protected

Spread force on centre to other points.

Spread force fc to f1, f2, f3. Each receives 1/3 of the force.

Definition at line 295 of file potential_base.cpp.

                                                   {
  for (int i = 0; i < 3; ++i) {
    double const f = fc[i] / 3;
    f1[i] += f;
    f2[i] += f;
    f3[i] += f;
  };
}

spreadWeightedForce()

void PotentialBase::spreadWeightedForce ( double const w1, double const w2, double const w3, double f1[], double f2[], double f3[], double const fc[] )

protected

Spread force on barycentre to atoms.

This is the inverse functions calculateWeightedCentre()

Precondition

Sum of weigths w1, w2, w3 must be 1.

Parameters

[in]	w1, w2, w3	weights.
[in,out]	f1, f2, f3	vector to spread fc to.
[in,out]	fc	force to spread.

Definition at line 312 of file potential_base.cpp.

                                                           {
  for (int i = 0; i < 3; ++i) {
    double const f = fc[i];
    f1[i] += f * w1;
    f2[i] += f * w2;
    f3[i] += f * w3;
  };
}

switching() [1/3]

template<int N, class R, class F>

void forcefields::PotentialBase::switching ( ChargeGroup< N, R, F > & g1, ChargeGroup< N, R, F > & g2, double & energy, double cutoff, double switchingWidth )

protected

Definition at line 163 of file potential_base.hpp.

                                                                  {
  double z[3], z1;
  distance(g1.centre_, g2.centre_, z, z1);
  assert(z1 <= cutoff);
  assert(z1 >= cutoff - switchingWidth);
  double x = (z1 - cutoff + switchingWidth) / switchingWidth;
  double S = (2 * x - 3) * x * x + 1;
  double dS = 6 * x * (x - 1);
  double const u = energy;
  for (int i = 0; i < N; ++i) {
    for (int k = 0; k < 3; ++k) {
      g1.force_[i][k] =
          g1.force_[i][k] * S - u * dS / switchingWidth * z[k] / z1 / N;
      g2.force_[i][k] =
          g2.force_[i][k] * S + u * dS / switchingWidth * z[k] / z1 / N;
    }
  };
  energy *= S;
}

switching() [2/3]

void PotentialBase::switching ( double const distance, double & force, double & energy )

protected

Smooth cutoff switch off.

Should be called when the distance between the two atoms is in the switching zone (see getSwitchinWidth() for more explanation).

Parameters

[in]	distance	Distance between two atoms.
[in,out]	force	On atoms. The input values should be the full force for the interaction between the two atoms. The function returns the corrected force.
[in,out]	energy	Potential energy. The input should be the full energy for the interaction between the two atoms. The functions returns the corrected energy.

Note

The function is for atom based cutoff only.

See also

getSwitchinWidth().

Definition at line 481 of file potential_base.cpp.

                                              {
  double x = (distance - cutoff_ + switchingWidth_) / switchingWidth_;
  if (x >= 1.0) { // if more: cutoff
    force = 0.0;
    energy = 0.0;
  }
  if (x > 0.0) { // if more: switching zone
    double S = (2 * x - 3) * x * x + 1;
    double dS = 6 * x * (x - 1);
    double const u = energy;
    force = force * S - u * dS / switchingWidth_;
    energy *= S;
  };
  // else: full interaction, leave force and energy as the are.
}

switching() [3/3]

void PotentialBase::switching ( double const r1[], double const r2[], double f1[], double f2[], double & energy )

protected

Smooth cutoff switch off.

Should be called when the distance between the two atoms is in the switching zone (see getSwitchinWidth() for more explanation).

Parameters

[in]	r1, r2	Positions of atoms.
[in,out]	f1, f2	Forces on atoms. The input values should be the full forces for the interaction between the two atoms. The function returns the corrected forces.
[in,out]	energy	Potential energy. The input should be the full energy for the interaction between the two atoms. The functions returns the corrected energy.

Note

The function is for atom based cutoff only.

See also

getSwitchinWidth().

Definition at line 511 of file potential_base.cpp.

                                                           {
  double z[3], z1;
  distance(r1, r2, z, z1);
  double x = (z1 - cutoff_ + switchingWidth_) / switchingWidth_;
  assert(x > 0.0);
  assert(x < 1.0);
  double S = (2 * x - 3) * x * x + 1;
  double dS = 6 * x * (x - 1);
  double const u = energy;
  for (int k = 0; k < 3; ++k) {
    f1[k] = f1[k] * S - u * dS / switchingWidth_ * z[k] / z1;
    f2[k] = f2[k] * S + u * dS / switchingWidth_ * z[k] / z1;
  };
  energy *= S;
}

unBreak0() [1/2]

double PotentialBase::unBreak0 ( double const r, double const ref, double const period )

staticprotected

Undo the separation of two atoms created by the periodic boundaries.

Example: two atoms H and O are at respectively x = -9.5 and x = 9.5 and the period along x is 20. The real distance between the two atoms is 1. We use unBreak() on the coordinates of these two atoms:

unBreak(rh, ro);

The coordinates rh of H is now 10.5 thus making O and H appears at a distance of 1 of each other.
Apply periodic boundaries condition to vector \(r - r_{ef} \).

Parameters

[in]	r	Returns \( \textrm{applyPeriodicity}(r - r_{ef}) +r_{ef} \)
[in]	ref
[in]	period

Definition at line 623 of file potential_base.cpp.

                                                    {
  return applyPeriodicity0(r - ref, period) + ref;
}

unBreak0() [2/2]

void PotentialBase::unBreak0 ( double r[], double const ref[], double const periods[] )

staticprotected

Undo the separation of two atoms created by the periodic boundaries.

See also

unBreak0()

Definition at line 629 of file potential_base.cpp.

                                                     {
  for (int a = 0; a < 3; ++a)
    r[a] = unBreak0(r[a], ref[a], periods[a]);
}

unBreak1() [1/2]

double PotentialBase::unBreak1 ( double const r, double const ref, int const axis )

protected

Undo the separation of two atoms created by the periodic boundaries.

See also

unBreak0()

Definition at line 641 of file potential_base.cpp.

                                               {
  return unBreak0(r, ref, periods_[axis]);
}

unBreak1() [2/2]

void PotentialBase::unBreak1 ( double r[], double const ref[] )

protected

Undo the separation of two atoms created by the periodic boundaries.

See also

unBreak0()

Definition at line 636 of file potential_base.cpp.

                                                           {
  unBreak0(r, ref, periods_);
}

Member Data Documentation

cutoff_

double forcefields::PotentialBase::cutoff_

protected

Definition at line 147 of file potential_base.hpp.

EPSILON_PT

const double PotentialBase::EPSILON_PT = 0.68165797577788501

staticprotected

Platinum Lennard-Jones.

Lennard-Jones parameters for platinum.

See also

computePt().

Definition at line 132 of file potential_base.hpp.

ONE_OVER_4_PI_EPSILON0

const double PotentialBase::ONE_OVER_4_PI_EPSILON0

staticprotected

Initial value:

                                               =
14.399644532010862

Definition at line 54 of file potential_base.hpp.

periods_

double forcefields::PotentialBase::periods_[3]

protected

Definition at line 148 of file potential_base.hpp.

SIGMA_PT

const double PotentialBase::SIGMA_PT = 2.54

staticprotected

Definition at line 133 of file potential_base.hpp.

switchingWidth_

double forcefields::PotentialBase::switchingWidth_

protected

Definition at line 149 of file potential_base.hpp.

---

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

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