Add Urey-Bradley listed interaction to nblib

diff --git a/api/nblib/molecules.h b/api/nblib/molecules.h
index 14c7f087c95032b23c08d97f3d838050ff963ce7..0ab4ec388425c03cee19326d103976453163ccb3 100644 (file)
--- a/api/nblib/molecules.h
+++ b/api/nblib/molecules.h
@@ -243,5 +243,39 @@ private:
     InteractionTuple interactionData_;
 };
 
+/*! \brief Helper function to add harmonic angle and harmonic bonds for Urey-Bradley term.
+ *
+ * Urey-Bradley consist of two harmonic terms:
+ *   1. Harmonic angle, connecting all three particles.
+ *   2. Harmonic correction to the distance between two non-central particles (particles 1 and 3)
+ * This function creates theese terms and adds them to the \c molecule as independent harmonic
+ * angle and harmonic bond.
+ *
+ * \todo This should be moved to another location (e.g. to TPR reader).
+ *
+ * \param[in,out] molecule   The molecule to add Urey-bradley to.
+ * \param[in]     particleI  First interacting particle.
+ * \param[in]     particleJ  Second (central) interacting particle.
+ * \param[in]     particleK  Third interacting particle.
+ * \param[in]     theta0     Equilibrium angle (in radians).
+ * \param[in]     kTheta     Force-constant for angle.
+ * \param[in]     r130       Equilibrium distance between particles 1 and 3.
+ * \param[in]     kUB        Force constant for bond correction term.
+ */
+static inline void addUreyBradleyInteraction(Molecule&                 molecule,
+                                             const ParticleIdentifier& particleI,
+                                             const ParticleIdentifier& particleJ,
+                                             const ParticleIdentifier& particleK,
+                                             const Radians             theta0,
+                                             const ForceConstant       kTheta,
+                                             const EquilConstant       r130,
+                                             const ForceConstant       kUB)
+{
+    HarmonicAngle    ubAngle(kTheta, theta0);
+    HarmonicBondType ubBond(kUB, r130);
+    molecule.addInteraction(particleI, particleJ, particleK, ubAngle);
+    molecule.addInteraction(particleI, particleK, ubBond);
+}
+
 } // namespace nblib
 #endif // NBLIB_MOLECULES_H

diff --git a/api/nblib/tests/molecules.cpp b/api/nblib/tests/molecules.cpp
index e7cde25810616a6aff64359a1eef6e74aa51c8d5..9f8af499a049e314b312f9373e49eff6402fdf02 100644 (file)
--- a/api/nblib/tests/molecules.cpp
+++ b/api/nblib/tests/molecules.cpp
@@ -212,6 +212,32 @@ TEST(NBlibTest, CanAddInteractions)
     EXPECT_EQ(pickType<HarmonicAngle>(interactionData).interactions_.size(), 1);
 }
 
+TEST(NBlibTest, CanAddUreyBradley)
+{
+    Molecule     molecule(MoleculeName("UreyBradleyTest"));
+    ParticleType O(ParticleTypeName("Ow"), Mass(1));
+    ParticleType H(ParticleTypeName("Hw"), Mass(1));
+    molecule.addParticle(ParticleName("O"), O);
+    molecule.addParticle(ParticleName("H1"), H);
+    molecule.addParticle(ParticleName("H2"), H);
+
+    addUreyBradleyInteraction(molecule,
+                              ParticleName("H1"),
+                              ParticleName("O"),
+                              ParticleName("H2"),
+                              Radians(1.82),
+                              ForceConstant(2.0),
+                              EquilConstant(0.15),
+                              ForceConstant(3.0));
+
+    const auto& interactionData = molecule.interactionData();
+
+    //! harmonic bonds
+    EXPECT_EQ(pickType<HarmonicBondType>(interactionData).interactions_.size(), 1);
+    //! angular interactions
+    EXPECT_EQ(pickType<HarmonicAngle>(interactionData).interactions_.size(), 1);
+}
+
 } // namespace
 } // namespace test
 } // namespace nblib