From af4703504734dd000872d1484e309220cecd8f08 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Mon, 6 Jul 2020 15:14:31 +0000 Subject: [PATCH] Move foreign potential energy accumulation The potential energy contributions to the foreign lambda Hamiltonian differences were summed in sum_dhdl() which meant that the foreign energy differences were not complete when do_force() returns. Now there are summed in accumulatePotentialEnergies() which is called instead of sum_epot(). Also consistently changed the lambda vector to ArrayRef in the force routines and made it const. --- src/gromacs/gmxlib/nonbonded/nb_kernel.h | 2 +- src/gromacs/listed_forces/listed_forces.cpp | 2 +- src/gromacs/mdlib/enerdata_utils.cpp | 112 +++++++++++------- src/gromacs/mdlib/enerdata_utils.h | 25 +++- src/gromacs/mdlib/force.h | 2 +- src/gromacs/mdlib/sim_util.cpp | 17 +-- src/gromacs/mdrun/md.cpp | 10 +- src/gromacs/mdrun/mimic.cpp | 2 +- src/gromacs/mdrun/minimize.cpp | 5 +- src/gromacs/mdrun/rerun.cpp | 7 -- src/gromacs/mdrun/shellfc.cpp | 4 +- .../modularsimulator/energyelement.cpp | 3 +- src/gromacs/nbnxm/kerneldispatch.cpp | 6 +- src/gromacs/nbnxm/nbnxm.h | 2 +- 14 files changed, 123 insertions(+), 76 deletions(-) diff --git a/src/gromacs/gmxlib/nonbonded/nb_kernel.h b/src/gromacs/gmxlib/nonbonded/nb_kernel.h index 2cf148730d..ff75e9159a 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_kernel.h +++ b/src/gromacs/gmxlib/nonbonded/nb_kernel.h @@ -53,7 +53,7 @@ typedef struct { int flags; const struct t_blocka* exclusions; - real* lambda; + const real* lambda; real* dvdl; /* pointers to tables */ diff --git a/src/gromacs/listed_forces/listed_forces.cpp b/src/gromacs/listed_forces/listed_forces.cpp index ea2c9a0a9e..e88f050fbd 100644 --- a/src/gromacs/listed_forces/listed_forces.cpp +++ b/src/gromacs/listed_forces/listed_forces.cpp @@ -718,7 +718,7 @@ void ListedForces::calculate(struct gmx_wallcycle* wcycle, 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++) { diff --git a/src/gromacs/mdlib/enerdata_utils.cpp b/src/gromacs/mdlib/enerdata_utils.cpp index c5e1bd33d8..4d9a758bb1 100644 --- a/src/gromacs/mdlib/enerdata_utils.cpp +++ b/src/gromacs/mdlib/enerdata_utils.cpp @@ -66,20 +66,20 @@ static real sum_v(int n, gmx::ArrayRef v) 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++) @@ -91,11 +91,15 @@ void sum_epot(gmx_grppairener_t* grpp, real* epot) } } -void sum_dhdl(gmx_enerdata_t* enerd, gmx::ArrayRef 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++) { @@ -108,6 +112,7 @@ void sum_dhdl(gmx_enerdata_t* enerd, gmx::ArrayRef lambda, const t_l if (fepvals.separate_dvdl[i]) { /* could this be done more readably/compactly? */ + int index; switch (i) { case (efptMASS): index = F_DKDL; break; @@ -134,7 +139,52 @@ void sum_dhdl(gmx_enerdata_t* enerd, gmx::ArrayRef lambda, const t_l } } } +} + +void accumulatePotentialEnergies(gmx_enerdata_t* enerd, gmx::ArrayRef 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 lambda, + const t_lambda& fepvals) +{ if (fepvals.separate_dvdl[efptBONDED]) { enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR]; @@ -151,40 +201,22 @@ void sum_dhdl(gmx_enerdata_t* enerd, gmx::ArrayRef lambda, const t_l 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]; } } diff --git a/src/gromacs/mdlib/enerdata_utils.h b/src/gromacs/mdlib/enerdata_utils.h index bfefd04bfe..5e0b066e0c 100644 --- a/src/gromacs/mdlib/enerdata_utils.h +++ b/src/gromacs/mdlib/enerdata_utils.h @@ -45,6 +45,7 @@ struct gmx_enerdata_t; struct gmx_grppairener_t; +struct t_fepvals; struct t_lambda; void reset_foreign_enerdata(gmx_enerdata_t* enerd); @@ -56,10 +57,26 @@ void reset_dvdl_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 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 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 lambda, + const t_lambda& fepvals); #endif diff --git a/src/gromacs/mdlib/force.h b/src/gromacs/mdlib/force.h index 81312671b9..3db3f96db4 100644 --- a/src/gromacs/mdlib/force.h +++ b/src/gromacs/mdlib/force.h @@ -91,7 +91,7 @@ void do_force(FILE* log, tensor vir_force, const t_mdatoms* mdatoms, gmx_enerdata_t* enerd, - gmx::ArrayRef lambda, + gmx::ArrayRef lambda, t_forcerec* fr, gmx::MdrunScheduleWorkload* runScheduleWork, gmx::VirtualSitesHandler* vsite, diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 9ace847885..5159f2438d 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -610,7 +610,7 @@ static void computeSpecialForces(FILE* fplog, const matrix box, gmx::ArrayRef x, const t_mdatoms* mdatoms, - real* lambda, + gmx::ArrayRef lambda, const StepWorkload& stepWork, gmx::ForceWithVirial* forceWithVirial, gmx_enerdata_t* enerd, @@ -632,7 +632,7 @@ static void computeSpecialForces(FILE* fplog, 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) { @@ -1011,7 +1011,7 @@ void do_force(FILE* fplog, tensor vir_force, const t_mdatoms* mdatoms, gmx_enerdata_t* enerd, - gmx::ArrayRef lambda, + gmx::ArrayRef lambda, t_forcerec* fr, gmx::MdrunScheduleWorkload* runScheduleWork, gmx::VirtualSitesHandler* vsite, @@ -1525,14 +1525,14 @@ void do_force(FILE* fplog, 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); } } @@ -1579,7 +1579,7 @@ void do_force(FILE* fplog, 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); @@ -1882,8 +1882,9 @@ void do_force(FILE* fplog, 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)) { diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 5f5f863e81..f891b5bce7 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -1123,10 +1123,10 @@ void gmx::LegacySimulator::do_md() { 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); } } @@ -1468,9 +1468,9 @@ void gmx::LegacySimulator::do_md() 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, diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index daf0143822..438a426859 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -506,7 +506,7 @@ void gmx::LegacySimulator::do_mimic() { /* 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 */ diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index 64882158e2..c67ff09cc2 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -901,7 +901,10 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres, 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)) { diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index 2c59d9635d..2d12115cf5 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -588,13 +588,6 @@ void gmx::LegacySimulator::do_rerun() 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)) { diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index 01ddc14ceb..c33a68d3c2 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -1035,7 +1035,7 @@ void relax_shell_flexcon(FILE* fplog, 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]); @@ -1109,7 +1109,7 @@ void relax_shell_flexcon(FILE* fplog, 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); diff --git a/src/gromacs/modularsimulator/energyelement.cpp b/src/gromacs/modularsimulator/energyelement.cpp index 5e8abff78e..eb6de7e482 100644 --- a/src/gromacs/modularsimulator/energyelement.cpp +++ b/src/gromacs/modularsimulator/energyelement.cpp @@ -238,7 +238,8 @@ void EnergyElement::doStep(Time time, bool isEnergyCalculationStep, bool isFreeE } if (freeEnergyPerturbationElement_) { - sum_dhdl(enerd_, freeEnergyPerturbationElement_->constLambdaView(), *inputrec_->fepvals); + accumulateKineticLambdaComponents(enerd_, freeEnergyPerturbationElement_->constLambdaView(), + *inputrec_->fepvals); dummyLegacyState_.fep_state = freeEnergyPerturbationElement_->currentFEPState(); } if (parrinelloRahmanBarostat_) diff --git a/src/gromacs/nbnxm/kerneldispatch.cpp b/src/gromacs/nbnxm/kerneldispatch.cpp index 71b9d87a17..ea47bb2090 100644 --- a/src/gromacs/nbnxm/kerneldispatch.cpp +++ b/src/gromacs/nbnxm/kerneldispatch.cpp @@ -460,7 +460,7 @@ void nonbonded_verlet_t::dispatchFreeEnergyKernel(gmx::InteractionLocality iLo gmx::ForceWithShiftForces* forceWithShiftForces, const t_mdatoms& mdatoms, t_lambda* fepvals, - real* lambda, + gmx::ArrayRef lambda, gmx_enerdata_t* enerd, const gmx::StepWorkload& stepWork, t_nrnb* nrnb) @@ -493,7 +493,7 @@ void nonbonded_verlet_t::dispatchFreeEnergyKernel(gmx::InteractionLocality iLo 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(); @@ -557,7 +557,7 @@ void nonbonded_verlet_t::dispatchFreeEnergyKernel(gmx::InteractionLocality iLo 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]; } diff --git a/src/gromacs/nbnxm/nbnxm.h b/src/gromacs/nbnxm/nbnxm.h index f1d7bed568..bcb3b34c33 100644 --- a/src/gromacs/nbnxm/nbnxm.h +++ b/src/gromacs/nbnxm/nbnxm.h @@ -331,7 +331,7 @@ public: gmx::ForceWithShiftForces* forceWithShiftForces, const t_mdatoms& mdatoms, t_lambda* fepvals, - real* lambda, + gmx::ArrayRef lambda, gmx_enerdata_t* enerd, const gmx::StepWorkload& stepWork, t_nrnb* nrnb); -- 2.22.0