new dump_modify refresh and compute displace/atom refresh commands for incremental dump files

doc/src/compute_displace_atom.txt

@@ -12,18 +12,23 @@ compute displace/atom command :h3
 
 compute ID group-ID displace/atom :pre
 
-ID, group-ID are documented in "compute"_compute.html command
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
-displace/atom = style name of this compute command :ul
+displace/atom = style name of this compute command :l
+zero or more keyword/arg pairs may be appended :l
+keyword = {refresh} :l
+  {replace} arg = per-atom variable ID :pre
+:ule
 
 [Examples:]
 
-compute 1 all displace/atom :pre
+compute 1 all displace/atom
+compute 1 all displace/atom refresh myVar :pre
 
 [Description:]
 
 Define a computation that calculates the current displacement of each
-atom in the group from its original coordinates, including all effects
+atom in the group from its original (reference) coordinates, including
-due to atoms passing thru periodic boundaries.
+all effects due to atoms passing thru periodic boundaries.
 
 A vector of four quantities per atom is calculated by this compute.
 The first 3 elements of the vector are the dx,dy,dz displacements.
@@ -49,6 +54,60 @@ This is so that the fix this compute creates to store per-atom
 quantities will also have the same ID, and thus be initialized
 correctly with time=0 atom coordinates from the restart file.
 
+:line
+
+The {refresh} option can be used in conjuction with the "dump_modify
+refresh" command to generate incremental dump files.
+
+The definition and motivation of an incremental dump file is as
+follows.  Instead of outputting all atoms at each snapshot (with some
+associated values), you may only wish to output the subset of atoms
+with a value that has changed in some way compared to the value the
+last time that atom was output.  In some scenarios this can result in
+a dramatically smaller dump file.  If desired, by post-processing the
+sequence of snapshots, the values for all atoms at all timesteps can
+be inferred.
+
+A concrete example using this compute, is a simulation of atom
+diffusion in a solid, represented as atoms on a lattice.  Diffusive
+hops are rare.  Imagine that when a hop occurs an atom moves more than
+a distance {Dhop}.  For any snapshot we only want to output atoms that
+have hopped since the last snapshot.  This can be accomplished with
+something like the following commands:
+
+variable        Dhop equal 0.6
+variable        check atom "c_dsp[4] > v_Dhop"
+compute         dsp all displace/atom refresh check
+dump            1 all custom 20 tmp.dump id type x y z
+dump_modify     1 append yes thresh c_dsp[4] > ${Dhop} refresh c_dsp :pre
+
+The "dump_modify thresh"_dump_modify.html command will cause only
+atoms that have displaced more than 0.6 Angstroms to be output on a
+given snapshot (assuming metal units).  The dump_modify {refresh}
+option triggers a call to this compute at the end of every dump.
+
+The argument for compute displace/atom refresh is the ID of an
+"atom-style variable"_variable.html which calculates a Boolean value
+(0 or 1) based on the same criterion used by dump_modify thresh.  What
+this compute will then do is evaluate the atom-style variable and for
+each atom that returns 1 (true), it will update the original
+(reference) coordinates of the atom, which are stored by this compute.
+
+The effect of these commands is that a particular atom will only be
+output in the dump file on shapshots following a diffusive hop.  It
+will not be output again until it performs another hop.
+
+Note that the first snapshot (geneated by the first run after the dump
+and this compute are defined), no atoms will be typically be output.
+That is because the initial displacement for all atoms is 0.0.  If an
+initial dump snapshot is desired, containing the initial reference
+position of all atoms, a command like this can be used before the
+first simulation is performed
+
+write_dump      all custom tmp.dump.initial id type x y z :pre
+
+:line
+
 [Output info:]
 
 This compute calculates a per-atom array with 4 columns, which can be
@@ -59,6 +118,10 @@ options.
 
 The per-atom array values will be in distance "units"_units.html.
 
+This compute supports the {refresh} option as explained above, for use
+in conjunction with "dump_modify refresh"_dump_modify.html to generate
+incremental dump files.
+
 [Restrictions:] none
 
 [Related commands:]

doc/src/dump.txt

@@ -120,9 +120,9 @@ which dump output is written can also be controlled by a variable.
 See the "dump_modify every"_dump_modify.html command.
 
 Only information for atoms in the specified group is dumped.  The
-"dump_modify thresh and region"_dump_modify.html commands can also
+"dump_modify thresh and region and refresh"_dump_modify.html commands
-alter what atoms are included.  Not all styles support all these
+can also alter what atoms are included.  Not all styles support all
-options; see details below.
+these options; see details below.
 
 As described below, the filename determines the kind of output (text
 or binary or gzipped, one big file or one per timestep, one big file

doc/src/dump_modify.txt

@@ -41,6 +41,7 @@ keyword = {append} or {at} or {buffer} or {element} or {every} or {fileper} or {
   {pbc} arg = {yes} or {no} = remap atoms via periodic boundary conditions
   {precision} arg = power-of-10 value from 10 to 1000000
   {region} arg = region-ID or "none"
+  {refresh} arg = c_ID = compute ID that supports a refresh operation
   {scale} arg = {yes} or {no}
   {sfactor} arg = coordinate scaling factor (> 0.0)
   {thermo} arg = {yes} or {no}
@@ -407,25 +408,67 @@ nanometer accuracy, e.g. for N = 1000, the coordinates are written to
 
 :line
 
-The {sfactor} and {tfactor} keywords only apply to the dump {xtc}
+The {refresh} keyword only applies to the dump {custom}, {cfg},
-style.  They allow customization of the unit conversion factors used
+{image}, and {movie} styles.  It allows an "incremental" dump file to
-when writing to XTC files.  By default they are initialized for
+be written, by refreshing a compute that is used as a threshold for
-whatever "units"_units.html style is being used, to write out
+determining which atoms are included in a dump snapshot.  The
-coordinates in nanometers and time in picoseconds.  I.e. for {real}
+specified {c_ID} gives the ID of the compute.  It is prefixed by "c_"
-units, LAMMPS defines {sfactor} = 0.1 and {tfactor} = 0.001, since the
+to indicate a compute, which is the only current option.  At some
-Angstroms and fmsec used by {real} units are 0.1 nm and 0.001 psec
+point, other options may be added, e.g. fixes or variables.
-respectively.  If you are using a units system with distance and time
-units far from nm and psec, you may wish to write XTC files with
-different units, since the compression algorithm used in XTC files is
-most effective when the typical magnitude of position data is between
-10.0 and 0.1.
 
-:line
+NOTE: This keyword can only be specified once for a dump.  Refreshes
+of multiple computes cannot yet be performed.
 
-The {thermo} keyword ({netcdf} only) triggers writing of "thermo"_thermo.html
+The definition and motivation of an incremental dump file is as
-information to the dump file alongside per-atom data. The data included in the
+follows.  Instead of outputting all atoms at each snapshot (with some
-dump file is identical to the data specified by
+associated values), you may only wish to output the subset of atoms
-"thermo_style"_thermo_style.html.
+with a value that has changed in some way compared to the value the
+last time that atom was output.  In some scenarios this can result in
+a dramatically smaller dump file.  If desired, by post-processing the
+sequence of snapshots, the values for all atoms at all timesteps can
+be inferred.
+
+A concrete example is a simulation of atom diffusion in a solid,
+represented as atoms on a lattice.  Diffusive hops are rare.  Imagine
+that when a hop occurs an atom moves more than a distance {Dhop}.  For
+any snapshot we only want to output atoms that have hopped since the
+last snapshot.  This can be accomplished with something the following
+commands:
+
+variable        Dhop equal 0.6
+variable        check atom "c_dsp[4] > v_Dhop"
+compute         dsp all displace/atom refresh check
+dump            1 all custom 20 tmp.dump id type x y z
+dump_modify     1 append yes thresh c_dsp[4] > ${Dhop} refresh c_dsp :pre
+
+The "compute displace/atom"_compute_displace_atom.html command
+calculates the displacement of each atom from its reference position.
+The "4" index is the scalar displacement; 1,2,3 are the xyz components
+of the displacement.  The "dump_modify thresh"_dump_modify.html
+command will cause only atoms that have displaced more than 0.6
+Angstroms to be output on a given snapshot (assuming metal units).
+However, note that when an atom is output, we also need to update the
+reference position for that atom to its new coordinates.  So that it
+will not be output in every snapshot thereafter.  That reference
+position is stored by "compute
+displace/atom"_compute_displace_atom.html.  So the dump_modify
+{refresh} option triggers a call to compute displace/atom at the end
+of every dump to perform that update.  The {refresh check} option
+shown as part of the "compute
+displace/atom"_compute_displace_atom.html command enables the compute
+to respond to the call from the dump command, and update the
+appropriate reference positions.  This is done be defining an
+"atom-style variable"_variable.html, {check} in this example, which
+calculates a Boolean value (0 or 1) for each atom, based on the same
+criterion used by dump_modify thresh.  See the "compute
+displace/atom"_compute_displace_atom.html command for more details.
+
+Note that only computes with a {refresh} option will work with
+dump_modify refresh.  See individual compute doc pages for details.
+Currently, only compute displace/atom supports this option.  Others
+may be added at some point.  If you use a compute that doesn't support
+refresh operations, LAMMPS will not complain; dump_modify refresh will
+simply do nothing.
 
 :line
 
@@ -449,6 +492,21 @@ value of {no} means they are written in absolute distance units
 
 :line
 
+The {sfactor} and {tfactor} keywords only apply to the dump {xtc}
+style.  They allow customization of the unit conversion factors used
+when writing to XTC files.  By default they are initialized for
+whatever "units"_units.html style is being used, to write out
+coordinates in nanometers and time in picoseconds.  I.e. for {real}
+units, LAMMPS defines {sfactor} = 0.1 and {tfactor} = 0.001, since the
+Angstroms and fmsec used by {real} units are 0.1 nm and 0.001 psec
+respectively.  If you are using a units system with distance and time
+units far from nm and psec, you may wish to write XTC files with
+different units, since the compression algorithm used in XTC files is
+most effective when the typical magnitude of position data is between
+10.0 and 0.1.
+
+:line
+
 The {sort} keyword determines whether lines of per-atom output in a
 snapshot are sorted or not.  A sort value of {off} means they will
 typically be written in indeterminate order, either in serial or
@@ -472,6 +530,13 @@ as well as memory, versus unsorted output.
 
 :line
 
+The {thermo} keyword only applies the dump {netcdf} style.  It
+triggers writing of "thermo"_thermo.html information to the dump file
+alongside per-atom data.  The values included in the dump file are
+identical to the values specified by "thermo_style"_thermo_style.html.
+
+:line
+
 The {thresh} keyword only applies to the dump {custom}, {cfg},
 {image}, and {movie} styles.  Multiple thresholds can be specified.
 Specifying {none} turns off all threshold criteria.  If thresholds are

src/compute.h

@@ -129,6 +129,8 @@ class Compute : protected Pointers {
   virtual void lock(class Fix *, bigint, bigint) {}
   virtual void unlock(class Fix *) {}
 
+  virtual void refresh() {}
+
   void addstep(bigint);
   int matchstep(bigint);
   void clearstep();

src/compute_displace_atom.cpp

@@ -21,6 +21,8 @@
 #include "modify.h"
 #include "fix.h"
 #include "fix_store.h"
+#include "input.h"
+#include "variable.h"
 #include "memory.h"
 #include "error.h"
 
@@ -32,12 +34,42 @@ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
   displace(NULL), id_fix(NULL)
 {
-  if (narg != 3) error->all(FLERR,"Illegal compute displace/atom command");
+  if (narg < 3) error->all(FLERR,"Illegal compute displace/atom command");
 
   peratom_flag = 1;
   size_peratom_cols = 4;
   create_attribute = 1;
 
+  // optional args
+
+  refreshflag = 0;
+  rvar = NULL;
+
+  int iarg = 3;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"refresh") == 0) {
+      if (iarg+2 > narg) 
+        error->all(FLERR,"Illegal compute displace/atom command");
+      refreshflag = 1;
+      delete [] rvar;
+      int n = strlen(arg[iarg+1]) + 1;
+      rvar = new char[n];
+      strcpy(rvar,arg[iarg+1]);
+      iarg += 2;
+    } else error->all(FLERR,"Illegal compute displace/atom command");
+  }
+
+  // error check
+
+  if (refreshflag) {
+    ivar = input->variable->find(rvar);
+    if (ivar < 0) 
+      error->all(FLERR,"Variable name for compute displace/atom does not exist");
+    if (input->variable->atomstyle(ivar) == 0)
+      error->all(FLERR,"Compute displace/atom variable "
+                 "is not atom-style variable");
+  }
+
   // create a new fix STORE style
   // id = compute-ID + COMPUTE_STORE, fix group = compute group
 
@@ -76,7 +108,8 @@ ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) :
 
   // per-atom displacement array
 
-  nmax = 0;
+  nmax = nvmax = 0;
+  varatom = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -89,6 +122,8 @@ ComputeDisplaceAtom::~ComputeDisplaceAtom()
 
   delete [] id_fix;
   memory->destroy(displace);
+  delete [] rvar;
+  memory->destroy(varatom);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -100,6 +135,12 @@ void ComputeDisplaceAtom::init()
   int ifix = modify->find_fix(id_fix);
   if (ifix < 0) error->all(FLERR,"Could not find compute displace/atom fix ID");
   fix = (FixStore *) modify->fix[ifix];
+
+  if (refreshflag) {
+    ivar = input->variable->find(rvar);
+    if (ivar < 0)
+      error->all(FLERR,"Variable name for compute displace/atom does not exist");
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -183,6 +224,30 @@ void ComputeDisplaceAtom::set_arrays(int i)
   xoriginal[i][2] = x[i][2];
 }
 
+/* ----------------------------------------------------------------------
+   reset per-atom storage values, based on atom-style variable evaluation
+   called by dump when dump_modify refresh is set
+------------------------------------------------------------------------- */
+
+void ComputeDisplaceAtom::refresh()
+{
+  if (atom->nmax > nvmax) {
+    nvmax = atom->nmax;
+    memory->destroy(varatom);
+    memory->create(varatom,nvmax,"displace/atom:varatom");
+  }
+
+  input->variable->compute_atom(ivar,igroup,varatom,1,0);
+
+  double **xoriginal = fix->astore;
+  double **x = atom->x;
+  imageint *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (varatom[i]) domain->unmap(x[i],image[i],xoriginal[i]);
+}
+
 /* ----------------------------------------------------------------------
    memory usage of local atom-based array
 ------------------------------------------------------------------------- */
@@ -190,5 +255,6 @@ void ComputeDisplaceAtom::set_arrays(int i)
 double ComputeDisplaceAtom::memory_usage()
 {
   double bytes = nmax*4 * sizeof(double);
+  bytes += nvmax * sizeof(double);
   return bytes;
 }

src/compute_displace_atom.h

@@ -31,6 +31,7 @@ class ComputeDisplaceAtom : public Compute {
   void init();
   void compute_peratom();
   void set_arrays(int);
+  void refresh();
   double memory_usage();
 
  private:
@@ -38,6 +39,10 @@ class ComputeDisplaceAtom : public Compute {
   double **displace;
   char *id_fix;
   class FixStore *fix;
+
+  int refreshflag,ivar,nvmax;    // refresh option is enabled
+  char *rvar;                    // for incremental dumps
+  double *varatom;
 };
 
 }

src/dump.cpp

@@ -22,11 +22,12 @@
 #include "domain.h"
 #include "group.h"
 #include "output.h"
-#include "memory.h"
-#include "error.h"
 #include "force.h"
 #include "modify.h"
 #include "fix.h"
+#include "compute.h"
+#include "memory.h"
+#include "error.h"
 
 using namespace LAMMPS_NS;
 
@@ -78,6 +79,9 @@ Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   format_bigint_user = NULL;
   format_column_user = NULL;
 
+  refreshflag = 0;
+  refresh = NULL;
+
   clearstep = 0;
   sort_flag = 0;
   append_flag = 0;
@@ -160,6 +164,8 @@ Dump::~Dump()
   delete [] format_int_user;
   delete [] format_bigint_user;
 
+  delete [] refresh;
+
   // format_column_user is deallocated by child classes that use it
 
   memory->destroy(buf);
@@ -277,6 +283,16 @@ void Dump::init()
     }
   }
 
+  // search for refresh compute specified by dump_modify refresh
+
+  if (refreshflag) {
+    int icompute;
+    for (icompute = 0; icompute < modify->ncompute; icompute++)
+      if (strcmp(refresh,modify->compute[icompute]->id) == 0) break;
+    if (icompute < modify->ncompute) irefresh = icompute;
+    else error->all(FLERR,"Dump could not find refresh compute ID");
+  }
+
   // preallocation for PBC copies if requested
 
   if (pbcflag && atom->nlocal > maxpbc) pbc_allocate();
@@ -478,6 +494,11 @@ void Dump::write()
     atom->image = imagehold;
   }
 
+  // trigger post-dump refresh by specified compute
+  // currently used for incremental dump files
+
+  if (refreshflag) modify->compute[irefresh]->refresh();
+
   // if file per timestep, close file if I am filewriter
 
   if (multifile) {

src/dump.h

@@ -78,6 +78,10 @@ class Dump : protected Pointers {
   int sortcolm1;             // sortcol - 1
   int sortorder;             // ASCEND or DESCEND
 
+  int refreshflag;           // 1 if dump_modify refresh specified
+  char *refresh;             // compute ID to invoke refresh() on
+  int irefresh;              // index of compute
+
   char boundstr[9];          // encoding of boundary flags
 
   char *format;              // format string for the file write

src/dump_custom.cpp

@@ -60,7 +60,8 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
   earg(NULL), vtype(NULL), vformat(NULL), columns(NULL), choose(NULL),
   dchoose(NULL), clist(NULL), field2index(NULL), argindex(NULL), id_compute(NULL),
   compute(NULL), id_fix(NULL), fix(NULL), id_variable(NULL), variable(NULL),
-  vbuf(NULL), id_custom(NULL), flag_custom(NULL), typenames(NULL), pack_choice(NULL)
+  vbuf(NULL), id_custom(NULL), flag_custom(NULL), typenames(NULL), 
+  pack_choice(NULL)
 {
   if (narg == 5) error->all(FLERR,"No dump custom arguments specified");
 
@@ -1673,6 +1674,19 @@ int DumpCustom::modify_param(int narg, char **arg)
     return ntypes+1;
   }
 
+  if (strcmp(arg[0],"refresh") == 0) {
+    if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
+    if (strncmp(arg[1],"c_",2) != 0) 
+      error->all(FLERR,"Illegal dump_modify command");
+    if (refreshflag) error->all(FLERR,"Dump modify can only have one refresh");
+
+    refreshflag = 1;
+    int n = strlen(arg[1]);
+    refresh = new char[n];
+    strcpy(refresh,&arg[1][2]);
+    return 2;
+  }
+
   if (strcmp(arg[0],"thresh") == 0) {
     if (narg < 2) error->all(FLERR,"Illegal dump_modify command");
     if (strcmp(arg[1],"none") == 0) {