src/gromacs/mdtypes/multipletimestepping.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020,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 #include "gmxpre.h"
  36
  37 #include "multipletimestepping.h"
  38
  39 #include <optional>
  40
  41 #include "gromacs/mdtypes/inputrec.h"
  42 #include "gromacs/mdtypes/pull_params.h"
  43 #include "gromacs/utility/arrayref.h"
  44 #include "gromacs/utility/gmxassert.h"
  45 #include "gromacs/utility/stringutil.h"
  46
  47 namespace gmx
  48 {
  49
  50 int nonbondedMtsFactor(const t_inputrec& ir)
  51 {
  52     GMX_RELEASE_ASSERT(!ir.useMts || ir.mtsLevels.size() == 2, "Only 2 MTS levels supported here");
  53
  54     if (ir.useMts && ir.mtsLevels[1].forceGroups[static_cast<int>(MtsForceGroups::Nonbonded)])
  55     {
  56         return ir.mtsLevels[1].stepFactor;
  57     }
  58     else
  59     {
  60         return 1;
  61     }
  62 }
  63
  64 std::vector<MtsLevel> setupMtsLevels(const GromppMtsOpts& mtsOpts, std::vector<std::string>* errorMessages)
  65 {
  66     std::vector<MtsLevel> mtsLevels;
  67
  68     if (mtsOpts.numLevels != 2)
  69     {
  70         if (errorMessages)
  71         {
  72             errorMessages->push_back("Only mts-levels = 2 is supported");
  73         }
  74     }
  75     else
  76     {
  77         mtsLevels.resize(2);
  78
  79         const std::vector<std::string> inputForceGroups = gmx::splitString(mtsOpts.level2Forces);
  80         auto&                          forceGroups      = mtsLevels[1].forceGroups;
  81         for (const auto& inputForceGroup : inputForceGroups)
  82         {
  83             bool found     = false;
  84             int  nameIndex = 0;
  85             for (const auto& forceGroupName : gmx::mtsForceGroupNames)
  86             {
  87                 if (gmx::equalCaseInsensitive(inputForceGroup, forceGroupName))
  88                 {
  89                     forceGroups.set(nameIndex);
  90                     found = true;
  91                 }
  92                 nameIndex++;
  93             }
  94             if (!found && errorMessages)
  95             {
  96                 errorMessages->push_back(
  97                         gmx::formatString("Unknown MTS force group '%s'", inputForceGroup.c_str()));
  98             }
  99         }
 100
 101         // Make the level 0 use the complement of the force groups of group 1
 102         mtsLevels[0].forceGroups = ~mtsLevels[1].forceGroups;
 103         mtsLevels[0].stepFactor  = 1;
 104
 105         mtsLevels[1].stepFactor = mtsOpts.level2Factor;
 106
 107         if (errorMessages && mtsLevels[1].stepFactor <= 1)
 108         {
 109             errorMessages->push_back("mts-factor should be larger than 1");
 110         }
 111     }
 112
 113     return mtsLevels;
 114 }
 115
 116 bool haveValidMtsSetup(const t_inputrec& ir)
 117 {
 118     return (ir.useMts && ir.mtsLevels.size() == 2 && ir.mtsLevels[1].stepFactor > 1);
 119 }
 120
 121 namespace
 122 {
 123
 124 //! Checks whether \p nstValue is a multiple of the largest MTS step, returns an error string for parameter \p param when this is not the case
 125 std::optional<std::string> checkMtsInterval(ArrayRef<const MtsLevel> mtsLevels, const char* param, const int nstValue)
 126 {
 127     GMX_RELEASE_ASSERT(mtsLevels.size() >= 2, "Need at least two levels for MTS");
 128
 129     const int mtsFactor = mtsLevels.back().stepFactor;
 130     if (nstValue % mtsFactor == 0)
 131     {
 132         return {};
 133     }
 134     else
 135     {
 136         return gmx::formatString(
 137                 "With MTS, %s = %d should be a multiple of mts-factor = %d", param, nstValue, mtsFactor);
 138     }
 139 }
 140
 141 } // namespace
 142
 143 std::vector<std::string> checkMtsRequirements(const t_inputrec& ir)
 144 {
 145     std::vector<std::string> errorMessages;
 146
 147     if (!ir.useMts)
 148     {
 149         return errorMessages;
 150     }
 151
 152     GMX_RELEASE_ASSERT(haveValidMtsSetup(ir), "MTS setup should be valid here");
 153
 154     ArrayRef<const MtsLevel> mtsLevels = ir.mtsLevels;
 155
 156     if (!(ir.eI == IntegrationAlgorithm::MD || ir.eI == IntegrationAlgorithm::SD1))
 157     {
 158         errorMessages.push_back(gmx::formatString(
 159                 "Multiple time stepping is only supported with integrators %s and %s",
 160                 enumValueToString(IntegrationAlgorithm::MD),
 161                 enumValueToString(IntegrationAlgorithm::SD1)));
 162     }
 163
 164     if ((EEL_FULL(ir.coulombtype) || EVDW_PME(ir.vdwtype))
 165         && forceGroupMtsLevel(ir.mtsLevels, MtsForceGroups::LongrangeNonbonded) == 0)
 166     {
 167         errorMessages.emplace_back(
 168                 "With long-range electrostatics and/or LJ treatment, the long-range part "
 169                 "has to be part of the mts-level2-forces");
 170     }
 171
 172     std::optional<std::string> mesg;
 173     if (ir.nstcalcenergy > 0)
 174     {
 175         if ((mesg = checkMtsInterval(mtsLevels, "nstcalcenergy", ir.nstcalcenergy)))
 176         {
 177             errorMessages.push_back(mesg.value());
 178         }
 179     }
 180     if ((mesg = checkMtsInterval(mtsLevels, "nstenergy", ir.nstenergy)))
 181     {
 182         errorMessages.push_back(mesg.value());
 183     }
 184     if ((mesg = checkMtsInterval(mtsLevels, "nstlog", ir.nstlog)))
 185     {
 186         errorMessages.push_back(mesg.value());
 187     }
 188     if ((mesg = checkMtsInterval(mtsLevels, "nstfout", ir.nstfout)))
 189     {
 190         errorMessages.push_back(mesg.value());
 191     }
 192     if (ir.efep != FreeEnergyPerturbationType::No)
 193     {
 194         if ((mesg = checkMtsInterval(mtsLevels, "nstdhdl", ir.fepvals->nstdhdl)))
 195         {
 196             errorMessages.push_back(mesg.value());
 197         }
 198     }
 199     if (mtsLevels.back().forceGroups[static_cast<int>(gmx::MtsForceGroups::Nonbonded)])
 200     {
 201         if ((mesg = checkMtsInterval(mtsLevels, "nstlist", ir.nstlist)))
 202         {
 203             errorMessages.push_back(mesg.value());
 204         }
 205     }
 206
 207     if (ir.bPull)
 208     {
 209         const int pullMtsLevel  = gmx::forceGroupMtsLevel(ir.mtsLevels, gmx::MtsForceGroups::Pull);
 210         const int mtsStepFactor = ir.mtsLevels[pullMtsLevel].stepFactor;
 211         if (ir.pull->nstxout % mtsStepFactor != 0)
 212         {
 213             errorMessages.emplace_back("pull-nstxout should be a multiple of mts-factor");
 214         }
 215         if (ir.pull->nstfout % mtsStepFactor != 0)
 216         {
 217             errorMessages.emplace_back("pull-nstfout should be a multiple of mts-factor");
 218         }
 219     }
 220
 221     return errorMessages;
 222 }
 223
 224 } // namespace gmx