src/testutils/simulationdatabase.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2016,2018,2019,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 /*! \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     // Capped alanine peptide in vacuo (no virtual sites)
 141     { "alanine_vacuo", { { { "constraints", "h-bonds" } }, { 1, 2, 3, 4, 5, 6, 7, 8, 9 } } },
 142     // Zwitterionic glycine in vacuo
 143     { "glycine_vacuo", { { { "constraints", "h-bonds" } }, { 1, 2, 3, 4, 5, 6, 7, 8, 9 } } },
 144     // Zwitterionic glycine in vacuo, without constraints
 145     { "glycine_no_constraints_vacuo", { { { "constraints", "none" } }, { 1, 2, 3, 4, 5, 6, 7, 8, 9 } } },
 146     // Simple mdrun tests of energy
 147     { "angles1", { {}, { 1, 2 } } },
 148     // Scaled water for NMA
 149     { "scaled-water", { {}, { 1, 2, 3, 4, 5, 6 } } },
 150     // Villin for NMA
 151     { "villin", { {}, { 1, 2, 3, 4, 5, 6 } } },
 152     // SPC-Dimer for NMA
 153     { "spc-dimer", { {}, { 1, 2, 3, 4, 5, 6 } } },
 154     // SW-Dimer for NMA
 155     { "sw-dimer", { { { "nstcalcenergy", "1" } }, { 1, 2, 3, 4, 5, 6 } } },
 156     // TIP5P for NMA
 157     { "one-tip5p", { {}, { 1, 2, 3, 4, 5, 6 } } },
 158     // ICE-Binding protein for NMA
 159     { "ice-binding", { {}, { 1, 2, 3, 4, 5, 6 } } },
 160     // Nonanol molecule in vacuo, topology suitable for testing FEP
 161     // on KE, angles, dihedral restraints, coulomb and vdw
 162     { "nonanol_vacuo",
 163       { { { "nsteps", "16" },
 164           { "compressibility", "5e-10" },
 165           { "tau-p", "1000" },
 166           { "constraints", "h-bonds" },
 167           { "free-energy", "yes" },
 168           { "sc-alpha", "0.5" },
 169           { "sc-r-power", "6" },
 170           { "mass-lambdas", "0.0 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0" },
 171           { "bonded-lambdas", "0.0 0.0 0.0 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0" },
 172           { "restraint-lambdas", "0.0 0.0 0.0 0.0 0.0 0.5 1.0 1.0 1.0 1.0 1.0" },
 173           { "vdw-lambdas", "0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 1.0 1.0 1.0" },
 174           { "coul-lambdas", "0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 1.0" },
 175           { ";couple-moltype", "nonanol" },
 176           { ";couple-lambda0", "none" },
 177           { ";couple-lambda1", "vdw-q" },
 178           { ";couple-intramol", "yes" } },
 179         { 1, 2, 3, 4, 5, 6, 8, 9 } } },
 180     // Artificial test system including all virtual site types
 181     { "vsite_test", { {}, { 1, 2, 3, 4, 5, 6, 7, 8, 9 } } },
 182 };
 183
 184 //! Helper typedef for mdp database
 185 using MdpDatabase = std::map<MdpParameterDatabase, MdpFieldValues>;
 186 //! Database of additional mdp options used for specific algorithms
 187 const MdpDatabase c_additionalMdpOptions{ { MdpParameterDatabase::Default, {} },
 188                                           { MdpParameterDatabase::Pull,
 189                                             { { "coulombtype", "reaction-field" },
 190                                               { "pull", "yes" },
 191                                               // Prev step reference is checkpointed - rest of pull is not cpt-dependent
 192                                               { "pull-pbc-ref-prev-step-com", "yes" },
 193                                               { "pull-ngroups", "2" },
 194                                               { "pull-group1-name", "FirstWaterMolecule" },
 195                                               { "pull-group2-name", "SecondWaterMolecule" },
 196                                               { "pull-ncoords", "1" },
 197                                               { "pull-coord1-type", "umbrella" },
 198                                               { "pull-coord1-geometry", "distance" },
 199                                               { "pull-coord1-groups", "1 2" },
 200                                               { "pull-coord1-init", "1" },
 201                                               { "pull-coord1-k", "10000" } } },
 202                                           { MdpParameterDatabase::Awh,
 203                                             { { "pull", "yes" },
 204                                               { "pull-ngroups", "5" },
 205                                               { "pull-ncoords", "2" },
 206                                               { "pull-group1-name", "C_&_r_1" },
 207                                               { "pull-group2-name", "N_&_r_2" },
 208                                               { "pull-group3-name", "CA" },
 209                                               { "pull-group4-name", "C_&_r_2" },
 210                                               { "pull-group5-name", "N_&_r_3" },
 211                                               { "pull-coord1-geometry", "dihedral" },
 212                                               { "pull-coord1-groups", "1 2 2 3 3 4" },
 213                                               { "pull-coord1-k", "4000" },
 214                                               { "pull-coord1-kB", "1000" },
 215                                               { "pull-coord2-geometry", "dihedral" },
 216                                               { "pull-coord2-groups", "2 3 3 4 4 5" },
 217                                               { "pull-coord2-k", "4000" },
 218                                               { "pull-coord2-kB", "1000" },
 219                                               { "pull-coord1-type", "external-potential" },
 220                                               { "pull-coord1-potential-provider", "awh" },
 221                                               { "pull-coord2-type", "external-potential" },
 222                                               { "pull-coord2-potential-provider", "awh" },
 223                                               { "awh", "yes" },
 224                                               { "awh-potential", "convolved" },
 225                                               { "awh-nstout", "4" },
 226                                               { "awh-nstsample", "4" },
 227                                               { "awh-nsamples-update", "1" },
 228                                               { "awh-share-multisim", "no" },
 229                                               { "awh-nbias", "2" },
 230                                               { "awh1-ndim", "1" },
 231                                               { "awh1-dim1-coord-index", "2" },
 232                                               { "awh1-dim1-start", "150" },
 233                                               { "awh1-dim1-end", "180" },
 234                                               { "awh1-dim1-force-constant", "4000" },
 235                                               { "awh1-dim1-diffusion", "0.1" },
 236                                               { "awh2-ndim", "1" },
 237                                               { "awh2-dim1-coord-index", "1" },
 238                                               { "awh2-dim1-start", "178" },
 239                                               { "awh2-dim1-end", "-178" },
 240                                               { "awh2-dim1-force-constant", "4000" },
 241                                               { "awh2-dim1-diffusion", "0.1" } } } };
 242
 243 /*! \brief Prepare default .mdp values
 244  *
 245  * Insert suitable .mdp defaults, so that \c mdpFileValueDatabase_g
 246  * does not need to specify repetitive values. This works because
 247  * std::map.insert() does not over-write elements that already exist.
 248  *
 249  * \todo ideally some of these default values should be the same as
 250  * grompp uses, and sourced from the same place, but that code is a
 251  * bit of a jungle until we transition to using IMdpOptions more.
 252  *
 253  * \throws  std::bad_alloc     if out of memory
 254  *          std::out_of_range  if \c simulationName is not in the database
 255  */
 256 MdpFieldValues prepareDefaultMdpFieldValues(const std::string& simulationName)
 257 {
 258     using MdpField = MdpFieldValues::value_type;
 259
 260     auto mdpFieldValues = mdpFileValueDatabase_g.at(simulationName).mdpFieldValues;
 261     mdpFieldValues.insert(MdpField("nsteps", "16"));
 262     mdpFieldValues.insert(MdpField("nstenergy", "4"));
 263     mdpFieldValues.insert(MdpField("nstxout", "4"));
 264     mdpFieldValues.insert(MdpField("nstvout", "4"));
 265     mdpFieldValues.insert(MdpField("nstfout", "4"));
 266     mdpFieldValues.insert(MdpField("nstxout-compressed", "0"));
 267     mdpFieldValues.insert(MdpField("nstdhdl", "4"));
 268     mdpFieldValues.insert(MdpField("comm-mode", "linear"));
 269     mdpFieldValues.insert(MdpField("nstcomm", "4"));
 270     mdpFieldValues.insert(MdpField("ref-t", "298"));
 271     mdpFieldValues.insert(MdpField("nsttcouple", "4"));
 272     mdpFieldValues.insert(MdpField("tau-p", "1"));
 273     mdpFieldValues.insert(MdpField("nstpcouple", "4"));
 274     mdpFieldValues.insert(MdpField("compressibility", "5e-5"));
 275     mdpFieldValues.insert(MdpField("constraints", "none"));
 276     mdpFieldValues.insert(MdpField("coulombtype", "Cut-off"));
 277     mdpFieldValues.insert(MdpField("rcoulomb", "0.7"));
 278     mdpFieldValues.insert(MdpField("vdwtype", "Cut-off"));
 279     mdpFieldValues.insert(MdpField("rvdw", "0.7"));
 280     mdpFieldValues.insert(MdpField("nstcalcenergy", "100"));
 281     mdpFieldValues.insert(MdpField("rlist", "-1"));
 282     mdpFieldValues.insert(MdpField("bd-fric", "1000"));
 283     mdpFieldValues.insert(MdpField("verlet-buffer-tolerance", "0.000001"));
 284     mdpFieldValues.insert(MdpField("nstlist", "8"));
 285     mdpFieldValues.insert(MdpField("ld-seed", "234262"));
 286     mdpFieldValues.insert(MdpField("tau-t", "1"));
 287     mdpFieldValues.insert(MdpField("tc-grps", "System"));
 288     mdpFieldValues.insert(MdpField("pcoupltype", "isotropic"));
 289     mdpFieldValues.insert(MdpField("ref-p", "1"));
 290     mdpFieldValues.insert(MdpField("constraint-algorithm", "lincs"));
 291     mdpFieldValues.insert(MdpField("lincs-order", "2"));
 292     mdpFieldValues.insert(MdpField("lincs-iter", "5"));
 293
 294     return mdpFieldValues;
 295 }
 296
 297 } // namespace
 298
 299 bool isNumberOfPpRanksSupported(const std::string& simulationName, int possibleNumberOfPpRanks)
 300 {
 301     const auto& possibleNumbers = mdpFileValueDatabase_g.at(simulationName).validPpRankCounts;
 302     return (std::find(std::begin(possibleNumbers), std::end(possibleNumbers), possibleNumberOfPpRanks)
 303             != std::end(possibleNumbers));
 304 }
 305
 306 std::string reportNumbersOfPpRanksSupported(const std::string& simulationName)
 307 {
 308     const auto& possibleNumbers = mdpFileValueDatabase_g.at(simulationName).validPpRankCounts;
 309     return formatAndJoin(
 310             std::begin(possibleNumbers), std::end(possibleNumbers), ",", StringFormatter("%d"));
 311 }
 312
 313 MdpFieldValues prepareMdpFieldValues(const std::string&   simulationName,
 314                                      const std::string&   integrator,
 315                                      const std::string&   tcoupl,
 316                                      const std::string&   pcoupl,
 317                                      MdpParameterDatabase additionalMdpParameters)
 318 {
 319     using MdpField = MdpFieldValues::value_type;
 320
 321     auto mdpFieldValues = prepareDefaultMdpFieldValues(simulationName);
 322     mdpFieldValues.insert(MdpField("integrator", integrator));
 323     mdpFieldValues.insert(MdpField("tcoupl", tcoupl));
 324     mdpFieldValues.insert(MdpField("pcoupl", pcoupl));
 325     for (const auto& mdpField : c_additionalMdpOptions.at(additionalMdpParameters))
 326     {
 327         // Here, we are overwriting default values - we assume the additional
 328         // parameters take precedence over the default parameters
 329         mdpFieldValues[mdpField.first] = mdpField.second;
 330     }
 331
 332     return mdpFieldValues;
 333 }
 334
 335 MdpFieldValues prepareMdpFieldValues(const char*          simulationName,
 336                                      const char*          integrator,
 337                                      const char*          tcoupl,
 338                                      const char*          pcoupl,
 339                                      MdpParameterDatabase additionalMdpParameters)
 340 {
 341     return prepareMdpFieldValues(
 342             std::string(simulationName), integrator, tcoupl, pcoupl, additionalMdpParameters);
 343 }
 344 std::string prepareMdpFileContents(const MdpFieldValues& mdpFieldValues)
 345 {
 346     /* Set up an .mdp file that permits a highly reproducible
 347      * simulation. The format string needs to be configured with
 348      * values for various .mdp fields to suit the kind of system
 349      * used and testing needed. It also
 350      * - writes frames from different kinds of steps: starting, ending, intermediate NS, intermediate non-NS
 351      * - has other steps between frame-writing steps
 352      * - has enough buffer that e.g. a rerun will compute the same potential energy even though it does NS every frame
 353      * - has very slow pressure coupling and weak compressibility (since otherwise the box will shrink too fast)
 354      * - can have arbitrary chunks of .mdp content appended to it (but it is up to grompp how it deals with duplicate fields)
 355      * - sets random seeds to known values
 356      * - uses cutoffs that fit inside boxes even after GPU mdrun scales rlist
 357      *
 358      * Note that forces computed in the absence of energy computations
 359      * generally follow a different code path from those computed with
 360      * energies. However a rerun always computes energies, so we don't
 361      * currently have a good way to compare forces at steps where
 362      * energies were not computed with those from rerun on the same
 363      * coordinates.
 364      */
 365     std::string optionString;
 366     for (auto const& [key, value] : mdpFieldValues)
 367     {
 368         optionString += formatString("%-24s = %s\n", key.c_str(), value.c_str());
 369     }
 370     return optionString;
 371 }
 372
 373 } // namespace test
 374 } // namespace gmx