src/testutils/simulationdatabase.cpp

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