src/gromacs/awh/tests/bias.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2017,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 #include "gmxpre.h"
  36
  37 #include "gromacs/awh/bias.h"
  38
  39 #include <cmath>
  40
  41 #include <memory>
  42 #include <tuple>
  43 #include <vector>
  44
  45 #include <gmock/gmock.h>
  46 #include <gtest/gtest.h>
  47
  48 #include "gromacs/awh/correlationgrid.h"
  49 #include "gromacs/awh/pointstate.h"
  50 #include "gromacs/mdtypes/awh-params.h"
  51 #include "gromacs/utility/stringutil.h"
  52
  53 #include "testutils/refdata.h"
  54 #include "testutils/testasserts.h"
  55
  56 namespace gmx
  57 {
  58
  59 namespace test
  60 {
  61
  62 /*! \internal \brief
  63  * Struct that gathers all input for setting up and using a Bias
  64  */
  65 struct AwhTestParameters
  66 {
  67     AwhTestParameters() = default;
  68     //!Move constructor
  69     AwhTestParameters(AwhTestParameters &&o) noexcept : beta(o.beta), awhDimParams(o.awhDimParams),
  70                                                         awhBiasParams(o.awhBiasParams), awhParams(o.awhParams),
  71                                                         dimParams(std::move(o.dimParams))
  72     {
  73         awhBiasParams.dimParams = &awhDimParams;
  74         awhParams.awhBiasParams = &awhBiasParams;
  75     }
  76     double                 beta;          //!< 1/(kB*T)
  77
  78     AwhDimParams           awhDimParams;  //!< Dimension parameters pointed to by \p awhBiasParams
  79     AwhBiasParams          awhBiasParams; //!< Bias parameters pointed to by \[ awhParams
  80     AwhParams              awhParams;     //!< AWH parameters, this is the struct to actually use
  81
  82     std::vector<DimParams> dimParams;     //!< Dimension parameters for setting up Bias
  83 };
  84
  85 //! Helper function to set up the C-style AWH parameters for the test
  86 static AwhTestParameters getAwhTestParameters(int eawhgrowth,
  87                                               int eawhpotential)
  88 {
  89     AwhTestParameters params;
  90
  91     params.beta = 0.4;
  92
  93     AwhDimParams &awhDimParams = params.awhDimParams;
  94
  95     awhDimParams.period           = 0;
  96     awhDimParams.diffusion        = 0.1;
  97     awhDimParams.origin           = 0.5;
  98     awhDimParams.end              = 1.5;
  99     awhDimParams.coordValueInit   = awhDimParams.origin;
 100     awhDimParams.coverDiameter    = 0;
 101
 102     AwhBiasParams &awhBiasParams = params.awhBiasParams;
 103
 104     awhBiasParams.ndim                 = 1;
 105     awhBiasParams.dimParams            = &awhDimParams;
 106     awhBiasParams.eTarget              = eawhtargetCONSTANT;
 107     awhBiasParams.targetBetaScaling    = 0;
 108     awhBiasParams.targetCutoff         = 0;
 109     awhBiasParams.eGrowth              = eawhgrowth;
 110     awhBiasParams.bUserData            = FALSE;
 111     awhBiasParams.errorInitial         = 0.5/params.beta;
 112     awhBiasParams.shareGroup           = 0;
 113     awhBiasParams.equilibrateHistogram = FALSE;
 114
 115     double                 convFactor  = 1;
 116     double                 k           = 1000;
 117     int64_t                seed        = 93471803;
 118
 119     params.dimParams.emplace_back(convFactor, k, params.beta);
 120
 121     AwhParams             &awhParams = params.awhParams;
 122
 123     awhParams.numBias                    = 1;
 124     awhParams.awhBiasParams              = &awhBiasParams;
 125     awhParams.seed                       = seed;
 126     awhParams.nstOut                     = 0;
 127     awhParams.nstSampleCoord             = 1;
 128     awhParams.numSamplesUpdateFreeEnergy = 10;
 129     awhParams.ePotential                 = eawhpotential;
 130     awhParams.shareBiasMultisim          = FALSE;
 131
 132     return params;
 133 }
 134
 135 //! Database of 21 test coordinates that represent a trajectory */
 136 const double g_coords[] = {
 137     0.62,
 138     0.70,
 139     0.68,
 140     0.80,
 141     0.93,
 142     0.87,
 143     1.16,
 144     1.14,
 145     0.95,
 146     0.89,
 147     0.91,
 148     0.86,
 149     0.88,
 150     0.79,
 151     0.75,
 152     0.82,
 153     0.74,
 154     0.70,
 155     0.68,
 156     0.71,
 157     0.73
 158 };
 159
 160 //! Convenience typedef: growth type enum, potential type enum, disable update skips
 161 typedef std::tuple<int, int, BiasParams::DisableUpdateSkips> BiasTestParameters;
 162
 163 /*! \brief Test fixture for testing Bias updates
 164  */
 165 class BiasTest : public ::testing::TestWithParam<BiasTestParameters>
 166 {
 167     public:
 168         //! Random seed for AWH MC sampling
 169         int64_t           seed_;
 170
 171         //! Coordinates representing a trajectory in time
 172         std::vector<double>   coordinates_;
 173         //! The awh Bias
 174         std::unique_ptr<Bias> bias_;
 175
 176         BiasTest() :
 177             coordinates_(std::begin(g_coords), std::end(g_coords))
 178         {
 179             /* We test all combinations of:
 180              *   eawhgrowth:
 181              *     eawhgrowthLINEAR:     final, normal update phase
 182              *     ewahgrowthEXP_LINEAR: intial phase, updated size is constant
 183              *   eawhpotential (should only affect the force output):
 184              *     eawhpotentialUMBRELLA:  MC on lambda (umbrella potential location)
 185              *     eawhpotentialCONVOLVED: MD on a convolved potential landscape
 186              *   disableUpdateSkips (should not affect the results):
 187              *     BiasParams::DisableUpdateSkips::yes: update the point state for every sample
 188              *     BiasParams::DisableUpdateSkips::no:  update the point state at an interval > 1 sample
 189              *
 190              * Note: It would be nice to explicitly check that eawhpotential
 191              *       and disableUpdateSkips do not affect the point state.
 192              *       But the reference data will also ensure this.
 193              */
 194             int eawhgrowth;
 195             int eawhpotential;
 196             BiasParams::DisableUpdateSkips disableUpdateSkips;
 197             std::tie(eawhgrowth, eawhpotential, disableUpdateSkips) = GetParam();
 198
 199             /* Set up a basic AWH setup with a single, 1D bias with parameters
 200              * such that we can measure the effects of different parameters.
 201              * The idea is to, among other things, have part of the interval
 202              * not covered by samples.
 203              */
 204             const AwhTestParameters params = getAwhTestParameters(eawhgrowth, eawhpotential);
 205
 206             seed_ = params.awhParams.seed;
 207
 208             double mdTimeStep = 0.1;
 209
 210             int    numSamples = coordinates_.size() - 1; // No sample taken at step 0
 211             GMX_RELEASE_ASSERT(numSamples % params.awhParams.numSamplesUpdateFreeEnergy == 0, "This test is intended to reproduce the situation when the might need to write output during a normal AWH run, therefore the number of samples should be a multiple of the free-energy update interval (but the test should also runs fine without this condition).");
 212
 213             bias_ = std::make_unique<Bias>(-1, params.awhParams, params.awhBiasParams, params.dimParams, params.beta, mdTimeStep, 1, "", Bias::ThisRankWillDoIO::No, disableUpdateSkips);
 214         }
 215 };
 216
 217 TEST_P(BiasTest, ForcesBiasPmf)
 218 {
 219     gmx::test::TestReferenceData     data;
 220     gmx::test::TestReferenceChecker  checker(data.rootChecker());
 221
 222     Bias                            &bias = *bias_;
 223
 224     /* Make strings with the properties we expect to be different in the tests.
 225      * These also helps to interpret the reference data.
 226      */
 227     std::vector<std::string>         props;
 228     props.push_back(formatString("stage:           %s", bias.state().inInitialStage() ? "initial" : "final"));
 229     props.push_back(formatString("convolve forces: %s", bias.params().convolveForce ? "yes" : "no"));
 230     props.push_back(formatString("skip updates:    %s", bias.params().skipUpdates() ? "yes" : "no"));
 231
 232     SCOPED_TRACE(gmx::formatString("%s, %s, %s", props[0].c_str(), props[1].c_str(), props[2].c_str()));
 233
 234     std::vector<double> force, pot, potJump;
 235
 236     double              coordMaxValue = 0;
 237     double              potentialJump = 0;
 238     int64_t             step          = 0;
 239     for (auto &coord : coordinates_)
 240     {
 241         coordMaxValue = std::max(coordMaxValue, std::abs(coord));
 242
 243         awh_dvec                    coordValue = { coord, 0, 0, 0 };
 244         double                      potential  = 0;
 245         gmx::ArrayRef<const double> biasForce  =
 246             bias.calcForceAndUpdateBias(coordValue,
 247                                         &potential, &potentialJump,
 248                                         nullptr, nullptr, step, step, seed_,
 249                                         nullptr);
 250
 251         force.push_back(biasForce[0]);
 252         pot.push_back(potential);
 253         potJump.push_back(potentialJump);
 254
 255         step++;
 256     }
 257
 258     /* When skipping updates, ensure all skipped updates are performed here.
 259      * This should result in the same bias state as at output in a normal run.
 260      */
 261     if (bias.params().skipUpdates())
 262     {
 263         bias.doSkippedUpdatesForAllPoints();
 264     }
 265
 266     std::vector<double> pointBias, logPmfsum;
 267     for (auto &point : bias.state().points())
 268     {
 269         pointBias.push_back(point.bias());
 270         logPmfsum.push_back(point.logPmfSum());
 271     }
 272
 273     /* The umbrella force is computed from the coordinate deviation.
 274      * In taking this deviation we lose a lot of precision, so we should
 275      * compare against k*max(coord) instead of the instantaneous force.
 276      */
 277     const double  kCoordMax = bias.dimParams()[0].k*coordMaxValue;
 278
 279     constexpr int ulpTol    = 10;
 280
 281     checker.checkSequence(props.begin(),     props.end(),     "Properties");
 282     checker.setDefaultTolerance(absoluteTolerance(kCoordMax*GMX_DOUBLE_EPS*ulpTol));
 283     checker.checkSequence(force.begin(),     force.end(),     "Force");
 284     checker.checkSequence(pot.begin(),       pot.end(),       "Potential");
 285     checker.checkSequence(potJump.begin(),   potJump.end(),   "PotentialJump");
 286     checker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpTol));
 287     checker.checkSequence(pointBias.begin(), pointBias.end(), "PointBias");
 288     checker.checkSequence(logPmfsum.begin(), logPmfsum.end(), "PointLogPmfsum");
 289 }
 290
 291 /* Scan initial/final phase, MC/convolved force and update skip (not) allowed
 292  * Both the convolving and skipping should not affect the bias and PMF.
 293  * It would be nice if the test would explicitly check for this.
 294  * Currently this is tested through identical reference data.
 295  */
 296 INSTANTIATE_TEST_CASE_P(WithParameters, BiasTest,
 297                             ::testing::Combine(
 298                                     ::testing::Values(eawhgrowthLINEAR, eawhgrowthEXP_LINEAR),
 299                                     ::testing::Values(eawhpotentialUMBRELLA, eawhpotentialCONVOLVED),
 300                                     ::testing::Values(BiasParams::DisableUpdateSkips::yes, BiasParams::DisableUpdateSkips::no)));
 301
 302 // Test that we detect coverings and exit the initial stage at the correct step
 303 TEST(BiasTest, DetectsCovering)
 304 {
 305     const AwhTestParameters params       = getAwhTestParameters(eawhgrowthEXP_LINEAR, eawhpotentialCONVOLVED);
 306     const AwhDimParams     &awhDimParams = params.awhParams.awhBiasParams[0].dimParams[0];
 307
 308     const double            mdTimeStep   = 0.1;
 309
 310     Bias                    bias(-1, params.awhParams, params.awhBiasParams, params.dimParams, params.beta, mdTimeStep, 1, "", Bias::ThisRankWillDoIO::No);
 311
 312     /* We use a trajectory of the sum of two sines to cover the reaction
 313      * coordinate range in a semi-realistic way. The period is 4*pi=12.57.
 314      * We get out of the initial stage after 4 coverings at step 300.
 315      */
 316     const int64_t     exitStepRef = 300;
 317     const double      midPoint    = 0.5*(awhDimParams.end + awhDimParams.origin);
 318     const double      halfWidth   = 0.5*(awhDimParams.end - awhDimParams.origin);
 319
 320     bool              inInitialStage = bias.state().inInitialStage();
 321     /* Normally this loop exits at exitStepRef, but we extend with failure */
 322     int64_t           step;
 323     for (step = 0; step <= 2*exitStepRef; step++)
 324     {
 325         double   t     = step*mdTimeStep;
 326         double   coord = midPoint + halfWidth*(0.5*std::sin(t) + 0.55*std::sin(1.5*t));
 327
 328         awh_dvec coordValue    = { coord, 0, 0, 0 };
 329         double   potential     = 0;
 330         double   potentialJump = 0;
 331         bias.calcForceAndUpdateBias(coordValue,
 332                                     &potential, &potentialJump,
 333                                     nullptr, nullptr,
 334                                     step, step, params.awhParams.seed, nullptr);
 335
 336         inInitialStage = bias.state().inInitialStage();
 337         if (!inInitialStage)
 338         {
 339             break;
 340         }
 341     }
 342
 343     EXPECT_EQ(false, inInitialStage);
 344     if (!inInitialStage)
 345     {
 346         EXPECT_EQ(exitStepRef, step);
 347     }
 348 }
 349
 350 }  // namespace test
 351 }  // namespace gmx