Merge branch release-2016
[alexxy/gromacs.git] / src / programs / mdrun / md.cpp
index 4c4d069ef4f0712d010b4ca666c4c7bd618e180a..7fc77b592cd43f0fac67f7ae2a2e898b91d1b7c0 100644 (file)
@@ -45,6 +45,7 @@
 #include <stdlib.h>
 
 #include <algorithm>
+#include <memory>
 
 #include "thread_mpi/threads.h"
 
@@ -55,7 +56,6 @@
 #include "gromacs/ewald/pme.h"
 #include "gromacs/ewald/pme-load-balancing.h"
 #include "gromacs/fileio/trxio.h"
-#include "gromacs/gmxlib/md_logging.h"
 #include "gromacs/gmxlib/network.h"
 #include "gromacs/gmxlib/nrnb.h"
 #include "gromacs/gpu_utils/gpu_utils.h"
@@ -76,6 +76,7 @@
 #include "gromacs/mdlib/mdebin.h"
 #include "gromacs/mdlib/mdoutf.h"
 #include "gromacs/mdlib/mdrun.h"
+#include "gromacs/mdlib/mdsetup.h"
 #include "gromacs/mdlib/nb_verlet.h"
 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
 #include "gromacs/mdlib/ns.h"
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/logger.h"
 #include "gromacs/utility/real.h"
 #include "gromacs/utility/smalloc.h"
 
@@ -153,7 +155,7 @@ static void checkNumberOfBondedInteractions(FILE *fplog, t_commrec *cr, int tota
     }
 }
 
-static void reset_all_counters(FILE *fplog, t_commrec *cr,
+static void reset_all_counters(FILE *fplog, const gmx::MDLogger &mdlog, t_commrec *cr,
                                gmx_int64_t step,
                                gmx_int64_t *step_rel, t_inputrec *ir,
                                gmx_wallcycle_t wcycle, t_nrnb *nrnb,
@@ -163,8 +165,9 @@ static void reset_all_counters(FILE *fplog, t_commrec *cr,
     char sbuf[STEPSTRSIZE];
 
     /* Reset all the counters related to performance over the run */
-    md_print_warn(cr, fplog, "step %s: resetting all time and cycle counters\n",
-                  gmx_step_str(step, sbuf));
+    GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
+            "step %s: resetting all time and cycle counters",
+            gmx_step_str(step, sbuf));
 
     if (use_GPU(nbv))
     {
@@ -188,7 +191,7 @@ static void reset_all_counters(FILE *fplog, t_commrec *cr,
 }
 
 /*! \libinternal
-    \copydoc integrator_t (FILE *fplog, t_commrec *cr,
+    \copydoc integrator_t (FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
                            int nfile, const t_filenm fnm[],
                            const gmx_output_env_t *oenv, gmx_bool bVerbose,
                            int nstglobalcomm,
@@ -207,7 +210,8 @@ static void reset_all_counters(FILE *fplog, t_commrec *cr,
                            unsigned long Flags,
                            gmx_walltime_accounting_t walltime_accounting)
  */
-double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
+double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
+                  int nfile, const t_filenm fnm[],
                   const gmx_output_env_t *oenv, gmx_bool bVerbose,
                   int nstglobalcomm,
                   gmx_vsite_t *vsite, gmx_constr_t constr,
@@ -215,6 +219,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                   gmx_mtop_t *top_global,
                   t_fcdata *fcd,
                   t_state *state_global,
+                  energyhistory_t *energyHistory,
                   t_mdatoms *mdatoms,
                   t_nrnb *nrnb, gmx_wallcycle_t wcycle,
                   gmx_edsam_t ed, t_forcerec *fr,
@@ -225,7 +230,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                   unsigned long Flags,
                   gmx_walltime_accounting_t walltime_accounting)
 {
-    gmx_mdoutf_t    outf = NULL;
+    gmx_mdoutf_t    outf = nullptr;
     gmx_int64_t     step, step_rel;
     double          elapsed_time;
     double          t, t0, lam0[efptNR];
@@ -243,18 +248,17 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     t_trxstatus      *status;
     rvec              mu_tot;
     t_vcm            *vcm;
-    matrix            pcoupl_mu, M;
+    matrix            parrinellorahmanMu, M;
     t_trxframe        rerun_fr;
-    gmx_repl_ex_t     repl_ex = NULL;
+    gmx_repl_ex_t     repl_ex = nullptr;
     int               nchkpt  = 1;
     gmx_localtop_t   *top;
-    t_mdebin         *mdebin   = NULL;
-    t_state          *state    = NULL;
+    t_mdebin         *mdebin   = nullptr;
     gmx_enerdata_t   *enerd;
-    rvec             *f = NULL;
+    PaddedRVecVector  f {};
     gmx_global_stat_t gstat;
-    gmx_update_t     *upd   = NULL;
-    t_graph          *graph = NULL;
+    gmx_update_t     *upd   = nullptr;
+    t_graph          *graph = nullptr;
     gmx_groups_t     *groups;
     gmx_ekindata_t   *ekind;
     gmx_shellfc_t    *shellfc;
@@ -262,7 +266,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     gmx_bool          bResetCountersHalfMaxH = FALSE;
     gmx_bool          bTemp, bPres, bTrotter;
     real              dvdl_constr;
-    rvec             *cbuf        = NULL;
+    rvec             *cbuf        = nullptr;
     int               cbuf_nalloc = 0;
     matrix            lastbox;
     int               lamnew  = 0;
@@ -278,7 +282,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
 
     /* PME load balancing data for GPU kernels */
-    pme_load_balancing_t *pme_loadbal      = NULL;
+    pme_load_balancing_t *pme_loadbal      = nullptr;
     gmx_bool              bPMETune         = FALSE;
     gmx_bool              bPMETunePrinting = FALSE;
 
@@ -331,7 +335,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         nstglobalcomm     = 1;
     }
 
-    nstglobalcomm   = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
+    nstglobalcomm   = check_nstglobalcomm(mdlog, nstglobalcomm, ir);
     bGStatEveryStep = (nstglobalcomm == 1);
 
     if (bRerunMD)
@@ -344,7 +348,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     {
         /* Initialize ion swapping code */
         init_swapcoords(fplog, bVerbose, ir, opt2fn_master("-swap", nfile, fnm, cr),
-                        top_global, state_global->x, state_global->box, &state_global->swapstate, cr, oenv, Flags);
+                        top_global, as_rvec_array(state_global->x.data()), state_global->box, &state_global->swapstate, cr, oenv, Flags);
     }
 
     /* Initial values */
@@ -360,13 +364,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     snew(enerd, 1);
     init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
                   enerd);
-    if (DOMAINDECOMP(cr))
-    {
-        f = NULL;
-    }
-    else
+    if (!DOMAINDECOMP(cr))
     {
-        snew(f, top_global->natoms);
+        f.resize(top_global->natoms + 1);
     }
 
     /* Kinetic energy data */
@@ -410,43 +410,31 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
     }
 
+    std::unique_ptr<t_state> stateInstance;
+    t_state *                state;
+
     if (DOMAINDECOMP(cr))
     {
         top = dd_init_local_top(top_global);
 
-        snew(state, 1);
+        stateInstance = std::unique_ptr<t_state>(new t_state);
+        state         = stateInstance.get();
         dd_init_local_state(cr->dd, state_global, state);
     }
     else
     {
-        top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
-
-        state    = serial_init_local_state(state_global);
+        /* Copy the pointer to the global state */
+        state = state_global;
 
-        atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
+        snew(top, 1);
+        mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
+                                  &graph, mdatoms, vsite, shellfc);
 
-        if (vsite)
-        {
-            set_vsite_top(vsite, top, mdatoms, cr);
-        }
-
-        if (ir->ePBC != epbcNONE && !fr->bMolPBC)
-        {
-            graph = mk_graph(fplog, &(top->idef), 0, top_global->natoms, FALSE, FALSE);
-        }
-
-        if (shellfc)
-        {
-            make_local_shells(cr, mdatoms, shellfc);
-        }
-
-        setup_bonded_threading(fr, &top->idef);
-
-        update_realloc(upd, state->nalloc);
+        update_realloc(upd, state->natoms);
     }
 
     /* Set up interactive MD (IMD) */
-    init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, state_global->x,
+    init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, as_rvec_array(state_global->x.data()),
              nfile, fnm, oenv, imdport, Flags);
 
     if (DOMAINDECOMP(cr))
@@ -456,9 +444,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                             state_global, top_global, ir,
                             state, &f, mdatoms, top, fr,
                             vsite, constr,
-                            nrnb, NULL, FALSE);
+                            nrnb, nullptr, FALSE);
         shouldCheckNumberOfBondedInteractions = true;
-        update_realloc(upd, state->nalloc);
+        update_realloc(upd, state->natoms);
     }
 
     update_mdatoms(mdatoms, state->lambda[efptMASS]);
@@ -467,7 +455,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
     if (ir->bExpanded)
     {
-        init_expanded_ensemble(startingFromCheckpoint, ir, &state->dfhist);
+        init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
     }
 
     if (MASTER(cr))
@@ -477,19 +465,18 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             /* Update mdebin with energy history if appending to output files */
             if (Flags & MD_APPENDFILES)
             {
-                restore_energyhistory_from_state(mdebin, state_global->enerhist);
+                restore_energyhistory_from_state(mdebin, energyHistory);
             }
             else
             {
                 /* We might have read an energy history from checkpoint,
                  * free the allocated memory and reset the counts.
                  */
-                done_energyhistory(state_global->enerhist);
-                init_energyhistory(state_global->enerhist);
+                *energyHistory = {};
             }
         }
         /* Set the initial energy history in state by updating once */
-        update_energyhistory(state_global->enerhist, mdebin);
+        update_energyhistory(energyHistory, mdebin);
     }
 
     /* Initialize constraints */
@@ -511,23 +498,31 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 !(Flags & MD_REPRODUCIBLE));
     if (bPMETune)
     {
-        pme_loadbal_init(&pme_loadbal, cr, fplog, ir, state->box,
+        pme_loadbal_init(&pme_loadbal, cr, mdlog, ir, state->box,
                          fr->ic, fr->pmedata, use_GPU(fr->nbv),
                          &bPMETunePrinting);
     }
 
     if (!ir->bContinuation && !bRerunMD)
     {
-        if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
+        if (state->flags & (1 << estV))
         {
-            /* Set the velocities of frozen particles to zero */
+            /* Set the velocities of vsites, shells and frozen atoms to zero */
             for (i = 0; i < mdatoms->homenr; i++)
             {
-                for (m = 0; m < DIM; m++)
+                if (mdatoms->ptype[i] == eptVSite ||
+                    mdatoms->ptype[i] == eptShell)
                 {
-                    if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+                    clear_rvec(state->v[i]);
+                }
+                else if (mdatoms->cFREEZE)
+                {
+                    for (m = 0; m < DIM; m++)
                     {
-                        state->v[i][m] = 0;
+                        if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+                        {
+                            state->v[i][m] = 0;
+                        }
                     }
                 }
             }
@@ -542,7 +537,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         if (vsite)
         {
             /* Construct the virtual sites for the initial configuration */
-            construct_vsites(vsite, state->x, ir->delta_t, NULL,
+            construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, nullptr,
                              top->idef.iparams, top->idef.il,
                              fr->ePBC, fr->bMolPBC, cr, state->box);
         }
@@ -582,7 +577,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
     bSumEkinhOld = FALSE;
     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
-                    NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+                    nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
                     constr, &nullSignaller, state->box,
                     &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags
                     | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
@@ -598,9 +593,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
            perhaps loses some logic?*/
 
         compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
-                        NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+                        nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
                         constr, &nullSignaller, state->box,
-                        NULL, &bSumEkinhOld,
+                        nullptr, &bSumEkinhOld,
                         cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
     }
 
@@ -776,15 +771,15 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     {
         if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
         {
-            md_print_info(cr, fplog,
-                          "Note: The number of steps is not consistent across multi simulations,\n"
-                          "but we are proceeding anyway!\n");
+            GMX_LOG(mdlog.warning).appendText(
+                    "Note: The number of steps is not consistent across multi simulations,\n"
+                    "but we are proceeding anyway!");
         }
         if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
         {
-            md_print_info(cr, fplog,
-                          "Note: The initial step is not consistent across multi simulations,\n"
-                          "but we are proceeding anyway!\n");
+            GMX_LOG(mdlog.warning).appendText(
+                    "Note: The initial step is not consistent across multi simulations,\n"
+                    "but we are proceeding anyway!");
         }
     }
 
@@ -800,8 +795,8 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         {
             /* PME grid + cut-off optimization with GPUs or PME nodes */
             pme_loadbal_do(pme_loadbal, cr,
-                           (bVerbose && MASTER(cr)) ? stderr : NULL,
-                           fplog,
+                           (bVerbose && MASTER(cr)) ? stderr : nullptr,
+                           fplog, mdlog,
                            ir, fr, state,
                            wcycle,
                            step, step_rel,
@@ -898,15 +893,15 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                     /* Following is necessary because the graph may get out of sync
                      * with the coordinates if we only have every N'th coordinate set
                      */
-                    mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
-                    shift_self(graph, state->box, state->x);
+                    mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
+                    shift_self(graph, state->box, as_rvec_array(state->x.data()));
                 }
-                construct_vsites(vsite, state->x, ir->delta_t, state->v,
+                construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
                                  top->idef.iparams, top->idef.il,
                                  fr->ePBC, fr->bMolPBC, cr, state->box);
                 if (graph)
                 {
-                    unshift_self(graph, state->box, state->x);
+                    unshift_self(graph, state->box, as_rvec_array(state->x.data()));
                 }
             }
         }
@@ -983,7 +978,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                                     nrnb, wcycle,
                                     do_verbose && !bPMETunePrinting);
                 shouldCheckNumberOfBondedInteractions = true;
-                update_realloc(upd, state->nalloc);
+                update_realloc(upd, state->natoms);
             }
         }
 
@@ -1004,7 +999,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
              * the full step kinetic energy and possibly for T-coupling.*/
             /* This may not be quite working correctly yet . . . . */
             compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
-                            wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+                            wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
                             constr, &nullSignaller, state->box,
                             &totalNumberOfBondedInteractions, &bSumEkinhOld,
                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
@@ -1077,10 +1072,10 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             relax_shell_flexcon(fplog, cr, bVerbose, step,
                                 ir, bNS, force_flags, top,
                                 constr, enerd, fcd,
-                                state, f, force_vir, mdatoms,
+                                state, &f, force_vir, mdatoms,
                                 nrnb, wcycle, graph, groups,
                                 shellfc, fr, bBornRadii, t, mu_tot,
-                                vsite, mdoutf_get_fp_field(outf));
+                                vsite);
         }
         else
         {
@@ -1090,17 +1085,17 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
              * Check comments in sim_util.c
              */
             do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
-                     state->box, state->x, &state->hist,
-                     f, force_vir, mdatoms, enerd, fcd,
+                     state->box, &state->x, &state->hist,
+                     &f, force_vir, mdatoms, enerd, fcd,
                      state->lambda, graph,
-                     fr, vsite, mu_tot, t, mdoutf_get_fp_field(outf), ed, bBornRadii,
+                     fr, vsite, mu_tot, t, ed, bBornRadii,
                      (bNS ? GMX_FORCE_NS : 0) | force_flags);
         }
 
         if (EI_VV(ir->eI) && !startingFromCheckpoint && !bRerunMD)
         /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
         {
-            rvec *vbuf = NULL;
+            rvec *vbuf = nullptr;
 
             wallcycle_start(wcycle, ewcUPDATE);
             if (ir->eI == eiVV && bInitStep)
@@ -1112,7 +1107,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                  * so that the input is actually the initial step.
                  */
                 snew(vbuf, state->natoms);
-                copy_rvecn(state->v, vbuf, 0, state->natoms); /* should make this better for parallelizing? */
+                copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
             }
             else
             {
@@ -1120,15 +1115,15 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
             }
 
-            update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+            update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
                           ekind, M, upd, etrtVELOCITY1,
                           cr, constr);
 
             if (!bRerunMD || rerun_fr.bV || bForceUpdate)         /* Why is rerun_fr.bV here?  Unclear. */
             {
                 wallcycle_stop(wcycle, ewcUPDATE);
-                update_constraints(fplog, step, NULL, ir, mdatoms,
-                                   state, fr->bMolPBC, graph, f,
+                update_constraints(fplog, step, nullptr, ir, mdatoms,
+                                   state, fr->bMolPBC, graph, &f,
                                    &top->idef, shake_vir,
                                    cr, nrnb, wcycle, upd, constr,
                                    TRUE, bCalcVir);
@@ -1138,7 +1133,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             {
                 /* Need to unshift here if a do_force has been
                    called in the previous step */
-                unshift_self(graph, state->box, state->x);
+                unshift_self(graph, state->box, as_rvec_array(state->x.data()));
             }
             /* if VV, compute the pressure and constraints */
             /* For VV2, we strictly only need this if using pressure
@@ -1198,7 +1193,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                     if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
                     {
                         /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
-                        enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE);
+                        enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
                         enerd->term[F_EKIN] = trace(ekind->ekin);
                     }
                 }
@@ -1209,9 +1204,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                      * the full step kinetic energy and possibly for T-coupling.*/
                     /* This may not be quite working correctly yet . . . . */
                     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
-                                    wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+                                    wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
                                     constr, &nullSignaller, state->box,
-                                    NULL, &bSumEkinhOld,
+                                    nullptr, &bSumEkinhOld,
                                     CGLO_GSTAT | CGLO_TEMPERATURE);
                     wallcycle_start(wcycle, ewcUPDATE);
                 }
@@ -1219,7 +1214,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             /* if it's the initial step, we performed this first step just to get the constraint virial */
             if (ir->eI == eiVV && bInitStep)
             {
-                copy_rvecn(vbuf, state->v, 0, state->natoms);
+                copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
                 sfree(vbuf);
             }
             wallcycle_stop(wcycle, ewcUPDATE);
@@ -1228,7 +1223,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         /* compute the conserved quantity */
         if (EI_VV(ir->eI))
         {
-            saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ);
+            saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
             if (ir->eI == eiVV)
             {
                 last_ekin = enerd->term[F_EKIN];
@@ -1253,9 +1248,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                statistics, but if performing simulated tempering, we
                do update the velocities and the tau_t. */
 
-            lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, &state->dfhist, step, state->v, mdatoms);
+            lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
             /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
-            copy_df_history(&state_global->dfhist, &state->dfhist);
+            copy_df_history(state_global->dfhist, state->dfhist);
         }
 
         /* Now we have the energies and forces corresponding to the
@@ -1263,13 +1258,14 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
          * the update.
          */
         do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
-                                 ir, state, state_global, top_global, fr,
-                                 outf, mdebin, ekind, f,
+                                 ir, state, state_global, energyHistory,
+                                 top_global, fr,
+                                 outf, mdebin, ekind, &f,
                                  &nchkpt,
                                  bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
                                  bSumEkinhOld);
         /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
-        bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
+        bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
 
         /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
         if (startingFromCheckpoint && bTrotter)
@@ -1360,8 +1356,8 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
             if (constr && bIfRandomize)
             {
-                update_constraints(fplog, step, NULL, ir, mdatoms,
-                                   state, fr->bMolPBC, graph, f,
+                update_constraints(fplog, step, nullptr, ir, mdatoms,
+                                   state, fr->bMolPBC, graph, &f,
                                    &top->idef, tmp_vir,
                                    cr, nrnb, wcycle, upd, constr,
                                    TRUE, bCalcVir);
@@ -1392,13 +1388,15 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             else
             {
                 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
-                update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep);
+                update_pcouple_before_coordinates(fplog, step, ir, state,
+                                                  parrinellorahmanMu, M,
+                                                  bInitStep);
             }
 
             if (EI_VV(ir->eI))
             {
                 /* velocity half-step update */
-                update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+                update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
                               ekind, M, upd, etrtVELOCITY2,
                               cr, constr);
             }
@@ -1416,15 +1414,15 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                     cbuf_nalloc = state->natoms;
                     srenew(cbuf, cbuf_nalloc);
                 }
-                copy_rvecn(state->x, cbuf, 0, state->natoms);
+                copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
             }
 
-            update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+            update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
                           ekind, M, upd, etrtPOSITION, cr, constr);
             wallcycle_stop(wcycle, ewcUPDATE);
 
             update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
-                               fr->bMolPBC, graph, f,
+                               fr->bMolPBC, graph, &f,
                                &top->idef, shake_vir,
                                cr, nrnb, wcycle, upd, constr,
                                FALSE, bCalcVir);
@@ -1436,27 +1434,27 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
                                 constr, &nullSignaller, lastbox,
-                                NULL, &bSumEkinhOld,
+                                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 again */
-                copy_rvecn(cbuf, state->x, 0, state->natoms);
+                copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
 
-                update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+                update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
                               ekind, M, upd, etrtPOSITION, cr, constr);
                 wallcycle_stop(wcycle, ewcUPDATE);
 
-                /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
+                /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
                 /* are the small terms in the shake_vir here due
                  * to numerical errors, or are they important
                  * physically? I'm thinking they are just errors, but not completely sure.
                  * For now, will call without actually constraining, constr=NULL*/
-                update_constraints(fplog, step, NULL, ir, mdatoms,
-                                   state, fr->bMolPBC, graph, f,
+                update_constraints(fplog, step, nullptr, ir, mdatoms,
+                                   state, fr->bMolPBC, graph, &f,
                                    &top->idef, tmp_vir,
-                                   cr, nrnb, wcycle, upd, NULL,
+                                   cr, nrnb, wcycle, upd, nullptr,
                                    FALSE, bCalcVir);
             }
             if (EI_VV(ir->eI))
@@ -1485,23 +1483,23 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         else if (graph)
         {
             /* Need to unshift here */
-            unshift_self(graph, state->box, state->x);
+            unshift_self(graph, state->box, as_rvec_array(state->x.data()));
         }
 
-        if (vsite != NULL)
+        if (vsite != nullptr)
         {
             wallcycle_start(wcycle, ewcVSITECONSTR);
-            if (graph != NULL)
+            if (graph != nullptr)
             {
-                shift_self(graph, state->box, state->x);
+                shift_self(graph, state->box, as_rvec_array(state->x.data()));
             }
-            construct_vsites(vsite, state->x, ir->delta_t, state->v,
+            construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
                              top->idef.iparams, top->idef.il,
                              fr->ePBC, fr->bMolPBC, cr, state->box);
 
-            if (graph != NULL)
+            if (graph != nullptr)
             {
-                unshift_self(graph, state->box, state->x);
+                unshift_self(graph, state->box, as_rvec_array(state->x.data()));
             }
             wallcycle_stop(wcycle, ewcVSITECONSTR);
         }
@@ -1560,8 +1558,10 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                Currently done every step so that dhdl is correct in the .edr */
             sum_dhdl(enerd, state->lambda, ir->fepvals);
         }
-        update_box(fplog, step, ir, mdatoms, state, f,
-                   pcoupl_mu, nrnb, upd);
+
+        update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
+                                         pres, parrinellorahmanMu,
+                                         state, nrnb, upd);
 
         /* ################# END UPDATE STEP 2 ################# */
         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
@@ -1590,7 +1590,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
         else
         {
-            enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir, state, &MassQ);
+            enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
         }
         /* #########  END PREPARING EDR OUTPUT  ###########  */
 
@@ -1600,8 +1600,8 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             if (fplog && do_log && bDoExpanded)
             {
                 /* only needed if doing expanded ensemble */
-                PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
-                                          &state_global->dfhist, state->fep_state, ir->nstlog, step);
+                PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
+                                          state_global->dfhist, state->fep_state, ir->nstlog, step);
             }
             if (bCalcEner)
             {
@@ -1619,7 +1619,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             gmx_bool do_dr  = do_per_step(step, ir->nstdisreout);
             gmx_bool do_or  = do_per_step(step, ir->nstorireout);
 
-            print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : NULL,
+            print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
                        step, t,
                        eprNORMAL, mdebin, fcd, groups, &(ir->opts));
 
@@ -1662,9 +1662,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             do_per_step(step, ir->swap->nstswap))
         {
             bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
-                                             bRerunMD ? rerun_fr.x   : state->x,
+                                             bRerunMD ? rerun_fr.x   : as_rvec_array(state->x.data()),
                                              bRerunMD ? rerun_fr.box : state->box,
-                                             top_global, MASTER(cr) && bVerbose, bRerunMD);
+                                             MASTER(cr) && bVerbose, bRerunMD);
 
             if (bNeedRepartition && DOMAINDECOMP(cr))
             {
@@ -1689,7 +1689,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                                 vsite, constr,
                                 nrnb, wcycle, FALSE);
             shouldCheckNumberOfBondedInteractions = true;
-            update_realloc(upd, state->nalloc);
+            update_realloc(upd, state->natoms);
         }
 
         bFirstStep             = FALSE;
@@ -1712,9 +1712,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
 
-        if ( (membed != NULL) && (!bLastStep) )
+        if ( (membed != nullptr) && (!bLastStep) )
         {
-            rescale_membed(step_rel, membed, state_global->x);
+            rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
         }
 
         if (bRerunMD)
@@ -1766,8 +1766,8 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                           "resetting counters later in the run, e.g. with gmx "
                           "mdrun -resetstep.", step);
             }
-            reset_all_counters(fplog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
-                               use_GPU(fr->nbv) ? fr->nbv : NULL);
+            reset_all_counters(fplog, mdlog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
+                               use_GPU(fr->nbv) ? fr->nbv : nullptr);
             wcycle_set_reset_counters(wcycle, -1);
             if (!(cr->duty & DUTY_PME))
             {
@@ -1815,11 +1815,11 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
     }
 
-    done_mdoutf(outf);
+    done_mdoutf(outf, ir);
 
     if (bPMETune)
     {
-        pme_loadbal_done(pme_loadbal, cr, fplog, use_GPU(fr->nbv));
+        pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
     }
 
     done_shellfc(fplog, shellfc, step_rel);