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