src/gromacs/awh/tests/bias.cpp

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