Remove unnecessary config.h includes
[alexxy/gromacs.git] / src / gromacs / mdlib / stat.c
index 88dbfef23d5e410620e159ccd2548fe6bf4e53f6..e0d2ee19c60aef83fb58f280b28ab7ef5711709b 100644 (file)
@@ -1,74 +1,67 @@
-/* -*- 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.
+ *
+ * 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.
+ *
+ * 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.
+ *
  * 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.
- * 
- * For more info, check our website at http://www.gromacs.org
- * 
- * And Hey:
- * GROningen Mixture of Alchemy and Childrens' Stories
+ * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
 
 #include <string.h>
 #include <stdio.h>
-#include "typedefs.h"
-#include "sysstuff.h"
-#include "gmx_fatal.h"
-#include "network.h"
-#include "txtdump.h"
-#include "names.h"
-#include "physics.h"
-#include "vec.h"
-#include "maths.h"
-#include "mvdata.h"
-#include "main.h"
-#include "force.h"
-#include "vcm.h"
-#include "smalloc.h"
-#include "futil.h"
-#include "network.h"
-#include "rbin.h"
-#include "tgroup.h"
-#include "xtcio.h"
-#include "gmxfio.h"
-#include "trnio.h"
-#include "statutil.h"
-#include "domdec.h"
-#include "partdec.h"
-#include "constr.h"
-#include "checkpoint.h"
-#include "xvgr.h"
-#include "md_support.h"
-#include "mdrun.h"
-#include "sim_util.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/vcm.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/rbin.h"
+#include "gromacs/legacyheaders/tgroup.h"
+#include "gromacs/fileio/xtcio.h"
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/trnio.h"
+#include "gromacs/legacyheaders/domdec.h"
+#include "gromacs/legacyheaders/constr.h"
+#include "gromacs/legacyheaders/checkpoint.h"
+#include "gromacs/legacyheaders/md_support.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/sim_util.h"
 
 typedef struct gmx_global_stat
 {
@@ -81,11 +74,11 @@ gmx_global_stat_t global_stat_init(t_inputrec *ir)
 {
     gmx_global_stat_t gs;
 
-    snew(gs,1);
-    
+    snew(gs, 1);
+
     gs->rb = mk_bin();
-    snew(gs->itc0,ir->opts.ngtc);
-    snew(gs->itc1,ir->opts.ngtc);
+    snew(gs->itc0, ir->opts.ngtc);
+    snew(gs->itc1, ir->opts.ngtc);
 
     return gs;
 }
@@ -99,12 +92,13 @@ void global_stat_destroy(gmx_global_stat_t gs)
 }
 
 static int filter_enerdterm(real *afrom, gmx_bool bToBuffer, real *ato,
-                            gmx_bool bTemp, gmx_bool bPres, gmx_bool bEner) {
-    int i,to,from;
+                            gmx_bool bTemp, gmx_bool bPres, gmx_bool bEner)
+{
+    int i, to, from;
 
     from = 0;
     to   = 0;
-    for (i=0;i<F_NRE;i++)
+    for (i = 0; i < F_NRE; i++)
     {
         if (bToBuffer)
         {
@@ -114,601 +108,316 @@ static int filter_enerdterm(real *afrom, gmx_bool bToBuffer, real *ato,
         {
             to = i;
         }
-        switch (i) {
-        case F_EKIN:
-        case F_TEMP:
-        case F_DKDL:
-            if (bTemp)
-            {
-                ato[to++] = afrom[from++];
-            }
-            break;
-        case F_PRES:    
-        case F_PDISPCORR:
-        case F_VTEMP:
-            if (bPres)
-            {
-                ato[to++] = afrom[from++];
-            }
-            break;
-        default:
-            if (bEner)
-            {
-                ato[to++] = afrom[from++];
-            }
-            break;
+        switch (i)
+        {
+            case F_EKIN:
+            case F_TEMP:
+            case F_DKDL:
+                if (bTemp)
+                {
+                    ato[to++] = afrom[from++];
+                }
+                break;
+            case F_PRES:
+            case F_PDISPCORR:
+                if (bPres)
+                {
+                    ato[to++] = afrom[from++];
+                }
+                break;
+            default:
+                if (bEner)
+                {
+                    ato[to++] = afrom[from++];
+                }
+                break;
         }
     }
 
     return to;
 }
 
-void global_stat(FILE *fplog,gmx_global_stat_t gs,
-                 t_commrec *cr,gmx_enerdata_t *enerd,
-                 tensor fvir,tensor svir,rvec mu_tot,
+void global_stat(FILE *fplog, gmx_global_stat_t gs,
+                 t_commrec *cr, gmx_enerdata_t *enerd,
+                 tensor fvir, tensor svir, rvec mu_tot,
                  t_inputrec *inputrec,
-                 gmx_ekindata_t *ekind,gmx_constr_t constr,
+                 gmx_ekindata_t *ekind, gmx_constr_t constr,
                  t_vcm *vcm,
-                 int nsig,real *sig,
-                 gmx_mtop_t *top_global, t_state *state_local, 
+                 int nsig, real *sig,
+                 gmx_mtop_t *top_global, t_state *state_local,
                  gmx_bool bSumEkinhOld, int flags)
 /* instead of current system, gmx_booleans for summing virial, kinetic energy, and other terms */
 {
-  t_bin  *rb;
-  int    *itc0,*itc1;
-  int    ie=0,ifv=0,isv=0,irmsd=0,imu=0;
-  int    idedl=0,idvdll=0,idvdlnl=0,iepl=0,icm=0,imass=0,ica=0,inb=0;
-  int    isig=-1;
-  int    icj=-1,ici=-1,icx=-1;
-  int    inn[egNR];
-  real   copyenerd[F_NRE];
-  int    nener,j;
-  real   *rmsd_data=NULL;
-  double nb;
-  gmx_bool   bVV,bTemp,bEner,bPres,bConstrVir,bEkinAveVel,bFirstIterate,bReadEkin;
-
-  bVV           = EI_VV(inputrec->eI);
-  bTemp         = flags & CGLO_TEMPERATURE;
-  bEner         = flags & CGLO_ENERGY;
-  bPres         = (flags & CGLO_PRESSURE); 
-  bConstrVir    = (flags & CGLO_CONSTRAINT);
-  bFirstIterate = (flags & CGLO_FIRSTITERATE);
-  bEkinAveVel   = (inputrec->eI==eiVV || (inputrec->eI==eiVVAK && bPres));
-  bReadEkin     = (flags & CGLO_READEKIN);
-
-  rb   = gs->rb;
-  itc0 = gs->itc0;
-  itc1 = gs->itc1;
-  
-
-  reset_bin(rb);
-  /* This routine copies all the data to be summed to one big buffer
-   * using the t_bin struct. 
-   */
-
-  /* First, we neeed to identify which enerd->term should be
-     communicated.  Temperature and pressure terms should only be
-     communicated and summed when they need to be, to avoid repeating
-     the sums and overcounting. */
-
-  nener = filter_enerdterm(enerd->term,TRUE,copyenerd,bTemp,bPres,bEner);
-  
-  /* First, the data that needs to be communicated with velocity verlet every time
-     This is just the constraint virial.*/
-  if (bConstrVir) {
-      isv = add_binr(rb,DIM*DIM,svir[0]);
-      where();
-  }
-  
-/* We need the force virial and the kinetic energy for the first time through with velocity verlet */
-  if (bTemp || !bVV)
-  {
-      if (ekind) 
-      {
-          for(j=0; (j<inputrec->opts.ngtc); j++) 
-          {
-              if (bSumEkinhOld) 
-              {
-                  itc0[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinh_old[0]);
-              }
-              if (bEkinAveVel && !bReadEkin) 
-              {
-                  itc1[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinf[0]);
-              } 
-              else if (!bReadEkin)
-              {
-                  itc1[j]=add_binr(rb,DIM*DIM,ekind->tcstat[j].ekinh[0]);
-              }
-          }
-          /* these probably need to be put into one of these categories */
-          where();
-          idedl = add_binr(rb,1,&(ekind->dekindl));
-          where();
-          ica   = add_binr(rb,1,&(ekind->cosacc.mvcos));
-          where();
-      }  
-  }      
-  where();
-  
-  if ((bPres || !bVV) && bFirstIterate)
-  {
-      ifv = add_binr(rb,DIM*DIM,fvir[0]);
-  }
-
-
-  if (bEner) 
-  { 
-      where();
-      if (bFirstIterate) 
-      {
-          ie  = add_binr(rb,nener,copyenerd);
-      }
-      where();
-      if (constr) 
-      {
-          rmsd_data = constr_rmsd_data(constr);
-          if (rmsd_data) 
-          {
-              irmsd = add_binr(rb,inputrec->eI==eiSD2 ? 3 : 2,rmsd_data);
-          }
-      } 
-      if (!NEED_MUTOT(*inputrec)) 
-      {
-          imu = add_binr(rb,DIM,mu_tot);
-          where();
-      }
-      
-      if (bFirstIterate) 
-      {
-          for(j=0; (j<egNR); j++)
-          {
-              inn[j]=add_binr(rb,enerd->grpp.nener,enerd->grpp.ener[j]);
-          }
-          where();
-          if (inputrec->efep != efepNO) 
-          {
-              idvdll  = add_bind(rb,efptNR,enerd->dvdl_lin);
-              idvdlnl = add_bind(rb,efptNR,enerd->dvdl_nonlin);
-              if (enerd->n_lambda > 0) 
-              {
-                  iepl = add_bind(rb,enerd->n_lambda,enerd->enerpart_lambda);
-              }
-          }
-      }
-  }
-
-  if (vcm)
-  {
-      icm   = add_binr(rb,DIM*vcm->nr,vcm->group_p[0]);
-      where();
-      imass = add_binr(rb,vcm->nr,vcm->group_mass);
-      where();
-      if (vcm->mode == ecmANGULAR)
-      {
-          icj   = add_binr(rb,DIM*vcm->nr,vcm->group_j[0]);
-          where();
-          icx   = add_binr(rb,DIM*vcm->nr,vcm->group_x[0]);
-          where();
-          ici   = add_binr(rb,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
-          where();
-      }
-  }
-
-  if (DOMAINDECOMP(cr)) 
-  {
-      nb = cr->dd->nbonded_local;
-      inb = add_bind(rb,1,&nb);
-      }
-  where();
-  if (nsig > 0) 
-  {
-      isig = add_binr(rb,nsig,sig);
-  }
-
-  /* Global sum it all */
-  if (debug)
-  {
-      fprintf(debug,"Summing %d energies\n",rb->maxreal);
-  }
-  sum_bin(rb,cr);
-  where();
-
-  /* Extract all the data locally */
-
-  if (bConstrVir) 
-  {
-      extract_binr(rb,isv ,DIM*DIM,svir[0]);
-  }
-
-  /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
-  if (bTemp || !bVV)
-  {
-      if (ekind) 
-      {
-          for(j=0; (j<inputrec->opts.ngtc); j++) 
-          {
-              if (bSumEkinhOld)
-              {
-                  extract_binr(rb,itc0[j],DIM*DIM,ekind->tcstat[j].ekinh_old[0]);
-              }
-              if (bEkinAveVel && !bReadEkin) {
-                  extract_binr(rb,itc1[j],DIM*DIM,ekind->tcstat[j].ekinf[0]);
-              }
-              else if (!bReadEkin)
-              {
-                  extract_binr(rb,itc1[j],DIM*DIM,ekind->tcstat[j].ekinh[0]);              
-              }
-          }
-          extract_binr(rb,idedl,1,&(ekind->dekindl));
-          extract_binr(rb,ica,1,&(ekind->cosacc.mvcos));
-          where();
-      }
-  }
-  if ((bPres || !bVV) && bFirstIterate)
-  {
-      extract_binr(rb,ifv ,DIM*DIM,fvir[0]);
-  }
-
-  if (bEner) 
-  {
-      if (bFirstIterate) 
-      {
-          extract_binr(rb,ie,nener,copyenerd);
-          if (rmsd_data) 
-          {
-              extract_binr(rb,irmsd,inputrec->eI==eiSD2 ? 3 : 2,rmsd_data);
-          }
-          if (!NEED_MUTOT(*inputrec))
-          {
-              extract_binr(rb,imu,DIM,mu_tot);
-          }
-
-          for(j=0; (j<egNR); j++)
-          {
-              extract_binr(rb,inn[j],enerd->grpp.nener,enerd->grpp.ener[j]);
-          }
-          if (inputrec->efep != efepNO) 
-          {
-              extract_bind(rb,idvdll ,efptNR,enerd->dvdl_lin);
-              extract_bind(rb,idvdlnl,efptNR,enerd->dvdl_nonlin);
-              if (enerd->n_lambda > 0) 
-              {
-                  extract_bind(rb,iepl,enerd->n_lambda,enerd->enerpart_lambda);
-              }
-          }
-          if (DOMAINDECOMP(cr)) 
-          {
-              extract_bind(rb,inb,1,&nb);
-              if ((int)(nb + 0.5) != cr->dd->nbonded_global) 
-              {
-                  dd_print_missing_interactions(fplog,cr,(int)(nb + 0.5),top_global,state_local);
-              }
-          }
-          where();
-
-          filter_enerdterm(copyenerd,FALSE,enerd->term,bTemp,bPres,bEner);    
-      }
-  }
-
-  if (vcm)
-  {
-      extract_binr(rb,icm,DIM*vcm->nr,vcm->group_p[0]);
-      where();
-      extract_binr(rb,imass,vcm->nr,vcm->group_mass);
-      where();
-      if (vcm->mode == ecmANGULAR)
-      {
-          extract_binr(rb,icj,DIM*vcm->nr,vcm->group_j[0]);
-          where();
-          extract_binr(rb,icx,DIM*vcm->nr,vcm->group_x[0]);
-          where();
-          extract_binr(rb,ici,DIM*DIM*vcm->nr,vcm->group_i[0][0]);
-          where();
-      }
-  }
-
-  if (nsig > 0) 
-  {
-      extract_binr(rb,isig,nsig,sig);
-  }
-  where();
-}
+    t_bin     *rb;
+    int       *itc0, *itc1;
+    int        ie    = 0, ifv = 0, isv = 0, irmsd = 0, imu = 0;
+    int        idedl = 0, idvdll = 0, idvdlnl = 0, iepl = 0, icm = 0, imass = 0, ica = 0, inb = 0;
+    int        isig  = -1;
+    int        icj   = -1, ici = -1, icx = -1;
+    int        inn[egNR];
+    real       copyenerd[F_NRE];
+    int        nener, j;
+    real      *rmsd_data = NULL;
+    double     nb;
+    gmx_bool   bVV, bTemp, bEner, bPres, bConstrVir, bEkinAveVel, bFirstIterate, bReadEkin;
+
+    bVV           = EI_VV(inputrec->eI);
+    bTemp         = flags & CGLO_TEMPERATURE;
+    bEner         = flags & CGLO_ENERGY;
+    bPres         = (flags & CGLO_PRESSURE);
+    bConstrVir    = (flags & CGLO_CONSTRAINT);
+    bFirstIterate = (flags & CGLO_FIRSTITERATE);
+    bEkinAveVel   = (inputrec->eI == eiVV || (inputrec->eI == eiVVAK && bPres));
+    bReadEkin     = (flags & CGLO_READEKIN);
+
+    rb   = gs->rb;
+    itc0 = gs->itc0;
+    itc1 = gs->itc1;
+
+
+    reset_bin(rb);
+    /* This routine copies all the data to be summed to one big buffer
+     * using the t_bin struct.
+     */
+
+    /* First, we neeed to identify which enerd->term should be
+       communicated.  Temperature and pressure terms should only be
+       communicated and summed when they need to be, to avoid repeating
+       the sums and overcounting. */
+
+    nener = filter_enerdterm(enerd->term, TRUE, copyenerd, bTemp, bPres, bEner);
+
+    /* First, the data that needs to be communicated with velocity verlet every time
+       This is just the constraint virial.*/
+    if (bConstrVir)
+    {
+        isv = add_binr(rb, DIM*DIM, svir[0]);
+        where();
+    }
 
-int do_per_step(gmx_large_int_t step,gmx_large_int_t nstep)
-{
-  if (nstep != 0) 
-    return ((step % nstep)==0); 
-  else 
-    return 0;
-}
+/* We need the force virial and the kinetic energy for the first time through with velocity verlet */
+    if (bTemp || !bVV)
+    {
+        if (ekind)
+        {
+            for (j = 0; (j < inputrec->opts.ngtc); j++)
+            {
+                if (bSumEkinhOld)
+                {
+                    itc0[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinh_old[0]);
+                }
+                if (bEkinAveVel && !bReadEkin)
+                {
+                    itc1[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinf[0]);
+                }
+                else if (!bReadEkin)
+                {
+                    itc1[j] = add_binr(rb, DIM*DIM, ekind->tcstat[j].ekinh[0]);
+                }
+            }
+            /* these probably need to be put into one of these categories */
+            where();
+            idedl = add_binr(rb, 1, &(ekind->dekindl));
+            where();
+            ica   = add_binr(rb, 1, &(ekind->cosacc.mvcos));
+            where();
+        }
+    }
+    where();
 
-static void moveit(t_commrec *cr,
-                  int left,int right,const char *s,rvec xx[])
-{
-  if (!xx) 
-    return;
+    if ((bPres || !bVV) && bFirstIterate)
+    {
+        ifv = add_binr(rb, DIM*DIM, fvir[0]);
+    }
 
-  move_rvecs(cr,FALSE,FALSE,left,right,
-            xx,NULL,(cr->nnodes-cr->npmenodes)-1,NULL);
-}
 
-gmx_mdoutf_t *init_mdoutf(int nfile,const t_filenm fnm[],int mdrun_flags,
-                          const t_commrec *cr,const t_inputrec *ir,
-                          const output_env_t oenv)
-{
-    gmx_mdoutf_t *of;
-    char filemode[3];
-    gmx_bool bAppendFiles;
-
-    snew(of,1);
-
-    of->fp_trn   = NULL;
-    of->fp_ene   = NULL;
-    of->fp_xtc   = NULL;
-    of->fp_dhdl  = NULL;
-    of->fp_field = NULL;
-    
-    of->eIntegrator     = ir->eI;
-    of->bExpanded       = ir->bExpanded;
-    of->elamstats       = ir->expandedvals->elamstats;
-    of->simulation_part = ir->simulation_part;
-
-    if (MASTER(cr))
+    if (bEner)
     {
-        bAppendFiles = (mdrun_flags & MD_APPENDFILES);
-
-        of->bKeepAndNumCPT = (mdrun_flags & MD_KEEPANDNUMCPT);
-
-        sprintf(filemode, bAppendFiles ? "a+" : "w+");  
-        
-        if ((EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
-#ifndef GMX_FAHCORE
-            &&
-            !(EI_DYNAMICS(ir->eI) &&
-              ir->nstxout == 0 &&
-              ir->nstvout == 0 &&
-              ir->nstfout == 0)
-#endif
-           )
+        where();
+        if (bFirstIterate)
         {
-            of->fp_trn = open_trn(ftp2fn(efTRN,nfile,fnm), filemode);
+            ie  = add_binr(rb, nener, copyenerd);
         }
-        if (EI_DYNAMICS(ir->eI) &&
-            ir->nstxtcout > 0)
+        where();
+        if (constr)
         {
-            of->fp_xtc = open_xtc(ftp2fn(efXTC,nfile,fnm), filemode);
-            of->xtc_prec = ir->xtcprec;
+            rmsd_data = constr_rmsd_data(constr);
+            if (rmsd_data)
+            {
+                irmsd = add_binr(rb, inputrec->eI == eiSD2 ? 3 : 2, rmsd_data);
+            }
         }
-        if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
+        if (!NEED_MUTOT(*inputrec))
         {
-            of->fp_ene = open_enx(ftp2fn(efEDR,nfile,fnm), filemode);
+            imu = add_binr(rb, DIM, mu_tot);
+            where();
         }
-        of->fn_cpt = opt2fn("-cpo",nfile,fnm);
-        
-        if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
-            (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
-            EI_DYNAMICS(ir->eI))
+
+        if (bFirstIterate)
         {
-            if (bAppendFiles)
-            {
-                of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl",nfile,fnm),filemode);
-            }
-            else
+            for (j = 0; (j < egNR); j++)
             {
-                of->fp_dhdl = open_dhdl(opt2fn("-dhdl",nfile,fnm),ir,oenv);
+                inn[j] = add_binr(rb, enerd->grpp.nener, enerd->grpp.ener[j]);
             }
-        }
-        
-        if (opt2bSet("-field",nfile,fnm) &&
-            (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
-        {
-            if (bAppendFiles)
+            where();
+            if (inputrec->efep != efepNO)
             {
-                of->fp_dhdl = gmx_fio_fopen(opt2fn("-field",nfile,fnm),
-                                            filemode);
-            }
-            else
-            {                            
-                of->fp_field = xvgropen(opt2fn("-field",nfile,fnm),
-                                        "Applied electric field","Time (ps)",
-                                        "E (V/nm)",oenv);
+                idvdll  = add_bind(rb, efptNR, enerd->dvdl_lin);
+                idvdlnl = add_bind(rb, efptNR, enerd->dvdl_nonlin);
+                if (enerd->n_lambda > 0)
+                {
+                    iepl = add_bind(rb, enerd->n_lambda, enerd->enerpart_lambda);
+                }
             }
         }
     }
 
-    return of;
-}
-
-void done_mdoutf(gmx_mdoutf_t *of)
-{
-    if (of->fp_ene != NULL)
+    if (vcm)
     {
-        close_enx(of->fp_ene);
+        icm   = add_binr(rb, DIM*vcm->nr, vcm->group_p[0]);
+        where();
+        imass = add_binr(rb, vcm->nr, vcm->group_mass);
+        where();
+        if (vcm->mode == ecmANGULAR)
+        {
+            icj   = add_binr(rb, DIM*vcm->nr, vcm->group_j[0]);
+            where();
+            icx   = add_binr(rb, DIM*vcm->nr, vcm->group_x[0]);
+            where();
+            ici   = add_binr(rb, DIM*DIM*vcm->nr, vcm->group_i[0][0]);
+            where();
+        }
     }
-    if (of->fp_xtc)
+
+    if (DOMAINDECOMP(cr))
     {
-        close_xtc(of->fp_xtc);
+        nb  = cr->dd->nbonded_local;
+        inb = add_bind(rb, 1, &nb);
     }
-    if (of->fp_trn)
+    where();
+    if (nsig > 0)
     {
-        close_trn(of->fp_trn);
+        isig = add_binr(rb, nsig, sig);
     }
-    if (of->fp_dhdl != NULL)
+
+    /* Global sum it all */
+    if (debug)
     {
-        gmx_fio_fclose(of->fp_dhdl);
+        fprintf(debug, "Summing %d energies\n", rb->maxreal);
     }
-    if (of->fp_field != NULL)
+    sum_bin(rb, cr);
+    where();
+
+    /* Extract all the data locally */
+
+    if (bConstrVir)
     {
-        gmx_fio_fclose(of->fp_field);
+        extract_binr(rb, isv, DIM*DIM, svir[0]);
     }
 
-    sfree(of);
-}
-
-void write_traj(FILE *fplog,t_commrec *cr,
-                gmx_mdoutf_t *of,
-                int mdof_flags,
-                gmx_mtop_t *top_global,
-                gmx_large_int_t step,double t,
-                t_state *state_local,t_state *state_global,
-                rvec *f_local,rvec *f_global,
-                int *n_xtc,rvec **x_xtc)
-{
-    int     i,j;
-    gmx_groups_t *groups;
-    rvec    *xxtc;
-    rvec *local_v;
-    rvec *global_v;
-    
-#define MX(xvf) moveit(cr,GMX_LEFT,GMX_RIGHT,#xvf,xvf)
-
-    /* MRS -- defining these variables is to manage the difference
-     * between half step and full step velocities, but there must be a better way . . . */
-
-    local_v  = state_local->v;
-    global_v = state_global->v;
-    
-    if (DOMAINDECOMP(cr))
+    /* We need the force virial and the kinetic energy for the first time through with velocity verlet */
+    if (bTemp || !bVV)
     {
-        if (mdof_flags & MDOF_CPT)
-        {
-            dd_collect_state(cr->dd,state_local,state_global);
-        }
-        else
+        if (ekind)
         {
-            if (mdof_flags & (MDOF_X | MDOF_XTC))
-            {
-                dd_collect_vec(cr->dd,state_local,state_local->x,
-                               state_global->x);
-            }
-            if (mdof_flags & MDOF_V)
+            for (j = 0; (j < inputrec->opts.ngtc); j++)
             {
-                dd_collect_vec(cr->dd,state_local,local_v,
-                               global_v);
+                if (bSumEkinhOld)
+                {
+                    extract_binr(rb, itc0[j], DIM*DIM, ekind->tcstat[j].ekinh_old[0]);
+                }
+                if (bEkinAveVel && !bReadEkin)
+                {
+                    extract_binr(rb, itc1[j], DIM*DIM, ekind->tcstat[j].ekinf[0]);
+                }
+                else if (!bReadEkin)
+                {
+                    extract_binr(rb, itc1[j], DIM*DIM, ekind->tcstat[j].ekinh[0]);
+                }
             }
-        }
-        if (mdof_flags & MDOF_F)
-        {
-            dd_collect_vec(cr->dd,state_local,f_local,f_global);
+            extract_binr(rb, idedl, 1, &(ekind->dekindl));
+            extract_binr(rb, ica, 1, &(ekind->cosacc.mvcos));
+            where();
         }
     }
-    else
+    if ((bPres || !bVV) && bFirstIterate)
     {
-        if (mdof_flags & MDOF_CPT)
-        {
-            /* All pointers in state_local are equal to state_global,
-             * but we need to copy the non-pointer entries.
-             */
-            state_global->lambda = state_local->lambda;
-            state_global->veta = state_local->veta;
-            state_global->vol0 = state_local->vol0;
-            copy_mat(state_local->box,state_global->box);
-            copy_mat(state_local->boxv,state_global->boxv);
-            copy_mat(state_local->svir_prev,state_global->svir_prev);
-            copy_mat(state_local->fvir_prev,state_global->fvir_prev);
-            copy_mat(state_local->pres_prev,state_global->pres_prev);
-        }
-        if (cr->nnodes > 1)
+        extract_binr(rb, ifv, DIM*DIM, fvir[0]);
+    }
+
+    if (bEner)
+    {
+        if (bFirstIterate)
         {
-            /* Particle decomposition, collect the data on the master node */
-            if (mdof_flags & MDOF_CPT)
+            extract_binr(rb, ie, nener, copyenerd);
+            if (rmsd_data)
             {
-                if (state_local->flags & (1<<estX))   MX(state_global->x);
-                if (state_local->flags & (1<<estV))   MX(state_global->v);
-                if (state_local->flags & (1<<estSDX)) MX(state_global->sd_X);
-                if (state_global->nrngi > 1) {
-                    if (state_local->flags & (1<<estLD_RNG)) {
-#ifdef GMX_MPI
-                        MPI_Gather(state_local->ld_rng ,
-                                   state_local->nrng*sizeof(state_local->ld_rng[0]),MPI_BYTE,
-                                   state_global->ld_rng,
-                                   state_local->nrng*sizeof(state_local->ld_rng[0]),MPI_BYTE,
-                                   MASTERRANK(cr),cr->mpi_comm_mygroup);
-#endif
-                    }
-                    if (state_local->flags & (1<<estLD_RNGI))
-                    {
-#ifdef GMX_MPI
-                        MPI_Gather(state_local->ld_rngi,
-                                   sizeof(state_local->ld_rngi[0]),MPI_BYTE,
-                                   state_global->ld_rngi,
-                                   sizeof(state_local->ld_rngi[0]),MPI_BYTE,
-                                   MASTERRANK(cr),cr->mpi_comm_mygroup);
-#endif
-                    }
-                }
+                extract_binr(rb, irmsd, inputrec->eI == eiSD2 ? 3 : 2, rmsd_data);
             }
-            else
+            if (!NEED_MUTOT(*inputrec))
             {
-                if (mdof_flags & (MDOF_X | MDOF_XTC)) MX(state_global->x);
-                if (mdof_flags & MDOF_V)              MX(global_v);
+                extract_binr(rb, imu, DIM, mu_tot);
             }
-            if (mdof_flags & MDOF_F) MX(f_global);
-         }
-     }
-
-     if (MASTER(cr))
-     {
-         if (mdof_flags & MDOF_CPT)
-         {
-             write_checkpoint(of->fn_cpt,of->bKeepAndNumCPT,
-                              fplog,cr,of->eIntegrator,of->simulation_part,
-                              of->bExpanded,of->elamstats,step,t,state_global);
-         }
-
-         if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))
-         {
-            fwrite_trn(of->fp_trn,step,t,state_local->lambda[efptFEP],
-                       state_local->box,top_global->natoms,
-                       (mdof_flags & MDOF_X) ? state_global->x : NULL,
-                       (mdof_flags & MDOF_V) ? global_v : NULL,
-                       (mdof_flags & MDOF_F) ? f_global : NULL);
-            if (gmx_fio_flush(of->fp_trn) != 0)
+
+            for (j = 0; (j < egNR); j++)
             {
-                gmx_file("Cannot write trajectory; maybe you are out of disk space?");
+                extract_binr(rb, inn[j], enerd->grpp.nener, enerd->grpp.ener[j]);
             }
-            gmx_fio_check_file_position(of->fp_trn);
-        }      
-        if (mdof_flags & MDOF_XTC) {
-            groups = &top_global->groups;
-            if (*n_xtc == -1)
+            if (inputrec->efep != efepNO)
             {
-                *n_xtc = 0;
-                for(i=0; (i<top_global->natoms); i++)
-                {
-                    if (ggrpnr(groups,egcXTC,i) == 0)
-                    {
-                        (*n_xtc)++;
-                    }
-                }
-                if (*n_xtc != top_global->natoms)
+                extract_bind(rb, idvdll, efptNR, enerd->dvdl_lin);
+                extract_bind(rb, idvdlnl, efptNR, enerd->dvdl_nonlin);
+                if (enerd->n_lambda > 0)
                 {
-                    snew(*x_xtc,*n_xtc);
+                    extract_bind(rb, iepl, enerd->n_lambda, enerd->enerpart_lambda);
                 }
             }
-            if (*n_xtc == top_global->natoms)
-            {
-                xxtc = state_global->x;
-            }
-            else
+            if (DOMAINDECOMP(cr))
             {
-                xxtc = *x_xtc;
-                j = 0;
-                for(i=0; (i<top_global->natoms); i++)
+                extract_bind(rb, inb, 1, &nb);
+                if ((int)(nb + 0.5) != cr->dd->nbonded_global)
                 {
-                    if (ggrpnr(groups,egcXTC,i) == 0)
-                    {
-                        copy_rvec(state_global->x[i],xxtc[j++]);
-                    }
+                    dd_print_missing_interactions(fplog, cr, (int)(nb + 0.5), top_global, state_local);
                 }
             }
-            if (write_xtc(of->fp_xtc,*n_xtc,step,t,
-                          state_local->box,xxtc,of->xtc_prec) == 0)
-            {
-                gmx_fatal(FARGS,"XTC error - maybe you are out of disk space?");
-            }
-            gmx_fio_check_file_position(of->fp_xtc);
+            where();
+
+            filter_enerdterm(copyenerd, FALSE, enerd->term, bTemp, bPres, bEner);
+        }
+    }
+
+    if (vcm)
+    {
+        extract_binr(rb, icm, DIM*vcm->nr, vcm->group_p[0]);
+        where();
+        extract_binr(rb, imass, vcm->nr, vcm->group_mass);
+        where();
+        if (vcm->mode == ecmANGULAR)
+        {
+            extract_binr(rb, icj, DIM*vcm->nr, vcm->group_j[0]);
+            where();
+            extract_binr(rb, icx, DIM*vcm->nr, vcm->group_x[0]);
+            where();
+            extract_binr(rb, ici, DIM*DIM*vcm->nr, vcm->group_i[0][0]);
+            where();
         }
     }
+
+    if (nsig > 0)
+    {
+        extract_binr(rb, isig, nsig, sig);
+    }
+    where();
 }
 
+int do_per_step(gmx_int64_t step, gmx_int64_t nstep)
+{
+    if (nstep != 0)
+    {
+        return ((step % nstep) == 0);
+    }
+    else
+    {
+        return 0;
+    }
+}