{
int flags;
const struct t_blocka* exclusions;
- real* lambda;
+ const real* lambda;
real* dvdl;
/* pointers to tables */
calc_listed_lambda(idef, threading_.get(), x, fr, pbc, forceBufferLambda_,
shiftForceBufferLambda_, &(enerd->foreign_grpp), enerd->foreign_term,
dvdl, nrnb, lam_i, md, &fcdata, global_atom_index);
- sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
+ sum_epot(enerd->foreign_grpp, enerd->foreign_term);
enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
for (int j = 0; j < efptNR; j++)
{
return t;
}
-void sum_epot(gmx_grppairener_t* grpp, real* epot)
+void sum_epot(const gmx_grppairener_t& grpp, real* epot)
{
int i;
/* Accumulate energies */
- epot[F_COUL_SR] = sum_v(grpp->nener, grpp->ener[egCOULSR]);
- epot[F_LJ] = sum_v(grpp->nener, grpp->ener[egLJSR]);
- epot[F_LJ14] = sum_v(grpp->nener, grpp->ener[egLJ14]);
- epot[F_COUL14] = sum_v(grpp->nener, grpp->ener[egCOUL14]);
+ epot[F_COUL_SR] = sum_v(grpp.nener, grpp.ener[egCOULSR]);
+ epot[F_LJ] = sum_v(grpp.nener, grpp.ener[egLJSR]);
+ epot[F_LJ14] = sum_v(grpp.nener, grpp.ener[egLJ14]);
+ epot[F_COUL14] = sum_v(grpp.nener, grpp.ener[egCOUL14]);
/* lattice part of LR doesnt belong to any group
* and has been added earlier
*/
- epot[F_BHAM] = sum_v(grpp->nener, grpp->ener[egBHAMSR]);
+ epot[F_BHAM] = sum_v(grpp.nener, grpp.ener[egBHAMSR]);
epot[F_EPOT] = 0;
for (i = 0; (i < F_EPOT); i++)
}
}
-void sum_dhdl(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda& fepvals)
+/* Adds computed dV/dlambda contributions to the requested outputs
+ *
+ * Also adds the dispersion correction dV/dlambda to the VdW term.
+ * Note that kinetic energy and constraint contributions are handled in sum_dkdl()
+ */
+static void sum_dvdl(gmx_enerdata_t* enerd, const t_lambda& fepvals)
{
- int index;
-
- enerd->dvdl_lin[efptVDW] += enerd->term[F_DVDL_VDW]; /* include dispersion correction */
+ // Add dispersion correction to the VdW component
+ enerd->dvdl_lin[efptVDW] += enerd->term[F_DVDL_VDW];
for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++)
{
if (fepvals.separate_dvdl[i])
{
/* could this be done more readably/compactly? */
+ int index;
switch (i)
{
case (efptMASS): index = F_DKDL; break;
}
}
}
+}
+
+void accumulatePotentialEnergies(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda* fepvals)
+{
+ sum_epot(enerd->grpp, enerd->term);
+
+ if (fepvals)
+ {
+ /* Note that sum_dvdl() adds the dispersion correction enerd->dvdl_lin[efptVDW],
+ * so sum_dvdl() should be called before computing the foreign lambda energy differences.
+ */
+ sum_dvdl(enerd, *fepvals);
+
+ /* 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++)
+ {
+ /* note we are iterating over fepvals here!
+ For the current lam, dlam = 0 automatically,
+ so we don't need to add anything to the
+ enerd->enerpart_lambda[0] */
+
+ /* we don't need to worry about dvdl_lin contributions to dE at
+ current lambda, because the contributions to the current
+ lambda are automatically zeroed */
+
+ double& enerpart_lambda = enerd->enerpart_lambda[i + 1];
+
+ 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];
+ enerpart_lambda += dlam * enerd->dvdl_lin[j];
+ }
+ }
+ }
+}
+
+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];
so we don't need to add anything to the
enerd->enerpart_lambda[0] */
- /* we don't need to worry about dvdl_lin contributions to dE at
- current lambda, because the contributions to the current
- lambda are automatically zeroed */
-
double& enerpart_lambda = enerd->enerpart_lambda[i + 1];
- 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];
-
- enerpart_lambda += dlam * enerd->dvdl_lin[j];
+ /* Note that potential energy terms have been added by sum_epot() -> sum_dvdl() */
- /* Constraints can not be evaluated at foreign lambdas, so we add
- * a linear extrapolation. This is an approximation, but usually
- * quite accurate since constraints change little between lambdas.
- */
- if ((j == efptBONDED && fepvals.separate_dvdl[efptBONDED])
- || (j == efptFEP && !fepvals.separate_dvdl[efptBONDED]))
- {
- enerpart_lambda += dlam * enerd->term[F_DVDL_CONSTR];
- }
+ /* Constraints can not be evaluated at foreign lambdas, so we add
+ * a linear extrapolation. This is an approximation, but usually
+ * quite accurate since constraints change little between lambdas.
+ */
+ const int lambdaIndex = (fepvals.separate_dvdl[efptBONDED] ? efptBONDED : efptFEP);
+ const double dlam = fepvals.all_lambda[lambdaIndex][i] - lambda[lambdaIndex];
+ enerpart_lambda += dlam * enerd->term[F_DVDL_CONSTR];
- if (j == efptMASS && !fepvals.separate_dvdl[j])
- {
- enerpart_lambda += dlam * enerd->term[F_DKDL];
- }
-
- if (debug)
- {
- fprintf(debug, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
- fepvals.all_lambda[j][i], efpt_names[j],
- enerpart_lambda - enerd->enerpart_lambda[0], dlam, enerd->dvdl_lin[j]);
- }
+ if (!fepvals.separate_dvdl[efptMASS])
+ {
+ const double dlam = fepvals.all_lambda[efptMASS][i] - lambda[efptMASS];
+ enerpart_lambda += dlam * enerd->term[F_DKDL];
}
}
struct gmx_enerdata_t;
struct gmx_grppairener_t;
+struct t_fepvals;
struct t_lambda;
void reset_foreign_enerdata(gmx_enerdata_t* enerd);
void reset_enerdata(gmx_enerdata_t* enerd);
/* Resets the energy data */
-void sum_epot(gmx_grppairener_t* grpp, real* epot);
-/* Locally sum the non-bonded potential energy terms */
+/*! \brief Sums energy group pair contributions into epot */
+void sum_epot(const gmx_grppairener_t& grpp, real* epot);
-void sum_dhdl(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda& fepvals);
-/* Sum the free energy contributions */
+/*! \brief Accumulates potential energy contributions to obtain final potential energies
+ *
+ * Accumulates energy group pair contributions into the output energy components
+ * and sums all potential energies into the total potential energy term.
+ * With free-energy also computes the foreign lambda potential energy differences.
+ *
+ * \param[in,out] enerd Energy data with components to sum and to accumulate into
+ * \param[in] lambda Vector of free-energy lambdas
+ * \param[in] fepvals Foreign lambda energy differences, only summed with !=nullptr
+ */
+void accumulatePotentialEnergies(gmx_enerdata_t* enerd,
+ gmx::ArrayRef<const real> lambda,
+ const t_lambda* fepvals);
+
+/*! \brief Accumulates kinetic and constraint contributions to dH/dlambda and foreign energies */
+void accumulateKineticLambdaComponents(gmx_enerdata_t* enerd,
+ gmx::ArrayRef<const real> lambda,
+ const t_lambda& fepvals);
#endif
tensor vir_force,
const t_mdatoms* mdatoms,
gmx_enerdata_t* enerd,
- gmx::ArrayRef<real> lambda,
+ gmx::ArrayRef<const real> lambda,
t_forcerec* fr,
gmx::MdrunScheduleWorkload* runScheduleWork,
gmx::VirtualSitesHandler* vsite,
const matrix box,
gmx::ArrayRef<const gmx::RVec> x,
const t_mdatoms* mdatoms,
- real* lambda,
+ gmx::ArrayRef<const real> lambda,
const StepWorkload& stepWork,
gmx::ForceWithVirial* forceWithVirial,
gmx_enerdata_t* enerd,
if (inputrec->bPull && pull_have_potential(pull_work))
{
pull_potential_wrapper(cr, inputrec, box, x, forceWithVirial, mdatoms, enerd, pull_work,
- lambda, t, wcycle);
+ lambda.data(), t, wcycle);
if (awh)
{
tensor vir_force,
const t_mdatoms* mdatoms,
gmx_enerdata_t* enerd,
- gmx::ArrayRef<real> lambda,
+ gmx::ArrayRef<const real> lambda,
t_forcerec* fr,
gmx::MdrunScheduleWorkload* runScheduleWork,
gmx::VirtualSitesHandler* vsite,
nbv->dispatchFreeEnergyKernel(InteractionLocality::Local, fr,
as_rvec_array(x.unpaddedArrayRef().data()),
&forceOut.forceWithShiftForces(), *mdatoms, inputrec->fepvals,
- lambda.data(), enerd, stepWork, nrnb);
+ lambda, enerd, stepWork, nrnb);
if (havePPDomainDecomposition(cr))
{
nbv->dispatchFreeEnergyKernel(InteractionLocality::NonLocal, fr,
as_rvec_array(x.unpaddedArrayRef().data()),
&forceOut.forceWithShiftForces(), *mdatoms,
- inputrec->fepvals, lambda.data(), enerd, stepWork, nrnb);
+ inputrec->fepvals, lambda, enerd, stepWork, nrnb);
}
}
wallcycle_stop(wcycle, ewcFORCE);
computeSpecialForces(fplog, cr, inputrec, awh, enforcedRotation, imdSession, pull_work, step, t,
- wcycle, fr->forceProviders, box, x.unpaddedArrayRef(), mdatoms, lambda.data(),
+ wcycle, fr->forceProviders, box, x.unpaddedArrayRef(), mdatoms, lambda,
stepWork, &forceOut.forceWithVirial(), enerd, ed, stepWork.doNeighborSearch);
if (stepWork.computeEnergy)
{
- /* Sum the potential energy terms from group contributions */
- sum_epot(&(enerd->grpp), enerd->term);
+ /* Compute the final potential energy terms */
+ accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals);
+ ;
if (!EI_TPI(inputrec->eI))
{
{
saved_conserved_quantity -= enerd->term[F_DISPCORR];
}
- /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
+ /* sum up the foreign kinetic energy and dK/dl terms for vv. currently done every step so that dhdl is correct in the .edr */
if (ir->efep != efepNO)
{
- sum_dhdl(enerd, state->lambda, *ir->fepvals);
+ accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
}
}
if (ir->efep != efepNO && !EI_VV(ir->eI))
{
- /* Sum up the foreign energy and dhdl terms for md and sd.
- Currently done every step so that dhdl is correct in the .edr */
- sum_dhdl(enerd, state->lambda, *ir->fepvals);
+ /* Sum up the foreign energy and dK/dl terms for md and sd.
+ Currently done every step so that dH/dl is correct in the .edr */
+ accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
}
update_pcouple_after_coordinates(fplog, step, ir, mdatoms, pres, force_vir, shake_vir,
{
/* Sum up the foreign energy and dhdl terms for md and sd.
Currently done every step so that dhdl is correct in the .edr */
- sum_dhdl(enerd, state->lambda, *ir->fepvals);
+ accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
}
/* Output stuff */
clear_mat(ekin);
enerd->term[F_PRES] = calc_pres(fr->pbcType, inputrec->nwall, ems->s.box, ekin, vir, pres);
- sum_dhdl(enerd, ems->s.lambda, *inputrec->fepvals);
+ if (inputrec->efep != efepNO)
+ {
+ accumulateKineticLambdaComponents(enerd, ems->s.lambda, *inputrec->fepvals);
+ }
if (EI_ENERGY_MINIMIZATION(inputrec->eI))
{
but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
- if (ir->efep != efepNO)
- {
- /* Sum up the foreign energy and dhdl terms for md and sd.
- Currently done every step so that dhdl is correct in the .edr */
- sum_dhdl(enerd, state->lambda, *ir->fepvals);
- }
-
/* Output stuff */
if (MASTER(cr))
{
sf_dir += md->massT[i] * norm2(shfc->acc_dir[i]);
}
}
- sum_epot(&(enerd->grpp), enerd->term);
+ accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals);
Epot[Min] = enerd->term[F_EPOT];
df[Min] = rms_force(cr, forceWithPadding[Min].paddedArrayRef(), shells, nflexcon, &sf_dir, &Epot[Min]);
wcycle, top, box, posWithPadding[Try], hist, forceWithPadding[Try], force_vir, md,
enerd, lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, shellfc_flags,
ddBalanceRegionHandler);
- sum_epot(&(enerd->grpp), enerd->term);
+ accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals);
if (gmx_debug_at)
{
pr_rvecs(debug, 0, "RELAX: force[Min]", as_rvec_array(force[Min].data()), homenr);
}
if (freeEnergyPerturbationElement_)
{
- sum_dhdl(enerd_, freeEnergyPerturbationElement_->constLambdaView(), *inputrec_->fepvals);
+ accumulateKineticLambdaComponents(enerd_, freeEnergyPerturbationElement_->constLambdaView(),
+ *inputrec_->fepvals);
dummyLegacyState_.fep_state = freeEnergyPerturbationElement_->currentFEPState();
}
if (parrinelloRahmanBarostat_)
gmx::ForceWithShiftForces* forceWithShiftForces,
const t_mdatoms& mdatoms,
t_lambda* fepvals,
- real* lambda,
+ gmx::ArrayRef<real const> lambda,
gmx_enerdata_t* enerd,
const gmx::StepWorkload& stepWork,
t_nrnb* nrnb)
nb_kernel_data_t kernel_data;
real dvdl_nb[efptNR] = { 0 };
kernel_data.flags = donb_flags;
- kernel_data.lambda = lambda;
+ kernel_data.lambda = lambda.data();
kernel_data.dvdl = dvdl_nb;
kernel_data.energygrp_elec = enerd->grpp.ener[egCOULSR].data();
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
- sum_epot(&(enerd->foreign_grpp), enerd->foreign_term);
+ sum_epot(enerd->foreign_grpp, enerd->foreign_term);
enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];
enerd->dhdlLambda[i] += dvdl_nb[efptVDW] + dvdl_nb[efptCOUL];
}
gmx::ForceWithShiftForces* forceWithShiftForces,
const t_mdatoms& mdatoms,
t_lambda* fepvals,
- real* lambda,
+ gmx::ArrayRef<const real> lambda,
gmx_enerdata_t* enerd,
const gmx::StepWorkload& stepWork,
t_nrnb* nrnb);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff