From: Oliver Fleetwood Date: Wed, 27 Oct 2021 14:37:29 +0000 (+0000) Subject: Add preprocessing and integration test for transformation coordinate X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?p=alexxy%2Fgromacs.git;a=commitdiff_plain;h=2acfc914079ae665c4f2ce16ec951565eaebbbff 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 2524976cb5..fed4eabb6d 100644 --- 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 de21db7179..586c95f99c 100644 --- 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 ab3003fbde..de7d8598eb 100644 --- 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 #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 48665c68e0..41c692a5b0 100644 --- a/src/programs/mdrun/tests/pull.cpp +++ b/src/programs/mdrun/tests/pull.cpp @@ -47,6 +47,7 @@ #include #include #include +#include #include @@ -71,14 +72,10 @@ namespace { /*! \brief Database of energy tolerances on the various systems. */ -std::unordered_map energyToleranceForSystem_g = { - { { "spc216", relativeToleranceAsFloatingPoint(1, 1e-4) } } -}; +std::unordered_map energyToleranceForSystem_g = { { { "spc216", 1e-4 } } }; /*! \brief Database of pressure tolerances on the various systems. */ -std::unordered_map pressureToleranceForSystem_g = { - { { "spc216", relativeToleranceAsFloatingPoint(1, 2e-4) } } -}; +std::unordered_map pressureToleranceForSystem_g = { { { "spc216", 2e-4 } } }; const std::unordered_map>> c_mdpPullParams = { { { "umbrella-3D", @@ -125,7 +122,26 @@ const std::unordered_map 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 index 0000000000..cc4df966ff --- /dev/null +++ b/src/programs/mdrun/tests/refdata/PullTest_PullIntegrationTest_WithinTolerances_3.xml @@ -0,0 +1,36 @@ + + + + + + + 7047.8843 + 6813.9585 + 6432.3613 + 5985.5601 + 5571.9185 + + + 1627.3685 + 1620.8447 + 1637.7889 + 1665.6293 + 1688.8979 + + + -7927.1426 + -7920.6416 + -7937.6411 + -7965.5195 + -7988.6768 + + + 0.031361848 + 0.039091855 + 0.047064677 + 0.055442076 + 0.064478263 + + + +