diff --git a/doc/src/info.txt b/doc/src/info.txt
index 0b761ccc1..43b430f74 100644
--- a/doc/src/info.txt
+++ b/doc/src/info.txt
@@ -1,111 +1,117 @@
 "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
 
 :link(lws,http://lammps.sandia.gov)
 :link(ld,Manual.html)
 :link(lc,Section_commands.html#comm)
 
 :line
 
 info command :h3
 
 [Syntax:]
 
 info args :pre
 
-args = one or more of the following keywords: {out}, {all}, {system}, {communication}, {computes}, {dumps}, {fixes}, {groups}, {regions}, {variables}, {styles}, {time}, or {configuration}
+args = one or more of the following keywords: {out}, {all}, {system}, {memory}, {communication}, {computes}, {dumps}, {fixes}, {groups}, {regions}, {variables}, {styles}, {time}, or {configuration}
      {out} values = {screen}, {log}, {append} filename, {overwrite} filename
      {styles} values = {all}, {angle}, {atom}, {bond}, {compute}, {command}, {dump}, {dihedral}, {fix}, {improper}, {integrate}, {kspace}, {minimize}, {pair}, {region} :ul
 
 [Examples:]
 
 info system
 info groups computes variables
 info all out log
 info all out append info.txt
 info styles all
 info styles atom :pre
 
 [Description:]
 
 Print out information about the current internal state of the running
 LAMMPS process. This can be helpful when debugging or validating
 complex input scripts.  Several output categories are available and
 one or more output category may be requested.
 
 The {out} flag controls where the output is sent. It can only be sent
 to one target. By default this is the screen, if it is active. The
 {log} argument selects the log file instead. With the {append} and
 {overwrite} option, followed by a filename, the output is written
 to that file, which is either appended to or overwritten, respectively.
 
 The {all} flag activates printing all categories listed below.
 
+The {configuration} category prints some information about the
+LAMMPS version as well as architecture and OS it is run on.
+
+The {memory} category prints some information about the current
+memory allocation of MPI rank 0 (this the amount of dynamically
+allocated memory reported by LAMMPS classes). Where supported,
+also some OS specific information about the size of the reserved
+memory pool size (this is where malloc() and the new operator
+request memory from) and the maximum resident set size is reported
+(this is the maximum amount of physical memory occupied so far).
+
 The {system} category prints a general system overview listing.  This
 includes the unit style, atom style, number of atoms, bonds, angles,
 dihedrals, and impropers and the number of the respective types, box
 dimensions and properties, force computing styles and more.
 
 The {communication} category prints a variety of information about
 communication and parallelization: the MPI library version level,
 the number of MPI ranks and OpenMP threads, the communication style
 and layout, the processor grid dimensions, ghost atom communication
 mode, cutoff, and related settings.
 
 The {computes} category prints a list of all currently defined
 computes, their IDs and styles and groups they operate on.
 
 The {dumps} category prints a list of all currently active dumps,
 their IDs, styles, filenames, groups, and dump frequencies.
 
 The {fixes} category prints a list of all currently defined fixes,
 their IDs and styles and groups they operate on.
 
 The {groups} category prints a list of all currently defined groups.
 
 The {regions} category prints a list of all currently defined regions,
 their IDs and styles and whether "inside" or "outside" atoms are
 selected.
 
 The {variables} category prints a list of all currently defined
 variables, their names, styles, definition and last computed value, if
 available.
 
 The {styles} category prints the list of styles available in the
 current LAMMPS binary. It supports one of the following options
 to control which category of styles is printed out:
 
 all
 angle
 atom
 bond
 compute
 command
 dump
 dihedral
 fix
 improper
 integrate
 kspace
 minimize
 pair
 region :ul
 
 The {time} category prints the accumulated CPU and wall time for the
 process that writes output (usually MPI rank 0).
 
-The {configuration} command prints some information about the LAMMPS
-version and architecture and OS it is run on. Where supported, also
-information about the memory consumption provided by the OS is
-reported.
-
 [Restrictions:] none
 
 [Related commands:]
 
 "print"_print.html
 
 [Default:]
 
 The {out} option has the default {screen}.
 
 The {styles} option has the default {all}.
diff --git a/src/info.cpp b/src/info.cpp
index e1b79951a..6a7f70b39 100644
--- a/src/info.cpp
+++ b/src/info.cpp
@@ -1,1088 +1,1107 @@
 /* ----------------------------------------------------------------------
   fputs
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    Contributing authors:  Axel Kohlmeyer (Temple U),
                           Richard Berger (Temple U)
 ------------------------------------------------------------------------- */
 
 #include <string.h>
 #include "info.h"
 #include "accelerator_kokkos.h"
 #include "atom.h"
 #include "comm.h"
 #include "compute.h"
 #include "domain.h"
 #include "dump.h"
 #include "fix.h"
 #include "force.h"
 #include "pair.h"
 #include "pair_hybrid.h"
 #include "group.h"
 #include "input.h"
 #include "modify.h"
 #include "neighbor.h"
 #include "output.h"
 #include "region.h"
 #include "universe.h"
 #include "variable.h"
 #include "update.h"
 #include "error.h"
 
 #include <time.h>
 #include <vector>
 #include <string>
 #include <algorithm>
 
 #ifdef _WIN32
 #define PSAPI_VERSION=1
 #include <windows.h>
 #include <stdint.h>
 #include <psapi.h>
 #else
 #include <sys/time.h>
 #include <sys/resource.h>
 #include <sys/utsname.h>
 #endif
 
 #if defined __linux
 #include <malloc.h>
 #endif
 
 namespace LAMMPS_NS {
 // same as in variable.cpp
 enum {INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,GETENV,
       SCALARFILE,ATOMFILE,FORMAT,EQUAL,ATOM,PYTHON};
 
 enum {COMPUTES=1<<0,
       DUMPS=1<<1,
       FIXES=1<<2,
       GROUPS=1<<3,
       REGIONS=1<<4,
       CONFIG=1<<5,
       TIME=1<<6,
-      VARIABLES=1<<7,
-      SYSTEM=1<<8,
-      COMM=1<<9,
-      ATOM_STYLES=1<<10,
-      INTEGRATE_STYLES=1<<11,
-      MINIMIZE_STYLES=1<<12,
-      PAIR_STYLES=1<<13,
-      BOND_STYLES=1<<14,
-      ANGLE_STYLES=1<<15,
-      DIHEDRAL_STYLES=1<<16,
-      IMPROPER_STYLES=1<<17,
-      KSPACE_STYLES=1<<18,
-      FIX_STYLES=1<<19,
-      COMPUTE_STYLES=1<<20,
-      REGION_STYLES=1<<21,
-      DUMP_STYLES=1<<22,
-      COMMAND_STYLES=1<<23,
+      MEMORY=1<<7,
+      VARIABLES=1<<8,
+      SYSTEM=1<<9,
+      COMM=1<<10,
+      ATOM_STYLES=1<<11,
+      INTEGRATE_STYLES=1<<12,
+      MINIMIZE_STYLES=1<<13,
+      PAIR_STYLES=1<<14,
+      BOND_STYLES=1<<15,
+      ANGLE_STYLES=1<<16,
+      DIHEDRAL_STYLES=1<<17,
+      IMPROPER_STYLES=1<<18,
+      KSPACE_STYLES=1<<19,
+      FIX_STYLES=1<<20,
+      COMPUTE_STYLES=1<<21,
+      REGION_STYLES=1<<22,
+      DUMP_STYLES=1<<23,
+      COMMAND_STYLES=1<<24,
       ALL=~0};
 
-static const int STYLES = ATOM_STYLES | INTEGRATE_STYLES | MINIMIZE_STYLES | PAIR_STYLES | BOND_STYLES | \
-                         ANGLE_STYLES | DIHEDRAL_STYLES | IMPROPER_STYLES | KSPACE_STYLES | FIX_STYLES | \
-                         COMPUTE_STYLES | REGION_STYLES | DUMP_STYLES | COMMAND_STYLES;
-
+static const int STYLES = ATOM_STYLES | INTEGRATE_STYLES | MINIMIZE_STYLES
+                        | PAIR_STYLES | BOND_STYLES | ANGLE_STYLES
+                        | DIHEDRAL_STYLES | IMPROPER_STYLES | KSPACE_STYLES
+                        | FIX_STYLES | COMPUTE_STYLES | REGION_STYLES
+                        | DUMP_STYLES | COMMAND_STYLES;
 }
 
-
 static const char *varstyles[] = {
   "index", "loop", "world", "universe", "uloop", "string", "getenv",
   "file", "atomfile", "format", "equal", "atom", "python", "(unknown)"};
 
 static const char *mapstyles[] = { "none", "array", "hash" };
 
 static const char *commstyles[] = { "brick", "tiled" };
 static const char *commlayout[] = { "uniform", "nonuniform", "irregular" };
 
 static const char bstyles[] = "pfsm";
 
 using namespace LAMMPS_NS;
 using namespace std;
 
 static void print_columns(FILE* fp, vector<string> & styles);
 
 /* ---------------------------------------------------------------------- */
 
 void Info::command(int narg, char **arg)
 {
   FILE *out=screen;
   int flags=0;
 
   if (comm->me != 0) return;
 
   // parse arguments
   int idx = 0;
 
   while (idx < narg) {
     if (strncmp(arg[idx],"all",3) == 0) {
       flags |= ALL;
       ++idx;
     } else if ((idx+1 < narg) && (strncmp(arg[idx],"out",3) == 0)
                && (strncmp(arg[idx+1],"screen",3) == 0)) {
       if ((out != screen) && (out != logfile)) fclose(out);
       out = screen;
       idx += 2;
     } else if ((idx+1 < narg) && (strncmp(arg[idx],"out",3) == 0)
                && (strncmp(arg[idx+1],"log",3) == 0)) {
       if ((out != screen) && (out != logfile)) fclose(out);
       out = logfile;
       idx += 2;
     } else if ((idx+2 < narg) && (strncmp(arg[idx],"out",3) == 0)
                && (strncmp(arg[idx+1],"append",3) == 0)) {
       if ((out != screen) && (out != logfile)) fclose(out);
       out = fopen(arg[idx+2],"a");
       idx += 3;
     } else if ((idx+2 < narg) && (strncmp(arg[idx],"out",3) == 0)
                && (strncmp(arg[idx+1],"overwrite",3) == 0)) {
       if ((out != screen) && (out != logfile)) fclose(out);
       out = fopen(arg[idx+2],"w");
       idx += 3;
     } else if (strncmp(arg[idx],"communication",5) == 0) {
       flags |= COMM;
       ++idx;
     } else if (strncmp(arg[idx],"computes",5) == 0) {
       flags |= COMPUTES;
       ++idx;
     } else if (strncmp(arg[idx],"dumps",5) == 0) {
       flags |= DUMPS;
       ++idx;
     } else if (strncmp(arg[idx],"fixes",5) == 0) {
       flags |= FIXES;
       ++idx;
     } else if (strncmp(arg[idx],"groups",3) == 0) {
       flags |= GROUPS;
       ++idx;
     } else if (strncmp(arg[idx],"regions",3) == 0) {
       flags |= REGIONS;
       ++idx;
     } else if (strncmp(arg[idx],"config",3) == 0) {
       flags |= CONFIG;
       ++idx;
     } else if (strncmp(arg[idx],"time",3) == 0) {
       flags |= TIME;
       ++idx;
+    } else if (strncmp(arg[idx],"memory",3) == 0) {
+      flags |= MEMORY;
+      ++idx;
     } else if (strncmp(arg[idx],"variables",3) == 0) {
       flags |= VARIABLES;
       ++idx;
     } else if (strncmp(arg[idx],"system",3) == 0) {
       flags |= SYSTEM;
       ++idx;
     } else if (strncmp(arg[idx],"styles",3) == 0) {
       if (idx+1 < narg) {
         ++idx;
         if (strncmp(arg[idx],"all",3) == 0) {
           flags |= STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"atom",3) == 0) {
           flags |= ATOM_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"integrate",3) == 0) {
           flags |= INTEGRATE_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"minimize",3) == 0) {
           flags |= MINIMIZE_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"pair",3) == 0) {
           flags |= PAIR_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"bond",3) == 0) {
           flags |= BOND_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"angle",3) == 0) {
           flags |= ANGLE_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"dihedral",3) == 0) {
           flags |= DIHEDRAL_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"improper",3) == 0) {
           flags |= IMPROPER_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"kspace",3) == 0) {
           flags |= KSPACE_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"fix",3) == 0) {
           flags |= FIX_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"compute",4) == 0) {
           flags |= COMPUTE_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"region",3) == 0) {
           flags |= REGION_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"dump",3) == 0) {
           flags |= DUMP_STYLES;
           ++idx;
         } else if (strncmp(arg[idx],"command",4) == 0) {
           flags |= COMMAND_STYLES;
           ++idx;
         } else {
           flags |= STYLES;
         }
       } else {
         flags |= STYLES;
         ++idx;
       }
     } else {
       error->warning(FLERR,"Ignoring unknown or incorrect info command flag");
       ++idx;
     }
   }
 
   if (out == NULL) return;
 
   fputs("\nInfo-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info\n",out);
   time_t now = time(NULL);
   fprintf(out,"Printed on %s\n",ctime(&now));
 
   if (flags & CONFIG) {
 
     fprintf(out,"\nLAMMPS version: %s / %s\n",
             universe->version, universe->num_ver);
     fprintf(out,"sizeof(smallint): %3d-bit\n",(int)sizeof(smallint)*8);
     fprintf(out,"sizeof(imageint): %3d-bit\n",(int)sizeof(imageint)*8);
     fprintf(out,"sizeof(tagint):   %3d-bit\n",(int)sizeof(tagint)*8);
     fprintf(out,"sizeof(bigint):   %3d-bit\n",(int)sizeof(bigint)*8);
 
 #if defined(_WIN32)
     DWORD fullversion,majorv,minorv,buildv=0;
 
     fullversion = GetVersion();
     majorv = (DWORD) (LOBYTE(LOWORD(fullversion)));
     minorv = (DWORD) (HIBYTE(LOWORD(fullversion)));
     if (fullversion < 0x80000000)
       buildv = (DWORD) (HIWORD(fullversion));
 
     SYSTEM_INFO si;
     GetSystemInfo(&si);
 
     const char *machine;
     switch (si.wProcessorArchitecture) {
     case PROCESSOR_ARCHITECTURE_AMD64:
       machine = (const char *) "x86_64";
       break;
     case PROCESSOR_ARCHITECTURE_ARM:
       machine = (const char *) "arm";
       break;
     case PROCESSOR_ARCHITECTURE_IA64:
       machine = (const char *) "ia64";
       break;
     case PROCESSOR_ARCHITECTURE_INTEL:
       machine = (const char *) "i386";
       break;
     default:
       machine = (const char *) "(unknown)";
     }
     fprintf(out,"\nOS information: Windows %d.%d (%d) on %s\n",
             majorv,minorv,buildv,machine);
 #else
     struct utsname ut;
     uname(&ut);
     fprintf(out,"\nOS information: %s %s on %s\n",
             ut.sysname, ut.release, ut.machine);
 #endif
-
-    fprintf(out,"\nMemory allocation information (MPI rank 0)\n");
+  }
+  
+  if (flags & MEMORY) {
+
+    fprintf(out,"\nMemory allocation information (MPI rank 0):\n\n");
+
+    bigint bytes = 0;
+    bytes += atom->memory_usage();
+    bytes += neighbor->memory_usage();
+    bytes += comm->memory_usage();
+    bytes += update->memory_usage();
+    bytes += force->memory_usage();
+    bytes += modify->memory_usage();
+    for (int i = 0; i < output->ndump; i++)
+      bytes += output->dump[i]->memory_usage();
+    double mbytes = bytes/1024.0/1024.0;
+    fprintf(out,"Total dynamically allocated memory: %.4g Mbyte\n",mbytes);
 
 #if defined(_WIN32)
     HANDLE phandle = GetCurrentProcess();
     PROCESS_MEMORY_COUNTERS_EX pmc;
     GetProcessMemoryInfo(phandle,(PROCESS_MEMORY_COUNTERS *)&pmc,sizeof(pmc));
-    fprintf(out,"Non-shared memory use: %.3g Mbyte\n",
+    fprintf(out,"Non-shared memory use: %.4g Mbyte\n",
             (double)pmc.PrivateUsage/1048576.0);
-    fprintf(out,"Maximum working set size: %.3g Mbyte\n",
+    fprintf(out,"Maximum working set size: %.4g Mbyte\n",
             (double)pmc.PeakWorkingSetSize/1048576.0);
 #else
 #if defined(__linux)
     struct mallinfo mi;
     mi = mallinfo();
-    fprintf(out,"Total dynamically allocated memory: %.3g Mbyte\n",
-            (double)mi.uordblks/1048576.0);
+    fprintf(out,"Current reserved memory pool size: %.4g Mbyte\n",
+            (double)mi.uordblks/1048576.0+(double)mi.hblkhd/1048576.0);
 #endif
     struct rusage ru;
     if (getrusage(RUSAGE_SELF, &ru) == 0) {
-      fprintf(out,"Maximum resident set size: %.3g Mbyte\n",
+      fprintf(out,"Maximum resident set size: %.4g Mbyte\n",
               (double)ru.ru_maxrss/1024.0);
     }
 #endif
   }
 
   if (flags & COMM) {
     int major,minor;
     MPI_Get_version(&major,&minor);
 
     fprintf(out,"\nCommunication information:\n");
     fprintf(out,"MPI library level: MPI v%d.%d\n",major,minor);
     fprintf(out,"Comm style = %s,  Comm layout = %s\n",
             commstyles[comm->style], commlayout[comm->layout]);
     fprintf(out,"Communicate velocities for ghost atoms = %s\n",
             comm->ghost_velocity ? "yes" : "no");
 
     if (comm->mode == 0) {
       fprintf(out,"Communication mode = single\n");
       fprintf(out,"Communication cutoff = %g\n",
               MAX(comm->cutghostuser,neighbor->cutneighmax));
     }
 
     if (comm->mode == 1) {
       fprintf(out,"Communication mode = multi\n");
       double cut;
       for (int i=1; i <= atom->ntypes && neighbor->cuttype; ++i) {
         cut = neighbor->cuttype[i];
         if (comm->cutusermulti) cut = MAX(cut,comm->cutusermulti[i]);
         fprintf(out,"Communication cutoff for type %d = %g\n", i, cut);
       }
     }
     fprintf(out,"Nprocs = %d,   Nthreads = %d\n",
             comm->nprocs, comm->nthreads);
     if (domain->box_exist)
       fprintf(out,"Processor grid = %d x %d x %d\n",comm->procgrid[0],
             comm->procgrid[1], comm->procgrid[2]);
   }
 
   if (flags & SYSTEM) {
     fprintf(out,"\nSystem information:\n");
     fprintf(out,"Units      = %s\n",update->unit_style);
     fprintf(out,"Atom style = %s\n", atom->atom_style);
     fprintf(out,"Atom map   = %s\n", mapstyles[atom->map_style]);
     if (atom->molecular > 0) {
       const char *msg;
       msg = (atom->molecular == 2) ? "template" : "standard";
       fprintf(out,"Molecule type = %s\n",msg);
     }
     fprintf(out,"Atoms = " BIGINT_FORMAT ",  types = %d,  style = %s\n",
             atom->natoms, atom->ntypes, force->pair_style);
 
     if (force->pair && strstr(force->pair_style,"hybrid")) {
       PairHybrid *hybrid = (PairHybrid *)force->pair;
       fprintf(out,"Hybrid sub-styles:");
       for (int i=0; i < hybrid->nstyles; ++i)
         fprintf(out," %s", hybrid->keywords[i]);
       fputc('\n',out);
     }
     if (atom->molecular > 0) {
       const char *msg;
       msg = force->bond_style ? force->bond_style : "none";
       fprintf(out,"Bonds = " BIGINT_FORMAT ",  types = %d,  style = %s\n",
               atom->nbonds, atom->nbondtypes, msg);
 
       msg = force->angle_style ? force->angle_style : "none";
       fprintf(out,"Angles = " BIGINT_FORMAT ",  types = %d,  style = %s\n",
               atom->nangles, atom->nangletypes, msg);
 
       msg = force->dihedral_style ? force->dihedral_style : "none";
       fprintf(out,"Dihedrals = " BIGINT_FORMAT ",  types = %d,  style = %s\n",
               atom->ndihedrals, atom->ndihedraltypes, msg);
 
       msg = force->improper_style ? force->improper_style : "none";
       fprintf(out,"Impropers = " BIGINT_FORMAT ",  types = %d,  style = %s\n",
               atom->nimpropers, atom->nimpropertypes, msg);
 
       const double * const special_lj   = force->special_lj;
       const double * const special_coul = force->special_coul;
 
       fprintf(out,"Special bond factors lj =   %-10g %-10g %-10g\n"
               "Special bond factors coul = %-10g %-10g %-10g\n",
               special_lj[1],special_lj[2],special_lj[3],
               special_coul[1],special_coul[2],special_coul[3]);
     }
 
     fprintf(out,"Kspace style = %s\n",
             force->kspace ? force->kspace_style : "none");
 
     if (domain->box_exist) {
       fprintf(out,"\nDimensions = %d\n",domain->dimension);
       fprintf(out,"%s box = %g x %g x %g\n",
               domain->triclinic ? "Triclinic" : "Orthogonal",
               domain->xprd, domain->yprd, domain->zprd);
       fprintf(out,"Boundaries = %c,%c %c,%c %c,%c\n",
               bstyles[domain->boundary[0][0]],bstyles[domain->boundary[0][1]],
               bstyles[domain->boundary[1][0]],bstyles[domain->boundary[1][1]],
               bstyles[domain->boundary[2][0]],bstyles[domain->boundary[2][1]]);
       fprintf(out,"xlo, xhi = %g, %g\n", domain->boxlo[0], domain->boxhi[0]);
       fprintf(out,"ylo, yhi = %g, %g\n", domain->boxlo[1], domain->boxhi[1]);
       fprintf(out,"zlo, zhi = %g, %g\n", domain->boxlo[2], domain->boxhi[2]);
       if (domain->triclinic)
           fprintf(out,"Xy, xz, yz = %g, %g, %g\n",
                   domain->xy, domain->xz, domain->yz);
     } else {
       fputs("\nBox has not yet been created\n",out);
     }
   }
 
   if (flags & GROUPS) {
     int ngroup = group->ngroup;
     char **names = group->names;
     int *dynamic = group->dynamic;
     fprintf(out,"\nGroup information:\n");
     for (int i=0; i < ngroup; ++i) {
       if (names[i])
         fprintf(out,"Group[%2d]: %s (%s)\n",
                 i, names[i], dynamic[i] ? "dynamic" : "static");
     }
   }
 
   if (flags & REGIONS) {
     int nreg = domain->nregion;
     Region **regs = domain->regions;
     fprintf(out,"\nRegion information:\n");
     for (int i=0; i < nreg; ++i) {
       fprintf(out,"Region[%3d]: %s,  style = %s,  side = %s\n",
               i, regs[i]->id, regs[i]->style,
               regs[i]->interior ? "in" : "out");
     }
   }
 
   if (flags & COMPUTES) {
     int ncompute = modify->ncompute;
     Compute **compute = modify->compute;
     char **names = group->names;
     fprintf(out,"\nCompute information:\n");
     for (int i=0; i < ncompute; ++i) {
       fprintf(out,"Compute[%3d]: %s,  style = %s,  group = %s\n",
               i, compute[i]->id, compute[i]->style,
               names[compute[i]->igroup]);
     }
   }
 
   if (flags & DUMPS) {
     int ndump = output->ndump;
     Dump **dump = output->dump;
     int *nevery = output->every_dump;           \
     char **vnames = output->var_dump;
     char **names = group->names;
     fprintf(out,"\nDump information:\n");
     for (int i=0; i < ndump; ++i) {
       fprintf(out,"Dump[%3d]: %s,  file = %s,  style = %s,  group = %s,  ",
               i, dump[i]->id, dump[i]->filename,
               dump[i]->style, names[dump[i]->igroup]);
       if (nevery[i]) {
         fprintf(out,"every = %d\n", nevery[i]);
       } else {
         fprintf(out,"every = %s\n", vnames[i]);
       }
     }
   }
 
   if (flags & FIXES) {
     int nfix = modify->nfix;
     Fix **fix = modify->fix;
     char **names = group->names;
     fprintf(out,"\nFix information:\n");
     for (int i=0; i < nfix; ++i) {
       fprintf(out,"Fix[%3d]: %s,  style = %s,  group = %s\n",
               i, fix[i]->id, fix[i]->style, names[fix[i]->igroup]);
     }
   }
 
   if (flags & VARIABLES) {
     int nvar = input->variable->nvar;
     int *style = input->variable->style;
     char **names = input->variable->names;
     char ***data = input->variable->data;
     fprintf(out,"\nVariable information:\n");
     for (int i=0; i < nvar; ++i) {
       int ndata = 1;
       fprintf(out,"Variable[%3d]: %-10s,  style = %-10s,  def =",
               i,names[i],varstyles[style[i]]);
       if ((style[i] != LOOP) && (style[i] != ULOOP))
         ndata = input->variable->num[i];
       for (int j=0; j < ndata; ++j)
         fprintf(out," %s",data[i][j]);
       fputs("\n",out);
     }
   }
 
   if (flags & TIME) {
     double wallclock = MPI_Wtime() - lmp->initclock;
     double cpuclock = 0.0;
 
 #if defined(_WIN32)
     // from MSD docs.
     FILETIME ct,et,kt,ut;
     union { FILETIME ft; uint64_t ui; } cpu;
     if (GetProcessTimes(GetCurrentProcess(),&ct,&et,&kt,&ut)) {
       cpu.ft = ut;
       cpuclock = cpu.ui * 0.0000001;
     }
 #else /* POSIX */
     struct rusage ru;
     if (getrusage(RUSAGE_SELF, &ru) == 0) {
       cpuclock  = (double) ru.ru_utime.tv_sec;
       cpuclock += (double) ru.ru_utime.tv_usec * 0.000001;
     }
 #endif /* ! _WIN32 */
 
     int cpuh,cpum,cpus,wallh,wallm,walls;
     cpus = fmod(cpuclock,60.0);
     cpuclock = (cpuclock - cpus) / 60.0;
     cpum = fmod(cpuclock,60.0);
     cpuh = (cpuclock - cpum) / 60.0;
     walls = fmod(wallclock,60.0);
     wallclock = (wallclock - walls) / 60.0;
     wallm = fmod(wallclock,60.0);
     wallh = (wallclock - wallm) / 60.0;
     fprintf(out,"\nTotal time information (MPI rank 0):\n"
             "  CPU time: %4d:%02d:%02d\n"
             " Wall time: %4d:%02d:%02d\n",
             cpuh,cpum,cpus,wallh,wallm,walls);
   }
 
   if (flags & STYLES) {
     available_styles(out, flags);
   }
 
   fputs("\nInfo-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info\n\n",out);
 
   // close output file pointer if opened locally thus forcing a hard sync.
   if ((out != screen) && (out != logfile))
     fclose(out);
 }
 
 
 void Info::available_styles(FILE * out, int flags)
 {
 
   fprintf(out,"\nStyles information:\n");
 
   if(flags & ATOM_STYLES)      atom_styles(out);
   if(flags & INTEGRATE_STYLES) integrate_styles(out);
   if(flags & MINIMIZE_STYLES)  minimize_styles(out);
   if(flags & PAIR_STYLES)      pair_styles(out);
   if(flags & BOND_STYLES)      bond_styles(out);
   if(flags & ANGLE_STYLES)     angle_styles(out);
   if(flags & DIHEDRAL_STYLES)  dihedral_styles(out);
   if(flags & IMPROPER_STYLES)  improper_styles(out);
   if(flags & KSPACE_STYLES)    kspace_styles(out);
   if(flags & FIX_STYLES)       fix_styles(out);
   if(flags & COMPUTE_STYLES)   compute_styles(out);
   if(flags & REGION_STYLES)    region_styles(out);
   if(flags & DUMP_STYLES)      dump_styles(out);
   if(flags & COMMAND_STYLES)   command_styles(out);
 }
 
 void Info::atom_styles(FILE * out)
 {
   fprintf(out, "\nAtom styles:\n");
 
   vector<string> styles;
 
   for(Atom::AtomVecCreatorMap::iterator it = atom->avec_map->begin(); it != atom->avec_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::integrate_styles(FILE * out)
 {
   fprintf(out, "\nIntegrate styles:\n");
 
   vector<string> styles;
 
   for(Update::IntegrateCreatorMap::iterator it = update->integrate_map->begin(); it != update->integrate_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::minimize_styles(FILE * out)
 {
   fprintf(out, "\nMinimize styles:\n");
 
   vector<string> styles;
 
   for(Update::MinimizeCreatorMap::iterator it = update->minimize_map->begin(); it != update->minimize_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::pair_styles(FILE * out)
 {
   fprintf(out, "\nPair styles:\n");
 
   vector<string> styles;
 
   for(Force::PairCreatorMap::iterator it = force->pair_map->begin(); it != force->pair_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::bond_styles(FILE * out)
 {
   fprintf(out, "\nBond styles:\n");
 
   vector<string> styles;
 
   for(Force::BondCreatorMap::iterator it = force->bond_map->begin(); it != force->bond_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::angle_styles(FILE * out)
 {
   fprintf(out, "\nAngle styles:\n");
 
   vector<string> styles;
 
   for(Force::AngleCreatorMap::iterator it = force->angle_map->begin(); it != force->angle_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::dihedral_styles(FILE * out)
 {
   fprintf(out, "\nDihedral styles:\n");
 
   vector<string> styles;
 
   for(Force::DihedralCreatorMap::iterator it = force->dihedral_map->begin(); it != force->dihedral_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::improper_styles(FILE * out)
 {
   fprintf(out, "\nImproper styles:\n");
 
   vector<string> styles;
 
   for(Force::ImproperCreatorMap::iterator it = force->improper_map->begin(); it != force->improper_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::kspace_styles(FILE * out)
 {
   fprintf(out, "\nKSpace styles:\n");
 
   vector<string> styles;
 
   for(Force::KSpaceCreatorMap::iterator it = force->kspace_map->begin(); it != force->kspace_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::fix_styles(FILE * out)
 {
   fprintf(out, "\nFix styles:\n");
 
   vector<string> styles;
 
   for(Modify::FixCreatorMap::iterator it = modify->fix_map->begin(); it != modify->fix_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::compute_styles(FILE * out)
 {
   fprintf(out, "\nCompute styles:\n");
 
   vector<string> styles;
 
   for(Modify::ComputeCreatorMap::iterator it = modify->compute_map->begin(); it != modify->compute_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::region_styles(FILE * out)
 {
   fprintf(out, "\nRegion styles:\n");
 
   vector<string> styles;
 
   for(Domain::RegionCreatorMap::iterator it = domain->region_map->begin(); it != domain->region_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::dump_styles(FILE * out)
 {
   fprintf(out, "\nDump styles:\n");
 
   vector<string> styles;
 
   for(Output::DumpCreatorMap::iterator it = output->dump_map->begin(); it != output->dump_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 void Info::command_styles(FILE * out)
 {
   fprintf(out, "\nCommand styles (add-on input script commands):\n");
 
   vector<string> styles;
 
   for(Input::CommandCreatorMap::iterator it = input->command_map->begin(); it != input->command_map->end(); ++it) {
     styles.push_back(it->first);
   }
 
   print_columns(out, styles);
   fprintf(out, "\n\n\n");
 }
 
 
 /* ---------------------------------------------------------------------- */
 
 // the is_active() function returns true if the selected style or name
 // in the selected category is currently in use.
 
 bool Info::is_active(const char *category, const char *name)
 {
   if ((category == NULL) || (name == NULL)) return false;
   const char *style = "none";
   const int len = strlen(name);
 
   if (strcmp(category,"package") == 0) {
     if (strcmp(name,"gpu") == 0) {
       return (modify->find_fix("package_gpu") >= 0) ? true : false;
     } else if (strcmp(name,"intel") == 0) {
       return (modify->find_fix("package_intel") >= 0) ? true : false;
     } else if (strcmp(name,"kokkos") == 0) {
       return (lmp->kokkos && lmp->kokkos->kokkos_exists) ? true : false;
     } else if (strcmp(name,"omp") == 0) {
       return (modify->find_fix("package_omp") >= 0) ? true : false;
     } else error->all(FLERR,"Unknown name for info package category");
 
   } else if (strcmp(category,"newton") == 0) {
     if (strcmp(name,"pair") == 0) return (force->newton_pair != 0);
     else if (strcmp(name,"bond") == 0) return (force->newton_bond != 0);
     else if (strcmp(name,"any") == 0) return (force->newton != 0);
     else error->all(FLERR,"Unknown name for info newton category");
 
   } else if (strcmp(category,"pair") == 0) {
     if (force->pair == NULL) return false;
     if (strcmp(name,"single") == 0) return (force->pair->single_enable != 0);
     else if (strcmp(name,"respa") == 0) return (force->pair->respa_enable != 0);
     else if (strcmp(name,"manybody") == 0) return (force->pair->manybody_flag != 0);
     else if (strcmp(name,"tail") == 0) return (force->pair->tail_flag != 0);
     else if (strcmp(name,"shift") == 0) return (force->pair->offset_flag != 0);
     else error->all(FLERR,"Unknown name for info pair category");
 
   } else if (strcmp(category,"comm_style") == 0) {
     style = commstyles[comm->style];
   } else if (strcmp(category,"min_style") == 0) {
     style = update->minimize_style;
   } else if (strcmp(category,"run_style") == 0) {
     style = update->integrate_style;
   } else if (strcmp(category,"atom_style") == 0) {
     style = atom->atom_style;
   } else if (strcmp(category,"pair_style") == 0) {
     style = force->pair_style;
   } else if (strcmp(category,"bond_style") == 0) {
     style = force->bond_style;
   } else if (strcmp(category,"angle_style") == 0) {
     style = force->angle_style;
   } else if (strcmp(category,"dihedral_style") == 0) {
     style = force->dihedral_style;
   } else if (strcmp(category,"improper_style") == 0) {
     style = force->improper_style;
   } else if (strcmp(category,"kspace_style") == 0) {
     style = force->kspace_style;
   } else error->all(FLERR,"Unknown category for info is_active()");
 
   int match = 0;
   if (strcmp(style,name) == 0) match = 1;
 
   if (!match && lmp->suffix_enable) {
     if (lmp->suffix) {
       char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix)];
       sprintf(name_w_suffix,"%s/%s",name,lmp->suffix);
       if (strcmp(style,name_w_suffix) == 0) match = 1;
       delete[] name_w_suffix;
     }
     if (!match && lmp->suffix2) {
       char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix2)];
       sprintf(name_w_suffix,"%s/%s",name,lmp->suffix2);
       if (strcmp(style,name_w_suffix) == 0) match = 1;
       delete[] name_w_suffix;
     }
   }
   return match ? true : false;
 }
 
 /* ---------------------------------------------------------------------- */
 
 // the is_available() function returns true if the selected style
 // or name in the selected category is available for use (but need
 // not be currently active).
 
 bool Info::is_available(const char *category, const char *name)
 {
   if ((category == NULL) || (name == NULL)) return false;
   const int len = strlen(name);
   int match = 0;
 
   if (strcmp(category,"command") == 0) {
     if (input->command_map->find(name) != input->command_map->end())
       match = 1;
 
   } else if (strcmp(category,"compute") == 0) {
     if (modify->compute_map->find(name) != modify->compute_map->end())
       match = 1;
 
     if (!match && lmp->suffix_enable) {
       if (lmp->suffix) {
         char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix)];
         sprintf(name_w_suffix,"%s/%s",name,lmp->suffix);
         if (modify->compute_map->find(name_w_suffix) != modify->compute_map->end())
           match = 1;
         delete[] name_w_suffix;
       }
       if (!match && lmp->suffix2) {
         char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix2)];
         sprintf(name_w_suffix,"%s/%s",name,lmp->suffix2);
         if (modify->compute_map->find(name_w_suffix) != modify->compute_map->end())
           match = 1;
         delete[] name_w_suffix;
       }
     }
   } else if (strcmp(category,"fix") == 0) {
     if (modify->fix_map->find(name) != modify->fix_map->end())
       match = 1;
 
     if (!match && lmp->suffix_enable) {
       if (lmp->suffix) {
         char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix)];
         sprintf(name_w_suffix,"%s/%s",name,lmp->suffix);
         if (modify->fix_map->find(name_w_suffix) != modify->fix_map->end())
           match = 1;
         delete[] name_w_suffix;
       }
       if (!match && lmp->suffix2) {
         char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix2)];
         sprintf(name_w_suffix,"%s/%s",name,lmp->suffix2);
         if (modify->fix_map->find(name_w_suffix) != modify->fix_map->end())
           match = 1;
         delete[] name_w_suffix;
       }
     }
   } else if (strcmp(category,"pair_style") == 0) {
     if (force->pair_map->find(name) != force->pair_map->end())
       match = 1;
 
     if (!match && lmp->suffix_enable) {
       if (lmp->suffix) {
         char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix)];
         sprintf(name_w_suffix,"%s/%s",name,lmp->suffix);
         if (force->pair_map->find(name_w_suffix) != force->pair_map->end())
           match = 1;
         delete[] name_w_suffix;
       }
       if (!match && lmp->suffix2) {
         char *name_w_suffix = new char [len + 2 + strlen(lmp->suffix2)];
         sprintf(name_w_suffix,"%s/%s",name,lmp->suffix2);
         if (force->pair_map->find(name_w_suffix) != force->pair_map->end())
           match = 1;
         delete[] name_w_suffix;
       }
     }
   } else if (strcmp(category,"feature") == 0) {
     if (strcmp(name,"gzip") == 0) {
       return has_gzip_support();
     } else if (strcmp(name,"png") == 0) {
       return has_png_support();
     } else if (strcmp(name,"jpeg") == 0) {
       return has_jpeg_support();
     } else if (strcmp(name,"ffmpeg") == 0) {
       return has_ffmpeg_support();
     } else if (strcmp(name,"exceptions") == 0) {
       return has_exceptions();
     }
   } else error->all(FLERR,"Unknown category for info is_available()");
 
   return match ? true : false;
 }
 
 /* ---------------------------------------------------------------------- */
 
 // the is_defined() function returns true if a particular ID of the
 // selected category (e.g. fix ID, group ID, region ID etc.) has been
 // defined and thus can be accessed. It does *NOT* check whether a
 // particular ID has a particular style.
 
 bool Info::is_defined(const char *category, const char *name)
 {
   if ((category == NULL) || (name == NULL)) return false;
 
   if (strcmp(category,"compute") == 0) {
     int ncompute = modify->ncompute;
     Compute **compute = modify->compute;
     for (int i=0; i < ncompute; ++i) {
       if (strcmp(compute[i]->id,name) == 0)
         return true;
     }
   } else if (strcmp(category,"dump") == 0) {
     int ndump = output->ndump;
     Dump **dump = output->dump;
     for (int i=0; i < ndump; ++i) {
       if (strcmp(dump[i]->id,name) == 0)
         return true;
     }
   } else if (strcmp(category,"fix") == 0) {
     int nfix = modify->nfix;
     Fix **fix = modify->fix;
     for (int i=0; i < nfix; ++i) {
       if (strcmp(fix[i]->id,name) == 0)
         return true;
     }
   } else if (strcmp(category,"group") == 0) {
     int ngroup = group->ngroup;
     char **names = group->names;
     for (int i=0; i < ngroup; ++i) {
       if (strcmp(names[i],name) == 0)
         return true;
     }
   } else if (strcmp(category,"region") == 0) {
     int nreg = domain->nregion;
     Region **regs = domain->regions;
     for (int i=0; i < nreg; ++i) {
       if (strcmp(regs[i]->id,name) == 0)
         return true;
     }
   } else if (strcmp(category,"variable") == 0) {
     int nvar = input->variable->nvar;
     char **names = input->variable->names;
     for (int i=0; i < nvar; ++i) {
       if (strcmp(names[i],name) == 0)
         return true;
     }
   } else error->all(FLERR,"Unknown category for info is_defined()");
 
   return false;
 }
 
 static void print_columns(FILE* fp, vector<string> & styles)
 {
   if (styles.size() == 0) {
     fprintf(fp, "\nNone");
     return;
   }
 
   std::sort(styles.begin(), styles.end());
 
   int pos = 80;
   for (int i = 0; i < styles.size(); ++i) {
 
     // skip "secret" styles
     if (isupper(styles[i][0])) continue;
 
     int len = styles[i].length();
     if (pos + len > 80) {
       fprintf(fp,"\n");
       pos = 0;
     }
 
     if (len < 16) {
       fprintf(fp,"%-16s",styles[i].c_str());
       pos += 16;
     } else if (len < 32) {
       fprintf(fp,"%-32s",styles[i].c_str());
       pos += 32;
     } else if (len < 48) {
       fprintf(fp,"%-48s",styles[i].c_str());
       pos += 48;
     } else if (len < 64) {
       fprintf(fp,"%-64s",styles[i].c_str());
       pos += 64;
     } else {
       fprintf(fp,"%-80s",styles[i].c_str());
       pos += 80;
     }
   }
 }
 
 bool Info::has_gzip_support() const {
 #ifdef LAMMPS_GZIP
   return true;
 #else
   return false;
 #endif
 }
 
 bool Info::has_png_support() const {
 #ifdef LAMMPS_PNG
   return true;
 #else
   return false;
 #endif
 }
 
 bool Info::has_jpeg_support() const {
 #ifdef LAMMPS_JPEG
   return true;
 #else
   return false;
 #endif
 }
 
 bool Info::has_ffmpeg_support() const {
 #ifdef LAMMPS_FFMPEG
   return true;
 #else
   return false;
 #endif
 }
 
 bool Info::has_exceptions() const {
 #ifdef LAMMPS_EXCEPTIONS
   return true;
 #else
   return false;
 #endif
 }
 
 /* ---------------------------------------------------------------------- */
 
 char **Info::get_variable_names(int &num) {
     num = input->variable->nvar;
     return input->variable->names;
 }
diff --git a/src/output.cpp b/src/output.cpp
index 89e5c8d9c..a2275b74b 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -1,835 +1,842 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include "output.h"
 #include "style_dump.h"
 #include "atom.h"
 #include "neighbor.h"
 #include "input.h"
 #include "variable.h"
 #include "comm.h"
 #include "update.h"
 #include "group.h"
 #include "domain.h"
 #include "thermo.h"
 #include "modify.h"
 #include "compute.h"
 #include "force.h"
 #include "dump.h"
 #include "write_restart.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 #define DELTA 1
 
 /* ----------------------------------------------------------------------
    initialize all output
 ------------------------------------------------------------------------- */
 
 Output::Output(LAMMPS *lmp) : Pointers(lmp)
 {
   // create default computes for temp,pressure,pe
 
   char **newarg = new char*[4];
   newarg[0] = (char *) "thermo_temp";
   newarg[1] = (char *) "all";
   newarg[2] = (char *) "temp";
   modify->add_compute(3,newarg,1);
 
   newarg[0] = (char *) "thermo_press";
   newarg[1] = (char *) "all";
   newarg[2] = (char *) "pressure";
   newarg[3] = (char *) "thermo_temp";
   modify->add_compute(4,newarg,1);
 
   newarg[0] = (char *) "thermo_pe";
   newarg[1] = (char *) "all";
   newarg[2] = (char *) "pe";
   modify->add_compute(3,newarg,1);
 
   delete [] newarg;
 
   // create default Thermo class
 
   newarg = new char*[1];
   newarg[0] = (char *) "one";
   thermo = new Thermo(lmp,1,newarg);
   delete [] newarg;
 
   thermo_every = 0;
   var_thermo = NULL;
 
   ndump = 0;
   max_dump = 0;
   every_dump = NULL;
   next_dump = NULL;
   last_dump = NULL;
   var_dump = NULL;
   ivar_dump = NULL;
   dump = NULL;
 
   restart_flag = restart_flag_single = restart_flag_double = 0;
   restart_every_single = restart_every_double = 0;
   last_restart = -1;
   restart1 = restart2a = restart2b = NULL;
   var_restart_single = var_restart_double = NULL;
   restart = NULL;
 
   dump_map = new DumpCreatorMap();
 
 #define DUMP_CLASS
 #define DumpStyle(key,Class) \
   (*dump_map)[#key] = &dump_creator<Class>;
 #include "style_dump.h"
 #undef DumpStyle
 #undef DUMP_CLASS
 }
 
 /* ----------------------------------------------------------------------
    free all memory
 ------------------------------------------------------------------------- */
 
 Output::~Output()
 {
   if (thermo) delete thermo;
   delete [] var_thermo;
 
   memory->destroy(every_dump);
   memory->destroy(next_dump);
   memory->destroy(last_dump);
   for (int i = 0; i < ndump; i++) delete [] var_dump[i];
   memory->sfree(var_dump);
   memory->destroy(ivar_dump);
   for (int i = 0; i < ndump; i++) delete dump[i];
   memory->sfree(dump);
 
   delete [] restart1;
   delete [] restart2a;
   delete [] restart2b;
   delete [] var_restart_single;
   delete [] var_restart_double;
   delete restart;
 
   delete dump_map;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void Output::init()
 {
   thermo->init();
   if (var_thermo) {
     ivar_thermo = input->variable->find(var_thermo);
     if (ivar_thermo < 0)
       error->all(FLERR,"Variable name for thermo every does not exist");
     if (!input->variable->equalstyle(ivar_thermo))
       error->all(FLERR,"Variable for thermo every is invalid style");
   }
 
   for (int i = 0; i < ndump; i++) dump[i]->init();
   for (int i = 0; i < ndump; i++)
     if (every_dump[i] == 0) {
       ivar_dump[i] = input->variable->find(var_dump[i]);
       if (ivar_dump[i] < 0)
         error->all(FLERR,"Variable name for dump every does not exist");
       if (!input->variable->equalstyle(ivar_dump[i]))
         error->all(FLERR,"Variable for dump every is invalid style");
     }
 
   if (restart_flag_single && restart_every_single == 0) {
     ivar_restart_single = input->variable->find(var_restart_single);
     if (ivar_restart_single < 0)
       error->all(FLERR,"Variable name for restart does not exist");
     if (!input->variable->equalstyle(ivar_restart_single))
       error->all(FLERR,"Variable for restart is invalid style");
   }
   if (restart_flag_double && restart_every_double == 0) {
     ivar_restart_double = input->variable->find(var_restart_double);
     if (ivar_restart_double < 0)
       error->all(FLERR,"Variable name for restart does not exist");
     if (!input->variable->equalstyle(ivar_restart_double))
       error->all(FLERR,"Variable for restart is invalid style");
   }
 }
 
 /* ----------------------------------------------------------------------
    perform output for setup of run/min
    do dump first, so memory_usage will include dump allocation
    do thermo last, so will print after memory_usage
    memflag = 0/1 for printing out memory usage
 ------------------------------------------------------------------------- */
 
 void Output::setup(int memflag)
 {
   bigint ntimestep = update->ntimestep;
 
   // perform dump at start of run only if:
   //   current timestep is multiple of every and last dump not >= this step
   //   this is first run after dump created and firstflag is set
   //   note that variable freq will not write unless triggered by firstflag
   // set next_dump to multiple of every or variable value
   // set next_dump_any to smallest next_dump
   // wrap dumps that invoke computes and variable eval with clear/add
   // if dump not written now, use addstep_compute_all() since don't know
   //   what computes the dump write would invoke
   // if no dumps, set next_dump_any to last+1 so will not influence next
 
   int writeflag;
 
   if (ndump && update->restrict_output == 0) {
     for (int idump = 0; idump < ndump; idump++) {
       if (dump[idump]->clearstep || every_dump[idump] == 0)
         modify->clearstep_compute();
       writeflag = 0;
       if (every_dump[idump] && ntimestep % every_dump[idump] == 0 &&
           last_dump[idump] != ntimestep) writeflag = 1;
       if (last_dump[idump] < 0 && dump[idump]->first_flag == 1) writeflag = 1;
 
       if (writeflag) {
         dump[idump]->write();
         last_dump[idump] = ntimestep;
       }
       if (every_dump[idump])
         next_dump[idump] =
           (ntimestep/every_dump[idump])*every_dump[idump] + every_dump[idump];
       else {
         bigint nextdump = static_cast<bigint>
           (input->variable->compute_equal(ivar_dump[idump]));
         if (nextdump <= ntimestep)
           error->all(FLERR,"Dump every variable returned a bad timestep");
         next_dump[idump] = nextdump;
       }
       if (dump[idump]->clearstep || every_dump[idump] == 0) {
         if (writeflag) modify->addstep_compute(next_dump[idump]);
         else modify->addstep_compute_all(next_dump[idump]);
       }
       if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]);
       else next_dump_any = next_dump[0];
     }
   } else next_dump_any = update->laststep + 1;
 
   // do not write restart files at start of run
   // set next_restart values to multiple of every or variable value
   // wrap variable eval with clear/add
   // if no restarts, set next_restart to last+1 so will not influence next
 
   if (restart_flag && update->restrict_output == 0) {
     if (restart_flag_single) {
       if (restart_every_single)
         next_restart_single =
           (ntimestep/restart_every_single)*restart_every_single +
           restart_every_single;
       else {
         bigint nextrestart = static_cast<bigint>
           (input->variable->compute_equal(ivar_restart_single));
         if (nextrestart <= ntimestep)
           error->all(FLERR,"Restart variable returned a bad timestep");
         next_restart_single = nextrestart;
       }
     } else next_restart_single = update->laststep + 1;
     if (restart_flag_double) {
       if (restart_every_double)
         next_restart_double =
           (ntimestep/restart_every_double)*restart_every_double +
           restart_every_double;
       else {
         bigint nextrestart = static_cast<bigint>
           (input->variable->compute_equal(ivar_restart_double));
         if (nextrestart <= ntimestep)
           error->all(FLERR,"Restart variable returned a bad timestep");
         next_restart_double = nextrestart;
       }
     } else next_restart_double = update->laststep + 1;
     next_restart = MIN(next_restart_single,next_restart_double);
   } else next_restart = update->laststep + 1;
 
   // print memory usage unless being called between multiple runs
 
   if (memflag) memory_usage();
 
   // set next_thermo to multiple of every or variable eval if var defined
   // insure thermo output on last step of run
   // thermo may invoke computes so wrap with clear/add
 
   modify->clearstep_compute();
 
   thermo->header();
   thermo->compute(0);
   last_thermo = ntimestep;
 
   if (var_thermo) {
     next_thermo = static_cast<bigint>
       (input->variable->compute_equal(ivar_thermo));
     if (next_thermo <= ntimestep)
       error->all(FLERR,"Thermo every variable returned a bad timestep");
   } else if (thermo_every) {
     next_thermo = (ntimestep/thermo_every)*thermo_every + thermo_every;
     next_thermo = MIN(next_thermo,update->laststep);
   } else next_thermo = update->laststep;
 
   modify->addstep_compute(next_thermo);
 
   // next = next timestep any output will be done
 
   next = MIN(next_dump_any,next_restart);
   next = MIN(next,next_thermo);
 }
 
 /* ----------------------------------------------------------------------
    perform all output for this timestep
    only perform output if next matches current step and last output doesn't
    do dump/restart before thermo so thermo CPU time will include them
 ------------------------------------------------------------------------- */
 
 void Output::write(bigint ntimestep)
 {
   // next_dump does not force output on last step of run
   // wrap dumps that invoke computes or eval of variable with clear/add
 
   if (next_dump_any == ntimestep) {
     for (int idump = 0; idump < ndump; idump++) {
       if (next_dump[idump] == ntimestep) {
         if (dump[idump]->clearstep || every_dump[idump] == 0)
           modify->clearstep_compute();
         if (last_dump[idump] != ntimestep) {
           dump[idump]->write();
           last_dump[idump] = ntimestep;
         }
         if (every_dump[idump]) next_dump[idump] += every_dump[idump];
         else {
           bigint nextdump = static_cast<bigint>
             (input->variable->compute_equal(ivar_dump[idump]));
           if (nextdump <= ntimestep)
             error->all(FLERR,"Dump every variable returned a bad timestep");
           next_dump[idump] = nextdump;
         }
         if (dump[idump]->clearstep || every_dump[idump] == 0)
           modify->addstep_compute(next_dump[idump]);
       }
       if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]);
       else next_dump_any = next_dump[0];
     }
   }
 
   // next_restart does not force output on last step of run
   // for toggle = 0, replace "*" with current timestep in restart filename
   // eval of variable may invoke computes so wrap with clear/add
 
   if (next_restart == ntimestep) {
     if (next_restart_single == ntimestep) {
       char *file = new char[strlen(restart1) + 16];
       char *ptr = strchr(restart1,'*');
       *ptr = '\0';
       sprintf(file,"%s" BIGINT_FORMAT "%s",restart1,ntimestep,ptr+1);
       *ptr = '*';
       if (last_restart != ntimestep) restart->write(file);
       delete [] file;
       if (restart_every_single) next_restart_single += restart_every_single;
       else {
         modify->clearstep_compute();
         bigint nextrestart = static_cast<bigint>
           (input->variable->compute_equal(ivar_restart_single));
         if (nextrestart <= ntimestep)
           error->all(FLERR,"Restart variable returned a bad timestep");
         next_restart_single = nextrestart;
         modify->addstep_compute(next_restart_single);
       }
     }
     if (next_restart_double == ntimestep) {
       if (last_restart != ntimestep) {
         if (restart_toggle == 0) {
           restart->write(restart2a);
           restart_toggle = 1;
         } else {
           restart->write(restart2b);
           restart_toggle = 0;
         }
       }
       if (restart_every_double) next_restart_double += restart_every_double;
       else {
         modify->clearstep_compute();
         bigint nextrestart = static_cast<bigint>
           (input->variable->compute_equal(ivar_restart_double));
         if (nextrestart <= ntimestep)
           error->all(FLERR,"Restart variable returned a bad timestep");
         next_restart_double = nextrestart;
         modify->addstep_compute(next_restart_double);
       }
     }
     last_restart = ntimestep;
     next_restart = MIN(next_restart_single,next_restart_double);
   }
 
   // insure next_thermo forces output on last step of run
   // thermo may invoke computes so wrap with clear/add
 
   if (next_thermo == ntimestep) {
     modify->clearstep_compute();
     if (last_thermo != ntimestep) thermo->compute(1);
     last_thermo = ntimestep;
     if (var_thermo) {
       next_thermo = static_cast<bigint>
         (input->variable->compute_equal(ivar_thermo));
       if (next_thermo <= ntimestep)
         error->all(FLERR,"Thermo every variable returned a bad timestep");
     } else if (thermo_every) next_thermo += thermo_every;
     else next_thermo = update->laststep;
     next_thermo = MIN(next_thermo,update->laststep);
     modify->addstep_compute(next_thermo);
   }
 
   // next = next timestep any output will be done
 
   next = MIN(next_dump_any,next_restart);
   next = MIN(next,next_thermo);
 }
 
 /* ----------------------------------------------------------------------
    force a snapshot to be written for all dumps
    called from PRD and TAD
 ------------------------------------------------------------------------- */
 
 void Output::write_dump(bigint ntimestep)
 {
   for (int idump = 0; idump < ndump; idump++) {
     dump[idump]->write();
     last_dump[idump] = ntimestep;
   }
 }
 
 /* ----------------------------------------------------------------------
    force restart file(s) to be written
    called from PRD and TAD
 ------------------------------------------------------------------------- */
 
 void Output::write_restart(bigint ntimestep)
 {
   if (restart_flag_single) {
     char *file = new char[strlen(restart1) + 16];
     char *ptr = strchr(restart1,'*');
     *ptr = '\0';
     sprintf(file,"%s" BIGINT_FORMAT "%s",restart1,ntimestep,ptr+1);
     *ptr = '*';
     restart->write(file);
     delete [] file;
   }
 
   if (restart_flag_double) {
     if (restart_toggle == 0) {
       restart->write(restart2a);
       restart_toggle = 1;
     } else {
       restart->write(restart2b);
       restart_toggle = 0;
     }
   }
 
   last_restart = ntimestep;
 }
 
 /* ----------------------------------------------------------------------
    timestep is being changed, called by update->reset_timestep()
    reset next timestep values for dumps, restart, thermo output
    reset to smallest value >= new timestep
    if next timestep set by variable evaluation,
      eval for ntimestep-1, so current ntimestep can be returned if needed
      no guarantee that variable can be evaluated for ntimestep-1
        if it depends on computes, but live with that rare case for now
 ------------------------------------------------------------------------- */
 
 void Output::reset_timestep(bigint ntimestep)
 {
   next_dump_any = MAXBIGINT;
   for (int idump = 0; idump < ndump; idump++) {
     if (every_dump[idump]) {
       next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump];
       if (next_dump[idump] < ntimestep) next_dump[idump] += every_dump[idump];
     } else {
       // ivar_dump may not be initialized
       if (ivar_dump[idump] < 0) {
         ivar_dump[idump] = input->variable->find(var_dump[idump]);
         if (ivar_dump[idump] < 0)
           error->all(FLERR,"Variable name for dump every does not exist");
         if (!input->variable->equalstyle(ivar_dump[idump]))
           error->all(FLERR,"Variable for dump every is invalid style");
       }
       modify->clearstep_compute();
       update->ntimestep--;
       bigint nextdump = static_cast<bigint>
         (input->variable->compute_equal(ivar_dump[idump]));
       if (nextdump < ntimestep)
         error->all(FLERR,"Dump every variable returned a bad timestep");
       update->ntimestep++;
       next_dump[idump] = nextdump;
       modify->addstep_compute(next_dump[idump]);
     }
     next_dump_any = MIN(next_dump_any,next_dump[idump]);
   }
 
   if (restart_flag_single) {
     if (restart_every_single) {
       next_restart_single =
         (ntimestep/restart_every_single)*restart_every_single;
       if (next_restart_single < ntimestep)
         next_restart_single += restart_every_single;
     } else {
       modify->clearstep_compute();
       update->ntimestep--;
       bigint nextrestart = static_cast<bigint>
         (input->variable->compute_equal(ivar_restart_single));
       if (nextrestart < ntimestep)
         error->all(FLERR,"Restart variable returned a bad timestep");
       update->ntimestep++;
       next_restart_single = nextrestart;
       modify->addstep_compute(next_restart_single);
     }
   } else next_restart_single = update->laststep + 1;
 
   if (restart_flag_double) {
     if (restart_every_double) {
       next_restart_double =
         (ntimestep/restart_every_double)*restart_every_double;
       if (next_restart_double < ntimestep)
         next_restart_double += restart_every_double;
     } else {
       modify->clearstep_compute();
       update->ntimestep--;
       bigint nextrestart = static_cast<bigint>
         (input->variable->compute_equal(ivar_restart_double));
       if (nextrestart < ntimestep)
         error->all(FLERR,"Restart variable returned a bad timestep");
       update->ntimestep++;
       next_restart_double = nextrestart;
       modify->addstep_compute(next_restart_double);
     }
   } else next_restart_double = update->laststep + 1;
 
   next_restart = MIN(next_restart_single,next_restart_double);
 
   if (var_thermo) {
     modify->clearstep_compute();
     update->ntimestep--;
     next_thermo = static_cast<bigint>
       (input->variable->compute_equal(ivar_thermo));
     if (next_thermo < ntimestep)
       error->all(FLERR,"Thermo_modify every variable returned a bad timestep");
     update->ntimestep++;
     next_thermo = MIN(next_thermo,update->laststep);
     modify->addstep_compute(next_thermo);
   } else if (thermo_every) {
     next_thermo = (ntimestep/thermo_every)*thermo_every;
     if (next_thermo < ntimestep) next_thermo += thermo_every;
     next_thermo = MIN(next_thermo,update->laststep);
   } else next_thermo = update->laststep;
 
   next = MIN(next_dump_any,next_restart);
   next = MIN(next,next_thermo);
 }
 
 /* ----------------------------------------------------------------------
    add a Dump to list of Dumps
 ------------------------------------------------------------------------- */
 
 void Output::add_dump(int narg, char **arg)
 {
   if (narg < 5) error->all(FLERR,"Illegal dump command");
 
   // error checks
 
   for (int idump = 0; idump < ndump; idump++)
     if (strcmp(arg[0],dump[idump]->id) == 0)
       error->all(FLERR,"Reuse of dump ID");
   int igroup = group->find(arg[1]);
   if (igroup == -1) error->all(FLERR,"Could not find dump group ID");
   if (force->inumeric(FLERR,arg[3]) <= 0)
     error->all(FLERR,"Invalid dump frequency");
 
   // extend Dump list if necessary
 
   if (ndump == max_dump) {
     max_dump += DELTA;
     dump = (Dump **)
       memory->srealloc(dump,max_dump*sizeof(Dump *),"output:dump");
     memory->grow(every_dump,max_dump,"output:every_dump");
     memory->grow(next_dump,max_dump,"output:next_dump");
     memory->grow(last_dump,max_dump,"output:last_dump");
     var_dump = (char **)
       memory->srealloc(var_dump,max_dump*sizeof(char *),"output:var_dump");
     memory->grow(ivar_dump,max_dump,"output:ivar_dump");
   }
 
   // initialize per-dump data to suitable default values
 
   every_dump[ndump] = 0;
   last_dump[ndump] = -1;
   var_dump[ndump] = NULL;
   ivar_dump[ndump] = -1;
 
   // create the Dump
 
   if (dump_map->find(arg[2]) != dump_map->end()) {
     DumpCreator dump_creator = (*dump_map)[arg[2]];
     dump[ndump] = dump_creator(lmp, narg, arg);
   }
   else error->all(FLERR,"Unknown dump style");
 
   every_dump[ndump] = force->inumeric(FLERR,arg[3]);
   if (every_dump[ndump] <= 0) error->all(FLERR,"Illegal dump command");
   last_dump[ndump] = -1;
   var_dump[ndump] = NULL;
   ndump++;
 }
 
 /* ----------------------------------------------------------------------
    one instance per dump style in style_dump.h
 ------------------------------------------------------------------------- */
 
 template <typename T>
 Dump *Output::dump_creator(LAMMPS *lmp, int narg, char ** arg)
 {
   return new T(lmp, narg, arg);
 }
 
 /* ----------------------------------------------------------------------
    modify parameters of a Dump
 ------------------------------------------------------------------------- */
 
 void Output::modify_dump(int narg, char **arg)
 {
   if (narg < 1) error->all(FLERR,"Illegal dump_modify command");
 
   // find which dump it is
 
   int idump;
   for (idump = 0; idump < ndump; idump++)
     if (strcmp(arg[0],dump[idump]->id) == 0) break;
   if (idump == ndump) error->all(FLERR,"Cound not find dump_modify ID");
 
   dump[idump]->modify_params(narg-1,&arg[1]);
 }
 
 /* ----------------------------------------------------------------------
    delete a Dump from list of Dumps
 ------------------------------------------------------------------------- */
 
 void Output::delete_dump(char *id)
 {
   // find which dump it is and delete it
 
   int idump;
   for (idump = 0; idump < ndump; idump++)
     if (strcmp(id,dump[idump]->id) == 0) break;
   if (idump == ndump) error->all(FLERR,"Could not find undump ID");
 
   delete dump[idump];
   delete [] var_dump[idump];
 
   // move other dumps down in list one slot
 
   for (int i = idump+1; i < ndump; i++) {
     dump[i-1] = dump[i];
     every_dump[i-1] = every_dump[i];
     next_dump[i-1] = next_dump[i];
     last_dump[i-1] = last_dump[i];
     var_dump[i-1] = var_dump[i];
     ivar_dump[i-1] = ivar_dump[i];
   }
   ndump--;
 }
 
 /* ----------------------------------------------------------------------
    set thermo output frequency from input script
 ------------------------------------------------------------------------- */
 
 void Output::set_thermo(int narg, char **arg)
 {
   if (narg != 1) error->all(FLERR,"Illegal thermo command");
 
   if (strstr(arg[0],"v_") == arg[0]) {
     delete [] var_thermo;
     int n = strlen(&arg[0][2]) + 1;
     var_thermo = new char[n];
     strcpy(var_thermo,&arg[0][2]);
   } else {
     thermo_every = force->inumeric(FLERR,arg[0]);
     if (thermo_every < 0) error->all(FLERR,"Illegal thermo command");
   }
 }
 
 /* ----------------------------------------------------------------------
    new Thermo style
 ------------------------------------------------------------------------- */
 
 void Output::create_thermo(int narg, char **arg)
 {
   if (narg < 1) error->all(FLERR,"Illegal thermo_style command");
 
   // don't allow this so that dipole style can safely allocate inertia vector
 
   if (domain->box_exist == 0)
     error->all(FLERR,"Thermo_style command before simulation box is defined");
 
   // warn if previous thermo had been modified via thermo_modify command
 
   if (thermo->modified && comm->me == 0)
     error->warning(FLERR,"New thermo_style command, "
                    "previous thermo_modify settings will be lost");
 
   // set thermo = NULL in case new Thermo throws an error
 
   delete thermo;
   thermo = NULL;
   thermo = new Thermo(lmp,narg,arg);
 }
 
 /* ----------------------------------------------------------------------
    setup restart capability for single or double output files
    if only one filename and it contains no "*", then append ".*"
 ------------------------------------------------------------------------- */
 
 void Output::create_restart(int narg, char **arg)
 {
   if (narg < 1) error->all(FLERR,"Illegal restart command");
 
   int every = 0;
   int varflag = 0;
 
   if (strstr(arg[0],"v_") == arg[0]) varflag = 1;
   else every = force->inumeric(FLERR,arg[0]);
 
   if (!varflag && every == 0) {
     if (narg != 1) error->all(FLERR,"Illegal restart command");
 
     restart_flag = restart_flag_single = restart_flag_double = 0;
     last_restart = -1;
 
     delete restart;
     restart = NULL;
     delete [] restart1;
     delete [] restart2a;
     delete [] restart2b;
     restart1 = restart2a = restart2b = NULL;
     delete [] var_restart_single;
     delete [] var_restart_double;
     var_restart_single = var_restart_double = NULL;
 
     return;
   }
 
   if (narg < 2) error->all(FLERR,"Illegal restart command");
 
   int nfile = 0;
   if (narg % 2 == 0) nfile = 1;
   else nfile = 2;
 
   if (nfile == 1) {
     restart_flag = restart_flag_single = 1;
 
     if (varflag) {
       delete [] var_restart_single;
       int n = strlen(&arg[0][2]) + 1;
       var_restart_single = new char[n];
       strcpy(var_restart_single,&arg[0][2]);
       restart_every_single = 0;
     } else restart_every_single = every;
 
     int n = strlen(arg[1]) + 3;
     delete [] restart1;
     restart1 = new char[n];
     strcpy(restart1,arg[1]);
     if (strchr(restart1,'*') == NULL) strcat(restart1,".*");
   }
 
   if (nfile == 2) {
     restart_flag = restart_flag_double = 1;
 
     if (varflag) {
       delete [] var_restart_double;
       int n = strlen(&arg[0][2]) + 1;
       var_restart_double = new char[n];
       strcpy(var_restart_double,&arg[0][2]);
       restart_every_double = 0;
     } else restart_every_double = every;
 
     delete [] restart2a;
     delete [] restart2b;
     restart_toggle = 0;
     int n = strlen(arg[1]) + 3;
     restart2a = new char[n];
     strcpy(restart2a,arg[1]);
     n = strlen(arg[2]) + 1;
     restart2b = new char[n];
     strcpy(restart2b,arg[2]);
   }
 
   // check for multiproc output and an MPI-IO filename
   // if 2 filenames, must be consistent
 
   int multiproc;
   if (strchr(arg[1],'%')) multiproc = comm->nprocs;
   else multiproc = 0;
   if (nfile == 2) {
     if (multiproc && !strchr(arg[2],'%'))
       error->all(FLERR,"Both restart files must use % or neither");
     if (!multiproc && strchr(arg[2],'%'))
       error->all(FLERR,"Both restart files must use % or neither");
   }
 
   int mpiioflag;
   if (strstr(arg[1],".mpi")) mpiioflag = 1;
   else mpiioflag = 0;
   if (nfile == 2) {
     if (mpiioflag && !strstr(arg[2],".mpi"))
       error->all(FLERR,"Both restart files must use MPI-IO or neither");
     if (!mpiioflag && strstr(arg[2],".mpi"))
       error->all(FLERR,"Both restart files must use MPI-IO or neither");
   }
 
   // setup output style and process optional args
 
   delete restart;
   restart = new WriteRestart(lmp);
   int iarg = nfile+1;
   restart->multiproc_options(multiproc,mpiioflag,narg-iarg,&arg[iarg]);
 }
 
 /* ----------------------------------------------------------------------
    sum and print memory usage
    result is only memory on proc 0, not averaged across procs
 ------------------------------------------------------------------------- */
-
 void Output::memory_usage()
 {
+
   bigint bytes = 0;
   bytes += atom->memory_usage();
   bytes += neighbor->memory_usage();
   bytes += comm->memory_usage();
   bytes += update->memory_usage();
   bytes += force->memory_usage();
   bytes += modify->memory_usage();
   for (int i = 0; i < ndump; i++) bytes += dump[i]->memory_usage();
 
   double mbytes = bytes/1024.0/1024.0;
+  double mbavg,mbmin,mbmax;
+  MPI_Reduce(&mbytes,&mbavg,1,MPI_DOUBLE,MPI_SUM,0,world);
+  MPI_Reduce(&mbytes,&mbmin,1,MPI_DOUBLE,MPI_MIN,0,world);
+  MPI_Reduce(&mbytes,&mbmax,1,MPI_DOUBLE,MPI_MAX,0,world);
+  mbavg /= comm->nprocs;
 
   if (comm->me == 0) {
     if (screen)
-      fprintf(screen,"Memory usage per processor = %g Mbytes\n",mbytes);
+      fprintf(screen,"Per MPI rank memory allocation (min/avg/max) = "
+              "%.4g | %.4g | %.4g Mbytes\n",mbmin,mbavg,mbmax);
     if (logfile)
-      fprintf(logfile,"Memory usage per processor = %g Mbytes\n",mbytes);
+      fprintf(logfile,"Per MPI rank memory allocation (min/avg/max) = "
+              "%.4g | %.4g | %.4g Mbytes\n",mbmin,mbavg,mbmax);
   }
 }