src/gromacs/mdlib/freeenergyparameters.cpp

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