Introduce AngleInteractionType
authorejjordan <ejjordan@kth.se>
Mon, 21 Jun 2021 14:34:29 +0000 (16:34 +0200)
committerMark Abraham <mark.j.abraham@gmail.com>
Thu, 24 Jun 2021 09:48:58 +0000 (09:48 +0000)
This adds a template from which subsequent angle types can be generated.

api/nblib/listed_forces/bondtypes.h
api/nblib/listed_forces/kernels.hpp
api/nblib/listed_forces/tests/calculator.cpp
api/nblib/listed_forces/tests/conversions.cpp
api/nblib/listed_forces/tests/helpers.cpp
api/nblib/listed_forces/tests/linear_chain_input.hpp
api/nblib/listed_forces/tests/typetests.cpp
api/nblib/samples/methane-water-integration.cpp
api/nblib/tests/molecules.cpp
api/nblib/tests/testsystems.cpp
api/nblib/tests/topology.cpp

index 7e281c045272aaa7446bac0ad5951d921836f5e4..8e9844f76ab855cb787d9efc4f2ded441d98c19a 100644 (file)
@@ -216,26 +216,42 @@ inline bool operator==(const MorseBondType& a, const MorseBondType& b)
 }
 
 
-/*! \brief Harmonic angle type
+/*! \brief Basic template for interactions with 2 parameters named forceConstant and equilAngle
+ *
+ * \tparam Phantom unused template parameter for type distinction
+ *
+ * Distinct angle types can be generated from this template with using declarations
+ * and declared, but undefined structs. For example:
+ * using HarmonicAngleType = AngleInteractionType<struct HarmonicAngleParameter>;
+ * HarmonicAngleParameter does not have to be defined.
  *
  * Note: the angle is always stored as radians internally
  */
-struct HarmonicAngle : public TwoParameterInteraction<struct HarmonicAngleTypeParameter>
+template<class Phantom>
+class AngleInteractionType : public TwoParameterInteraction<Phantom>
 {
-    HarmonicAngle() = default;
+public:
+    AngleInteractionType() = default;
     //! \brief construct from angle given in radians
-    HarmonicAngle(Radians angle, ForceConstant f) :
-        TwoParameterInteraction<struct HarmonicAngleTypeParameter>{ f, angle }
+    AngleInteractionType(ForceConstant f, Radians angle) :
+        TwoParameterInteraction<Phantom>{ f, angle }
     {
     }
 
     //! \brief construct from angle given in degrees
-    HarmonicAngle(Degrees angle, ForceConstant f) :
-        TwoParameterInteraction<struct HarmonicAngleTypeParameter>{ f, angle * DEG2RAD }
+    AngleInteractionType(ForceConstant f, Degrees angle) :
+        TwoParameterInteraction<Phantom>{ f, angle * DEG2RAD }
     {
     }
 };
 
+/*! \brief Harmonic angle type
+ *
+ * It represents the interaction of the form
+ * V(theta, forceConstant, equilAngle) = 0.5 * forceConstant * (theta - equilAngle)^2
+ */
+using HarmonicAngle = AngleInteractionType<struct HarmonicAngleParameter>;
+
 /*! \brief Proper Dihedral Implementation
  */
 class ProperDihedral
index 5b98c1c32f842eabbfbbcb9d4d146e020a021b2f..ba4b1abfad84339a673cce635970c7f06ac63a31 100644 (file)
@@ -533,13 +533,13 @@ inline void spreadTwoCenterForces(const T bondForce,
  * @param force_k
  */
 template <class T>
-inline void spreadThreeCenterForces(T cos_theta,
-                             T force,
-                             const gmx::RVec& r_ij,
-                             const gmx::RVec& r_kj,
-                             gmx::RVec* force_i,
-                             gmx::RVec* force_j,
-                             gmx::RVec* force_k)
+inline void spreadThreeCenterForces(T                          cos_theta,
+                                    T                          force,
+                                    const gmx::BasicVector<T>& r_ij,
+                                    const gmx::BasicVector<T>& r_kj,
+                                    gmx::BasicVector<T>*       force_i,
+                                    gmx::BasicVector<T>*       force_j,
+                                    gmx::BasicVector<T>*       force_k)
 {
     T cos_theta2 = cos_theta * cos_theta;
     if (cos_theta2 < 1)
@@ -556,9 +556,9 @@ inline void spreadThreeCenterForces(T cos_theta,
         T cii = sth * nrij_1 * nrij_1; /*   2          */
         T ckk = sth * nrkj_1 * nrkj_1; /*   2          */
 
-        gmx::RVec f_i{0, 0, 0};
-        gmx::RVec f_j{0, 0, 0};
-        gmx::RVec f_k{0, 0, 0};
+        gmx::BasicVector<T> f_i{0, 0, 0};
+        gmx::BasicVector<T> f_j{0, 0, 0};
+        gmx::BasicVector<T> f_k{0, 0, 0};
         for (int m = 0; m < dimSize; m++)
         { /*  39               */
             f_i[m] = -(cik * r_kj[m] - cii * r_ij[m]);
index 4a56cca2731b5da01b69fc110e82d5d45ab979d3..447ab2385d46d8acab094066bb112fc0c0e39e0c 100644 (file)
@@ -104,7 +104,7 @@ protected:
         // one bond between atoms 0-1 with bond1 parameters and another between atoms 1-2 with bond2 parameters
         std::vector<InteractionIndex<HarmonicBondType>> bondIndices{ { 0, 1, 0 }, { 1, 2, 1 } };
 
-        HarmonicAngle                                angle(Degrees(108.53), 397.5);
+        HarmonicAngle                                angle(397.5, Degrees(108.53));
         std::vector<HarmonicAngle>                   angles{ angle };
         std::vector<InteractionIndex<HarmonicAngle>> angleIndices{ { 0, 1, 2, 0 } };
 
index 29dafc5099b0deeb40bb27b99cba28200bc38fa9..b77aa65eb2da337bb887a0978dcc8c2b37c30983 100644 (file)
@@ -65,8 +65,8 @@ ListedInteractionData someBondsAndAngles()
     std::vector<HarmonicBondType> bonds{ bond1, bond2 };
     pickType<HarmonicBondType>(interactions).parameters = bonds;
 
-    HarmonicAngle              angle1(Degrees(100), 100);
-    HarmonicAngle              angle2(Degrees(101), 200);
+    HarmonicAngle              angle1(100, Degrees(100));
+    HarmonicAngle              angle2(200, Degrees(101));
     std::vector<HarmonicAngle> angles{ angle1, angle2 };
     pickType<HarmonicAngle>(interactions).parameters = angles;
 
index 22fde45d78cc318c6d784074f95ea94e76f86db9..cb91699c8feb7f9c063e0efe7733064d92aa6b09 100644 (file)
@@ -62,7 +62,7 @@ TEST(NBlibTest, CanSplitListedWork)
 {
     ListedInteractionData interactions;
 
-    HarmonicAngle    angle(Degrees(1), 1);
+    HarmonicAngle    angle(1, Degrees(1));
     HarmonicBondType bond(1, 1);
 
     int largestIndex = 20;
index 05dc6480d6ca4f790d597c3670d6b88f2cc1b8e4..4bc72d4401b0e1acf174b5394c9a8cfe8cb71b3d 100644 (file)
@@ -63,7 +63,7 @@ public:
         std::vector<HarmonicBondType> bonds{ bond1, bond2 };
         pickType<HarmonicBondType>(interactions).parameters = bonds;
 
-        HarmonicAngle              angle(Degrees(179.9), 397.5);
+        HarmonicAngle              angle(397.5, Degrees(179.9));
         std::vector<HarmonicAngle> angles{ angle };
         pickType<HarmonicAngle>(interactions).parameters = angles;
 
index b67f3b6f540f92ec324be68d0d26b76e2ec92b82..2fe6a0c3976f129b9bddc2b4a3692e6e7cd41409 100644 (file)
@@ -67,7 +67,7 @@ std::vector<std::vector<HarmonicBondType>> c_InputHarmonicBond = { { HarmonicBon
 
 // Parameters for harmonic angles
 std::vector<InteractionIndex<HarmonicAngle>> c_HarmonicAngleIndices{ { 0, 1, 2, 0 }, { 1, 2, 3, 0 } };
-std::vector<std::vector<HarmonicAngle>> c_InputHarmonicAngle = { { HarmonicAngle(Degrees(100), 50.0) } };
+std::vector<std::vector<HarmonicAngle>> c_InputHarmonicAngle = { { HarmonicAngle(50.0, Degrees(100)) } };
 
 //! Function types for testing dihedrals. Add new terms at the end.
 std::vector<std::vector<ProperDihedral>> c_InputDihs = { { { ProperDihedral(Degrees(-105.0), 15.0, 2) } } /*, { ImproperDihedral(100.0, 50.0) }*/ };
index 8012de656b24d31d72ef7cbcc360eae4dc44f582..127433d232f245d51194e76a69ab27ddcd92c5bd 100644 (file)
@@ -96,8 +96,8 @@ int main()
     HarmonicBondType ohHarmonicBond(1, 1);
     HarmonicBondType hcHarmonicBond(2, 1);
 
-    HarmonicAngle hohAngle(Degrees(120), 1);
-    HarmonicAngle hchAngle(Degrees(109.5), 1);
+    HarmonicAngle hohAngle(1, Degrees(120));
+    HarmonicAngle hchAngle(1, Degrees(109.5));
 
     // add harmonic bonds for water
     water.addInteraction(ParticleName("O"), ParticleName("H1"), ohHarmonicBond);
index d4fb155d0e09de59e494f8ea1f1e832087d44c5e..e7cde25810616a6aff64359a1eef6e74aa51c8d5 100644 (file)
@@ -195,7 +195,7 @@ TEST(NBlibTest, CanAddInteractions)
 
     HarmonicBondType hb(1, 2);
     CubicBondType    cub(1, 2, 3);
-    HarmonicAngle    ang(Degrees(1), 1);
+    HarmonicAngle    ang(1, Degrees(1));
 
     molecule.addInteraction(ParticleName("O"), ParticleName("H1"), hb);
     molecule.addInteraction(ParticleName("O"), ParticleName("H2"), hb);
index db9ca65893e5fb3fd3631470fafb7e66fafe99fe..d0096792855d3d1cc673a01023ebbb6d7339e671 100644 (file)
@@ -166,7 +166,7 @@ MethanolMoleculeBuilder::MethanolMoleculeBuilder() : methanol_(MoleculeName("MeO
     HarmonicBondType ometBond(1.1, 1.2);
     methanol_.addInteraction(ParticleName("O2"), ParticleName("Me1"), ometBond);
 
-    HarmonicAngle ochAngle(Degrees(108.52), 397.5);
+    HarmonicAngle ochAngle(397.5, Degrees(108.52));
     methanol_.addInteraction(ParticleName("O2"), ParticleName("Me1"), ParticleName("H3"), ochAngle);
 }
 
index 84622a4bfb234bb8a71927a4d0960177696694a0..3e8d73e1162c11103ac18df78c8b6faddbfefc8a 100644 (file)
@@ -360,7 +360,7 @@ TEST(NBlibTest, TopologyListedInteractionsMultipleTypes)
     Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
 
     CubicBondType testBond(1., 1., 1.);
-    HarmonicAngle testAngle(Degrees(1), 1);
+    HarmonicAngle testAngle(1, Degrees(1));
 
     water.addInteraction(ParticleName("H1"), ParticleName("H2"), testBond);
     water.addInteraction(ParticleName("H1"), ParticleName("Oxygen"), ParticleName("H2"), testAngle);
@@ -406,7 +406,7 @@ TEST(NBlibTest, TopologyListedInteractionsMultipleTypes)
     EXPECT_EQ(cubicBondsReference, cubicBonds.parameters);
     EXPECT_EQ(cubicIndicesReference, cubicBonds.indices);
 
-    HarmonicAngle                                methanolAngle(Degrees(108.52), 397.5);
+    HarmonicAngle                                methanolAngle(397.5, Degrees(108.52));
     std::vector<HarmonicAngle>                   angleReference{ testAngle, methanolAngle };
     std::vector<InteractionIndex<HarmonicAngle>> angleIndicesReference{
         { std::min(H1, H2), Ow, std::max(H1, H2), 0 }, { std::min(MeH1, MeO1), Me1, std::max(MeO1, MeH1), 1 }