/src/variable.cpp

Large files files are truncated, but you can click here to view the full file

/* ----------------------------------------------------------------------
  LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under 
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "ctype.h"
#include "unistd.h"
#include "variable.h"
#include "universe.h"
#include "atom.h"
#include "update.h"
#include "group.h"
#include "domain.h"
#include "region.h"
#include "modify.h"
#include "compute.h"
#include "fix.h"
#include "output.h"
#include "thermo.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;

#define VARDELTA 4
#define MAXLEVEL 4

#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)

#define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)

enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,EQUAL,ATOM};
enum{ARG,OP};

// customize by adding a function

enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
       NOT,EQ,NE,LT,LE,GT,GE,AND,OR,
       SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
       RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,
       VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,GRMASK,
       VALUE,ATOMARRAY,TYPEARRAY,INTARRAY};

// customize by adding a special function

enum{SUM,XMIN,XMAX,AVE,TRAP};

#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
#define INVOKED_ARRAY 4
#define INVOKED_PERATOM 8

#define BIG 1.0e20

/* ---------------------------------------------------------------------- */

Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
{
  MPI_Comm_rank(world,&me);

  nvar = maxvar = 0;
  names = NULL;
  style = NULL;
  num = NULL;
  which = NULL;
  pad = NULL;
  data = NULL;

  randomequal = NULL;
  randomatom = NULL;

  precedence[DONE] = 0;
  precedence[OR] = 1;
  precedence[AND] = 2;
  precedence[EQ] = precedence[NE] = 3;
  precedence[LT] = precedence[LE] = precedence[GT] = precedence[GE] = 4;
  precedence[ADD] = precedence[SUBTRACT] = 5;
  precedence[MULTIPLY] = precedence[DIVIDE] = 6;
  precedence[CARAT] = 7;
  precedence[UNARY] = precedence[NOT] = 8;

  PI = 4.0*atan(1.0);
}

/* ---------------------------------------------------------------------- */

Variable::~Variable()
{
  for (int i = 0; i < nvar; i++) {
    delete [] names[i];
    if (style[i] == LOOP || style[i] == ULOOP) delete [] data[i][0];
    else for (int j = 0; j < num[i]; j++) delete [] data[i][j];
    delete [] data[i];
  }
  memory->sfree(names);
  memory->destroy(style);
  memory->destroy(num);
  memory->destroy(which);
  memory->destroy(pad);
  memory->sfree(data);

  delete randomequal;
  delete randomatom;
}

/* ----------------------------------------------------------------------
   called by variable command in input script
------------------------------------------------------------------------- */

void Variable::set(int narg, char **arg)
{
  if (narg < 2) error->all("Illegal variable command");

  // DELETE
  // doesn't matter if variable no longer exists

  if (strcmp(arg[1],"delete") == 0) {
    if (narg != 2) error->all("Illegal variable command");
    if (find(arg[0]) >= 0) remove(find(arg[0]));
    return;

  // INDEX
  // num = listed args, which = 1st value, data = copied args

  } else if (strcmp(arg[1],"index") == 0) {
    if (narg < 3) error->all("Illegal variable command");
    if (find(arg[0]) >= 0) return;
    if (nvar == maxvar) extend();
    style[nvar] = INDEX;
    num[nvar] = narg - 2;
    which[nvar] = 0;
    pad[nvar] = 0;
    data[nvar] = new char*[num[nvar]];
    copy(num[nvar],&arg[2],data[nvar]);

  // LOOP
  // 1 arg + pad: num = N, which = 1st value, data = single string
  // 2 args + pad: num = N2, which = N1, data = single string

  } else if (strcmp(arg[1],"loop") == 0) {
    if (find(arg[0]) >= 0) return;
    if (nvar == maxvar) extend();
    style[nvar] = LOOP;
    int nfirst,nlast;
    if (narg == 3 || (narg == 4 && strcmp(arg[3],"pad") == 0)) {
      nfirst = 1;
      nlast = atoi(arg[2]);
      if (nlast <= 0) error->all("Illegal variable command");
      if (narg == 4 && strcmp(arg[3],"pad") == 0) {
	char digits[12];
	sprintf(digits,"%d",nlast);
	pad[nvar] = strlen(digits);
      } else pad[nvar] = 0;
    } else if (narg == 4 || (narg == 5 && strcmp(arg[4],"pad") == 0)) {
      nfirst = atoi(arg[2]);
      nlast = atoi(arg[3]);
      if (nfirst > nlast || nlast <= 0) error->all("Illegal variable command");
      if (narg == 5 && strcmp(arg[4],"pad") == 0) {
	char digits[12];
	sprintf(digits,"%d",nlast);
	pad[nvar] = strlen(digits);
      } else pad[nvar] = 0;
    } else error->all("Illegal variable command");
    num[nvar] = nlast;
    which[nvar] = nfirst-1;
    data[nvar] = new char*[1];
    data[nvar][0] = NULL;

  // WORLD
  // num = listed args, which = partition this proc is in, data = copied args
  // error check that num = # of worlds in universe

  } else if (strcmp(arg[1],"world") == 0) {
    if (narg < 3) error->all("Illegal variable command");
    if (find(arg[0]) >= 0) return;
    if (nvar == maxvar) extend();
    style[nvar] = WORLD;
    num[nvar] = narg - 2;
    if (num[nvar] != universe->nworlds)
      error->all("World variable count doesn't match # of partitions");
    which[nvar] = universe->iworld;
    pad[nvar] = 0;
    data[nvar] = new char*[num[nvar]];
    copy(num[nvar],&arg[2],data[nvar]);

  // UNIVERSE and ULOOP
  // for UNIVERSE: num = listed args, data = copied args
  // for ULOOP: num = N, data = single string
  // which = partition this proc is in
  // universe proc 0 creates lock file
  // error check that all other universe/uloop variables are same length

  } else if (strcmp(arg[1],"universe") == 0 || strcmp(arg[1],"uloop") == 0) {
    if (strcmp(arg[1],"universe") == 0) {
      if (narg < 3) error->all("Illegal variable command");
      if (find(arg[0]) >= 0) return;
      if (nvar == maxvar) extend();
      style[nvar] = UNIVERSE;
      num[nvar] = narg - 2;
      data[nvar] = new char*[num[nvar]];
      copy(num[nvar],&arg[2],data[nvar]);
    } else if (strcmp(arg[1],"uloop") == 0) {
      if (narg < 3 || narg > 4 || (narg == 4 && strcmp(arg[3],"pad") != 0))
	error->all("Illegal variable command");
      if (find(arg[0]) >= 0) return;
      if (nvar == maxvar) extend();
      style[nvar] = ULOOP;
      num[nvar] = atoi(arg[2]);
      data[nvar] = new char*[1];
      data[nvar][0] = NULL;
      if (narg == 4) {
	char digits[12];
	sprintf(digits,"%d",num[nvar]);
	pad[nvar] = strlen(digits);
      } else pad[nvar] = 0;
    }

    if (num[nvar] < universe->nworlds)
      error->all("Universe/uloop variable count < # of partitions");
    which[nvar] = universe->iworld;

    if (universe->me == 0) {
      FILE *fp = fopen("tmp.lammps.variable","w");
      fprintf(fp,"%d\n",universe->nworlds);
      fclose(fp);
    }

    for (int jvar = 0; jvar < nvar; jvar++)
      if (num[jvar] && (style[jvar] == UNIVERSE || style[jvar] == ULOOP) && 
	  num[nvar] != num[jvar])
	error->all("All universe/uloop variables must have same # of values");

  // STRING
  // remove pre-existing var if also style STRING (allows it to be reset)
  // num = 1, which = 1st value
  // data = 1 value, string to eval

  } else if (strcmp(arg[1],"string") == 0) {
    if (narg != 3) error->all("Illegal variable command");
    if (find(arg[0]) >= 0) {
      if (style[find(arg[0])] != STRING)
	error->all("Cannot redefine variable as a different style");
      remove(find(arg[0]));
    }
    if (nvar == maxvar) extend();
    style[nvar] = STRING;
    num[nvar] = 1;
    which[nvar] = 0;
    pad[nvar] = 0;
    data[nvar] = new char*[num[nvar]];
    copy(1,&arg[2],data[nvar]);
    
  // EQUAL
  // remove pre-existing var if also style EQUAL (allows it to be reset)
  // num = 2, which = 1st value
  // data = 2 values, 1st is string to eval, 2nd is filled on retrieval

  } else if (strcmp(arg[1],"equal") == 0) {
    if (narg != 3) error->all("Illegal variable command");
    if (find(arg[0]) >= 0) {
      if (style[find(arg[0])] != EQUAL)
	error->all("Cannot redefine variable as a different style");
      remove(find(arg[0]));
    }
    if (nvar == maxvar) extend();
    style[nvar] = EQUAL;
    num[nvar] = 2;
    which[nvar] = 0;
    pad[nvar] = 0;
    data[nvar] = new char*[num[nvar]];
    copy(1,&arg[2],data[nvar]);
    data[nvar][1] = NULL;
    
  // ATOM
  // remove pre-existing var if also style ATOM (allows it to be reset)
  // num = 1, which = 1st value
  // data = 1 value, string to eval

  } else if (strcmp(arg[1],"atom") == 0) {
    if (narg != 3) error->all("Illegal variable command");
    if (find(arg[0]) >= 0) {
      if (style[find(arg[0])] != ATOM)
	error->all("Cannot redefine variable as a different style");
      remove(find(arg[0]));
    }
    if (nvar == maxvar) extend();
    style[nvar] = ATOM;
    num[nvar] = 1;
    which[nvar] = 0;
    pad[nvar] = 0;
    data[nvar] = new char*[num[nvar]];
    copy(1,&arg[2],data[nvar]);
    
  } else error->all("Illegal variable command");

  // set name of variable
  // must come at end, since STRING/EQUAL/ATOM reset may have removed name
  // name must be all alphanumeric chars or underscores

  int n = strlen(arg[0]) + 1;
  names[nvar] = new char[n];
  strcpy(names[nvar],arg[0]);

  for (int i = 0; i < n-1; i++)
    if (!isalnum(names[nvar][i]) && names[nvar][i] != '_')
      error->all("Variable name must be alphanumeric or "
		 "underscore characters");
  nvar++;
}

/* ----------------------------------------------------------------------
   INDEX variable created by command-line argument
   make it INDEX rather than STRING so cannot be re-defined in input script
------------------------------------------------------------------------- */

void Variable::set(char *name, int narg, char **arg)
{
  char **newarg = new char*[2+narg];
  newarg[0] = name;
  newarg[1] = (char *) "index";
  for (int i = 0; i < narg; i++) newarg[2+i] = arg[i];
  set(2+narg,newarg);
  delete [] newarg;
}

/* ----------------------------------------------------------------------
   increment variable(s)
   return 0 if OK if successfully incremented
   return 1 if any variable is exhausted, free the variable to allow re-use
------------------------------------------------------------------------- */

int Variable::next(int narg, char **arg)
{
  int ivar;

  if (narg == 0) error->all("Illegal next command");

  // check that variables exist and are all the same style
  // exception: UNIVERSE and ULOOP variables can be mixed in same next command

  for (int iarg = 0; iarg < narg; iarg++) {
    ivar = find(arg[iarg]);
    if (ivar == -1) error->all("Invalid variable in next command");
    if (style[ivar] == ULOOP && style[find(arg[0])] == UNIVERSE) continue;
    else if (style[ivar] == UNIVERSE && style[find(arg[0])] == ULOOP) continue;
    else if (style[ivar] != style[find(arg[0])])
      error->all("All variables in next command must be same style");
  }

  // invalid styles STRING or EQUAL or WORLD or ATOM

  int istyle = style[find(arg[0])];
  if (istyle == STRING || istyle == EQUAL || istyle == WORLD || istyle == ATOM)
    error->all("Invalid variable style with next command");

  // increment all variables in list
  // if any variable is exhausted, set flag = 1 and remove var to allow re-use

  int flag = 0;

  if (istyle == INDEX || istyle == LOOP) {
    for (int iarg = 0; iarg < narg; iarg++) {
      ivar = find(arg[iarg]);
      which[ivar]++;
      if (which[ivar] >= num[ivar]) {
	flag = 1;
	remove(ivar);
      }
    }

  } else if (istyle == UNIVERSE || istyle == ULOOP) {

    // wait until lock file can be created and owned by proc 0 of this world
    // read next available index and Bcast it within my world
    // set all variables in list to nextindex

    int nextindex;
    if (me == 0) {
      while (1) {
	if (!rename("tmp.lammps.variable","tmp.lammps.variable.lock")) break;
	usleep(100000);
      }
      FILE *fp = fopen("tmp.lammps.variable.lock","r");
      fscanf(fp,"%d",&nextindex);
      fclose(fp);
      fp = fopen("tmp.lammps.variable.lock","w");
      fprintf(fp,"%d\n",nextindex+1);
      fclose(fp);
      rename("tmp.lammps.variable.lock","tmp.lammps.variable");
      if (universe->uscreen)
	fprintf(universe->uscreen,
		"Increment via next: value %d on partition %d\n",
		nextindex+1,universe->iworld);
      if (universe->ulogfile)
	fprintf(universe->ulogfile,
		"Increment via next: value %d on partition %d\n",
		nextindex+1,universe->iworld);
    }
    MPI_Bcast(&nextindex,1,MPI_INT,0,world);

    for (int iarg = 0; iarg < narg; iarg++) {
      ivar = find(arg[iarg]);
      which[ivar] = nextindex;
      if (which[ivar] >= num[ivar]) {
	flag = 1;
	remove(ivar);
      }
    }
  }

  return flag;
}

/* ----------------------------------------------------------------------
   return ptr to the data text associated with a variable
   if INDEX or WORLD or UNIVERSE or STRING var, return ptr to stored string
   if LOOP or ULOOP var, write int to data[0] and return ptr to string
   if EQUAL var, evaluate variable and put result in str
   if ATOM var, return NULL
   return NULL if no variable or which is bad, caller must respond
------------------------------------------------------------------------- */

char *Variable::retrieve(char *name)
{
  int ivar = find(name);
  if (ivar == -1) return NULL;
  if (which[ivar] >= num[ivar]) return NULL;

  char *str;
  if (style[ivar] == INDEX || style[ivar] == WORLD || 
      style[ivar] == UNIVERSE || style[ivar] == STRING) {
    str = data[ivar][which[ivar]];
  } else if (style[ivar] == LOOP || style[ivar] == ULOOP) {
    char result[16];
    if (pad[ivar] == 0) sprintf(result,"%d",which[ivar]+1);
    else {
      char padstr[16];
      sprintf(padstr,"%%0%dd",pad[ivar]);
      sprintf(result,padstr,which[ivar]+1);
    }
    int n = strlen(result) + 1;
    delete [] data[ivar][0];
    data[ivar][0] = new char[n];
    strcpy(data[ivar][0],result);
    str = data[ivar][0];
  } else if (style[ivar] == EQUAL) {
    char result[32];
    double answer = evaluate(data[ivar][0],NULL);
    sprintf(result,"%.10g",answer);
    int n = strlen(result) + 1;
    if (data[ivar][1]) delete [] data[ivar][1];
    data[ivar][1] = new char[n];
    strcpy(data[ivar][1],result);
    str = data[ivar][1];
  } else if (style[ivar] == ATOM) return NULL;

  return str;
}

/* ----------------------------------------------------------------------
   return result of equal-style variable evaluation
------------------------------------------------------------------------- */

double Variable::compute_equal(int ivar)
{
  return evaluate(data[ivar][0],NULL);
}

/* ----------------------------------------------------------------------
   compute result of atom-style variable evaluation
   only computed for atoms in igroup, else result is 0.0
   answers are placed every stride locations into result
   if sumflag, add variable values to existing result
------------------------------------------------------------------------- */

void Variable::compute_atom(int ivar, int igroup,
			    double *result, int stride, int sumflag)
{
  Tree *tree;
  double tmp = evaluate(data[ivar][0],&tree);
  tmp = collapse_tree(tree);

  int groupbit = group->bitmask[igroup];
  int *mask = atom->mask;
  int nlocal = atom->nlocal;

  if (sumflag == 0) {
    int m = 0;
    for (int i = 0; i < nlocal; i++) {
      if (mask[i] && groupbit) result[m] = eval_tree(tree,i);
      else result[m] = 0.0;
      m += stride;
    }

  } else {
    int m = 0;
    for (int i = 0; i < nlocal; i++) {
      if (mask[i] && groupbit) result[m] += eval_tree(tree,i);
      m += stride;
    }
  }

  free_tree(tree);
}

/* ----------------------------------------------------------------------
   search for name in list of variables names
   return index or -1 if not found
------------------------------------------------------------------------- */
  
int Variable::find(char *name)
{
  for (int i = 0; i < nvar; i++)
    if (strcmp(name,names[i]) == 0) return i;
  return -1;
}

/* ----------------------------------------------------------------------
   return 1 if variable is EQUAL style, 0 if not
------------------------------------------------------------------------- */
  
int Variable::equalstyle(int ivar)
{
  if (style[ivar] == EQUAL) return 1;
  return 0;
}

/* ----------------------------------------------------------------------
   return 1 if variable is ATOM style, 0 if not
------------------------------------------------------------------------- */
  
int Variable::atomstyle(int ivar)
{
  if (style[ivar] == ATOM) return 1;
  return 0;
}

/* ----------------------------------------------------------------------
   remove Nth variable from list and compact list
------------------------------------------------------------------------- */
  
void Variable::remove(int n)
{
  delete [] names[n];
  if (style[n] == LOOP || style[n] == ULOOP) delete [] data[n][0];
  else for (int i = 0; i < num[n]; i++) delete [] data[n][i];
  delete [] data[n];

  for (int i = n+1; i < nvar; i++) {
    names[i-1] = names[i];
    style[i-1] = style[i];
    num[i-1] = num[i];
    which[i-1] = which[i];
    pad[i-1] = pad[i];
    data[i-1] = data[i];
  }
  nvar--;
}

/* ----------------------------------------------------------------------
  make space in arrays for new variable
------------------------------------------------------------------------- */

void Variable::extend()
{
  maxvar += VARDELTA;
  names = (char **)
    memory->srealloc(names,maxvar*sizeof(char *),"var:names");
  memory->grow(style,maxvar,"var:style");
  memory->grow(num,maxvar,"var:num");
  memory->grow(which,maxvar,"var:which");
  memory->grow(pad,maxvar,"var:pad");
  data = (char ***) 
    memory->srealloc(data,maxvar*sizeof(char **),"var:data");
}

/* ----------------------------------------------------------------------
   copy narg strings from **from to **to, and allocate space for them
------------------------------------------------------------------------- */
  
void Variable::copy(int narg, char **from, char **to)
{
  int n;
  for (int i = 0; i < narg; i++) {
    n = strlen(from[i]) + 1;
    to[i] = new char[n];
    strcpy(to[i],from[i]);
  }
}

/* ----------------------------------------------------------------------
   recursive evaluation of a string str
   str is an equal-style or atom-style formula containing one or more items:
     number = 0.0, -5.45, 2.8e-4, ...
     constant = PI
     thermo keyword = ke, vol, atoms, ...
     math operation = (),-x,x+y,x-y,x*y,x/y,x^y,
                      x==y,x!=y,x<y,x<=y,x>y,x>=y,x&&y,x||y,
                      sqrt(x),exp(x),ln(x),log(x),
		      sin(x),cos(x),tan(x),asin(x),atan2(y,x),...
     group function = count(group), mass(group), xcm(group,x), ...
     special function = sum(x),min(x), ...
     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]
   equal-style variables passes in tree = NULL:
     evaluate the formula, return result as a double
   atom-style variable passes in tree = non-NULL:
     parse the formula but do not evaluate it
     create a parse tree and return it
------------------------------------------------------------------------- */

double Variable::evaluate(char *str, Tree **tree)
{
  int op,opprevious;
  double value1,value2;
  char onechar;
  char *ptr;

  double argstack[MAXLEVEL];
  Tree *treestack[MAXLEVEL];
  int opstack[MAXLEVEL];
  int nargstack = 0;
  int ntreestack = 0;
  int nopstack = 0;

  int i = 0;
  int expect = ARG;

  while (1) {
    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");
      expect = OP;

      char *contents;
      i = find_matching_paren(str,i,contents);
      i++;

      // evaluate contents and push on stack

      if (tree) {
	Tree *newtree;
	double tmp = evaluate(contents,&newtree);
	treestack[ntreestack++] = newtree;
      } else argstack[nargstack++] = evaluate(contents,NULL);

      delete [] contents;

    // ----------------
    // number: push value onto stack
    // ----------------

    } else if (isdigit(onechar) || onechar == '.') {
      if (expect == OP) error->all("Invalid syntax in variable formula");
      expect = OP;

      // istop = end of number, including scientific notation

      int istart = i;
      while (isdigit(str[i]) || str[i] == '.') i++;
      if (str[i] == 'e' || str[i] == 'E') {
	i++;
	if (str[i] == '+' || str[i] == '-') i++;
	while (isdigit(str[i])) i++;
      }
      int istop = i - 1;

      int n = istop - istart + 1;
      char *number = new char[n+1];
      strncpy(number,&str[istart],n);
      number[n] = '\0';

      if (tree) {
        Tree *newtree = new Tree();
	newtree->type = VALUE;
	newtree->value = atof(number);
	newtree->left = newtree->middle = newtree->right = NULL;
	treestack[ntreestack++] = newtree;
      } else argstack[nargstack++] = atof(number);

      delete [] number;

    // ----------------
    // letter: c_ID, c_ID[], c_ID[][], f_ID, f_ID[], f_ID[][],
    //         v_name, v_name[], exp(), xcm(,), x, x[], PI, vol
    // ----------------

    } else if (isalpha(onechar)) {
      if (expect == OP) error->all("Invalid syntax in variable formula");
      expect = OP;

      // istop = end of word
      // word = all alphanumeric or underscore

      int istart = i;
      while (isalnum(str[i]) || str[i] == '_') i++;
      int istop = i-1;

      int n = istop - istart + 1;
      char *word = new char[n+1];
      strncpy(word,&str[istart],n);
      word[n] = '\0';

      // ----------------
      // compute
      // ----------------

      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]);

	int icompute = modify->find_compute(id);
	if (icompute < 0) error->all("Invalid compute ID in variable formula");
	Compute *compute = modify->compute[icompute];
	delete [] id;

	// parse zero or one or two trailing brackets
	// point i beyond last bracket
	// nbracket = # of bracket pairs
	// index1,index2 = int inside each bracket pair

	int nbracket,index1,index2;
	if (str[i] != '[') nbracket = 0;
	else {
	  nbracket = 1;
	  ptr = &str[i];
	  index1 = int_between_brackets(ptr);
	  i = ptr-str+1;
	  if (str[i] == '[') {
	    nbracket = 2;
	    ptr = &str[i];
	    index2 = int_between_brackets(ptr);
	    i = ptr-str+1;
	  }
	}

        // c_ID = scalar from global scalar

	if (nbracket == 0 && compute->scalar_flag) {

	  if (update->whichflag == 0) {
	    if (compute->invoked_scalar != update->ntimestep)
	      error->all("Compute used in variable between runs "
			 "is not current");
	  } else if (!(compute->invoked_flag & INVOKED_SCALAR)) {
	    compute->compute_scalar();
	    compute->invoked_flag |= INVOKED_SCALAR;
	  }

	  value1 = compute->scalar;
	  if (tree) {
	    Tree *newtree = new Tree();
	    newtree->type = VALUE;
	    newtree->value = value1;
	    newtree->left = newtree->middle = newtree->right = NULL;
	    treestack[ntreestack++] = newtree;
	  } else argstack[nargstack++] = value1;

        // c_ID[i] = scalar from global vector

	} else if (nbracket == 1 && compute->vector_flag) {

	  if (index1 > compute->size_vector)
	    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 "
			 "is not current");
	  } else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
	    compute->compute_vector();
	    compute->invoked_flag |= INVOKED_VECTOR;
	  }

	  value1 = compute->vector[index1-1];
	  if (tree) {
	    Tree *newtree = new Tree();
	    newtree->type = VALUE;
	    newtree->value = value1;
	    newtree->left = newtree->middle = newtree->right = NULL;
	    treestack[ntreestack++] = newtree;
	  } else argstack[nargstack++] = value1;

        // c_ID[i][j] = scalar from global array

	} else if (nbracket == 2 && compute->array_flag) {

	  if (index1 > compute->size_array_rows)
	    error->all("Variable formula compute array "
		       "is accessed out-of-range");
	  if (index2 > compute->size_array_cols)
	    error->all("Variable formula compute array "
		       "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->middle = 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 (update->whichflag == 0) {
	    if (compute->invoked_peratom != update->ntimestep)
	      error->all("Compute used in variable between runs "
			 "is not current");
	  } else if (!(compute->invoked_flag & INVOKED_PERATOM)) {
	    compute->compute_peratom();
	    compute->invoked_flag |= INVOKED_PERATOM;
	  }

	  peratom2global(1,NULL,compute->vector_atom,1,index1,
			 tree,treestack,ntreestack,argstack,nargstack);

        // c_ID[i][j] = scalar from per-atom array

	} else if (nbracket == 2 && compute->peratom_flag &&
		   compute->size_peratom_cols > 0) {

	  if (index2 > compute->size_peratom_cols)
	    error->all("Variable formula compute array "
		       "is accessed out-of-range");
	  if (update->whichflag == 0) {
	    if (compute->invoked_peratom != update->ntimestep)
	      error->all("Compute used in variable between runs "
			 "is not current");
	  } else if (!(compute->invoked_flag & INVOKED_PERATOM)) {
	    compute->compute_peratom();
	    compute->invoked_flag |= INVOKED_PERATOM;
	  }

	  peratom2global(1,NULL,&compute->array_atom[0][index2-1],
			 compute->size_peratom_cols,index1,
			 tree,treestack,ntreestack,argstack,nargstack);

        // c_ID = vector from per-atom vector

	} 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 (update->whichflag == 0) {
	    if (compute->invoked_peratom != update->ntimestep)
	      error->all("Compute used in variable between runs "
			 "is not current");
	  } else if (!(compute->invoked_flag & INVOKED_PERATOM)) {
	    compute->compute_peratom();
	    compute->invoked_flag |= INVOKED_PERATOM;
	  }

	  Tree *newtree = new Tree();
	  newtree->type = ATOMARRAY;
	  newtree->array = compute->vector_atom;
	  newtree->nstride = 1;
	  newtree->left = newtree->middle = newtree->right = NULL;
	  treestack[ntreestack++] = newtree;

        // c_ID[i] = vector from per-atom array

	} 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 (index1 > compute->size_peratom_cols)
	    error->all("Variable formula compute array "
		       "is accessed out-of-range");
	  if (update->whichflag == 0) {
	    if (compute->invoked_peratom != update->ntimestep)
	      error->all("Compute used in variable between runs "
			 "is not current");
	  } else if (!(compute->invoked_flag & INVOKED_PERATOM)) {
	    compute->compute_peratom();
	    compute->invoked_flag |= INVOKED_PERATOM;
	  }

	  Tree *newtree = new Tree();
	  newtree->type = ATOMARRAY;
	  newtree->array = &compute->array_atom[0][index1-1];
	  newtree->nstride = compute->size_peratom_cols;
	  newtree->left = newtree->middle = newtree->right = NULL;
	  treestack[ntreestack++] = newtree;

	} else error->all("Mismatched compute in variable formula");

      // ----------------
      // fix
      // ----------------

      } 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]);

	int ifix = modify->find_fix(id);
	if (ifix < 0) error->all("Invalid fix ID in variable formula");
	Fix *fix = modify->fix[ifix];
	delete [] id;

	// parse zero or one or two trailing brackets
	// point i beyond last bracket
	// nbracket = # of bracket pairs
	// index1,index2 = int inside each bracket pair

	int nbracket,index1,index2;
	if (str[i] != '[') nbracket = 0;
	else {
	  nbracket = 1;
	  ptr = &str[i];
	  index1 = int_between_brackets(ptr);
	  i = ptr-str+1;
	  if (str[i] == '[') {
	    nbracket = 2;
	    ptr = &str[i];
	    index2 = int_between_brackets(ptr);
	    i = ptr-str+1;
	  }
	}

        // 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();
	    newtree->type = VALUE;
	    newtree->value = value1;
	    newtree->left = newtree->middle = newtree->right = NULL;
	    treestack[ntreestack++] = newtree;
	  } else argstack[nargstack++] = value1;

        // f_ID[i] = scalar from global vector

	} else if (nbracket == 1 && fix->vector_flag) {

	  if (index1 > fix->size_vector)
	    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();
	    newtree->type = VALUE;
	    newtree->value = value1;
	    newtree->left = newtree->middle = newtree->right = NULL;
	    treestack[ntreestack++] = newtree;
	  } else argstack[nargstack++] = value1;

        // f_ID[i][j] = scalar from global array

	} else if (nbracket == 2 && fix->array_flag) {

	  if (index1 > fix->size_array_rows)
	    error->all("Variable formula fix array is accessed out-of-range");
	  if (index2 > fix->size_array_cols)
	    error->all("Variable formula fix array 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->middle = 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 (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,
			 tree,treestack,ntreestack,argstack,nargstack);

        // 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 array 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);

        // f_ID = vector from per-atom vector

	} 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 (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->middle = newtree->right = NULL;
	  treestack[ntreestack++] = newtree;

        // f_ID[i] = vector from per-atom array

	} 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 (index1 > fix->size_peratom_cols)
	    error->all("Variable formula fix array is accessed out-of-range");
	  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][index1-1];
	  newtree->nstride = fix->size_peratom_cols;
	  newtree->left = newtree->middle = newtree->right = NULL;
	  treestack[ntreestack++] = newtree;


	} else error->all("Mismatched fix in variable formula");

      // ----------------
      // variable
      // ----------------

      } else if (strncmp(word,"v_",2) == 0) {
	n = strlen(word) - 2 + 1;
	char *id = new char[n];
	strcpy(id,&word[2]);

	int ivar = find(id);
	if (ivar < 0) error->all("Invalid variable name in variable formula");

	// parse zero or one trailing brackets
	// point i beyond last bracket
	// nbracket = # of bracket pairs
	// index = int inside bracket

	int nbracket,index;
	if (str[i] != '[') nbracket = 0;
	else {
	  nbracket = 1;
	  ptr = &str[i];
	  index = int_between_brackets(ptr);
	  i = ptr-str+1;
	}

        // v_name = scalar from non atom-style global scalar

	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);
	    newtree->left = newtree->middle = newtree->right = NULL;
	    treestack[ntreestack++] = newtree;
	  } else argstack[nargstack++] = atof(var);

        // v_name = vector from atom-style per-atom vector

	} else if (nbracket == 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;

        // v_name[N] = scalar from atom-style per-atom vector
	// compute the per-atom variable in result
	// use peratom2global to extract single value from result

	} else if (nbracket && style[ivar] == ATOM) {

	  double *result;
	  memory->create(result,atom->nlocal,"variable:result");
	  compute_atom(ivar,0,result,1,0);
	  peratom2global(1,NULL,result,1,index,
			 tree,treestack,ntreestack,argstack,nargstack);
	  memory->destroy(result);

	} else error->all("Mismatched variable in variable formula");

	delete [] id;

      // ----------------
      // math/group/special function or atom value/vector or
      // constant or thermo keyword
      // ----------------

      } else {

	// ----------------
	// math or group or special function
	// ----------------

	if (str[i] == '(') {
	  char *contents;
	  i = find_matching_paren(str,i,contents);
	  i++;

	  if (math_function(word,contents,tree,
			    treestack,ntreestack,argstack,nargstack));
	  else if (group_function(word,contents,tree,
				  treestack,ntreestack,argstack,nargstack));
	  else if (special_function(word,contents,tree,
				    treestack,ntreestack,argstack,nargstack));
	  else error->all("Invalid math/group/special function "
			  "in variable formula");
	  delete [] contents;

	// ----------------
	// atom value
	// ----------------

	} else if (str[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;
 
	  peratom2global(0,word,NULL,0,id,
			 tree,treestack,ntreestack,argstack,nargstack);

	// ----------------
	// 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);

	// ----------------
	// constant
	// ----------------

	} else if (is_constant(word)) {
	  value1 = constant(word);
	  if (tree) {
	    Tree *newtree = new Tree();
	    newtree->type = VALUE;
	    newtree->value = value1;
	    newtree->left = newtree->middle = newtree->right = NULL;
	    treestack[ntreestack++] = newtree;
	  } else argstack[nargstack++] = value1;

	// ----------------
	// thermo keyword
	// ----------------

	} else {
	  if (domain->box_exist == 0)
	    error->all("Variable evaluation before simulation box is defined");
 
	  int flag = output->thermo->evaluate_keyword(word,&value1);
	  if (flag) error->all("Invalid thermo keyword in variable formula");
	  if (tree) {
	    Tree *newtree = new Tree();
	    newtree->type = VALUE;
	    newtree->value = value1;
	    newtree->left = newtree->middle = newtree->right = NULL;
	    treestack[ntreestack++] = newtree;
	  } else argstack[nargstack++] = value1;
	}
      }

      delete [] word;

    // ----------------
    // math operator, including end-of-string
    // ----------------

    } else if (strchr("+-*/^<>=!&|\0",onechar)) {
      if (onechar == '+') op = ADD;
      else if (onechar == '-') op = SUBTRACT;
      else if (onechar == '*') op = MULTIPLY;
      else if (onechar == '/') op = DIVIDE;
      else if (onechar == '^') op = CARAT;
      else if (onechar == '=') {
	if (str[i+1] != '=') error->all("Invalid syntax in variable formula");
	op = EQ;
	i++;
      } else if (onechar == '!') {
	if (str[i+1] == '=') {
	  op = NE;
	  i++;
	} else op = NOT;
      } else if (onechar == '<') {
	if (str[i+1] != '=') op = LT;
	else {
	  op = LE;
	  i++;
	}
      } else if (onechar == '>') {
	if (str[i+1] != '=') op = GT;
	else {
	  op = GE;
	  i++;
	}
      } else if (onechar == '&') {
	if (str[i+1] != '&') error->all("Invalid syntax in variable formula");
	op = AND;
	i++;
      } else if (onechar == '|') {
	if (str[i+1] != '|') error->all("Invalid syntax in variable formula");
	op = OR;
	i++;
      } else op = DONE;

      i++;

      if (op == SUBTRACT && expect == ARG) {
	opstack[nopstack++] = UNARY;
	continue;
      }
      if (op == NOT && expect == ARG) {
	opstack[nopstack++] = op;
	continue;
      }

      if (expect == ARG) error->all("Invalid syntax in variable formula");
      expect = ARG;

      // evaluate stack as deep as possible while respecting precedence
      // before pushing current op onto stack

      while (nopstack && precedence[opstack[nopstack-1]] >= precedence[op]) {
	opprevious = opstack[--nopstack];

	if (tree) {
	  Tree *newtree = new Tree();
	  newtree->type = opprevious;
	  if (opprevious == UNARY) {
	    newtree->left = treestack[--ntreestack];
	    newtree->middle = newtree->right = NULL;
	  } else {
	    newtree->right = treestack[--ntreestack];
	    newtree->middle = NULL;
	    newtree->left = treestack[--ntreestack];
	  }
	  treestack[ntreestack++] = newtree;

	} else {
	  value2 = argstack[--nargstack];
	  if (opprevious != UNARY && opprevious != NOT)
	    value1 = argstack[--nargstack];

	  if (opprevious == ADD)
	    argstack[nargstack++] = value1 + value2;
	  else if (opprevious == SUBTRACT)
	    argstack[nargstack++] = value1 - value2;
	  else if (opprevious == MULTIPLY)
	    argstack[nargstack++] = value1 * value2;
	  else if (opprevious == DIVIDE) {
	    if (value2 == 0.0) error->all("Divide by 0 in variable formula");
	    argstack[nargstack++] = value1 / value2;
	  } else if (opprevious == CARAT) {
	    if (value2 == 0.0) error->all("Power by 0 in variable formula");
	    argstack[nargstack++] = pow(value1,value2);
	  } else if (opprevious == UNARY) {
	    argstack[nargstack++] = -value2;
	  } else if (opprevious == NOT) {
	    if (value2 == 0.0) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  } else if (opprevious == EQ) {
	    if (value1 == value2) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  } else if (opprevious == NE) {
	    if (value1 != value2) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  } else if (opprevious == LT) {
	    if (value1 < value2) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  } else if (opprevious == LE) {
	    if (value1 <= value2) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  } else if (opprevious == GT) {
	    if (value1 > value2) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  } else if (opprevious == GE) {
	    if (value1 >= value2) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  } else if (opprevious == AND) {
	    if (value1 != 0.0 && value2 != 0.0) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  } else if (opprevious == OR) {
	    if (value1 != 0.0 || value2 != 0.0) argstack[nargstack++] = 1.0;
	    else argstack[nargstack++] = 0.0;
	  }
	}
      }

      // if end-of-string, break out of entire formula evaluation loop

      if (op == DONE) break;

      // push current operation onto stack

      opstack[nopstack++] = op;

    } else error->all("Invalid syntax in variable formula");
  }

  if (nopstack) error->all("Invalid syntax in variable formula");

  // for atom-style variable, return remaining tree
  // for equal-style variable, return remaining arg

  if (tree) {
    if (ntreestack != 1) error->all("Invalid syntax in variable formula");
    *tree = treestack[0];
    return 0.0;
  } else {
    if (nargstack != 1) error->all("Invalid syntax in variable formula");
    return argstack[0];
  }
}

/* ----------------------------------------------------------------------
   one-time collapse of an atom-style variable parse tree
   tree was created by one-time parsing of formula string via evaulate()
   only keep tree nodes that depend on ATOMARRAY, TYPEARRAY, INTARRAY
   remainder is converted to single VALUE
   this enables optimal eval_tree loop over atoms
   customize by adding a function:
     sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan(),
     atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
     ramp(x,y),stagger(x,y),logfreq(x,y,z),
     vdisplace(x,y),swiggle(x,y,z),cwiggle(x,y,z),
     gmask(x),rmask(x),grmask(x,y)
---------------------------------------------------------------------- */

double Variable::collapse_tree(Tree *tree)
{
  double arg1,arg2,arg3;

  if (tree->type == VALUE) return tree->value;
  if (tree->type == ATOMARRAY) return 0.0;
  if (tree->type == TYPEARRAY) return 0.0;
  if (tree->type == INTARRAY) return 0.0;

  if (tree->type == ADD) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    tree->value = arg1 + arg2;
    return tree->value;
  }

  if (tree->type == SUBTRACT) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    tree->value = arg1 - arg2;
    return tree->value;
  }

  if (tree->type == MULTIPLY) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    tree->value = arg1 * arg2;
    return tree->value;
  }

  if (tree->type == DIVIDE) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg2 == 0.0) error->one("Divide by 0 in variable formula");
    tree->value = arg1 / arg2;
    return tree->value;
  }

  if (tree->type == CARAT) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg2 == 0.0) error->one("Power by 0 in variable formula");
    tree->value = pow(arg1,arg2);
    return tree->value;
  }

  if (tree->type == UNARY) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    tree->value = -arg1;
    return tree->value;
  }

  if (tree->type == NOT) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 == 0.0) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == EQ) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 == arg2) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == NE) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 != arg2) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == LT) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 < arg2) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == LE) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 <= arg2) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == GT) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 > arg2) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == GE) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 >= arg2) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == AND) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 != 0.0 && arg2 != 0.0) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == OR) {
    arg1 = collapse_tree(tree->left);
    arg2 = collapse_tree(tree->right);
    if (tree->left->type != VALUE || tree->right->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 != 0.0 || arg2 != 0.0) tree->value = 1.0;
    else tree->value = 0.0;
    return tree->value;
  }

  if (tree->type == SQRT) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 < 0.0) error->one("Sqrt of negative value in variable formula");
    tree->value = sqrt(arg1);
    return tree->value;
  }

  if (tree->type == EXP) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    tree->value = exp(arg1);
    return tree->value;
  }

  if (tree->type == LN) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 <= 0.0) 
      error->one("Log of zero/negative value in variable formula");
    tree->value = log(arg1);
    return tree->value;
  }

  if (tree->type == LOG) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 <= 0.0) 
      error->one("Log of zero/negative value in variable formula");
    tree->value = log10(arg1);
    return tree->value;
  }

  if (tree->type == SIN) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    tree->value = sin(arg1);
    return tree->value;
  }

  if (tree->type == COS) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    tree->value = cos(arg1);
    return tree->value;
  }

  if (tree->type == TAN) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    tree->value = tan(arg1);
    return tree->value;
  }

  if (tree->type == ASIN) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
    tree->type = VALUE;
    if (arg1 < -1.0 || arg1 > 1.0)
      error->one("Arcsin of invalid value in variable formula");
    tree->value = asin(arg1);
    return tree->value;
  }

  if (tree->type == ACOS) {
    arg1 = collapse_tree(tree->left);
    if (tree->left->type != VALUE) return 0.0;
…