Merge release-5-0 into master
[alexxy/gromacs.git] / src / gromacs / mdlib / minimize.c
index 10e4ae43fc38261baeb3068805bc2b008dbc2e42..5255d544d146f0fa3cea26efe922ef5f8665a302 100644 (file)
@@ -1,84 +1,84 @@
-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
  *
- *
- *                This source code is part of
- *
- *                 G   R   O   M   A   C   S
- *
- *          GROningen MAchine for Chemical Simulations
- *
- *                        VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
  * of the License, or (at your option) any later version.
  *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
  *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  *
- * For more info, check our website at http://www.gromacs.org
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
  *
- * And Hey:
- * GROwing Monsters And Cloning Shrimps
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "config.h"
 
+#include <math.h>
 #include <string.h>
 #include <time.h>
-#include <math.h>
-#include "sysstuff.h"
-#include "string2.h"
-#include "network.h"
-#include "confio.h"
-#include "smalloc.h"
-#include "nrnb.h"
-#include "main.h"
-#include "force.h"
-#include "macros.h"
-#include "random.h"
-#include "names.h"
-#include "gmx_fatal.h"
-#include "txtdump.h"
-#include "typedefs.h"
-#include "update.h"
-#include "constr.h"
-#include "vec.h"
-#include "statutil.h"
-#include "tgroup.h"
-#include "mdebin.h"
-#include "vsite.h"
-#include "force.h"
-#include "mdrun.h"
-#include "md_support.h"
-#include "domdec.h"
-#include "partdec.h"
-#include "trnio.h"
-#include "mdatoms.h"
-#include "ns.h"
-#include "gmx_wallcycle.h"
-#include "mtop_util.h"
-#include "gmxfio.h"
-#include "pme.h"
-#include "bondf.h"
-#include "gmx_omp_nthreads.h"
-
-
-#include "gromacs/linearalgebra/mtxio.h"
+
+#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/mtxio.h"
+#include "gromacs/fileio/trajectory_writing.h"
+#include "gromacs/imd/imd.h"
+#include "gromacs/legacyheaders/bonded-threading.h"
+#include "gromacs/legacyheaders/constr.h"
+#include "gromacs/legacyheaders/domdec.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/md_logging.h"
+#include "gromacs/legacyheaders/md_support.h"
+#include "gromacs/legacyheaders/mdatoms.h"
+#include "gromacs/legacyheaders/mdebin.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/ns.h"
+#include "gromacs/legacyheaders/pme.h"
+#include "gromacs/legacyheaders/sim_util.h"
+#include "gromacs/legacyheaders/tgroup.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/update.h"
+#include "gromacs/legacyheaders/vsite.h"
+#include "gromacs/legacyheaders/types/commrec.h"
 #include "gromacs/linearalgebra/sparsematrix.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
 
 typedef struct {
     t_state  s;
@@ -101,25 +101,22 @@ static em_state_t *init_em_state()
     return ems;
 }
 
-static void print_em_start(FILE *fplog, t_commrec *cr, gmx_runtime_t *runtime,
-                           gmx_wallcycle_t wcycle,
-                           const char *name)
+static void print_em_start(FILE                     *fplog,
+                           t_commrec                *cr,
+                           gmx_walltime_accounting_t walltime_accounting,
+                           gmx_wallcycle_t           wcycle,
+                           const char               *name)
 {
-    char buf[STRLEN];
-
-    runtime_start(runtime);
-
-    sprintf(buf, "Started %s", name);
-    print_date_and_time(fplog, cr->nodeid, buf, NULL);
-
+    walltime_accounting_start(walltime_accounting);
     wallcycle_start(wcycle, ewcRUN);
+    print_start(fplog, cr, walltime_accounting, name);
 }
-static void em_time_end(FILE *fplog, t_commrec *cr, gmx_runtime_t *runtime,
-                        gmx_wallcycle_t wcycle)
+static void em_time_end(gmx_walltime_accounting_t walltime_accounting,
+                        gmx_wallcycle_t           wcycle)
 {
     wallcycle_stop(wcycle, ewcRUN);
 
-    runtime_end(runtime);
+    walltime_accounting_end(walltime_accounting);
 }
 
 static void sp_header(FILE *out, const char *minimizer, real ftol, int nsteps)
@@ -136,25 +133,25 @@ static void warn_step(FILE *fp, real ftol, gmx_bool bLastStep, gmx_bool bConstra
     if (bLastStep)
     {
         sprintf(buffer,
-                "\nEnergy minimization reached the maximum number"
-                "of steps before the forces reached the requested"
+                "\nEnergy minimization reached the maximum number "
+                "of steps before the forces reached the requested "
                 "precision Fmax < %g.\n", ftol);
     }
     else
     {
         sprintf(buffer,
-                "\nEnergy minimization has stopped, but the forces have"
-                "not converged to the requested precision Fmax < %g (which"
-                "may not be possible for your system). It stopped"
-                "because the algorithm tried to make a new step whose size"
-                "was too small, or there was no change in the energy since"
-                "last step. Either way, we regard the minimization as"
-                "converged to within the available machine precision,"
+                "\nEnergy minimization has stopped, but the forces have "
+                "not converged to the requested precision Fmax < %g (which "
+                "may not be possible for your system). It stopped "
+                "because the algorithm tried to make a new step whose size "
+                "was too small, or there was no change in the energy since "
+                "last step. Either way, we regard the minimization as "
+                "converged to within the available machine precision, "
                 "given your starting configuration and EM parameters.\n%s%s",
                 ftol,
                 sizeof(real) < sizeof(double) ?
-                "\nDouble precision normally gives you higher accuracy, but"
-                "this is often not needed for preparing to run molecular"
+                "\nDouble precision normally gives you higher accuracy, but "
+                "this is often not needed for preparing to run molecular "
                 "dynamics.\n" :
                 "",
                 bConstrain ?
@@ -168,7 +165,7 @@ static void warn_step(FILE *fp, real ftol, gmx_bool bLastStep, gmx_bool bConstra
 
 
 static void print_converged(FILE *fp, const char *alg, real ftol,
-                            gmx_large_int_t count, gmx_bool bDone, gmx_large_int_t nsteps,
+                            gmx_int64_t count, gmx_bool bDone, gmx_int64_t nsteps,
                             real epot, real fmax, int nfmax, real fnorm)
 {
     char buf[STEPSTRSIZE];
@@ -216,8 +213,8 @@ static void get_f_norm_max(t_commrec *cr,
     fmax2  = 0;
     la_max = -1;
     gf     = 0;
-    start  = mdatoms->start;
-    end    = mdatoms->homenr + start;
+    start  = 0;
+    end    = mdatoms->homenr;
     if (mdatoms->cFREEZE)
     {
         for (i = start; i < end; i++)
@@ -312,9 +309,11 @@ void init_em(FILE *fplog, const char *title,
              t_graph **graph, t_mdatoms *mdatoms, gmx_global_stat_t *gstat,
              gmx_vsite_t *vsite, gmx_constr_t constr,
              int nfile, const t_filenm fnm[],
-             gmx_mdoutf_t **outf, t_mdebin **mdebin)
+             gmx_mdoutf_t *outf, t_mdebin **mdebin,
+             int imdport, unsigned long gmx_unused Flags,
+             gmx_wallcycle_t wcycle)
 {
-    int  start, homenr, i;
+    int  i;
     real dvdl_constr;
 
     if (fplog)
@@ -329,6 +328,10 @@ void init_em(FILE *fplog, const char *title,
 
     init_nrnb(nrnb);
 
+    /* Interactive molecular dynamics */
+    init_IMD(ir, cr, top_global, fplog, 1, state_global->x,
+             nfile, fnm, NULL, imdport, Flags);
+
     if (DOMAINDECOMP(cr))
     {
         *top = dd_init_local_top(top_global);
@@ -369,22 +372,12 @@ void init_em(FILE *fplog, const char *title,
         }
         copy_mat(state_global->box, ems->s.box);
 
-        if (PAR(cr) && ir->eI != eiNM)
-        {
-            /* Initialize the particle decomposition and split the topology */
-            *top = split_system(fplog, top_global, ir, cr);
-
-            pd_cg_range(cr, &fr->cg0, &fr->hcg);
-        }
-        else
-        {
-            *top = gmx_mtop_generate_local_top(top_global, ir);
-        }
+        *top      = gmx_mtop_generate_local_top(top_global, ir);
         *f_global = *f;
 
-        forcerec_set_excl_load(fr, *top, cr);
+        forcerec_set_excl_load(fr, *top);
 
-        init_bonded_thread_force_reduction(fr, &(*top)->idef);
+        setup_bonded_threading(fr, &(*top)->idef);
 
         if (ir->ePBC != epbcNONE && !fr->bMolPBC)
         {
@@ -395,17 +388,7 @@ void init_em(FILE *fplog, const char *title,
             *graph = NULL;
         }
 
-        if (PARTDECOMP(cr))
-        {
-            pd_at_range(cr, &start, &homenr);
-            homenr -= start;
-        }
-        else
-        {
-            start  = 0;
-            homenr = top_global->natoms;
-        }
-        atoms2md(top_global, ir, 0, NULL, start, homenr, mdatoms);
+        atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
         update_mdatoms(mdatoms, state_global->lambda[efptFEP]);
 
         if (vsite)
@@ -433,7 +416,7 @@ void init_em(FILE *fplog, const char *title,
             /* Constrain the starting coordinates */
             dvdl_constr = 0;
             constrain(PAR(cr) ? NULL : fplog, TRUE, TRUE, constr, &(*top)->idef,
-                      ir, NULL, cr, -1, 0, mdatoms,
+                      ir, NULL, cr, -1, 0, 1.0, mdatoms,
                       ems->s.x, ems->s.x, NULL, fr->bMolPBC, ems->s.box,
                       ems->s.lambda[efptFEP], &dvdl_constr,
                       NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
@@ -445,7 +428,7 @@ void init_em(FILE *fplog, const char *title,
         *gstat = global_stat_init(ir);
     }
 
-    *outf = init_mdoutf(nfile, fnm, 0, cr, ir, NULL);
+    *outf = init_mdoutf(fplog, nfile, fnm, 0, cr, ir, top_global, NULL, wcycle);
 
     snew(*enerd, 1);
     init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
@@ -454,15 +437,16 @@ void init_em(FILE *fplog, const char *title,
     if (mdebin != NULL)
     {
         /* Init bin for energy stuff */
-        *mdebin = init_mdebin((*outf)->fp_ene, top_global, ir, NULL);
+        *mdebin = init_mdebin(mdoutf_get_fp_ene(*outf), top_global, ir, NULL);
     }
 
     clear_rvec(mu_tot);
     calc_shifts(ems->s.box, fr->shift_vec);
 }
 
-static void finish_em(FILE *fplog, t_commrec *cr, gmx_mdoutf_t *outf,
-                      gmx_runtime_t *runtime, gmx_wallcycle_t wcycle)
+static void finish_em(t_commrec *cr, gmx_mdoutf_t outf,
+                      gmx_walltime_accounting_t walltime_accounting,
+                      gmx_wallcycle_t wcycle)
 {
     if (!(cr->duty & DUTY_PME))
     {
@@ -472,7 +456,7 @@ static void finish_em(FILE *fplog, t_commrec *cr, gmx_mdoutf_t *outf,
 
     done_mdoutf(outf);
 
-    em_time_end(fplog, cr, runtime, wcycle);
+    em_time_end(walltime_accounting, wcycle);
 }
 
 static void swap_em_state(em_state_t *ems1, em_state_t *ems2)
@@ -495,16 +479,24 @@ static void copy_em_coords(em_state_t *ems, t_state *state)
 }
 
 static void write_em_traj(FILE *fplog, t_commrec *cr,
-                          gmx_mdoutf_t *outf,
+                          gmx_mdoutf_t outf,
                           gmx_bool bX, gmx_bool bF, const char *confout,
                           gmx_mtop_t *top_global,
-                          t_inputrec *ir, gmx_large_int_t step,
+                          t_inputrec *ir, gmx_int64_t step,
                           em_state_t *state,
                           t_state *state_global, rvec *f_global)
 {
-    int mdof_flags;
+    int      mdof_flags;
+    gmx_bool bIMDout = FALSE;
+
 
-    if ((bX || bF || confout != NULL) && !DOMAINDECOMP(cr))
+    /* Shall we do IMD output? */
+    if (ir->bIMD)
+    {
+        bIMDout = do_per_step(step, IMD_get_step(ir->imd->setup));
+    }
+
+    if ((bX || bF || bIMDout || confout != NULL) && !DOMAINDECOMP(cr))
     {
         copy_em_coords(state, state_global);
         f_global = state->f;
@@ -519,9 +511,16 @@ static void write_em_traj(FILE *fplog, t_commrec *cr,
     {
         mdof_flags |= MDOF_F;
     }
-    write_traj(fplog, cr, outf, mdof_flags,
-               top_global, step, (double)step,
-               &state->s, state_global, state->f, f_global, NULL, NULL);
+
+    /* If we want IMD output, set appropriate MDOF flag */
+    if (ir->bIMD)
+    {
+        mdof_flags |= MDOF_IMD;
+    }
+
+    mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
+                                     top_global, step, (double)step,
+                                     &state->s, state_global, state->f, f_global);
 
     if (confout != NULL && MASTER(cr))
     {
@@ -543,7 +542,7 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
                        em_state_t *ems1, real a, rvec *f, em_state_t *ems2,
                        gmx_constr_t constr, gmx_localtop_t *top,
                        t_nrnb *nrnb, gmx_wallcycle_t wcycle,
-                       gmx_large_int_t count)
+                       gmx_int64_t count)
 
 {
     t_state *s1, *s2;
@@ -551,6 +550,7 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
     int      start, end;
     rvec    *x1, *x2;
     real     dvdl_constr;
+    int      nthreads gmx_unused;
 
     s1 = &ems1->s;
     s2 = &ems2->s;
@@ -582,13 +582,14 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
     }
     copy_mat(s1->box, s2->box);
 
-    start = md->start;
-    end   = md->start + md->homenr;
+    start = 0;
+    end   = md->homenr;
 
     x1 = s1->x;
     x2 = s2->x;
 
-#pragma omp parallel num_threads(gmx_omp_nthreads_get(emntUpdate))
+    nthreads = gmx_omp_nthreads_get(emntUpdate);
+#pragma omp parallel num_threads(nthreads)
     {
         int gf, i, m;
 
@@ -650,7 +651,7 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
         wallcycle_start(wcycle, ewcCONSTR);
         dvdl_constr = 0;
         constrain(NULL, TRUE, TRUE, constr, &top->idef,
-                  ir, NULL, cr, count, 0, md,
+                  ir, NULL, cr, count, 0, 1.0, md,
                   s1->x, s2->x, NULL, bMolPBC, s2->box,
                   s2->lambda[efptBONDED], &dvdl_constr,
                   NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
@@ -676,8 +677,8 @@ static void em_dd_partition_system(FILE *fplog, int step, t_commrec *cr,
     wallcycle_stop(wcycle, ewcDOMDEC);
 }
 
-static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
-                            t_state *state_global, gmx_mtop_t *top_global,
+static void evaluate_energy(FILE *fplog, t_commrec *cr,
+                            gmx_mtop_t *top_global,
                             em_state_t *ems, gmx_localtop_t *top,
                             t_inputrec *inputrec,
                             t_nrnb *nrnb, gmx_wallcycle_t wcycle,
@@ -687,7 +688,7 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
                             t_graph *graph, t_mdatoms *mdatoms,
                             t_forcerec *fr, rvec mu_tot,
                             gmx_enerdata_t *enerd, tensor vir, tensor pres,
-                            gmx_large_int_t count, gmx_bool bFirst)
+                            gmx_int64_t count, gmx_bool bFirst)
 {
     real     t;
     gmx_bool bNS;
@@ -702,7 +703,7 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
     if (bFirst ||
         (DOMAINDECOMP(cr) && ems->s.ddp_count < cr->dd->ddp_count))
     {
-        /* This the first state or an old state used before the last ns */
+        /* This is the first state or an old state used before the last ns */
         bNS = TRUE;
     }
     else
@@ -725,20 +726,17 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
 
     if (vsite)
     {
-        construct_vsites(fplog, vsite, ems->s.x, nrnb, 1, NULL,
+        construct_vsites(vsite, ems->s.x, 1, NULL,
                          top->idef.iparams, top->idef.il,
-                         fr->ePBC, fr->bMolPBC, graph, cr, ems->s.box);
+                         fr->ePBC, fr->bMolPBC, cr, ems->s.box);
     }
 
-    if (DOMAINDECOMP(cr))
+    if (DOMAINDECOMP(cr) && bNS)
     {
-        if (bNS)
-        {
-            /* Repartition the domain decomposition */
-            em_dd_partition_system(fplog, count, cr, top_global, inputrec,
-                                   ems, top, mdatoms, fr, vsite, constr,
-                                   nrnb, wcycle);
-        }
+        /* Repartition the domain decomposition */
+        em_dd_partition_system(fplog, count, cr, top_global, inputrec,
+                               ems, top, mdatoms, fr, vsite, constr,
+                               nrnb, wcycle);
     }
 
     /* Calc force & energy on new trial position  */
@@ -746,7 +744,7 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
      * We do not unshift, so molecules are always whole in congrad.c
      */
     do_force(fplog, cr, inputrec,
-             count, nrnb, wcycle, top, top_global, &top_global->groups,
+             count, nrnb, wcycle, top, &top_global->groups,
              ems->s.box, ems->s.x, &ems->s.hist,
              ems->f, force_vir, mdatoms, enerd, fcd,
              ems->s.lambda, graph, fr, vsite, mu_tot, t, NULL, NULL, TRUE,
@@ -775,7 +773,7 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
     }
 
     /* Calculate long range corrections to pressure and energy */
-    calc_dispcorr(fplog, inputrec, fr, count, top_global->natoms, ems->s.box, ems->s.lambda[efptVDW],
+    calc_dispcorr(inputrec, fr, top_global->natoms, ems->s.box, ems->s.lambda[efptVDW],
                   pres, force_vir, &prescorr, &enercorr, &dvdlcorr);
     enerd->term[F_DISPCORR] = enercorr;
     enerd->term[F_EPOT]    += enercorr;
@@ -790,14 +788,10 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
         wallcycle_start(wcycle, ewcCONSTR);
         dvdl_constr = 0;
         constrain(NULL, FALSE, FALSE, constr, &top->idef,
-                  inputrec, NULL, cr, count, 0, mdatoms,
+                  inputrec, NULL, cr, count, 0, 1.0, mdatoms,
                   ems->s.x, ems->f, ems->f, fr->bMolPBC, ems->s.box,
                   ems->s.lambda[efptBONDED], &dvdl_constr,
                   NULL, &shake_vir, nrnb, econqForceDispl, FALSE, 0, 0);
-        if (fr->bSepDVDL && fplog)
-        {
-            fprintf(fplog, sepdvdlformat, "Constraints", t, dvdl_constr);
-        }
         enerd->term[F_DVDL_CONSTR] += dvdl_constr;
         m_add(force_vir, shake_vir, vir);
         wallcycle_stop(wcycle, ewcCONSTR);
@@ -920,7 +914,7 @@ static real pr_beta(t_commrec *cr, t_grpopts *opts, t_mdatoms *mdatoms,
         /* This part of code can be incorrect with DD,
          * since the atom ordering in s_b and s_min might differ.
          */
-        for (i = mdatoms->start; i < mdatoms->start+mdatoms->homenr; i++)
+        for (i = 0; i < mdatoms->homenr; i++)
         {
             if (mdatoms->cFREEZE)
             {
@@ -950,23 +944,24 @@ static real pr_beta(t_commrec *cr, t_grpopts *opts, t_mdatoms *mdatoms,
 
 double do_cg(FILE *fplog, t_commrec *cr,
              int nfile, const t_filenm fnm[],
-             const output_env_t oenv, gmx_bool bVerbose, gmx_bool bCompact,
-             int nstglobalcomm,
+             const output_env_t gmx_unused oenv, gmx_bool bVerbose, gmx_bool gmx_unused bCompact,
+             int gmx_unused nstglobalcomm,
              gmx_vsite_t *vsite, gmx_constr_t constr,
-             int stepout,
+             int gmx_unused stepout,
              t_inputrec *inputrec,
              gmx_mtop_t *top_global, t_fcdata *fcd,
              t_state *state_global,
              t_mdatoms *mdatoms,
              t_nrnb *nrnb, gmx_wallcycle_t wcycle,
-             gmx_edsam_t ed,
+             gmx_edsam_t gmx_unused ed,
              t_forcerec *fr,
-             int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
-             gmx_membed_t membed,
-             real cpt_period, real max_hours,
-             const char *deviceOptions,
-             unsigned long Flags,
-             gmx_runtime_t *runtime)
+             int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+             gmx_membed_t gmx_unused membed,
+             real gmx_unused cpt_period, real gmx_unused max_hours,
+             const char gmx_unused *deviceOptions,
+             int imdport,
+             unsigned long gmx_unused Flags,
+             gmx_walltime_accounting_t walltime_accounting)
 {
     const char       *CG = "Polak-Ribiere Conjugate Gradients";
 
@@ -989,7 +984,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
     gmx_bool          do_log = FALSE, do_ene = FALSE, do_x, do_f;
     tensor            vir, pres;
     int               number_steps, neval = 0, nstcg = inputrec->nstcgsteep;
-    gmx_mdoutf_t     *outf;
+    gmx_mdoutf_t      outf;
     int               i, m, gf, step, nminstep;
     real              terminate = 0;
 
@@ -1004,10 +999,10 @@ double do_cg(FILE *fplog, t_commrec *cr,
     init_em(fplog, CG, cr, inputrec,
             state_global, top_global, s_min, &top, &f, &f_global,
             nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
-            nfile, fnm, &outf, &mdebin);
+            nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
 
     /* Print to log file */
-    print_em_start(fplog, cr, runtime, wcycle, CG);
+    print_em_start(fplog, cr, walltime_accounting, wcycle, CG);
 
     /* Max number of steps */
     number_steps = inputrec->nsteps;
@@ -1025,8 +1020,8 @@ double do_cg(FILE *fplog, t_commrec *cr,
     /* do_force always puts the charge groups in the box and shifts again
      * We do not unshift, so molecules are always whole in congrad.c
      */
-    evaluate_energy(fplog, bVerbose, cr,
-                    state_global, top_global, s_min, top,
+    evaluate_energy(fplog, cr,
+                    top_global, s_min, top,
                     inputrec, nrnb, wcycle, gstat,
                     vsite, constr, fcd, graph, mdatoms, fr,
                     mu_tot, enerd, vir, pres, -1, TRUE);
@@ -1040,7 +1035,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
                    NULL, NULL, vir, pres, NULL, mu_tot, constr);
 
         print_ebin_header(fplog, step, step, s_min->s.lambda[efptFEP]);
-        print_ebin(outf->fp_ene, TRUE, FALSE, FALSE, fplog, step, step, eprNORMAL,
+        print_ebin(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, step, eprNORMAL,
                    TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
     }
     where();
@@ -1080,7 +1075,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
         sf  = s_min->f;
         gpa = 0;
         gf  = 0;
-        for (i = mdatoms->start; i < mdatoms->start+mdatoms->homenr; i++)
+        for (i = 0; i < mdatoms->homenr; i++)
         {
             if (mdatoms->cFREEZE)
             {
@@ -1133,7 +1128,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
          * relative change in coordinate is smaller than precision
          */
         minstep = 0;
-        for (i = mdatoms->start; i < mdatoms->start+mdatoms->homenr; i++)
+        for (i = 0; i < mdatoms->homenr; i++)
         {
             for (m = 0; m < DIM; m++)
             {
@@ -1202,8 +1197,8 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
         neval++;
         /* Calculate energy for the trial step */
-        evaluate_energy(fplog, bVerbose, cr,
-                        state_global, top_global, s_c, top,
+        evaluate_energy(fplog, cr,
+                        top_global, s_c, top,
                         inputrec, nrnb, wcycle, gstat,
                         vsite, constr, fcd, graph, mdatoms, fr,
                         mu_tot, enerd, vir, pres, -1, FALSE);
@@ -1212,7 +1207,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
         p   = s_c->s.cg_p;
         sf  = s_c->f;
         gpc = 0;
-        for (i = mdatoms->start; i < mdatoms->start+mdatoms->homenr; i++)
+        for (i = 0; i < mdatoms->homenr; i++)
         {
             for (m = 0; m < DIM; m++)
             {
@@ -1310,8 +1305,8 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
                 neval++;
                 /* Calculate energy for the trial step */
-                evaluate_energy(fplog, bVerbose, cr,
-                                state_global, top_global, s_b, top,
+                evaluate_energy(fplog, cr,
+                                top_global, s_b, top,
                                 inputrec, nrnb, wcycle, gstat,
                                 vsite, constr, fcd, graph, mdatoms, fr,
                                 mu_tot, enerd, vir, pres, -1, FALSE);
@@ -1322,7 +1317,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
                 p   = s_b->s.cg_p;
                 sf  = s_b->f;
                 gpb = 0;
-                for (i = mdatoms->start; i < mdatoms->start+mdatoms->homenr; i++)
+                for (i = 0; i < mdatoms->homenr; i++)
                 {
                     for (m = 0; m < DIM; m++)
                     {
@@ -1471,15 +1466,25 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
             do_log = do_per_step(step, inputrec->nstlog);
             do_ene = do_per_step(step, inputrec->nstenergy);
+
+            /* Prepare IMD energy record, if bIMD is TRUE. */
+            IMD_fill_energy_record(inputrec->bIMD, inputrec->imd, enerd, step, TRUE);
+
             if (do_log)
             {
                 print_ebin_header(fplog, step, step, s_min->s.lambda[efptFEP]);
             }
-            print_ebin(outf->fp_ene, do_ene, FALSE, FALSE,
+            print_ebin(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
                        do_log ? fplog : NULL, step, step, eprNORMAL,
                        TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
         }
 
+        /* Send energies and positions to the IMD client if bIMD is TRUE. */
+        if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+        {
+            IMD_send_positions(inputrec->imd);
+        }
+
         /* Stop when the maximum force lies below tolerance.
          * If we have reached machine precision, converged is already set to true.
          */
@@ -1487,6 +1492,9 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
     } /* End of the loop */
 
+    /* IMD cleanup, if bIMD is TRUE. */
+    IMD_finalize(inputrec->bIMD, inputrec->imd);
+
     if (converged)
     {
         step--; /* we never took that last step in this case */
@@ -1515,7 +1523,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
         if (!do_ene || !do_log)
         {
             /* Write final energy file entries */
-            print_ebin(outf->fp_ene, !do_ene, FALSE, FALSE,
+            print_ebin(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
                        !do_log ? fplog : NULL, step, step, eprNORMAL,
                        TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
         }
@@ -1553,10 +1561,10 @@ double do_cg(FILE *fplog, t_commrec *cr,
         fprintf(fplog, "\nPerformed %d energy evaluations in total.\n", neval);
     }
 
-    finish_em(fplog, cr, outf, runtime, wcycle);
+    finish_em(cr, outf, walltime_accounting, wcycle);
 
     /* To print the actual number of steps we needed somewhere */
-    runtime->nsteps_done = step;
+    walltime_accounting_set_nsteps_done(walltime_accounting, step);
 
     return 0;
 } /* That's all folks */
@@ -1564,23 +1572,24 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
 double do_lbfgs(FILE *fplog, t_commrec *cr,
                 int nfile, const t_filenm fnm[],
-                const output_env_t oenv, gmx_bool bVerbose, gmx_bool bCompact,
-                int nstglobalcomm,
+                const output_env_t gmx_unused oenv, gmx_bool bVerbose, gmx_bool gmx_unused bCompact,
+                int gmx_unused nstglobalcomm,
                 gmx_vsite_t *vsite, gmx_constr_t constr,
-                int stepout,
+                int gmx_unused stepout,
                 t_inputrec *inputrec,
                 gmx_mtop_t *top_global, t_fcdata *fcd,
                 t_state *state,
                 t_mdatoms *mdatoms,
                 t_nrnb *nrnb, gmx_wallcycle_t wcycle,
-                gmx_edsam_t ed,
+                gmx_edsam_t gmx_unused ed,
                 t_forcerec *fr,
-                int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
-                gmx_membed_t membed,
-                real cpt_period, real max_hours,
-                const char *deviceOptions,
-                unsigned long Flags,
-                gmx_runtime_t *runtime)
+                int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+                gmx_membed_t gmx_unused membed,
+                real gmx_unused cpt_period, real gmx_unused max_hours,
+                const char gmx_unused *deviceOptions,
+                int imdport,
+                unsigned long gmx_unused Flags,
+                gmx_walltime_accounting_t walltime_accounting)
 {
     static const char *LBFGS = "Low-Memory BFGS Minimizer";
     em_state_t         ems;
@@ -1604,7 +1613,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
     gmx_bool           do_log, do_ene, do_x, do_f, foundlower, *frozen;
     tensor             vir, pres;
     int                start, end, number_steps;
-    gmx_mdoutf_t      *outf;
+    gmx_mdoutf_t       outf;
     int                i, k, m, n, nfmax, gf, step;
     int                mdof_flags;
     /* not used */
@@ -1662,7 +1671,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
     init_em(fplog, LBFGS, cr, inputrec,
             state, top_global, &ems, &top, &f, &f_global,
             nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
-            nfile, fnm, &outf, &mdebin);
+            nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
     /* Do_lbfgs is not completely updated like do_steep and do_cg,
      * so we free some memory again.
      */
@@ -1672,11 +1681,11 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
     xx = (real *)state->x;
     ff = (real *)f;
 
-    start = mdatoms->start;
-    end   = mdatoms->homenr + start;
+    start = 0;
+    end   = mdatoms->homenr;
 
     /* Print to log file */
-    print_em_start(fplog, cr, runtime, wcycle, LBFGS);
+    print_em_start(fplog, cr, walltime_accounting, wcycle, LBFGS);
 
     do_log = do_ene = do_x = do_f = TRUE;
 
@@ -1707,9 +1716,9 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
 
     if (vsite)
     {
-        construct_vsites(fplog, vsite, state->x, nrnb, 1, NULL,
+        construct_vsites(vsite, state->x, 1, NULL,
                          top->idef.iparams, top->idef.il,
-                         fr->ePBC, fr->bMolPBC, graph, cr, state->box);
+                         fr->ePBC, fr->bMolPBC, cr, state->box);
     }
 
     /* Call the force routine and some auxiliary (neighboursearching etc.) */
@@ -1719,8 +1728,8 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
     neval++;
     ems.s.x = state->x;
     ems.f   = f;
-    evaluate_energy(fplog, bVerbose, cr,
-                    state, top_global, &ems, top,
+    evaluate_energy(fplog, cr,
+                    top_global, &ems, top,
                     inputrec, nrnb, wcycle, gstat,
                     vsite, constr, fcd, graph, mdatoms, fr,
                     mu_tot, enerd, vir, pres, -1, TRUE);
@@ -1734,7 +1743,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
                    NULL, NULL, vir, pres, NULL, mu_tot, constr);
 
         print_ebin_header(fplog, step, step, state->lambda[efptFEP]);
-        print_ebin(outf->fp_ene, TRUE, FALSE, FALSE, fplog, step, step, eprNORMAL,
+        print_ebin(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, step, eprNORMAL,
                    TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
     }
     where();
@@ -1779,7 +1788,6 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
     }
 
     stepsize  = 1.0/fnorm;
-    converged = FALSE;
 
     /* Start the loop over BFGS steps.
      * Each successful step is counted, and we continue until
@@ -1808,8 +1816,13 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
             mdof_flags |= MDOF_F;
         }
 
-        write_traj(fplog, cr, outf, mdof_flags,
-                   top_global, step, (real)step, state, state, f, f, NULL, NULL);
+        if (inputrec->bIMD)
+        {
+            mdof_flags |= MDOF_IMD;
+        }
+
+        mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
+                                         top_global, step, (real)step, state, state, f, f);
 
         /* Do the linesearching in the direction dx[point][0..(n-1)] */
 
@@ -1912,8 +1925,8 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
         /* Calculate energy for the trial step */
         ems.s.x = (rvec *)xc;
         ems.f   = (rvec *)fc;
-        evaluate_energy(fplog, bVerbose, cr,
-                        state, top_global, &ems, top,
+        evaluate_energy(fplog, cr,
+                        top_global, &ems, top,
                         inputrec, nrnb, wcycle, gstat,
                         vsite, constr, fcd, graph, mdatoms, fr,
                         mu_tot, enerd, vir, pres, step, FALSE);
@@ -2010,8 +2023,8 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
                 /* Calculate energy for the trial step */
                 ems.s.x = (rvec *)xb;
                 ems.f   = (rvec *)fb;
-                evaluate_energy(fplog, bVerbose, cr,
-                                state, top_global, &ems, top,
+                evaluate_energy(fplog, cr,
+                                top_global, &ems, top,
                                 inputrec, nrnb, wcycle, gstat,
                                 vsite, constr, fcd, graph, mdatoms, fr,
                                 mu_tot, enerd, vir, pres, step, FALSE);
@@ -2264,11 +2277,17 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
             {
                 print_ebin_header(fplog, step, step, state->lambda[efptFEP]);
             }
-            print_ebin(outf->fp_ene, do_ene, FALSE, FALSE,
+            print_ebin(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
                        do_log ? fplog : NULL, step, step, eprNORMAL,
                        TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
         }
 
+        /* Send x and E to IMD client, if bIMD is TRUE. */
+        if (do_IMD(inputrec->bIMD, step, cr, TRUE, state->box, state->x, inputrec, 0, wcycle) && MASTER(cr))
+        {
+            IMD_send_positions(inputrec->imd);
+        }
+
         /* Stop when the maximum force lies below tolerance.
          * If we have reached machine precision, converged is already set to true.
          */
@@ -2277,6 +2296,9 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
 
     } /* End of the loop */
 
+    /* IMD cleanup, if bIMD is TRUE. */
+    IMD_finalize(inputrec->bIMD, inputrec->imd);
+
     if (converged)
     {
         step--; /* we never took that last step in this case */
@@ -2301,7 +2323,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
     }
     if (!do_ene || !do_log) /* Write final energy file entries */
     {
-        print_ebin(outf->fp_ene, !do_ene, FALSE, FALSE,
+        print_ebin(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
                    !do_log ? fplog : NULL, step, step, eprNORMAL,
                    TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
     }
@@ -2335,10 +2357,10 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
         fprintf(fplog, "\nPerformed %d energy evaluations in total.\n", neval);
     }
 
-    finish_em(fplog, cr, outf, runtime, wcycle);
+    finish_em(cr, outf, walltime_accounting, wcycle);
 
     /* To print the actual number of steps we needed somewhere */
-    runtime->nsteps_done = step;
+    walltime_accounting_set_nsteps_done(walltime_accounting, step);
 
     return 0;
 } /* That's all folks */
@@ -2346,23 +2368,24 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
 
 double do_steep(FILE *fplog, t_commrec *cr,
                 int nfile, const t_filenm fnm[],
-                const output_env_t oenv, gmx_bool bVerbose, gmx_bool bCompact,
-                int nstglobalcomm,
+                const output_env_t gmx_unused oenv, gmx_bool bVerbose, gmx_bool gmx_unused bCompact,
+                int gmx_unused nstglobalcomm,
                 gmx_vsite_t *vsite, gmx_constr_t constr,
-                int stepout,
+                int gmx_unused stepout,
                 t_inputrec *inputrec,
                 gmx_mtop_t *top_global, t_fcdata *fcd,
                 t_state *state_global,
                 t_mdatoms *mdatoms,
                 t_nrnb *nrnb, gmx_wallcycle_t wcycle,
-                gmx_edsam_t ed,
+                gmx_edsam_t gmx_unused  ed,
                 t_forcerec *fr,
-                int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
-                gmx_membed_t membed,
-                real cpt_period, real max_hours,
-                const char *deviceOptions,
-                unsigned long Flags,
-                gmx_runtime_t *runtime)
+                int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+                gmx_membed_t gmx_unused membed,
+                real gmx_unused cpt_period, real gmx_unused max_hours,
+                const char  gmx_unused *deviceOptions,
+                int imdport,
+                unsigned long gmx_unused Flags,
+                gmx_walltime_accounting_t walltime_accounting)
 {
     const char       *SD = "Steepest Descents";
     em_state_t       *s_min, *s_try;
@@ -2374,7 +2397,7 @@ double do_steep(FILE *fplog, t_commrec *cr,
     t_graph          *graph;
     real              stepsize, constepsize;
     real              ustep, fnormn;
-    gmx_mdoutf_t     *outf;
+    gmx_mdoutf_t      outf;
     t_mdebin         *mdebin;
     gmx_bool          bDone, bAbort, do_x, do_f;
     tensor            vir, pres;
@@ -2392,10 +2415,10 @@ double do_steep(FILE *fplog, t_commrec *cr,
     init_em(fplog, SD, cr, inputrec,
             state_global, top_global, s_try, &top, &f, &f_global,
             nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
-            nfile, fnm, &outf, &mdebin);
+            nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
 
     /* Print to log file  */
-    print_em_start(fplog, cr, runtime, wcycle, SD);
+    print_em_start(fplog, cr, walltime_accounting, wcycle, SD);
 
     /* Set variables for stepsize (in nm). This is the largest
      * step that we are going to make in any direction.
@@ -2437,8 +2460,8 @@ double do_steep(FILE *fplog, t_commrec *cr,
                        constr, top, nrnb, wcycle, count);
         }
 
-        evaluate_energy(fplog, bVerbose, cr,
-                        state_global, top_global, s_try, top,
+        evaluate_energy(fplog, cr,
+                        top_global, s_try, top,
                         inputrec, nrnb, wcycle, gstat,
                         vsite, constr, fcd, graph, mdatoms, fr,
                         mu_tot, enerd, vir, pres, count, count == 0);
@@ -2469,7 +2492,11 @@ double do_steep(FILE *fplog, t_commrec *cr,
                 upd_mdebin(mdebin, FALSE, FALSE, (double)count,
                            mdatoms->tmass, enerd, &s_try->s, inputrec->fepvals, inputrec->expandedvals,
                            s_try->s.box, NULL, NULL, vir, pres, NULL, mu_tot, constr);
-                print_ebin(outf->fp_ene, TRUE,
+
+                /* Prepare IMD energy record, if bIMD is TRUE. */
+                IMD_fill_energy_record(inputrec->bIMD, inputrec->imd, enerd, count, TRUE);
+
+                print_ebin(mdoutf_get_fp_ene(outf), TRUE,
                            do_per_step(steps_accepted, inputrec->nstdisreout),
                            do_per_step(steps_accepted, inputrec->nstorireout),
                            fplog, count, count, eprNORMAL, TRUE,
@@ -2538,10 +2565,19 @@ double do_steep(FILE *fplog, t_commrec *cr,
             bAbort = TRUE;
         }
 
+        /* Send IMD energies and positions, if bIMD is TRUE. */
+        if (do_IMD(inputrec->bIMD, count, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+        {
+            IMD_send_positions(inputrec->imd);
+        }
+
         count++;
     } /* End of the loop  */
 
-    /* Print some shit...  */
+    /* IMD cleanup, if bIMD is TRUE. */
+    IMD_finalize(inputrec->bIMD, inputrec->imd);
+
+    /* Print some data...  */
     if (MASTER(cr))
     {
         fprintf(stderr, "\nwriting lowest energy coordinates.\n");
@@ -2560,12 +2596,12 @@ double do_steep(FILE *fplog, t_commrec *cr,
                         s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
     }
 
-    finish_em(fplog, cr, outf, runtime, wcycle);
+    finish_em(cr, outf, walltime_accounting, wcycle);
 
     /* To print the actual number of steps we needed somewhere */
     inputrec->nsteps = count;
 
-    runtime->nsteps_done = count;
+    walltime_accounting_set_nsteps_done(walltime_accounting, count);
 
     return 0;
 } /* That's all folks */
@@ -2573,26 +2609,27 @@ double do_steep(FILE *fplog, t_commrec *cr,
 
 double do_nm(FILE *fplog, t_commrec *cr,
              int nfile, const t_filenm fnm[],
-             const output_env_t oenv, gmx_bool bVerbose, gmx_bool bCompact,
-             int nstglobalcomm,
+             const output_env_t gmx_unused oenv, gmx_bool bVerbose, gmx_bool gmx_unused  bCompact,
+             int gmx_unused nstglobalcomm,
              gmx_vsite_t *vsite, gmx_constr_t constr,
-             int stepout,
+             int gmx_unused stepout,
              t_inputrec *inputrec,
              gmx_mtop_t *top_global, t_fcdata *fcd,
              t_state *state_global,
              t_mdatoms *mdatoms,
              t_nrnb *nrnb, gmx_wallcycle_t wcycle,
-             gmx_edsam_t ed,
+             gmx_edsam_t  gmx_unused ed,
              t_forcerec *fr,
-             int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
-             gmx_membed_t membed,
-             real cpt_period, real max_hours,
-             const char *deviceOptions,
-             unsigned long Flags,
-             gmx_runtime_t *runtime)
+             int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+             gmx_membed_t gmx_unused membed,
+             real gmx_unused cpt_period, real gmx_unused max_hours,
+             const char gmx_unused *deviceOptions,
+             int imdport,
+             unsigned long gmx_unused Flags,
+             gmx_walltime_accounting_t walltime_accounting)
 {
     const char          *NM = "Normal Mode Analysis";
-    gmx_mdoutf_t        *outf;
+    gmx_mdoutf_t         outf;
     int                  natoms, atom, d;
     int                  nnodes, node;
     rvec                *f_global;
@@ -2607,7 +2644,7 @@ double do_nm(FILE *fplog, t_commrec *cr,
     rvec                 mu_tot;
     rvec                *fneg, *dfdx;
     gmx_bool             bSparse; /* use sparse matrix storage format */
-    size_t               sz;
+    size_t               sz = 0;
     gmx_sparsematrix_t * sparse_matrix           = NULL;
     real           *     full_matrix             = NULL;
     em_state_t       *   state_work;
@@ -2630,7 +2667,7 @@ double do_nm(FILE *fplog, t_commrec *cr,
             state_global, top_global, state_work, &top,
             &f, &f_global,
             nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
-            nfile, fnm, &outf, NULL);
+            nfile, fnm, &outf, NULL, imdport, Flags, wcycle);
 
     natoms = top_global->natoms;
     snew(fneg, natoms);
@@ -2655,32 +2692,35 @@ double do_nm(FILE *fplog, t_commrec *cr,
      */
     if (EEL_FULL(fr->eeltype) || fr->rlist == 0.0)
     {
-        fprintf(stderr, "Non-cutoff electrostatics used, forcing full Hessian format.\n");
+        md_print_info(cr, fplog, "Non-cutoff electrostatics used, forcing full Hessian format.\n");
         bSparse = FALSE;
     }
     else if (top_global->natoms < 1000)
     {
-        fprintf(stderr, "Small system size (N=%d), using full Hessian format.\n", top_global->natoms);
+        md_print_info(cr, fplog, "Small system size (N=%d), using full Hessian format.\n", top_global->natoms);
         bSparse = FALSE;
     }
     else
     {
-        fprintf(stderr, "Using compressed symmetric sparse Hessian format.\n");
+        md_print_info(cr, fplog, "Using compressed symmetric sparse Hessian format.\n");
         bSparse = TRUE;
     }
 
-    sz = DIM*top_global->natoms;
+    if (MASTER(cr))
+    {
+        sz = DIM*top_global->natoms;
 
-    fprintf(stderr, "Allocating Hessian memory...\n\n");
+        fprintf(stderr, "Allocating Hessian memory...\n\n");
 
-    if (bSparse)
-    {
-        sparse_matrix = gmx_sparsematrix_init(sz);
-        sparse_matrix->compressed_symmetric = TRUE;
-    }
-    else
-    {
-        snew(full_matrix, sz*sz);
+        if (bSparse)
+        {
+            sparse_matrix = gmx_sparsematrix_init(sz);
+            sparse_matrix->compressed_symmetric = TRUE;
+        }
+        else
+        {
+            snew(full_matrix, sz*sz);
+        }
     }
 
     /* Initial values */
@@ -2694,7 +2734,7 @@ double do_nm(FILE *fplog, t_commrec *cr,
     where();
 
     /* Write start time and temperature */
-    print_em_start(fplog, cr, runtime, wcycle, NM);
+    print_em_start(fplog, cr, walltime_accounting, wcycle, NM);
 
     /* fudge nr of steps to nr of atoms */
     inputrec->nsteps = natoms*2;
@@ -2709,8 +2749,8 @@ double do_nm(FILE *fplog, t_commrec *cr,
 
     /* Make evaluate_energy do a single node force calculation */
     cr->nnodes = 1;
-    evaluate_energy(fplog, bVerbose, cr,
-                    state_global, top_global, state_work, top,
+    evaluate_energy(fplog, cr,
+                    top_global, state_work, top,
                     inputrec, nrnb, wcycle, gstat,
                     vsite, constr, fcd, graph, mdatoms, fr,
                     mu_tot, enerd, vir, pres, -1, TRUE);
@@ -2719,16 +2759,14 @@ double do_nm(FILE *fplog, t_commrec *cr,
     /* if forces are not small, warn user */
     get_state_f_norm_max(cr, &(inputrec->opts), mdatoms, state_work);
 
-    if (MASTER(cr))
+    md_print_info(cr, fplog, "Maximum force:%12.5e\n", state_work->fmax);
+    if (state_work->fmax > 1.0e-3)
     {
-        fprintf(stderr, "Maximum force:%12.5e\n", state_work->fmax);
-        if (state_work->fmax > 1.0e-3)
-        {
-            fprintf(stderr, "Maximum force probably not small enough to");
-            fprintf(stderr, " ensure that you are in an \nenergy well. ");
-            fprintf(stderr, "Be aware that negative eigenvalues may occur");
-            fprintf(stderr, " when the\nresulting matrix is diagonalized.\n");
-        }
+        md_print_info(cr, fplog,
+                      "The force is probably not small enough to "
+                      "ensure that you are at a minimum.\n"
+                      "Be aware that negative eigenvalues may occur\n"
+                      "when the resulting matrix is diagonalized.\n\n");
     }
 
     /***********************************************************
@@ -2752,8 +2790,8 @@ double do_nm(FILE *fplog, t_commrec *cr,
 
             /* Make evaluate_energy do a single node force calculation */
             cr->nnodes = 1;
-            evaluate_energy(fplog, bVerbose, cr,
-                            state_global, top_global, state_work, top,
+            evaluate_energy(fplog, cr,
+                            top_global, state_work, top,
                             inputrec, nrnb, wcycle, gstat,
                             vsite, constr, fcd, graph, mdatoms, fr,
                             mu_tot, enerd, vir, pres, atom*2, FALSE);
@@ -2765,8 +2803,8 @@ double do_nm(FILE *fplog, t_commrec *cr,
 
             state_work->s.x[atom][d] = x_min + der_range;
 
-            evaluate_energy(fplog, bVerbose, cr,
-                            state_global, top_global, state_work, top,
+            evaluate_energy(fplog, cr,
+                            top_global, state_work, top,
                             inputrec, nrnb, wcycle, gstat,
                             vsite, constr, fcd, graph, mdatoms, fr,
                             mu_tot, enerd, vir, pres, atom*2+1, FALSE);
@@ -2855,9 +2893,9 @@ double do_nm(FILE *fplog, t_commrec *cr,
         gmx_mtxio_write(ftp2fn(efMTX, nfile, fnm), sz, sz, full_matrix, sparse_matrix);
     }
 
-    finish_em(fplog, cr, outf, runtime, wcycle);
+    finish_em(cr, outf, walltime_accounting, wcycle);
 
-    runtime->nsteps_done = natoms*2;
+    walltime_accounting_set_nsteps_done(walltime_accounting, natoms*2);
 
     return 0;
 }