#include <array>
+#include "nblib/exception.h"
#include "nblib/particletype.h"
#include "nblib/util/util.hpp"
/*! \brief FENE bond type
*
* It represents the interaction of the form
- * V(r, forceConstant, equilConstant) = - 0.5 * forceConstant * equilDistance^2 * log( 1 - (r / equilConstant)^2)
+ * V(r, forceConstant, equilDistance) = - 0.5 * forceConstant * equilDistance^2 * log( 1 - (r / equilConstant)^2)
*/
using FENEBondType = TwoParameterInteraction<struct FENEBondTypeParameter>;
/*! \brief Half-attractive quartic bond type
*
* It represents the interaction of the form
- * V(r, forceConstant, equilConstant) = 0.5 * forceConstant * (r - equilConstant)^4
+ * V(r, forceConstant, equilConstant) = 0.5 * forceConstant * (r - equilDistance)^4
*/
using HalfAttractiveQuarticBondType =
TwoParameterInteraction<struct HalfAttractiveQuarticBondTypeParameter>;
*
* It represents the interaction of the form
* V(r, quadraticForceConstant, cubicForceConstant, equilConstant) = quadraticForceConstant * (r -
- * equilConstant)^2 + quadraticForceConstant * cubicForceConstant * (r - equilDistance)
+ * equilDistance)^2 + quadraticForceConstant * cubicForceConstant * (r - equilConstant)
*/
-struct CubicBondType
+class CubicBondType
{
+public:
CubicBondType() = default;
CubicBondType(ForceConstant fq, ForceConstant fc, EquilConstant d) :
quadraticForceConstant_(fq), cubicForceConstant_(fc), equilDistance_(d)
== std::tie(b.forceConstant(), b.exponent(), b.equilDistance());
}
+/*! \brief Non-Bonded Pair Interaction Type
+ *
+ * It represents the interaction of the form
+ * of LJ interactions, but occur between atoms
+ * of the same bonded chain
+ * V(r, c_6, c_12) = c_12*(r^-12) - c_6*(r^-6)
+ */
+class PairLJType
+{
+public:
+ PairLJType() = default;
+ PairLJType(C6 c6, C12 c12) : c6_(c6), c12_(c12) {}
+
+ [[nodiscard]] const C6& c6() const { return c6_; }
+ [[nodiscard]] const C12& c12() const { return c12_; }
+
+private:
+ C6 c6_;
+ C12 c12_;
+};
+
+inline bool operator<(const PairLJType& a, const PairLJType& b)
+{
+ return std::tie(a.c6(), a.c12()) < std::tie(b.c6(), b.c12());
+}
+
+inline bool operator==(const PairLJType& a, const PairLJType& b)
+{
+ return std::tie(a.c6(), a.c12()) == std::tie(b.c6(), b.c12());
+}
/*! \brief Basic template for interactions with 2 parameters named forceConstant and equilAngle
*
*/
using HarmonicAngle = AngleInteractionType<struct HarmonicAngleParameter>;
+/*! \brief linear angle type
+ *
+ * It represents the interaction of the form
+ * V(theta, forceConstant, a) = 0.5 * forceConstant * (dr)^2
+ * where dr = - a * r_ij - (1 - a) * r_kj
+ */
+using LinearAngle = TwoParameterInteraction<struct LinearAngleParameter>;
+
+/*! \brief Basic template for angle types that use the cosines of their equilibrium angles
+ * in their potential expression
+ */
+template<class Phantom>
+class CosineParamAngle : public TwoParameterInteraction<Phantom>
+{
+public:
+ CosineParamAngle() = default;
+ //! \brief construct from angle given in radians
+ CosineParamAngle(ForceConstant f, Radians angle) :
+ TwoParameterInteraction<Phantom>{ f, std::cos(angle) }
+ {
+ }
+
+ //! \brief construct from angle given in degrees
+ CosineParamAngle(ForceConstant f, Degrees angle) :
+ TwoParameterInteraction<Phantom>{ f, std::cos(angle * DEG2RAD) }
+ {
+ }
+};
+
+/*! \brief G96 or Cosine-based angle type
+ *
+ * This represents the interaction of the form
+ * V(cos(theta), forceConstant, cos(equilAngle)) = 0.5 * forceConstant * (cos(theta) - cos(equilAngle))^2
+ */
+using G96Angle = CosineParamAngle<struct G96AngleParameter>;
+
+/*! \brief Restricted angle type
+ *
+ * This represents the interaction of the form
+ * V(cos(theta), forceConstant, cos(equilAngle)) =
+ * 0.5 * forceConstant * (cos(theta) - cos(equilAngle))^2 / (sin(theta))^2
+ */
+using RestrictedAngle = CosineParamAngle<struct RestrictedAngleParameter>;
+
+/*! \brief Quartic angle type
+ *
+ * It represents the interaction of the form of a fourth order polynomial
+ * V(theta, forceConstant, equilAngle) = sum[i = 0 -> 4](forceConstant_i * (theta - equilAngle)^i
+ */
+class QuarticAngle
+{
+public:
+ QuarticAngle() = default;
+ //! \brief construct from given angle in radians
+ QuarticAngle(ForceConstant f0, ForceConstant f1, ForceConstant f2, ForceConstant f3, ForceConstant f4, Radians angle) :
+ forceConstants_{ f0, f1, f2, f3, f4 }, equilConstant_(angle)
+ {
+ }
+
+ //! \brief construct from given angle in degrees
+ QuarticAngle(ForceConstant f0, ForceConstant f1, ForceConstant f2, ForceConstant f3, ForceConstant f4, Degrees angle) :
+ forceConstants_{ f0, f1, f2, f3, f4 }, equilConstant_(angle * DEG2RAD)
+ {
+ }
+
+ [[nodiscard]] const std::array<ForceConstant, 5>& forceConstants() const
+ {
+ return forceConstants_;
+ }
+
+ ForceConstant forceConstant(int order) const
+ {
+ switch (order)
+ {
+ case 0: return forceConstants_[0];
+ case 1: return forceConstants_[1];
+ case 2: return forceConstants_[2];
+ case 3: return forceConstants_[3];
+ case 4: return forceConstants_[4];
+ default:
+ throw InputException(
+ "Please enter a value between 0-4 for the Quartic Angle force constants");
+ }
+ }
+
+ Radians equilConstant() const { return equilConstant_; }
+
+private:
+ std::array<ForceConstant, 5> forceConstants_;
+ Radians equilConstant_;
+};
+
+inline bool operator<(const QuarticAngle& a, const QuarticAngle& b)
+{
+ return (a.forceConstants() < b.forceConstants()) && (a.equilConstant() < b.equilConstant());
+}
+
+inline bool operator==(const QuarticAngle& a, const QuarticAngle& b)
+{
+ return (a.forceConstants() == b.forceConstants()) && (a.equilConstant() == b.equilConstant());
+}
+
+
+/*! \brief Cross bond-bond interaction type
+ */
+class CrossBondBond
+{
+public:
+ CrossBondBond() = default;
+ CrossBondBond(ForceConstant f, EquilConstant r0ij, EquilConstant r0kj) :
+ forceConstant_(f), r0ij_(r0ij), r0kj_(r0kj)
+ {
+ }
+
+ [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
+ [[nodiscard]] const EquilConstant& equilDistanceIJ() const { return r0ij_; }
+ [[nodiscard]] const EquilConstant& equilDistanceKJ() const { return r0kj_; }
+
+private:
+ ForceConstant forceConstant_;
+ EquilConstant r0ij_;
+ EquilConstant r0kj_;
+};
+
+inline bool operator<(const CrossBondBond& a, const CrossBondBond& b)
+{
+ return std::tie(a.forceConstant(), a.equilDistanceIJ(), a.equilDistanceKJ())
+ < std::tie(b.forceConstant(), b.equilDistanceIJ(), b.equilDistanceKJ());
+}
+
+inline bool operator==(const CrossBondBond& a, const CrossBondBond& b)
+{
+ return std::tie(a.forceConstant(), a.equilDistanceIJ(), a.equilDistanceKJ())
+ == std::tie(b.forceConstant(), b.equilDistanceIJ(), b.equilDistanceKJ());
+}
+
+/*! \brief Cross bond-angle interaction type
+ */
+class CrossBondAngle
+{
+public:
+ CrossBondAngle() = default;
+ CrossBondAngle(ForceConstant f, EquilConstant r0ij, EquilConstant r0kj, EquilConstant r0ik) :
+ forceConstant_(f), r0ij_(r0ij), r0kj_(r0kj), r0ik_(r0ik)
+ {
+ }
+
+ [[nodiscard]] const ForceConstant& forceConstant() const { return forceConstant_; }
+ [[nodiscard]] const EquilConstant& equilDistanceIJ() const { return r0ij_; }
+ [[nodiscard]] const EquilConstant& equilDistanceKJ() const { return r0kj_; }
+ [[nodiscard]] const EquilConstant& equilDistanceIK() const { return r0ik_; }
+
+private:
+ ForceConstant forceConstant_;
+ EquilConstant r0ij_;
+ EquilConstant r0kj_;
+ EquilConstant r0ik_;
+};
+
+inline bool operator<(const CrossBondAngle& a, const CrossBondAngle& b)
+{
+ return std::tie(a.forceConstant(), a.equilDistanceIJ(), a.equilDistanceKJ(), a.equilDistanceIK())
+ < std::tie(b.forceConstant(), b.equilDistanceIJ(), b.equilDistanceKJ(), b.equilDistanceIK());
+}
+
+inline bool operator==(const CrossBondAngle& a, const CrossBondAngle& b)
+{
+ return std::tie(a.forceConstant(), a.equilDistanceIJ(), a.equilDistanceKJ(), a.equilDistanceIK())
+ == std::tie(b.forceConstant(), b.equilDistanceIJ(), b.equilDistanceKJ(), b.equilDistanceIK());
+}
+
/*! \brief Proper Dihedral Implementation
*/
class ProperDihedral
/*! \brief Improper Dihedral Implementation
*/
-struct ImproperDihedral : public TwoParameterInteraction<struct ImproperDihdedralParameter>
+class ImproperDihedral : public TwoParameterInteraction<struct ImproperDihdedralParameter>
{
+public:
ImproperDihedral() = default;
ImproperDihedral(Radians phi, ForceConstant f) :
TwoParameterInteraction<struct ImproperDihdedralParameter>{ f, phi }
return energy;
}
+template<class BasicVector, class ShiftForce>
+inline NBLIB_ALWAYS_INLINE
+auto computeThreeCenter(const LinearAngle& parameters, const BasicVector& rij, const BasicVector& rkj,
+ const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
+ ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
+{
+ using ValueType = BasicVectorValueType_t<BasicVector>;
+
+ ValueType b = parameters.equilConstant() - 1;
+ auto dr_vec = b * rkj - parameters.equilConstant() * rij;
+ ValueType dr = norm(dr_vec);
+
+ auto [ci, ck, energy] = threeCenterKernel(dr, parameters);
+
+ BasicVector fi_loc = ci * dr_vec;
+ BasicVector fk_loc = ck * dr_vec;
+ BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
+ *fi += fi_loc;
+ *fj += fj_loc;
+ *fk += fk_loc;
+
+ addShiftForce(fi_loc, shf_ij);
+ addShiftForce(fj_loc, shf_c);
+ addShiftForce(fk_loc, shf_kj);
+
+ return energy;
+}
+
+template<class BasicVector, class ShiftForce>
+inline NBLIB_ALWAYS_INLINE
+auto computeThreeCenter(const CrossBondBond& parameters, const BasicVector& rij, const BasicVector& rkj,
+ const BasicVector& /* rik */, BasicVector* fi, BasicVector* fj, BasicVector* fk,
+ ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
+{
+ using ValueType = BasicVectorValueType_t<BasicVector>;
+ // 28 flops from the norm() calls
+ auto [ci, ck, energy] = threeCenterKernel(norm(rij), norm(rkj), parameters);
+
+ BasicVector fi_loc = ci * rij;
+ BasicVector fk_loc = ck * rkj;
+ BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
+ *fi += fi_loc;
+ *fj += fj_loc;
+ *fk += fk_loc;
+
+ addShiftForce(fi_loc, shf_ij);
+ addShiftForce(fj_loc, shf_c);
+ addShiftForce(fk_loc, shf_kj);
+
+ return energy;
+}
+
+template<class BasicVector, class ShiftForce>
+inline NBLIB_ALWAYS_INLINE
+auto computeThreeCenter(const CrossBondAngle& parameters, const BasicVector& rij, const BasicVector& rkj,
+ const BasicVector& rik, BasicVector* fi, BasicVector* fj, BasicVector* fk,
+ ShiftForce* shf_ij, ShiftForce* shf_kj, ShiftForce* shf_c)
+{
+ using ValueType = BasicVectorValueType_t<BasicVector>;
+ // 42 flops from the norm() calls
+ auto [ci, cj, ck, energy] = crossBondAngleKernel(norm(rij), norm(rkj), norm(rik), parameters);
+
+ BasicVector fi_loc = ci * rij + ck * rik;
+ BasicVector fk_loc = cj * rkj - ck * rik;
+ BasicVector fj_loc = ValueType(-1.0) * (fi_loc + fk_loc);
+ *fi += fi_loc;
+ *fj += fj_loc;
+ *fk += fk_loc;
+
+ addShiftForce(fi_loc, shf_ij);
+ addShiftForce(fj_loc, shf_c);
+ addShiftForce(fk_loc, shf_kj);
+
+ return energy;
+}
+
/*! \brief Calculate three-center interactions
*
* \tparam Force buffer type
*/
using SupportedTwoCenterTypes =
- TypeList<HarmonicBondType, G96BondType, CubicBondType, FENEBondType, HalfAttractiveQuarticBondType>;
-using SupportedThreeCenterTypes = TypeList<HarmonicAngle>;
+ TypeList<HarmonicBondType, G96BondType, CubicBondType, MorseBondType, FENEBondType, HalfAttractiveQuarticBondType, PairLJType>;
+using SupportedThreeCenterTypes =
+ TypeList<HarmonicAngle, G96Angle, QuarticAngle, RestrictedAngle, CrossBondBond, CrossBondAngle, LinearAngle>;
using SupportedFourCenterTypes = TypeList<ProperDihedral, ImproperDihedral, RyckaertBellemanDihedral>;
using SupportedFiveCenterTypes = TypeList<Default5Center>;
{
using type = InteractionType;
- // tuple format: <particleID i, particleID j, ..., InteractionInstanceIndex>
- std::vector<InteractionIndex<InteractionType>> indices;
// vector of unique TwoCenterType instances
std::vector<InteractionType> parameters;
+ // tuple format: <particleID i, particleID j, ..., InteractionInstanceIndex>
+ std::vector<InteractionIndex<InteractionType>> indices;
};
using TwoCenterInteraction = Reduce<std::variant, SupportedTwoCenterTypes>;
return morseScalarForce(bond.forceConstant(), bond.exponent(), bond.equilDistance(), dr);
}
+/*! \brief kernel to calculate the scalar part for the 1-4 LJ non-bonded forces
+ *
+ * \param c6 C6 parameter of LJ potential
+ * \param c12 C12 parameter of LJ potential
+ * \param r distance between the atoms
+ *
+ * \return tuple<force, potential energy>
+ */
+template <class T>
+inline std::tuple<T, T> pairLJScalarForce(C6 c6, C12 c12, T r)
+{
+ T rinv = 1./r; /* 1 */
+ T rinv2 = rinv * rinv; /* 1 */
+ T rinv6 = rinv2 * rinv2 * rinv2; /* 2 */
+
+ T epot = rinv6*(c12*rinv6 - c6); /* 3 */
+
+ T c6_ = 6.*c6; /* 1 */
+ T c12_ = 12.*c12; /* 1 */
+
+ T force = rinv6*(c12_*rinv6 - c6_)*rinv; /* 4 */
+
+ return std::make_tuple(force, epot);
+
+ /* That was 13 flops */
+}
+
+//! Abstraction layer for different 2-center bonds. 1-4 LJ pair interactions case
+template <class T>
+inline auto bondKernel(T dr, const PairLJType& bond)
+{
+ return pairLJScalarForce(bond.c6(), bond.c12(), dr);
+}
+
/*! \brief kernel to calculate the scalar part for the FENE pontential bond forces
* for lambda = 0
//! Three-center interaction type kernels
+/*! \brief kernel to calculate the scalar part for linear angle forces
+ * for lambda = 0
+ *
+ * \param k force constant
+ * \param a contribution of rij vector
+ * \param dr weighted distance metric
+ *
+ * \return tuple<force, potential energy>
+ */
+template <class T>
+inline std::tuple<T, T, T> linearAnglesScalar(T k, T a, T dr)
+{
+ T b = T(1.0) - a;
+
+ T kdr = k * dr;
+ T epot = 0.5 * kdr * dr;
+
+ T ci = a * k;
+ T ck = b * k;
+
+ return std::make_tuple(ci, ck, epot);
+
+ /* That was 5 flops */
+}
+
+template <class T>
+inline auto threeCenterKernel(T dr, const LinearAngle& angle)
+{
+ return linearAnglesScalar(angle.forceConstant(), angle.equilConstant(), dr);
+}
//! Harmonic Angle
template <class T>
return harmonicScalarForce(angle.forceConstant(), angle.equilConstant(), dr);
}
+//! Cosine based (GROMOS-96) Angle
+template <class T>
+inline auto threeCenterKernel(T dr, const G96Angle& angle)
+{
+ auto costheta = std::cos(dr);
+ auto feTuple = g96ScalarForce(angle.forceConstant(), angle.equilConstant(), costheta);
+
+ // The above kernel call effectively computes the derivative of the potential with respect to
+ // cos(theta). However, we need the derivative with respect to theta. We use this extra -sin(theta)
+ // factor to account for this before the forces are spread between the particles.
+
+ std::get<0>(feTuple) *= -std::sqrt(1 - costheta*costheta);
+ return feTuple;
+}
+
+/*! \brief kernel to calculate the scalar part for cross bond-bond forces
+ * for lambda = 0
+ *
+ * \param k force constant
+ * \param r0ij equilibrium distance between particles i & j
+ * \param r0kj equilibrium distance between particles k & j
+ * \param rij input bond length between particles i & j
+ * \param rkj input bond length between particles k & j
+ *
+ * \return tuple<force scalar i, force scalar k, potential energy>
+ */
+
+template <class T>
+inline std::tuple<T, T, T> crossBondBondScalarForce(T k, T r0ij, T r0kj, T rij, T rkj)
+{
+ T si = rij - r0ij;
+ T sk = rkj - r0kj;
+
+ T epot = k * si * sk;
+
+ T ci = -k * sk / rij;
+ T ck = -k * si / rkj;
+
+ return std::make_tuple(ci, ck, epot);
+
+ /* That was 8 flops */
+}
+
+//! Cross bond-bond interaction
+template <class T>
+inline auto threeCenterKernel(T drij, T drkj, const CrossBondBond& crossBondBond)
+{
+ return crossBondBondScalarForce(crossBondBond.forceConstant(), crossBondBond.equilDistanceIJ(), crossBondBond.equilDistanceKJ(), drij, drkj);
+}
+
+/*! \brief kernel to calculate the scalar part for cross bond-angle forces
+ * for lambda = 0
+ *
+ * \param k force constant
+ * \param r0ij equilibrium distance between particles i & j
+ * \param r0kj equilibrium distance between particles k & j
+ * \param r0ik equilibrium distance between particles i & k
+ * \param rij input bond length between particles i & j
+ * \param rkj input bond length between particles k & j
+ * \param rik input bond length between particles i & k
+ *
+ * \return tuple<force, potential energy>
+ */
+
+template <class T>
+inline std::tuple<T, T, T, T> crossBondAngleScalarForce(T k, T r0ij, T r0kj, T r0ik, T rij, T rkj, T rik)
+{
+ T sij = rij - r0ij;
+ T skj = rkj - r0kj;
+ T sik = rik - r0ik;
+
+ T epot = k * sik * (sij + skj);
+
+ T ci = -k * sik / rij;
+ T cj = -k * sik / rkj;
+ T ck = -k * (sij + skj) / rik;
+
+ return std::make_tuple(ci, cj, ck, epot);
+
+ /* That was 13 flops */
+}
+
+//! Cross bond-bond interaction
+template <class T>
+inline auto crossBondAngleKernel(T drij, T drkj, T drik, const CrossBondAngle& crossBondAngle)
+{
+ return crossBondAngleScalarForce(crossBondAngle.forceConstant(), crossBondAngle.equilDistanceIJ(), crossBondAngle.equilDistanceKJ(), crossBondAngle.equilDistanceIK(), drij, drkj, drik);
+}
+
+//! Quartic Angle
+template <class T>
+inline auto threeCenterKernel(T dr, const QuarticAngle& angle)
+{
+ T dt = dr - angle.equilConstant(); /* 1 */
+
+ T force = 0;
+ T energy = angle.forceConstant(0);
+ T dtp = 1.0;
+ for (auto j = 1; j <= 4; j++)
+ { /* 24 */
+ T c = angle.forceConstant(j);
+ force -= j * c * dtp; /* 3 */
+ dtp *= dt; /* 1 */
+ energy += c * dtp; /* 2 */
+ }
+
+ /* TOTAL 25 */
+ return std::make_tuple(force, energy);
+}
+
+//! \brief Restricted Angle potential. Returns scalar force and energy
+template <class T>
+inline auto threeCenterKernel(T theta, const RestrictedAngle& angle)
+{
+ T costheta = std::cos(theta);
+ auto [force, ePot] = harmonicScalarForce(angle.forceConstant(), angle.equilConstant(), costheta);
+
+ // The above kernel call effectively computes the derivative of the potential with respect to
+ // cos(theta). However, we need the derivative with respect to theta.
+ // This introduces the extra (cos(theta)*cos(eqAngle) - 1)/(sin(theta)^3 factor for the force
+ // The call also computes the potential energy without the sin(theta)^-2 factor
+
+ T sintheta2 = (1 - costheta*costheta);
+ T sintheta = std::sqrt(sintheta2);
+ force *= (costheta*angle.equilConstant() - 1)/(sintheta2*sintheta);
+ ePot /= sintheta2;
+ return std::make_tuple(force, ePot);
+}
//! \brief Computes and returns the proper dihedral force
template <class T>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">2.5625292086196829</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">3</Int>
+ <Vector>
+ <Real Name="X">20.584931942381417</Real>
+ <Real Name="Y">32.355853840934138</Real>
+ <Real Name="Z">-40.148174152833434</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-1.5714394048923204</Real>
+ <Real Name="Y">1.7263018958837222</Real>
+ <Real Name="Z">0.65914251711356542</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-19.013492537489096</Real>
+ <Real Name="Y">-34.08215573681786</Real>
+ <Real Name="Z">39.489031635719869</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">15.930047836493104</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">3</Int>
+ <Vector>
+ <Real Name="X">16.257587884629597</Real>
+ <Real Name="Y">2.253527033513016</Real>
+ <Real Name="Z">-18.35014870146313</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-18.352335074096068</Real>
+ <Real Name="Y">23.83559705347821</Real>
+ <Real Name="Z">5.5912340019856792</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">2.0947471894664718</Real>
+ <Real Name="Y">-26.089124086991227</Real>
+ <Real Name="Z">12.758914699477451</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">0.878581203916453</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">3</Int>
+ <Vector>
+ <Real Name="X">19.022908985367618</Real>
+ <Real Name="Y">-52.884817112836259</Real>
+ <Real Name="Z">10.359003227565086</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-61.096971498538466</Real>
+ <Real Name="Y">67.108390759454124</Real>
+ <Real Name="Z">25.632702402756305</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">42.074062513170844</Real>
+ <Real Name="Y">-14.223573646617863</Real>
+ <Real Name="Z">-35.991705630321391</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">0.20264300000000018</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">3</Int>
+ <Vector>
+ <Real Name="X">-0.94000000000000127</Real>
+ <Real Name="Y">1.532</Real>
+ <Real Name="Z">0.10800000000000157</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">2.3500000000000032</Real>
+ <Real Name="Y">-3.8300000000000001</Real>
+ <Real Name="Z">-0.2700000000000039</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-1.4100000000000019</Real>
+ <Real Name="Y">2.298</Real>
+ <Real Name="Z">0.16200000000000236</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">3.7653963456510402</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">3</Int>
+ <Vector>
+ <Real Name="X">28.524143507898081</Real>
+ <Real Name="Y">-79.298813545069052</Real>
+ <Real Name="Z">15.532939514620491</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-91.612633181538157</Real>
+ <Real Name="Y">100.62653246562743</Real>
+ <Real Name="Z">38.435282552939967</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">63.088489673640076</Real>
+ <Real Name="Y">-21.327718920558382</Real>
+ <Real Name="Z">-53.968222067560454</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">1.0103314126941512</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">3</Int>
+ <Vector>
+ <Real Name="X">23.578508969297374</Real>
+ <Real Name="Y">-65.549655712163073</Real>
+ <Real Name="Z">12.839773911655692</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">-75.728454127772068</Real>
+ <Real Name="Y">83.179486095108572</Real>
+ <Real Name="Z">31.771213539191535</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">52.149945158474708</Real>
+ <Real Name="Y">-17.629830382945496</Real>
+ <Real Name="Z">-44.610987450847226</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Real Name="Epot">0.0018852207068626682</Real>
+ <Sequence Name="forces">
+ <Int Name="Length">2</Int>
+ <Vector>
+ <Real Name="X">-0.84131268781639024</Real>
+ <Real Name="Y">-0.11661760029138095</Real>
+ <Real Name="Z">0.94960045951553052</Real>
+ </Vector>
+ <Vector>
+ <Real Name="X">0.84131268781639024</Real>
+ <Real Name="Y">0.11661760029138095</Real>
+ <Real Name="Z">-0.94960045951553052</Real>
+ </Vector>
+ </Sequence>
+</ReferenceData>
//! Function types for testing cubic bonds
static const std::vector<CubicBondType> c_InputCubicBonds = { { CubicBondType(50.0, 2.0, 0.16) } };
+//! Function types for testing Morse bonds
+static const std::vector<MorseBondType> c_InputMorseBonds = { { MorseBondType(30.0, 2.7, 0.15) } };
+
//! Function types for testing FENE bonds
static const std::vector<FENEBondType> c_InputFeneBonds = { { FENEBondType(5.0, 0.4) } };
//! Function types for testing Harmonic angles
static const std::vector<HarmonicAngle> c_InputHarmonicAngles = { { HarmonicAngle(50.0, Degrees(100)) } };
+//! Function types for testing Linear angles
+static const std::vector<LinearAngle> c_InputLinearAngles = { { LinearAngle(50.0, 0.4) } };
+
+//! Function types for testing G96 angles
+static const std::vector<G96Angle> c_InputG96Angles = { { G96Angle(50.0, Degrees(100)) } };
+
+//! Function types for testing Restricted angles
+static const std::vector<RestrictedAngle> c_InputRestrictedAngles = { { RestrictedAngle(50.0, Degrees(100)) } };
+
+//! Function types for testing Quartic angles
+static const std::vector<QuarticAngle> c_InputQuarticAngles = {
+ { QuarticAngle(1.1, 2.3, 4.6, 7.8, 9.2, Degrees(87)) }
+};
+
+//! Function types for testing cross bond-bond interaction
+static const std::vector<CrossBondBond> c_InputCrossBondBond = { { CrossBondBond(45.0, 0.8, 0.7) } };
+
+//! Function types for testing cross bond-angle interaction
+static const std::vector<CrossBondAngle> c_InputCrossBondAngle = { { CrossBondAngle(45.0, 0.8, 0.7, 0.3) } };
+
//! Function types for testing dihedrals
static const std::vector<ProperDihedral> c_InputDihs = { { ProperDihedral(Degrees(-105.0), 15.0, 2) } };
checkForcesAndEnergiesWithRefData(c_InputDihs, c_coordinatesForDihTests);
}
+TEST(ThreeCenter, ListedForcesG96AngleTest)
+{
+ checkForcesAndEnergiesWithRefData(c_InputG96Angles, c_coordinatesForAngleTests);
+}
+
TEST(ThreeCenter, ListedForcesHarmonicAngleTest)
{
checkForcesAndEnergiesWithRefData(c_InputHarmonicAngles, c_coordinatesForAngleTests);
}
+TEST(ThreeCenter, ListedForcesLinearAngleTest)
+{
+ checkForcesAndEnergiesWithRefData(c_InputLinearAngles, c_coordinatesForAngleTests);
+}
+
+TEST(ThreeCenter, ListedForcesCrossBondBondTest)
+{
+ checkForcesAndEnergiesWithRefData(c_InputCrossBondBond, c_coordinatesForAngleTests);
+}
+
+TEST(ThreeCenter, ListedForcesCrossBondAngleTest)
+{
+ checkForcesAndEnergiesWithRefData(c_InputCrossBondAngle, c_coordinatesForAngleTests);
+}
+
+TEST(ThreeCenter, ListedForcesQuarticAngleTest)
+{
+ checkForcesAndEnergiesWithRefData(c_InputQuarticAngles, c_coordinatesForAngleTests);
+}
+
+TEST(ThreeCenter, ListedForcesRestrictedAngleTest)
+{
+ checkForcesAndEnergiesWithRefData(c_InputRestrictedAngles, c_coordinatesForAngleTests);
+}
+
TEST(TwoCenter, ListedForcesHarmonicBondTest)
{
checkForcesAndEnergiesWithRefData(c_InputHarmonicBonds, c_coordinatesForBondTests);
checkForcesAndEnergiesWithRefData(c_InputCubicBonds, c_coordinatesForBondTests);
}
+TEST(TwoCenter, ListedForcesMorseBondTest)
+{
+ checkForcesAndEnergiesWithRefData(c_InputMorseBonds, c_coordinatesForBondTests);
+}
+
TEST(TwoCenter, ListedForcesFeneBondTest)
{
checkForcesAndEnergiesWithRefData(c_InputFeneBonds, c_coordinatesForBondTests);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff