old_start_lambda = fr->block[i].sub[0].dval[3];
delta_lambda = fr->block[i].sub[0].dval[4];
- if (delta_lambda > 0)
+ if (delta_lambda != 0)
{
gmx_fatal(FARGS, "Lambda values not constant in %s: can't apply BAR method", fn);
}
/* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
if (ir->efep != efepNO)
{
- if (fep->delta_lambda > 0)
+ if (fep->delta_lambda != 0)
{
ir->efep = efepSLOWGROWTH;
}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief Implements routines in freeenergyparameters.h .
+ *
+ * \author Christian Blau <blau@kth.se>
+ */
+
+#include "gmxpre.h"
+
+#include "freeenergyparameters.h"
+
+#include "gromacs/mdtypes/inputrec.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+std::array<real, efptNR> lambdasAtState(const int stateIndex, double** const lambdaArray, const int lambdaArrayExtent)
+{
+ std::array<real, efptNR> lambda;
+ // set lambda from an fep state index from stateIndex, if stateIndex was defined (> -1)
+ if (stateIndex >= 0 && stateIndex < lambdaArrayExtent)
+ {
+ for (int i = 0; i < efptNR; i++)
+ {
+ lambda[i] = lambdaArray[i][stateIndex];
+ }
+ }
+ return lambda;
+}
+
+/*! \brief Evaluates where in the lambda arrays we are at currently.
+ *
+ * \param[in] step the current simulation step
+ * \param[in] deltaLambdaPerStep the change of the overall controlling lambda parameter per step
+ * \param[in] initialFEPStateIndex the FEP state at the start of the simulation. -1 if not set.
+ * \param[in] initialLambda the lambda at the start of the simulation . -1 if not set.
+ * \param[in] lambdaArrayExtent how many lambda values we have
+ * \returns a number that reflects where in the lambda arrays we are at the moment
+ */
+double currentGlobalLambda(const int64_t step,
+ const double deltaLambdaPerStep,
+ const int initialFEPStateIndex,
+ const double initialLambda,
+ const int lambdaArrayExtent)
+{
+ const real fracSimulationLambda = step * deltaLambdaPerStep;
+
+ // Set initial lambda value for the simulation either from initialFEPStateIndex or,
+ // if not set, from the initial lambda.
+ double initialGlobalLambda = 0;
+ if (initialFEPStateIndex > -1)
+ {
+ if (lambdaArrayExtent > 1)
+ {
+ initialGlobalLambda = static_cast<double>(initialFEPStateIndex) / (lambdaArrayExtent - 1);
+ }
+ }
+ else
+ {
+ if (initialLambda > -1)
+ {
+ initialGlobalLambda = initialLambda;
+ }
+ }
+
+ return initialGlobalLambda + fracSimulationLambda;
+}
+
+/*! \brief Returns an array of lambda values from linear interpolation of a lambda value matrix.
+ *
+ * \note If there is nothing to interpolate, fills the array with the global current lambda.
+ * \note Returns array boundary values if currentGlobalLambda <0 or >1 .
+ *
+ * \param[in] currentGlobalLambda determines at which position in the lambda array to interpolate
+ * \param[in] lambdaArray array of lambda values
+ * \param[in] lambdaArrayExtent number of lambda values
+ */
+std::array<real, efptNR> interpolatedLambdas(const double currentGlobalLambda,
+ double** const lambdaArray,
+ const int lambdaArrayExtent)
+{
+ std::array<real, efptNR> lambda;
+ // when there is no lambda value array, set all lambdas to steps * deltaLambdaPerStep
+ if (lambdaArrayExtent <= 0)
+ {
+ std::fill(std::begin(lambda), std::end(lambda), currentGlobalLambda);
+ return lambda;
+ }
+
+ // if we run over the boundary of the lambda array, return the boundary array values
+ if (currentGlobalLambda <= 0)
+ {
+ for (int i = 0; i < efptNR; i++)
+ {
+ lambda[i] = lambdaArray[i][0];
+ }
+ return lambda;
+ }
+ if (currentGlobalLambda >= 1)
+ {
+ for (int i = 0; i < efptNR; i++)
+ {
+ lambda[i] = lambdaArray[i][lambdaArrayExtent - 1];
+ }
+ return lambda;
+ }
+
+ // find out between which two value lambda array elements to interpolate
+ const int fepStateLeft = static_cast<int>(std::floor(currentGlobalLambda * (lambdaArrayExtent - 1)));
+ const int fepStateRight = fepStateLeft + 1;
+ // interpolate between this state and the next
+ const double fracBetween = currentGlobalLambda * (lambdaArrayExtent - 1) - fepStateLeft;
+ for (int i = 0; i < efptNR; i++)
+ {
+ lambda[i] = lambdaArray[i][fepStateLeft]
+ + fracBetween * (lambdaArray[i][fepStateRight] - lambdaArray[i][fepStateLeft]);
+ }
+ return lambda;
+}
+
+} // namespace
+
+std::array<real, efptNR> currentLambdas(const int64_t step, const t_lambda& fepvals, const int currentLambdaState)
+{
+ if (fepvals.delta_lambda == 0)
+ {
+ if (currentLambdaState > -1)
+ {
+ return lambdasAtState(currentLambdaState, fepvals.all_lambda, fepvals.n_lambda);
+ }
+
+ if (fepvals.init_fep_state > -1)
+ {
+ return lambdasAtState(fepvals.init_fep_state, fepvals.all_lambda, fepvals.n_lambda);
+ }
+
+ std::array<real, efptNR> lambdas;
+ std::fill(std::begin(lambdas), std::end(lambdas), fepvals.init_lambda);
+ return lambdas;
+ }
+ const double globalLambda = currentGlobalLambda(step, fepvals.delta_lambda, fepvals.init_fep_state,
+ fepvals.init_lambda, fepvals.n_lambda);
+ return interpolatedLambdas(globalLambda, fepvals.all_lambda, fepvals.n_lambda);
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief Handling of free energy parameters
+ *
+ * \author Christian Blau <blau@kth.se>
+ * \ingroup module_mdlib
+ */
+
+#ifndef GMX_MDLIB_FREEENERGYPARAMETERS_H
+#define GMX_MDLIB_FREEENERGYPARAMETERS_H
+
+#include <array>
+
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/utility/real.h"
+
+struct t_lambda;
+
+namespace gmx
+{
+
+/*! \brief Evaluate the current lambdas
+ *
+ * \param[in] step the current simulation step
+ * \param[in] fepvals describing the lambda setup
+ * \param[in] currentLambdaState the lambda state to use to set the lambdas, -1 if not set
+ * \returns the current lambda-value array
+ */
+std::array<real, efptNR> currentLambdas(int64_t step, const t_lambda& fepvals, int currentLambdaState);
+
+} // namespace gmx
+
+#endif
}
}
-void setCurrentLambdasRerun(int64_t step,
- const t_lambda* fepvals,
- const t_trxframe* rerun_fr,
- const double* lam0,
- t_state* globalState)
-{
- GMX_RELEASE_ASSERT(globalState != nullptr,
- "setCurrentLambdasGlobalRerun should be called with a valid state object");
-
- if (rerun_fr->bLambda)
- {
- if (fepvals->delta_lambda == 0)
- {
- globalState->lambda[efptFEP] = rerun_fr->lambda;
- }
- else
- {
- /* find out between which two value of lambda we should be */
- real frac = step * fepvals->delta_lambda;
- int fep_state = static_cast<int>(std::floor(frac * fepvals->n_lambda));
- /* interpolate between this state and the next */
- /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
- frac = frac * fepvals->n_lambda - fep_state;
- for (int i = 0; i < efptNR; i++)
- {
- globalState->lambda[i] =
- lam0[i] + (fepvals->all_lambda[i][fep_state])
- + frac * (fepvals->all_lambda[i][fep_state + 1] - fepvals->all_lambda[i][fep_state]);
- }
- }
- }
- else if (rerun_fr->bFepState)
- {
- globalState->fep_state = rerun_fr->fep_state;
- for (int i = 0; i < efptNR; i++)
- {
- globalState->lambda[i] = fepvals->all_lambda[i][globalState->fep_state];
- }
- }
-}
-
-void setCurrentLambdasLocal(const int64_t step,
- const t_lambda* fepvals,
- const double* lam0,
- gmx::ArrayRef<real> lambda,
- const int currentFEPState)
-/* find the current lambdas. If rerunning, we either read in a state, or a lambda value,
- requiring different logic. */
-{
- if (fepvals->delta_lambda != 0)
- {
- /* find out between which two value of lambda we should be */
- real frac = step * fepvals->delta_lambda;
- if (fepvals->n_lambda > 0)
- {
- int fep_state = static_cast<int>(std::floor(frac * fepvals->n_lambda));
- /* interpolate between this state and the next */
- /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
- frac = frac * fepvals->n_lambda - fep_state;
- for (int i = 0; i < efptNR; i++)
- {
- lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state])
- + frac * (fepvals->all_lambda[i][fep_state + 1] - fepvals->all_lambda[i][fep_state]);
- }
- }
- else
- {
- for (int i = 0; i < efptNR; i++)
- {
- lambda[i] = lam0[i] + frac;
- }
- }
- }
- else
- {
- /* if < 0, fep_state was never defined, and we should not set lambda from the state */
- if (currentFEPState > -1)
- {
- for (int i = 0; i < efptNR; i++)
- {
- lambda[i] = fepvals->all_lambda[i][currentFEPState];
- }
- }
- }
-}
-
static void min_zero(int* n, int i)
{
if (i > 0 && (*n == 0 || i < *n))
#define GMX_MDLIB_MD_SUPPORT_H
#include "gromacs/mdlib/vcm.h"
+#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/timing/wallcycle.h"
struct gmx_ekindata_t;
struct t_forcerec;
struct t_grpopts;
struct t_inputrec;
-struct t_lambda;
struct t_nrnb;
class t_state;
struct t_trxframe;
//! \brief Allocate and initialize node-local state entries
void set_state_entries(t_state* state, const t_inputrec* ir, bool useModularSimulator);
-/* Set the lambda values in the global state from a frame read with rerun */
-void setCurrentLambdasRerun(int64_t step,
- const t_lambda* fepvals,
- const t_trxframe* rerun_fr,
- const double* lam0,
- t_state* globalState);
-
-/* Set the lambda values at each step of mdrun when they change */
-void setCurrentLambdasLocal(int64_t step,
- const t_lambda* fepvals,
- const double* lam0,
- gmx::ArrayRef<real> lambda,
- int currentFEPState);
-
-
/* Compute global variables during integration
*
* Coordinates x are needed for kinetic energy calculation with cosine accelation
constrtestrunners.cpp
ebin.cpp
energyoutput.cpp
+ freeenergyparameters.cpp
leapfrog.cpp
leapfrogtestdata.cpp
leapfrogtestrunners.cpp
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2020, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ *
+ * \brief Tests routines in freeenergyparameters.h .
+ *
+ * \author Christian Blau <blau@kth.se>
+ */
+
+#include "gmxpre.h"
+
+#include "gromacs/mdlib/freeenergyparameters.h"
+
+#include <gtest/gtest.h>
+
+#include "gromacs/mdtypes/inputrec.h"
+
+#include "testutils/testasserts.h"
+#include "testutils/testmatchers.h"
+
+
+namespace gmx
+{
+namespace test
+{
+namespace
+{
+
+
+/*! \brief Parameters that will vary from test to test.
+ */
+struct FreeEnergyParameterTestParameters
+{
+ //! current state of lambda in the simulation, -1 if not set
+ int currentLambdaState = -1;
+ //! Fractional value of lambda to start from, -1 if not set
+ double initLambda = -1;
+ //! The initial number of the state, -1 if not set
+ int initFepState = -1;
+ //! Change of lambda per time step (fraction of (0.1)
+ double deltaLambda = 0;
+ //! number of lambda entries
+ int nLambda = 0;
+ //! the current simulation step
+ int64_t step = 0;
+ //! the expected lambda at the current simulation step
+ std::array<real, efptNR> expectedLambdas = { -1, -1, -1, -1, -1, -1, -1 };
+};
+
+
+/*! \brief Sets of parameters on which to run the tests.
+ */
+const FreeEnergyParameterTestParameters freeEnergyParameterSets[] = {
+ // no parameters set at all
+ { -1, -1, -1, 0, 1, 0, { -1, -1, -1, -1, -1, -1, -1 } },
+ // setting current lambda state to 0, no other variables set, using eftpNR * [0.8] as lambda state matrix
+ { 0, -1, -1, 0, 1, 1, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ // setting current lambda state to 1, no other variables set, using eftpNR * [0.2,0.8] as lambda state matrix
+ { 1, -1, -1, 0, 2, 1, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ // test that non-zero deltaLambda trumps current lambda state
+ // setting current lambda state to 1, using deltaLambda 0.1 and setp 0 and eftpNR * [0.2,0.8] as lambda state matrix
+ { 1, -1, -1, 0.1, 2, 0, { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2 } },
+ // test that interpolating lambda values starting
+ // from lambda = 0, deltaLambda = 0.1, step = 10 results in values at end-state
+ { 1, 0, -1, 0.1, 2, 10, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ // interpolating half the way
+ { 1, 0, -1, 0.1, 2, 5, { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 } },
+ // starting from end-state 1 and move lambda half-way with negative deltaLambda = -0.1
+ { -1, -1, 1, -0.1, 2, 5, { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 } },
+ // starting from end-state 1 and move one step backwards with negative deltaLambda = -0.1
+ { -1, -1, 1, -0.1, 2, 1, { 0.74, 0.74, 0.74, 0.74, 0.74, 0.74, 0.74 } },
+ // three lambda states, the last two equal ([0.2,0.8,0.8]), move forward with deltaLambda = 0.1
+ { -1, -1, 0, 0.1, 3, 0, { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2 } },
+ { -1, -1, 0, 0.1, 3, 3, { 0.56, 0.56, 0.56, 0.56, 0.56, 0.56, 0.56 } },
+ { -1, -1, 0, 0.1, 3, 7, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ { -1, -1, 0, 0.1, 3, 8, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ { -1, -1, 0, 0.1, 3, 10, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ // three lambda states, the last two equal ([0.2,0.8,0.8]), move backwards
+ { -1, -1, 2, -0.1, 3, 1, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ { -1, -1, 2, -0.1, 3, 2, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ { -1, -1, 2, -0.1, 3, 3, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ { -1, -1, 2, -0.1, 3, 5, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ { -1, -1, 2, -0.1, 3, 7, { 0.56, 0.56, 0.56, 0.56, 0.56, 0.56, 0.56 } },
+ { -1, -1, 2, -0.1, 3, 8, { 0.44, 0.44, 0.44, 0.44, 0.44, 0.44, 0.44 } },
+ { -1, -1, 2, -0.1, 3, 10, { 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2 } },
+ // three lambda states, the last two equal ([0.2,0.8,0.8]), start in middle state, move backwards
+ { -1, -1, 1, -0.1, 3, 0, { 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8 } },
+ { -1, -1, 1, -0.1, 3, 2, { 0.56, 0.56, 0.56, 0.56, 0.56, 0.56, 0.56 } },
+ { -1, -1, 1, -0.1, 3, 3, { 0.44, 0.44, 0.44, 0.44, 0.44, 0.44, 0.44 } },
+};
+
+/*! \brief Test for setting free energy parameters.
+ */
+class FreeEnergyParameterTest : public ::testing::TestWithParam<FreeEnergyParameterTestParameters>
+{
+public:
+ //! Fills in the required FEP values for testing from the test parameters
+ t_lambda getFepVals()
+ {
+ t_lambda fepvals;
+ fepvals.init_fep_state = GetParam().initFepState;
+ fepvals.init_lambda = GetParam().initLambda;
+ fepvals.delta_lambda = GetParam().deltaLambda;
+ fepvals.all_lambda = getLambdaMatrix(GetParam().nLambda);
+ fepvals.n_lambda = GetParam().nLambda;
+ return fepvals;
+ }
+
+ /*! \brief construct a lambda matrix
+ *
+ * \param[in] nLambda
+ * \returns nLambda * eftpNR matrix with pre-defined values
+ */
+ double** getLambdaMatrix(int nLambda)
+ {
+ for (int i = 0; i < efptNR; ++i)
+ {
+ allLambda_[i] = defaultLambdaArrayForTest_[nLambda].data();
+ }
+ return allLambda_.data();
+ }
+
+private:
+ //! Construction aide for double ** matrix without snew
+ std::array<double*, efptNR> allLambda_;
+ //! a set of default lambda arrays for different lengths
+ std::vector<std::vector<double>> defaultLambdaArrayForTest_ = { {}, { 0.8 }, { 0.2, 0.8 }, { 0.2, 0.8, 0.8 } };
+};
+
+TEST_P(FreeEnergyParameterTest, CorrectLambdas)
+{
+ EXPECT_THAT(GetParam().expectedLambdas,
+ Pointwise(RealEq(defaultRealTolerance()),
+ currentLambdas(GetParam().step, getFepVals(), GetParam().currentLambdaState)));
+}
+
+INSTANTIATE_TEST_CASE_P(WithParameters, FreeEnergyParameterTest, ::testing::ValuesIn(freeEnergyParameterSets));
+
+} // namespace
+
+} // namespace test
+
+} // namespace gmx
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/freeenergyparameters.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
// will go away eventually.
t_inputrec* ir = inputrec;
int64_t step, step_rel;
- double t, t0 = ir->init_t, lam0[efptNR];
+ double t, t0 = ir->init_t;
gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
gmx_bool bNS = FALSE, bNStList, bStopCM, bFirstStep, bInitStep, bLastStep = FALSE;
gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
int* fep_state = MASTER(cr) ? &state_global->fep_state : nullptr;
gmx::ArrayRef<real> lambda = MASTER(cr) ? state_global->lambda : gmx::ArrayRef<real>();
- initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda, lam0);
+ initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda);
Update upd(*ir, deform);
const bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
if (ir->efep != efepNO || ir->bSimTemp)
{
/* find and set the current lambdas */
- setCurrentLambdasLocal(step, ir->fepvals, lam0, state->lambda, state->fep_state);
+ state->lambda = currentLambdas(step, *(ir->fepvals), state->fep_state);
bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/freeenergyparameters.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
{
t_inputrec* ir = inputrec;
int64_t step, step_rel;
- double t, lam0[efptNR];
+ double t;
bool isLastStep = false;
bool doFreeEnergyPerturbation = false;
unsigned int force_flags;
nonConstGlobalTopology->intermolecularExclusionGroup = genQmmmIndices(*top_global);
}
- initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda, lam0);
+ initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda);
const bool simulationsShareState = false;
gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider,
if (ir->efep != efepNO)
{
- setCurrentLambdasLocal(step, ir->fepvals, lam0, state->lambda, state->fep_state);
+ state->lambda = currentLambdas(step, *(ir->fepvals), state_global->fep_state);
}
if (MASTER(cr))
}
int* fep_state = MASTER(cr) ? &state_global->fep_state : nullptr;
gmx::ArrayRef<real> lambda = MASTER(cr) ? state_global->lambda : gmx::ArrayRef<real>();
- initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda, nullptr);
+ initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda);
if (ir->eI == eiNM)
{
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/freeenergyparameters.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
// will go away eventually.
t_inputrec* ir = inputrec;
int64_t step, step_rel;
- double t, lam0[efptNR];
+ double t;
bool isLastStep = false;
bool doFreeEnergyPerturbation = false;
unsigned int force_flags;
}
int* fep_state = MASTER(cr) ? &state_global->fep_state : nullptr;
gmx::ArrayRef<real> lambda = MASTER(cr) ? state_global->lambda : gmx::ArrayRef<real>();
- initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda, lam0);
+ initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda);
const bool simulationsShareState = false;
gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider,
mdModulesNotifier, ir, top_global, oenv, wcycle,
if (ir->efep != efepNO && MASTER(cr))
{
- setCurrentLambdasRerun(step, ir->fepvals, &rerun_fr, lam0, state_global);
+ if (rerun_fr.bLambda)
+ {
+ ir->fepvals->init_lambda = rerun_fr.lambda;
+ }
+ else
+ {
+ if (rerun_fr.bFepState)
+ {
+ state->fep_state = rerun_fr.fep_state;
+ }
+ }
+
+ state_global->lambda = currentLambdas(step, *(ir->fepvals), state->fep_state);
}
if (MASTER(cr))
}
}
-void initialize_lambdas(FILE* fplog,
- const t_inputrec& ir,
- bool isMaster,
- int* fep_state,
- gmx::ArrayRef<real> lambda,
- double* lam0)
+void initialize_lambdas(FILE* fplog, const t_inputrec& ir, bool isMaster, int* fep_state, gmx::ArrayRef<real> lambda)
{
/* TODO: Clean up initialization of fep_state and lambda in
t_state. This function works, but could probably use a logic
{
lambda[i] = thisLambda;
}
- if (lam0 != nullptr)
- {
- lam0[i] = thisLambda;
- }
}
if (ir.bSimTemp)
{
void printLambdaStateToLog(FILE* fplog, gmx::ArrayRef<real> lambda, bool isInitialOutput);
-/*! \brief Fills fep_state, lambda, and lam0 if needed
+/*! \brief Fills fep_state and lambda if needed
*
- * If FEP or simulated tempering is in use:
- *
- * fills non-null lam0 with the initial lambda values, and
- * on master rank fills fep_state and lambda.
+ * If FEP or simulated tempering is in use, fills fep_state
+ * and lambda on master rank.
*
* Reports the initial lambda state to the log file. */
-void initialize_lambdas(FILE* fplog,
- const t_inputrec& ir,
- bool isMaster,
- int* fep_state,
- gmx::ArrayRef<real> lambda,
- double* lam0);
+void initialize_lambdas(FILE* fplog, const t_inputrec& ir, bool isMaster, int* fep_state, gmx::ArrayRef<real> lambda);
#endif
#include "freeenergyperturbationdata.h"
+#include "gromacs/mdlib/freeenergyparameters.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdtypes/inputrec.h"
FreeEnergyPerturbationData::FreeEnergyPerturbationData(FILE* fplog, const t_inputrec* inputrec, MDAtoms* mdAtoms) :
element_(std::make_unique<Element>(this, inputrec->fepvals->delta_lambda)),
lambda_(),
- lambda0_(),
currentFEPState_(0),
fplog_(fplog),
inputrec_(inputrec),
mdAtoms_(mdAtoms)
{
lambda_.fill(0);
- lambda0_.fill(0);
- initialize_lambdas(fplog_, *inputrec_, true, ¤tFEPState_, lambda_, lambda0_.data());
+ initialize_lambdas(fplog_, *inputrec_, true, ¤tFEPState_, lambda_);
update_mdatoms(mdAtoms_->mdatoms(), lambda_[efptMASS]);
}
void FreeEnergyPerturbationData::updateLambdas(Step step)
{
// at beginning of step (if lambdas change...)
- setCurrentLambdasLocal(step, inputrec_->fepvals, lambda0_.data(), lambda_, currentFEPState_);
+ lambda_ = currentLambdas(step, *(inputrec_->fepvals), currentFEPState_);
update_mdatoms(mdAtoms_->mdatoms(), lambda_[efptMASS]);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff