Expand signatures for nblib listed forces calculator
[alexxy/gromacs.git] / api / nblib / listed_forces / tests / calculator.cpp
index 7987f350af47f4fa12b0d069ff05e3e89851cbd6..1c7bff1b4e4b07ae9510665b0b983d35f7da3de0 100644 (file)
@@ -83,7 +83,9 @@ void compareVectors(const TestSeq&                    forces,
         for (int m = 0; m < dimSize; ++m)
         {
             EXPECT_FLOAT_DOUBLE_EQ_TOL(
-                    forces[i][m], refForcesFloat[i][m], refForcesDouble[i][m],
+                    forces[i][m],
+                    refForcesFloat[i][m],
+                    refForcesDouble[i][m],
                     // Todo: why does the tolerance need to be so low?
                     gmx::test::relativeToleranceAsFloatingPoint(refForcesDouble[i][m], 5e-5));
         }
@@ -102,48 +104,22 @@ 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 } };
 
-        HarmonicAngleType                                angle(Degrees(108.53), 397.5);
-        std::vector<HarmonicAngleType>                   angles{ angle };
-        std::vector<InteractionIndex<HarmonicAngleType>> angleIndices{ { 0, 1, 2, 0 } };
+        HarmonicAngle                                angle(397.5, Degrees(108.53));
+        std::vector<HarmonicAngle>                   angles{ angle };
+        std::vector<InteractionIndex<HarmonicAngle>> angleIndices{ { 0, 1, 2, 0 } };
 
         pickType<HarmonicBondType>(interactions).indices    = bondIndices;
         pickType<HarmonicBondType>(interactions).parameters = bonds;
 
-        pickType<HarmonicAngleType>(interactions).indices    = angleIndices;
-        pickType<HarmonicAngleType>(interactions).parameters = angles;
+        pickType<HarmonicAngle>(interactions).indices    = angleIndices;
+        pickType<HarmonicAngle>(interactions).parameters = angles;
 
         // initial position for the methanol atoms from the spc-water example
         x = std::vector<gmx::RVec>{ { 1.97, 1.46, 1.209 }, { 1.978, 1.415, 1.082 }, { 1.905, 1.46, 1.03 } };
         forces = std::vector<gmx::RVec>(3, gmx::RVec{ 0, 0, 0 });
 
-        refBondForcesFloat =
-                std::valarray<gmx::BasicVector<float>>{ { -22.8980637, 128.801575, 363.505951 },
-                                                        { -43.2698593, -88.0130997, -410.639252 },
-                                                        { 66.167923, -40.788475, 47.1333084 } };
-        refAngleForcesFloat =
-                std::valarray<gmx::BasicVector<float>>{ { 54.7276611, -40.1688995, 17.6805191 },
-                                                        { -81.8118973, 86.1988525, 60.1752243 },
-                                                        { 27.0842342, -46.0299492, -77.8557434 } };
-
-        refBondForcesDouble = std::valarray<gmx::BasicVector<double>>{
-            { -22.89764839974935, 128.79927224858977, 363.50016834602064 },
-            { -43.24622441913251, -88.025652017772231, -410.61635172385434 },
-            { 66.14387281888186, -40.773620230817542, 47.116183377833721 }
-        };
-        refAngleForcesDouble = std::valarray<gmx::BasicVector<double>>{
-            { 54.726206806506234, -40.167809526198099, 17.680008528590257 },
-            { -81.809781666748606, 86.196545126117257, 60.173723525141448 },
-            { 27.083574860242372, -46.028735599919159, -77.853732053731704 }
-        };
-
-        refBondEnergyFloat  = 0.2113433;
-        refAngleEnergyFloat = 0.112774156;
-
-        refBondEnergyDouble  = 0.2113273434867636;
-        refAngleEnergyDouble = 0.11276812148357591;
-
         box.reset(new Box(3, 3, 3));
-        pbc.reset(new PbcHolder(*box));
+        pbc.reset(new PbcHolder(PbcType::Xyz, *box));
     }
 
     std::vector<gmx::RVec> x;
@@ -153,50 +129,56 @@ protected:
 
     std::shared_ptr<Box>       box;
     std::shared_ptr<PbcHolder> pbc;
-
-    // reference values
-    std::valarray<gmx::BasicVector<float>>  refBondForcesFloat, refAngleForcesFloat;
-    std::valarray<gmx::BasicVector<double>> refBondForcesDouble, refAngleForcesDouble;
-
-    float  refBondEnergyFloat, refAngleEnergyFloat;
-    double refBondEnergyDouble, refAngleEnergyDouble;
 };
 
-TEST_F(ListedExampleData, DISABLED_ComputeHarmonicBondForces)
+TEST_F(ListedExampleData, ComputeHarmonicBondForces)
 {
-    auto indices = pickType<HarmonicBondType>(interactions).indices;
-    auto bonds   = pickType<HarmonicBondType>(interactions).parameters;
-    real energy  = computeForces(indices, bonds, x, &forces, *pbc);
+    gmx::ArrayRef<const InteractionIndex<HarmonicBondType>> indices =
+            pickType<HarmonicBondType>(interactions).indices;
+    gmx::ArrayRef<const HarmonicBondType> bonds = pickType<HarmonicBondType>(interactions).parameters;
+    computeForces(indices, bonds, x, &forces, *pbc);
+
+    RefDataChecker vector3DTest(1e-3);
+    vector3DTest.testArrays<Vec3>(forces, "Bond forces");
+}
 
-    EXPECT_FLOAT_DOUBLE_EQ_TOL(energy, refBondEnergyFloat, refBondEnergyDouble,
-                               gmx::test::relativeToleranceAsFloatingPoint(refBondEnergyDouble, 1e-5));
+TEST_F(ListedExampleData, ComputeHarmonicBondEnergies)
+{
+    gmx::ArrayRef<const InteractionIndex<HarmonicBondType>> indices =
+            pickType<HarmonicBondType>(interactions).indices;
+    gmx::ArrayRef<const HarmonicBondType> bonds = pickType<HarmonicBondType>(interactions).parameters;
+    real                                  energy = computeForces(indices, bonds, x, &forces, *pbc);
 
-    compareVectors(forces, refBondForcesFloat, refBondForcesDouble);
+    RefDataChecker vector3DTest(1e-4);
+    vector3DTest.testReal(energy, "Bond energy");
 }
 
 TEST_F(ListedExampleData, ComputeHarmonicAngleForces)
 {
-    auto indices = pickType<HarmonicAngleType>(interactions).indices;
-    auto angles  = pickType<HarmonicAngleType>(interactions).parameters;
-    real energy  = computeForces(indices, angles, x, &forces, *pbc);
+    gmx::ArrayRef<const InteractionIndex<HarmonicAngle>> indices =
+            pickType<HarmonicAngle>(interactions).indices;
+    gmx::ArrayRef<const HarmonicAngle> angles = pickType<HarmonicAngle>(interactions).parameters;
+    computeForces(indices, angles, x, &forces, *pbc);
 
-    EXPECT_FLOAT_DOUBLE_EQ_TOL(energy, refAngleEnergyFloat, refAngleEnergyDouble,
-                               gmx::test::relativeToleranceAsFloatingPoint(refAngleEnergyDouble, 1e-5));
+    RefDataChecker vector3DTest(1e-4);
+    vector3DTest.testArrays<Vec3>(forces, "Angle forces");
+}
 
-    compareVectors(forces, refAngleForcesFloat, refAngleForcesDouble);
+TEST_F(ListedExampleData, CanReduceForces)
+{
+    reduceListedForces(interactions, x, &forces, *pbc);
+
+    RefDataChecker vector3DTest(1e-2);
+    vector3DTest.testArrays<Vec3>(forces, "Reduced forces");
 }
 
-TEST_F(ListedExampleData, DISABLED_CanReduceForces)
+TEST_F(ListedExampleData, CanReduceEnergies)
 {
     auto energies    = reduceListedForces(interactions, x, &forces, *pbc);
     real totalEnergy = std::accumulate(begin(energies), end(energies), 0.0);
 
-    EXPECT_FLOAT_DOUBLE_EQ_TOL(totalEnergy, refBondEnergyFloat + refAngleEnergyFloat,
-                               refBondEnergyDouble + refAngleEnergyDouble,
-                               gmx::test::relativeToleranceAsFloatingPoint(refBondEnergyDouble, 1e-5));
-
-    compareVectors(forces, refBondForcesFloat + refAngleForcesFloat,
-                   refBondForcesDouble + refAngleForcesDouble);
+    RefDataChecker vector3DTest(1e-4);
+    vector3DTest.testReal(totalEnergy, "Reduced energy");
 }
 
 
@@ -205,7 +187,8 @@ void compareArray(const ListedForceCalculator::EnergyType& energies,
 {
     for (size_t i = 0; i < energies.size(); ++i)
     {
-        EXPECT_REAL_EQ_TOL(energies[i], refEnergies[i],
+        EXPECT_REAL_EQ_TOL(energies[i],
+                           refEnergies[i],
                            gmx::test::relativeToleranceAsFloatingPoint(refEnergies[i], 1e-5));
     }
 }
@@ -223,7 +206,6 @@ protected:
         interactions = data.interactions;
         box          = data.box;
         refForces    = data.forces;
-        // pbc.reset(new PbcHolder(*box));
 
         refEnergies = reduceListedForces(interactions, x, &refForces, NoPbc{});
     }
@@ -241,7 +223,6 @@ protected:
     std::vector<gmx::RVec> x;
     ListedInteractionData  interactions;
     std::shared_ptr<Box>   box;
-    // std::shared_ptr<PbcHolder> pbc;
 
 private:
     std::vector<gmx::RVec>            refForces;