Add preprocessing and integration test for transformation coordinate

We have general checks for preprocessing in readir.cpp. This ensures that all relevant mdp options are validated properly.
The mathematical expression, however, is evaluated in the grompp integration test.
Then we have the mdrun integration test, which is based on the existing pull test for
a 3D umbrealla potential. Basically, we use the exact same setup and expect the same output, but apply the force on a transformation coordinate instead of a regular coordinate, only with a slightly lower numerical accuracy.

diff --git a/src/gromacs/gmxpreprocess/readpull.cpp b/src/gromacs/gmxpreprocess/readpull.cpp
index 2524976cb583bf942e49fdc43b1e51d74ade2a99..fed4eabb6dd329f15a9afdc14d76ed17d1b3ee44 100644 (file)
--- a/src/gromacs/gmxpreprocess/readpull.cpp
+++ b/src/gromacs/gmxpreprocess/readpull.cpp
@@ -134,23 +134,23 @@ static void process_pull_dim(char* dim_buf, ivec dim, const t_pull_coord* pcrd)
     }
 }
 
-static void initTransformationPullCoord(t_pull_coord* pcrd, const pull_params_t& pull)
+static void initTransformationPullCoord(t_pull_coord* pcrd, const pull_params_t& pull, warninp_t wi)
 {
     const int coord_index_for_output = pull.coord.size() + 1;
     if (pcrd->eType == PullingAlgorithm::Constraint)
     {
-        GMX_THROW(gmx::InvalidInputError(gmx::formatString(
-                "pull-coord%d can not have type 'constraint' and geometry 'transformation'",
-                coord_index_for_output)));
+        warning_error(
+                wi,
+                gmx::formatString(
+                        "pull-coord%d cannot have type 'constraint' and geometry 'transformation'",
+                        coord_index_for_output));
     }
 
     /*Validate the mathematical expression to epullgTRANSFORMATION*/
     if (pcrd->expression.empty())
     {
-        GMX_THROW(gmx::InvalidInputError(
-                gmx::formatString("pull-coord%d-expression not set for pull coordinate of geometry "
-                                  "'transformation'",
-                                  coord_index_for_output)));
+        warning_error(
+                wi, gmx::formatString("pull-coord%d-expression not set for pull coordinate of geometry 'transformation'", coord_index_for_output));
     }
     else if (pcrd->expression[0] == '"' || pcrd->expression[0] == '\'')
     {
@@ -160,17 +160,19 @@ static void initTransformationPullCoord(t_pull_coord* pcrd, const pull_params_t&
     }
     if (pcrd->dx == 0)
     {
-        GMX_THROW(gmx::InvalidInputError(gmx::formatString(
-                "pull-coord%d-dx cannot be set to zero for pull coordinate of geometry "
-                "'transformation'",
-                coord_index_for_output)));
+        warning_error(
+                wi,
+                gmx::formatString(
+                        "pull-coord%d-dx cannot be set to zero for pull coordinate of geometry "
+                        "'transformation'",
+                        coord_index_for_output));
     }
     /* make sure that the kappa of all previous pull coords is 0*/
     int previousCoordOutputIndex = 0;
     for (const auto& previousPcrd : pull.coord)
     {
         previousCoordOutputIndex++;
-        // See if the previous variable is used by the transformatino coord
+        // See if the previous variable is used by the transformation coord
         // Note that a simple std::string::find won't work since we don't want x1 to match x11 etc.
         std::string previousPcrdName = gmx::formatString("x%d(\\D|$)", previousCoordOutputIndex);
         std::regex  rx(previousPcrdName);
@@ -187,11 +189,25 @@ static void initTransformationPullCoord(t_pull_coord* pcrd, const pull_params_t&
 
         if (previousPcrd.eType == PullingAlgorithm::Constraint)
         {
-            GMX_THROW(gmx::InvalidInputError(
-                    gmx::formatString("pull-coord%d can not use pull-coord%d in the "
-                                      "transformation since this is a constraint",
+            warning_error(wi,
+                          gmx::formatString("pull-coord%d can not use pull-coord%d in the "
+                                            "transformation since this is a "
+                                            "constraint",
+                                            coord_index_for_output,
+                                            previousCoordOutputIndex));
+        }
+        else if (previousPcrd.k != 0 && pcrd->k != 0)
+        {
+            warning_note(
+                    wi,
+                    gmx::formatString("pull-coord%d has a non-zero force constant and is also "
+                                      "referenced in pull-coord%d-expression. "
+                                      "Make sure that this is intended. "
+                                      "In most use cases, the pull coordinates referenced by a "
+                                      "transformation coordinate should have their force constant "
+                                      "set to zero.",
                                       coord_index_for_output,
-                                      previousCoordOutputIndex)));
+                                      previousCoordOutputIndex));
         }
     }
 }
@@ -365,7 +381,7 @@ static void init_pull_coord(t_pull_coord*        pcrd,
     }
     if (pcrd->eGeom == PullGroupGeometry::Transformation)
     {
-        initTransformationPullCoord(pcrd, pull);
+        initTransformationPullCoord(pcrd, pull, wi);
     }
 
     for (m = 0; m < DIM; m++)

diff --git a/src/gromacs/gmxpreprocess/tests/readir.cpp b/src/gromacs/gmxpreprocess/tests/readir.cpp
index de21db7179f34cf84051af8ef9743f8dfa00fe58..586c95f99ccd5991d23124007678a93494c83769 100644 (file)
--- a/src/gromacs/gmxpreprocess/tests/readir.cpp
+++ b/src/gromacs/gmxpreprocess/tests/readir.cpp
@@ -86,8 +86,9 @@ public:
     //! Tells whether warnings and/or errors are expected from inputrec parsing and checking, and whether we should compare the output
     enum class TestBehavior
     {
-        NoErrorAndCompareOutput,   //!< Expect no warnings/error and compare output
-        ErrorAndCompareOutput,     //!< Expect at least one warning/error and compare output
+        NoErrorAndCompareOutput,      //!< Expect no warnings/error and compare output
+        NoErrorAndDoNotCompareOutput, //!< Expect no warnings/error and do not compare output
+        ErrorAndCompareOutput,        //!< Expect at least one warning/error and compare output
         ErrorAndDoNotCompareOutput //!< Expect at least one warning/error and do not compare output
     };
 
@@ -99,8 +100,10 @@ public:
     void runTest(const std::string& inputMdpFileContents,
                  const TestBehavior testBehavior = TestBehavior::NoErrorAndCompareOutput)
     {
-        const bool expectError   = testBehavior != TestBehavior::NoErrorAndCompareOutput;
-        const bool compareOutput = testBehavior != TestBehavior::ErrorAndDoNotCompareOutput;
+        const bool expectError = testBehavior == TestBehavior::ErrorAndCompareOutput
+                                 || testBehavior == TestBehavior::ErrorAndDoNotCompareOutput;
+        const bool compareOutput = testBehavior == TestBehavior::ErrorAndCompareOutput
+                                   || testBehavior == TestBehavior::NoErrorAndCompareOutput;
 
         std::string inputMdpFilename = fileManager_.getTemporaryFilePath("input.mdp");
         std::string outputMdpFilename;
@@ -108,7 +111,6 @@ public:
         {
             outputMdpFilename = fileManager_.getTemporaryFilePath("output.mdp");
         }
-
         TextWriter::writeFileFromString(inputMdpFilename, inputMdpFileContents);
 
         get_ir(inputMdpFilename.c_str(),
@@ -281,5 +283,73 @@ TEST_F(GetIrTest, AcceptsMimic)
     runTest(joinStrings(inputMdpFile, "\n"));
 }
 
+#if HAVE_MUPARSER
+
+TEST_F(GetIrTest, AcceptsTransformationCoord)
+{
+    const char* inputMdpFile[] = {
+        "pull = yes",
+        "pull-ngroups = 2",
+        "pull-ncoords = 2",
+        "pull-coord1-geometry = distance",
+        "pull-coord1-groups = 1 2",
+        "pull-coord1-k = 1",
+        "pull-coord2-geometry = transformation",
+        "pull-coord2-expression = 1/(1+x1)",
+        "pull-coord2-k = 10",
+    };
+    runTest(joinStrings(inputMdpFile, "\n"), TestBehavior::NoErrorAndDoNotCompareOutput);
+}
+
+TEST_F(GetIrTest, InvalidTransformationCoordWithConstraint)
+{
+    const char* inputMdpFile[] = {
+        "pull = yes",
+        "pull-ncoords = 1",
+        "pull-coord1-geometry = transformation",
+        "pull-coord1-type = constraint", // INVALID
+        "pull-coord1-expression = 10",
+    };
+    runTest(joinStrings(inputMdpFile, "\n"), TestBehavior::ErrorAndDoNotCompareOutput);
+}
+
+TEST_F(GetIrTest, InvalidPullCoordWithConstraintInTransformationExpression)
+{
+    const char* inputMdpFile[] = {
+        "pull = yes",
+        "pull-ngroups = 2",
+        "pull-ncoords = 2",
+        "pull-coord1-geometry = distance",
+        "pull-coord1-type = constraint", // INVALID
+        "pull-coord1-groups = 1 2",
+        "pull-coord2-geometry = transformation",
+        "pull-coord2-expression = x1",
+    };
+    runTest(joinStrings(inputMdpFile, "\n"), TestBehavior::ErrorAndDoNotCompareOutput);
+}
+
+TEST_F(GetIrTest, InvalidTransformationCoordDxValue)
+{
+    const char* inputMdpFile[] = {
+        "pull = yes",
+        "pull-ncoords = 1",
+        "pull-coord1-geometry = transformation",
+        "pull-coord1-expression = 10",
+        "pull-coord1-dx = 0", // INVALID
+    };
+    runTest(joinStrings(inputMdpFile, "\n"), TestBehavior::ErrorAndDoNotCompareOutput);
+}
+
+TEST_F(GetIrTest, MissingTransformationCoordExpression)
+{
+    const char* inputMdpFile[] = {
+        "pull = yes",
+        "pull-ncoords = 1",
+        "pull-coord1-geometry = transformation",
+    };
+    runTest(joinStrings(inputMdpFile, "\n"), TestBehavior::ErrorAndDoNotCompareOutput);
+}
+#endif // HAVE_MUPARSER
+
 } // namespace test
 } // namespace gmx

diff --git a/src/programs/mdrun/tests/grompp.cpp b/src/programs/mdrun/tests/grompp.cpp
index ab3003fbde0869530655f8cd21d4a2819bf216d3..de7d8598eba7b196b5468ebfed72b32b682bbffe 100644 (file)
--- a/src/programs/mdrun/tests/grompp.cpp
+++ b/src/programs/mdrun/tests/grompp.cpp
@@ -46,10 +46,15 @@
  */
 #include "gmxpre.h"
 
+#include "gromacs/utility/stringutil.h"
+#include "gromacs/gmxpreprocess/readir.h"
+
 #include <gtest/gtest.h>
 
 #include "moduletest.h"
 
+#include "testutils/testexceptions.h"
+
 namespace
 {
 
@@ -96,5 +101,43 @@ TEST_F(GromppTest, SimulatedAnnealingWorksWithMultipleGroups)
     runTest();
 }
 
+#if HAVE_MUPARSER
+
+TEST_F(GromppTest, ValidTransformationCoord)
+{
+    const char* inputMdpFile[] = {
+        "pull = yes",
+        "pull-ncoords = 2",
+        "pull-ngroups = 2",
+        "pull-group1-name = SOL",
+        "pull-group2-name = Methanol",
+        "pull-coord1-geometry = distance",
+        "pull-coord1-groups = 1 2",
+        "pull-coord2-geometry = transformation",
+        "pull-coord2-expression = x1", // Valid expression
+    };
+    runner_.useStringAsMdpFile(gmx::joinStrings(inputMdpFile, "\n"));
+    runTest();
+}
+
+TEST_F(GromppTest, InvalidTransformationCoord)
+{
+    const char* inputMdpFile[] = {
+        "pull = yes",
+        "pull-ncoords = 2",
+        "pull-ngroups = 2",
+        "pull-group1-name = SOL",
+        "pull-group2-name = Methanol",
+        "pull-coord1-geometry = distance",
+        "pull-coord1-groups = 1 2",
+        "pull-coord2-geometry = transformation",
+        "pull-coord2-expression = x2", // Invalid expression -> evaluation should fail
+    };
+    runner_.useStringAsMdpFile(gmx::joinStrings(inputMdpFile, "\n"));
+    ASSERT_THROW(runTest(), gmx::InconsistentInputError);
+    done_inputrec_strings(); // This allows grompp to be called again in another test
+}
+#endif // HAVE_MUPARSER
+
 
 } // namespace

diff --git a/src/programs/mdrun/tests/pull.cpp b/src/programs/mdrun/tests/pull.cpp
index 48665c68e0ff2024cd5d952e8dba8ff92a0c5323..41c692a5b0f26d2369e170a4e59e6354117b062d 100644 (file)
--- a/src/programs/mdrun/tests/pull.cpp
+++ b/src/programs/mdrun/tests/pull.cpp
@@ -47,6 +47,7 @@
 #include <string>
 #include <unordered_map>
 #include <vector>
+#include <regex>
 
 #include <gtest/gtest.h>
 
@@ -71,14 +72,10 @@ namespace
 {
 
 /*! \brief Database of energy tolerances on the various systems. */
-std::unordered_map<std::string, FloatingPointTolerance> energyToleranceForSystem_g = {
-    { { "spc216", relativeToleranceAsFloatingPoint(1, 1e-4) } }
-};
+std::unordered_map<std::string, double> energyToleranceForSystem_g = { { { "spc216", 1e-4 } } };
 
 /*! \brief Database of pressure tolerances on the various systems. */
-std::unordered_map<std::string, FloatingPointTolerance> pressureToleranceForSystem_g = {
-    { { "spc216", relativeToleranceAsFloatingPoint(1, 2e-4) } }
-};
+std::unordered_map<std::string, double> pressureToleranceForSystem_g = { { { "spc216", 2e-4 } } };
 
 const std::unordered_map<std::string, std::vector<std::pair<std::string, std::string>>> c_mdpPullParams = {
     { { "umbrella-3D",
@@ -125,7 +122,26 @@ const std::unordered_map<std::string, std::vector<std::pair<std::string, std::st
           { "pull-coord2-geometry", "distance" },
           { "pull-coord2-dim", "Y Y Y" },
           { "pull-coord2-init", "0.4" },
-          { "pull-coord2-k", "100" } } } }
+          { "pull-coord2-k", "100" } } },
+      { "transformation-coord-umbrella-3D", // should yield the same force on the first pull coordinate as in the umbrella-3D test above
+        {
+                { "pull-ngroups", "2" },
+                { "pull-ncoords", "2" },
+                { "pull-nstxout", "0" },
+                { "pull-nstfout", "0" },
+                { "pull-group1-name", "r_1" },
+                { "pull-group2-name", "r_2" },
+                { "pull-coord1-groups", "1 2" },
+                { "pull-coord1-type", "umbrella" },
+                { "pull-coord1-geometry", "distance" },
+                { "pull-coord1-dim", "Y Y Y" },
+                { "pull-coord1-k", "0" }, // set to zero since the force comes from the transformation coordinate
+                { "pull-coord2-geometry", "transformation" },
+                { "pull-coord2-type", "umbrella" },
+                { "pull-coord2-expression", "2*(0.6 - x1)" },
+                { "pull-coord2-init", "0" }, // -> pull-coord1-init = 0.6
+                { "pull-coord2-k", "25" }    // -> pull-coord1-k = 2^2*25 = 100
+        } } }
 };
 
 //! Helper type
@@ -155,6 +171,14 @@ void addBasicMdpValues(MdpFieldValues* mdpFieldValues)
     (*mdpFieldValues)["vdwtype"]       = "Cut-off";
 }
 
+bool isTransformationPullSetup(const std::string& pullSetupName)
+{
+    const auto* pattern{ "transformation" };
+    // std::regex_constants::icase - TO IGNORE CASE.
+    auto rx = std::regex{ pattern, std::regex_constants::icase };
+    return std::regex_search(pullSetupName, rx);
+}
+
 TEST_P(PullIntegrationTest, WithinTolerances)
 {
     auto params         = GetParam();
@@ -195,13 +219,23 @@ TEST_P(PullIntegrationTest, WithinTolerances)
     }
     // Do mdrun
     {
+        // conver the tolerance to relative floating point tolerance
+        auto energyTolerance   = energyToleranceForSystem_g.at(simulationName);
+        auto pressureTolerance = pressureToleranceForSystem_g.at(simulationName);
+        if (isTransformationPullSetup(pullSetup)) // need to increase the tolerance a bit due to numerical evaluations
+        {
+            energyTolerance *= 10;
+            pressureTolerance *= 10;
+        }
+        auto relativeEnergyTolerance   = relativeToleranceAsFloatingPoint(1, energyTolerance);
+        auto relativePressureTolerance = relativeToleranceAsFloatingPoint(1, pressureTolerance);
         CommandLine mdrunCaller;
         ASSERT_EQ(0, runner_.callMdrun(mdrunCaller));
         EnergyTermsToCompare energyTermsToCompare{ {
-                { interaction_function[F_COM_PULL].longname, energyToleranceForSystem_g.at(simulationName) },
-                { interaction_function[F_EPOT].longname, energyToleranceForSystem_g.at(simulationName) },
-                { interaction_function[F_EKIN].longname, energyToleranceForSystem_g.at(simulationName) },
-                { interaction_function[F_PRES].longname, pressureToleranceForSystem_g.at(simulationName) },
+                { interaction_function[F_COM_PULL].longname, relativeEnergyTolerance },
+                { interaction_function[F_EPOT].longname, relativeEnergyTolerance },
+                { interaction_function[F_EKIN].longname, relativeEnergyTolerance },
+                { interaction_function[F_PRES].longname, relativePressureTolerance },
         } };
         TestReferenceData    refData;
         auto                 checker = refData.rootChecker()
@@ -216,7 +250,8 @@ INSTANTIATE_TEST_SUITE_P(PullTest,
                          ::testing::Combine(::testing::Values("spc216"),
                                             ::testing::Values("umbrella-3D",
                                                               "umbrella-2D",
-                                                              "constraint-flatbottom")));
+                                                              "constraint-flatbottom",
+                                                              "transformation-coord-umbrella-3D")));
 
 } // namespace
 } // namespace test

diff --git a/src/programs/mdrun/tests/refdata/PullTest_PullIntegrationTest_WithinTolerances_3.xml b/src/programs/mdrun/tests/refdata/PullTest_PullIntegrationTest_WithinTolerances_3.xml
new file mode 100644 (file)
index 0000000..cc4df96
--- /dev/null
+++ b/src/programs/mdrun/tests/refdata/PullTest_PullIntegrationTest_WithinTolerances_3.xml
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <Simulation Name="spc216">
+    <PullSetup Name="transformation-coord-umbrella-3D">
+      <Energy Name="Pressure">
+        <Real Name="Time 0.000000 Step 0 in frame 0">7047.8843</Real>
+        <Real Name="Time 0.005000 Step 5 in frame 1">6813.9585</Real>
+        <Real Name="Time 0.010000 Step 10 in frame 2">6432.3613</Real>
+        <Real Name="Time 0.015000 Step 15 in frame 3">5985.5601</Real>
+        <Real Name="Time 0.020000 Step 20 in frame 4">5571.9185</Real>
+      </Energy>
+      <Energy Name="Kinetic En.">
+        <Real Name="Time 0.000000 Step 0 in frame 0">1627.3685</Real>
+        <Real Name="Time 0.005000 Step 5 in frame 1">1620.8447</Real>
+        <Real Name="Time 0.010000 Step 10 in frame 2">1637.7889</Real>
+        <Real Name="Time 0.015000 Step 15 in frame 3">1665.6293</Real>
+        <Real Name="Time 0.020000 Step 20 in frame 4">1688.8979</Real>
+      </Energy>
+      <Energy Name="Potential">
+        <Real Name="Time 0.000000 Step 0 in frame 0">-7927.1426</Real>
+        <Real Name="Time 0.005000 Step 5 in frame 1">-7920.6416</Real>
+        <Real Name="Time 0.010000 Step 10 in frame 2">-7937.6411</Real>
+        <Real Name="Time 0.015000 Step 15 in frame 3">-7965.5195</Real>
+        <Real Name="Time 0.020000 Step 20 in frame 4">-7988.6768</Real>
+      </Energy>
+      <Energy Name="COM Pull En.">
+        <Real Name="Time 0.000000 Step 0 in frame 0">0.031361848</Real>
+        <Real Name="Time 0.005000 Step 5 in frame 1">0.039091855</Real>
+        <Real Name="Time 0.010000 Step 10 in frame 2">0.047064677</Real>
+        <Real Name="Time 0.015000 Step 15 in frame 3">0.055442076</Real>
+        <Real Name="Time 0.020000 Step 20 in frame 4">0.064478263</Real>
+      </Energy>
+    </PullSetup>
+  </Simulation>
+</ReferenceData>