src/testutils/simulationdatabase.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2016,2018,2019,2020, 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 /*! \internal \file
  36  * \brief
  37  * Implements declarations from in simulationdatabase.h
  38  *
  39  * \author Mark Abraham <mark.j.abraham@gmail.com>
  40  * \ingroup module_testutils
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "testutils/simulationdatabase.h"
  45
  46 #include <algorithm>
  47 #include <map>
  48 #include <string>
  49
  50 #include "gromacs/utility/stringutil.h"
  51
  52 #include "testutils/testasserts.h"
  53
  54 namespace gmx
  55 {
  56 namespace test
  57 {
  58
  59 namespace
  60 {
  61
  62 struct DatabaseEntry
  63 {
  64     MdpFieldValues   mdpFieldValues;
  65     std::vector<int> validPpRankCounts;
  66 };
  67
  68 //! Helper typedef
  69 using MdpFileValues = std::map<std::string, DatabaseEntry>;
  70
  71 //! Database of .mdp strings that supports prepareDefaultMdpValues()
  72 const MdpFileValues mdpFileValueDatabase_g{
  73     // Simple system with 12 argon atoms, fairly widely separated
  74     { "argon12",
  75       { { { "ref-t", "80" }, { "compressibility", "5e-10" }, { "tau-p", "1000" } }, { 1, 2, 3, 4 } } },
  76     // Simple system with 5 water molecules, fairly widely separated
  77     { "tip3p5", { { { "compressibility", "5e-10" }, { "tau-p", "1000" } }, { 1, 2, 3, 4, 5, 6, 8, 9 } } },
  78     // Simple system with 5832 argon atoms, suitable for normal pressure coupling
  79     { "argon5832",
  80       { { { "ref-t", "80" } },
  81         { // TODO This test case is not currently used, so we
  82           // have not tested which rank counts work.
  83           1,
  84           2,
  85           3,
  86           4,
  87           5,
  88           6,
  89           7,
  90           8,
  91           9 } } },
  92     // Simple system with 2 nearby water molecules
  93     { "spc2",
  94       { {},
  95         { // TODO This test case is not currently used, so we
  96           // have not tested which rank counts work.
  97           1,
  98           2,
  99           3,
 100           4,
 101           5,
 102           6,
 103           7,
 104           8,
 105           9 } } },
 106     // Simple system with 216 water molecules, condensed phase
 107     { "spc216",
 108       { {},
 109         {
 110                 // TODO This test case is not currently used, so we
 111                 // have not tested which rank counts work.
 112                 1,
 113                 2,
 114                 3,
 115                 4,
 116                 5,
 117                 6,
 118                 7,
 119                 8,
 120                 9 // TODO tpi test
 121         } } },
 122     // Capped alanine peptide in vacuo with virtual sites
 123     { "alanine_vsite_vacuo",
 124       { { { "constraints", "all-bonds" }, { "compressibility", "5e-10" }, { "tau-p", "1000" } },
 125         { 1, 2, 3, 4, 6, 9 } } },
 126     // Capped alanine peptide in aqueous condensed phase, with virtual sites
 127     { "alanine_vsite_solvated",
 128       { { { "constraints", "all-bonds" }, { "compressibility", "5e-10" }, { "tau-p", "1000" } },
 129         { // TODO This test case is not currently used, so we
 130           // have not tested which rank counts work.
 131           1,
 132           2,
 133           3,
 134           4,
 135           5,
 136           6,
 137           7,
 138           8,
 139           9 } } },
 140     // Zwitterionic glycine in vacuo
 141     { "glycine_vacuo", { { { "constraints", "h-bonds" } }, { 1, 2, 3, 4, 5, 6, 7, 8, 9 } } },
 142     // Zwitterionic glycine in vacuo, without constraints
 143     { "glycine_no_constraints_vacuo", { { { "constraints", "none" } }, { 1, 2, 3, 4, 5, 6, 7, 8, 9 } } },
 144     // Simple mdrun tests of energy
 145     { "angles1", { {}, { 1, 2 } } },
 146     // Scaled water for NMA
 147     { "scaled-water", { {}, { 1, 2, 3, 4, 5, 6 } } },
 148     // Villin for NMA
 149     { "villin", { {}, { 1, 2, 3, 4, 5, 6 } } },
 150     // SPC-Dimer for NMA
 151     { "spc-dimer", { {}, { 1, 2, 3, 4, 5, 6 } } },
 152     // SW-Dimer for NMA
 153     { "sw-dimer", { { { "nstcalcenergy", "1" } }, { 1, 2, 3, 4, 5, 6 } } },
 154     // TIP5P for NMA
 155     { "one-tip5p", { {}, { 1, 2, 3, 4, 5, 6 } } },
 156     // ICE-Binding protein for NMA
 157     { "ice-binding", { {}, { 1, 2, 3, 4, 5, 6 } } },
 158     // Nonanol molecule in vacuo, topology suitable for testing FEP
 159     // on KE, angles, dihedral restraints, coulomb and vdw
 160     { "nonanol_vacuo",
 161       { { { "nsteps", "16" },
 162           { "compressibility", "5e-10" },
 163           { "tau-p", "1000" },
 164           { "constraints", "h-bonds" },
 165           { "other",
 166             R"(free-energy         = yes
 167                                   sc-alpha            = 0.5
 168                                   sc-r-power          = 6
 169                                   mass-lambdas        = 0.0 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
 170                                   bonded-lambdas      = 0.0 0.0 0.0 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
 171                                   restraint-lambdas   = 0.0 0.0 0.0 0.0 0.0 0.5 1.0 1.0 1.0 1.0 1.0
 172                                   vdw-lambdas         = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 1.0 1.0 1.0
 173                                   coul-lambdas        = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 1.0
 174                                   ;couple-moltype      = nonanol
 175                                   ;couple-lambda0      = none
 176                                   ;couple-lambda1      = vdw-q
 177                                   ;couple-intramol     = yes)" } },
 178         { 1, 2, 3, 4, 5, 6, 8, 9 } } }
 179 };
 180
 181 /*! \brief Prepare default .mdp values
 182  *
 183  * Insert suitable .mdp defaults, so that \c mdpFileValueDatabase_g
 184  * does not need to specify repetitive values. This works because
 185  * std::map.insert() does not over-write elements that already exist.
 186  *
 187  * \todo ideally some of these default values should be the same as
 188  * grompp uses, and sourced from the same place, but that code is a
 189  * bit of a jungle until we transition to using IMdpOptions more.
 190  *
 191  * \throws  std::bad_alloc     if out of memory
 192  *          std::out_of_range  if \c simulationName is not in the database
 193  *
 194  * Note: Any mdp options that are not added here cannot be used
 195  */
 196 MdpFieldValues prepareDefaultMdpFieldValues(const std::string& simulationName)
 197 {
 198     using MdpField = MdpFieldValues::value_type;
 199
 200     auto mdpFieldValues = mdpFileValueDatabase_g.at(simulationName).mdpFieldValues;
 201     mdpFieldValues.insert(MdpField("nsteps", "16"));
 202     mdpFieldValues.insert(MdpField("nstenergy", "4"));
 203     mdpFieldValues.insert(MdpField("nstxout", "4"));
 204     mdpFieldValues.insert(MdpField("nstvout", "4"));
 205     mdpFieldValues.insert(MdpField("nstfout", "4"));
 206     mdpFieldValues.insert(MdpField("nstxout-compressed", "0"));
 207     mdpFieldValues.insert(MdpField("nstdhdl", "4"));
 208     mdpFieldValues.insert(MdpField("comm-mode", "linear"));
 209     mdpFieldValues.insert(MdpField("nstcomm", "4"));
 210     mdpFieldValues.insert(MdpField("ref-t", "298"));
 211     mdpFieldValues.insert(MdpField("nsttcouple", "4"));
 212     mdpFieldValues.insert(MdpField("tau-p", "1"));
 213     mdpFieldValues.insert(MdpField("nstpcouple", "4"));
 214     mdpFieldValues.insert(MdpField("compressibility", "5e-5"));
 215     mdpFieldValues.insert(MdpField("constraints", "none"));
 216     mdpFieldValues.insert(MdpField("other", ""));
 217     mdpFieldValues.insert(MdpField("coulombtype", "Cut-off"));
 218     mdpFieldValues.insert(MdpField("rcoulomb", "0.7"));
 219     mdpFieldValues.insert(MdpField("vdwtype", "Cut-off"));
 220     mdpFieldValues.insert(MdpField("rvdw", "0.7"));
 221     mdpFieldValues.insert(MdpField("nstcalcenergy", "100"));
 222
 223     return mdpFieldValues;
 224 }
 225
 226 } // namespace
 227
 228 bool isNumberOfPpRanksSupported(const std::string& simulationName, int possibleNumberOfPpRanks)
 229 {
 230     const auto& possibleNumbers = mdpFileValueDatabase_g.at(simulationName).validPpRankCounts;
 231     return (std::find(std::begin(possibleNumbers), std::end(possibleNumbers), possibleNumberOfPpRanks)
 232             != std::end(possibleNumbers));
 233 }
 234
 235 std::string reportNumbersOfPpRanksSupported(const std::string& simulationName)
 236 {
 237     const auto& possibleNumbers = mdpFileValueDatabase_g.at(simulationName).validPpRankCounts;
 238     return formatAndJoin(
 239             std::begin(possibleNumbers), std::end(possibleNumbers), ",", StringFormatter("%d"));
 240 }
 241
 242 MdpFieldValues prepareMdpFieldValues(const std::string& simulationName,
 243                                      const std::string& integrator,
 244                                      const std::string& tcoupl,
 245                                      const std::string& pcoupl)
 246 {
 247     using MdpField = MdpFieldValues::value_type;
 248
 249     auto mdpFieldValues = prepareDefaultMdpFieldValues(simulationName);
 250     mdpFieldValues.insert(MdpField("integrator", integrator));
 251     mdpFieldValues.insert(MdpField("tcoupl", tcoupl));
 252     mdpFieldValues.insert(MdpField("pcoupl", pcoupl));
 253     return mdpFieldValues;
 254 }
 255
 256 MdpFieldValues prepareMdpFieldValues(const char* simulationName,
 257                                      const char* integrator,
 258                                      const char* tcoupl,
 259                                      const char* pcoupl)
 260 {
 261     return prepareMdpFieldValues(std::string(simulationName), integrator, tcoupl, pcoupl);
 262 }
 263 std::string prepareMdpFileContents(const MdpFieldValues& mdpFieldValues)
 264 {
 265     /* Set up an .mdp file that permits a highly reproducible
 266      * simulation. The format string needs to be configured with
 267      * values for various .mdp fields to suit the kind of system
 268      * used and testing needed. It also
 269      * - writes frames from different kinds of steps: starting, ending, intermediate NS, intermediate non-NS
 270      * - has other steps between frame-writing steps
 271      * - has enough buffer that e.g. a rerun will compute the same potential energy even though it does NS every frame
 272      * - has very slow pressure coupling and weak compressibility (since otherwise the box will shrink too fast)
 273      * - can have arbitrary chunks of .mdp content appended to it (but it is up to grompp how it deals with duplicate fields)
 274      * - sets random seeds to known values
 275      * - uses cutoffs that fit inside boxes even after GPU mdrun scales rlist
 276      *
 277      * Note that forces computed in the absence of energy computations
 278      * generally follow a different code path from those computed with
 279      * energies. However a rerun always computes energies, so we don't
 280      * currently have a good way to compare forces at steps where
 281      * energies were not computed with those from rerun on the same
 282      * coordinates.
 283      *
 284      * Note: Any mdp options that are not printed here cannot be used
 285      */
 286     return formatString(
 287             R"(coulombtype             = %s
 288                            rcoulomb                = %s
 289                            vdwtype                 = %s
 290                            rvdw                    = %s
 291                            rlist                   = -1
 292                            bd-fric                 = 1000
 293                            verlet-buffer-tolerance = 0.000001
 294                            nsteps                  = %s
 295                            nstenergy               = %s
 296                            nstxout                 = %s
 297                            nstvout                 = %s
 298                            nstfout                 = %s
 299                            nstxout-compressed      = %s
 300                            nstdhdl                 = %s
 301                            nstlist                 = 8
 302                            integrator              = %s
 303                            ld-seed                 = 234262
 304                            tcoupl                  = %s
 305                            nsttcouple              = %s
 306                            ref-t                   = %s
 307                            tau-t                   = 1
 308                            tc-grps                 = System
 309                            pcoupl                  = %s
 310                            nstpcouple              = %s
 311                            pcoupltype              = isotropic
 312                            ref-p                   = 1
 313                            tau-p                   = %s
 314                            compressibility         = %s
 315                            constraints             = %s
 316                            constraint-algorithm    = lincs
 317                            lincs-order             = 2
 318                            lincs-iter              = 5
 319                            nstcalcenergy           = %s
 320                            comm-mode               = %s
 321                            nstcomm                 = %s
 322                            %s)",
 323             mdpFieldValues.at("coulombtype").c_str(),
 324             mdpFieldValues.at("rcoulomb").c_str(),
 325             mdpFieldValues.at("vdwtype").c_str(),
 326             mdpFieldValues.at("rvdw").c_str(),
 327             mdpFieldValues.at("nsteps").c_str(),
 328             mdpFieldValues.at("nstenergy").c_str(),
 329             mdpFieldValues.at("nstxout").c_str(),
 330             mdpFieldValues.at("nstvout").c_str(),
 331             mdpFieldValues.at("nstfout").c_str(),
 332             mdpFieldValues.at("nstxout-compressed").c_str(),
 333             mdpFieldValues.at("nstdhdl").c_str(),
 334             mdpFieldValues.at("integrator").c_str(),
 335             mdpFieldValues.at("tcoupl").c_str(),
 336             mdpFieldValues.at("nsttcouple").c_str(),
 337             mdpFieldValues.at("ref-t").c_str(),
 338             mdpFieldValues.at("pcoupl").c_str(),
 339             mdpFieldValues.at("nstpcouple").c_str(),
 340             mdpFieldValues.at("tau-p").c_str(),
 341             mdpFieldValues.at("compressibility").c_str(),
 342             mdpFieldValues.at("constraints").c_str(),
 343             mdpFieldValues.at("nstcalcenergy").c_str(),
 344             mdpFieldValues.at("comm-mode").c_str(),
 345             mdpFieldValues.at("nstcomm").c_str(),
 346             mdpFieldValues.at("other").c_str());
 347 }
 348
 349 } // namespace test
 350 } // namespace gmx