Move dispersion correction call to do_force()
authorBerk Hess <hess@kth.se>
Fri, 26 Jun 2020 08:30:32 +0000 (08:30 +0000)
committerMagnus Lundborg <magnus.lundborg@scilifelab.se>
Fri, 26 Jun 2020 08:30:32 +0000 (08:30 +0000)
The call to dispersion correction was in compute_globals(), which
should only compute globals, not compute energies. Moving the call
to do_force() simplifies the logic, as correcting the virial before
computing the pressure removes the need to correct the pressure.

src/gromacs/mdlib/dispersioncorrection.h
src/gromacs/mdlib/enerdata_utils.cpp
src/gromacs/mdlib/md_support.cpp
src/gromacs/mdlib/md_support.h
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdrun/md.cpp
src/gromacs/mdrun/mimic.cpp
src/gromacs/mdrun/rerun.cpp
src/gromacs/modularsimulator/computeglobalselement.cpp

index 1e4b223f17fd920cfc62dcac0ad81a89e14c55f2..be656e2130c4769ff0963a28d5a3127acf8ce959 100644 (file)
@@ -109,18 +109,6 @@ public:
             }
         }
 
-        /*! \brief Correct the pressure tensor for the missing dispersion
-         *
-         * \param[in,out] pressureTensor  The pressure tensor to correct
-         */
-        void correctPressure(tensor pressureTensor) const
-        {
-            for (int m = 0; m < DIM; m++)
-            {
-                pressureTensor[m][m] += pressure;
-            }
-        }
-
         real virial   = 0; //!< Scalar correction to the virial
         real pressure = 0; //!< Scalar correction to the pressure
         real energy   = 0; //!< Correction to the energy
index 4ae6e8f22e4f86acb87386dfa270c75a460cf3c5..c5e1bd33d880f355b11a48ad15a9691e62506aa3 100644 (file)
@@ -232,21 +232,22 @@ void reset_enerdata(gmx_enerdata_t* enerd)
     {
         for (j = 0; (j < enerd->grpp.nener); j++)
         {
-            enerd->grpp.ener[i][j] = 0.0;
+            enerd->grpp.ener[i][j] = 0.0_real;
         }
     }
 
     /* Normal potential energy components */
     for (i = 0; (i <= F_EPOT); i++)
     {
-        enerd->term[i] = 0.0;
+        enerd->term[i] = 0.0_real;
     }
-    enerd->term[F_DVDL]           = 0.0;
-    enerd->term[F_DVDL_COUL]      = 0.0;
-    enerd->term[F_DVDL_VDW]       = 0.0;
-    enerd->term[F_DVDL_BONDED]    = 0.0;
-    enerd->term[F_DVDL_RESTRAINT] = 0.0;
-    enerd->term[F_DKDL]           = 0.0;
+    enerd->term[F_PDISPCORR]      = 0.0_real;
+    enerd->term[F_DVDL]           = 0.0_real;
+    enerd->term[F_DVDL_COUL]      = 0.0_real;
+    enerd->term[F_DVDL_VDW]       = 0.0_real;
+    enerd->term[F_DVDL_BONDED]    = 0.0_real;
+    enerd->term[F_DVDL_RESTRAINT] = 0.0_real;
+    enerd->term[F_DKDL]           = 0.0_real;
     std::fill(enerd->enerpart_lambda.begin(), enerd->enerpart_lambda.end(), 0);
     std::fill(enerd->dhdlLambda.begin(), enerd->dhdlLambda.end(), 0);
     /* reset foreign energy data and dvdl - separate functions since they are also called elsewhere */
index 8803c5a99b974abb058ae4f8b717028446d67da7..efad98cc8a98b5e53252c49e1ef77d5ffc10f748 100644 (file)
@@ -50,7 +50,6 @@
 #include "gromacs/gmxlib/nrnb.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/coupling.h"
-#include "gromacs/mdlib/dispersioncorrection.h"
 #include "gromacs/mdlib/gmx_omp_nthreads.h"
 #include "gromacs/mdlib/simulationsignal.h"
 #include "gromacs/mdlib/stat.h"
@@ -376,7 +375,6 @@ void compute_globals(gmx_global_stat*               gstat,
                      gmx::ArrayRef<const gmx::RVec> x,
                      gmx::ArrayRef<const gmx::RVec> v,
                      const matrix                   box,
-                     real                           vdwLambda,
                      const t_mdatoms*               mdatoms,
                      t_nrnb*                        nrnb,
                      t_vcm*                         vcm,
@@ -504,32 +502,6 @@ void compute_globals(gmx_global_stat*               gstat,
 
         enerd->term[F_PRES] = calc_pres(fr->pbcType, ir->nwall, lastbox, ekind->ekin, total_vir, pres);
     }
-
-    /* ##########  Long range energy information ###### */
-    if ((bEner || bPres) && fr->dispersionCorrection)
-    {
-        /* Calculate long range corrections to pressure and energy */
-        /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
-           and computes enerd->term[F_DISPCORR].  Also modifies the
-           total_vir and pres tensors */
-
-        const DispersionCorrection::Correction correction =
-                fr->dispersionCorrection->calculate(lastbox, vdwLambda);
-
-        if (bEner)
-        {
-            enerd->term[F_DISPCORR] = correction.energy;
-            enerd->term[F_EPOT] += correction.energy;
-            enerd->term[F_DVDL_VDW] += correction.dvdl;
-        }
-        if (bPres)
-        {
-            correction.correctVirial(total_vir);
-            correction.correctPressure(pres);
-            enerd->term[F_PDISPCORR] = correction.pressure;
-            enerd->term[F_PRES] += correction.pressure;
-        }
-    }
 }
 
 void setCurrentLambdasRerun(int64_t           step,
index 333cde82ec8c90f7c280afb0fe830ad874fadbc5..87368d5e6c916cfb660d4a368cd7e2bc713ed831 100644 (file)
@@ -139,7 +139,6 @@ void compute_globals(gmx_global_stat*               gstat,
                      gmx::ArrayRef<const gmx::RVec> x,
                      gmx::ArrayRef<const gmx::RVec> v,
                      const matrix                   box,
-                     real                           vdwLambda,
                      const t_mdatoms*               mdatoms,
                      t_nrnb*                        nrnb,
                      t_vcm*                         vcm,
index 725cdd16ce5c384147239102074d51904cd80a55..9ace847885b5c5bb15784b0c734a1e42f09d664b 100644 (file)
@@ -73,6 +73,7 @@
 #include "gromacs/mdlib/calcmu.h"
 #include "gromacs/mdlib/calcvir.h"
 #include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/dispersioncorrection.h"
 #include "gromacs/mdlib/enerdata_utils.h"
 #include "gromacs/mdlib/force.h"
 #include "gromacs/mdlib/force_flags.h"
@@ -1843,6 +1844,25 @@ void do_force(FILE*                               fplog,
                                         vir_force, *mdatoms, *fr, vsite, stepWork);
     }
 
+    // VdW dispersion correction, only computed on master rank to avoid double counting
+    if ((stepWork.computeEnergy || stepWork.computeVirial) && fr->dispersionCorrection && MASTER(cr))
+    {
+        // Calculate long range corrections to pressure and energy
+        const DispersionCorrection::Correction correction =
+                fr->dispersionCorrection->calculate(box, lambda[efptVDW]);
+
+        if (stepWork.computeEnergy)
+        {
+            enerd->term[F_DISPCORR] = correction.energy;
+            enerd->term[F_DVDL_VDW] += correction.dvdl;
+        }
+        if (stepWork.computeVirial)
+        {
+            correction.correctVirial(vir_force);
+            enerd->term[F_PDISPCORR] = correction.pressure;
+        }
+    }
+
     // TODO refactor this and unify with above GPU PME-PP / GPU update path call to the same function
     if (PAR(cr) && !thisRankHasDuty(cr, DUTY_PME) && !simulationWork.useGpuPmePpCommunication
         && !simulationWork.useGpuUpdate)
index 4ace36f8938c635c788c6d01a7a3ae46021638cf..9ee0697c84d6c171340871cf7e4c45752d64edf2 100644 (file)
@@ -578,9 +578,9 @@ void gmx::LegacySimulator::do_md()
             cglo_flags_iteration &= ~CGLO_TEMPERATURE;
         }
         compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                        makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms,
-                        nrnb, &vcm, nullptr, enerd, force_vir, shake_vir, total_vir, pres, constr,
-                        &nullSignaller, state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                        makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, nullptr,
+                        enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
+                        state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
                         cglo_flags_iteration
                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
                                                                          : 0));
@@ -607,9 +607,9 @@ void gmx::LegacySimulator::do_md()
            perhaps loses some logic?*/
 
         compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                        makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms,
-                        nrnb, &vcm, nullptr, enerd, force_vir, shake_vir, total_vir, pres, constr,
-                        &nullSignaller, state->box, nullptr, &bSumEkinhOld, cglo_flags & ~CGLO_PRESSURE);
+                        makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, nullptr,
+                        enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
+                        state->box, nullptr, &bSumEkinhOld, cglo_flags & ~CGLO_PRESSURE);
     }
 
     /* Calculate the initial half step temperature, and save the ekinh_old */
@@ -897,9 +897,9 @@ void gmx::LegacySimulator::do_md()
              * the full step kinetic energy and possibly for T-coupling.*/
             /* This may not be quite working correctly yet . . . . */
             compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                            makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms,
-                            nrnb, &vcm, wcycle, enerd, nullptr, nullptr, nullptr, nullptr, constr,
-                            &nullSignaller, state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                            makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
+                            enerd, nullptr, nullptr, nullptr, nullptr, constr, &nullSignaller,
+                            state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
             checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global,
                                             &top, makeConstArrayRef(state->x), state->box,
@@ -1034,17 +1034,16 @@ void gmx::LegacySimulator::do_md()
             if (bGStat || do_per_step(step - 1, nstglobalcomm))
             {
                 wallcycle_stop(wcycle, ewcUPDATE);
-                compute_globals(
-                        gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                        makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms,
-                        nrnb, &vcm, wcycle, enerd, force_vir, shake_vir, total_vir, pres, constr,
-                        &nullSignaller, state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
-                        (bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
-                                | (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
-                                | (bPres ? CGLO_CONSTRAINT : 0) | (bStopCM ? CGLO_STOPCM : 0)
-                                | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
-                                                                         : 0)
-                                | CGLO_SCALEEKIN);
+                compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+                                makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
+                                enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
+                                state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                                (bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
+                                        | (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
+                                        | (bPres ? CGLO_CONSTRAINT : 0) | (bStopCM ? CGLO_STOPCM : 0)
+                                        | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
+                                                                                 : 0)
+                                        | CGLO_SCALEEKIN);
                 /* explanation of above:
                    a) We compute Ekin at the full time step
                    if 1) we are using the AveVel Ekin, and it's not the
@@ -1097,10 +1096,9 @@ void gmx::LegacySimulator::do_md()
                      * the full step kinetic energy and possibly for T-coupling.*/
                     /* This may not be quite working correctly yet . . . . */
                     compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                                    makeConstArrayRef(state->v), state->box, state->lambda[efptVDW],
-                                    mdatoms, nrnb, &vcm, wcycle, enerd, nullptr, nullptr, nullptr,
-                                    nullptr, constr, &nullSignaller, state->box, nullptr,
-                                    &bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
+                                    makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
+                                    enerd, nullptr, nullptr, nullptr, nullptr, constr, &nullSignaller,
+                                    state->box, nullptr, &bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
                     wallcycle_start(wcycle, ewcUPDATE);
                 }
             }
@@ -1347,10 +1345,9 @@ void gmx::LegacySimulator::do_md()
             /* erase F_EKIN and F_TEMP here? */
             /* just compute the kinetic energy at the half step to perform a trotter step */
             compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                            makeConstArrayRef(state->v), state->box, state->lambda[efptVDW],
-                            mdatoms, nrnb, &vcm, wcycle, enerd, force_vir, shake_vir, total_vir,
-                            pres, constr, &nullSignaller, lastbox, nullptr, &bSumEkinhOld,
-                            (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE);
+                            makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle, enerd,
+                            force_vir, shake_vir, total_vir, pres, constr, &nullSignaller, lastbox,
+                            nullptr, &bSumEkinhOld, (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE);
             wallcycle_start(wcycle, ewcUPDATE);
             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
             /* now we know the scaling, we can compute the positions again */
@@ -1425,17 +1422,16 @@ void gmx::LegacySimulator::do_md()
                 bool                doIntraSimSignal = true;
                 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
 
-                compute_globals(
-                        gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                        makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms,
-                        nrnb, &vcm, wcycle, enerd, force_vir, shake_vir, total_vir, pres, constr,
-                        &signaller, lastbox, &totalNumberOfBondedInteractions, &bSumEkinhOld,
-                        (bGStat ? CGLO_GSTAT : 0) | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
-                                | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
-                                | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
-                                | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0) | CGLO_CONSTRAINT
-                                | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
-                                                                         : 0));
+                compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+                                makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm,
+                                wcycle, enerd, force_vir, shake_vir, total_vir, pres, constr,
+                                &signaller, lastbox, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                                (bGStat ? CGLO_GSTAT : 0) | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
+                                        | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
+                                        | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
+                                        | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0) | CGLO_CONSTRAINT
+                                        | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
+                                                                                 : 0));
                 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
                                                 top_global, &top, makeConstArrayRef(state->x),
                                                 state->box, &shouldCheckNumberOfBondedInteractions);
index ef5a04f3d52b26faa98a56812501bdba460c4ad3..daf014382202636d6f7017ff884389f0e9de4ee6 100644 (file)
@@ -295,9 +295,9 @@ void gmx::LegacySimulator::do_mimic()
         bool   bSumEkinhOld = false;
         t_vcm* vcm          = nullptr;
         compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                        makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms, nrnb,
-                        vcm, nullptr, enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
-                        state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags);
+                        makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, nullptr, enerd,
+                        force_vir, shake_vir, total_vir, pres, constr, &nullSignaller, state->box,
+                        &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags);
     }
     checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global, &top,
                                     makeConstArrayRef(state->x), state->box,
@@ -464,9 +464,9 @@ void gmx::LegacySimulator::do_mimic()
             SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
 
             compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                            makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms,
-                            nrnb, vcm, wcycle, enerd, nullptr, nullptr, nullptr, nullptr, constr,
-                            &signaller, state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                            makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle,
+                            enerd, nullptr, nullptr, nullptr, nullptr, constr, &signaller,
+                            state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
                             CGLO_GSTAT | CGLO_ENERGY
                                     | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
                                                                              : 0));
index 164c5082384cd6d192a94af153d8249f01c73194..2c59d9635dd0aabff1bd7bc851ecc2d737b194ad 100644 (file)
@@ -347,9 +347,9 @@ void gmx::LegacySimulator::do_rerun()
         bool   bSumEkinhOld = false;
         t_vcm* vcm          = nullptr;
         compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                        makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms, nrnb,
-                        vcm, nullptr, enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
-                        state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags);
+                        makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, nullptr, enerd,
+                        force_vir, shake_vir, total_vir, pres, constr, &nullSignaller, state->box,
+                        &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags);
     }
     checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global, &top,
                                     makeConstArrayRef(state->x), state->box,
@@ -572,9 +572,9 @@ void gmx::LegacySimulator::do_rerun()
             SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
 
             compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                            makeConstArrayRef(state->v), state->box, state->lambda[efptVDW], mdatoms,
-                            nrnb, vcm, wcycle, enerd, force_vir, shake_vir, total_vir, pres, constr,
-                            &signaller, state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                            makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle,
+                            enerd, force_vir, shake_vir, total_vir, pres, constr, &signaller,
+                            state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
                             CGLO_GSTAT | CGLO_ENERGY
                                     | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
                                                                              : 0));
index 6a3bcc55e745deb8c8011e00d25ccf2223b3d2a6..233e6fca83780692c44abc4d20e3305133055dc9 100644 (file)
@@ -280,13 +280,9 @@ void ComputeGlobalsElement<algorithm>::compute(gmx::Step            step,
     auto lastbox = useLastBox ? statePropagatorData_->constPreviousBox()
                               : statePropagatorData_->constBox();
 
-    const real vdwLambda = freeEnergyPerturbationElement_
-                                   ? freeEnergyPerturbationElement_->constLambdaView()[efptVDW]
-                                   : 0;
-
     compute_globals(
-            gstat_, cr_, inputrec_, fr_, energyElement_->ekindata(), x, v, box, vdwLambda,
-            mdAtoms_->mdatoms(), nrnb_, &vcm_, step != -1 ? wcycle_ : nullptr, energyElement_->enerdata(),
+            gstat_, cr_, inputrec_, fr_, energyElement_->ekindata(), x, v, box, mdAtoms_->mdatoms(),
+            nrnb_, &vcm_, step != -1 ? wcycle_ : nullptr, energyElement_->enerdata(),
             energyElement_->forceVirial(step), energyElement_->constraintVirial(step),
             energyElement_->totalVirial(step), energyElement_->pressure(step), constr_, signaller,
             lastbox, &totalNumberOfBondedInteractions_, energyElement_->needToSumEkinhOld(),