Add lambda state to dhdl.xvg with AWH fep
[alexxy/gromacs.git] / src / gromacs / mdrun / mimic.cpp
index a434de0126b78830281683731e4fc81f90615ab6..92d183f08328fb497629d4f8d73c2118a0432a33 100644 (file)
@@ -57,6 +57,7 @@
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/domdec/domdec_network.h"
 #include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/domdec/localtopologychecker.h"
 #include "gromacs/domdec/mdsetup.h"
 #include "gromacs/domdec/partition.h"
 #include "gromacs/essentialdynamics/edsam.h"
 #include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/mdtypes/mdrunoptions.h"
 #include "gromacs/mdtypes/observableshistory.h"
+#include "gromacs/mdtypes/observablesreducer.h"
 #include "gromacs/mdtypes/simulation_workload.h"
 #include "gromacs/mdtypes/state.h"
 #include "gromacs/mimic/communicator.h"
@@ -142,7 +144,6 @@ using gmx::SimulationSignaller;
 void gmx::LegacySimulator::do_mimic()
 {
     const t_inputrec* ir = inputrec;
-    int64_t           step, step_rel;
     double            t;
     bool              isLastStep               = false;
     bool              doFreeEnergyPerturbation = false;
@@ -206,16 +207,18 @@ void gmx::LegacySimulator::do_mimic()
     int        nstglobalcomm = 1;
     const bool bNS           = true;
 
+    ObservablesReducer observablesReducer = observablesReducerBuilder->build();
+
     if (MASTER(cr))
     {
         MimicCommunicator::init();
-        auto nonConstGlobalTopology = const_cast<gmx_mtop_t*>(&top_global);
+        auto* nonConstGlobalTopology = const_cast<gmx_mtop_t*>(&top_global);
         MimicCommunicator::sendInitData(nonConstGlobalTopology, state_global->x);
         // TODO: Avoid changing inputrec (#3854)
         auto* nonConstInputrec   = const_cast<t_inputrec*>(inputrec);
         nonConstInputrec->nsteps = MimicCommunicator::getStepNumber();
     }
-    if (DOMAINDECOMP(cr))
+    if (haveDDAtomOrdering(*cr))
     {
         // TODO: Avoid changing inputrec (#3854)
         auto* nonConstInputrec = const_cast<t_inputrec*>(inputrec);
@@ -224,12 +227,15 @@ void gmx::LegacySimulator::do_mimic()
 
     const SimulationGroups* groups = &top_global.groups;
     {
-        auto nonConstGlobalTopology                          = const_cast<gmx_mtop_t*>(&top_global);
+        auto* nonConstGlobalTopology                         = const_cast<gmx_mtop_t*>(&top_global);
         nonConstGlobalTopology->intermolecularExclusionGroup = genQmmmIndices(top_global);
     }
 
     initialize_lambdas(fplog,
-                       *ir,
+                       ir->efep,
+                       ir->bSimTemp,
+                       *ir->fepvals,
+                       ir->simtempvals->temperatures,
                        gmx::arrayRefFromArray(ir->opts.ref_t, ir->opts.ngtc),
                        MASTER(cr),
                        &state_global->fep_state,
@@ -242,7 +248,7 @@ void gmx::LegacySimulator::do_mimic()
                                    mdrunOptions,
                                    cr,
                                    outputProvider,
-                                   mdModulesNotifier,
+                                   mdModulesNotifiers,
                                    ir,
                                    top_global,
                                    oenv,
@@ -258,7 +264,7 @@ void gmx::LegacySimulator::do_mimic()
                                    true,
                                    StartingBehavior::NewSimulation,
                                    simulationsShareState,
-                                   mdModulesNotifier);
+                                   mdModulesNotifiers);
 
     gstat = global_stat_init(ir);
 
@@ -267,7 +273,7 @@ void gmx::LegacySimulator::do_mimic()
                                  top_global,
                                  constr ? constr->numFlexibleConstraints() : 0,
                                  ir->nstcalcenergy,
-                                 DOMAINDECOMP(cr),
+                                 haveDDAtomOrdering(*cr),
                                  runScheduleWork->simulationWork.useGpuPme);
 
     {
@@ -278,16 +284,9 @@ void gmx::LegacySimulator::do_mimic()
         }
     }
 
-    // Local state only becomes valid now.
-    std::unique_ptr<t_state> stateInstance;
-    t_state*                 state;
-
-    gmx_localtop_t top(top_global.ffparams);
-
-    if (DOMAINDECOMP(cr))
+    if (haveDDAtomOrdering(*cr))
     {
-        stateInstance = std::make_unique<t_state>();
-        state         = stateInstance.get();
+        // Local state only becomes valid now.
         dd_init_local_state(*cr->dd, state_global, state);
 
         /* Distribute the charge groups over the nodes from the master node */
@@ -305,7 +304,7 @@ void gmx::LegacySimulator::do_mimic()
                             state,
                             &f,
                             mdAtoms,
-                            &top,
+                            top,
                             fr,
                             vsite,
                             constr,
@@ -316,13 +315,10 @@ void gmx::LegacySimulator::do_mimic()
     else
     {
         state_change_natoms(state_global, state_global->natoms);
-        /* Copy the pointer to the global state */
-        state = state_global;
-
-        mdAlgorithmsSetupAtomData(cr, *ir, top_global, &top, fr, &f, mdAtoms, constr, vsite, shellfc);
+        mdAlgorithmsSetupAtomData(cr, *ir, top_global, top, fr, &f, mdAtoms, constr, vsite, shellfc);
     }
 
-    auto mdatoms = mdAtoms->mdatoms();
+    auto* mdatoms = mdAtoms->mdatoms();
 
     // NOTE: The global state is no longer used at this point.
     // But state_global is still used as temporary storage space for writing
@@ -330,18 +326,18 @@ void gmx::LegacySimulator::do_mimic()
     // (Global topology should persist.)
 
     update_mdatoms(mdatoms, state->lambda[FreeEnergyPerturbationCouplingType::Mass]);
+    fr->longRangeNonbondeds->updateAfterPartition(*mdatoms);
 
     if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl != 0)
     {
         doFreeEnergyPerturbation = true;
     }
 
+    int64_t step     = ir->init_step;
+    int64_t step_rel = 0;
+
     {
-        int cglo_flags = CGLO_GSTAT;
-        if (DOMAINDECOMP(cr) && shouldCheckNumberOfBondedInteractions(*cr->dd))
-        {
-            cglo_flags |= CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS;
-        }
+        int    cglo_flags   = CGLO_GSTAT;
         bool   bSumEkinhOld = false;
         t_vcm* vcm          = nullptr;
         compute_globals(gstat,
@@ -361,16 +357,14 @@ void gmx::LegacySimulator::do_mimic()
                         shake_vir,
                         total_vir,
                         pres,
-                        gmx::ArrayRef<real>{},
                         &nullSignaller,
                         state->box,
                         &bSumEkinhOld,
-                        cglo_flags);
-        if (DOMAINDECOMP(cr))
-        {
-            checkNumberOfBondedInteractions(
-                    mdlog, cr, top_global, &top, makeConstArrayRef(state->x), state->box);
-        }
+                        cglo_flags,
+                        step,
+                        &observablesReducer);
+        // Clean up after pre-step use of compute_globals()
+        observablesReducer.markAsReadyToReduce();
     }
 
     if (MASTER(cr))
@@ -387,7 +381,7 @@ void gmx::LegacySimulator::do_mimic()
     }
 
     walltime_accounting_start_time(walltime_accounting);
-    wallcycle_start(wcycle, ewcRUN);
+    wallcycle_start(wcycle, WallCycleCounter::Run);
     print_start(fplog, cr, walltime_accounting, "mdrun");
 
     /***********************************************************
@@ -411,9 +405,6 @@ void gmx::LegacySimulator::do_mimic()
                     "MiMiC does not report kinetic energy, total energy, temperature, virial and "
                     "pressure.");
 
-    step     = ir->init_step;
-    step_rel = 0;
-
     auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
             compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]),
             false,
@@ -438,13 +429,13 @@ void gmx::LegacySimulator::do_mimic()
     while (!isLastStep)
     {
         isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel == ir->nsteps));
-        wallcycle_start(wcycle, ewcSTEP);
+        wallcycle_start(wcycle, WallCycleCounter::Step);
 
         t = step;
 
         if (MASTER(cr))
         {
-            MimicCommunicator::getCoords(&state_global->x, state_global->natoms);
+            MimicCommunicator::getCoords(state_global->x, state_global->natoms);
         }
 
         if (ir->efep != FreeEnergyPerturbationType::No)
@@ -455,7 +446,7 @@ void gmx::LegacySimulator::do_mimic()
         if (MASTER(cr))
         {
             const bool constructVsites = ((vsite != nullptr) && mdrunOptions.rerunConstructVsites);
-            if (constructVsites && DOMAINDECOMP(cr))
+            if (constructVsites && haveDDAtomOrdering(*cr))
             {
                 gmx_fatal(FARGS,
                           "Vsite recalculation with -rerun is not implemented with domain "
@@ -464,13 +455,13 @@ void gmx::LegacySimulator::do_mimic()
             }
             if (constructVsites)
             {
-                wallcycle_start(wcycle, ewcVSITECONSTR);
+                wallcycle_start(wcycle, WallCycleCounter::VsiteConstr);
                 vsite->construct(state->x, state->v, state->box, VSiteOperation::PositionsAndVelocities);
-                wallcycle_stop(wcycle, ewcVSITECONSTR);
+                wallcycle_stop(wcycle, WallCycleCounter::VsiteConstr);
             }
         }
 
-        if (DOMAINDECOMP(cr))
+        if (haveDDAtomOrdering(*cr))
         {
             /* Repartition the domain decomposition */
             const bool bMasterState = true;
@@ -488,7 +479,7 @@ void gmx::LegacySimulator::do_mimic()
                                 state,
                                 &f,
                                 mdAtoms,
-                                &top,
+                                top,
                                 fr,
                                 vsite,
                                 constr,
@@ -507,6 +498,8 @@ void gmx::LegacySimulator::do_mimic()
             update_mdatoms(mdatoms, state->lambda[FreeEnergyPerturbationCouplingType::Mass]);
         }
 
+        fr->longRangeNonbondeds->updateAfterPartition(*mdatoms);
+
         force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES
                        | GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 is solved
                        GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0));
@@ -525,7 +518,7 @@ void gmx::LegacySimulator::do_mimic()
                                 pull_work,
                                 bNS,
                                 force_flags,
-                                &top,
+                                top,
                                 constr,
                                 enerd,
                                 state->natoms,
@@ -537,6 +530,7 @@ void gmx::LegacySimulator::do_mimic()
                                 &f.view(),
                                 force_vir,
                                 *mdatoms,
+                                fr->longRangeNonbondeds.get(),
                                 nrnb,
                                 wcycle,
                                 shellfc,
@@ -567,7 +561,7 @@ void gmx::LegacySimulator::do_mimic()
                      step,
                      nrnb,
                      wcycle,
-                     &top,
+                     top,
                      state->box,
                      state->x.arrayRefWithPadding(),
                      &state->hist,
@@ -582,6 +576,7 @@ void gmx::LegacySimulator::do_mimic()
                      mu_tot,
                      t,
                      ed,
+                     fr->longRangeNonbondeds.get(),
                      GMX_FORCE_NS | force_flags,
                      ddBalanceRegionHandler);
         }
@@ -627,10 +622,6 @@ void gmx::LegacySimulator::do_mimic()
             SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
 
             int cglo_flags = CGLO_GSTAT | CGLO_ENERGY;
-            if (DOMAINDECOMP(cr) && shouldCheckNumberOfBondedInteractions(*cr->dd))
-            {
-                cglo_flags |= CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS;
-            }
             compute_globals(gstat,
                             cr,
                             ir,
@@ -648,23 +639,19 @@ void gmx::LegacySimulator::do_mimic()
                             nullptr,
                             nullptr,
                             nullptr,
-                            constr != nullptr ? constr->rmsdData() : gmx::ArrayRef<real>{},
                             &signaller,
                             state->box,
                             &bSumEkinhOld,
-                            cglo_flags);
-            if (DOMAINDECOMP(cr))
-            {
-                checkNumberOfBondedInteractions(
-                        mdlog, cr, top_global, &top, makeConstArrayRef(state->x), state->box);
-            }
+                            cglo_flags,
+                            step,
+                            &observablesReducer);
         }
 
         {
             gmx::HostVector<gmx::RVec>     fglobal(top_global.natoms);
             gmx::ArrayRef<gmx::RVec>       ftemp;
             gmx::ArrayRef<const gmx::RVec> flocal = f.view().force();
-            if (DOMAINDECOMP(cr))
+            if (haveDDAtomOrdering(*cr))
             {
                 ftemp = gmx::makeArrayRef(fglobal);
                 dd_collect_vec(cr->dd, state->ddp_count, state->ddp_count_cg_gl, state->cg_gl, flocal, ftemp);
@@ -704,7 +691,6 @@ void gmx::LegacySimulator::do_mimic()
                                              mdatoms->tmass,
                                              enerd,
                                              ir->fepvals.get(),
-                                             ir->expandedvals.get(),
                                              state->box,
                                              PTCouplingArrays({ state->boxv,
                                                                 state->nosehoover_xi,
@@ -712,8 +698,6 @@ void gmx::LegacySimulator::do_mimic()
                                                                 state->nhpres_xi,
                                                                 state->nhpres_vxi }),
                                              state->fep_state,
-                                             shake_vir,
-                                             force_vir,
                                              total_vir,
                                              pres,
                                              ekind,
@@ -756,8 +740,8 @@ void gmx::LegacySimulator::do_mimic()
             print_time(stderr, walltime_accounting, step, ir, cr);
         }
 
-        cycles = wallcycle_stop(wcycle, ewcSTEP);
-        if (DOMAINDECOMP(cr) && wcycle)
+        cycles = wallcycle_stop(wcycle, WallCycleCounter::Step);
+        if (haveDDAtomOrdering(*cr) && wcycle)
         {
             dd_cycles_add(cr->dd, cycles, ddCyclStep);
         }
@@ -765,6 +749,7 @@ void gmx::LegacySimulator::do_mimic()
         /* increase the MD step number */
         step++;
         step_rel++;
+        observablesReducer.markAsReadyToReduce();
     }
     /* End of main MD loop */