src/gromacs/mdlib/freeenergyparameters.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020, 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  *
  37  * \brief Implements routines in freeenergyparameters.h .
  38  *
  39  * \author Christian Blau <blau@kth.se>
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "freeenergyparameters.h"
  45
  46 #include "gromacs/mdtypes/inputrec.h"
  47
  48 namespace gmx
  49 {
  50
  51 namespace
  52 {
  53
  54 std::array<real, efptNR> lambdasAtState(const int stateIndex, double** const lambdaArray, const int lambdaArrayExtent)
  55 {
  56     std::array<real, efptNR> lambda;
  57     // set lambda from an fep state index from stateIndex, if stateIndex was defined (> -1)
  58     if (stateIndex >= 0 && stateIndex < lambdaArrayExtent)
  59     {
  60         for (int i = 0; i < efptNR; i++)
  61         {
  62             lambda[i] = lambdaArray[i][stateIndex];
  63         }
  64     }
  65     return lambda;
  66 }
  67
  68 /*! \brief Evaluates where in the lambda arrays we are at currently.
  69  *
  70  * \param[in] step the current simulation step
  71  * \param[in] deltaLambdaPerStep the change of the overall controlling lambda parameter per step
  72  * \param[in] initialFEPStateIndex the FEP state at the start of the simulation. -1 if not set.
  73  * \param[in] initialLambda the lambda at the start of the simulation . -1 if not set.
  74  * \param[in] lambdaArrayExtent how many lambda values we have
  75  * \returns a number that reflects where in the lambda arrays we are at the moment
  76  */
  77 double currentGlobalLambda(const int64_t step,
  78                            const double  deltaLambdaPerStep,
  79                            const int     initialFEPStateIndex,
  80                            const double  initialLambda,
  81                            const int     lambdaArrayExtent)
  82 {
  83     const real fracSimulationLambda = step * deltaLambdaPerStep;
  84
  85     // Set initial lambda value for the simulation either from initialFEPStateIndex or,
  86     // if not set, from the initial lambda.
  87     double initialGlobalLambda = 0;
  88     if (initialFEPStateIndex > -1)
  89     {
  90         if (lambdaArrayExtent > 1)
  91         {
  92             initialGlobalLambda = static_cast<double>(initialFEPStateIndex) / (lambdaArrayExtent - 1);
  93         }
  94     }
  95     else
  96     {
  97         if (initialLambda > -1)
  98         {
  99             initialGlobalLambda = initialLambda;
 100         }
 101     }
 102
 103     return initialGlobalLambda + fracSimulationLambda;
 104 }
 105
 106 /*! \brief Returns an array of lambda values from linear interpolation of a lambda value matrix.
 107  *
 108  * \note If there is nothing to interpolate, fills the array with the global current lambda.
 109  * \note Returns array boundary values if currentGlobalLambda <0 or >1 .
 110  *
 111  * \param[in] currentGlobalLambda determines at which position in the lambda array to interpolate
 112  * \param[in] lambdaArray array of lambda values
 113  * \param[in] lambdaArrayExtent number of lambda values
 114  */
 115 std::array<real, efptNR> interpolatedLambdas(const double   currentGlobalLambda,
 116                                              double** const lambdaArray,
 117                                              const int      lambdaArrayExtent)
 118 {
 119     std::array<real, efptNR> lambda;
 120     // when there is no lambda value array, set all lambdas to steps * deltaLambdaPerStep
 121     if (lambdaArrayExtent <= 0)
 122     {
 123         std::fill(std::begin(lambda), std::end(lambda), currentGlobalLambda);
 124         return lambda;
 125     }
 126
 127     // if we run over the boundary of the lambda array, return the boundary array values
 128     if (currentGlobalLambda <= 0)
 129     {
 130         for (int i = 0; i < efptNR; i++)
 131         {
 132             lambda[i] = lambdaArray[i][0];
 133         }
 134         return lambda;
 135     }
 136     if (currentGlobalLambda >= 1)
 137     {
 138         for (int i = 0; i < efptNR; i++)
 139         {
 140             lambda[i] = lambdaArray[i][lambdaArrayExtent - 1];
 141         }
 142         return lambda;
 143     }
 144
 145     // find out between which two value lambda array elements to interpolate
 146     const int fepStateLeft = static_cast<int>(std::floor(currentGlobalLambda * (lambdaArrayExtent - 1)));
 147     const int fepStateRight = fepStateLeft + 1;
 148     // interpolate between this state and the next
 149     const double fracBetween = currentGlobalLambda * (lambdaArrayExtent - 1) - fepStateLeft;
 150     for (int i = 0; i < efptNR; i++)
 151     {
 152         lambda[i] = lambdaArray[i][fepStateLeft]
 153                     + fracBetween * (lambdaArray[i][fepStateRight] - lambdaArray[i][fepStateLeft]);
 154     }
 155     return lambda;
 156 }
 157
 158 } // namespace
 159
 160 std::array<real, efptNR> currentLambdas(const int64_t step, const t_lambda& fepvals, const int currentLambdaState)
 161 {
 162     if (fepvals.delta_lambda == 0)
 163     {
 164         if (currentLambdaState > -1)
 165         {
 166             return lambdasAtState(currentLambdaState, fepvals.all_lambda, fepvals.n_lambda);
 167         }
 168
 169         if (fepvals.init_fep_state > -1)
 170         {
 171             return lambdasAtState(fepvals.init_fep_state, fepvals.all_lambda, fepvals.n_lambda);
 172         }
 173
 174         std::array<real, efptNR> lambdas;
 175         std::fill(std::begin(lambdas), std::end(lambdas), fepvals.init_lambda);
 176         return lambdas;
 177     }
 178     const double globalLambda = currentGlobalLambda(step, fepvals.delta_lambda, fepvals.init_fep_state,
 179                                                     fepvals.init_lambda, fepvals.n_lambda);
 180     return interpolatedLambdas(globalLambda, fepvals.all_lambda, fepvals.n_lambda);
 181 }
 182
 183 } // namespace gmx