Store dHdL for all neighbor lambdas

This is required for calculating forces along a lambda dimension
in AWH.

Change-Id: Ic5247a53998b1e733858419d20338b5419dd7ae4

diff --git a/src/gromacs/listed_forces/listed_forces.cpp b/src/gromacs/listed_forces/listed_forces.cpp
index ceaa6cdc3f0090f41f5a5765505109cc2ee3cd44..2a850054af4e1f1c2cb142af03e3967fbbd09b29 100644 (file)
--- 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<real>   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 d9c1df340bd29075851e9d2c4e711da335b3256a..5ad7d35dc00093d80aa3827daffd6951ee6b1a54 100644 (file)
--- 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<typename>
+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<real>   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 36482d071239191eca072af68550b25b83a073c9..30449fe0e29a39282ceb4bf62d4cb314edcc604c 100644 (file)
--- 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<false>(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 43afbccb43f69afb0ce0ca3197a05bc7a6eb15ce..4ae6e8f22e4f86acb87386dfa270c75a460cf3c5 100644 (file)
--- 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<const real> 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 e1e9d3231ffc28c98429c108f414c3f676b04726..bfefd04bfee9bab5fe2ba312896048a7b1266417 100644 (file)
--- 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 e5ecf8c0c3bcd56139dad36a960c8be01475fb33..acaf09c158c531b1b1b3c60c4d49465eba5edf1a 100644 (file)
--- 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 fb3a048e72466a8034d0eeed3aeb2ece39bf24ab..9eb763a25d1a47b826353e1950c422999c3cd49e 100644 (file)
--- 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<double>(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 585ad0a80b24431d08da501fd4d652fc317d8a35..df856edbd33d639b322110534746e437b5812e16 100644 (file)
--- 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 5180f4e3534a85d68ad51d9d23b67799ec364aba..9866d7740b9252c19cad779945323ac8e2770283 100644 (file)
--- 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<double> enerpart_lambda; /* Partial Hamiltonian for lambda and flambda[], includes at least all perturbed terms */
+    std::vector<double> 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 6cd3991b6b40df90ad551973d801e4b8c3812ac4..71b9d87a17831c3193d7c64c0437eb25ec38ad9a 100644 (file)
--- 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);