diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp
index 5920002dc1c332314eed17c16047ff70eebbbffd..79c55fc24c10a1c1149c254091989d19167d38f2 100755
--- a/src/ASPHERE/compute_temp_asphere.cpp
+++ b/src/ASPHERE/compute_temp_asphere.cpp
@@ -31,8 +31,6 @@
using namespace LAMMPS_NS;
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
-
/* ---------------------------------------------------------------------- */
ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
diff --git a/src/compute.h b/src/compute.h
index bf3aae4708a4facc54870e578d0e2cdc17fac849..83a6882b57ffeac0c15d387086beea254fbc5358 100644
--- a/src/compute.h
+++ b/src/compute.h
@@ -25,7 +25,7 @@ class Compute : protected Pointers {
double scalar; // computed global scalar
double *vector; // computed global vector
- double *array; // computed global array
+ double **array; // computed global array
double *vector_atom; // computed per-atom vector
double **array_atom; // computed per-atom array
double *vector_local; // computed local vector
@@ -48,7 +48,7 @@ class Compute : protected Pointers {
int extscalar; // 0/1 if global scalar is intensive/extensive
int extvector; // 0/1/-1 if global vector is all int/ext/extlist
int *extlist; // list of 0/1 int/ext for each vec component
- int extarray; // 0/1 if global array is intensive/extensive
+ int extarray; // 0/1 if global array is all intensive/extensive
int tempflag; // 1 if Compute can be used as temperature
// must have both compute_scalar, compute_vector
diff --git a/src/compute_heat_flux.cpp b/src/compute_heat_flux.cpp
index 738509903cf78c41ccdb4733a811eca289d04855..e038020524b309b21a95ba071c4b14881fcf4331 100644
--- a/src/compute_heat_flux.cpp
+++ b/src/compute_heat_flux.cpp
@@ -35,7 +35,7 @@
using namespace LAMMPS_NS;
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+#define INVOKED_PERATOM 8
/* ---------------------------------------------------------------------- */
@@ -67,6 +67,7 @@ ComputeHeatFlux::ComputeHeatFlux(LAMMPS *lmp, int narg, char **arg) :
ComputeHeatFlux::~ComputeHeatFlux()
{
+ delete [] id_atomPE;
delete [] vector;
}
diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp
index a827198b1a42b15cb9e7b5d4450d608675d6d381..e8c2221df2eb3d52141b52111b68ccb7952055a0 100644
--- a/src/compute_pressure.cpp
+++ b/src/compute_pressure.cpp
@@ -31,8 +31,6 @@
using namespace LAMMPS_NS;
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
-
/* ---------------------------------------------------------------------- */
ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp
index cdd159b6c9d73f132f81c1cc043010efbc281185..5cda51396b6442958c4e154a8259a67880047b57 100644
--- a/src/compute_reduce.cpp
+++ b/src/compute_reduce.cpp
@@ -21,6 +21,7 @@
#include "fix.h"
#include "force.h"
#include "comm.h"
+#include "group.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
@@ -28,9 +29,14 @@
using namespace LAMMPS_NS;
-enum{SUM,MINN,MAXX};
+enum{SUM,MINN,MAXX,AVE};
enum{X,V,F,COMPUTE,FIX,VARIABLE};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+enum{GLOBAL,PERATOM,LOCAL};
+
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+#define INVOKED_PERATOM 8
+#define INVOKED_LOCAL 16
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
@@ -55,6 +61,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg],"sum") == 0) mode = SUM;
else if (strcmp(arg[iarg],"min") == 0) mode = MINN;
else if (strcmp(arg[iarg],"max") == 0) mode = MAXX;
+ else if (strcmp(arg[iarg],"ave") == 0) mode = AVE;
else error->all("Illegal compute reduce command");
iarg++;
@@ -62,6 +69,7 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
which = new int[narg-4];
argindex = new int[narg-4];
+ flavor = new int[narg-4];
ids = new char*[narg-4];
value2index = new int[narg-4];
nvalues = 0;
@@ -136,26 +144,88 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
int icompute = modify->find_compute(ids[i]);
if (icompute < 0)
error->all("Compute ID for compute reduce does not exist");
- if (modify->compute[icompute]->peratom_flag == 0)
- error->all("Compute reduce compute does not "
- "calculate per-atom values");
- if (argindex[i] == 0 &&
- modify->compute[icompute]->size_peratom_cols != 0)
- error->all("Compute reduce compute does not "
- "calculate a per-atom vector");
- if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
- error->all("Compute reduce compute does not "
- "calculate a per-atom array");
+ if (modify->compute[icompute]->vector_flag ||
+ modify->compute[icompute]->array_flag) {
+ flavor[i] = GLOBAL;
+ if (argindex[i] == 0 &&
+ modify->compute[icompute]->vector_flag == 0)
+ error->all("Compute reduce compute does not "
+ "calculate a global vector");
+ if (argindex[i] && modify->compute[icompute]->array_flag == 0)
+ error->all("Compute reduce compute does not "
+ "calculate a global array");
+ if (argindex[i] &&
+ argindex[i] > modify->compute[icompute]->size_array_cols)
+ error->all("Compute reduce compute array is accessed out-of-range");
+ } else if (modify->compute[icompute]->peratom_flag) {
+ flavor[i] = PERATOM;
+ if (argindex[i] == 0 &&
+ modify->compute[icompute]->size_peratom_cols != 0)
+ error->all("Compute reduce compute does not "
+ "calculate a per-atom vector");
+ if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
+ error->all("Compute reduce compute does not "
+ "calculate a per-atom array");
+ if (argindex[i] &&
+ argindex[i] > modify->compute[icompute]->size_peratom_cols)
+ error->all("Compute reduce compute array is accessed out-of-range");
+ } else if (modify->compute[icompute]->local_flag) {
+ flavor[i] = LOCAL;
+ if (argindex[i] == 0 &&
+ modify->compute[icompute]->size_local_cols != 0)
+ error->all("Compute reduce compute does not "
+ "calculate a local vector");
+ if (argindex[i] && modify->compute[icompute]->size_local_cols == 0)
+ error->all("Compute reduce compute does not "
+ "calculate a local array");
+ if (argindex[i] &&
+ argindex[i] > modify->compute[icompute]->size_local_cols)
+ error->all("Compute reduce compute array is accessed out-of-range");
+ }
+
} else if (which[i] == FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all("Fix ID for compute reduce does not exist");
- if (modify->fix[ifix]->peratom_flag == 0)
- error->all("Compute reduce fix does not calculate per-atom values");
- if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
- error->all("Compute reduce fix does not calculate a per-atom vector");
- if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
- error->all("Compute reduce fix does not calculate a per-atom array");
+ if (modify->fix[ifix]->vector_flag ||
+ modify->fix[ifix]->array_flag) {
+ flavor[i] = GLOBAL;
+ if (argindex[i] == 0 &&
+ modify->fix[ifix]->vector_flag == 0)
+ error->all("Compute reduce fix does not "
+ "calculate a global vector");
+ if (argindex[i] && modify->fix[ifix]->array_flag == 0)
+ error->all("Compute reduce fix does not "
+ "calculate a global array");
+ if (argindex[i] &&
+ argindex[i] > modify->fix[ifix]->size_array_cols)
+ error->all("Compute reduce fix array is accessed out-of-range");
+ } else if (modify->fix[ifix]->peratom_flag) {
+ flavor[i] = PERATOM;
+ if (argindex[i] == 0 &&
+ modify->fix[ifix]->size_peratom_cols != 0)
+ error->all("Compute reduce fix does not "
+ "calculate a per-atom vector");
+ if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
+ error->all("Compute reduce fix does not "
+ "calculate a per-atom array");
+ if (argindex[i] &&
+ argindex[i] > modify->fix[ifix]->size_peratom_cols)
+ error->all("Compute reduce fix array is accessed out-of-range");
+ } else if (modify->fix[ifix]->local_flag) {
+ flavor[i] = LOCAL;
+ if (argindex[i] == 0 &&
+ modify->fix[ifix]->size_local_cols != 0)
+ error->all("Compute reduce fix does not "
+ "calculate a local vector");
+ if (argindex[i] && modify->fix[ifix]->size_local_cols == 0)
+ error->all("Compute reduce fix does not "
+ "calculate a local array");
+ if (argindex[i] &&
+ argindex[i] > modify->fix[ifix]->size_local_cols)
+ error->all("Compute reduce fix array is accessed out-of-range");
+ }
+
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
@@ -192,6 +262,7 @@ ComputeReduce::~ComputeReduce()
{
delete [] which;
delete [] argindex;
+ delete [] flavor;
for (int m = 0; m < nvalues; m++) delete [] ids[m];
delete [] ids;
delete [] value2index;
@@ -238,12 +309,17 @@ double ComputeReduce::compute_scalar()
invoked_scalar = update->ntimestep;
double one = compute_one(0);
+
if (mode == SUM)
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
else if (mode == MINN)
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MIN,world);
else if (mode == MAXX)
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MAX,world);
+ else if (mode == AVE) {
+ MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
+ scalar /= count(0);
+ }
return scalar;
}
@@ -254,12 +330,17 @@ void ComputeReduce::compute_vector()
invoked_vector = update->ntimestep;
for (int m = 0; m < nvalues; m++) onevec[m] = compute_one(m);
+
if (mode == SUM)
MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
else if (mode == MINN)
MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_MIN,world);
else if (mode == MAXX)
MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_MAX,world);
+ else if (mode == AVE) {
+ MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_MAX,world);
+ for (int m = 0; m < nvalues; m++) vector[m] /= count(m);
+ }
}
/* ---------------------------------------------------------------------- */
@@ -269,15 +350,15 @@ double ComputeReduce::compute_one(int m)
int i;
// invoke the appropriate attribute,compute,fix,variable
- // compute scalar quantity by summing over atom scalars
- // only include atoms in group
-
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
+ // compute scalar quantity by summing over atom properties
+ // only include atoms in group for atom properties and per-atom quantities
int n = value2index[m];
int j = argindex[m];
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
double one;
if (mode == SUM) one = 0.0;
else if (mode == MINN) one = BIG;
@@ -300,37 +381,114 @@ double ComputeReduce::compute_one(int m)
} else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[n];
- if (!(compute->invoked_flag & INVOKED_PERATOM)) {
- compute->compute_peratom();
- compute->invoked_flag |= INVOKED_PERATOM;
- }
-
- if (j == 0) {
- double *compute_vector = compute->vector_atom;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) combine(one,compute_vector[i]);
- } else {
- double **compute_array = compute->array_atom;
- int jm1 = j - 1;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) combine(one,compute_array[i][jm1]);
+
+ if (flavor[m] == GLOBAL) {
+ if (j == 0) {
+ if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+ compute->compute_vector();
+ compute->invoked_flag |= INVOKED_VECTOR;
+ }
+ double *compute_vector = compute->vector;
+ int n = compute->size_vector;
+ for (i = 0; i < n; i++)
+ combine(one,compute_vector[i]);
+ } else {
+ if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+ compute->compute_array();
+ compute->invoked_flag |= INVOKED_ARRAY;
+ }
+ double **compute_array = compute->array;
+ int n = compute->size_array_rows;
+ int jm1 = j - 1;
+ for (i = 0; i < n; i++)
+ combine(one,compute_array[i][jm1]);
+ }
+
+ } else if (flavor[m] == PERATOM) {
+ if (!(compute->invoked_flag & INVOKED_PERATOM)) {
+ compute->compute_peratom();
+ compute->invoked_flag |= INVOKED_PERATOM;
+ }
+
+ if (j == 0) {
+ double *compute_vector = compute->vector_atom;
+ int n = nlocal;
+ for (i = 0; i < n; i++)
+ if (mask[i] & groupbit) combine(one,compute_vector[i]);
+ } else {
+ double **compute_array = compute->array_atom;
+ int n = nlocal;
+ int jm1 = j - 1;
+ for (i = 0; i < n; i++)
+ if (mask[i] & groupbit) combine(one,compute_array[i][jm1]);
+ }
+
+ } else if (flavor[m] == LOCAL) {
+ if (!(compute->invoked_flag & INVOKED_LOCAL)) {
+ compute->compute_local();
+ compute->invoked_flag |= INVOKED_LOCAL;
+ }
+
+ if (j == 0) {
+ double *compute_vector = compute->vector_local;
+ int n = compute->size_local_rows;
+ for (i = 0; i < n; i++)
+ combine(one,compute_vector[i]);
+ } else {
+ double **compute_array = compute->array_local;
+ int n = compute->size_local_rows;
+ int jm1 = j - 1;
+ for (i = 0; i < n; i++)
+ combine(one,compute_array[i][jm1]);
+ }
}
- // access fix fields, check if frequency is a match
+ // access fix fields, check if fix frequency is a match
} else if (which[m] == FIX) {
if (update->ntimestep % modify->fix[n]->peratom_freq)
error->all("Fix used in compute reduce not computed at compatible time");
-
- if (j == 0) {
- double *fix_vector = modify->fix[n]->vector_atom;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) combine(one,fix_vector[i]);
- } else {
- double **fix_array = modify->fix[n]->array_atom;
- int jm1 = j - 1;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) combine(one,fix_array[i][jm1]);
+ Fix *fix = modify->fix[n];
+
+ if (flavor[m] == GLOBAL) {
+ if (j == 0) {
+ int n = fix->size_vector;
+ for (i = 0; i < n; i++)
+ combine(one,fix->compute_vector(i));
+ } else {
+ int n = fix->size_array_rows;
+ int jm1 = j - 1;
+ for (i = 0; i < nlocal; i++)
+ combine(one,fix->compute_array(i,jm1));
+ }
+
+ } else if (flavor[m] == PERATOM) {
+ if (j == 0) {
+ double *fix_vector = fix->vector_atom;
+ int n = nlocal;
+ for (i = 0; i < n; i++)
+ if (mask[i] & groupbit) combine(one,fix_vector[i]);
+ } else {
+ double **fix_array = fix->array_atom;
+ int n = nlocal;
+ int jm1 = j - 1;
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) combine(one,fix_array[i][jm1]);
+ }
+
+ } else if (flavor[m] == LOCAL) {
+ if (j == 0) {
+ double *fix_vector = fix->vector_local;
+ int n = fix->size_local_rows;
+ for (i = 0; i < n; i++)
+ combine(one,fix_vector[i]);
+ } else {
+ double **fix_array = fix->array_local;
+ int n = fix->size_local_rows;
+ int jm1 = j - 1;
+ for (i = 0; i < n; i++)
+ combine(one,fix_array[i][jm1]);
+ }
}
// evaluate atom-style variable
@@ -351,13 +509,57 @@ double ComputeReduce::compute_one(int m)
return one;
}
+/* ---------------------------------------------------------------------- */
+
+double ComputeReduce::count(int m)
+{
+ int n = value2index[m];
+ int j = argindex[m];
+
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ if (which[m] == X || which[m] == V || which[m] == F)
+ return group->count(igroup);
+ else if (which[m] == COMPUTE) {
+ Compute *compute = modify->compute[n];
+ if (flavor[m] == GLOBAL) {
+ if (j == 0) return(1.0*compute->size_vector);
+ else return(1.0*compute->size_array_rows);
+ } else if (flavor[m] == PERATOM) {
+ return group->count(igroup);
+ } else if (flavor[m] == LOCAL) {
+ double ncount = compute->size_local_rows;
+ double ncountall;
+ MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
+ return ncountall;
+ }
+ } else if (which[m] == FIX) {
+ Fix *fix = modify->fix[n];
+ if (flavor[m] == GLOBAL) {
+ if (j == 0) return(1.0*fix->size_vector);
+ else return(1.0*fix->size_array_rows);
+ } else if (flavor[m] == PERATOM) {
+ return group->count(igroup);
+ } else if (flavor[m] == LOCAL) {
+ double ncount = fix->size_local_rows;
+ double ncountall;
+ MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
+ return ncountall;
+ }
+ } else if (which[m] == VARIABLE)
+ return group->count(igroup);
+
+ return 0.0;
+}
+
/* ----------------------------------------------------------------------
combine two values according to reduction mode
------------------------------------------------------------------------- */
void ComputeReduce::combine(double &one, double two)
{
- if (mode == SUM) one += two;
+ if (mode == SUM || mode == AVE) one += two;
else if (mode == MINN) one = MIN(one,two);
else if (mode == MAXX) one = MAX(one,two);
}
diff --git a/src/compute_reduce.h b/src/compute_reduce.h
index 82fe48ffb6d7cc4db7489139f325ed08e0eb7435..633cd6ab6985a6d1e8cc382e8883a1061111a6cf 100644
--- a/src/compute_reduce.h
+++ b/src/compute_reduce.h
@@ -29,7 +29,7 @@ class ComputeReduce : public Compute {
protected:
int mode,nvalues,iregion;
- int *which,*argindex,*value2index;
+ int *which,*argindex,*flavor,*value2index;
char **ids;
double *onevec;
@@ -37,6 +37,7 @@ class ComputeReduce : public Compute {
double *varatom;
virtual double compute_one(int);
+ virtual double count(int);
void combine(double &, double);
};
diff --git a/src/compute_reduce_region.cpp b/src/compute_reduce_region.cpp
index 9e46aab4fefdd89f04850d19c09794284db2fb79..3976714e87f70eda1b66860f22aafbeaa99e46a6 100644
--- a/src/compute_reduce_region.cpp
+++ b/src/compute_reduce_region.cpp
@@ -18,6 +18,7 @@
#include "update.h"
#include "modify.h"
#include "domain.h"
+#include "group.h"
#include "region.h"
#include "fix.h"
#include "input.h"
@@ -29,7 +30,12 @@ using namespace LAMMPS_NS;
enum{SUM,MINN,MAXX};
enum{X,V,F,COMPUTE,FIX,VARIABLE};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+enum{GLOBAL,PERATOM,LOCAL};
+
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+#define INVOKED_PERATOM 8
+#define INVOKED_LOCAL 16
#define BIG 1.0e20
@@ -81,41 +87,118 @@ double ComputeReduceRegion::compute_one(int m)
} else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[n];
- if (!(compute->invoked_flag & INVOKED_PERATOM)) {
- compute->compute_peratom();
- compute->invoked_flag |= INVOKED_PERATOM;
- }
+
+ if (flavor[m] == GLOBAL) {
+ if (j == 0) {
+ if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+ compute->compute_vector();
+ compute->invoked_flag |= INVOKED_VECTOR;
+ }
+ double *compute_vector = compute->vector;
+ int n = compute->size_vector;
+ for (i = 0; i < n; i++)
+ combine(one,compute_vector[i]);
+ } else {
+ if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+ compute->compute_array();
+ compute->invoked_flag |= INVOKED_ARRAY;
+ }
+ double **compute_array = compute->array;
+ int n = compute->size_array_rows;
+ int jm1 = j - 1;
+ for (i = 0; i < n; i++)
+ combine(one,compute_array[i][jm1]);
+ }
+
+ } else if (flavor[m] == PERATOM) {
+ if (!(compute->invoked_flag & INVOKED_PERATOM)) {
+ compute->compute_peratom();
+ compute->invoked_flag |= INVOKED_PERATOM;
+ }
+
+ if (j == 0) {
+ double *compute_vector = compute->vector_atom;
+ int n = nlocal;
+ for (i = 0; i < n; i++)
+ if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+ combine(one,compute_vector[i]);
+ } else {
+ double **compute_array = compute->array_atom;
+ int n = nlocal;
+ int jm1 = j - 1;
+ for (i = 0; i < n; i++)
+ if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+ combine(one,compute_array[i][jm1]);
+ }
- if (j == 0) {
- double *compute_vector = compute->vector_atom;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+ } else if (flavor[m] == LOCAL) {
+ if (!(compute->invoked_flag & INVOKED_LOCAL)) {
+ compute->compute_local();
+ compute->invoked_flag |= INVOKED_LOCAL;
+ }
+
+ if (j == 0) {
+ double *compute_vector = compute->vector_local;
+ int n = compute->size_local_rows;
+ for (i = 0; i < n; i++)
combine(one,compute_vector[i]);
- } else {
- double **compute_array = compute->array_atom;
- int jm1 = j - 1;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+ } else {
+ double **compute_array = compute->array_local;
+ int n = compute->size_local_rows;
+ int jm1 = j - 1;
+ for (i = 0; i < n; i++)
combine(one,compute_array[i][jm1]);
+ }
}
- // access fix fields, check if frequency is a match
+ // check if fix frequency is a match
} else if (which[m] == FIX) {
if (update->ntimestep % modify->fix[n]->peratom_freq)
error->all("Fix used in compute reduce not computed at compatible time");
+ Fix *fix = modify->fix[n];
+
+ if (flavor[m] == GLOBAL) {
+ if (j == 0) {
+ int n = fix->size_vector;
+ for (i = 0; i < n; i++)
+ combine(one,fix->compute_vector(i));
+ } else {
+ int n = fix->size_array_rows;
+ int jm1 = j - 1;
+ for (i = 0; i < nlocal; i++)
+ combine(one,fix->compute_array(i,jm1));
+ }
+
+ } else if (flavor[m] == PERATOM) {
+ if (j == 0) {
+ double *fix_vector = fix->vector_atom;
+ int n = nlocal;
+ for (i = 0; i < n; i++)
+ if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+ combine(one,fix_vector[i]);
+ } else {
+ double **fix_array = fix->array_atom;
+ int n = nlocal;
+ int jm1 = j - 1;
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+ combine(one,fix_array[i][jm1]);
+ }
- if (j == 0) {
- double *fix_vector = modify->fix[n]->vector_atom;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+ } else if (flavor[m] == LOCAL) {
+ if (j == 0) {
+ double *fix_vector = fix->vector_local;
+ int n = fix->size_local_rows;
+ for (i = 0; i < n; i++)
combine(one,fix_vector[i]);
- } else {
- double **fix_array = modify->fix[n]->array_atom;
- int jm1 = j - 1;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
+ } else {
+ double **fix_array = fix->array_local;
+ int n = fix->size_local_rows;
+ int jm1 = j - 1;
+ for (i = 0; i < n; i++)
combine(one,fix_array[i][jm1]);
+ }
}
// evaluate atom-style variable
@@ -136,3 +219,47 @@ double ComputeReduceRegion::compute_one(int m)
return one;
}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeReduceRegion::count(int m)
+{
+ int n = value2index[m];
+ int j = argindex[m];
+
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ if (which[m] == X || which[m] == V || which[m] == F)
+ return group->count(igroup,iregion);
+ else if (which[m] == COMPUTE) {
+ Compute *compute = modify->compute[n];
+ if (flavor[m] == GLOBAL) {
+ if (j == 0) return(1.0*compute->size_vector);
+ else return(1.0*compute->size_array_rows);
+ } else if (flavor[m] == PERATOM) {
+ return group->count(igroup,iregion);
+ } else if (flavor[m] == LOCAL) {
+ double ncount = compute->size_local_rows;
+ double ncountall;
+ MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
+ return ncountall;
+ }
+ } else if (which[m] == FIX) {
+ Fix *fix = modify->fix[n];
+ if (flavor[m] == GLOBAL) {
+ if (j == 0) return(1.0*fix->size_vector);
+ else return(1.0*fix->size_array_rows);
+ } else if (flavor[m] == PERATOM) {
+ return group->count(igroup,iregion);
+ } else if (flavor[m] == LOCAL) {
+ double ncount = fix->size_local_rows;
+ double ncountall;
+ MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
+ return ncountall;
+ }
+ } else if (which[m] == VARIABLE)
+ return group->count(igroup,iregion);
+
+ return 0.0;
+}
diff --git a/src/compute_reduce_region.h b/src/compute_reduce_region.h
index 6aaf1741e58dabba3dd2b3af446ab6e480657c8b..5ceb834d4dbb685380b5d6016d730eed28fd036f 100644
--- a/src/compute_reduce_region.h
+++ b/src/compute_reduce_region.h
@@ -25,6 +25,7 @@ class ComputeReduceRegion : public ComputeReduce {
private:
double compute_one(int);
+ double count(int);
};
}
diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp
index 43409fe79d7d899e9bb209b3fc981bf948952389..bd8aa03362b076d0d08f17b1186bfd09ee95379e 100644
--- a/src/compute_temp_sphere.cpp
+++ b/src/compute_temp_sphere.cpp
@@ -26,8 +26,6 @@
using namespace LAMMPS_NS;
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
-
#define INERTIA 0.4 // moment of inertia for sphere
/* ---------------------------------------------------------------------- */
diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp
index 1b0fb20060be8402f8882af5e208f7ea27df2c80..07b97fb0550ea7d45365f3de340fc53df86ca1d3 100644
--- a/src/dump_custom.cpp
+++ b/src/dump_custom.cpp
@@ -40,7 +40,8 @@ enum{ID,MOL,TYPE,MASS,
COMPUTE,FIX,VARIABLE};
enum{LT,LE,GT,GE,EQ,NEQ};
enum{INT,DOUBLE};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define INVOKED_PERATOM 8
/* ---------------------------------------------------------------------- */
@@ -960,7 +961,7 @@ void DumpCustom::parse_fields(int narg, char **arg)
vtype[i] = DOUBLE;
// compute value = c_ID
- // if no trailing [], then arg is set to 0, else arg is between []
+ // if no trailing [], then arg is set to 0, else arg is int between []
} else if (strncmp(arg[iarg],"c_",2) == 0) {
pack_choice[i] = &DumpCustom::pack_compute;
@@ -981,14 +982,14 @@ void DumpCustom::parse_fields(int narg, char **arg)
n = modify->find_compute(suffix);
if (n < 0) error->all("Could not find dump custom compute ID");
if (modify->compute[n]->peratom_flag == 0)
- error->all("Dump custom compute ID does not compute per-atom info");
+ error->all("Dump custom compute does not compute per-atom info");
if (argindex[i] == 0 && modify->compute[n]->size_peratom_cols > 0)
- error->all("Dump custom compute ID does not compute per-atom vector");
+ error->all("Dump custom compute does not calculate per-atom vector");
if (argindex[i] > 0 && modify->compute[n]->size_peratom_cols == 0)
- error->all("Dump custom compute ID does not compute per-atom array");
+ error->all("Dump custom compute does not calculate per-atom array");
if (argindex[i] > 0 &&
argindex[i] > modify->compute[n]->size_peratom_cols)
- error->all("Dump custom compute ID vector is not large enough");
+ error->all("Dump custom compute vector is accessed out-of-range");
field2index[i] = add_compute(suffix);
delete [] suffix;
@@ -1015,14 +1016,14 @@ void DumpCustom::parse_fields(int narg, char **arg)
n = modify->find_fix(suffix);
if (n < 0) error->all("Could not find dump custom fix ID");
if (modify->fix[n]->peratom_flag == 0)
- error->all("Dump custom fix ID does not compute per-atom info");
+ error->all("Dump custom fix does not compute per-atom info");
if (argindex[i] == 0 && modify->fix[n]->size_peratom_cols > 0)
- error->all("Dump custom fix ID does not compute per-atom vector");
+ error->all("Dump custom fix does not compute per-atom vector");
if (argindex[i] > 0 && modify->fix[n]->size_peratom_cols == 0)
- error->all("Dump custom fix ID does not compute per-atom array");
+ error->all("Dump custom fix does not compute per-atom array");
if (argindex[i] > 0 &&
argindex[i] > modify->fix[n]->size_peratom_cols)
- error->all("Dump custom fix ID vector is not large enough");
+ error->all("Dump custom fix vector is accessed out-of-range");
field2index[i] = add_fix(suffix);
delete [] suffix;
@@ -1398,7 +1399,6 @@ void DumpCustom::pack_compute(int n)
double *vector = compute[field2index[n]]->vector_atom;
double **array = compute[field2index[n]]->array_atom;
int index = argindex[n];
-
int nlocal = atom->nlocal;
if (index == 0) {
@@ -1424,7 +1424,6 @@ void DumpCustom::pack_fix(int n)
double *vector = fix[field2index[n]]->vector_atom;
double **array = fix[field2index[n]]->array_atom;
int index = argindex[n];
-
int nlocal = atom->nlocal;
if (index == 0) {
@@ -1448,7 +1447,6 @@ void DumpCustom::pack_fix(int n)
void DumpCustom::pack_variable(int n)
{
double *vector = vbuf[field2index[n]];
-
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
diff --git a/src/fix.h b/src/fix.h
index b8bd3a525a76877675c1e4da5df2f399bbd1a20e..9379d87d40fd55116fadcd070d20ecb321b770bf 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -136,7 +136,7 @@ class Fix : protected Pointers {
virtual double compute_scalar() {return 0.0;}
virtual double compute_vector(int) {return 0.0;}
- virtual double compute_array(int) {return 0.0;}
+ virtual double compute_array(int,int) {return 0.0;}
virtual int dof(int) {return 0;}
virtual void deform(int) {}
diff --git a/src/fix_ave_atom.cpp b/src/fix_ave_atom.cpp
index 3a637884e7e122747bc85761029a3c615f63f4ab..acc1ff8dc92130f806ddcbd4fd907112594349e9 100644
--- a/src/fix_ave_atom.cpp
+++ b/src/fix_ave_atom.cpp
@@ -27,7 +27,8 @@
using namespace LAMMPS_NS;
enum{X,XU,V,F,COMPUTE,FIX,VARIABLE};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define INVOKED_PERATOM 8
/* ---------------------------------------------------------------------- */
@@ -125,6 +126,7 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
}
// setup and error check
+ // for fix inputs, check that fix frequency is acceptable
if (nevery <= 0) error->all("Illegal fix ave/atom command");
if (peratom_freq < nevery || peratom_freq % nevery ||
@@ -147,19 +149,23 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
"calculate a per-atom array");
if (argindex[i] &&
argindex[i] > modify->compute[icompute]->size_peratom_cols)
- error->all("Fix ave/atom compute vector is accessed out-of-range");
+ error->all("Fix ave/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("Fix ID for fix ave/atom does not exist");
if (modify->fix[ifix]->peratom_flag == 0)
error->all("Fix ave/atom fix does not calculate per-atom values");
- if (argindex[i] && modify->fix[ifix]->size_peratom_cols != 0)
+ if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
error->all("Fix ave/atom fix does not calculate a per-atom vector");
if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
error->all("Fix ave/atom fix does not calculate a per-atom array");
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
- error->all("Fix ave/atom fix vector is accessed out-of-range");
+ error->all("Fix ave/atom fix array is accessed out-of-range");
+ if (nevery % modify->fix[ifix]->peratom_freq)
+ error->all("Fix for fix ave/atom not computed at compatible time");
+
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
@@ -240,7 +246,6 @@ int FixAveAtom::setmask()
void FixAveAtom::init()
{
// set indices and check validity of all computes,fixes,variables
- // check that fix frequency is acceptable
for (int m = 0; m < nvalues; m++) {
if (which[m] == COMPUTE) {
@@ -255,9 +260,6 @@ void FixAveAtom::init()
error->all("Fix ID for fix ave/atom does not exist");
value2index[m] = ifix;
- if (nevery % modify->fix[ifix]->peratom_freq)
- error->all("Fix for fix ave/atom not computed at compatible time");
-
} else if (which[m] == VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp
index 06ba4f8d86d740a8659e17ad57055e0f507817f6..2af3684477b0d0864b47e16e25fae5473bdd717e 100644
--- a/src/fix_ave_spatial.cpp
+++ b/src/fix_ave_spatial.cpp
@@ -36,7 +36,8 @@ enum{X,V,F,DENSITY_NUMBER,DENSITY_MASS,COMPUTE,FIX,VARIABLE};
enum{SAMPLE,ALL};
enum{BOX,LATTICE,REDUCED};
enum{ONE,RUNNING,WINDOW};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define INVOKED_PERATOM 8
#define BIG 1000000000
@@ -154,6 +155,9 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
scaleflag = LATTICE;
fp = NULL;
ave = ONE;
+ char *title1 = NULL;
+ char *title2 = NULL;
+ char *title3 = NULL;
while (iarg < narg) {
if (strcmp(arg[iarg],"norm") == 0) {
@@ -193,6 +197,27 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
}
iarg += 2;
if (ave == WINDOW) iarg++;
+ } else if (strcmp(arg[iarg],"title1") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
+ delete [] title1;
+ int n = strlen(arg[iarg+1]) + 1;
+ title1 = new char[n];
+ strcpy(title1,arg[iarg+1]);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"title2") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
+ delete [] title2;
+ int n = strlen(arg[iarg+1]) + 1;
+ title2 = new char[n];
+ strcpy(title2,arg[iarg+1]);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"title3") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
+ delete [] title3;
+ int n = strlen(arg[iarg+1]) + 1;
+ title3 = new char[n];
+ strcpy(title3,arg[iarg+1]);
+ iarg += 2;
} else error->all("Illegal fix ave/spatial command");
}
@@ -243,27 +268,38 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
}
}
- // print header into file
+ // print file comment lines
if (fp && me == 0) {
- fprintf(fp,"# Spatial-averaged data for fix %s and group %s\n",id,arg[1]);
- fprintf(fp,"# TimeStep Number-of-layers\n");
- fprintf(fp,"# Layer Coordinate Natoms");
- for (int i = 0; i < nvalues; i++)
- if (which[i] == COMPUTE) fprintf(fp," c_%s",ids[i]);
- else if (which[i] == FIX) fprintf(fp," f_%s",ids[i]);
- else if (which[i] == VARIABLE) fprintf(fp," v_%s",ids[i]);
- else fprintf(fp," %s",arg[9+i]);
- fprintf(fp,"\n");
+ if (title1) fprintf(fp,"%s\n",title1);
+ else fprintf(fp,"# Spatial-averaged data for fix %s and group %s\n",
+ id,arg[1]);
+ if (title2) fprintf(fp,"%s\n",title2);
+ else fprintf(fp,"# Timestep Number-of-layers\n");
+ if (title3) fprintf(fp,"%s\n",title3);
+ else {
+ fprintf(fp,"# Layer Coord Ncount");
+ for (int i = 0; i < nvalues; i++) {
+ if (which[i] == COMPUTE) fprintf(fp," c_%s",ids[i]);
+ else if (which[i] == FIX) fprintf(fp," f_%s",ids[i]);
+ else if (which[i] == VARIABLE) fprintf(fp," v_%s",ids[i]);
+ else fprintf(fp," %s",arg[9+i]);
+ }
+ fprintf(fp,"\n");
+ }
}
-
- // this fix produces a global vector
- // set size_vector to BIG since compute_vector() checks bounds on-the-fly
- vector_flag = 1;
- size_vector = BIG;
+ delete [] title1;
+ delete [] title2;
+ delete [] title3;
+
+ // this fix produces a global array
+
+ array_flag = 1;
+ size_local_rows = BIG;
+ size_local_cols = nvalues+2;
global_freq = nfreq;
- extvector = 0;
+ extarray = 0;
// setup scaling
@@ -803,18 +839,18 @@ void FixAveSpatial::end_of_step()
}
/* ----------------------------------------------------------------------
- return Nth vector value
- since values_sum is 2d array, map N into ilayer and ivalue
- if ilayer exceeds current layers, return 0.0 instead of generate an error
+ return I,J array value
+ if I exceeds current layers, return 0.0 instead of generating an error
+ 1st column = bin coord, 2nd column = count, remaining columns = Nvalues
------------------------------------------------------------------------- */
-double FixAveSpatial::compute_vector(int n)
+double FixAveSpatial::compute_array(int i, int j)
{
- int ivalue = n % nvalues;
- int ilayer = n / nvalues;
- if (ilayer >= nlayers) return 0.0;
if (values_total == NULL) return 0.0;
- return values_total[ilayer][ivalue]/norm;
+ if (i >= nlayers) return 0.0;
+ if (j == 0) return coord[i];
+ if (j == 1) return count_total[i]/norm;
+ return values_total[i][j]/norm;
}
/* ----------------------------------------------------------------------
diff --git a/src/fix_ave_spatial.h b/src/fix_ave_spatial.h
index 0a4da9876e414389f9a5c69a5e6e612ac457c05d..cc11d55d40261cbe3260eeeb001224005da88722 100644
--- a/src/fix_ave_spatial.h
+++ b/src/fix_ave_spatial.h
@@ -27,7 +27,7 @@ class FixAveSpatial : public Fix {
void init();
void setup(int);
void end_of_step();
- double compute_vector(int);
+ double compute_array(int,int);
double memory_usage();
private:
@@ -35,6 +35,7 @@ class FixAveSpatial : public Fix {
int nrepeat,nfreq,nvalid,irepeat;
int dim,originflag,normflag;
double origin,delta;
+ char *tstring,*sstring;
int *which,*argindex,*value2index;
char **ids;
FILE *fp;
diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp
index d042f042d18b30ad3d39249f6cfd6fc0939e4200..43824525be7cc5099ef6588c203bb07111c02f6a 100644
--- a/src/fix_ave_time.cpp
+++ b/src/fix_ave_time.cpp
@@ -31,7 +31,10 @@ using namespace LAMMPS_NS;
enum{COMPUTE,FIX,VARIABLE};
enum{ONE,RUNNING,WINDOW};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
/* ---------------------------------------------------------------------- */
@@ -129,6 +132,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
}
// setup and error check
+ // for fix inputs, check that fix frequency is acceptable
if (nevery <= 0) error->all("Illegal fix ave/time command");
if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
@@ -145,6 +149,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
error->all("Fix ave/time compute does not calculate a vector");
if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector)
error->all("Fix ave/time compute vector is accessed out-of-range");
+
} else if (which[i] == FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
@@ -155,6 +160,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
error->all("Fix ave/time fix does not calculate a vector");
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector)
error->all("Fix ave/time fix vector is accessed out-of-range");
+ if (nevery % modify->fix[ifix]->global_freq)
+ error->all("Fix for fix ave/time not computed at compatible time");
+
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
@@ -285,8 +293,7 @@ int FixAveTime::setmask()
void FixAveTime::init()
{
- // set indices and check validity of all computes,fixes,variables
- // check that fix frequency is acceptable
+ // set current indices for all computes,fixes,variables
for (int i = 0; i < nvalues; i++) {
if (which[i] == COMPUTE) {
@@ -301,9 +308,6 @@ void FixAveTime::init()
error->all("Fix ID for fix ave/time does not exist");
value2index[i] = ifix;
- if (nevery % modify->fix[ifix]->global_freq)
- error->all("Fix for fix ave/time not computed at compatible time");
-
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
@@ -440,7 +444,6 @@ void FixAveTime::end_of_step()
/* ----------------------------------------------------------------------
return scalar value
- could be ONE, RUNNING, or WINDOW value
------------------------------------------------------------------------- */
double FixAveTime::compute_scalar()
@@ -450,12 +453,21 @@ double FixAveTime::compute_scalar()
}
/* ----------------------------------------------------------------------
- return Nth vector value
- could be ONE, RUNNING, or WINDOW value
+ return Ith vector value
+------------------------------------------------------------------------- */
+
+double FixAveTime::compute_vector(int i)
+{
+ if (norm) return vector_total[i]/norm;
+ return 0.0;
+}
+
+/* ----------------------------------------------------------------------
+ return I,J array value
------------------------------------------------------------------------- */
-double FixAveTime::compute_vector(int n)
+double FixAveTime::compute_array(int i, int j)
{
- if (norm) return vector_total[n]/norm;
+ if (norm) return array_total[i][j]/norm;
return 0.0;
}
diff --git a/src/fix_ave_time.h b/src/fix_ave_time.h
index f4461583570dcf8578ccff314f6e5aac68f59a92..0d4b562cdcfcce5def5106582e278f940706c693 100644
--- a/src/fix_ave_time.h
+++ b/src/fix_ave_time.h
@@ -29,6 +29,7 @@ class FixAveTime : public Fix {
void end_of_step();
double compute_scalar();
double compute_vector(int);
+ double compute_array(int,int);
private:
int me,nvalues;
@@ -43,6 +44,7 @@ class FixAveTime : public Fix {
int norm,iwindow,window_limit;
double *vector_total;
double **vector_list;
+ double **array_total;
};
}
diff --git a/src/style.h b/src/style.h
index d6dd6003aa2cc46ae6d7da73d0c2e410f05b35f8..f124a47f49a574284d735230fd03aed7be96a3d3 100644
--- a/src/style.h
+++ b/src/style.h
@@ -90,6 +90,7 @@ CommandStyle(write_restart,WriteRestart)
#include "compute_pe.h"
#include "compute_pe_atom.h"
#include "compute_pressure.h"
+#include "compute_rdf.h"
#include "compute_reduce.h"
#include "compute_reduce_region.h"
#include "compute_erotate_sphere.h"
@@ -121,6 +122,7 @@ ComputeStyle(pe/atom,ComputePEAtom)
ComputeStyle(pressure,ComputePressure)
ComputeStyle(reduce,ComputeReduce)
ComputeStyle(reduce/region,ComputeReduceRegion)
+ComputeStyle(rdf,ComputeRDF)
ComputeStyle(erotate/sphere,ComputeERotateSphere)
ComputeStyle(stress/atom,ComputeStressAtom)
ComputeStyle(temp,ComputeTemp)
diff --git a/src/thermo.cpp b/src/thermo.cpp
index ddd26643e163864630e24c13a24c18c7d78ea959..d14220452c3c17973be0c5aabb0def725a0ce3d7 100644
--- a/src/thermo.cpp
+++ b/src/thermo.cpp
@@ -53,7 +53,11 @@ using namespace LAMMPS_NS;
enum{IGNORE,WARN,ERROR}; // same as write_restart.cpp
enum{ONELINE,MULTILINE};
enum{INT,FLOAT};
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
+enum{SCALAR,VECTOR,ARRAY};
+
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
#define MAXLINE 1024
#define DELTA 8
@@ -271,16 +275,21 @@ void Thermo::compute(int flag)
// which = 0 is global scalar, which = 1 is global vector
for (i = 0; i < ncompute; i++)
- if (compute_which[i] == 0) {
+ if (compute_which[i] == SCALAR) {
if (!(computes[i]->invoked_flag & INVOKED_SCALAR)) {
computes[i]->compute_scalar();
computes[i]->invoked_flag |= INVOKED_SCALAR;
}
- } else {
+ } else if (compute_which[i] == VECTOR) {
if (!(computes[i]->invoked_flag & INVOKED_VECTOR)) {
computes[i]->compute_vector();
computes[i]->invoked_flag |= INVOKED_VECTOR;
}
+ } else if (compute_which[i] == ARRAY) {
+ if (!(computes[i]->invoked_flag & INVOKED_ARRAY)) {
+ computes[i]->compute_array();
+ computes[i]->invoked_flag |= INVOKED_ARRAY;
+ }
}
// if lineflag = MULTILINE, prepend step/cpu header line
@@ -502,7 +511,8 @@ void Thermo::allocate()
for (int i = 0; i < n; i++) format_user[i] = NULL;
field2index = new int[n];
- argindex = new int[n];
+ argindex1 = new int[n];
+ argindex2 = new int[n];
// factor of 3 is max number of computes a single field can add
@@ -539,7 +549,8 @@ void Thermo::deallocate()
delete [] format_user;
delete [] field2index;
- delete [] argindex;
+ delete [] argindex1;
+ delete [] argindex2;
for (int i = 0; i < ncompute; i++) delete [] id_compute[i];
delete [] id_compute;
@@ -580,56 +591,56 @@ void Thermo::parse_fields(char *str)
} else if (strcmp(word,"temp") == 0) {
addfield("Temp",&Thermo::compute_temp,FLOAT);
- index_temp = add_compute(id_temp,0);
+ index_temp = add_compute(id_temp,SCALAR);
} else if (strcmp(word,"press") == 0) {
addfield("Press",&Thermo::compute_press,FLOAT);
- index_press_scalar = add_compute(id_press,0);
+ index_press_scalar = add_compute(id_press,SCALAR);
} else if (strcmp(word,"pe") == 0) {
addfield("PotEng",&Thermo::compute_pe,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"ke") == 0) {
addfield("KinEng",&Thermo::compute_ke,FLOAT);
- index_temp = add_compute(id_temp,0);
+ index_temp = add_compute(id_temp,SCALAR);
} else if (strcmp(word,"etotal") == 0) {
addfield("TotEng",&Thermo::compute_etotal,FLOAT);
- index_temp = add_compute(id_temp,0);
- index_pe = add_compute(id_pe,0);
+ index_temp = add_compute(id_temp,SCALAR);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"enthalpy") == 0) {
addfield("Enthalpy",&Thermo::compute_enthalpy,FLOAT);
- index_temp = add_compute(id_temp,0);
- index_press_scalar = add_compute(id_press,0);
- index_pe = add_compute(id_pe,0);
+ index_temp = add_compute(id_temp,SCALAR);
+ index_press_scalar = add_compute(id_press,SCALAR);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"evdwl") == 0) {
addfield("E_vdwl",&Thermo::compute_evdwl,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"ecoul") == 0) {
addfield("E_coul",&Thermo::compute_ecoul,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"epair") == 0) {
addfield("E_pair",&Thermo::compute_epair,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"ebond") == 0) {
addfield("E_bond",&Thermo::compute_ebond,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"eangle") == 0) {
addfield("E_angle",&Thermo::compute_eangle,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"edihed") == 0) {
addfield("E_dihed",&Thermo::compute_edihed,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"eimp") == 0) {
addfield("E_impro",&Thermo::compute_eimp,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"emol") == 0) {
addfield("E_mol",&Thermo::compute_emol,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"elong") == 0) {
addfield("E_long",&Thermo::compute_elong,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"etail") == 0) {
addfield("E_tail",&Thermo::compute_etail,FLOAT);
- index_pe = add_compute(id_pe,0);
+ index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"vol") == 0) {
addfield("Volume",&Thermo::compute_vol,FLOAT);
@@ -662,25 +673,25 @@ void Thermo::parse_fields(char *str)
} else if (strcmp(word,"pxx") == 0) {
addfield("Pxx",&Thermo::compute_pxx,FLOAT);
- index_press_vector = add_compute(id_press,1);
+ index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pyy") == 0) {
addfield("Pyy",&Thermo::compute_pyy,FLOAT);
- index_press_vector = add_compute(id_press,1);
+ index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pzz") == 0) {
addfield("Pzz",&Thermo::compute_pzz,FLOAT);
- index_press_vector = add_compute(id_press,1);
+ index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pxy") == 0) {
addfield("Pxy",&Thermo::compute_pxy,FLOAT);
- index_press_vector = add_compute(id_press,1);
+ index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pxz") == 0) {
addfield("Pxz",&Thermo::compute_pxz,FLOAT);
- index_press_vector = add_compute(id_press,1);
+ index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pyz") == 0) {
addfield("Pyz",&Thermo::compute_pyz,FLOAT);
- index_press_vector = add_compute(id_press,1);
+ index_press_vector = add_compute(id_press,VECTOR);
// compute value = c_ID, fix value = f_ID, variable value = v_ID
- // if no trailing [], then arg is set to 0, else arg is between []
+ // count trailing [] and store int arguments
// copy = at most 8 chars of ID to pass to addfield
} else if ((strncmp(word,"c_",2) == 0) || (strncmp(word,"f_",2) == 0) ||
@@ -693,38 +704,66 @@ void Thermo::parse_fields(char *str)
strncpy(copy,id,8);
copy[8] = '\0';
+ // parse zero or one or two trailing brackets from ID
+ // argindex1,argindex2 = int inside each bracket pair, 0 if no bracket
+
char *ptr = strchr(id,'[');
- if (ptr) {
- if (id[strlen(id)-1] != ']')
- error->all("Invalid keyword in thermo_style custom command");
- argindex[nfield] = atoi(ptr+1);
+ if (ptr == NULL) argindex1[nfield] = 0;
+ else {
*ptr = '\0';
- } else argindex[nfield] = 0;
+ argindex1[nfield] = input->variable->int_between_brackets(ptr);
+ ptr++;
+ if (*ptr == '[') {
+ argindex2[nfield] = input->variable->int_between_brackets(ptr);
+ ptr++;
+ } else argindex2[nfield] = 0;
+ }
if (word[0] == 'c') {
n = modify->find_compute(id);
if (n < 0) error->all("Could not find thermo custom compute ID");
- if (argindex[nfield] == 0 && modify->compute[n]->scalar_flag == 0)
- error->all("Thermo compute ID does not compute scalar info");
- if (argindex[nfield] > 0 && modify->compute[n]->vector_flag == 0)
- error->all("Thermo compute ID does not compute vector info");
- if (argindex[nfield] > 0 &&
- argindex[nfield] > modify->compute[n]->size_vector)
- error->all("Thermo compute ID vector is not large enough");
-
- field2index[nfield] = add_compute(id,MIN(argindex[nfield],1));
+ if (argindex1[nfield] == 0 && modify->compute[n]->scalar_flag == 0)
+ error->all("Thermo compute does not compute scalar");
+ if (argindex1[nfield] > 0 && argindex2[nfield] == 0) {
+ if (modify->compute[n]->vector_flag == 0)
+ error->all("Thermo compute does not compute vector");
+ if (argindex1[nfield] > modify->compute[n]->size_vector)
+ error->all("Thermo compute vector is accessed out-of-range");
+ }
+ if (argindex1[nfield] > 0 && argindex2[nfield] > 0) {
+ if (modify->compute[n]->array_flag == 0)
+ error->all("Thermo compute does not compute array");
+ if (argindex1[nfield] > modify->compute[n]->size_array_rows ||
+ argindex2[nfield] > modify->compute[n]->size_array_cols)
+ error->all("Thermo compute array is accessed out-of-range");
+ }
+
+ if (argindex1[nfield] == 0)
+ field2index[nfield] = add_compute(id,SCALAR);
+ else if (argindex2[nfield] == 0)
+ field2index[nfield] = add_compute(id,VECTOR);
+ else
+ field2index[nfield] = add_compute(id,ARRAY);
addfield(copy,&Thermo::compute_compute,FLOAT);
} else if (word[0] == 'f') {
n = modify->find_fix(id);
if (n < 0) error->all("Could not find thermo custom fix ID");
- if (argindex[nfield] == 0 && modify->fix[n]->scalar_flag == 0)
- error->all("Thermo fix ID does not compute scalar info");
- if (argindex[nfield] > 0 && modify->fix[n]->vector_flag == 0)
- error->all("Thermo fix ID does not compute vector info");
- if (argindex[nfield] > 0 &&
- argindex[nfield] > modify->fix[n]->size_vector)
- error->all("Thermo fix ID vector is not large enough");
+ if (argindex1[nfield] == 0 && modify->fix[n]->scalar_flag == 0)
+ error->all("Thermo fix does not compute scalar");
+ if (argindex1[nfield] > 0 && argindex2[nfield] == 0) {
+ if (modify->fix[n]->vector_flag == 0)
+ error->all("Thermo fix does not compute vector");
+ if (argindex1[nfield] > modify->fix[n]->size_vector)
+ error->all("Thermo fix vector is accessed out-of-range");
+ }
+ if (argindex1[nfield] > 0 && argindex2[nfield] > 0) {
+ if (modify->fix[n]->array_flag == 0)
+ error->all("Thermo fix does not compute array");
+ if (argindex1[nfield] > modify->fix[n]->size_array_rows ||
+ argindex2[nfield] > modify->fix[n]->size_array_cols)
+ error->all("Thermo fix array is accessed out-of-range");
+ }
field2index[nfield] = add_fix(id);
addfield(copy,&Thermo::compute_fix,FLOAT);
@@ -734,6 +773,8 @@ void Thermo::parse_fields(char *str)
if (n < 0) error->all("Could not find thermo custom variable name");
if (input->variable->equalstyle(n) == 0)
error->all("Thermo custom variable is not equal-style variable");
+ if (argindex1[nfield])
+ error->all("Thermo custom variable cannot be indexed");
field2index[nfield] = add_variable(id);
addfield(copy,&Thermo::compute_variable,FLOAT);
@@ -1150,17 +1191,22 @@ int Thermo::evaluate_keyword(char *word, double *answer)
void Thermo::compute_compute()
{
- Compute *compute = computes[field2index[ifield]];
- if (argindex[ifield] == 0) {
+ int m = field2index[ifield];
+ Compute *compute = computes[m];
+
+ if (compute_which[m] == SCALAR) {
dvalue = compute->scalar;
if (normflag && compute->extscalar) dvalue /= natoms;
- } else {
- dvalue = compute->vector[argindex[ifield]-1];
+ } else if (compute_which[m] == VECTOR) {
+ dvalue = compute->vector[argindex1[ifield]-1];
if (normflag) {
if (compute->extvector == 0) return;
else if (compute->extvector == 1) dvalue /= natoms;
- else if (compute->extlist[argindex[ifield]-1]) dvalue /= natoms;
+ else if (compute->extlist[argindex1[ifield]-1]) dvalue /= natoms;
}
+ } else {
+ dvalue = compute->array[argindex1[ifield]-1][argindex2[ifield]-1];
+ if (normflag && compute->extarray) dvalue /= natoms;
}
}
@@ -1168,17 +1214,22 @@ void Thermo::compute_compute()
void Thermo::compute_fix()
{
- Fix *fix = fixes[field2index[ifield]];
- if (argindex[ifield] == 0) {
+ int m = field2index[ifield];
+ Fix *fix = fixes[m];
+
+ if (argindex1[ifield] == 0) {
dvalue = fix->compute_scalar();
if (normflag && fix->extscalar) dvalue /= natoms;
- } else {
- dvalue = fix->compute_vector(argindex[ifield]-1);
+ } else if (argindex2[ifield] == 0) {
+ dvalue = fix->compute_vector(argindex1[ifield]-1);
if (normflag) {
if (fix->extvector == 0) return;
else if (fix->extvector == 1) dvalue /= natoms;
- else if (fix->extlist[argindex[ifield]-1]) dvalue /= natoms;
+ else if (fix->extlist[argindex1[ifield]-1]) dvalue /= natoms;
}
+ } else {
+ dvalue = fix->compute_array(argindex1[ifield]-1,argindex2[ifield]-1);
+ if (normflag && fix->extarray) dvalue /= natoms;
}
}
diff --git a/src/thermo.h b/src/thermo.h
index 4c47a978ecdac80070624fa3c798920df5e5ad4e..06bb96416ed567eca575df7e8acb3def1ec583e0 100644
--- a/src/thermo.h
+++ b/src/thermo.h
@@ -62,8 +62,8 @@ class Thermo : protected Pointers {
double dvalue,natoms; // dvalue = double value to print
int ifield; // which field in thermo output is being computed
int *field2index; // which compute,fix,variable calcs this field
- int *argindex; // index into compute,fix scalar,vector
-
+ int *argindex1; // indices into compute,fix scalar,vector
+ int *argindex2;
// data for keyword-specific Compute objects
// index = where they are in computes list
// id = ID of Compute objects
diff --git a/src/variable.cpp b/src/variable.cpp
index 3c6edb55a2d3152123267944af08b4b33c571a17..28b1c63b351d8a3d57975a36dd828e79e1a8e727 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -37,13 +37,17 @@ using namespace LAMMPS_NS;
#define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
-enum{DUMMY0,INVOKED_SCALAR,INVOKED_VECTOR,DUMMMY3,INVOKED_PERATOM};
enum{INDEX,LOOP,EQUAL,WORLD,UNIVERSE,ULOOP,ATOM};
enum{ARG,OP};
enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN,
CEIL,FLOOR,ROUND,VALUE,ATOMARRAY,TYPEARRAY,INTARRAY};
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+#define INVOKED_ARRAY 4
+#define INVOKED_PERATOM 8
+
/* ---------------------------------------------------------------------- */
Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
@@ -524,11 +528,11 @@ void Variable::copy(int narg, char **from, char **to)
sqrt(x),exp(x),ln(x),log(x),
sin(x),cos(x),tan(x),asin(x),acos(x),atan(x)
group function = count(group), mass(group), xcm(group,x), ...
- atom value = x[123], y[3], vx[34], ...
- atom vector = x[], y[], vx[], ...
- compute = c_ID, c_ID[2], c_ID[], c_ID[][2], c_ID[N], c_ID[N][2]
- fix = f_ID, f_ID[2], f_ID[], f_ID[][2], f_ID[N], f_ID[N][2]
- variable = v_name, v_name[], v_name[N]
+ atom value = x[i], y[i], vx[i], ...
+ atom vector = x, y, vx, ...
+ compute = c_ID, c_ID[i], c_ID[i][j]
+ fix = f_ID, f_ID[i], f_ID[i][j]
+ variable = v_name, v_name[i]
if tree ptr passed in, return ptr to parse tree created for formula
if no tree ptr passed in, return value of string as double
------------------------------------------------------------------------- */
@@ -538,6 +542,7 @@ double Variable::evaluate(char *str, Tree **tree)
int op,opprevious;
double value1,value2;
char onechar;
+ char *ptr;
double argstack[MAXLEVEL];
Tree *treestack[MAXLEVEL];
@@ -613,7 +618,8 @@ double Variable::evaluate(char *str, Tree **tree)
delete [] number;
// ----------------
- // letter: c_ID, f_ID, v_name, exp(), xcm(,), x[234], x[], vol
+ // letter: c_ID, c_ID[], c_ID[][], f_ID, f_ID[], f_ID[][],
+ // v_name, v_name[], exp(), xcm(,), x, x[], vol
// ----------------
} else if (islower(onechar)) {
@@ -637,6 +643,9 @@ double Variable::evaluate(char *str, Tree **tree)
// ----------------
if (strncmp(word,"c_",2) == 0) {
+ if (domain->box_exist == 0)
+ error->all("Variable evaluation before simulation box is defined");
+
n = strlen(word) - 2 + 1;
char *id = new char[n];
strcpy(id,&word[2]);
@@ -646,28 +655,27 @@ double Variable::evaluate(char *str, Tree **tree)
Compute *compute = modify->compute[icompute];
delete [] id;
- if (domain->box_exist == 0)
- error->all("Variable evaluation before simulation box is defined");
-
// parse zero or one or two trailing brackets
// point i beyond last bracket
// nbracket = # of bracket pairs
- // index1,index2 = int inside each bracket pair, 0 if empty
+ // index1,index2 = int inside each bracket pair
int nbracket,index1,index2;
if (str[i] != '[') nbracket = 0;
else {
nbracket = 1;
- i = int_between_brackets(str,i,index1,1);
- i++;
+ ptr = &str[i];
+ index1 = int_between_brackets(ptr);
+ i = ptr-str+1;
if (str[i] == '[') {
nbracket = 2;
- i = int_between_brackets(str,i,index2,0);
- i++;
+ ptr = &str[i];
+ index2 = int_between_brackets(ptr);
+ i = ptr-str+1;
}
}
- // c_ID = global scalar
+ // c_ID = scalar from global scalar
if (nbracket == 0 && compute->scalar_flag) {
@@ -689,12 +697,13 @@ double Variable::evaluate(char *str, Tree **tree)
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
- // c_ID[2] = global vector
+ // c_ID[i] = scalar from global vector
- } else if (nbracket == 1 && index1 > 0 && compute->vector_flag) {
+ } else if (nbracket == 1 && compute->vector_flag) {
if (index1 > compute->size_vector)
- error->all("Compute vector in variable formula is too small");
+ error->all("Variable formula compute vector "
+ "is accessed out-of-range");
if (update->whichflag == 0) {
if (compute->invoked_vector != update->ntimestep)
error->all("Compute used in variable between runs "
@@ -713,13 +722,39 @@ double Variable::evaluate(char *str, Tree **tree)
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
- // c_ID[] = per-atom vector
+ // c_ID[i][j] = scalar from global array
- } else if (nbracket == 1 && index1 == 0 &&
- compute->peratom_flag && compute->size_peratom_cols == 0) {
+ } else if (nbracket == 2 && compute->array_flag) {
+
+ if (index1 > compute->size_array_rows)
+ error->all("Variable formula compute vector "
+ "is accessed out-of-range");
+ if (index2 > compute->size_array_cols)
+ error->all("Variable formula compute vector "
+ "is accessed out-of-range");
+ if (update->whichflag == 0) {
+ if (compute->invoked_array != update->ntimestep)
+ error->all("Compute used in variable between runs "
+ "is not current");
+ } else if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+ compute->compute_array();
+ compute->invoked_flag |= INVOKED_ARRAY;
+ }
+
+ value1 = compute->array[index1-1][index2-1];
+ if (tree) {
+ Tree *newtree = new Tree();
+ newtree->type = VALUE;
+ newtree->value = value1;
+ newtree->left = newtree->right = NULL;
+ treestack[ntreestack++] = newtree;
+ } else argstack[nargstack++] = value1;
+
+ // c_ID[i] = scalar from per-atom vector
+
+ } else if (nbracket == 1 && compute->peratom_flag &&
+ compute->size_peratom_cols == 0) {
- if (tree == NULL)
- error->all("Per-atom compute in equal-style variable formula");
if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep)
error->all("Compute used in variable between runs "
@@ -729,18 +764,17 @@ double Variable::evaluate(char *str, Tree **tree)
compute->invoked_flag |= INVOKED_PERATOM;
}
- Tree *newtree = new Tree();
- newtree->type = ATOMARRAY;
- newtree->array = compute->vector_atom;
- newtree->nstride = 1;
- newtree->left = newtree->right = NULL;
- treestack[ntreestack++] = newtree;
+ peratom2global(1,NULL,compute->vector_atom,1,index1,
+ tree,treestack,ntreestack,argstack,nargstack);
- // c_ID[N] = global value from per-atom vector
+ // c_ID[i][j] = scalar from per-atom array
- } else if (nbracket == 1 && index1 > 0 &&
- compute->peratom_flag && compute->size_peratom_cols == 0) {
+ } else if (nbracket == 2 && compute->peratom_flag &&
+ compute->size_peratom_cols > 0) {
+ if (index2 > compute->size_peratom_cols)
+ error->all("Variable formula compute vector "
+ "is accessed out-of-range");
if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep)
error->all("Compute used in variable between runs "
@@ -750,18 +784,17 @@ double Variable::evaluate(char *str, Tree **tree)
compute->invoked_flag |= INVOKED_PERATOM;
}
- peratom2global(1,NULL,compute->vector_atom,1,index1,
+ peratom2global(1,NULL,&compute->array_atom[0][index2-1],
+ compute->size_peratom_cols,index1,
tree,treestack,ntreestack,argstack,nargstack);
- // c_ID[][2] = per-atom array
+ // c_ID = vector from per-atom vector
- } else if (nbracket == 2 && index1 == 0 && index2 > 0 &&
- compute->peratom_flag) {
+ } else if (nbracket == 0 && compute->peratom_flag &&
+ compute->size_peratom_cols == 0) {
if (tree == NULL)
error->all("Per-atom compute in equal-style variable formula");
- if (index2 > compute->size_peratom_cols)
- error->all("Compute vector in variable formula is too small");
if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep)
error->all("Compute used in variable between runs "
@@ -773,18 +806,21 @@ double Variable::evaluate(char *str, Tree **tree)
Tree *newtree = new Tree();
newtree->type = ATOMARRAY;
- newtree->array = &compute->array_atom[0][index2-1];
- newtree->nstride = compute->size_peratom_cols;
+ newtree->array = compute->vector_atom;
+ newtree->nstride = 1;
newtree->left = newtree->right = NULL;
treestack[ntreestack++] = newtree;
- // c_ID[N][2] = global value from per-atom array
+ // c_ID[i] = vector from per-atom array
- } else if (nbracket == 2 && index1 > 0 && index2 > 0 &&
- compute->peratom_flag) {
+ } else if (nbracket == 1 && compute->peratom_flag &&
+ compute->size_peratom_cols > 0) {
- if (index2 > compute->size_peratom_cols)
- error->all("Compute vector in variable formula is too small");
+ if (tree == NULL)
+ error->all("Per-atom compute in equal-style variable formula");
+ if (index1 > compute->size_peratom_cols)
+ error->all("Variable formula compute vector "
+ "is accessed out-of-range");
if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep)
error->all("Compute used in variable between runs "
@@ -794,9 +830,12 @@ double Variable::evaluate(char *str, Tree **tree)
compute->invoked_flag |= INVOKED_PERATOM;
}
- peratom2global(1,NULL,&compute->array_atom[0][index2-1],
- compute->size_peratom_cols,index1,
- tree,treestack,ntreestack,argstack,nargstack);
+ Tree *newtree = new Tree();
+ newtree->type = ATOMARRAY;
+ newtree->array = &compute->array_atom[0][index1-1];
+ newtree->nstride = compute->size_peratom_cols;
+ newtree->left = newtree->right = NULL;
+ treestack[ntreestack++] = newtree;
} else error->all("Mismatched compute in variable formula");
@@ -805,6 +844,9 @@ double Variable::evaluate(char *str, Tree **tree)
// ----------------
} else if (strncmp(word,"f_",2) == 0) {
+ if (domain->box_exist == 0)
+ error->all("Variable evaluation before simulation box is defined");
+
n = strlen(word) - 2 + 1;
char *id = new char[n];
strcpy(id,&word[2]);
@@ -814,33 +856,33 @@ double Variable::evaluate(char *str, Tree **tree)
Fix *fix = modify->fix[ifix];
delete [] id;
- if (domain->box_exist == 0)
- error->all("Variable evaluation before simulation box is defined");
-
// parse zero or one or two trailing brackets
// point i beyond last bracket
// nbracket = # of bracket pairs
- // index1,index2 = int inside each bracket pair, 0 if empty
+ // index1,index2 = int inside each bracket pair
int nbracket,index1,index2;
if (str[i] != '[') nbracket = 0;
else {
nbracket = 1;
- i = int_between_brackets(str,i,index1,1);
- i++;
+ ptr = &str[i];
+ index1 = int_between_brackets(ptr);
+ i = ptr-str+1;
if (str[i] == '[') {
nbracket = 2;
- i = int_between_brackets(str,i,index2,0);
- i++;
+ ptr = &str[i];
+ index2 = int_between_brackets(ptr);
+ i = ptr-str+1;
}
}
- // f_ID = global scalar
+ // f_ID = scalar from global scalar
if (nbracket == 0 && fix->scalar_flag) {
if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
error->all("Fix in variable not computed at compatible time");
+
value1 = fix->compute_scalar();
if (tree) {
Tree *newtree = new Tree();
@@ -850,14 +892,15 @@ double Variable::evaluate(char *str, Tree **tree)
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
- // f_ID[2] = global vector
+ // f_ID[i] = scalar from global vector
- } else if (nbracket == 1 && index1 > 0 && fix->vector_flag) {
+ } else if (nbracket == 1 && fix->vector_flag) {
if (index1 > fix->size_vector)
- error->all("Fix vector in variable formula is too small");
+ error->all("Variable formula fix vector is accessed out-of-range");
if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
error->all("Fix in variable not computed at compatible time");
+
value1 = fix->compute_vector(index1-1);
if (tree) {
Tree *newtree = new Tree();
@@ -867,66 +910,91 @@ double Variable::evaluate(char *str, Tree **tree)
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = value1;
- // f_ID[] = per-atom vector
+ // f_ID[i][j] = scalar from global array
- } else if (nbracket == 1 && index1 == 0 &&
- fix->peratom_flag && fix->size_peratom_cols == 0) {
+ } else if (nbracket == 2 && fix->array_flag) {
+
+ if (index1 > fix->size_array_rows)
+ error->all("Variable formula fix vector is accessed out-of-range");
+ if (index2 > fix->size_array_cols)
+ error->all("Variable formula fix vector is accessed out-of-range");
+ if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+ error->all("Fix in variable not computed at compatible time");
+
+ value1 = fix->compute_array(index1-1,index2-1);
+ if (tree) {
+ Tree *newtree = new Tree();
+ newtree->type = VALUE;
+ newtree->value = value1;
+ newtree->left = newtree->right = NULL;
+ treestack[ntreestack++] = newtree;
+ } else argstack[nargstack++] = value1;
+
+ // f_ID[i] = scalar from per-atom vector
+
+ } else if (nbracket == 1 && fix->peratom_flag &&
+ fix->size_peratom_cols == 0) {
- if (tree == NULL)
- error->all("Per-atom fix in equal-style variable formula");
if (update->whichflag > 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
- Tree *newtree = new Tree();
- newtree->type = ATOMARRAY;
- newtree->array = fix->vector_atom;
- newtree->nstride = 1;
- newtree->left = newtree->right = NULL;
- treestack[ntreestack++] = newtree;
- // f_ID[N] = global value from per-atom vector
+ peratom2global(1,NULL,fix->vector_atom,1,index1,
+ tree,treestack,ntreestack,argstack,nargstack);
- } else if (nbracket == 1 && index1 > 0 &&
- fix->peratom_flag && fix->size_peratom_cols == 0) {
+ // f_ID[i][j] = scalar from per-atom array
+ } else if (nbracket == 2 && fix->peratom_flag &&
+ fix->size_peratom_cols > 0) {
+
+ if (index2 > fix->size_peratom_cols)
+ error->all("Variable formula fix vector is accessed out-of-range");
if (update->whichflag > 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
- peratom2global(1,NULL,fix->vector_atom,1,index1,
+
+ peratom2global(1,NULL,&fix->array_atom[0][index2-1],
+ fix->size_peratom_cols,index1,
tree,treestack,ntreestack,argstack,nargstack);
- // f_ID[][2] = per-atom array
+ // f_ID = vector from per-atom vector
- } else if (nbracket == 2 && index1 == 0 && index2 > 0 &&
- fix->peratom_flag) {
+ } else if (nbracket == 0 && fix->peratom_flag &&
+ fix->size_peratom_cols == 0) {
if (tree == NULL)
error->all("Per-atom fix in equal-style variable formula");
- if (index2 > fix->size_peratom_cols)
- error->all("Fix vector in variable formula is too small");
if (update->whichflag > 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
+
Tree *newtree = new Tree();
newtree->type = ATOMARRAY;
- newtree->array = &fix->array_atom[0][index2-1];
- newtree->nstride = fix->size_peratom_cols;
+ newtree->array = fix->vector_atom;
+ newtree->nstride = 1;
newtree->left = newtree->right = NULL;
treestack[ntreestack++] = newtree;
- // f_ID[N][2] = global value from per-atom array
+ // f_ID[i] = vector from per-atom array
- } else if (nbracket == 2 && index1 > 0 && index2 > 0 &&
- fix->peratom_flag) {
+ } else if (nbracket == 1 && fix->peratom_flag &&
+ fix->size_peratom_cols > 0) {
- if (index2 > fix->size_peratom_cols)
- error->all("Fix vector in variable formula is too small");
+ if (tree == NULL)
+ error->all("Per-atom fix in equal-style variable formula");
+ if (index1 > fix->size_peratom_cols)
+ error->all("Variable formula fix vector is accessed out-of-range");
if (update->whichflag > 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
- peratom2global(1,NULL,&fix->array_atom[0][index2-1],
- fix->size_peratom_cols,index1,
- tree,treestack,ntreestack,argstack,nargstack);
+
+ Tree *newtree = new Tree();
+ newtree->type = ATOMARRAY;
+ newtree->array = &fix->array_atom[0][index1-1];
+ newtree->nstride = fix->size_peratom_cols;
+ newtree->left = newtree->right = NULL;
+ treestack[ntreestack++] = newtree;
+
} else error->all("Mismatched fix in variable formula");
@@ -945,14 +1013,15 @@ double Variable::evaluate(char *str, Tree **tree)
// parse zero or one trailing brackets
// point i beyond last bracket
// nbracket = # of bracket pairs
- // index = int inside bracket, 0 if empty
+ // index = int inside bracket
int nbracket,index;
if (str[i] != '[') nbracket = 0;
else {
nbracket = 1;
- i = int_between_brackets(str,i,index,1);
- i++;
+ ptr = &str[i];
+ index = int_between_brackets(ptr);
+ i = ptr-str+1;
}
// v_name = non atom-style variable = global value
@@ -970,9 +1039,9 @@ double Variable::evaluate(char *str, Tree **tree)
treestack[ntreestack++] = newtree;
} else argstack[nargstack++] = atof(var);
- // v_name[] = atom-style variable
+ // v_name = atom-style variable
- } else if (nbracket && index == 0 && style[ivar] == ATOM) {
+ } else if (nbracket == 0 && style[ivar] == ATOM) {
if (tree == NULL)
error->all("Atom-style variable in equal-style variable formula");
@@ -984,7 +1053,7 @@ double Variable::evaluate(char *str, Tree **tree)
// compute the per-atom variable in result
// use peratom2global to extract single value from result
- } else if (nbracket && index > 0 && style[ivar] == ATOM) {
+ } else if (nbracket && style[ivar] == ATOM) {
double *result = (double *)
memory->smalloc(atom->nlocal*sizeof(double),"variable:result");
@@ -1003,7 +1072,9 @@ double Variable::evaluate(char *str, Tree **tree)
} else {
+ // ----------------
// math or group function
+ // ----------------
if (str[i] == '(') {
char *contents;
@@ -1019,25 +1090,29 @@ double Variable::evaluate(char *str, Tree **tree)
delete [] contents;
// ----------------
- // atom value or vector
+ // atom value
// ----------------
} else if (str[i] == '[') {
- int id;
- i = int_between_brackets(str,i,id,1);
- i++;
-
if (domain->box_exist == 0)
error->all("Variable evaluation before simulation box is defined");
+
+ ptr = &str[i];
+ int id = int_between_brackets(ptr);
+ i = ptr-str+1;
- // ID between brackets exists: atom value
- // empty brackets: atom vector
+ peratom2global(0,word,NULL,0,id,
+ tree,treestack,ntreestack,argstack,nargstack);
- if (id > 0)
- peratom2global(0,word,NULL,0,id,
- tree,treestack,ntreestack,argstack,nargstack);
- else
- atom_vector(word,tree,treestack,ntreestack);
+ // ----------------
+ // atom vector
+ // ----------------
+
+ } else if (is_atom_vector(word)) {
+ if (domain->box_exist == 0)
+ error->all("Variable evaluation before simulation box is defined");
+
+ atom_vector(word,tree,treestack,ntreestack);
// ----------------
// thermo keyword
@@ -1270,37 +1345,33 @@ int Variable::find_matching_paren(char *str, int i,char *&contents)
}
/* ----------------------------------------------------------------------
- find int between brackets and set index to its value
- if emptyflag, then brackets can be empty (index = 0)
- else if not emptyflag and brackets empty, is an error
- else contents of brackets must be positive int
- i = left bracket
- return loc of right bracket
+ find int between brackets and return it
+ ptr initially points to left bracket
+ return it pointing to right bracket
+ error if no right bracket or brackets are empty
+ error if any between-bracket chars are non-digits or value == 0
------------------------------------------------------------------------- */
-int Variable::int_between_brackets(char *str, int i, int &index, int emptyflag)
+int Variable::int_between_brackets(char *&ptr)
{
- // istop = matching ']'
+ char *start = ++ptr;
- int istart = i;
- while (!str[i] || str[i] != ']') i++;
- if (!str[i]) error->all("Invalid syntax in variable formula");
- int istop = i;
-
- int n = istop - istart - 1;
- char *arg = new char[n+1];
- strncpy(arg,&str[istart+1],n);
- arg[n] = '\0';
-
- if (n == 0 && emptyflag) index = 0;
- else if (n == 0) error->all("Empty brackets in variable formula");
- else {
- index = atoi(arg);
- if (index <= 0) error->all("Invalid index in variable formula");
+ while (*ptr && *ptr != ']') {
+ if (!isdigit(*ptr))
+ error->all("Non digit character between brackets in input command");
+ ptr++;
}
- delete [] arg;
-
- return istop;
+
+ if (*ptr != ']') error->all("Mismatched brackets in input command");
+ if (ptr == start) error->all("Empty brackets in input command");
+
+ *ptr = '\0';
+ int index = atoi(start);
+ *ptr = ']';
+
+ if (index == 0)
+ error->all("Index between input command brackets must be positive");
+ return index;
}
/* ----------------------------------------------------------------------
@@ -1645,6 +1716,28 @@ void Variable::peratom2global(int flag, char *word,
} else argstack[nargstack++] = value;
}
+/* ----------------------------------------------------------------------
+ check if word matches an atom vector
+ return 1 if yes, else 0
+ customize by adding an atom vector: mass,type,x,y,z,vx,vy,vz,fx,fy,fz
+------------------------------------------------------------------------- */
+
+int Variable::is_atom_vector(char *word)
+{
+ if (strcmp(word,"mass") == 0) return 1;
+ if (strcmp(word,"type") == 0) return 1;
+ if (strcmp(word,"x") == 0) return 1;
+ if (strcmp(word,"y") == 0) return 1;
+ if (strcmp(word,"z") == 0) return 1;
+ if (strcmp(word,"vx") == 0) return 1;
+ if (strcmp(word,"vy") == 0) return 1;
+ if (strcmp(word,"vz") == 0) return 1;
+ if (strcmp(word,"fx") == 0) return 1;
+ if (strcmp(word,"fy") == 0) return 1;
+ if (strcmp(word,"fz") == 0) return 1;
+ return 0;
+}
+
/* ----------------------------------------------------------------------
process an atom vector in formula
push result onto tree
@@ -1686,6 +1779,4 @@ void Variable::atom_vector(char *word, Tree **tree,
else if (strcmp(word,"fx") == 0) newtree->array = &atom->f[0][0];
else if (strcmp(word,"fy") == 0) newtree->array = &atom->f[0][1];
else if (strcmp(word,"fz") == 0) newtree->array = &atom->f[0][2];
-
- else error->all("Invalid atom vector in variable formula");
}
diff --git a/src/variable.h b/src/variable.h
index 457b759ac0b4ea31283ab30d121e1f0a7ba292fc..6afe3bf8a32a4e16ce62acbe13c0b7443c7be217 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -31,6 +31,7 @@ class Variable : protected Pointers {
char *retrieve(char *);
double compute_equal(int);
void compute_atom(int, int, double *, int, int);
+ int int_between_brackets(char *&);
private:
int me;
@@ -59,12 +60,12 @@ class Variable : protected Pointers {
double eval_tree(Tree *, int);
void free_tree(Tree *);
int find_matching_paren(char *, int, char *&);
- int int_between_brackets(char *, int, int &, int);
int math_function(char *, char *, Tree **, Tree **, int &, double *, int &);
int group_function(char *, char *, Tree **, Tree **, int &, double *, int &);
int region_function(char *);
void peratom2global(int, char *, double *, int, int,
Tree **, Tree **, int &, double *, int &);
+ int is_atom_vector(char *);
void atom_vector(char *, Tree **, Tree **, int &);
};