Add lambda state to dhdl.xvg with AWH fep

[alexxy/gromacs.git] / src / gromacs / mdrun / minimize.cpp

diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp
index 497b3accd9af174738c961cf74c5b17fa1b50881..eec144faf42604fa6ae6a41a0607abc4bf580253 100644 (file)
--- a/src/gromacs/mdrun/minimize.cpp
+++ b/src/gromacs/mdrun/minimize.cpp
@@ -315,7 +315,7 @@ static void get_f_norm_max(const t_commrec*               cr,
         }
     }
 
-    if (la_max >= 0 && DOMAINDECOMP(cr))
+    if (la_max >= 0 && haveDDAtomOrdering(*cr))
     {
         a_max = cr->dd->globalAtomIndices[la_max];
     }
@@ -414,7 +414,7 @@ static void init_em(FILE*                fplog,
                                       top_global,
                                       constr ? constr->numFlexibleConstraints() : 0,
                                       ir->nstcalcenergy,
-                                      DOMAINDECOMP(cr),
+                                      haveDDAtomOrdering(*cr),
                                       thisRankHasDuty(cr, DUTY_PME));
     }
     else
@@ -432,8 +432,9 @@ static void init_em(FILE*                fplog,
         }
     }
 
-    if (DOMAINDECOMP(cr))
+    if (haveDDAtomOrdering(*cr))
     {
+        // Local state only becomes valid now.
         dd_init_local_state(*cr->dd, state_global, &ems->s);
 
         /* Distribute the charge groups over the nodes from the master node */
@@ -595,7 +596,7 @@ static void write_em_traj(FILE*               fplog,
 
     if (confout != nullptr)
     {
-        if (DOMAINDECOMP(cr))
+        if (haveDDAtomOrdering(*cr))
         {
             /* If bX=true, x was collected to state_global in the call above */
             if (!bX)
@@ -613,7 +614,7 @@ static void write_em_traj(FILE*               fplog,
 
         if (MASTER(cr))
         {
-            if (ir->pbcType != PbcType::No && !ir->bPeriodicMols && DOMAINDECOMP(cr))
+            if (ir->pbcType != PbcType::No && !ir->bPeriodicMols && haveDDAtomOrdering(*cr))
             {
                 /* Make molecules whole only for confout writing */
                 do_pbc_mtop(ir->pbcType, state->s.box, &top_global, state_global->x.rvec_array());
@@ -654,7 +655,7 @@ static bool do_em_step(const t_commrec*                          cr,
     s1 = &ems1->s;
     s2 = &ems2->s;
 
-    if (DOMAINDECOMP(cr) && s1->ddp_count != cr->dd->ddp_count)
+    if (haveDDAtomOrdering(*cr) && s1->ddp_count != cr->dd->ddp_count)
     {
         gmx_incons("state mismatch in do_em_step");
     }
@@ -666,7 +667,7 @@ static bool do_em_step(const t_commrec*                          cr,
         state_change_natoms(s2, s1->natoms);
         ems2->f.resize(s2->natoms);
     }
-    if (DOMAINDECOMP(cr) && s2->cg_gl.size() != s1->cg_gl.size())
+    if (haveDDAtomOrdering(*cr) && s2->cg_gl.size() != s1->cg_gl.size())
     {
         s2->cg_gl.resize(s1->cg_gl.size());
     }
@@ -724,7 +725,7 @@ static bool do_em_step(const t_commrec*                          cr,
             }
         }
 
-        if (DOMAINDECOMP(cr))
+        if (haveDDAtomOrdering(*cr))
         {
             /* OpenMP does not supported unsigned loop variables */
 #pragma omp for schedule(static) nowait
@@ -735,7 +736,7 @@ static bool do_em_step(const t_commrec*                          cr,
         }
     }
 
-    if (DOMAINDECOMP(cr))
+    if (haveDDAtomOrdering(*cr))
     {
         s2->ddp_count       = s1->ddp_count;
         s2->ddp_count_cg_gl = s1->ddp_count_cg_gl;
@@ -966,7 +967,7 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres,
     // Compute the buffer size of the pair list
     const real bufferSize = inputrec->rlist - std::max(inputrec->rcoulomb, inputrec->rvdw);
 
-    if (bFirst || bufferSize <= 0 || (DOMAINDECOMP(cr) && ems->s.ddp_count != ddpCountPairSearch))
+    if (bFirst || bufferSize <= 0 || (haveDDAtomOrdering(*cr) && ems->s.ddp_count != ddpCountPairSearch))
     {
         /* This is the first state or an old state used before the last ns */
         bNS = TRUE;
@@ -980,7 +981,7 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres,
               > bufferSize;
     }
 
-    if (DOMAINDECOMP(cr) && bNS)
+    if (haveDDAtomOrdering(*cr) && bNS)
     {
         /* Repartition the domain decomposition */
         em_dd_partition_system(
@@ -1049,7 +1050,6 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres,
                     shake_vir,
                     *inputrec,
                     nullptr,
-                    gmx::ArrayRef<real>{},
                     nullptr,
                     std::vector<real>(1, terminate),
                     FALSE,
@@ -1201,7 +1201,7 @@ static real pr_beta(const t_commrec*  cr,
      * and might have to sum it in parallel runs.
      */
 
-    if (!DOMAINDECOMP(cr)
+    if (!haveDDAtomOrdering(*cr)
         || (s_min->s.ddp_count == cr->dd->ddp_count && s_b->s.ddp_count == cr->dd->ddp_count))
     {
         auto fm = s_min->f.view().force();
@@ -1246,7 +1246,6 @@ void LegacySimulator::do_cg()
 {
     const char* CG = "Polak-Ribiere Conjugate Gradients";
 
-    gmx_localtop_t    top(top_global.ffparams);
     gmx_global_stat_t gstat;
     double            tmp, minstep;
     real              stepsize;
@@ -1306,7 +1305,7 @@ void LegacySimulator::do_cg()
             state_global,
             top_global,
             s_min,
-            &top,
+            top,
             nrnb,
             fr,
             mdAtoms,
@@ -1359,7 +1358,7 @@ void LegacySimulator::do_cg()
                                      cr,
                                      ms,
                                      top_global,
-                                     &top,
+                                     top,
                                      inputrec,
                                      imdSession,
                                      pull_work,
@@ -1380,6 +1379,7 @@ void LegacySimulator::do_cg()
      * We do not unshift, so molecules are always whole in congrad.c
      */
     energyEvaluator.run(s_min, mu_tot, vir, pres, -1, TRUE, step);
+    observablesReducer.markAsReadyToReduce();
 
     if (MASTER(cr))
     {
@@ -1391,7 +1391,6 @@ void LegacySimulator::do_cg()
                                          mdatoms->tmass,
                                          enerd,
                                          nullptr,
-                                         nullptr,
                                          nullBox,
                                          PTCouplingArrays(),
                                          0,
@@ -1547,7 +1546,7 @@ void LegacySimulator::do_cg()
         a         = 0.0;
         c         = a + stepsize; /* reference position along line is zero */
 
-        if (DOMAINDECOMP(cr) && s_min->s.ddp_count < cr->dd->ddp_count)
+        if (haveDDAtomOrdering(*cr) && s_min->s.ddp_count < cr->dd->ddp_count)
         {
             em_dd_partition_system(fplog,
                                    mdlog,
@@ -1558,7 +1557,7 @@ void LegacySimulator::do_cg()
                                    imdSession,
                                    pull_work,
                                    s_min,
-                                   &top,
+                                   top,
                                    mdAtoms,
                                    fr,
                                    vsite,
@@ -1573,6 +1572,7 @@ void LegacySimulator::do_cg()
         neval++;
         /* Calculate energy for the trial step */
         energyEvaluator.run(s_c, mu_tot, vir, pres, -1, FALSE, step);
+        observablesReducer.markAsReadyToReduce();
 
         /* Calc derivative along line */
         const rvec*                    pc  = s_c->s.cg_p.rvec_array();
@@ -1661,7 +1661,7 @@ void LegacySimulator::do_cg()
                     b = 0.5 * (a + c);
                 }
 
-                if (DOMAINDECOMP(cr) && s_min->s.ddp_count != cr->dd->ddp_count)
+                if (haveDDAtomOrdering(*cr) && s_min->s.ddp_count != cr->dd->ddp_count)
                 {
                     /* Reload the old state */
                     em_dd_partition_system(fplog,
@@ -1673,7 +1673,7 @@ void LegacySimulator::do_cg()
                                            imdSession,
                                            pull_work,
                                            s_min,
-                                           &top,
+                                           top,
                                            mdAtoms,
                                            fr,
                                            vsite,
@@ -1688,6 +1688,7 @@ void LegacySimulator::do_cg()
                 neval++;
                 /* Calculate energy for the trial step */
                 energyEvaluator.run(s_b, mu_tot, vir, pres, -1, FALSE, step);
+                observablesReducer.markAsReadyToReduce();
 
                 /* p does not change within a step, but since the domain decomposition
                  * might change, we have to use cg_p of s_b here.
@@ -1840,7 +1841,6 @@ void LegacySimulator::do_cg()
                                              mdatoms->tmass,
                                              enerd,
                                              nullptr,
-                                             nullptr,
                                              nullBox,
                                              PTCouplingArrays(),
                                              0,
@@ -1880,7 +1880,7 @@ void LegacySimulator::do_cg()
          * If we have reached machine precision, converged is already set to true.
          */
         converged = converged || (s_min->fmax < inputrec->em_tol);
-
+        observablesReducer.markAsReadyToReduce();
     } /* End of the loop */
 
     if (converged)
@@ -1964,7 +1964,6 @@ void LegacySimulator::do_lbfgs()
 {
     static const char* LBFGS = "Low-Memory BFGS Minimizer";
     em_state_t         ems;
-    gmx_localtop_t     top(top_global.ffparams);
     gmx_global_stat_t  gstat;
     auto*              mdatoms = mdAtoms->mdatoms();
 
@@ -1976,6 +1975,10 @@ void LegacySimulator::do_lbfgs()
                     "be available in a different form in a future version of GROMACS, "
                     "e.g. gmx minimize and an .mdp option.");
 
+    if (haveDDAtomOrdering(*cr))
+    {
+        gmx_fatal(FARGS, "L_BFGS is currently not supported");
+    }
     if (PAR(cr))
     {
         gmx_fatal(FARGS, "L-BFGS minimization only supports a single rank");
@@ -2024,7 +2027,7 @@ void LegacySimulator::do_lbfgs()
             state_global,
             top_global,
             &ems,
-            &top,
+            top,
             nrnb,
             fr,
             mdAtoms,
@@ -2115,7 +2118,7 @@ void LegacySimulator::do_lbfgs()
                                      cr,
                                      ms,
                                      top_global,
-                                     &top,
+                                     top,
                                      inputrec,
                                      imdSession,
                                      pull_work,
@@ -2144,7 +2147,6 @@ void LegacySimulator::do_lbfgs()
                                          mdatoms->tmass,
                                          enerd,
                                          nullptr,
-                                         nullptr,
                                          nullBox,
                                          PTCouplingArrays(),
                                          0,
@@ -2666,7 +2668,6 @@ void LegacySimulator::do_lbfgs()
                                              mdatoms->tmass,
                                              enerd,
                                              nullptr,
-                                             nullptr,
                                              nullBox,
                                              PTCouplingArrays(),
                                              0,
@@ -2709,7 +2710,7 @@ void LegacySimulator::do_lbfgs()
          * If we have reached machine precision, converged is already set to true.
          */
         converged = converged || (ems.fmax < inputrec->em_tol);
-
+        observablesReducer.markAsReadyToReduce();
     } /* End of the loop */
 
     if (converged)
@@ -2781,7 +2782,6 @@ void LegacySimulator::do_lbfgs()
 void LegacySimulator::do_steep()
 {
     const char*       SD = "Steepest Descents";
-    gmx_localtop_t    top(top_global.ffparams);
     gmx_global_stat_t gstat;
     real              stepsize;
     real              ustep;
@@ -2819,7 +2819,7 @@ void LegacySimulator::do_steep()
             state_global,
             top_global,
             s_try,
-            &top,
+            top,
             nrnb,
             fr,
             mdAtoms,
@@ -2878,7 +2878,7 @@ void LegacySimulator::do_steep()
                                      cr,
                                      ms,
                                      top_global,
-                                     &top,
+                                     top,
                                      inputrec,
                                      imdSession,
                                      pull_work,
@@ -2960,7 +2960,6 @@ void LegacySimulator::do_steep()
                                                  mdatoms->tmass,
                                                  enerd,
                                                  nullptr,
-                                                 nullptr,
                                                  nullBox,
                                                  PTCouplingArrays(),
                                                  0,
@@ -3012,7 +3011,7 @@ void LegacySimulator::do_steep()
             /* If energy is not smaller make the step smaller...  */
             ustep *= 0.5;
 
-            if (DOMAINDECOMP(cr) && s_min->s.ddp_count != cr->dd->ddp_count)
+            if (haveDDAtomOrdering(*cr) && s_min->s.ddp_count != cr->dd->ddp_count)
             {
                 /* Reload the old state */
                 em_dd_partition_system(fplog,
@@ -3024,7 +3023,7 @@ void LegacySimulator::do_steep()
                                        imdSession,
                                        pull_work,
                                        s_min,
-                                       &top,
+                                       top,
                                        mdAtoms,
                                        fr,
                                        vsite,
@@ -3070,6 +3069,7 @@ void LegacySimulator::do_steep()
         }
 
         count++;
+        observablesReducer.markAsReadyToReduce();
     } /* End of the loop  */
 
     /* Print some data...  */
@@ -3114,7 +3114,6 @@ void LegacySimulator::do_nm()
 {
     const char*         NM = "Normal Mode Analysis";
     int                 nnodes;
-    gmx_localtop_t      top(top_global.ffparams);
     gmx_global_stat_t   gstat;
     tensor              vir, pres;
     rvec                mu_tot = { 0 };
@@ -3164,7 +3163,7 @@ void LegacySimulator::do_nm()
             state_global,
             top_global,
             &state_work,
-            &top,
+            top,
             nrnb,
             fr,
             mdAtoms,
@@ -3270,7 +3269,7 @@ void LegacySimulator::do_nm()
                                      cr,
                                      ms,
                                      top_global,
-                                     &top,
+                                     top,
                                      inputrec,
                                      imdSession,
                                      pull_work,
@@ -3352,7 +3351,7 @@ void LegacySimulator::do_nm()
                                         pull_work,
                                         bNS,
                                         force_flags,
-                                        &top,
+                                        top,
                                         constr,
                                         enerd,
                                         state_work.s.natoms,