From beb5a30f6750b7592cef773aa47b1ed5520fa9a2 Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Wed, 30 Nov 2016 14:01:27 -0700
Subject: [PATCH] new compute global/atom command, also bug fix for descending
dump sorts
---
doc/src/Section_commands.txt | 1 +
doc/src/compute_chunk_atom.txt | 3 +-
doc/src/compute_global_atom.txt | 220 +++++++++++++
doc/src/lammps.book | 1 +
src/compute.cpp | 3 +-
src/compute_bond_local.cpp | 2 +-
src/compute_global_atom.cpp | 540 ++++++++++++++++++++++++++++++++
src/compute_global_atom.h | 143 +++++++++
src/fix_ave_time.cpp | 6 +-
9 files changed, 914 insertions(+), 5 deletions(-)
create mode 100644 doc/src/compute_global_atom.txt
create mode 100644 src/compute_global_atom.cpp
create mode 100644 src/compute_global_atom.h
diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt
index 4707d25c4f..7cf0b3df90 100644
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@@ -767,6 +767,7 @@ KOKKOS, o = USER-OMP, t = OPT.
"erotate/sphere"_compute_erotate_sphere.html,
"erotate/sphere/atom"_compute_erotate_sphere_atom.html,
"event/displace"_compute_event_displace.html,
+"global/atom"_compute_global_atom.html,
"group/group"_compute_group_group.html,
"gyration"_compute_gyration.html,
"gyration/chunk"_compute_gyration_chunk.html,
diff --git a/doc/src/compute_chunk_atom.txt b/doc/src/compute_chunk_atom.txt
index 60516fe421..e6083ad615 100644
--- a/doc/src/compute_chunk_atom.txt
+++ b/doc/src/compute_chunk_atom.txt
@@ -641,7 +641,8 @@ the restarted simulation begins.
[Related commands:]
-"fix ave/chunk"_fix_ave_chunk.html
+"fix ave/chunk"_fix_ave_chunk.html,
+"compute global/atom"_compute_global_atom.html
[Default:]
diff --git a/doc/src/compute_global_atom.txt b/doc/src/compute_global_atom.txt
new file mode 100644
index 0000000000..f62efcff2e
--- /dev/null
+++ b/doc/src/compute_global_atom.txt
@@ -0,0 +1,220 @@
+"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
+
+compute global/atom command :h3
+
+[Syntax:]
+
+compute ID group-ID style index input1 input2 ... :pre
+
+ID, group-ID are documented in "compute"_compute.html command :ulb,l
+global/atom = style name of this compute command :l
+index = c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_name :l
+ c_ID = per-atom vector calculated by a compute with ID
+ c_ID\[I\] = Ith column of per-atom array calculated by a compute with ID
+ f_ID = per-atom vector calculated by a fix with ID
+ f_ID\[I\] = Ith column of per-atom array calculated by a fix with ID
+ v_name = per-atom vector calculated by an atom-style variable with name :pre
+one or more inputs can be listed :l
+input = c_ID, c_ID\[N\], f_ID, f_ID\[N\], v_name :l
+ c_ID = global vector calculated by a compute with ID
+ c_ID\[I\] = Ith column of global array calculated by a compute with ID, I can include wildcard (see below)
+ f_ID = global vector calculated by a fix with ID
+ f_ID\[I\] = Ith column of global array calculated by a fix with ID, I can include wildcard (see below)
+ v_name = global vector calculated by a vector-style variable with name :pre
+:ule
+
+[Examples:]
+
+compute 1 all global/atom c_chunk c_com\[1\\] c_com\[2\\] c_com\[3\\]
+compute 1 all global/atom c_chunk c_com\[*\\] :pre
+
+[Description:]
+
+Define a calculation that assigns global values to each atom from
+vectors or arrays of global values. The specified {index} parameter
+is used to determine which global value is assigned to each atom.
+
+The {index} parameter must reference a per-atom vector or array from a
+"compute"_compute.html or "fix"_fix.html or the evaluation of an
+atom-style "variable"_variable.html. Each {input} value must
+reference a global vector or array from a "compute"_compute.html or
+"fix"_fix.html or the evaluation of an vector-style
+"variable"_variable.html. Details are given below.
+
+The {index} value for an atom is used as a index I (from 1 to N) into
+the vector associated with each of the input values. The Ith value
+from the input vector becomes one output value for that atom. If the
+atom is not in the specified group, or the index I < 1 or I > M, where
+M is the actual length of the input vector, then an output value of
+0.0 is assigned to the atom.
+
+An example of how this command is useful, is in the context of
+"chunks" which are static or dyanmic subsets of atoms. The "compute
+chunk/atom"_compute_chunk_atom.html command assigns unique chunk IDs
+to each atom. It's output can be used as the {index} parameter for
+this command. Various other computes with "chunk" in their style
+name, such as "compute com/chunk"_compute_com_chunk.html or "compute
+msd/chunk"_compute_msd_chunk.html, calculate properties for each
+chunk. The output of these commands are global vectors or arrays,
+with one or more values per chunk, and can be used as input values for
+this command. This command will then assign the global chunk value to
+each atom in the chunk, producing a per-atom vector or per-atom array
+as output. The per-atom values can then be output to a dump file or
+used by any command that uses per-atom values from a compute as input,
+as discussed in "Section 6.15"_Section_howto.html#howto_15.
+
+As a concrete example, these commands will calculate the displacement
+of each atom from the center-of-mass of the molecule it is in, and
+dump those values to a dump file. In this case, each molecule is a
+chunk.
+
+compute cc1 all chunk/atom molecule
+compute myChunk all com/chunk cc1
+compute prop all property/atom xu yu zu
+compute glob all global/atom c_cc1 c_myChunk\[*\]
+variable dx atom c_prop\[1\]-c_glob\[1\]
+variable dy atom c_prop\[2\]-c_glob\[2\]
+variable dz atom c_prop\[3\]-c_glob\[3\]
+variable dist atom sqrt(v_dx*v_dx+v_dy*v_dy+v_dz*v_dz)
+dump 1 all custom 100 tmp.dump id xu yu zu c_glob\[1\] c_glob\[2\] c_glob\[3\] &
+ v_dx v_dy v_dz v_dist
+dump_modify 1 sort id :pre
+
+You can add these commands to the bench/in.chain script to see how
+they work.
+
+:line
+
+Note that for input values from a compute or fix, the bracketed index
+I can be specified using a wildcard asterisk with the index to
+effectively specify multiple values. This takes the form "*" or "*n"
+or "n*" or "m*n". If N = the size of the vector (for {mode} = scalar)
+or the number of columns in the array (for {mode} = vector), then an
+asterisk with no numeric values means all indices from 1 to N. A
+leading asterisk means all indices from 1 to n (inclusive). A
+trailing asterisk means all indices from n to N (inclusive). A middle
+asterisk means all indices from m to n (inclusive).
+
+Using a wildcard is the same as if the individual columns of the array
+had been listed one by one. E.g. these 2 compute global/atom commands
+are equivalent, since the "compute com/chunk"_compute_com_chunk.html
+command creates a global array with 3 columns:
+
+compute cc1 all chunk/atom molecule
+compute com all com/chunk cc1
+compute 1 all global/atom c_cc1 c_com\[1\] c_com\[2\] c_com\[3\]
+compute 1 all global/atom c_cc1 c_com\[*\] :pre
+
+:line
+
+This section explains the {index} parameter. Note that it must
+reference per-atom values, as contrasted with the {input} values which
+must reference global values.
+
+Note that all of these options generate floating point values. When
+they are used as an index into the specified input vectors, they
+simple rounded down to convert the value to integer indices. The
+final values should range from 1 to N (inclusive), since they are used
+to access values from N-length vectors.
+
+If {index} begins with "c_", a compute ID must follow which has been
+previously defined in the input script. The compute must generate
+per-atom quantities. See the individual "compute"_compute.html doc
+page for details. If no bracketed integer is appended, the per-atom
+vector calculated by the compute is used. If a bracketed integer is
+appended, the Ith column of the per-atom array calculated by the
+compute is used. Users can also write code for their own compute
+styles and "add them to LAMMPS"_Section_modify.html. See the
+discussion above for how I can be specified with a wildcard asterisk
+to effectively specify multiple values.
+
+If {index} begins with "f_", a fix ID must follow which has been
+previously defined in the input script. The Fix must generate
+per-atom quantities. See the individual "fix"_fix.html doc page for
+details. Note that some fixes only produce their values on certain
+timesteps, which must be compatible with when compute global/atom
+references the values, else an error results. If no bracketed integer
+is appended, the per-atom vector calculated by the fix is used. If a
+bracketed integer is appended, the Ith column of the per-atom array
+calculated by the fix is used. Users can also write code for their
+own fix style and "add them to LAMMPS"_Section_modify.html. See the
+discussion above for how I can be specified with a wildcard asterisk
+to effectively specify multiple values.
+
+If {index} begins with "v_", a variable name must follow which has
+been previously defined in the input script. It must be an
+"atom-style variable"_variable.html. Atom-style variables can
+reference thermodynamic keywords and various per-atom attributes, or
+invoke other computes, fixes, or variables when they are evaluated, so
+this is a very general means of generating per-atom quantities to use
+as {index}.
+
+:line
+
+This section explains the kinds of {input} values that can be used.
+Note that inputs reference global values, as contrasted with the
+{index} parameter which must reference per-atom values.
+
+If a value begins with "c_", a compute ID must follow which has been
+previously defined in the input script. The compute must generate a
+global vector or array. See the individual "compute"_compute.html doc
+page for details. If no bracketed integer is appended, the vector
+calculated by the compute is used. If a bracketed integer is
+appended, the Ith column of the array calculated by the compute is
+used. Users can also write code for their own compute styles and "add
+them to LAMMPS"_Section_modify.html. See the discussion above for how
+I can be specified with a wildcard asterisk to effectively specify
+multiple values.
+
+If a value begins with "f_", a fix ID must follow which has been
+previously defined in the input script. The fix must generate a
+global vector or array. See the individual "fix"_fix.html doc page
+for details. Note that some fixes only produce their values on
+certain timesteps, which must be compatible with when compute
+global/atom references the values, else an error results. If no
+bracketed integer is appended, the vector calculated by the fix is
+used. If a bracketed integer is appended, the Ith column of the array
+calculated by the fix is used. Users can also write code for their
+own fix style and "add them to LAMMPS"_Section_modify.html. See the
+discussion above for how I can be specified with a wildcard asterisk
+to effectively specify multiple values.
+
+If a value begins with "v_", a variable name must follow which has
+been previously defined in the input script. It must be a
+"vector-style variable"_variable.html. Vector-style variables can
+reference thermodynamic keywords and various other attributes of
+atoms, or invoke other computes, fixes, or variables when they are
+evaluated, so this is a very general means of generating a vector of
+global quantities which the {index} parameter will reference for
+assignement of global values to atoms.
+
+:line
+
+[Output info:]
+
+If a single input is specified this compute produces a per-atom
+vector. If multiple inputs are specified, this compute produces a
+per-atom array values, where the number of columns is equal to the
+number of inputs specified. These values can be used by any command
+that uses per-atom vector or array values from a compute as input.
+See "Section 6.15"_Section_howto.html#howto_15 for an overview of
+LAMMPS output options.
+
+The per-atom vector or array values will be in whatever units the
+corresponsing input values are in.
+
+[Restrictions:] none
+
+[Related commands:]
+
+"compute"_compute.html, "fix"_fix.html, "variable"_variable.html,
+"compute chunk/atom"_compute_chunk_atom.html, "compute
+reduce"_compute_reduce.html
+
+[Default:] none
diff --git a/doc/src/lammps.book b/doc/src/lammps.book
index b83b2021c4..ec689adf15 100644
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@@ -311,6 +311,7 @@ compute_erotate_sphere.html
compute_erotate_sphere_atom.html
compute_event_displace.html
compute_fep.html
+compute_global_atom.html
compute_group_group.html
compute_gyration.html
compute_gyration_chunk.html
diff --git a/src/compute.cpp b/src/compute.cpp
index d306b0b34f..88ae70c067 100644
--- a/src/compute.cpp
+++ b/src/compute.cpp
@@ -40,7 +40,8 @@ int Compute::instance_total = 0;
Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp),
id(NULL), style(NULL),
- vector(NULL), array(NULL), vector_atom(NULL), array_atom(NULL), vector_local(NULL), array_local(NULL),
+ vector(NULL), array(NULL), vector_atom(NULL),
+ array_atom(NULL), vector_local(NULL), array_local(NULL),
tlist(NULL), vbiasall(NULL)
{
instance_me = instance_total++;
diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp
index b866be6991..424059aa86 100644
--- a/src/compute_bond_local.cpp
+++ b/src/compute_bond_local.cpp
@@ -70,7 +70,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
// set velflag if compute any quantities based on velocities
singleflag = 0;
- ghostvelflag = 0;
+ velflag = 0;
for (int i = 0; i < nvalues; i++) {
if (bstyle[i] == ENGPOT || bstyle[i] == FORCE) singleflag = 1;
if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS ||
diff --git a/src/compute_global_atom.cpp b/src/compute_global_atom.cpp
new file mode 100644
index 0000000000..b786043a34
--- /dev/null
+++ b/src/compute_global_atom.cpp
@@ -0,0 +1,540 @@
+/* ----------------------------------------------------------------------
+ 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 <string.h>
+#include <stdlib.h>
+#include "compute_global_atom.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "modify.h"
+#include "fix.h"
+#include "force.h"
+#include "comm.h"
+#include "group.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{COMPUTE,FIX,VARIABLE};
+enum{VECTOR,ARRAY};
+
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+#define INVOKED_PERATOM 8
+
+#define BIG 1.0e20
+
+/* ---------------------------------------------------------------------- */
+
+ComputeGlobalAtom::ComputeGlobalAtom(LAMMPS *lmp, int narg, char **arg) :
+ Compute(lmp, narg, arg),
+ idref(NULL), which(NULL), argindex(NULL), value2index(NULL), ids(NULL),
+ indices(NULL), varatom(NULL), vecglobal(NULL)
+{
+ if (narg < 5) error->all(FLERR,"Illegal compute global/atom command");
+
+ peratom_flag = 1;
+
+ // process index arg
+
+ int iarg = 3;
+
+ if (strncmp(arg[iarg],"c_",2) == 0 ||
+ strncmp(arg[iarg],"f_",2) == 0 ||
+ strncmp(arg[iarg],"v_",2) == 0) {
+ if (arg[iarg][0] == 'c') whichref = COMPUTE;
+ else if (arg[iarg][0] == 'f') whichref = FIX;
+ else if (arg[iarg][0] == 'v') whichref = VARIABLE;
+
+ int n = strlen(arg[iarg]);
+ char *suffix = new char[n];
+ strcpy(suffix,&arg[iarg][2]);
+
+ char *ptr = strchr(suffix,'[');
+ if (ptr) {
+ if (suffix[strlen(suffix)-1] != ']')
+ error->all(FLERR,"Illegal compute global/atom command");
+ indexref = atoi(ptr+1);
+ *ptr = '\0';
+ } else indexref = 0;
+
+ n = strlen(suffix) + 1;
+ idref = new char[n];
+ strcpy(idref,suffix);
+ delete [] suffix;
+ } else error->all(FLERR,"Illegal compute global/atom command");
+
+ iarg++;
+
+ // expand args if any have wildcard character "*"
+
+ int expand = 0;
+ char **earg;
+ int nargnew = input->expand_args(narg-iarg,&arg[iarg],1,earg);
+
+ if (earg != &arg[iarg]) expand = 1;
+ arg = earg;
+
+ // parse values until one isn't recognized
+
+ which = new int[nargnew];
+ argindex = new int[nargnew];
+ ids = new char*[nargnew];
+ value2index = new int[nargnew];
+ nvalues = 0;
+
+ iarg = 0;
+ while (iarg < nargnew) {
+ ids[nvalues] = NULL;
+ if (strncmp(arg[iarg],"c_",2) == 0 ||
+ strncmp(arg[iarg],"f_",2) == 0 ||
+ strncmp(arg[iarg],"v_",2) == 0) {
+ if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
+ else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
+ else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
+
+ int n = strlen(arg[iarg]);
+ char *suffix = new char[n];
+ strcpy(suffix,&arg[iarg][2]);
+
+ char *ptr = strchr(suffix,'[');
+ if (ptr) {
+ if (suffix[strlen(suffix)-1] != ']')
+ error->all(FLERR,"Illegal compute global/atom command");
+ argindex[nvalues] = atoi(ptr+1);
+ *ptr = '\0';
+ } else argindex[nvalues] = 0;
+
+ n = strlen(suffix) + 1;
+ ids[nvalues] = new char[n];
+ strcpy(ids[nvalues],suffix);
+ nvalues++;
+ delete [] suffix;
+
+ } else error->all(FLERR,"Illegal compute global/atom command");
+
+ iarg++;
+ }
+
+ // if wildcard expansion occurred, free earg memory from expand_args()
+
+ if (expand) {
+ for (int i = 0; i < nargnew; i++) delete [] earg[i];
+ memory->sfree(earg);
+ }
+
+ // setup and error check both index arg and values
+
+ if (whichref == COMPUTE) {
+ int icompute = modify->find_compute(idref);
+ if (icompute < 0)
+ error->all(FLERR,"Compute ID for compute global/atom does not exist");
+
+ if (!modify->compute[icompute]->peratom_flag)
+ error->all(FLERR,"Compute global/atom compute does not "
+ "calculate a per-atom vector or array");
+ if (indexref == 0 &&
+ modify->compute[icompute]->size_peratom_cols != 0)
+ error->all(FLERR,"Compute global/atom compute does not "
+ "calculate a per-atom vector");
+ if (indexref && modify->compute[icompute]->size_peratom_cols == 0)
+ error->all(FLERR,"Compute global/atom compute does not "
+ "calculate a per-atom array");
+ if (indexref && indexref > modify->compute[icompute]->size_peratom_cols)
+ error->all(FLERR,
+ "Compute global/atom compute array is accessed out-of-range");
+
+ } else if (whichref == FIX) {
+ int ifix = modify->find_fix(idref);
+ if (ifix < 0)
+ error->all(FLERR,"Fix ID for compute global/atom does not exist");
+ if (!modify->fix[ifix]->peratom_flag)
+ error->all(FLERR,"Compute global/atom fix does not "
+ "calculate a per-atom vector or array");
+ if (indexref == 0 &&
+ modify->fix[ifix]->size_peratom_cols != 0)
+ error->all(FLERR,"Compute global/atom fix does not "
+ "calculate a per-atom vector");
+ if (indexref && modify->fix[ifix]->size_peratom_cols == 0)
+ error->all(FLERR,"Compute global/atom fix does not "
+ "calculate a per-atom array");
+ if (indexref && indexref > modify->fix[ifix]->size_peratom_cols)
+ error->all(FLERR,
+ "Compute global/atom fix array is accessed out-of-range");
+
+ } else if (whichref == VARIABLE) {
+ int ivariable = input->variable->find(idref);
+ if (ivariable < 0)
+ error->all(FLERR,"Variable name for compute global/atom does not exist");
+ if (input->variable->atomstyle(ivariable) == 0)
+ error->all(FLERR,"Compute global/atom variable is not "
+ "atom-style variable");
+ }
+
+ for (int i = 0; i < nvalues; i++) {
+ if (which[i] == COMPUTE) {
+ int icompute = modify->find_compute(ids[i]);
+ if (icompute < 0)
+ error->all(FLERR,"Compute ID for compute global/atom does not exist");
+ if (argindex[i] == 0) {
+ if (!modify->compute[icompute]->vector_flag)
+ error->all(FLERR,"Compute global/atom compute does not "
+ "calculate a global vector");
+ } else {
+ if (!modify->compute[icompute]->array_flag)
+ error->all(FLERR,"Compute global/atom compute does not "
+ "calculate a global array");
+ if (argindex[i] > modify->compute[icompute]->size_array_cols)
+ error->all(FLERR,"Compute global/atom compute array is "
+ "accessed out-of-range");
+ }
+
+ } else if (which[i] == FIX) {
+ int ifix = modify->find_fix(ids[i]);
+ if (ifix < 0)
+ error->all(FLERR,"Fix ID for compute global/atom does not exist");
+ if (argindex[i] == 0) {
+ if (!modify->fix[ifix]->vector_flag)
+ error->all(FLERR,"Compute global/atom fix does not "
+ "calculate a global vector");
+ } else {
+ if (!modify->fix[ifix]->array_flag)
+ error->all(FLERR,"Compute global/atom fix does not "
+ "calculate a global array");
+ if (argindex[i] > modify->fix[ifix]->size_array_cols)
+ error->all(FLERR,"Compute global/atom fix array is "
+ "accessed out-of-range");
+ }
+
+ } else if (which[i] == VARIABLE) {
+ int ivariable = input->variable->find(ids[i]);
+ if (ivariable < 0)
+ error->all(FLERR,"Variable name for compute global/atom "
+ "does not exist");
+ if (input->variable->vectorstyle(ivariable) == 0)
+ error->all(FLERR,"Compute global/atom variable is not "
+ "vector-style variable");
+ }
+ }
+
+ // this compute produces either a peratom vector or array
+
+ if (nvalues == 1) size_peratom_cols = 0;
+ else size_peratom_cols = nvalues;
+
+ nmax = maxvector = 0;
+ indices = NULL;
+ varatom = NULL;
+ vecglobal = NULL;
+
+ vector_atom = NULL;
+ array_atom = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeGlobalAtom::~ComputeGlobalAtom()
+{
+ delete [] idref;
+
+ delete [] which;
+ delete [] argindex;
+ for (int m = 0; m < nvalues; m++) delete [] ids[m];
+ delete [] ids;
+ delete [] value2index;
+
+ memory->destroy(indices);
+ memory->destroy(varatom);
+ memory->destroy(vecglobal);
+ memory->destroy(vector_atom);
+ memory->destroy(array_atom);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeGlobalAtom::init()
+{
+ // set indices of all computes,fixes,variables
+
+ if (whichref == COMPUTE) {
+ int icompute = modify->find_compute(idref);
+ if (icompute < 0)
+ error->all(FLERR,"Compute ID for compute global/atom does not exist");
+ ref2index = icompute;
+ } else if (whichref == FIX) {
+ int ifix = modify->find_fix(idref);
+ if (ifix < 0)
+ error->all(FLERR,"Fix ID for compute global/atom does not exist");
+ ref2index = ifix;
+ } else if (whichref == VARIABLE) {
+ int ivariable = input->variable->find(idref);
+ if (ivariable < 0)
+ error->all(FLERR,"Variable name for compute global/atom does not exist");
+ ref2index = ivariable;
+ }
+
+ for (int m = 0; m < nvalues; m++) {
+ if (which[m] == COMPUTE) {
+ int icompute = modify->find_compute(ids[m]);
+ if (icompute < 0)
+ error->all(FLERR,"Compute ID for compute global/atom does not exist");
+ value2index[m] = icompute;
+
+ } else if (which[m] == FIX) {
+ int ifix = modify->find_fix(ids[m]);
+ if (ifix < 0)
+ error->all(FLERR,"Fix ID for compute global/atom does not exist");
+ value2index[m] = ifix;
+
+ } else if (which[m] == VARIABLE) {
+ int ivariable = input->variable->find(ids[m]);
+ if (ivariable < 0)
+ error->all(FLERR,"Variable name for compute global/atom "
+ "does not exist");
+ value2index[m] = ivariable;
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeGlobalAtom::compute_peratom()
+{
+ int i,j;
+
+ invoked_peratom = update->ntimestep;
+
+ // grow indices and output vector or array if necessary
+
+ if (atom->nmax > nmax) {
+ nmax = atom->nmax;
+ memory->destroy(indices);
+ memory->create(indices,nmax,"global/atom:indices");
+ if (nvalues == 1) {
+ memory->destroy(vector_atom);
+ memory->create(vector_atom,nmax,"global/atom:vector_atom");
+ } else {
+ memory->destroy(array_atom);
+ memory->create(array_atom,nmax,nvalues,"global/atom:array_atom");
+ }
+ }
+
+ // setup current peratom indices
+ // integer indices are rounded down from double values
+ // indices are decremented from 1 to N -> 0 to N-1
+
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ if (whichref == COMPUTE) {
+ Compute *compute = modify->compute[ref2index];
+
+ if (!(compute->invoked_flag & INVOKED_PERATOM)) {
+ compute->compute_peratom();
+ compute->invoked_flag |= INVOKED_PERATOM;
+ }
+
+ if (indexref == 0) {
+ double *compute_vector = compute->vector_atom;
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ indices[i] = static_cast<int> (compute_vector[i]) - 1;
+ } else {
+ double **compute_array = compute->array_atom;
+ int im1 = indexref - 1;
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ indices[i] = static_cast<int> (compute_array[i][im1]) - 1;
+ }
+
+ } else if (whichref == FIX) {
+ if (update->ntimestep % modify->fix[ref2index]->peratom_freq)
+ error->all(FLERR,"Fix used in compute global/atom not "
+ "computed at compatible time");
+ Fix *fix = modify->fix[ref2index];
+
+ if (indexref == 0) {
+ double *fix_vector = fix->vector_atom;
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ indices[i] = static_cast<int> (fix_vector[i]) - 1;
+ } else {
+ double **fix_array = fix->array_atom;
+ int im1 = indexref - 1;
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ indices[i] = static_cast<int> (fix_array[i][im1]) - 1;
+ }
+
+ } else if (whichref == VARIABLE) {
+ if (atom->nmax > nmax) {
+ nmax = atom->nmax;
+ memory->destroy(varatom);
+ memory->create(varatom,nmax,"global/atom:varatom");
+ }
+
+ input->variable->compute_atom(ref2index,igroup,varatom,1,0);
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ indices[i] = static_cast<int> (varatom[i]) - 1;
+ }
+
+ // loop over values to fill output vector or array
+
+ for (int m = 0; m < nvalues; m++) {
+
+ // output = vector
+
+ if (argindex[m] == 0) {
+ int vmax;
+ double *source;
+
+ if (which[m] == COMPUTE) {
+ Compute *compute = modify->compute[value2index[m]];
+
+ if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+ compute->compute_vector();
+ compute->invoked_flag |= INVOKED_VECTOR;
+ }
+
+ source = compute->vector;
+ vmax = compute->size_vector;
+
+ } else if (which[m] == FIX) {
+ if (update->ntimestep % modify->fix[value2index[m]]->peratom_freq)
+ error->all(FLERR,"Fix used in compute global/atom not "
+ "computed at compatible time");
+ Fix *fix = modify->fix[value2index[m]];
+ vmax = fix->size_vector;
+
+ if (vmax > maxvector) {
+ maxvector = vmax;
+ memory->destroy(vecglobal);
+ memory->create(vecglobal,maxvector,"global/atom:vecglobal");
+ }
+
+ for (i = 0; i < vmax; i++)
+ vecglobal[i] = fix->compute_vector(i);
+
+ source = vecglobal;
+
+ } else if (which[m] == VARIABLE) {
+ vmax = input->variable->compute_vector(value2index[m],&source);
+ }
+
+ if (nvalues == 1) {
+ for (i = 0; i < nlocal; i++) {
+ vector_atom[i] = 0.0;
+ if (mask[i] & groupbit) {
+ j = indices[i];
+ if (j >= 0 && j < vmax) vector_atom[i] = source[j];
+ }
+ }
+ } else {
+ for (i = 0; i < nlocal; i++) {
+ array_atom[i][m] = 0.0;
+ if (mask[i] & groupbit) {
+ j = indices[i];
+ if (j >= 0 && j < vmax) array_atom[i][m] = source[j];
+ }
+ }
+ }
+
+ // output = array
+
+ } else {
+ int vmax;
+ double *source;
+ int col = argindex[m] - 1;
+
+ if (which[m] == COMPUTE) {
+ Compute *compute = modify->compute[value2index[m]];
+
+ if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+ compute->compute_array();
+ compute->invoked_flag |= INVOKED_ARRAY;
+ }
+
+ double **compute_array = compute->array;
+ vmax = compute->size_array_rows;
+
+ if (vmax > maxvector) {
+ maxvector = vmax;
+ memory->destroy(vecglobal);
+ memory->create(vecglobal,maxvector,"global/atom:vecglobal");
+ }
+
+ for (i = 0; i < vmax; i++)
+ vecglobal[i] = compute_array[i][col];
+
+ source = vecglobal;
+
+ } else if (which[m] == FIX) {
+ if (update->ntimestep % modify->fix[value2index[m]]->peratom_freq)
+ error->all(FLERR,"Fix used in compute global/atom not "
+ "computed at compatible time");
+ Fix *fix = modify->fix[value2index[m]];
+ vmax = fix->size_array_rows;
+
+ if (vmax > maxvector) {
+ maxvector = vmax;
+ memory->destroy(vecglobal);
+ memory->create(vecglobal,maxvector,"global/atom:vecglobal");
+ }
+
+ for (i = 0; i < vmax; i++)
+ vecglobal[i] = fix->compute_array(i,col);
+
+ source = vecglobal;
+
+ } else if (which[m] == VARIABLE) {
+ vmax = input->variable->compute_vector(value2index[m],&source);
+ }
+
+ if (nvalues == 1) {
+ for (i = 0; i < nlocal; i++) {
+ vector_atom[i] = 0.0;
+ if (mask[i] & groupbit) {
+ j = indices[i];
+ if (j >= 0 && j < vmax) vector_atom[i] = source[j];
+ }
+ }
+ } else {
+ for (i = 0; i < nlocal; i++) {
+ array_atom[i][m] = 0.0;
+ if (mask[i] & groupbit) {
+ j = indices[i];
+ if (j >= 0 && j < vmax) array_atom[i][m] = source[j];
+ }
+ }
+ }
+ }
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double ComputeGlobalAtom::memory_usage()
+{
+ double bytes = nmax*nvalues * sizeof(double);
+ bytes += nmax * sizeof(int); // indices
+ if (varatom) bytes += nmax * sizeof(double); // varatom
+ bytes += maxvector * sizeof(double); // vecglobal
+ return bytes;
+}
diff --git a/src/compute_global_atom.h b/src/compute_global_atom.h
new file mode 100644
index 0000000000..439720f712
--- /dev/null
+++ b/src/compute_global_atom.h
@@ -0,0 +1,143 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+
+ComputeStyle(global/atom,ComputeGlobalAtom)
+
+#else
+
+#ifndef LMP_COMPUTE_GLOBAL_ATOM_H
+#define LMP_COMPUTE_GLOBAL_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeGlobalAtom : public Compute {
+ public:
+ ComputeGlobalAtom(class LAMMPS *, int, char **);
+ virtual ~ComputeGlobalAtom();
+ void init();
+ void compute_peratom();
+ double memory_usage();
+
+ protected:
+ int whichref,indexref,ref2index;
+ char *idref;
+
+ int nvalues;
+ int *which,*argindex,*value2index;
+ char **ids;
+
+ int nmax,maxvector;
+ int *indices;
+ double *varatom;
+ double *vecglobal;
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory. Check the input script syntax and compare to the
+documentation for the command. You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Region ID for compute reduce/region does not exist
+
+Self-explanatory.
+
+E: Compute reduce replace requires min or max mode
+
+Self-explanatory.
+
+E: Invalid replace values in compute reduce
+
+Self-explanatory.
+
+E: Compute ID for compute reduce does not exist
+
+Self-explanatory.
+
+E: Compute reduce compute does not calculate a per-atom vector
+
+Self-explanatory.
+
+E: Compute reduce compute does not calculate a per-atom array
+
+Self-explanatory.
+
+E: Compute reduce compute array is accessed out-of-range
+
+An index for the array is out of bounds.
+
+E: Compute reduce compute does not calculate a local vector
+
+Self-explanatory.
+
+E: Compute reduce compute does not calculate a local array
+
+Self-explanatory.
+
+E: Compute reduce compute calculates global values
+
+A compute that calculates peratom or local values is required.
+
+E: Fix ID for compute reduce does not exist
+
+Self-explanatory.
+
+E: Compute reduce fix does not calculate a per-atom vector
+
+Self-explanatory.
+
+E: Compute reduce fix does not calculate a per-atom array
+
+Self-explanatory.
+
+E: Compute reduce fix array is accessed out-of-range
+
+An index for the array is out of bounds.
+
+E: Compute reduce fix does not calculate a local vector
+
+Self-explanatory.
+
+E: Compute reduce fix does not calculate a local array
+
+Self-explanatory.
+
+E: Compute reduce fix calculates global values
+
+A fix that calculates peratom or local values is required.
+
+E: Variable name for compute reduce does not exist
+
+Self-explanatory.
+
+E: Compute reduce variable is not atom-style variable
+
+Self-explanatory.
+
+E: Fix used in compute reduce not computed at compatible time
+
+Fixes generate their values on specific timesteps. Compute reduce is
+requesting a value on a non-allowed timestep.
+
+*/
diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp
index ff47af03a8..56a2458aad 100644
--- a/src/fix_ave_time.cpp
+++ b/src/fix_ave_time.cpp
@@ -43,8 +43,10 @@ enum{SCALAR,VECTOR};
FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
- nvalues(0), which(NULL), argindex(NULL), value2index(NULL), offcol(NULL), varlen(NULL), ids(NULL),
- fp(NULL), offlist(NULL), format(NULL), format_user(NULL), vector(NULL), vector_total(NULL), vector_list(NULL),
+ nvalues(0), which(NULL), argindex(NULL), value2index(NULL),
+ offcol(NULL), varlen(NULL), ids(NULL),
+ fp(NULL), offlist(NULL), format(NULL), format_user(NULL),
+ vector(NULL), vector_total(NULL), vector_list(NULL),
column(NULL), array(NULL), array_total(NULL), array_list(NULL)
{
if (narg < 7) error->all(FLERR,"Illegal fix ave/time command");
--
GitLab