*/
class ConstraintsTest : public ::testing::TestWithParam<ConstraintsTestParameters>
{
- public:
- //! PBC setups
- std::unordered_map <std::string, t_pbc> pbcs_;
- //! Algorithms (SHAKE and LINCS)
- std::unordered_map <std::string, void(*)(ConstraintsTestData *testData, t_pbc pbc)> algorithms_;
-
- /*! \brief Test setup function.
- *
- * Setting up the pbcs and algorithms. Note, that corresponding string keywords
- * have to be explicitly added at the end of this file when the tests are called.
- *
- */
- void SetUp() override
- {
+public:
+ //! PBC setups
+ std::unordered_map<std::string, t_pbc> pbcs_;
+ //! Algorithms (SHAKE and LINCS)
+ std::unordered_map<std::string, void (*)(ConstraintsTestData* testData, t_pbc pbc)> algorithms_;
+
+ /*! \brief Test setup function.
+ *
+ * Setting up the pbcs and algorithms. Note, that corresponding string keywords
+ * have to be explicitly added at the end of this file when the tests are called.
+ *
+ */
+ void SetUp() override
+ {
- //
- // PBC initialization
- //
- t_pbc pbc;
-
- // Infinitely small box
- matrix boxNone = { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} };
- set_pbc(&pbc, epbcNONE, boxNone);
- pbcs_["PBCNone"] = pbc;
-
- // Rectangular box
- matrix boxXyz = { {10.0, 0.0, 0.0}, {0.0, 20.0, 0.0}, {0.0, 0.0, 15.0} };
- set_pbc(&pbc, epbcXYZ, boxXyz);
- pbcs_["PBCXYZ"] = pbc;
-
- //
- // Algorithms
- //
- // SHAKE
- algorithms_["SHAKE"] = applyShake;
- // LINCS
- algorithms_["LINCS"] = applyLincs;
- // LINCS using CUDA (will only be called if CUDA is available)
- algorithms_["LINCS_CUDA"] = applyLincsCuda;
- }
+ //
+ // PBC initialization
+ //
+ t_pbc pbc;
+
+ // Infinitely small box
+ matrix boxNone = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } };
+ set_pbc(&pbc, epbcNONE, boxNone);
+ pbcs_["PBCNone"] = pbc;
+
+ // Rectangular box
+ matrix boxXyz = { { 10.0, 0.0, 0.0 }, { 0.0, 20.0, 0.0 }, { 0.0, 0.0, 15.0 } };
+ set_pbc(&pbc, epbcXYZ, boxXyz);
+ pbcs_["PBCXYZ"] = pbc;
+
+ //
+ // Algorithms
+ //
+ // SHAKE
+ algorithms_["SHAKE"] = applyShake;
+ // LINCS
+ algorithms_["LINCS"] = applyLincs;
+ // LINCS using CUDA (will only be called if CUDA is available)
+ algorithms_["LINCS_CUDA"] = applyLincsCuda;
+ }
- /*! \brief
- * The test on the final length of constrained bonds.
- *
- * Goes through all the constraints and checks if the final length of all the constraints is equal
- * to the target length with provided tolerance.
- *
- * \param[in] tolerance Allowed tolerance in final lengths.
- * \param[in] testData Test data structure.
- * \param[in] pbc Periodic boundary data.
- */
- void checkConstrainsLength(FloatingPointTolerance tolerance, const ConstraintsTestData &testData, t_pbc pbc)
- {
+ /*! \brief
+ * The test on the final length of constrained bonds.
+ *
+ * Goes through all the constraints and checks if the final length of all the constraints is
+ * equal to the target length with provided tolerance.
+ *
+ * \param[in] tolerance Allowed tolerance in final lengths.
+ * \param[in] testData Test data structure.
+ * \param[in] pbc Periodic boundary data.
+ */
+ void checkConstrainsLength(FloatingPointTolerance tolerance, const ConstraintsTestData& testData, t_pbc pbc)
+ {
- // Test if all the constraints are satisfied
- for (index c = 0; c < ssize(testData.constraints_)/3; c++)
+ // Test if all the constraints are satisfied
+ for (index c = 0; c < ssize(testData.constraints_) / 3; c++)
+ {
+ real r0 = testData.constraintsR0_.at(testData.constraints_.at(3 * c));
+ int i = testData.constraints_.at(3 * c + 1);
+ int j = testData.constraints_.at(3 * c + 2);
+ RVec xij0, xij1;
+ real d0, d1;
+ if (pbc.ePBC == epbcXYZ)
{
- real r0 = testData.constraintsR0_.at(testData.constraints_.at(3*c));
- int i = testData.constraints_.at(3*c + 1);
- int j = testData.constraints_.at(3*c + 2);
- RVec xij0, xij1;
- real d0, d1;
- if (pbc.ePBC == epbcXYZ)
- {
- pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
- pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
- }
- else
- {
- rvec_sub(testData.x_[i], testData.x_[j], xij0);
- rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
- }
- d0 = norm(xij0);
- d1 = norm(xij1);
- EXPECT_REAL_EQ_TOL(r0, d1, tolerance) << gmx::formatString(
- "rij = %f, which is not equal to r0 = %f for constraint #%zd, between atoms %d and %d"
- " (before constraining rij was %f).", d1, r0, c, i, j, d0);
+ pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
+ pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
}
- }
-
- /*! \brief
- * The test on the final length of constrained bonds.
- *
- * Goes through all the constraints and checks if the direction of constraint has not changed
- * by the algorithm (i.e. the constraints algorithm arrived to the solution that is closest
- * to the initial system conformation).
- *
- * \param[in] testData Test data structure.
- * \param[in] pbc Periodic boundary data.
- */
- void checkConstrainsDirection(const ConstraintsTestData &testData, t_pbc pbc)
- {
-
- for (index c = 0; c < ssize(testData.constraints_)/3; c++)
+ else
{
- int i = testData.constraints_.at(3*c + 1);
- int j = testData.constraints_.at(3*c + 2);
- RVec xij0, xij1;
- if (pbc.ePBC == epbcXYZ)
- {
- pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
- pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
- }
- else
- {
- rvec_sub(testData.x_[i], testData.x_[j], xij0);
- rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
- }
-
- real dot = xij0.dot(xij1);
-
- EXPECT_GE(dot, 0.0) << gmx::formatString(
- "The constraint %zd changed direction. Constraining algorithm might have returned the wrong root "
- "of the constraints equation.", c);
-
+ rvec_sub(testData.x_[i], testData.x_[j], xij0);
+ rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
}
+ d0 = norm(xij0);
+ d1 = norm(xij1);
+ EXPECT_REAL_EQ_TOL(r0, d1, tolerance) << gmx::formatString(
+ "rij = %f, which is not equal to r0 = %f for constraint #%zd, between atoms %d "
+ "and %d"
+ " (before constraining rij was %f).",
+ d1, r0, c, i, j, d0);
}
+ }
- /*! \brief
- * The test on the coordinates of the center of the mass (COM) of the system.
- *
- * Checks if the center of mass has not been shifted by the constraints. Note,
- * that this test does not take into account the periodic boundary conditions.
- * Hence it will not work should the constraints decide to move atoms across
- * PBC borders.
- *
- * \param[in] tolerance Allowed tolerance in COM coordinates.
- * \param[in] testData Test data structure.
- */
- void checkCOMCoordinates(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
- {
+ /*! \brief
+ * The test on the final length of constrained bonds.
+ *
+ * Goes through all the constraints and checks if the direction of constraint has not changed
+ * by the algorithm (i.e. the constraints algorithm arrived to the solution that is closest
+ * to the initial system conformation).
+ *
+ * \param[in] testData Test data structure.
+ * \param[in] pbc Periodic boundary data.
+ */
+ void checkConstrainsDirection(const ConstraintsTestData& testData, t_pbc pbc)
+ {
- RVec comPrime0({0.0, 0.0, 0.0});
- RVec comPrime({0.0, 0.0, 0.0});
- for (int i = 0; i < testData.numAtoms_; i++)
+ for (index c = 0; c < ssize(testData.constraints_) / 3; c++)
+ {
+ int i = testData.constraints_.at(3 * c + 1);
+ int j = testData.constraints_.at(3 * c + 2);
+ RVec xij0, xij1;
+ if (pbc.ePBC == epbcXYZ)
{
- comPrime0 += testData.masses_[i]*testData.xPrime0_[i];
- comPrime += testData.masses_[i]*testData.xPrime_[i];
+ pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
+ pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
+ }
+ else
+ {
+ rvec_sub(testData.x_[i], testData.x_[j], xij0);
+ rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
}
- comPrime0 /= testData.numAtoms_;
- comPrime /= testData.numAtoms_;
-
- EXPECT_REAL_EQ_TOL(comPrime[XX], comPrime0[XX], tolerance)
- << "Center of mass was shifted by constraints in x-direction.";
- EXPECT_REAL_EQ_TOL(comPrime[YY], comPrime0[YY], tolerance)
- << "Center of mass was shifted by constraints in y-direction.";
- EXPECT_REAL_EQ_TOL(comPrime[ZZ], comPrime0[ZZ], tolerance)
- << "Center of mass was shifted by constraints in z-direction.";
+ real dot = xij0.dot(xij1);
+ EXPECT_GE(dot, 0.0) << gmx::formatString(
+ "The constraint %zd changed direction. Constraining algorithm might have "
+ "returned the wrong root "
+ "of the constraints equation.",
+ c);
}
+ }
- /*! \brief
- * The test on the velocity of the center of the mass (COM) of the system.
- *
- * Checks if the velocity of the center of mass has not changed.
- *
- * \param[in] tolerance Allowed tolerance in COM velocity components.
- * \param[in] testData Test data structure.
- */
- void checkCOMVelocity(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
+ /*! \brief
+ * The test on the coordinates of the center of the mass (COM) of the system.
+ *
+ * Checks if the center of mass has not been shifted by the constraints. Note,
+ * that this test does not take into account the periodic boundary conditions.
+ * Hence it will not work should the constraints decide to move atoms across
+ * PBC borders.
+ *
+ * \param[in] tolerance Allowed tolerance in COM coordinates.
+ * \param[in] testData Test data structure.
+ */
+ void checkCOMCoordinates(FloatingPointTolerance tolerance, const ConstraintsTestData& testData)
+ {
+
+ RVec comPrime0({ 0.0, 0.0, 0.0 });
+ RVec comPrime({ 0.0, 0.0, 0.0 });
+ for (int i = 0; i < testData.numAtoms_; i++)
{
+ comPrime0 += testData.masses_[i] * testData.xPrime0_[i];
+ comPrime += testData.masses_[i] * testData.xPrime_[i];
+ }
- RVec comV0({0.0, 0.0, 0.0});
- RVec comV({0.0, 0.0, 0.0});
- for (int i = 0; i < testData.numAtoms_; i++)
- {
- comV0 += testData.masses_[i]*testData.v0_[i];
- comV += testData.masses_[i]*testData.v_[i];
- }
- comV0 /= testData.numAtoms_;
- comV /= testData.numAtoms_;
-
- EXPECT_REAL_EQ_TOL(comV[XX], comV0[XX], tolerance)
- << "Velocity of the center of mass in x-direction has been changed by constraints.";
- EXPECT_REAL_EQ_TOL(comV[YY], comV0[YY], tolerance)
- << "Velocity of the center of mass in y-direction has been changed by constraints.";
- EXPECT_REAL_EQ_TOL(comV[ZZ], comV0[ZZ], tolerance)
- << "Velocity of the center of mass in z-direction has been changed by constraints.";
+ comPrime0 /= testData.numAtoms_;
+ comPrime /= testData.numAtoms_;
+
+ EXPECT_REAL_EQ_TOL(comPrime[XX], comPrime0[XX], tolerance)
+ << "Center of mass was shifted by constraints in x-direction.";
+ EXPECT_REAL_EQ_TOL(comPrime[YY], comPrime0[YY], tolerance)
+ << "Center of mass was shifted by constraints in y-direction.";
+ EXPECT_REAL_EQ_TOL(comPrime[ZZ], comPrime0[ZZ], tolerance)
+ << "Center of mass was shifted by constraints in z-direction.";
+ }
+
+ /*! \brief
+ * The test on the velocity of the center of the mass (COM) of the system.
+ *
+ * Checks if the velocity of the center of mass has not changed.
+ *
+ * \param[in] tolerance Allowed tolerance in COM velocity components.
+ * \param[in] testData Test data structure.
+ */
+ void checkCOMVelocity(FloatingPointTolerance tolerance, const ConstraintsTestData& testData)
+ {
+
+ RVec comV0({ 0.0, 0.0, 0.0 });
+ RVec comV({ 0.0, 0.0, 0.0 });
+ for (int i = 0; i < testData.numAtoms_; i++)
+ {
+ comV0 += testData.masses_[i] * testData.v0_[i];
+ comV += testData.masses_[i] * testData.v_[i];
}
+ comV0 /= testData.numAtoms_;
+ comV /= testData.numAtoms_;
+
+ EXPECT_REAL_EQ_TOL(comV[XX], comV0[XX], tolerance)
+ << "Velocity of the center of mass in x-direction has been changed by constraints.";
+ EXPECT_REAL_EQ_TOL(comV[YY], comV0[YY], tolerance)
+ << "Velocity of the center of mass in y-direction has been changed by constraints.";
+ EXPECT_REAL_EQ_TOL(comV[ZZ], comV0[ZZ], tolerance)
+ << "Velocity of the center of mass in z-direction has been changed by constraints.";
+ }
- /*! \brief
- * The test of virial tensor.
- *
- * Checks if the values in the scaled virial tensor are equal to pre-computed values.
- *
- * \param[in] tolerance Tolerance for the tensor values.
- * \param[in] testData Test data structure.
- */
- void checkVirialTensor(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
+ /*! \brief
+ * The test of virial tensor.
+ *
+ * Checks if the values in the scaled virial tensor are equal to pre-computed values.
+ *
+ * \param[in] tolerance Tolerance for the tensor values.
+ * \param[in] testData Test data structure.
+ */
+ void checkVirialTensor(FloatingPointTolerance tolerance, const ConstraintsTestData& testData)
+ {
+ for (int i = 0; i < DIM; i++)
{
- for (int i = 0; i < DIM; i++)
+ for (int j = 0; j < DIM; j++)
{
- for (int j = 0; j < DIM; j++)
- {
- EXPECT_REAL_EQ_TOL(testData.virialScaledRef_[i][j], testData.virialScaled_[i][j],
- tolerance) << gmx::formatString(
- "Values in virial tensor at [%d][%d] are not within the tolerance from reference value.", i, j);
- }
+ EXPECT_REAL_EQ_TOL(testData.virialScaledRef_[i][j], testData.virialScaled_[i][j], tolerance)
+ << gmx::formatString(
+ "Values in virial tensor at [%d][%d] are not within the "
+ "tolerance from reference value.",
+ i, j);
}
}
+ }
};
-TEST_P(ConstraintsTest, SingleConstraint){
- std::string title = "one constraint (e.g. OH)";
- int numAtoms = 2;
-
- std::vector<real> masses = {1.0, 12.0};
- std::vector<int> constraints = {0, 0, 1};
- std::vector<real> constraintsR0 = {0.1};
-
- real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
-
- std::vector<RVec> x = {{ 0.0, oneTenthOverSqrtTwo, 0.0 },
- { oneTenthOverSqrtTwo, 0.0, 0.0 }};
- std::vector<RVec> xPrime = {{ 0.01, 0.08, 0.01 },
- { 0.06, 0.01, -0.01 }};
- std::vector<RVec> v = {{ 1.0, 2.0, 3.0 },
- { 3.0, 2.0, 1.0 }};
-
- tensor virialScaledRef = {{-5.58e-04, 5.58e-04, 0.00e+00 },
- { 5.58e-04, -5.58e-04, 0.00e+00 },
- { 0.00e+00, 0.00e+00, 0.00e+00 }};
-
- real shakeTolerance = 0.0001;
- gmx_bool shakeUseSOR = false;
-
- int lincsNIter = 1;
- int lincslincsExpansionOrder = 4;
- real lincsWarnAngle = 30.0;
-
- std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
- (title, numAtoms, masses,
- constraints, constraintsR0,
- true, virialScaledRef,
- false, 0,
- real(0.0), real(0.001),
- x, xPrime, v,
- shakeTolerance, shakeUseSOR,
- lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
+TEST_P(ConstraintsTest, SingleConstraint)
+{
+ std::string title = "one constraint (e.g. OH)";
+ int numAtoms = 2;
+
+ std::vector<real> masses = { 1.0, 12.0 };
+ std::vector<int> constraints = { 0, 0, 1 };
+ std::vector<real> constraintsR0 = { 0.1 };
+
+ real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
+
+ std::vector<RVec> x = { { 0.0, oneTenthOverSqrtTwo, 0.0 }, { oneTenthOverSqrtTwo, 0.0, 0.0 } };
+ std::vector<RVec> xPrime = { { 0.01, 0.08, 0.01 }, { 0.06, 0.01, -0.01 } };
+ std::vector<RVec> v = { { 1.0, 2.0, 3.0 }, { 3.0, 2.0, 1.0 } };
+
+ tensor virialScaledRef = { { -5.58e-04, 5.58e-04, 0.00e+00 },
+ { 5.58e-04, -5.58e-04, 0.00e+00 },
+ { 0.00e+00, 0.00e+00, 0.00e+00 } };
+
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
+
+ int lincsNIter = 1;
+ int lincslincsExpansionOrder = 4;
+ real lincsWarnAngle = 30.0;
+
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>(
+ title, numAtoms, masses, constraints, constraintsR0, true, virialScaledRef, false, 0,
+ real(0.0), real(0.001), x, xPrime, v, shakeTolerance, shakeUseSOR, lincsNIter,
+ lincslincsExpansionOrder, lincsWarnAngle);
std::string pbcName;
std::string algorithmName;
std::tie(pbcName, algorithmName) = GetParam();
- t_pbc pbc = pbcs_.at(pbcName);
+ t_pbc pbc = pbcs_.at(pbcName);
// Apply constraints
algorithms_.at(algorithmName)(testData.get(), pbc);
checkCOMVelocity(absoluteTolerance(0.0001), *testData);
checkVirialTensor(absoluteTolerance(0.0001), *testData);
-
}
-TEST_P(ConstraintsTest, TwoDisjointConstraints){
+TEST_P(ConstraintsTest, TwoDisjointConstraints)
+{
- std::string title = "two disjoint constraints";
- int numAtoms = 4;
- std::vector<real> masses = {0.5, 1.0/3.0, 0.25, 1.0};
- std::vector<int> constraints = {0, 0, 1, 1, 2, 3};
- std::vector<real> constraintsR0 = {2.0, 1.0};
+ std::string title = "two disjoint constraints";
+ int numAtoms = 4;
+ std::vector<real> masses = { 0.5, 1.0 / 3.0, 0.25, 1.0 };
+ std::vector<int> constraints = { 0, 0, 1, 1, 2, 3 };
+ std::vector<real> constraintsR0 = { 2.0, 1.0 };
- std::vector<RVec> x = {{ 2.50, -3.10, 15.70 },
- { 0.51, -3.02, 15.55 },
- { -0.50, -3.00, 15.20 },
- { -1.51, -2.95, 15.05 }};
+ std::vector<RVec> x = {
+ { 2.50, -3.10, 15.70 }, { 0.51, -3.02, 15.55 }, { -0.50, -3.00, 15.20 }, { -1.51, -2.95, 15.05 }
+ };
- std::vector<RVec> xPrime = {{ 2.50, -3.10, 15.70 },
- { 0.51, -3.02, 15.55 },
- { -0.50, -3.00, 15.20 },
- { -1.51, -2.95, 15.05 }};
+ std::vector<RVec> xPrime = {
+ { 2.50, -3.10, 15.70 }, { 0.51, -3.02, 15.55 }, { -0.50, -3.00, 15.20 }, { -1.51, -2.95, 15.05 }
+ };
- std::vector<RVec> v = {{ 0.0, 1.0, 0.0 },
- { 1.0, 0.0, 0.0 },
- { 0.0, 0.0, 1.0 },
- { 0.0, 0.0, 0.0 }};
+ std::vector<RVec> v = { { 0.0, 1.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 0.0, 0.0, 1.0 }, { 0.0, 0.0, 0.0 } };
- tensor virialScaledRef = {{ 3.3e-03, -1.7e-04, 5.6e-04 },
- {-1.7e-04, 8.9e-06, -2.8e-05 },
- { 5.6e-04, -2.8e-05, 8.9e-05 }};
+ tensor virialScaledRef = { { 3.3e-03, -1.7e-04, 5.6e-04 },
+ { -1.7e-04, 8.9e-06, -2.8e-05 },
+ { 5.6e-04, -2.8e-05, 8.9e-05 } };
- real shakeTolerance = 0.0001;
- gmx_bool shakeUseSOR = false;
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
- int lincsNIter = 1;
- int lincslincsExpansionOrder = 4;
- real lincsWarnAngle = 30.0;
+ int lincsNIter = 1;
+ int lincslincsExpansionOrder = 4;
+ real lincsWarnAngle = 30.0;
- std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
- (title, numAtoms, masses,
- constraints, constraintsR0,
- true, virialScaledRef,
- false, 0,
- real(0.0), real(0.001),
- x, xPrime, v,
- shakeTolerance, shakeUseSOR,
- lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>(
+ title, numAtoms, masses, constraints, constraintsR0, true, virialScaledRef, false, 0,
+ real(0.0), real(0.001), x, xPrime, v, shakeTolerance, shakeUseSOR, lincsNIter,
+ lincslincsExpansionOrder, lincsWarnAngle);
std::string pbcName;
std::string algorithmName;
std::tie(pbcName, algorithmName) = GetParam();
- t_pbc pbc = pbcs_.at(pbcName);
+ t_pbc pbc = pbcs_.at(pbcName);
// Apply constraints
algorithms_.at(algorithmName)(testData.get(), pbc);
checkCOMVelocity(absoluteTolerance(0.0001), *testData);
checkVirialTensor(absoluteTolerance(0.0001), *testData);
-
}
-TEST_P(ConstraintsTest, ThreeSequentialConstraints){
+TEST_P(ConstraintsTest, ThreeSequentialConstraints)
+{
- std::string title = "three atoms, connected longitudinally (e.g. CH2)";
- int numAtoms = 3;
- std::vector<real> masses = {1.0, 12.0, 16.0 };
- std::vector<int> constraints = {0, 0, 1, 1, 1, 2};
- std::vector<real> constraintsR0 = {0.1, 0.2};
+ std::string title = "three atoms, connected longitudinally (e.g. CH2)";
+ int numAtoms = 3;
+ std::vector<real> masses = { 1.0, 12.0, 16.0 };
+ std::vector<int> constraints = { 0, 0, 1, 1, 1, 2 };
+ std::vector<real> constraintsR0 = { 0.1, 0.2 };
- real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
- real twoTenthsOverSqrtThree = 0.2_real / std::sqrt(3.0_real);
+ real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
+ real twoTenthsOverSqrtThree = 0.2_real / std::sqrt(3.0_real);
- std::vector<RVec> x = {{ oneTenthOverSqrtTwo, oneTenthOverSqrtTwo, 0.0 },
- { 0.0, 0.0, 0.0 },
- { twoTenthsOverSqrtThree, twoTenthsOverSqrtThree, twoTenthsOverSqrtThree }};
+ std::vector<RVec> x = { { oneTenthOverSqrtTwo, oneTenthOverSqrtTwo, 0.0 },
+ { 0.0, 0.0, 0.0 },
+ { twoTenthsOverSqrtThree, twoTenthsOverSqrtThree, twoTenthsOverSqrtThree } };
- std::vector<RVec> xPrime = {{ 0.08, 0.07, 0.01 },
- { -0.02, 0.01, -0.02 },
- { 0.10, 0.12, 0.11 }};
+ std::vector<RVec> xPrime = { { 0.08, 0.07, 0.01 }, { -0.02, 0.01, -0.02 }, { 0.10, 0.12, 0.11 } };
- std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
- { 0.0, 1.0, 0.0 },
- { 0.0, 0.0, 1.0 }};
+ std::vector<RVec> v = { { 1.0, 0.0, 0.0 }, { 0.0, 1.0, 0.0 }, { 0.0, 0.0, 1.0 } };
- tensor virialScaledRef = {{ 4.14e-03, 4.14e-03, 3.31e-03},
- { 4.14e-03, 4.14e-03, 3.31e-03},
- { 3.31e-03, 3.31e-03, 3.31e-03}};
+ tensor virialScaledRef = { { 4.14e-03, 4.14e-03, 3.31e-03 },
+ { 4.14e-03, 4.14e-03, 3.31e-03 },
+ { 3.31e-03, 3.31e-03, 3.31e-03 } };
- real shakeTolerance = 0.0001;
- gmx_bool shakeUseSOR = false;
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
- int lincsNIter = 1;
- int lincslincsExpansionOrder = 4;
- real lincsWarnAngle = 30.0;
+ int lincsNIter = 1;
+ int lincslincsExpansionOrder = 4;
+ real lincsWarnAngle = 30.0;
- std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
- (title, numAtoms, masses,
- constraints, constraintsR0,
- true, virialScaledRef,
- false, 0,
- real(0.0), real(0.001),
- x, xPrime, v,
- shakeTolerance, shakeUseSOR,
- lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>(
+ title, numAtoms, masses, constraints, constraintsR0, true, virialScaledRef, false, 0,
+ real(0.0), real(0.001), x, xPrime, v, shakeTolerance, shakeUseSOR, lincsNIter,
+ lincslincsExpansionOrder, lincsWarnAngle);
std::string pbcName;
std::string algorithmName;
std::tie(pbcName, algorithmName) = GetParam();
- t_pbc pbc = pbcs_.at(pbcName);
+ t_pbc pbc = pbcs_.at(pbcName);
// Apply constraints
algorithms_.at(algorithmName)(testData.get(), pbc);
checkCOMVelocity(absoluteTolerance(0.0001), *testData);
checkVirialTensor(absoluteTolerance(0.0001), *testData);
-
}
-TEST_P(ConstraintsTest, ThreeConstraintsWithCentralAtom){
+TEST_P(ConstraintsTest, ThreeConstraintsWithCentralAtom)
+{
- std::string title = "three atoms, connected to the central atom (e.g. CH3)";
- int numAtoms = 4;
- std::vector<real> masses = {12.0, 1.0, 1.0, 1.0 };
- std::vector<int> constraints = {0, 0, 1, 0, 0, 2, 0, 0, 3};
- std::vector<real> constraintsR0 = {0.1};
+ std::string title = "three atoms, connected to the central atom (e.g. CH3)";
+ int numAtoms = 4;
+ std::vector<real> masses = { 12.0, 1.0, 1.0, 1.0 };
+ std::vector<int> constraints = { 0, 0, 1, 0, 0, 2, 0, 0, 3 };
+ std::vector<real> constraintsR0 = { 0.1 };
- std::vector<RVec> x = {{ 0.00, 0.00, 0.00 },
- { 0.10, 0.00, 0.00 },
- { 0.00, -0.10, 0.00 },
- { 0.00, 0.00, 0.10 }};
+ std::vector<RVec> x = {
+ { 0.00, 0.00, 0.00 }, { 0.10, 0.00, 0.00 }, { 0.00, -0.10, 0.00 }, { 0.00, 0.00, 0.10 }
+ };
- std::vector<RVec> xPrime = {{ 0.004, 0.009, -0.010 },
- { 0.110, -0.006, 0.003 },
- {-0.007, -0.102, -0.007 },
- {-0.005, 0.011, 0.102 }};
+ std::vector<RVec> xPrime = { { 0.004, 0.009, -0.010 },
+ { 0.110, -0.006, 0.003 },
+ { -0.007, -0.102, -0.007 },
+ { -0.005, 0.011, 0.102 } };
- std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
- { 1.0, 0.0, 0.0 },
- { 1.0, 0.0, 0.0 },
- { 1.0, 0.0, 0.0 }};
+ std::vector<RVec> v = { { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 } };
- tensor virialScaledRef = {{7.14e-04, 0.00e+00, 0.00e+00},
- {0.00e+00, 1.08e-03, 0.00e+00},
- {0.00e+00, 0.00e+00, 1.15e-03}};
+ tensor virialScaledRef = { { 7.14e-04, 0.00e+00, 0.00e+00 },
+ { 0.00e+00, 1.08e-03, 0.00e+00 },
+ { 0.00e+00, 0.00e+00, 1.15e-03 } };
- real shakeTolerance = 0.0001;
- gmx_bool shakeUseSOR = false;
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
- int lincsNIter = 1;
- int lincslincsExpansionOrder = 4;
- real lincsWarnAngle = 30.0;
+ int lincsNIter = 1;
+ int lincslincsExpansionOrder = 4;
+ real lincsWarnAngle = 30.0;
- std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
- (title, numAtoms, masses,
- constraints, constraintsR0,
- true, virialScaledRef,
- false, 0,
- real(0.0), real(0.001),
- x, xPrime, v,
- shakeTolerance, shakeUseSOR,
- lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>(
+ title, numAtoms, masses, constraints, constraintsR0, true, virialScaledRef, false, 0,
+ real(0.0), real(0.001), x, xPrime, v, shakeTolerance, shakeUseSOR, lincsNIter,
+ lincslincsExpansionOrder, lincsWarnAngle);
std::string pbcName;
std::string algorithmName;
std::tie(pbcName, algorithmName) = GetParam();
- t_pbc pbc = pbcs_.at(pbcName);
+ t_pbc pbc = pbcs_.at(pbcName);
// Apply constraints
algorithms_.at(algorithmName)(testData.get(), pbc);
checkVirialTensor(absoluteTolerance(0.0001), *testData);
}
-TEST_P(ConstraintsTest, FourSequentialConstraints){
+TEST_P(ConstraintsTest, FourSequentialConstraints)
+{
- std::string title = "four atoms, connected longitudinally";
- int numAtoms = 4;
- std::vector<real> masses = {0.5, 1.0/3.0, 0.25, 1.0};
- std::vector<int> constraints = {0, 0, 1, 1, 1, 2, 2, 2, 3};
- std::vector<real> constraintsR0 = {2.0, 1.0, 1.0};
+ std::string title = "four atoms, connected longitudinally";
+ int numAtoms = 4;
+ std::vector<real> masses = { 0.5, 1.0 / 3.0, 0.25, 1.0 };
+ std::vector<int> constraints = { 0, 0, 1, 1, 1, 2, 2, 2, 3 };
+ std::vector<real> constraintsR0 = { 2.0, 1.0, 1.0 };
- std::vector<RVec> x = {{ 2.50, -3.10, 15.70 },
- { 0.51, -3.02, 15.55 },
- { -0.50, -3.00, 15.20 },
- { -1.51, -2.95, 15.05 }};
+ std::vector<RVec> x = {
+ { 2.50, -3.10, 15.70 }, { 0.51, -3.02, 15.55 }, { -0.50, -3.00, 15.20 }, { -1.51, -2.95, 15.05 }
+ };
- std::vector<RVec> xPrime = {{ 2.50, -3.10, 15.70 },
- { 0.51, -3.02, 15.55 },
- { -0.50, -3.00, 15.20 },
- { -1.51, -2.95, 15.05 }};
+ std::vector<RVec> xPrime = {
+ { 2.50, -3.10, 15.70 }, { 0.51, -3.02, 15.55 }, { -0.50, -3.00, 15.20 }, { -1.51, -2.95, 15.05 }
+ };
- std::vector<RVec> v = {{ 0.0, 0.0, 2.0 },
- { 0.0, 0.0, 3.0 },
- { 0.0, 0.0, -4.0 },
- { 0.0, 0.0, -1.0 }};
+ std::vector<RVec> v = { { 0.0, 0.0, 2.0 }, { 0.0, 0.0, 3.0 }, { 0.0, 0.0, -4.0 }, { 0.0, 0.0, -1.0 } };
- tensor virialScaledRef = {{ 1.15e-01, -4.20e-03, 2.12e-02},
- {-4.20e-03, 1.70e-04, -6.41e-04},
- { 2.12e-02, -6.41e-04, 5.45e-03}};
+ tensor virialScaledRef = { { 1.15e-01, -4.20e-03, 2.12e-02 },
+ { -4.20e-03, 1.70e-04, -6.41e-04 },
+ { 2.12e-02, -6.41e-04, 5.45e-03 } };
- real shakeTolerance = 0.0001;
- gmx_bool shakeUseSOR = false;
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
- int lincsNIter = 4;
- int lincslincsExpansionOrder = 8;
- real lincsWarnAngle = 30.0;
+ int lincsNIter = 4;
+ int lincslincsExpansionOrder = 8;
+ real lincsWarnAngle = 30.0;
- std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
- (title, numAtoms, masses,
- constraints, constraintsR0,
- true, virialScaledRef,
- false, 0,
- real(0.0), real(0.001),
- x, xPrime, v,
- shakeTolerance, shakeUseSOR,
- lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>(
+ title, numAtoms, masses, constraints, constraintsR0, true, virialScaledRef, false, 0,
+ real(0.0), real(0.001), x, xPrime, v, shakeTolerance, shakeUseSOR, lincsNIter,
+ lincslincsExpansionOrder, lincsWarnAngle);
std::string pbcName;
std::string algorithmName;
std::tie(pbcName, algorithmName) = GetParam();
- t_pbc pbc = pbcs_.at(pbcName);
+ t_pbc pbc = pbcs_.at(pbcName);
// Apply constraints
algorithms_.at(algorithmName)(testData.get(), pbc);
checkCOMVelocity(absoluteTolerance(0.0001), *testData);
checkVirialTensor(absoluteTolerance(0.01), *testData);
-
}
-TEST_P(ConstraintsTest, TriangleOfConstraints){
+TEST_P(ConstraintsTest, TriangleOfConstraints)
+{
- std::string title = "basic triangle (tree atoms, connected to each other)";
- int numAtoms = 3;
- std::vector<real> masses = {1.0, 1.0, 1.0};
- std::vector<int> constraints = {0, 0, 1, 2, 0, 2, 1, 1, 2};
- std::vector<real> constraintsR0 = {0.1, 0.1, 0.1};
+ std::string title = "basic triangle (tree atoms, connected to each other)";
+ int numAtoms = 3;
+ std::vector<real> masses = { 1.0, 1.0, 1.0 };
+ std::vector<int> constraints = { 0, 0, 1, 2, 0, 2, 1, 1, 2 };
+ std::vector<real> constraintsR0 = { 0.1, 0.1, 0.1 };
- real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
+ real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
- std::vector<RVec> x = {{ oneTenthOverSqrtTwo, 0.0, 0.0 },
- { 0.0, oneTenthOverSqrtTwo, 0.0 },
- { 0.0, 0.0, oneTenthOverSqrtTwo }};
+ std::vector<RVec> x = { { oneTenthOverSqrtTwo, 0.0, 0.0 },
+ { 0.0, oneTenthOverSqrtTwo, 0.0 },
+ { 0.0, 0.0, oneTenthOverSqrtTwo } };
- std::vector<RVec> xPrime = {{ 0.09, -0.02, 0.01 },
- { -0.02, 0.10, -0.02 },
- { 0.03, -0.01, 0.07 }};
+ std::vector<RVec> xPrime = { { 0.09, -0.02, 0.01 }, { -0.02, 0.10, -0.02 }, { 0.03, -0.01, 0.07 } };
- std::vector<RVec> v = {{ 1.0, 1.0, 1.0 },
- { -2.0, -2.0, -2.0 },
- { 1.0, 1.0, 1.0 }};
+ std::vector<RVec> v = { { 1.0, 1.0, 1.0 }, { -2.0, -2.0, -2.0 }, { 1.0, 1.0, 1.0 } };
- tensor virialScaledRef = {{ 6.00e-04, -1.61e-03, 1.01e-03},
- {-1.61e-03, 2.53e-03, -9.25e-04},
- { 1.01e-03, -9.25e-04, -8.05e-05}};
+ tensor virialScaledRef = { { 6.00e-04, -1.61e-03, 1.01e-03 },
+ { -1.61e-03, 2.53e-03, -9.25e-04 },
+ { 1.01e-03, -9.25e-04, -8.05e-05 } };
- real shakeTolerance = 0.0001;
- gmx_bool shakeUseSOR = false;
+ real shakeTolerance = 0.0001;
+ gmx_bool shakeUseSOR = false;
- int lincsNIter = 1;
- int lincslincsExpansionOrder = 4;
- real lincsWarnAngle = 30.0;
+ int lincsNIter = 1;
+ int lincslincsExpansionOrder = 4;
+ real lincsWarnAngle = 30.0;
- std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
- (title, numAtoms, masses,
- constraints, constraintsR0,
- true, virialScaledRef,
- false, 0,
- real(0.0), real(0.001),
- x, xPrime, v,
- shakeTolerance, shakeUseSOR,
- lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
+ std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>(
+ title, numAtoms, masses, constraints, constraintsR0, true, virialScaledRef, false, 0,
+ real(0.0), real(0.001), x, xPrime, v, shakeTolerance, shakeUseSOR, lincsNIter,
+ lincslincsExpansionOrder, lincsWarnAngle);
std::string pbcName;
std::string runnerName;
std::tie(pbcName, runnerName) = GetParam();
- t_pbc pbc = pbcs_.at(pbcName);
+ t_pbc pbc = pbcs_.at(pbcName);
// Apply constraints
algorithms_.at(runnerName)(testData.get(), pbc);
checkCOMVelocity(absoluteTolerance(0.0001), *testData);
checkVirialTensor(absoluteTolerance(0.00001), *testData);
-
}
-INSTANTIATE_TEST_CASE_P(WithParameters, ConstraintsTest,
- ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
- ::testing::ValuesIn(getRunnersNames())));
+INSTANTIATE_TEST_CASE_P(WithParameters,
+ ConstraintsTest,
+ ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
+ ::testing::ValuesIn(getRunnersNames())));
} // namespace
} // namespace test
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff