{
}
-std::pair<std::vector<double>, std::vector<double>> ForeignLambdaTerms::getTerms(t_commrec* cr)
+std::pair<std::vector<double>, std::vector<double>> ForeignLambdaTerms::getTerms(t_commrec* cr) const
{
+ GMX_RELEASE_ASSERT(finalizedPotentialContributions_,
+ "The object needs to be finalized before calling getTerms");
+
std::vector<double> data(2 * numLambdas_);
for (int i = 0; i < numLambdas_; i++)
{
{
std::fill(energies_.begin(), energies_.end(), 0.0);
std::fill(dhdl_.begin(), dhdl_.end(), 0.0);
+ finalizedPotentialContributions_ = false;
}
gmx_enerdata_t::gmx_enerdata_t(int numEnergyGroups, int numFepLambdas) :
}
}
-void accumulatePotentialEnergies(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda* fepvals)
+void ForeignLambdaTerms::addConstantDhdl(const double dhdl)
{
- sum_epot(enerd->grpp, enerd->term);
+ for (double& foreignDhdl : dhdl_)
+ {
+ foreignDhdl += dhdl;
+ }
+}
- if (fepvals == nullptr)
+void ForeignLambdaTerms::finalizePotentialContributions(gmx::ArrayRef<const double> dvdlLinear,
+ gmx::ArrayRef<const real> lambda,
+ const t_lambda& fepvals)
+{
+ if (finalizedPotentialContributions_)
{
return;
}
- // Accumulate the foreign lambda terms
-
double dvdl_lin = 0;
for (int i = 0; i < efptNR; i++)
{
- dvdl_lin += enerd->dvdl_lin[i];
+ dvdl_lin += dvdlLinear[i];
}
- enerd->foreignLambdaTerms.addConstantDhdl(dvdl_lin);
-
- set_dvdl_output(enerd, *fepvals);
+ addConstantDhdl(dvdl_lin);
/* Sum the foreign lambda energy difference contributions.
* Note that here we only add the potential energy components.
* The constraint and kinetic energy components are add after integration
* by sum_dhdl().
*/
- for (int i = 0; i < fepvals->n_lambda; i++)
+ for (int i = 0; i < fepvals.n_lambda; i++)
{
/* note we are iterating over fepvals here!
For the current lam, dlam = 0 automatically,
for (gmx::index j = 0; j < lambda.ssize(); j++)
{
/* Note that this loop is over all dhdl components, not just the separated ones */
- const double dlam = fepvals->all_lambda[j][i] - lambda[j];
+ const double dlam = fepvals.all_lambda[j][i] - lambda[j];
- enerpart_lambda += dlam * enerd->dvdl_lin[j];
+ enerpart_lambda += dlam * dvdlLinear[j];
}
- enerd->foreignLambdaTerms.accumulate(1 + i, enerpart_lambda, 0);
+ accumulate(1 + i, enerpart_lambda, 0);
}
+
+ finalizedPotentialContributions_ = true;
}
-void accumulateKineticLambdaComponents(gmx_enerdata_t* enerd,
- gmx::ArrayRef<const real> lambda,
- const t_lambda& fepvals)
+void accumulatePotentialEnergies(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda* fepvals)
{
- if (fepvals.separate_dvdl[efptBONDED])
- {
- enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
- }
- else
+ sum_epot(enerd->grpp, enerd->term);
+
+ if (fepvals)
{
- enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
+ set_dvdl_output(enerd, *fepvals);
+
+ enerd->foreignLambdaTerms.finalizePotentialContributions(enerd->dvdl_lin, lambda, *fepvals);
}
+}
+
+void ForeignLambdaTerms::accumulateKinetic(int listIndex, double energy, double dhdl)
+{
+ energies_[listIndex] += energy;
+ dhdl_[listIndex] += dhdl;
+}
+
+void ForeignLambdaTerms::finalizeKineticContributions(gmx::ArrayRef<const real> energyTerms,
+ const double dhdlMass,
+ gmx::ArrayRef<const real> lambda,
+ const t_lambda& fepvals)
+{
+ // Add perturbed mass contributions
+ addConstantDhdl(dhdlMass);
- // Treat current lambda, which only needs a dV/dl contribution
- enerd->foreignLambdaTerms.accumulate(0, 0.0, enerd->term[F_DVDL_CONSTR]);
+ // Treat current lambda, the deltaH contribution is 0 as delta-lambda=0 for the current lambda
+ accumulateKinetic(0, 0.0, energyTerms[F_DVDL_CONSTR]);
if (!fepvals.separate_dvdl[efptMASS])
{
- enerd->foreignLambdaTerms.accumulate(0, 0.0, enerd->term[F_DKDL]);
+ accumulateKinetic(0, 0.0, energyTerms[F_DKDL]);
}
for (int i = 0; i < fepvals.n_lambda; i++)
*/
const int lambdaIndex = (fepvals.separate_dvdl[efptBONDED] ? efptBONDED : efptFEP);
const double dlam = fepvals.all_lambda[lambdaIndex][i] - lambda[lambdaIndex];
- enerd->foreignLambdaTerms.accumulate(1 + i, dlam * enerd->term[F_DVDL_CONSTR],
- enerd->term[F_DVDL_CONSTR]);
+ accumulateKinetic(1 + i, dlam * energyTerms[F_DVDL_CONSTR], energyTerms[F_DVDL_CONSTR]);
if (!fepvals.separate_dvdl[efptMASS])
{
const double dlam = fepvals.all_lambda[efptMASS][i] - lambda[efptMASS];
- enerd->foreignLambdaTerms.accumulate(1 + i, dlam * enerd->term[F_DKDL], enerd->term[F_DKDL]);
+ accumulateKinetic(1 + i, dlam * energyTerms[F_DKDL], energyTerms[F_DKDL]);
}
}
+}
+
+void accumulateKineticLambdaComponents(gmx_enerdata_t* enerd,
+ gmx::ArrayRef<const real> lambda,
+ const t_lambda& fepvals)
+{
+ if (fepvals.separate_dvdl[efptBONDED])
+ {
+ enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
+ }
+ else
+ {
+ enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
+ }
+
+ enerd->foreignLambdaTerms.finalizeKineticContributions(enerd->term, enerd->dvdl_lin[efptMASS],
+ lambda, fepvals);
/* The constrain contribution is now included in other terms, so clear it */
enerd->term[F_DVDL_CONSTR] = 0;
}
-
void reset_foreign_enerdata(gmx_enerdata_t* enerd)
{
int i, j;
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/topology/idef.h"
#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/real.h"
struct t_commrec;
+struct t_lambda;
// The non-bonded energy terms accumulated for energy group pairs
enum
/*! \brief Returns a list of partial energies, the part which depends on lambda),
* current lambda in entry 0, foreign lambda i in entry 1+i
+ *
+ * Note: the potential terms needs to be finalized before calling this method.
*/
- gmx::ArrayRef<double> energies() { return energies_; }
+ gmx::ArrayRef<double> energies()
+ {
+ GMX_ASSERT(finalizedPotentialContributions_, "Should be finalized");
+ return energies_;
+ }
/*! \brief Returns a list of partial energies, the part which depends on lambda),
* current lambda in entry 0, foreign lambda i in entry 1+i
+ *
+ * Note: the potential terms needs to be finalized before calling this method.
*/
- gmx::ArrayRef<const double> energies() const { return energies_; }
+ gmx::ArrayRef<const double> energies() const
+ {
+ GMX_ASSERT(finalizedPotentialContributions_, "Should be finalized");
+ return energies_;
+ }
- /*! \brief Adds an energy and dH/dl constribution to lambda list index \p listIndex
+ /*! \brief Adds an energy and dV/dl constribution to lambda list index \p listIndex
*
* This should only be used for terms with non-linear dependence on lambda
* The value passed as listIndex should be 0 for the current lambda
* and 1+i for foreign lambda index i.
*/
- void accumulate(int listIndex, double energy, double dhdl)
+ void accumulate(int listIndex, double energy, double dvdl)
{
+ GMX_ASSERT(!finalizedPotentialContributions_,
+ "Can only accumulate with an unfinalized object");
+
energies_[listIndex] += energy;
- dhdl_[listIndex] += dhdl;
+ dhdl_[listIndex] += dvdl;
}
- /*! \brief Add a dH/dl contribution that does not depend on lambda to all foreign dH/dl terms
+ /*! \brief Finalizes the potential (non-kinetic) terms
*
- * Note: this should not be called directly for energy terms that depend linearly on lambda,
- * as those are added automatically through the accumulated dvdl_lin term in gmx_enerdata_t.
+ * Note: This can be called multiple times during the same force calculations
+ * without affecting the results.
+ *
+ * \param[in] dvdlLinear List of dV/dlambda contributions of size efptNR with depend linearly on lambda
+ * \param[in] lambda Lambda values for the efptNR contribution types
+ * \param[in] fepvals Free-energy parameters
*/
- void addConstantDhdl(double dhdl)
- {
- for (double& foreignDhdl : dhdl_)
- {
- foreignDhdl += dhdl;
- }
- }
+ void finalizePotentialContributions(gmx::ArrayRef<const double> dvdlLinear,
+ gmx::ArrayRef<const real> lambda,
+ const t_lambda& fepvals);
+
+ /* !\brief Accumulates the kinetic and constraint free-energy contributions
+ *
+ * \param[in] energyTerms List of energy terms, pass \p term in \p gmx_enerdata_t
+ * \param[in] dhdlMass The mass dependent contribution to dH/dlambda
+ * \param[in] lambda Lambda values for the efptNR contribution types
+ * \param[in] fepvals Free-energy parameters
+ */
+ void finalizeKineticContributions(gmx::ArrayRef<const real> energyTerms,
+ double dhdlMass,
+ gmx::ArrayRef<const real> lambda,
+ const t_lambda& fepvals);
/*! \brief Returns a pair of lists of deltaH and dH/dlambda
*
* Both lists are of size numLambdas() and are indexed with the lambda index.
- * The returned lists are valid until the next call to this method.
+ *
+ * Note: should only be called after the object has been finalized by a call to
+ * accumulateLinearPotentialComponents() (is asserted).
*
* \param[in] cr Communication record, used to reduce the terms when !=nullptr
*/
- std::pair<std::vector<double>, std::vector<double>> getTerms(t_commrec* cr);
+ std::pair<std::vector<double>, std::vector<double>> getTerms(t_commrec* cr) const;
//! Sets all terms to 0
void zeroAllTerms();
private:
+ //! As accumulate(), but for kinetic contributions
+ void accumulateKinetic(int listIndex, double energy, double dhdl);
+
+ //! Add a dH/dl contribution that does not depend on lambda to all foreign dH/dl terms
+ void addConstantDhdl(double dhdl);
+
//! The number of foreign lambdas
int numLambdas_;
//! Storage for foreign lambda energies
std::vector<double> energies_;
//! Storage for foreign lambda dH/dlambda
std::vector<double> dhdl_;
+ //! Tells whether all potential energy contributions have been accumulated
+ bool finalizedPotentialContributions_ = false;
};
//! Struct for accumulating all potential energy terms and some kinetic energy terms
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff