From 2aef3a0e9670f32459dd5100dd23bb484cfcc663 Mon Sep 17 00:00:00 2001
From: "Steven J. Plimpton" <sjplimp@singsing.sandia.gov>
Date: Tue, 27 Mar 2018 14:37:04 -0600
Subject: [PATCH] new dump_modify refresh and compute displace/atom refresh
commands for incremental dump files
---
doc/src/compute_displace_atom.txt | 73 +++++++++++++++++++--
doc/src/dump.txt | 6 +-
doc/src/dump_modify.txt | 103 ++++++++++++++++++++++++------
src/compute.h | 2 +
src/compute_displace_atom.cpp | 70 +++++++++++++++++++-
src/compute_displace_atom.h | 5 ++
src/dump.cpp | 25 +++++++-
src/dump.h | 4 ++
src/dump_custom.cpp | 16 ++++-
9 files changed, 272 insertions(+), 32 deletions(-)
diff --git a/doc/src/compute_displace_atom.txt b/doc/src/compute_displace_atom.txt
index decf0bc8b9..0dbea1dcd6 100644
--- a/doc/src/compute_displace_atom.txt
+++ b/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
-displace/atom = style name of this compute command :ul
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+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
-due to atoms passing thru periodic boundaries.
+atom in the group from its original (reference) coordinates, including
+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:]
diff --git a/doc/src/dump.txt b/doc/src/dump.txt
index 1e6acc67e6..42ddbc47dc 100644
--- a/doc/src/dump.txt
+++ b/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
-alter what atoms are included. Not all styles support all these
-options; see details below.
+"dump_modify thresh and region and refresh"_dump_modify.html commands
+can also alter what atoms are included. Not all styles support all
+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
diff --git a/doc/src/dump_modify.txt b/doc/src/dump_modify.txt
index db727c2d4f..fddf88875f 100644
--- a/doc/src/dump_modify.txt
+++ b/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}
-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 {thermo} keyword ({netcdf} only) triggers writing of "thermo"_thermo.html
-information to the dump file alongside per-atom data. The data included in the
-dump file is identical to the data specified by
-"thermo_style"_thermo_style.html.
+The {refresh} keyword only applies to the dump {custom}, {cfg},
+{image}, and {movie} styles. It allows an "incremental" dump file to
+be written, by refreshing a compute that is used as a threshold for
+determining which atoms are included in a dump snapshot. The
+specified {c_ID} gives the ID of the compute. It is prefixed by "c_"
+to indicate a compute, which is the only current option. At some
+point, other options may be added, e.g. fixes or variables.
+
+NOTE: This keyword can only be specified once for a dump. Refreshes
+of multiple computes cannot yet be performed.
+
+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 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
diff --git a/src/compute.h b/src/compute.h
index f04ffebd61..a023834368 100644
--- a/src/compute.h
+++ b/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();
diff --git a/src/compute_displace_atom.cpp b/src/compute_displace_atom.cpp
index 9b9e165699..a5c32c53cc 100644
--- a/src/compute_displace_atom.cpp
+++ b/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;
}
diff --git a/src/compute_displace_atom.h b/src/compute_displace_atom.h
index 5b62f747e1..05c729daa4 100644
--- a/src/compute_displace_atom.h
+++ b/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;
};
}
diff --git a/src/dump.cpp b/src/dump.cpp
index ddd958c25c..6ae75a8824 100644
--- a/src/dump.cpp
+++ b/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) {
diff --git a/src/dump.h b/src/dump.h
index a5582ce5a5..5eabf25e08 100644
--- a/src/dump.h
+++ b/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
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index 9892387cbe..34ab218204 100644
--- a/src/dump_custom.cpp
+++ b/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) {
--
GitLab