src/programs/mdrun/tests/freezegroups.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  * \brief End-to-end tests checking sanity of results of simulations
  38  *        containing freeze groups
  39  *
  40  * \author Pascal Merz <pascal.merz@me.com>
  41  * \ingroup module_mdrun_integration_tests
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "config.h"
  46
  47 #include "gromacs/topology/ifunc.h"
  48 #include "gromacs/utility/stringutil.h"
  49
  50 #include "testutils/mpitest.h"
  51 #include "testutils/simulationdatabase.h"
  52 #include "testutils/testmatchers.h"
  53
  54 #include "moduletest.h"
  55 #include "simulatorcomparison.h"
  56 #include "trajectoryreader.h"
  57
  58 namespace gmx::test
  59 {
  60 namespace
  61 {
  62 /*! \brief Test fixture checking sanity of freeze group results
  63  *
  64  * This tests the sanity of simulation results containing fully and partially
  65  * frozen atoms. For fully frozen atoms, it checks that their reported position
  66  * is identical for all steps, and that their velocity is zero. For partially
  67  * frozen atoms (for simplicity only in z-direction), it checks that their
  68  * position is identical in the frozen dimension for all steps, and that their
  69  * velocity is zero in the frozen dimension.
  70  */
  71 using FreezeGroupTestParams = std::tuple<std::string, std::string, std::string>;
  72 class FreezeGroupTest : public MdrunTestFixture, public ::testing::WithParamInterface<FreezeGroupTestParams>
  73 {
  74 public:
  75     //! Check that the frozen positions don't change and velocities are zero
  76     static void checkFreezeGroups(const std::string&           trajectoryName,
  77                                   ArrayRef<const unsigned int> fullyFrozenAtoms,
  78                                   ArrayRef<const unsigned int> partiallyFrozenAtomsDimZ,
  79                                   const TrajectoryTolerances&  tolerances)
  80     {
  81         auto [fullyFrozenPositions, fullyFrozenVelocities] =
  82                 getFrozenPositionsAndVelocities(trajectoryName, fullyFrozenAtoms);
  83         auto [partiallyFrozenPositions, partiallyFrozenVelocities] =
  84                 getFrozenPositionsAndVelocities(trajectoryName, partiallyFrozenAtomsDimZ);
  85         GMX_RELEASE_ASSERT(fullyFrozenPositions.size() == fullyFrozenVelocities.size(),
  86                            "Position and velocity trajectory don't have the same length.");
  87         GMX_RELEASE_ASSERT(partiallyFrozenPositions.size() == partiallyFrozenVelocities.size(),
  88                            "Position and velocity trajectory don't have the same length.");
  89         GMX_RELEASE_ASSERT(fullyFrozenPositions.size() == partiallyFrozenPositions.size(),
  90                            "Fully and partially frozen trajectory don't have the same length.");
  91         const auto trajectorySize = fullyFrozenPositions.size();
  92
  93         for (auto frameIdx = decltype(trajectorySize){ 0 }; frameIdx < trajectorySize; frameIdx++)
  94         {
  95             SCOPED_TRACE(formatString("Checking frame %lu", frameIdx + 1));
  96             if (frameIdx > 0)
  97             {
  98                 checkFullyFrozenPositions(
  99                         fullyFrozenPositions[frameIdx], fullyFrozenPositions[frameIdx - 1], tolerances);
 100                 checkZDimFrozenPositions(partiallyFrozenPositions[frameIdx],
 101                                          partiallyFrozenPositions[frameIdx - 1],
 102                                          tolerances);
 103             }
 104             checkFullyFrozenVelocities(fullyFrozenVelocities[frameIdx], tolerances);
 105             checkZDimFrozenVelocities(partiallyFrozenVelocities[frameIdx], tolerances);
 106         }
 107     }
 108
 109     //! Check that fully frozen frame velocities are zero
 110     static void checkFullyFrozenVelocities(ArrayRef<const RVec>        velocities,
 111                                            const TrajectoryTolerances& tolerances)
 112     {
 113         SCOPED_TRACE("Checking fully frozen velocity frame");
 114         std::vector<RVec> zeroVelocities(velocities.size(), RVec{ 0, 0, 0 });
 115         EXPECT_THAT(zeroVelocities, Pointwise(RVecEq(tolerances.velocities), velocities));
 116     }
 117     //! Check that z-dimension frozen frame velocities are zero
 118     static void checkZDimFrozenVelocities(ArrayRef<const RVec>        velocities,
 119                                           const TrajectoryTolerances& tolerances)
 120     {
 121         SCOPED_TRACE("Checking z-dimension frozen velocity frame");
 122         std::vector<real> zVelocities;
 123         for (const auto& v : velocities)
 124         {
 125             zVelocities.emplace_back(v[ZZ]);
 126         }
 127         std::vector<real> zeroVelocities(zVelocities.size(), 0);
 128         EXPECT_THAT(zeroVelocities, Pointwise(RealEq(tolerances.velocities), zVelocities));
 129     }
 130     //! Check that fully frozen frame positions are static
 131     static void checkFullyFrozenPositions(ArrayRef<const RVec>        positions,
 132                                           ArrayRef<const RVec>        previousPositions,
 133                                           const TrajectoryTolerances& tolerances)
 134     {
 135         SCOPED_TRACE("Checking fully frozen position frame");
 136         EXPECT_THAT(previousPositions, Pointwise(RVecEq(tolerances.coordinates), positions));
 137     }
 138     //! Check that z-dimension frozen frame positions are zero
 139     static void checkZDimFrozenPositions(ArrayRef<const RVec>        positions,
 140                                          ArrayRef<const RVec>        previousPositions,
 141                                          const TrajectoryTolerances& tolerances)
 142     {
 143         SCOPED_TRACE("Checking z-dimension frozen position frame");
 144         std::vector<real> zPositions;
 145         for (const auto& p : positions)
 146         {
 147             zPositions.emplace_back(p[ZZ]);
 148         }
 149         std::vector<real> zPrevPositions;
 150         for (const auto& p : previousPositions)
 151         {
 152             zPrevPositions.emplace_back(p[ZZ]);
 153         }
 154         EXPECT_THAT(zPrevPositions, Pointwise(RealEq(tolerances.coordinates), zPositions));
 155     }
 156
 157     static std::tuple<std::vector<std::vector<RVec>>, std::vector<std::vector<RVec>>>
 158     getFrozenPositionsAndVelocities(const std::string& trajectoryName, ArrayRef<const unsigned int> frozenAtoms)
 159     {
 160         std::vector<std::vector<RVec>> positions;
 161         std::vector<std::vector<RVec>> velocities;
 162
 163         TrajectoryFrameReader trajectoryFrameReader(trajectoryName);
 164         while (trajectoryFrameReader.readNextFrame())
 165         {
 166             const auto frame = trajectoryFrameReader.frame();
 167             positions.emplace_back();
 168             velocities.emplace_back();
 169             for (const auto& index : frozenAtoms)
 170             {
 171                 positions.back().emplace_back(frame.x().at(index));
 172                 velocities.back().emplace_back(frame.v().at(index));
 173             }
 174         }
 175
 176         return { std::move(positions), std::move(velocities) };
 177     }
 178 };
 179
 180 TEST_P(FreezeGroupTest, WithinTolerances)
 181 {
 182     const auto& params         = GetParam();
 183     const auto& integrator     = std::get<0>(params);
 184     const auto& tcoupling      = std::get<1>(params);
 185     const auto& pcoupling      = std::get<2>(params);
 186     const auto& simulationName = "alanine_vacuo";
 187
 188     constexpr std::array<unsigned int, 5>  backbone   = { 4, 6, 8, 14, 16 };
 189     constexpr std::array<unsigned int, 13> sideChainH = { 0,  1,  2,  3,  9,  10, 11,
 190                                                           12, 13, 18, 19, 20, 21 };
 191
 192     if (integrator == "md-vv" && pcoupling == "parrinello-rahman")
 193     {
 194         // Parrinello-Rahman is not implemented in md-vv
 195         return;
 196     }
 197
 198     // Prepare mdp input
 199     auto mdpFieldValues = prepareMdpFieldValues(simulationName, integrator, tcoupling, pcoupling);
 200     mdpFieldValues["nsteps"]      = "8";
 201     mdpFieldValues["nstxout"]     = "4";
 202     mdpFieldValues["nstvout"]     = "4";
 203     mdpFieldValues["freezegrps"]  = "Backbone SideChain";
 204     mdpFieldValues["freezedim"]   = "Y Y Y N N Y";
 205     mdpFieldValues["constraints"] = "all-bonds";
 206
 207     // Run grompp
 208     runner_.useTopGroAndNdxFromDatabase(simulationName);
 209     runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
 210     // Allow one warning for COMM removal + partially frozen atoms
 211     runGrompp(&runner_, { SimulationOptionTuple("-maxwarn", "1") });
 212     // Run mdrun
 213     runMdrun(&runner_);
 214
 215     // Check frozen atoms
 216     checkFreezeGroups(runner_.fullPrecisionTrajectoryFileName_,
 217                       backbone,
 218                       sideChainH,
 219                       TrajectoryComparison::s_defaultTrajectoryTolerances);
 220 }
 221
 222 INSTANTIATE_TEST_SUITE_P(FreezeWorks,
 223                          FreezeGroupTest,
 224                          ::testing::Combine(::testing::Values("md", "md-vv", "sd", "bd"),
 225                                             ::testing::Values("no"),
 226                                             ::testing::Values("no")));
 227 } // namespace
 228 } // namespace gmx::test