From 439bff03a472589f9708b93818a57b2602dbb318 Mon Sep 17 00:00:00 2001 From: Magnus Lundborg Date: Sun, 23 Feb 2020 17:57:56 +0100 Subject: [PATCH] Store dHdL for all neighbor lambdas This is required for calculating forces along a lambda dimension in AWH. Change-Id: Ic5247a53998b1e733858419d20338b5419dd7ae4 --- src/gromacs/listed_forces/listed_forces.cpp | 12 ++++++++--- src/gromacs/listed_forces/listed_forces.h | 3 +++ src/gromacs/listed_forces/position_restraints.cpp | 10 +++++---- src/gromacs/mdlib/enerdata_utils.cpp | 25 +++++++++++++++++------ src/gromacs/mdlib/enerdata_utils.h | 3 +++ src/gromacs/mdlib/force.cpp | 10 +++++++++ src/gromacs/mdlib/md_support.cpp | 5 +++++ src/gromacs/mdlib/sim_util.cpp | 9 ++++++++ src/gromacs/mdtypes/enerdata.h | 3 ++- src/gromacs/nbnxm/kerneldispatch.cpp | 11 +++++++++- 10 files changed, 76 insertions(+), 15 deletions(-) diff --git a/src/gromacs/listed_forces/listed_forces.cpp b/src/gromacs/listed_forces/listed_forces.cpp index ceaa6cdc3f..2a850054af 100644 --- a/src/gromacs/listed_forces/listed_forces.cpp +++ b/src/gromacs/listed_forces/listed_forces.cpp @@ -586,6 +586,7 @@ void calc_listed_lambda(const t_idef* idef, const struct t_graph* g, gmx_grppairener_t* grpp, real* epot, + gmx::ArrayRef dvdl, t_nrnb* nrnb, const real* lambda, const t_mdatoms* md, @@ -593,7 +594,6 @@ void calc_listed_lambda(const t_idef* idef, int* global_atom_index) { real v; - real dvdl_dum[efptNR] = { 0 }; rvec4* f; rvec* fshift; const t_pbc* pbc_null; @@ -637,7 +637,7 @@ void calc_listed_lambda(const t_idef* idef, gmx::StepWorkload tempFlags; tempFlags.computeEnergy = true; v = calc_one_bond(0, ftype, idef, iatoms, iatoms.ssize(), workDivision, x, f, - fshift, fr, pbc_null, g, grpp, nrnb, lambda, dvdl_dum, md, fcd, + fshift, fr, pbc_null, g, grpp, nrnb, lambda, dvdl.data(), md, fcd, tempFlags, global_atom_index); epot[ftype] += v; } @@ -688,6 +688,7 @@ void do_force_listed(struct gmx_wallcycle* wcycle, */ if (fepvals->n_lambda > 0 && stepWork.computeDhdl) { + real dvdl[efptNR] = { 0 }; posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr); if (idef->ilsort != ilsortNO_FE) @@ -707,9 +708,14 @@ void do_force_listed(struct gmx_wallcycle* wcycle, lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i - 1]); } calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), - enerd->foreign_term, nrnb, lam_i, md, fcd, global_atom_index); + enerd->foreign_term, dvdl, nrnb, lam_i, md, fcd, global_atom_index); sum_epot(&(enerd->foreign_grpp), enerd->foreign_term); enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT]; + for (int j = 0; j < efptNR; j++) + { + enerd->dhdlLambda[i] += dvdl[j]; + dvdl[j] = 0; + } } wallcycle_sub_stop(wcycle, ewcsLISTED_FEP); } diff --git a/src/gromacs/listed_forces/listed_forces.h b/src/gromacs/listed_forces/listed_forces.h index d9c1df340b..5ad7d35dc0 100644 --- a/src/gromacs/listed_forces/listed_forces.h +++ b/src/gromacs/listed_forces/listed_forces.h @@ -89,6 +89,8 @@ namespace gmx { class ForceOutputs; class StepWorkload; +template +class ArrayRef; } // namespace gmx //! Type of CPU function to compute a bonded interaction. @@ -143,6 +145,7 @@ void calc_listed_lambda(const t_idef* idef, const struct t_graph* g, gmx_grppairener_t* grpp, real* epot, + gmx::ArrayRef dvdl, t_nrnb* nrnb, const real* lambda, const t_mdatoms* md, diff --git a/src/gromacs/listed_forces/position_restraints.cpp b/src/gromacs/listed_forces/position_restraints.cpp index 36482d0712..30449fe0e2 100644 --- a/src/gromacs/listed_forces/position_restraints.cpp +++ b/src/gromacs/listed_forces/position_restraints.cpp @@ -443,13 +443,15 @@ void posres_wrapper_lambda(struct gmx_wallcycle* wcycle, wallcycle_sub_start_nocount(wcycle, ewcsRESTRAINTS); for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++) { - real dvdl_dum = 0, lambda_dum; + real dvdl = 0; - lambda_dum = (i == 0 ? lambda[efptRESTRAINT] : fepvals->all_lambda[efptRESTRAINT][i - 1]); + const real lambda_dum = + (i == 0 ? lambda[efptRESTRAINT] : fepvals->all_lambda[efptRESTRAINT][i - 1]); v = posres(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x, - nullptr, fr->pbcType == PbcType::No ? nullptr : pbc, lambda_dum, - &dvdl_dum, fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB); + nullptr, fr->pbcType == PbcType::No ? nullptr : pbc, lambda_dum, &dvdl, + fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB); enerd->enerpart_lambda[i] += v; + enerd->dhdlLambda[i] += dvdl; } wallcycle_sub_stop(wcycle, ewcsRESTRAINTS); } diff --git a/src/gromacs/mdlib/enerdata_utils.cpp b/src/gromacs/mdlib/enerdata_utils.cpp index 43afbccb43..4ae6e8f22e 100644 --- a/src/gromacs/mdlib/enerdata_utils.cpp +++ b/src/gromacs/mdlib/enerdata_utils.cpp @@ -47,6 +47,7 @@ gmx_enerdata_t::gmx_enerdata_t(int numEnergyGroups, int numFepLambdas) : grpp(numEnergyGroups), enerpart_lambda(numFepLambdas == 0 ? 0 : numFepLambdas + 1), + dhdlLambda(numFepLambdas == 0 ? 0 : numFepLambdas + 1), foreign_grpp(numEnergyGroups) { } @@ -95,6 +96,12 @@ void sum_dhdl(gmx_enerdata_t* enerd, gmx::ArrayRef lambda, const t_l int index; enerd->dvdl_lin[efptVDW] += enerd->term[F_DVDL_VDW]; /* include dispersion correction */ + + for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++) + { + enerd->dhdlLambda[i] += enerd->term[F_DVDL_VDW]; + } + enerd->term[F_DVDL] = 0.0; for (int i = 0; i < efptNR; i++) { @@ -207,6 +214,15 @@ void reset_foreign_enerdata(gmx_enerdata_t* enerd) } } +void reset_dvdl_enerdata(gmx_enerdata_t* enerd) +{ + for (int i = 0; i < efptNR; i++) + { + enerd->dvdl_lin[i] = 0.0; + enerd->dvdl_nonlin[i] = 0.0; + } +} + void reset_enerdata(gmx_enerdata_t* enerd) { int i, j; @@ -219,11 +235,6 @@ void reset_enerdata(gmx_enerdata_t* enerd) enerd->grpp.ener[i][j] = 0.0; } } - for (i = 0; i < efptNR; i++) - { - enerd->dvdl_lin[i] = 0.0; - enerd->dvdl_nonlin[i] = 0.0; - } /* Normal potential energy components */ for (i = 0; (i <= F_EPOT); i++) @@ -237,6 +248,8 @@ void reset_enerdata(gmx_enerdata_t* enerd) enerd->term[F_DVDL_RESTRAINT] = 0.0; enerd->term[F_DKDL] = 0.0; std::fill(enerd->enerpart_lambda.begin(), enerd->enerpart_lambda.end(), 0); - /* reset foreign energy data - separate function since we also call it elsewhere */ + std::fill(enerd->dhdlLambda.begin(), enerd->dhdlLambda.end(), 0); + /* reset foreign energy data and dvdl - separate functions since they are also called elsewhere */ reset_foreign_enerdata(enerd); + reset_dvdl_enerdata(enerd); } diff --git a/src/gromacs/mdlib/enerdata_utils.h b/src/gromacs/mdlib/enerdata_utils.h index e1e9d3231f..bfefd04bfe 100644 --- a/src/gromacs/mdlib/enerdata_utils.h +++ b/src/gromacs/mdlib/enerdata_utils.h @@ -50,6 +50,9 @@ struct t_lambda; void reset_foreign_enerdata(gmx_enerdata_t* enerd); /* Resets only the foreign energy data */ +void reset_dvdl_enerdata(gmx_enerdata_t* enerd); +/* Resets only the dvdl energy data */ + void reset_enerdata(gmx_enerdata_t* enerd); /* Resets the energy data */ diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp index e5ecf8c0c3..acaf09c158 100644 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@ -139,6 +139,11 @@ void do_force_lowlevel(t_forcerec* fr, real dvdl_walls = do_walls(*ir, *fr, box, *md, x, &forceWithVirial, lambda[efptVDW], enerd->grpp.ener[egLJSR].data(), nrnb); enerd->dvdl_lin[efptVDW] += dvdl_walls; + + for (auto& dhdl : enerd->dhdlLambda) + { + dhdl += dvdl_walls; + } } /* Shift the coordinates. Must be done before listed forces and PPPM, @@ -350,6 +355,11 @@ void do_force_lowlevel(t_forcerec* fr, enerd->term[F_COUL_RECIP] = Vlr_q + ewaldOutput.Vcorr_q; enerd->term[F_LJ_RECIP] = Vlr_lj + ewaldOutput.Vcorr_lj; + for (auto& dhdl : enerd->dhdlLambda) + { + dhdl += ewaldOutput.dvdl[efptVDW] + ewaldOutput.dvdl[efptCOUL]; + } + if (debug) { fprintf(debug, "Vlr_q = %g, Vcorr_q = %g, Vlr_corr_q = %g\n", Vlr_q, diff --git a/src/gromacs/mdlib/md_support.cpp b/src/gromacs/mdlib/md_support.cpp index fb3a048e72..9eb763a25d 100644 --- a/src/gromacs/mdlib/md_support.cpp +++ b/src/gromacs/mdlib/md_support.cpp @@ -271,6 +271,11 @@ void compute_globals(gmx_global_stat* gstat, enerd->dvdl_lin[efptMASS] = static_cast(dvdl_ekin); enerd->term[F_EKIN] = trace(ekind->ekin); + + for (auto& dhdl : enerd->dhdlLambda) + { + dhdl += enerd->dvdl_lin[efptMASS]; + } } /* ########## Now pressure ############## */ diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 585ad0a80b..df856edbd3 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -198,6 +198,10 @@ static void pull_potential_wrapper(const t_commrec* cr, enerd->term[F_COM_PULL] += pull_potential(pull_work, mdatoms, &pbc, cr, t, lambda[efptRESTRAINT], as_rvec_array(x.data()), force, &dvdl); enerd->dvdl_lin[efptRESTRAINT] += dvdl; + for (auto& dhdl : enerd->dhdlLambda) + { + dhdl += dvdl; + } wallcycle_stop(wcycle, ewcPULLPOT); } @@ -228,6 +232,11 @@ static void pme_receive_force_ener(t_forcerec* fr, enerd->dvdl_lin[efptCOUL] += dvdl_q; enerd->dvdl_lin[efptVDW] += dvdl_lj; + for (auto& dhdl : enerd->dhdlLambda) + { + dhdl += dvdl_q + dvdl_lj; + } + if (wcycle) { dd_cycles_add(cr->dd, cycles_seppme, ddCyclPME); diff --git a/src/gromacs/mdtypes/enerdata.h b/src/gromacs/mdtypes/enerdata.h index 5180f4e353..9866d7740b 100644 --- a/src/gromacs/mdtypes/enerdata.h +++ b/src/gromacs/mdtypes/enerdata.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,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. @@ -81,6 +81,7 @@ struct gmx_enerdata_t int fep_state = 0; /*current fep state -- just for printing */ std::vector enerpart_lambda; /* Partial Hamiltonian for lambda and flambda[], includes at least all perturbed terms */ + std::vector dhdlLambda; /* dHdL at all neighboring lambda points (the current lambda point also at index 0). */ real foreign_term[F_NRE] = { 0 }; /* alternate array for storing foreign lambda energies */ struct gmx_grppairener_t foreign_grpp; /* alternate array for storing foreign lambda energies */ }; diff --git a/src/gromacs/nbnxm/kerneldispatch.cpp b/src/gromacs/nbnxm/kerneldispatch.cpp index 6cd3991b6b..71b9d87a17 100644 --- a/src/gromacs/nbnxm/kerneldispatch.cpp +++ b/src/gromacs/nbnxm/kerneldispatch.cpp @@ -534,12 +534,13 @@ void nonbonded_verlet_t::dispatchFreeEnergyKernel(gmx::InteractionLocality iLo kernel_data.flags = (donb_flags & ~(GMX_NONBONDED_DO_FORCE | GMX_NONBONDED_DO_SHIFTFORCE)) | GMX_NONBONDED_DO_FOREIGNLAMBDA; kernel_data.lambda = lam_i; + kernel_data.dvdl = dvdl_nb; kernel_data.energygrp_elec = enerd->foreign_grpp.ener[egCOULSR].data(); kernel_data.energygrp_vdw = enerd->foreign_grpp.ener[egLJSR].data(); - /* Note that we add to kernel_data.dvdl, but ignore the result */ for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++) { + std::fill(std::begin(dvdl_nb), std::end(dvdl_nb), 0); for (int j = 0; j < efptNR; j++) { lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i - 1]); @@ -558,6 +559,14 @@ void nonbonded_verlet_t::dispatchFreeEnergyKernel(gmx::InteractionLocality iLo 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]; + } + } + else + { + for (size_t i = 0; i < enerd->enerpart_lambda.size(); i++) + { + enerd->dhdlLambda[i] += dvdl_nb[efptVDW] + dvdl_nb[efptCOUL]; } } wallcycle_sub_stop(wcycle_, ewcsNONBONDED_FEP); -- 2.11.4.GIT