This adds a template from which subsequent angle types can be generated.
}
-/*! \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
* @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)
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]);
// 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 } };
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;
{
ListedInteractionData interactions;
- HarmonicAngle angle(Degrees(1), 1);
+ HarmonicAngle angle(1, Degrees(1));
HarmonicBondType bond(1, 1);
int largestIndex = 20;
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;
// 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) }*/ };
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);
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);
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);
}
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);
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 }