From 1022815b1b063e387d61fcf6219c6e5146a95f4d Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 8 Jan 2008 23:13:53 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1351 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- doc/variable.html | 124 +++++++------ doc/variable.txt | 120 ++++++------ src/variable.cpp | 452 ++++++++++++++++++++++++++++++---------------- src/variable.h | 3 +- 4 files changed, 435 insertions(+), 264 deletions(-) diff --git a/doc/variable.html b/doc/variable.html index d31356eb1d..c59bc68e1e 100644 --- a/doc/variable.html +++ b/doc/variable.html @@ -25,7 +25,7 @@ universe args = one or more strings uloop args = N = integer size of loop equal or atom args = one formula containing numbers, thermo keywords, math operations, group functions, atom values and vectors, compute and fix references, other variables - numbers = 0.0, -5.4, 2.8e-4, etc + numbers = 0.0, 100, -5.4, 2.8e-4, etc thermo keywords = vol, ke, press, etc from thermo_style math operations = (), -x, x+y, x-y, x*y, x/y, x^y, sqrt(x), exp(x), ln(x) group functions = count(group), mass(group), charge(group), @@ -35,9 +35,9 @@ vx[N], vy[N], vz[N], fx[N], fy[N], fz[N] atom vector = mass[], x[], y[], z[], vx[], vy[], vz[], fx[], fy[], fz[] - compute references = c_ID, c_ID[N] - fix references = f_ID, f_ID[N] - other variables = v_abc, v_x, etc + compute references = c_ID, c_ID[2], c_ID[N], c_ID[N][2], c_ID[], c_ID[][2] + fix references = f_ID, f_ID[2], f_ID[N], f_ID[N][2], f_ID[], f_ID[][2] + other variables = v_abc, v_abc[N], v_abc[] @@ -210,24 +210,25 @@ operations, group functions, atom values, atom vectors, compute references, fix references, and references to other variables.

- + - - - + +
Number 0.2, 1.0e20, -15.4, etc
Number 0.2, 100, 1.0e20, -15.4, etc
Thermo keywords vol, pe, ebond, etc
Math operations (), -x, x+y, x-y, x*y, x/y, x^y, sqrt(x), exp(x), ln(x)
Group functions count(ID), mass(ID), charge(ID), xcm(ID,dim), vcm(ID,dim), fcm(ID,dim), bound(ID,dir), gyration(ID)
Atom values mass[N], x[N], y[N], z[N], vx[N], vy[N], vz[N], fx[N], fy[N], fz[N]
Atom vectors mass[], x[], y[], z[], vx[], vy[], vz[], fx[], fy[], fz[]
Compute references c_ID, c_ID[N]
Fix references f_ID, f_ID[N]
Other variables v_abc, v_x, etc +
Compute references c_ID, c_ID[2], c_ID[N], c_ID[N][2], c_ID[], c_ID[][2]
Fix references f_ID, f_ID[2], f_ID[N], f_ID[N][2], f_ID[], f_ID[][2]
Other variables v_abc, v_abc[N], v_abc[]
-

Note that a formula for equal-style variables cannot use any formula -element that produces per-atom values. This includes atom vectors, a -compute that produces per-atom values, a fix that produces per-atoms -values, or an atom-style variable. A formula for an atom-style -variables can use these elements, as well as elements that produce a -global scalar or vector, e.g. a compute or fix that produces global -quantities. +

Note that formula elements that contain empty brackets, such as an +atom vector, produce per-atom values. All other formula elements +produce a global value. +

+

A formula for equal-style variables cannot use any formula element +that produces per-atom values. A formula for an atom-style variable +can use formula elements that produce either global values or per-atom +values.

The thermo keywords allowed in a formula are those defined by the "thermo_style custom" command. Since many thermodyanmic quantities @@ -262,50 +263,65 @@ gyration command for a definition of the formula. desired atom-ID, e.g. x[243]., which means use the x coordinate of the atom with ID=243.

-

Atom vectors take no argument. They generate one value per atom, so -that a reference like x[] means the x-coord of each atom will be -used when evaluating the variable. +

Atom vectors use empty brackets, i.e. they take no argument. They +generate one value per atom, so that a reference like x[] means the +x-coord of each atom will be used when evaluating the variable.

-

Compute references access global scalar or vector quantities or -per-atom scalar or vector quantities calculated by a +

Compute references access one or more quantities calculated by a compute. The ID in the reference should be replaced by the actual ID of the compute defined elsewhere in the input script. -See the compute command for details. Any specific -compute will either generate global or per-atom values, so there -is no amiguity as to which is used. +See the doc pages for individual computes to see which ones calculate +global versus per-atom quantities. If the compute reference contains +empty brackets, then per-atom values calculated by the compute are +accessed. Otherwise a single value (global or per-atom) calculated by +the compute is accessed.

-

If c_ID is used as a keyword, then the scalar quantity calculated by -the compute is used (global or per-atom). If c_ID[N] is used, -then one component of the vector quantity calculated by the compute is -used (global or per-atom). N should be an integer from 1-M, where M -is the length of the vector calculated by the compute. +

The different kinds of compute references are as follows. M is a +positive integer <= the number of vector values calculated by the +compute. N is a global atom ID (positive integer).

-

Fix references access global scalar or vector quantities or per-atom -scalar or vector quantities calculated by a fix. The ID in -the reference should be replaced by the actual ID of the fix defined -elsewhere in the input script. See the doc pages for individual fixes -to see which ones compute what kind of quantities. Any specific fix -compute will either generate global or per-atom values, so there is no -amiguity as to which is used. Some fixes only generate quantities on -certain timesteps. If a variable attempts to access the fix on -non-allowed timesteps, an error is generated. For example, the fix -ave/time command may only generate averaged -quantities every 100 steps. See the fix command for -details. +

+ + + + + +
c_ID scalar value of a global compute
c_ID[2] vector component of a global compute
c_ID[N] single atom's scalar value of a per-atom compute
c_ID[N][M] single atom's vector component of a per-atom compute
c_ID[] per-atom scalar from a per-atom compute
c_ID[][M] per-atom vector component from a per-atom compute +
+ +

Fix references access one or more quantities calculated by a +fix. The ID in the reference should be replaced by +the actual ID of the fix defined elsewhere in the input script. +See the doc pages for individual computes to see which ones calculate +global versus per-atom quantities. If the compute reference contains +empty brackets, then per-atom values calculated by the compute are +accessed. Otherwise a single value (global or per-atom) calculated by +the compute is accessed.

-

If f_ID is used as a keyword, then the scalar quantity calculated by -the fix is used (global or per-atom). If f_ID[N] is used, then -one component of the vector quantity calculated by the fix is used -(global or per-atom). N should be an integer from 1-M, where M is the -length of the vector calculated by the fix. +

Note that some fixes only generate quantities on certain timesteps. +If a variable attempts to access the fix on non-allowed timesteps, an +error is generated. For example, the fix ave/time +command may only generate averaged quantities every 100 steps. See +the doc pages for individual fix commands for details. +

+

The different kinds of fix references are exactly the same as the +compute references listed in the above table, where "c_" is replaced +by "f_", and the word "compute" is replaced by "fix".

The current values of other variables can be accessed by prepending a "v_" to the variable name. This will cause that variable to be -evaulated. Note that equal style variables generate single values -and atom style variables generate one value per atom. Thus an -atom-style variable cannot be referenced by the formula for an -equal-style variable, but the converse is allowed. +evaulated. Atom-style variables generate per-atom values; all other +styles of variables generate a single scalar value.

+

The different kinds of variable references are as follows. N is a +global atom ID (positive integer). +

+
+ + +
v_ID scalar value of a non atom-style variable
v_ID[N] single atom's scalar value from an atom-style variable
v_ID[] per-atom value from an atom-style variable +
+

IMPORTANT NOTE: If you define variables in circular manner like this: 

variable a equal v_b
@@ -359,11 +375,11 @@ the first case, but would change dynamically for the second case.
 
 
Restrictions:
 

-
The use of atom values in equal or atom style variables requires
-the atom style to use a global mapping in order to look up the vector
-indices.  By default, only atom styles with molecular information
-create global maps.  The atom_modify map command
-can override the default.
+
Indexing any formula element by global atom ID, such as an atom value,
+requires the atom style to use a global mapping in order to look up
+the vector indices.  By default, only atom styles with molecular
+information create global maps.  The atom_modify
+map command can override the default.
 

 
All universe- and uloop-style variables defined in an input script
 must have the same number of values.
diff --git a/doc/variable.txt b/doc/variable.txt
index 097128a7f2..112e29eda6 100644
--- a/doc/variable.txt
+++ b/doc/variable.txt
@@ -20,7 +20,7 @@ style = {index} or {loop} or {world} or {universe} or {uloop} or {equal} or {ato
   {universe} args = one or more strings
   {uloop} args = N = integer size of loop
   {equal} or {atom} args = one formula containing numbers, thermo keywords, math operations, group functions, atom values and vectors, compute and fix references, other variables
-    numbers = 0.0, -5.4, 2.8e-4, etc
+    numbers = 0.0, 100, -5.4, 2.8e-4, etc
     thermo keywords = vol, ke, press, etc from "thermo_style"_thermo_style.html
     math operations = (), -x, x+y, x-y, x*y, x/y, x^y, sqrt(x), exp(x), ln(x)
     group functions = count(group), mass(group), charge(group),
@@ -30,9 +30,9 @@ style = {index} or {loop} or {world} or {universe} or {uloop} or {equal} or {ato
                  vx\[N\], vy\[N\], vz\[N\], fx\[N\], fy\[N\], fz\[N\]
     atom vector = mass\[\], x\[\], y\[\], z\[\],
                   vx\[\], vy\[\], vz\[\], fx\[\], fy\[\], fz\[\]
-    compute references = c_ID, c_ID\[N\]
-    fix references = f_ID, f_ID\[N\]
-    other variables = v_abc, v_x, etc :pre
+    compute references = c_ID, c_ID\[2\], c_ID\[N\], c_ID\[N\]\[2\], c_ID\[\], c_ID\[\]\[2\]
+    fix references = f_ID, f_ID\[2\], f_ID\[N\], f_ID\[N\]\[2\], f_ID\[\], f_ID\[\]\[2\]
+    other variables = v_abc, v_abc\[N\], v_abc\[\] :pre
 :ule
 
 [Examples:]
@@ -203,7 +203,7 @@ Specifically, an formula can contain numbers, thermo keywords, math
 operations, group functions, atom values, atom vectors, compute
 references, fix references, and references to other variables.
 
-Number: 0.2, 1.0e20, -15.4, etc
+Number: 0.2, 100, 1.0e20, -15.4, etc
 Thermo keywords: vol, pe, ebond, etc
 Math operations: (), -x, x+y, x-y, x*y, x/y, x^y, sqrt(x), exp(x), ln(x)
 Group functions: count(ID), mass(ID), charge(ID), xcm(ID,dim), \
@@ -212,17 +212,18 @@ Atom values: mass\[N\], x\[N\], y\[N\], z\[N\], \
              vx\[N\], vy\[N\], vz\[N\], fx\[N\], fy\[N\], fz\[N\]
 Atom vectors: mass\[\], x\[\], y\[\], z\[\], \
               vx\[\], vy\[\], vz\[\], fx\[\], fy\[\], fz\[\]
-Compute references: c_ID, c_ID\[N\]
-Fix references: f_ID, f_ID\[N\]
-Other variables: v_abc, v_x, etc :tb(s=:)
+Compute references: c_ID, c_ID\[2\], c_ID\[N\], c_ID\[N\]\[2\], c_ID\[\], c_ID\[\]\[2\]
+Fix references: f_ID, f_ID\[2\], f_ID\[N\], f_ID\[N\]\[2\], f_ID\[\], f_ID\[\]\[2\]
+Other variables: v_abc, v_abc\[N\], v_abc\[\] :tb(s=:)
 
-Note that a formula for equal-style variables cannot use any formula
-element that produces per-atom values.  This includes atom vectors, a
-compute that produces per-atom values, a fix that produces per-atoms
-values, or an atom-style variable.  A formula for an atom-style
-variables can use these elements, as well as elements that produce a
-global scalar or vector, e.g. a compute or fix that produces global
-quantities.
+Note that formula elements that contain empty brackets, such as an
+atom vector, produce per-atom values.  All other formula elements
+produce a global value.
+
+A formula for equal-style variables cannot use any formula element
+that produces per-atom values.  A formula for an atom-style variable
+can use formula elements that produce either global values or per-atom
+values.
 
 The thermo keywords allowed in a formula are those defined by the
 "thermo_style custom" command.  Since many thermodyanmic quantities
@@ -257,49 +258,60 @@ Atom values take a single integer argument from 1-N, which is the
 desired atom-ID, e.g. x\[243\]., which means use the x coordinate of
 the atom with ID=243.
 
-Atom vectors take no argument.  They generate one value per atom, so
-that a reference like x\[\] means the x-coord of each atom will be
-used when evaluating the variable.
+Atom vectors use empty brackets, i.e. they take no argument.  They
+generate one value per atom, so that a reference like x\[\] means the
+x-coord of each atom will be used when evaluating the variable.
 
-Compute references access global scalar or vector quantities or
-per-atom scalar or vector quantities calculated by a
+Compute references access one or more quantities calculated by a
 "compute"_compute.html.  The ID in the reference should be replaced by
 the actual ID of the compute defined elsewhere in the input script.
-See the "compute"_compute.html command for details.  Any specific
-compute will either generate global or per-atom values, so there
-is no amiguity as to which is used.
+See the doc pages for individual computes to see which ones calculate
+global versus per-atom quantities.  If the compute reference contains
+empty brackets, then per-atom values calculated by the compute are
+accessed.  Otherwise a single value (global or per-atom) calculated by
+the compute is accessed.
 
-If {c_ID} is used as a keyword, then the scalar quantity calculated by
-the compute is used (global or per-atom).  If {c_ID\[N\]} is used,
-then one component of the vector quantity calculated by the compute is
-used (global or per-atom).  N should be an integer from 1-M, where M
-is the length of the vector calculated by the compute.
+The different kinds of compute references are as follows.  M is a
+positive integer <= the number of vector values calculated by the
+compute.  N is a global atom ID (positive integer).
 
-Fix references access global scalar or vector quantities or per-atom
-scalar or vector quantities calculated by a "fix"_fix.html.  The ID in
-the reference should be replaced by the actual ID of the fix defined
-elsewhere in the input script.  See the doc pages for individual fixes
-to see which ones compute what kind of quantities.  Any specific fix
-compute will either generate global or per-atom values, so there is no
-amiguity as to which is used.  Some fixes only generate quantities on
-certain timesteps.  If a variable attempts to access the fix on
-non-allowed timesteps, an error is generated.  For example, the "fix
-ave/time"_fix_ave_time.html command may only generate averaged
-quantities every 100 steps.  See the "fix"_fix.html command for
-details.
+c_ID: scalar value of a global compute
+c_ID\[2\]: vector component of a global compute
+c_ID\[N\]: single atom's scalar value of a per-atom compute
+c_ID\[N\]\[M\]: single atom's vector component of a per-atom compute
+c_ID\[\]: per-atom scalar from a per-atom compute
+c_ID\[\]\[M\]: per-atom vector component from a per-atom compute :tb(s=:)
 
-If {f_ID} is used as a keyword, then the scalar quantity calculated by
-the fix is used (global or per-atom).  If {f_ID\[N\]} is used, then
-one component of the vector quantity calculated by the fix is used
-(global or per-atom).  N should be an integer from 1-M, where M is the
-length of the vector calculated by the fix.
+Fix references access one or more quantities calculated by a
+"fix"_fix.html.  The ID in the reference should be replaced by
+the actual ID of the fix defined elsewhere in the input script.
+See the doc pages for individual computes to see which ones calculate
+global versus per-atom quantities.  If the compute reference contains
+empty brackets, then per-atom values calculated by the compute are
+accessed.  Otherwise a single value (global or per-atom) calculated by
+the compute is accessed.
+
+Note that some fixes only generate quantities on certain timesteps.
+If a variable attempts to access the fix on non-allowed timesteps, an
+error is generated.  For example, the "fix ave/time"_fix_ave_time.html
+command may only generate averaged quantities every 100 steps.  See
+the doc pages for individual fix commands for details.
+
+The different kinds of fix references are exactly the same as the
+compute references listed in the above table, where "c_" is replaced
+by "f_", and the word "compute" is replaced by "fix".
 
 The current values of other variables can be accessed by prepending a
 "v_" to the variable name.  This will cause that variable to be
-evaulated.  Note that {equal} style variables generate single values
-and {atom} style variables generate one value per atom.  Thus an
-atom-style variable cannot be referenced by the formula for an
-equal-style variable, but the converse is allowed.
+evaulated.  Atom-style variables generate per-atom values; all other
+styles of variables generate a single scalar value.
+
+The different kinds of variable references are as follows.  N is a
+global atom ID (positive integer).
+
+v_ID: scalar value of a non atom-style variable
+v_ID\[N\]: single atom's scalar value from an atom-style variable
+v_ID\[\]: per-atom value from an atom-style variable :tb(s=:)
 
 IMPORTANT NOTE: If you define variables in circular manner like this:
 
@@ -354,11 +366,11 @@ the first case, but would change dynamically for the second case.
 
 [Restrictions:]
 
-The use of atom values in {equal} or {atom} style variables requires
-the atom style to use a global mapping in order to look up the vector
-indices.  By default, only atom styles with molecular information
-create global maps.  The "atom_modify map"_atom_modify.html command
-can override the default.
+Indexing any formula element by global atom ID, such as an atom value,
+requires the atom style to use a global mapping in order to look up
+the vector indices.  By default, only atom styles with molecular
+information create global maps.  The "atom_modify
+map"_atom_modify.html command can override the default.
 
 All {universe}- and {uloop}-style variables defined in an input script
 must have the same number of values.
diff --git a/src/variable.cpp b/src/variable.cpp
index 14d1d73716..f315ceb387 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -494,7 +494,7 @@ void Variable::copy(int narg, char **from, char **to)
 }
 
 /* ----------------------------------------------------------------------
-   recursive evaluation of a "string"
+   recursive evaluation of a string str
    string is a equal-style or atom-style formula containing one or more items:
      number = 0.0, -5.45, 2.8e-4, ...
      thermo keyword = ke, vol, atoms, ...
@@ -502,9 +502,9 @@ void Variable::copy(int narg, char **from, char **to)
      group function = count(group), mass(group), xcm(group,x), ...
      atom value = x[123], y[3], vx[34], ...
      atom vector = x[], y[], vx[], ...
-     compute = c_mytemp, c_thermo_pressure[3], ...
-     fix = f_indent, f_setforce[3], ...
-     variable = v_a, v_myvar, ...
+     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]
    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
 ------------------------------------------------------------------------- */
@@ -513,6 +513,7 @@ double Variable::evaluate(char *str, Tree **tree)
 {
   int op,opprevious;
   double value1,value2;
+  char onechar;
 
   double argstack[MAXLEVEL];
   Tree *treestack[MAXLEVEL];
@@ -525,13 +526,15 @@ double Variable::evaluate(char *str, Tree **tree)
   int expect = ARG;
 
   while (1) {
-    char onechar = str[i];
+    onechar = str[i];
 
     // whitespace: just skip
 
     if (isspace(onechar)) i++;
 
+    // ----------------
     // parentheses: recursively evaluate contents of parens
+    // ----------------
 
     else if (onechar == '(') {
       if (expect == OP) error->all("Invalid syntax in variable formula");
@@ -551,7 +554,9 @@ double Variable::evaluate(char *str, Tree **tree)
 
       delete [] contents;
 
+    // ----------------
     // number: push value onto stack
+    // ----------------
 
     } else if (isdigit(onechar) || onechar == '.') {
       if (expect == OP) error->all("Invalid syntax in variable formula");
@@ -578,13 +583,16 @@ double Variable::evaluate(char *str, Tree **tree)
 
       delete [] number;
 
+    // ----------------
     // letter: c_ID, f_ID, v_name, exp(), xcm(,), x[234], x[], vol
+    // ----------------
 
     } else if (islower(onechar)) {
       if (expect == OP) error->all("Invalid syntax in variable formula");
       expect = OP;
 
-      // istop = end of word, assumed to be alphanumeric or underscore
+      // istop = end of word
+      // word = all alphanumeric or underscore
 
       int istart = i;
       while (isalnum(str[i]) || str[i] == '_') i++;
@@ -595,7 +603,9 @@ double Variable::evaluate(char *str, Tree **tree)
       strncpy(word,&str[istart],n);
       word[n] = '\0';
 
+      // ----------------
       // compute
+      // ----------------
 
       if (strncmp(word,"c_",2) == 0) {
 	n = strlen(word) - 2 + 1;
@@ -605,78 +615,119 @@ double Variable::evaluate(char *str, Tree **tree)
 	if (update->first_update == 0)
 	  error->all("Compute in variable formula before initial run");
 
-	int icompute;
-	for (icompute = 0; icompute < modify->ncompute; icompute++)
-	  if (strcmp(id,modify->compute[icompute]->id) == 0) break;
-	if (icompute == modify->ncompute)
-	  error->all("Invalid compute ID in variable formula");
+	int icompute = modify->find_compute(id);
+	if (icompute < 0) error->all("Invalid compute ID in variable formula");
 	Compute *compute = modify->compute[icompute];
 
 	delete [] id;
 
-	int index = 0;
-	if (str[i] == '[') {
-	  i = int_between_brackets(str,i,index,0);
+	// 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
+
+	int nbracket,index1,index2;
+	if (str[i] != '[') nbracket = 0;
+	else {
+	  nbracket = 1;
+	  i = int_between_brackets(str,i,index1,1);
 	  i++;
+	  if (str[i] == '[') {
+	    nbracket = 2;
+	    i = int_between_brackets(str,i,index2,0);
+	    i++;
+	  }
 	}
 
-	// global or per-atom scalar compute
+        // c_ID = global scalar
 
-	if (index == 0) {
-	  if (compute->scalar_flag) {
-	    if (!(compute->invoked & INVOKED_SCALAR))
-	      value1 = compute->compute_scalar();
-	    else value1 = compute->scalar;
-	    if (tree) {
-	      Tree *newtree = new Tree();
-	      newtree->type = VALUE;
-	      newtree->value = value1;
-	      treestack[ntreestack++] = newtree;
-	    } else argstack[nargstack++] = value1;
+	if (nbracket == 0 && compute->scalar_flag) {
 
-	  } else if (compute->peratom_flag && compute->size_peratom == 0) {
-	    if (tree == NULL)
-	      error->all("Per-atom compute in equal-style variable formula");
-	    if (!(compute->invoked & INVOKED_PERATOM))
-	      compute->compute_peratom();
+	  if (!(compute->invoked & INVOKED_SCALAR))
+	    value1 = compute->compute_scalar();
+	  else value1 = compute->scalar;
+	  if (tree) {
 	    Tree *newtree = new Tree();
-	    newtree->type = ATOMARRAY;
-	    newtree->array = compute->scalar_atom;
-	    newtree->nstride = 1;
+	    newtree->type = VALUE;
+	    newtree->value = value1;
 	    treestack[ntreestack++] = newtree;
-	  } else error->all("Mismatched compute in variable formula");
+	  } else argstack[nargstack++] = value1;
 
-	// global or per-atom vector compute
+        // c_ID[2] = global vector
 
-	} else {
-	  if (compute->vector_flag) {
-	    if (index > compute->size_vector)
+	} else if (nbracket == 1 && index1 > 0 && compute->vector_flag) {
+
+	  if (index1 > compute->size_vector)
 	      error->all("Compute vector in variable formula is too small");
-	    if (!(compute->invoked & INVOKED_VECTOR)) compute->compute_vector();
-	    value1 = compute->vector[index-1];
-	    if (tree) {
-	      Tree *newtree = new Tree();
-	      newtree->type = VALUE;
-	      newtree->value = value1;
-	      treestack[ntreestack++] = newtree;
-	    } else argstack[nargstack++] = value1;
-
-	  } else if (compute->peratom_flag && compute->size_peratom) {
-	    if (tree == NULL)
-	      error->all("Per-atom compute in equal-style variable formula");
-	    if (index > compute->size_peratom)
-	      error->all("Compute vector in variable formula is too small");
-	    if (!(compute->invoked & INVOKED_PERATOM))
-	      compute->compute_peratom();
+	  if (!(compute->invoked & INVOKED_VECTOR)) compute->compute_vector();
+	  value1 = compute->vector[index1-1];
+	  if (tree) {
 	    Tree *newtree = new Tree();
-	    newtree->type = ATOMARRAY;
-	    newtree->array = &compute->vector_atom[0][index-1];
-	    newtree->nstride = compute->size_peratom;
+	    newtree->type = VALUE;
+	    newtree->value = value1;
 	    treestack[ntreestack++] = newtree;
-	  } else error->all("Mismatched compute in variable formula");
-	}
+	  } else argstack[nargstack++] = value1;
 
+        // c_ID[] = per-atom scalar
+
+	} else if (nbracket == 1 && index1 == 0 && 
+		   compute->peratom_flag && compute->size_peratom == 0) {
+
+	  if (tree == NULL)
+	    error->all("Per-atom compute in equal-style variable formula");
+	  if (!(compute->invoked & INVOKED_PERATOM))
+	    compute->compute_peratom();
+	  Tree *newtree = new Tree();
+	  newtree->type = ATOMARRAY;
+	  newtree->array = compute->scalar_atom;
+	  newtree->nstride = 1;
+	  treestack[ntreestack++] = newtree;
+
+        // c_ID[N] = global value from per-atom scalar
+
+	} else if (nbracket == 1 && index1 > 0 && 
+		   compute->peratom_flag && compute->size_peratom == 0) {
+
+	  if (!(compute->invoked & INVOKED_PERATOM))
+	    compute->compute_peratom();
+	  peratom2global(1,NULL,compute->scalar_atom,1,index1,
+			 tree,treestack,ntreestack,argstack,nargstack);
+
+        // c_ID[][2] = per-atom vector
+
+	} else if (nbracket == 2 && index1 == 0 && index2 > 0 &&
+		   compute->peratom_flag) {
+
+	  if (tree == NULL)
+	    error->all("Per-atom compute in equal-style variable formula");
+	  if (index2 > compute->size_peratom)
+	    error->all("Compute vector in variable formula is too small");
+	  if (!(compute->invoked & INVOKED_PERATOM))
+	    compute->compute_peratom();
+	  Tree *newtree = new Tree();
+	  newtree->type = ATOMARRAY;
+	  newtree->array = &compute->vector_atom[0][index2-1];
+	  newtree->nstride = compute->size_peratom;
+	  treestack[ntreestack++] = newtree;
+
+        // c_ID[N][2] = global value from per-atom vector
+
+	} else if (nbracket == 2 && index1 > 0 && index2 > 0 &&
+		   compute->peratom_flag) {
+
+	  if (index2 > compute->size_peratom)
+	    error->all("Compute vector in variable formula is too small");
+	  if (!(compute->invoked & INVOKED_PERATOM))
+	    compute->compute_peratom();
+	  peratom2global(1,NULL,&compute->vector_atom[0][index2-1],
+			 compute->size_peratom,index1,
+			 tree,treestack,ntreestack,argstack,nargstack);
+
+	} else error->all("Mismatched compute in variable formula");
+
+      // ----------------
       // fix
+      // ----------------
 
       } else if (strncmp(word,"f_",2) == 0) {
 	n = strlen(word) - 2 + 1;
@@ -686,77 +737,120 @@ double Variable::evaluate(char *str, Tree **tree)
 	if (update->first_update == 0)
 	  error->all("Fix in variable formula before initial run");
 
-	int ifix;
-	for (ifix = 0; ifix < modify->nfix; ifix++)
-	  if (strcmp(id,modify->fix[ifix]->id) == 0) break;
-	if (ifix == modify->nfix)
-	  error->all("Invalid fix ID in variable formula");
+	int ifix = modify->find_fix(id);
+	if (ifix < 0) error->all("Invalid fix ID in variable formula");
 	Fix *fix = modify->fix[ifix];
 
 	delete [] id;
 
-	int index = 0;
-	if (str[i] == '[') {
-	  i = int_between_brackets(str,i,index,0);
+	// 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
+
+	int nbracket,index1,index2;
+	if (str[i] != '[') nbracket = 0;
+	else {
+	  nbracket = 1;
+	  i = int_between_brackets(str,i,index1,1);
 	  i++;
+	  if (str[i] == '[') {
+	    nbracket = 2;
+	    i = int_between_brackets(str,i,index2,0);
+	    i++;
+	  }
 	}
 
-	// global or per-atom scalar fix
+        // f_ID = global scalar
 
-	if (index == 0) {
-	  if (fix->scalar_flag) {
-	    if (update->ntimestep % fix->scalar_vector_freq)
-	      error->all("Fix in variable not computed at compatible time");
-	    value1 = fix->compute_scalar();
-	    if (tree) {
-	      Tree *newtree = new Tree();
-	      newtree->type = VALUE;
-	      newtree->value = value1;
-	      treestack[ntreestack++] = newtree;
-	    } else argstack[nargstack++] = value1;
+	if (nbracket == 0 && fix->scalar_flag) {
 
-	  } else if (fix->peratom_flag && fix->size_peratom == 0) {
-	    if (tree == NULL)
-	      error->all("Per-atom fix in equal-style variable formula");
-	    if (update->ntimestep % fix->peratom_freq)
-	      error->all("Fix in variable not computed at compatible time");
+	  if (update->ntimestep % fix->scalar_vector_freq)
+	    error->all("Fix in variable not computed at compatible time");
+	  value1 = fix->compute_scalar();
+	  if (tree) {
 	    Tree *newtree = new Tree();
-	    newtree->type = ATOMARRAY;
-	    newtree->array = fix->scalar_atom;
-	    newtree->nstride = 1;
+	    newtree->type = VALUE;
+	    newtree->value = value1;
 	    treestack[ntreestack++] = newtree;
-	  } else error->all("Mismatched fix in variable formula");
+	  } else argstack[nargstack++] = value1;
 
-	// global or per-atom vector fix
+        // f_ID[2] = global vector
 
-	} else {
-	  if (fix->vector_flag) {
-	    if (index > fix->size_vector)
+	} else if (nbracket == 1 && index1 > 0 && fix->vector_flag) {
+
+	  if (index1 > fix->size_vector)
 	      error->all("Fix vector in variable formula is too small");
-	    value1 = fix->compute_vector(index-1);
-	    if (tree) {
-	      Tree *newtree = new Tree();
-	      newtree->type = VALUE;
-	      newtree->value = value1;
-	      treestack[ntreestack++] = newtree;
-	    } else argstack[nargstack++] = value1;
-
-	  } else if (fix->peratom_flag && fix->size_peratom) {
-	    if (tree == NULL)
-	      error->all("Per-atom fix in equal-style variable formula");
-	    if (index > fix->size_peratom)
-	      error->all("Fix vector in variable formula is too small");
-	    if (update->ntimestep % fix->peratom_freq)
-	      error->all("Fix in variable not computed at compatible time");
+	  if (update->ntimestep % fix->scalar_vector_freq)
+	    error->all("Fix in variable not computed at compatible time");
+	  value1 = fix->compute_vector(index1-1);
+	  if (tree) {
 	    Tree *newtree = new Tree();
-	    newtree->type = ATOMARRAY;
-	    newtree->array = &fix->vector_atom[0][index-1];
-	    newtree->nstride = fix->size_peratom;
+	    newtree->type = VALUE;
+	    newtree->value = value1;
 	    treestack[ntreestack++] = newtree;
-	  } else error->all("Mismatched fix in variable formula");
-	}
+	  } else argstack[nargstack++] = value1;
 
+        // f_ID[] = per-atom scalar
+
+	} else if (nbracket == 1 && index1 == 0 && 
+		   fix->peratom_flag && fix->size_peratom == 0) {
+
+	  if (tree == NULL)
+	    error->all("Per-atom fix in equal-style variable formula");
+	  if (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->scalar_atom;
+	  newtree->nstride = 1;
+	  treestack[ntreestack++] = newtree;
+
+        // f_ID[N] = global value from per-atom scalar
+
+	} else if (nbracket == 1 && index1 > 0 && 
+		   fix->peratom_flag && fix->size_peratom == 0) {
+
+	  if (update->ntimestep % fix->peratom_freq)
+	    error->all("Fix in variable not computed at compatible time");
+	  peratom2global(1,NULL,fix->scalar_atom,1,index1,
+			 tree,treestack,ntreestack,argstack,nargstack);
+
+        // f_ID[][2] = per-atom vector
+
+	} else if (nbracket == 2 && index1 == 0 && index2 > 0 &&
+		   fix->peratom_flag) {
+
+	  if (tree == NULL)
+	    error->all("Per-atom fix in equal-style variable formula");
+	  if (index2 > fix->size_peratom)
+	    error->all("Fix vector in variable formula is too small");
+	  if (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[0][index2-1];
+	  newtree->nstride = fix->size_peratom;
+	  treestack[ntreestack++] = newtree;
+
+        // f_ID[N][2] = global value from per-atom vector
+
+	} else if (nbracket == 2 && index1 > 0 && index2 > 0 &&
+		   fix->peratom_flag) {
+
+	  if (index2 > fix->size_peratom)
+	    error->all("Fix vector in variable formula is too small");
+	  if (update->ntimestep % fix->peratom_freq)
+	    error->all("Fix in variable not computed at compatible time");
+	  peratom2global(1,NULL,&fix->vector_atom[0][index2-1],
+			 fix->size_peratom,index1,
+			 tree,treestack,ntreestack,argstack,nargstack);
+
+	} else error->all("Mismatched fix in variable formula");
+
+      // ----------------
       // variable
+      // ----------------
 
       } else if (strncmp(word,"v_",2) == 0) {
 	n = strlen(word) - 2 + 1;
@@ -764,29 +858,65 @@ double Variable::evaluate(char *str, Tree **tree)
 	strcpy(id,&word[2]);
 
 	int ivar = find(id);
-	if (ivar == -1) error->all("Invalid variable name in variable formula");
-	char *var = retrieve(id);
-	if (var == NULL && style[ivar] != ATOM)
-	  error->all("Invalid variable evaluation in variable formula");
+	if (ivar < 0) error->all("Invalid variable name in variable formula");
 
-	if (tree) {
-	  Tree *newtree;
-	  if (style[ivar] != ATOM) {
-	    newtree = new Tree();
+	// parse zero or one trailing brackets
+	// point i beyond last bracket
+	// nbracket = # of bracket pairs
+	// index = int inside bracket, 0 if empty
+
+	int nbracket,index;
+	if (str[i] != '[') nbracket = 0;
+	else {
+	  nbracket = 1;
+	  i = int_between_brackets(str,i,index,1);
+	  i++;
+	}
+
+        // v_name = non atom-style variable = global value
+
+	if (nbracket == 0 && style[ivar] != ATOM) {
+
+	  char *var = retrieve(id);
+	  if (var == NULL)
+	    error->all("Invalid variable evaluation in variable formula");
+	  if (tree) {
+	    Tree *newtree = new Tree();
 	    newtree->type = VALUE;
 	    newtree->value = atof(var);
-	  } else double tmp = evaluate(data[ivar][0],&newtree);
+	    treestack[ntreestack++] = newtree;
+	  } else argstack[nargstack++] = atof(var);
+
+        // v_name[] = atom-style variable
+
+	} else if (nbracket && index == 0 && style[ivar] == ATOM) {
+
+	  if (tree == NULL)
+	    error->all("Atom-style variable in equal-style variable formula");
+	  Tree *newtree;
+	  double tmp = evaluate(data[ivar][0],&newtree);
 	  treestack[ntreestack++] = newtree;
 
-	} else {
-	  if (style[ivar] == ATOM)
-	    error->all("Atom-style variable in equal-style variable formula");
-	  argstack[nargstack++] = atof(var);
-	}
+        // v_name[N] = global value from atom-style variable
+	// compute the per-atom variable in result
+	// use peratom2global to extract single value from result
+
+	} else if (nbracket && index > 0 && style[ivar] == ATOM) {
+
+	  double *result = (double *)
+	    memory->smalloc(atom->nlocal*sizeof(double),"variable:result");
+	  compute_atom(ivar,0,result,1,0);
+	  peratom2global(1,NULL,result,1,index,
+			 tree,treestack,ntreestack,argstack,nargstack);
+	  memory->sfree(result);
+
+	} else error->all("Mismatched variable in variable formula");
 
 	delete [] id;
 
+      // ----------------
       // math/group function or atom value/vector or thermo keyword
+      // ----------------
 
       } else {
 
@@ -805,22 +935,27 @@ double Variable::evaluate(char *str, Tree **tree)
 
 	  delete [] contents;
 
+	// ----------------
 	// atom value or vector
+	// ----------------
 
 	} else if (str[i] == '[') {
 	  int id;
 	  i = int_between_brackets(str,i,id,1);
 	  i++;
 
-	  // ID between brackets exists, atom value
-	  // empty brackets, atom vector
+	  // ID between brackets exists: atom value
+	  // empty brackets: atom vector
 
 	  if (id > 0)
-	    atom_value(word,id,tree,treestack,ntreestack,argstack,nargstack);
+	    peratom2global(0,word,NULL,0,id,
+			   tree,treestack,ntreestack,argstack,nargstack);
 	  else
 	    atom_vector(word,tree,treestack,ntreestack);
 
+	// ----------------
 	// thermo keyword
+	// ----------------
 
 	} else {
 	  if (update->first_update == 0)
@@ -838,7 +973,9 @@ double Variable::evaluate(char *str, Tree **tree)
 
       delete [] word;
 
+    // ----------------
     // math operator, including end-of-string
+    // ----------------
 
     } else if (strchr("+-*/^\0",onechar)) {
       if (onechar == '+') op = ADD;
@@ -1001,8 +1138,8 @@ int Variable::find_matching_paren(char *str, int i,char *&contents)
 
 /* ----------------------------------------------------------------------
    find int between brackets and set index to its value
-   if emptyflag and brackets empty, set index to 0
-   else if brackets empty, is an error
+   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
@@ -1211,40 +1348,45 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
 }
 
 /* ----------------------------------------------------------------------
-   process an atom value in formula
+   extract a global value from a per-atom quantity in a formula
+   flag = 0 -> word is an atom vector
+   flag = 1 -> vector is a per-atom compute or fix quantity with nstride
+   id = positive global ID of atom, converted to local index
    push result onto tree or arg stack
-   word = atom vector
-   index = positive global ID of atom
    customize by adding an atom vector: mass,x,y,z,vx,vy,vz,fx,fy,fz
 ------------------------------------------------------------------------- */
 
-void Variable::atom_value(char *word, int id, Tree **tree,
-			  Tree **treestack, int &ntreestack,
-			  double *argstack, int &nargstack)
+void Variable::peratom2global(int flag, char *word,
+			      double *vector, int nstride, int id,
+			      Tree **tree, Tree **treestack, int &ntreestack,
+			      double *argstack, int &nargstack)
 {
   if (atom->map_style == 0)
-    error->all("Atom value in variable formula without atom map");
+    error->all("Indexed per-atom vector in variable formula without atom map");
 
   int index = atom->map(id);
 
   double mine;
   if (index >= 0 && index < atom->nlocal) {
-	      
-    if (strcmp(word,"mass") == 0) {
-      if (atom->mass) mine = atom->mass[atom->type[index]];
-      else mine = atom->rmass[index];
-    }
-    else if (strcmp(word,"x") == 0) mine = atom->x[index][0];
-    else if (strcmp(word,"y") == 0) mine = atom->x[index][1];
-    else if (strcmp(word,"z") == 0) mine = atom->x[index][2];
-    else if (strcmp(word,"vx") == 0) mine = atom->v[index][0];
-    else if (strcmp(word,"vy") == 0) mine = atom->v[index][1];
-    else if (strcmp(word,"vz") == 0) mine = atom->v[index][2];
-    else if (strcmp(word,"fx") == 0) mine = atom->f[index][0];
-    else if (strcmp(word,"fy") == 0) mine = atom->f[index][1];
-    else if (strcmp(word,"fz") == 0) mine = atom->f[index][2];
-    
-    else error->one("Invalid atom vector in variable formula");
+
+    if (flag == 0) {
+      if (strcmp(word,"mass") == 0) {
+	if (atom->mass) mine = atom->mass[atom->type[index]];
+	else mine = atom->rmass[index];
+      }
+      else if (strcmp(word,"x") == 0) mine = atom->x[index][0];
+      else if (strcmp(word,"y") == 0) mine = atom->x[index][1];
+      else if (strcmp(word,"z") == 0) mine = atom->x[index][2];
+      else if (strcmp(word,"vx") == 0) mine = atom->v[index][0];
+      else if (strcmp(word,"vy") == 0) mine = atom->v[index][1];
+      else if (strcmp(word,"vz") == 0) mine = atom->v[index][2];
+      else if (strcmp(word,"fx") == 0) mine = atom->f[index][0];
+      else if (strcmp(word,"fy") == 0) mine = atom->f[index][1];
+      else if (strcmp(word,"fz") == 0) mine = atom->f[index][2];
+      
+      else error->one("Invalid atom vector in variable formula");
+
+    } else mine = vector[index*nstride];
     
   } else mine = 0.0;
 
diff --git a/src/variable.h b/src/variable.h
index 332d0879ad..8b0a98f533 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -60,7 +60,8 @@ class Variable : protected Pointers {
   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 &);
-  void atom_value(char *, int, Tree **, Tree **, int &, double *, int &);
+  void peratom2global(int, char *, double *, int, int,
+		      Tree **, Tree **, int &, double *, int &);
   void atom_vector(char *, Tree **, Tree **, int &);
 };