From: Joe Jordan Date: Thu, 26 Aug 2021 07:14:30 +0000 (+0000) Subject: Add several more listed forces interaction types to nblib X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=c965670dea5627acc23ff89f0690d794c71f9cb7;p=alexxy%2Fgromacs.git Add several more listed forces interaction types to nblib --- diff --git a/api/nblib/listed_forces/bondtypes.h b/api/nblib/listed_forces/bondtypes.h index 51baf103d6..5a69f787e5 100644 --- a/api/nblib/listed_forces/bondtypes.h +++ b/api/nblib/listed_forces/bondtypes.h @@ -54,6 +54,7 @@ #include +#include "nblib/exception.h" #include "nblib/particletype.h" #include "nblib/util/util.hpp" @@ -135,7 +136,7 @@ public: /*! \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; @@ -143,7 +144,7 @@ using FENEBondType = TwoParameterInteraction; /*! \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; @@ -153,10 +154,11 @@ using HalfAttractiveQuarticBondType = * * 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) @@ -224,6 +226,36 @@ inline bool operator==(const MorseBondType& a, const MorseBondType& b) == 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 * @@ -261,6 +293,177 @@ public: */ using HarmonicAngle = AngleInteractionType; +/*! \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; + +/*! \brief Basic template for angle types that use the cosines of their equilibrium angles + * in their potential expression + */ +template +class CosineParamAngle : public TwoParameterInteraction +{ +public: + CosineParamAngle() = default; + //! \brief construct from angle given in radians + CosineParamAngle(ForceConstant f, Radians angle) : + TwoParameterInteraction{ f, std::cos(angle) } + { + } + + //! \brief construct from angle given in degrees + CosineParamAngle(ForceConstant f, Degrees angle) : + TwoParameterInteraction{ 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; + +/*! \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; + +/*! \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& 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 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 @@ -303,8 +506,9 @@ inline bool operator==(const ProperDihedral& a, const ProperDihedral& b) /*! \brief Improper Dihedral Implementation */ -struct ImproperDihedral : public TwoParameterInteraction +class ImproperDihedral : public TwoParameterInteraction { +public: ImproperDihedral() = default; ImproperDihedral(Radians phi, ForceConstant f) : TwoParameterInteraction{ f, phi } diff --git a/api/nblib/listed_forces/dataflow.hpp b/api/nblib/listed_forces/dataflow.hpp index 1d20bf704a..14b002b85c 100644 --- a/api/nblib/listed_forces/dataflow.hpp +++ b/api/nblib/listed_forces/dataflow.hpp @@ -200,6 +200,82 @@ auto computeThreeCenter(const ThreeCenterType& parameters, const BasicVector& ri return energy; } +template +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; + + 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 +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; + // 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 +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; + // 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 diff --git a/api/nblib/listed_forces/definitions.h b/api/nblib/listed_forces/definitions.h index beb75a2a63..38bf00f80b 100644 --- a/api/nblib/listed_forces/definitions.h +++ b/api/nblib/listed_forces/definitions.h @@ -66,8 +66,9 @@ namespace nblib */ using SupportedTwoCenterTypes = - TypeList; -using SupportedThreeCenterTypes = TypeList; + TypeList; +using SupportedThreeCenterTypes = + TypeList; using SupportedFourCenterTypes = TypeList; using SupportedFiveCenterTypes = TypeList; @@ -131,10 +132,10 @@ struct ListedTypeData { using type = InteractionType; - // tuple format: - std::vector> indices; // vector of unique TwoCenterType instances std::vector parameters; + // tuple format: + std::vector> indices; }; using TwoCenterInteraction = Reduce; diff --git a/api/nblib/listed_forces/kernels.hpp b/api/nblib/listed_forces/kernels.hpp index b9064ce211..b2a28a812f 100644 --- a/api/nblib/listed_forces/kernels.hpp +++ b/api/nblib/listed_forces/kernels.hpp @@ -256,6 +256,40 @@ inline auto bondKernel(T dr, const MorseBondType& bond) 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 + */ +template +inline std::tuple 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 +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 @@ -397,6 +431,36 @@ inline auto bondKernel(T dr, const HalfAttractiveQuarticBondType& bond) //! 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 + */ +template +inline std::tuple 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 +inline auto threeCenterKernel(T dr, const LinearAngle& angle) +{ + return linearAnglesScalar(angle.forceConstant(), angle.equilConstant(), dr); +} //! Harmonic Angle template @@ -405,6 +469,134 @@ inline auto threeCenterKernel(T dr, const HarmonicAngle& angle) return harmonicScalarForce(angle.forceConstant(), angle.equilConstant(), dr); } +//! Cosine based (GROMOS-96) Angle +template +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 + */ + +template +inline std::tuple 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 +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 + */ + +template +inline std::tuple 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 +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 +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 +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 diff --git a/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesCrossBondAngleTest.xml b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesCrossBondAngleTest.xml new file mode 100644 index 0000000000..3e85dd72e5 --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesCrossBondAngleTest.xml @@ -0,0 +1,23 @@ + + + + 2.5625292086196829 + + 3 + + 20.584931942381417 + 32.355853840934138 + -40.148174152833434 + + + -1.5714394048923204 + 1.7263018958837222 + 0.65914251711356542 + + + -19.013492537489096 + -34.08215573681786 + 39.489031635719869 + + + diff --git a/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesCrossBondBondTest.xml b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesCrossBondBondTest.xml new file mode 100644 index 0000000000..88286b4ca2 --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesCrossBondBondTest.xml @@ -0,0 +1,23 @@ + + + + 15.930047836493104 + + 3 + + 16.257587884629597 + 2.253527033513016 + -18.35014870146313 + + + -18.352335074096068 + 23.83559705347821 + 5.5912340019856792 + + + 2.0947471894664718 + -26.089124086991227 + 12.758914699477451 + + + diff --git a/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesG96AngleTest.xml b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesG96AngleTest.xml new file mode 100644 index 0000000000..7c0ddec72e --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesG96AngleTest.xml @@ -0,0 +1,23 @@ + + + + 0.878581203916453 + + 3 + + 19.022908985367618 + -52.884817112836259 + 10.359003227565086 + + + -61.096971498538466 + 67.108390759454124 + 25.632702402756305 + + + 42.074062513170844 + -14.223573646617863 + -35.991705630321391 + + + diff --git a/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesLinearAngleTest.xml b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesLinearAngleTest.xml new file mode 100644 index 0000000000..e9f59bba20 --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesLinearAngleTest.xml @@ -0,0 +1,23 @@ + + + + 0.20264300000000018 + + 3 + + -0.94000000000000127 + 1.532 + 0.10800000000000157 + + + 2.3500000000000032 + -3.8300000000000001 + -0.2700000000000039 + + + -1.4100000000000019 + 2.298 + 0.16200000000000236 + + + diff --git a/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesQuarticAngleTest.xml b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesQuarticAngleTest.xml new file mode 100644 index 0000000000..1e299bea2e --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesQuarticAngleTest.xml @@ -0,0 +1,23 @@ + + + + 3.7653963456510402 + + 3 + + 28.524143507898081 + -79.298813545069052 + 15.532939514620491 + + + -91.612633181538157 + 100.62653246562743 + 38.435282552939967 + + + 63.088489673640076 + -21.327718920558382 + -53.968222067560454 + + + diff --git a/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesRestrictedAngleTest.xml b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesRestrictedAngleTest.xml new file mode 100644 index 0000000000..7594a57dc8 --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/ThreeCenter_ListedForcesRestrictedAngleTest.xml @@ -0,0 +1,23 @@ + + + + 1.0103314126941512 + + 3 + + 23.578508969297374 + -65.549655712163073 + 12.839773911655692 + + + -75.728454127772068 + 83.179486095108572 + 31.771213539191535 + + + 52.149945158474708 + -17.629830382945496 + -44.610987450847226 + + + diff --git a/api/nblib/listed_forces/tests/refdata/TwoCenter_ListedForcesMorseBondTest.xml b/api/nblib/listed_forces/tests/refdata/TwoCenter_ListedForcesMorseBondTest.xml new file mode 100644 index 0000000000..2f09194b2c --- /dev/null +++ b/api/nblib/listed_forces/tests/refdata/TwoCenter_ListedForcesMorseBondTest.xml @@ -0,0 +1,18 @@ + + + + 0.0018852207068626682 + + 2 + + -0.84131268781639024 + -0.11661760029138095 + 0.94960045951553052 + + + 0.84131268781639024 + 0.11661760029138095 + -0.94960045951553052 + + + diff --git a/api/nblib/listed_forces/tests/typetests.cpp b/api/nblib/listed_forces/tests/typetests.cpp index af26692cb6..4eb1684cfa 100644 --- a/api/nblib/listed_forces/tests/typetests.cpp +++ b/api/nblib/listed_forces/tests/typetests.cpp @@ -76,12 +76,35 @@ static const std::vector c_InputG96Bonds = { { G96BondType(50.0, 0. //! Function types for testing cubic bonds static const std::vector c_InputCubicBonds = { { CubicBondType(50.0, 2.0, 0.16) } }; +//! Function types for testing Morse bonds +static const std::vector c_InputMorseBonds = { { MorseBondType(30.0, 2.7, 0.15) } }; + //! Function types for testing FENE bonds static const std::vector c_InputFeneBonds = { { FENEBondType(5.0, 0.4) } }; //! Function types for testing Harmonic angles static const std::vector c_InputHarmonicAngles = { { HarmonicAngle(50.0, Degrees(100)) } }; +//! Function types for testing Linear angles +static const std::vector c_InputLinearAngles = { { LinearAngle(50.0, 0.4) } }; + +//! Function types for testing G96 angles +static const std::vector c_InputG96Angles = { { G96Angle(50.0, Degrees(100)) } }; + +//! Function types for testing Restricted angles +static const std::vector c_InputRestrictedAngles = { { RestrictedAngle(50.0, Degrees(100)) } }; + +//! Function types for testing Quartic angles +static const std::vector 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 c_InputCrossBondBond = { { CrossBondBond(45.0, 0.8, 0.7) } }; + +//! Function types for testing cross bond-angle interaction +static const std::vector c_InputCrossBondAngle = { { CrossBondAngle(45.0, 0.8, 0.7, 0.3) } }; + //! Function types for testing dihedrals static const std::vector c_InputDihs = { { ProperDihedral(Degrees(-105.0), 15.0, 2) } }; @@ -108,11 +131,41 @@ TEST(FourCenter, ListedForcesProperDihedralTest) 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); @@ -128,6 +181,11 @@ TEST(TwoCenter, ListedForcesCubicBondTest) checkForcesAndEnergiesWithRefData(c_InputCubicBonds, c_coordinatesForBondTests); } +TEST(TwoCenter, ListedForcesMorseBondTest) +{ + checkForcesAndEnergiesWithRefData(c_InputMorseBonds, c_coordinatesForBondTests); +} + TEST(TwoCenter, ListedForcesFeneBondTest) { checkForcesAndEnergiesWithRefData(c_InputFeneBonds, c_coordinatesForBondTests);