add optional nmax keyword to fix vector to allow using it as a sliding window

This commit is contained in:
Axel Kohlmeyer
2023-07-22 00:08:17 -04:00
parent a4a206e601
commit a48f4597a2
3 changed files with 55 additions and 25 deletions

View File

@ -8,11 +8,12 @@ Syntax
.. code-block:: LAMMPS .. code-block:: LAMMPS
fix ID group-ID vector Nevery value1 value2 ... fix ID group-ID vector Nevery [nmax <length>] value1 value2 ...
* ID, group-ID are documented in :doc:`fix <fix>` command * ID, group-ID are documented in :doc:`fix <fix>` command
* vector = style name of this fix command * vector = style name of this fix command
* Nevery = use input values every this many timesteps * Nevery = use input values every this many timesteps
* nmax <length> = set maximal length of vector to <length> (optional)
* one or more input values can be listed * one or more input values can be listed
* value = c_ID, c_ID[N], f_ID, f_ID[N], v_name * value = c_ID, c_ID[N], f_ID, f_ID[N], v_name
@ -32,6 +33,7 @@ Examples
fix 1 all vector 100 c_myTemp fix 1 all vector 100 c_myTemp
fix 1 all vector 5 c_myTemp v_integral fix 1 all vector 5 c_myTemp v_integral
fix 1 all vector 50 nmax 200 c_myTemp
Description Description
""""""""""" """""""""""
@ -43,6 +45,10 @@ values, they are stored as rows in a global array, whose number of
rows is growing. The resulting vector or array can be used by other rows is growing. The resulting vector or array can be used by other
:doc:`output commands <Howto_output>`. :doc:`output commands <Howto_output>`.
The optional *nmax* keyword can be used to restrict the length of the
vector to the given *length* value. Once the vector is filled, the
oldest entries will be discarded when new entries are added.
One way to to use this command is to accumulate a vector that is One way to to use this command is to accumulate a vector that is
time-integrated using the :doc:`variable trap() <variable>` function. time-integrated using the :doc:`variable trap() <variable>` function.
For example the velocity auto-correlation function (VACF) can be For example the velocity auto-correlation function (VACF) can be
@ -94,11 +100,12 @@ calculated by the compute is used.
Note that there is a :doc:`compute reduce <compute_reduce>` command Note that there is a :doc:`compute reduce <compute_reduce>` command
which can sum per-atom quantities into a global scalar or vector which which can sum per-atom quantities into a global scalar or vector which
can thus be accessed by fix vector. Or it can be a compute defined can thus be accessed by fix vector. Or it can be a compute defined not
not in your input script, but by :doc:`thermodynamic output <thermo_style>` or other fixes such as :doc:`fix nvt <fix_nh>` in your input script, but by :doc:`thermodynamic output <thermo_style>`
or :doc:`fix temp/rescale <fix_temp_rescale>`. See the doc pages for or other fixes such as :doc:`fix nvt <fix_nh>` or :doc:`fix temp/rescale
these commands which give the IDs of these computes. Users can also <fix_temp_rescale>`. See the doc pages for these commands which give
write code for their own compute styles and :doc:`add them to LAMMPS <Modify>`. the IDs of these computes. Users can also write code for their own
compute styles and :doc:`add them to LAMMPS <Modify>`.
If a value begins with "f\_", a fix ID must follow which has been If a value begins with "f\_", a fix ID must follow which has been
previously defined in the input script. If no bracketed term is previously defined in the input script. If no bracketed term is
@ -108,7 +115,8 @@ calculated by the fix is used.
Note that some fixes only produce their values on certain timesteps, Note that some fixes only produce their values on certain timesteps,
which must be compatible with *Nevery*, else an error will result. which must be compatible with *Nevery*, else an error will result.
Users can also write code for their own fix styles and :doc:`add them to LAMMPS <Modify>`. Users can also write code for their own fix styles and :doc:`add them to
LAMMPS <Modify>`.
If a value begins with "v\_", a variable name must follow which has If a value begins with "v\_", a variable name must follow which has
been previously defined in the input script. An equal-style or been previously defined in the input script. An equal-style or
@ -126,8 +134,9 @@ quantities to be stored by fix vector.
Restart, fix_modify, output, run start/stop, minimize info Restart, fix_modify, output, run start/stop, minimize info
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""""""""""""""""""""""""""""""
No information about this fix is written to :doc:`binary restart files <restart>`. None of the :doc:`fix_modify <fix_modify>` options No information about this fix is written to :doc:`binary restart files
are relevant to this fix. <restart>`. None of the :doc:`fix_modify <fix_modify>` options are
relevant to this fix.
This fix produces a global vector or global array which can be This fix produces a global vector or global array which can be
accessed by various :doc:`output commands <Howto_output>`. The values accessed by various :doc:`output commands <Howto_output>`. The values
@ -144,15 +153,15 @@ the vector are "intensive" or "extensive". If the fix produces an
array, then all elements in the array must be the same, either array, then all elements in the array must be the same, either
"intensive" or "extensive". If a compute or fix provides the value "intensive" or "extensive". If a compute or fix provides the value
stored, then the compute or fix determines whether the value is stored, then the compute or fix determines whether the value is
intensive or extensive; see the page for that compute or fix for intensive or extensive; see the page for that compute or fix for further
further info. Values produced by a variable are treated as intensive. info. Values produced by a variable are treated as intensive.
This fix can allocate storage for stored values accumulated over This fix can allocate storage for stored values accumulated over
multiple runs, using the *start* and *stop* keywords of the multiple runs, using the *start* and *stop* keywords of the :doc:`run
:doc:`run <run>` command. See the :doc:`run <run>` command for details of <run>` command. See the :doc:`run <run>` command for details of how to
how to do this. If using the :doc:`run pre no <run>` command option, do this. If using the :doc:`run pre no <run>` command option, this is
this is required to allow the fix to allocate sufficient storage for required to allow the fix to allocate sufficient storage for stored
stored values. values.
This fix is not invoked during :doc:`energy minimization <minimize>`. This fix is not invoked during :doc:`energy minimization <minimize>`.

View File

@ -35,10 +35,20 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
nevery = utils::inumeric(FLERR, arg[3], false, lmp); nevery = utils::inumeric(FLERR, arg[3], false, lmp);
if (nevery <= 0) error->all(FLERR, "Invalid fix vector every argument: {}", nevery); if (nevery <= 0) error->all(FLERR, "Invalid fix vector every argument: {}", nevery);
nmaxval = MAXSMALLINT;
nindex = 0;
int iarg = 4;
if (strcmp(arg[iarg], "nmax") == 0) {
nmaxval = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
if (nmaxval < 1) error->all(FLERR, "Invalid nmax value");
iarg += 2;
}
// parse values // parse values
values.clear(); values.clear();
for (int iarg = 4; iarg < narg; iarg++) { while (iarg < narg) {
ArgInfo argi(arg[iarg]); ArgInfo argi(arg[iarg]);
value_t val; value_t val;
@ -51,6 +61,7 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR, "Invalid fix vector argument: {}", arg[iarg]); error->all(FLERR, "Invalid fix vector argument: {}", arg[iarg]);
values.push_back(val); values.push_back(val);
++iarg;
} }
// setup and error check // setup and error check
@ -132,7 +143,7 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) :
vector = nullptr; vector = nullptr;
array = nullptr; array = nullptr;
ncount = ncountmax = 0; ncount = ncountmax = nindex = 0;
if (values.size() == 1) if (values.size() == 1)
size_vector = 0; size_vector = 0;
else else
@ -199,6 +210,7 @@ void FixVector::init()
bigint finalstep = update->endstep / nevery * nevery; bigint finalstep = update->endstep / nevery * nevery;
if (finalstep > update->endstep) finalstep -= nevery; if (finalstep > update->endstep) finalstep -= nevery;
ncountmax = (finalstep - initialstep) / nevery + 1; ncountmax = (finalstep - initialstep) / nevery + 1;
if (ncountmax > nmaxval) ncountmax = nmaxval;
if (values.size() == 1) if (values.size() == 1)
memory->grow(vector, ncountmax, "vector:vector"); memory->grow(vector, ncountmax, "vector:vector");
else else
@ -221,16 +233,18 @@ void FixVector::end_of_step()
// skip if not step which requires doing something // skip if not step which requires doing something
if (update->ntimestep != nextstep) return; if (update->ntimestep != nextstep) return;
if (ncount == ncountmax) error->all(FLERR, "Overflow of allocated fix vector storage");
// wrap around when vector/array is full
nindex = ncount % ncountmax;
// accumulate results of computes,fixes,variables to local copy // accumulate results of computes,fixes,variables to local copy
// compute/fix/variable may invoke computes so wrap with clear/add // compute/fix/variable may invoke computes so wrap with clear/add
double *result; double *result;
if (values.size() == 1) if (values.size() == 1)
result = &vector[ncount]; result = &vector[nindex];
else else
result = array[ncount]; result = array[nindex];
modify->clearstep_compute(); modify->clearstep_compute();
@ -290,9 +304,9 @@ void FixVector::end_of_step()
ncount++; ncount++;
if (values.size() == 1) if (values.size() == 1)
size_vector++; size_vector = MIN(size_vector + 1, ncountmax);
else else
size_array_rows++; size_array_rows = MIN(size_array_rows + 1, ncountmax);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -301,7 +315,9 @@ void FixVector::end_of_step()
double FixVector::compute_vector(int i) double FixVector::compute_vector(int i)
{ {
return vector[i]; int idx = i;
if (ncount >= ncountmax) idx = (i + ncount) % ncountmax;
return vector[idx];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -310,5 +326,7 @@ double FixVector::compute_vector(int i)
double FixVector::compute_array(int i, int j) double FixVector::compute_array(int i, int j)
{ {
return array[i][j]; int idx = i;
if (ncount >= ncountmax) idx = (i + ncount) % ncountmax;
return array[idx][j];
} }

View File

@ -52,6 +52,9 @@ class FixVector : public Fix {
int ncount; // # of values currently in growing vector or array int ncount; // # of values currently in growing vector or array
int ncountmax; // max # of values vector/array can hold int ncountmax; // max # of values vector/array can hold
int nmaxval; // maximum allowed number of values
int nindex; // start index of data, may wrap around
double *vector; double *vector;
double **array; double **array;
}; };