git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@322 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/KSPACE/ewald.h

@@ -50,7 +50,8 @@ class Ewald : public KSpace {
   void slabcorr(int);
 };
 
-#endif
 }
 
+#endif
+
 

src/compute.h

@@ -23,15 +23,15 @@ class Compute : protected Pointers {
   char *id,*style;
   int igroup,groupbit;
 
-  double scalar;            // computed scalar
+  double scalar;            // computed global scalar
-  double *vector;           // computed vector
+  double *vector;           // computed global vector
   double *scalar_atom;      // computed per-atom scalar
   double **vector_atom;     // computed per-atom vector
 
   int scalar_flag;          // 0/1 if compute_scalar() function exists
   int vector_flag;          // 0/1 if compute_vector() function exists
   int peratom_flag;         // 0/1 if compute_peratom() function exists
-  int size_vector;          // N = size of vector
+  int size_vector;          // N = size of global vector
   int size_peratom;         // 0 = just scalar_atom, N = size of vector_atom
 
   int extensive;      // 0/1 if scalar,vector are intensive/extensive values

src/fix.cpp

@@ -37,6 +37,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   restart_peratom = 0;
   force_reneighbor = 0;
   thermo_energy = 0;
+  pressure_every = 0;
 
   comm_forward = comm_reverse = 0;
   neigh_half_once = neigh_half_every = 0;
@@ -55,6 +56,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   POST_FORCE_RESPA = 256;
   FINAL_INTEGRATE_RESPA = 512;
   MIN_POST_FORCE = 1024;
+  MIN_ENERGY = 2048;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix.h

@@ -29,6 +29,7 @@ class Fix : protected Pointers {
   int next_reneighbor;           // next timestep to force a reneighboring
   int thermo_energy;             // 1 if ThEng enabled via fix_modify, 0 if not
   int nevery;                    // how often to call an end_of_step fix
+  int pressure_every;            // how often fix needs virial computed
 
   int comm_forward;              // size of forward communication (0 if none)
   int comm_reverse;              // size of reverse communication (0 if none)
@@ -42,7 +43,7 @@ class Fix : protected Pointers {
   int INITIAL_INTEGRATE,PRE_EXCHANGE,PRE_NEIGHBOR;    // mask settings
   int POST_FORCE,FINAL_INTEGRATE,END_OF_STEP,THERMO_ENERGY;
   int INITIAL_INTEGRATE_RESPA,POST_FORCE_RESPA,FINAL_INTEGRATE_RESPA;
-  int MIN_POST_FORCE;
+  int MIN_POST_FORCE,MIN_ENERGY;
 
   Fix(class LAMMPS *, int, char **);
   virtual ~Fix();
@@ -76,6 +77,9 @@ class Fix : protected Pointers {
   virtual void final_integrate_respa(int) {}
 
   virtual void min_post_force(int) {}
+  virtual double min_energy(double *, double *) {return 0.0;}
+  virtual int min_dof() {return 0;}
+  virtual void min_xinitial(double *) {}
 
   virtual int pack_comm(int, int *, double *, int *) {return 0;}
   virtual void unpack_comm(int, int, double *) {}

src/fix_nph.cpp

@@ -44,6 +44,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
   if (narg < 4) error->all("Illegal fix nph command");
 
   restart_global = 1;
+  pressure_every = 1;
 
   double p_period[3];
   if (strcmp(arg[3],"xyz") == 0) {

src/fix_npt.cpp

@@ -44,6 +44,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
   if (narg < 7) error->all("Illegal fix npt command");
 
   restart_global = 1;
+  pressure_every = 1;
 
   t_start = atof(arg[3]);
   t_stop = atof(arg[4]);

src/modify.cpp

@@ -47,6 +47,7 @@ using namespace LAMMPS_NS;
 #define POST_FORCE_RESPA        256
 #define FINAL_INTEGRATE_RESPA   512
 #define MIN_POST_FORCE         1024
+#define MIN_ENERGY             2048
 
 /* ---------------------------------------------------------------------- */
 
@@ -57,7 +58,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
   n_pre_exchange = n_pre_neighbor = 0;
   n_post_force = n_final_integrate = n_end_of_step = n_thermo_energy = 0;
   n_initial_integrate_respa = n_post_force_respa = n_final_integrate_respa = 0;
-  n_min_post_force = 0;
+  n_min_post_force = n_min_energy = 0;
 
   fix = NULL;
   fmask = NULL;
@@ -67,7 +68,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
   list_thermo_energy = NULL;
   list_initial_integrate_respa = list_post_force_respa = NULL;
   list_final_integrate_respa = NULL;
-  list_min_post_force = NULL;
+  list_min_post_force = list_min_energy = NULL;
 
   end_of_step_every = NULL;
 
@@ -103,6 +104,7 @@ Modify::~Modify()
   delete [] list_post_force_respa;
   delete [] list_final_integrate_respa;
   delete [] list_min_post_force;
+  delete [] list_min_energy;
 
   delete [] end_of_step_every;
 
@@ -149,6 +151,7 @@ void Modify::init()
 	    n_final_integrate_respa,list_final_integrate_respa);
 
   list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
+  list_init(MIN_ENERGY,n_min_energy,list_min_energy);
 
   // init each compute
 
@@ -283,6 +286,50 @@ void Modify::min_post_force(int vflag)
     fix[list_min_post_force[i]]->min_post_force(vflag);
 }
 
+/* ----------------------------------------------------------------------
+   minimizer energy, force evaluation only for relevant fixes
+   return energy and forces on extra degrees of freedom
+------------------------------------------------------------------------- */
+
+double Modify::min_energy(double *xextra, double *fextra)
+{
+  int ifix,index;
+
+  index = 0;
+  double energy_extra = 0.0;
+  for (int i = 0; i < n_min_energy; i++) {
+    ifix = list_min_energy[i];
+    energy_extra += fix[ifix]->min_energy(&xextra[index],&fextra[index]);
+    index += fix[ifix]->min_dof();
+  }
+  return energy_extra;
+}
+
+/* ----------------------------------------------------------------------
+   minimizer extra degrees of freedom from relevant fixes
+------------------------------------------------------------------------- */
+
+int Modify::min_dof()
+{
+  int ndof = 0;
+  for (int i = 0; i < n_min_energy; i++)
+    ndof += fix[list_min_energy[i]]->min_dof();
+  return ndof;
+}
+
+/* ----------------------------------------------------------------------
+   minimizer initial xextra values only from relevant fixes
+------------------------------------------------------------------------- */
+
+void Modify::min_xinitial(double *xextra)
+{
+  int index = 0;
+  for (int i = 0; i < n_min_energy; i++) {
+    fix[list_min_energy[i]]->min_xinitial(&xextra[index]);
+    index += fix[list_min_energy[i]]->min_dof();
+  }
+}
+
 /* ----------------------------------------------------------------------
    add a new fix or replace one with same ID
 ------------------------------------------------------------------------- */

src/modify.h

@@ -25,7 +25,7 @@ class Modify : protected Pointers {
   int n_initial_integrate,n_pre_decide,n_pre_exchange,n_pre_neighbor;
   int n_post_force,n_final_integrate,n_end_of_step,n_thermo_energy;
   int n_initial_integrate_respa,n_post_force_respa,n_final_integrate_respa;
-  int n_min_post_force;
+  int n_min_post_force,n_min_energy;
   int nfix_restart_peratom;
 
   class Fix **fix;           // list of fixes
@@ -52,6 +52,9 @@ class Modify : protected Pointers {
   void final_integrate_respa(int);
 
   void min_post_force(int);
+  double min_energy(double *, double *);
+  int min_dof();
+  void min_xinitial(double *);
 
   void add_fix(int, char **);
   void modify_fix(int, char **);
@@ -77,7 +80,7 @@ class Modify : protected Pointers {
   int *list_thermo_energy;
   int *list_initial_integrate_respa,*list_post_force_respa;
   int *list_final_integrate_respa;
-  int *list_min_post_force;
+  int *list_min_post_force,*list_min_energy;
 
   int *end_of_step_every;
 

src/respa.cpp

@@ -40,6 +40,9 @@
 
 using namespace LAMMPS_NS;
 
+#define MIN(A,B) ((A) < (B)) ? (A) : (B)
+#define MAX(A,B) ((A) > (B)) ? (A) : (B)
+
 /* ---------------------------------------------------------------------- */
 
 Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
@@ -273,25 +276,30 @@ void Respa::init()
     if (force->pair && force->pair->respa_enable == 0)
       error->all("Pair style does not support rRESPA inner/middle/outer");
 
-  // set flags for how virial should be computed when needed
+  // setup virial computations for timestepping
-  // pressure_flag is 1 if NPT,NPH
+  // virial_style = 1 (explicit) since never computed implicitly like Verlet
-  // virial_every is how virial should be computed every timestep
+  // virial_every = 1 if computed every timestep (NPT,NPH)
-  //   0 = not computed, 1 = computed explicity by pair
+  // fix arrays store info on fixes that need virial computed occasionally
-  // virial_thermo is how virial should be computed on thermo timesteps
-  //   1 = computed explicity by pair
-  // unlike Verlet, virial is never computed implicitly
 
-  int pressure_flag = 0;
+  virial_style = 1;
-  for (int i = 0; i < modify->nfix; i++) {
+
-    if (strcmp(modify->fix[i]->style,"npt") == 0) pressure_flag = 1;
+  virial_every = 0;
-    if (strcmp(modify->fix[i]->style,"nph") == 0) pressure_flag = 1;
+  nfix_virial = 0;
+  for (int i = 0; i < modify->nfix; i++)
+    if (modify->fix[i]->pressure_every == 1) virial_every = 1;
+    else if (modify->fix[i]->pressure_every > 1) nfix_virial++;
+
+  if (nfix_virial) {
+    delete [] fix_virial_every;
+    delete [] next_fix_virial;
+    fix_virial_every = new int[nfix_virial];
+    next_fix_virial = new int[nfix_virial];
+    nfix_virial = 0;
+    for (int i = 0; i < modify->nfix; i++)
+      if (modify->fix[i]->pressure_every > 1)
+	fix_virial_every[nfix_virial++] = modify->fix[i]->pressure_every;
   }
 
-  if (pressure_flag) virial_every = 1;
-  else virial_every = 0;
-
-  virial_thermo = 1;
-
   // step[] = timestep for each level
 
   step[nlevels-1] = update->dt;
@@ -338,8 +346,8 @@ void Respa::setup()
 
   // compute all forces
 
-  int eflag = 1;
+  eflag = 1;
-  int vflag = virial_thermo;
+  vflag = virial_style;
 
   for (int ilevel = 0; ilevel < nlevels; ilevel++) {
     force_clear(newton[ilevel]);
@@ -370,6 +378,21 @@ void Respa::setup()
   modify->setup();
   sum_flevel_f();
   output->setup(1);
+
+  // setup virial computations for timestepping
+
+  int ntimestep = update->ntimestep;
+  next_virial = 0;
+  if (virial_every) next_virial = ntimestep + 1;
+  else {
+    for (int ivirial = 0; ivirial < nfix_virial; ivirial++) {
+      next_fix_virial[ivirial] = 
+	(ntimestep/fix_virial_every[ivirial])*fix_virial_every[ivirial] + 
+	fix_virial_every[ivirial];
+      if (ivirial) next_virial = MIN(next_virial,next_fix_virial[ivirial]);
+      else next_virial = next_fix_virial[0];
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -378,22 +401,28 @@ void Respa::setup()
 
 void Respa::iterate(int n)
 {
+  int ntimestep;
+
   for (int i = 0; i < n; i++) {
 
-    update->ntimestep++;
+    ntimestep = ++update->ntimestep;
 
-    eflag = 0;
+    // eflag/vflag = 0/1 for energy/virial computation
-    vflag = virial_every;
+
-    if (output->next_thermo == update->ntimestep) {
+    if (ntimestep == output->next_thermo) eflag = 1;
-      eflag = 1;
+    else eflag = 0;
-      vflag = virial_thermo;
+
-    }
+    if (ntimestep == output->next_thermo || ntimestep == next_virial) {
+      vflag = virial_style;
+      if (virial_every) next_virial++;
+      else next_virial = fix_virial(ntimestep);
+    } else vflag = 0;
 
     recurse(nlevels-1);
 
     if (modify->n_end_of_step) modify->end_of_step();
 
-    if (output->next == update->ntimestep) {
+    if (ntimestep == output->next) {
       timer->stamp();
       sum_flevel_f();
       output->write(update->ntimestep);
@@ -582,3 +611,20 @@ void Respa::sum_flevel_f()
     }
   }
 }
+
+/* ----------------------------------------------------------------------
+   return next timestep virial should be computed
+   based on one or more fixes that need virial computed periodically
+------------------------------------------------------------------------- */
+
+int Respa::fix_virial(int ntimestep)
+{
+  int next;
+  for (int ivirial = 0; ivirial < nfix_virial; ivirial++) {
+    if (ntimestep == next_fix_virial[ivirial])
+      next_fix_virial[ivirial] += fix_virial_every[ivirial];
+    if (ivirial) next = MIN(next,next_fix_virial[ivirial]);
+    else next = next_fix_virial[0];
+  }
+  return next;
+}

src/respa.h

@@ -44,16 +44,21 @@ class Respa : public Integrate {
   void copy_flevel_f(int);
 
  private:
-  int virial_every;       // what vflag should be on every timestep (0,1)
+  int eflag,vflag;                  // flags for energy/virial computation
-  int virial_thermo;      // what vflag should be on thermo steps (1)
+  int virial_style;                 // compute virial explicitly (not implicit)
+  int virial_every;                 // 1 if virial computed every step
+  int next_virial;                  // next timestep to compute virial
+  int nfix_virial;                  // # of fixes that need virial occasionally
+  int *fix_virial_every;            // frequency they require it
+  int *next_fix_virial;             // next timestep they need it
 
   int *newton;                      // newton flag at each level
-  int eflag,vflag;              // flags for energy/virial computation
   class FixRespa *fix_respa;        // Fix to store the force level array
 
   void recurse(int);
   void force_clear(int);
   void sum_flevel_f();
+  int fix_virial(int);
 };
 
 }

src/style.h

@@ -84,6 +84,7 @@ CommandStyle(write_restart,WriteRestart)
 #include "compute_temp_partial.h"
 #include "compute_temp_ramp.h"
 #include "compute_temp_region.h"
+#include "compute_variable_atom.h"
 #endif
 
 #ifdef ComputeClass
@@ -99,6 +100,7 @@ ComputeStyle(temp,ComputeTemp)
 ComputeStyle(temp/partial,ComputeTempPartial)
 ComputeStyle(temp/ramp,ComputeTempRamp)
 ComputeStyle(temp/region,ComputeTempRegion)
+ComputeStyle(variable/atom,ComputeVariableAtom)
 #endif
 
 #ifdef DihedralInclude
@@ -124,6 +126,9 @@ DumpStyle(xyz,DumpXYZ)
 #ifdef FixInclude
 #include "fix_add_force.h"
 #include "fix_ave_force.h"
+#include "fix_ave_spatial.h"
+#include "fix_ave_time.h"
+  //#include "fix_box_relax.h"
 #include "fix_com.h"
 #include "fix_drag.h"
 #include "fix_deposit.h"
@@ -168,6 +173,9 @@ DumpStyle(xyz,DumpXYZ)
 #ifdef FixClass
 FixStyle(addforce,FixAddForce)
 FixStyle(aveforce,FixAveForce)
+FixStyle(ave/spatial,FixAveSpatial)
+FixStyle(ave/time,FixAveTime)
+  //FixStyle(box/relax,FixBoxRelax)
 FixStyle(com,FixCOM)
 FixStyle(drag,FixDrag)
 FixStyle(deposit,FixDeposit)

src/variable.cpp

@@ -33,7 +33,8 @@ using namespace LAMMPS_NS;
 
 #define VARDELTA 4
 
-enum{INDEX,LOOP,EQUAL,WORLD,UNIVERSE,ULOOP};
+enum{INDEX,LOOP,EQUAL,WORLD,UNIVERSE,ULOOP,ATOM};
+enum{VALUE,ATOMARRAY,TYPEARRAY,ADD,SUB,MULT,DIV,NEG,POW,EXP,LN,SQRT};
 
 /* ---------------------------------------------------------------------- */
 
@@ -73,9 +74,9 @@ void Variable::set(int narg, char **arg)
 {
   if (narg < 3) error->all("Illegal variable command");
 
-  // if var already exists, just skip (except EQUAL vars)
+  // if var already exists, just skip
 
-  if (find(arg[0]) >= 0 && strcmp(arg[1],"equal") != 0) return;
+  if (find(arg[0]) >= 0) return;
       
   // make space for new variable
 
@@ -90,6 +91,11 @@ void Variable::set(int narg, char **arg)
       memory->srealloc(data,maxvar*sizeof(char **),"var:data");
   }
 
+  // set name of variable
+
+  names[nvar] = new char[strlen(arg[0]) + 1];
+  strcpy(names[nvar],arg[0]);
+
   // INDEX
   // num = listed args, index = 1st value, data = copied args
 
@@ -111,16 +117,10 @@ void Variable::set(int narg, char **arg)
     for (int i = 0; i < num[nvar]; i++) data[nvar][i] = NULL;
     
   // EQUAL
-  // remove pre-existing var if also style EQUAL (allows it to be reset)
   // num = 2, index = 1st value
   // data = 2 values, 1st is string to eval, 2nd is filled on retrieval
 
   } else if (strcmp(arg[1],"equal") == 0) {
-    if (find(arg[0]) >= 0) {
-      if (style[find(arg[0])] != EQUAL)
-	error->all("Cannot redefine variable as a different style");
-      remove(find(arg[0]));
-    }
     style[nvar] = EQUAL;
     num[nvar] = 2;
     index[nvar] = 0;
@@ -187,17 +187,24 @@ void Variable::set(int narg, char **arg)
 		arg[0],index[nvar]+1,universe->iworld);
     }
 
+  // ATOM
+  // num = 1, index = 1st value
+  // data = 1 value, string to eval
+
+  } else if (strcmp(arg[1],"atom") == 0) {
+    style[nvar] = ATOM;
+    num[nvar] = 1;
+    index[nvar] = 0;
+    data[nvar] = new char*[num[nvar]];
+    copy(1,&arg[2],data[nvar]);
+    
   } else error->all("Illegal variable command");
 
-  // set variable name (after possible EQUAL remove)
-
-  names[nvar] = new char[strlen(arg[0]) + 1];
-  strcpy(names[nvar],arg[0]);
   nvar++;
 }
 
 /* ----------------------------------------------------------------------
-   single-value INDEX variable created by command-line argument
+   single-value EQUAL variable created by command-line argument
 ------------------------------------------------------------------------- */
 
 void Variable::set(char *name, char *value)
@@ -205,25 +212,12 @@ void Variable::set(char *name, char *value)
   int ivar = find(name);
   if (ivar >= 0) error->all("Command-line variable already exists");
 
-  if (nvar == maxvar) {
+  char **newarg = new char*[3];
-    maxvar += VARDELTA;
+  newarg[0] = name;
-    names = (char **)
+  newarg[1] = "equal";
-      memory->srealloc(names,maxvar*sizeof(char *),"var:names");
+  newarg[2] = value;
-    style = (int *) memory->srealloc(style,maxvar*sizeof(int),"var:style");
+  set(3,newarg);
-    num = (int *) memory->srealloc(num,maxvar*sizeof(int),"var:num");
+  delete [] newarg;
-    index = (int *) memory->srealloc(index,maxvar*sizeof(int),"var:index");
-    data = (char ***) 
-      memory->srealloc(data,maxvar*sizeof(char **),"var:data");
-  }
-
-  names[nvar] = new char[strlen(name) + 1];
-  strcpy(names[nvar],name);
-  style[nvar] = INDEX;
-  num[nvar] = 1;
-  index[nvar] = 0;
-  data[nvar] = new char*[num[nvar]];
-  copy(1,&value,data[nvar]);
-  nvar++;
 }
 
 /* ----------------------------------------------------------------------
@@ -247,10 +241,10 @@ int Variable::next(int narg, char **arg)
       error->all("All variables in next command must be same style");
   }
 
-  // check for invalid styles EQUAL or WORLD
+  // invalid styles EQUAL or WORLD or ATOM
 
   int istyle = style[find(arg[0])];
-  if (istyle == EQUAL || istyle == WORLD)
+  if (istyle == EQUAL || istyle == WORLD || istyle == ATOM)
     error->all("Invalid variable style with next command");
 
   // increment all variables in list
@@ -336,14 +330,16 @@ char *Variable::retrieve(char *name)
     delete [] value;
     str = data[ivar][0];
   } else if (style[ivar] == EQUAL) {
-    char *value = evaluate(data[ivar][0]);
+    char result[32];
-    int n = strlen(value) + 1;
+    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],value);
+    strcpy(data[ivar][1],result);
-    delete [] value;
     str = data[ivar][1];
-  }
+  } else if (style[ivar] == ATOM) return NULL;
+
   return str;
 }
 
@@ -404,34 +400,34 @@ void Variable::copy(int narg, char **from, char **to)
      group function = mass(group), xcm(group,x), ...
      atom vector = x[123], y[3], vx[34], ...
      compute vector = c_mytemp[0], c_thermo_press[3], ...
+   numbers start with a digit or "." or "-" (no parens or brackets)
+   keywords start with a lowercase letter (no parens or brackets)
    functions contain ()
      can have 1 or 2 args, each of which can be a "string"
    vectors contain []
      single arg must be integer
-     for atom vectors, it is global ID of atom
+     for atom vectors, arg is global ID of atom
      for compute vectors, 0 is the scalar value, 1-N are vector values
-   keywords start with a lowercase letter (no parens or brackets)
+   see lists of valid functions & vectors below
-   numbers start with a digit or "." or "-" (no parens or brackets)
+   return answer = value of string
-   see lists of valid functions, vectors, keywords below
-   when string is evaluated, put result in a newly allocated string
-   return the address of result string (will be freed by caller)
 ------------------------------------------------------------------------- */
 
-char *Variable::evaluate(char *str)
+double Variable::evaluate(char *str, Tree *tree)
 {
-  // allocate a new string for the eventual result
+  double answer = 0.0;
-
+  if (tree) {
-  char *result = new char[32];
+    tree->type = VALUE;
-  double answer;
+    tree->left = tree->right = NULL;
+  }
 
   // string is a "function" since contains ()
   // grab one or two args, separated by ","
   // evaulate args recursively
-  // then evaluate math or group function
+  // if tree = NULL, evaluate math or group function
+  // else store as leaf in tree
 
   if (strchr(str,'(')) {
-    if (str[strlen(str)-1] != ')')
+    if (str[strlen(str)-1] != ')') error->all("Cannot evaluate variable");
-      error->all("Cannot evaluate variable equal command");
 
     char *ptr = strchr(str,'(');
     int n = ptr - str;
@@ -442,8 +438,7 @@ char *Variable::evaluate(char *str)
     char *comma = ++ptr;
     int level = 0;
     while (1) {
-      if (*comma == '\0')
+      if (*comma == '\0') error->all("Cannot evaluate variable");
-	error->all("Cannot evaluate variable equal command");
       else if (*comma == ',' && level == 0) break;
       else if (*comma == ')' && level == 0) break;
       else if (*comma == '(') level++;
@@ -467,106 +462,132 @@ char *Variable::evaluate(char *str)
     } else arg2 = NULL;
     
     double value1,value2;
-    char *strarg1 = NULL;
-    char *strarg2 = NULL;
 
     // customize by adding math function to this list and to if statement
     //   add(x,y),sub(x,y),mult(x,y),div(x,y),neg(x),
     //   pow(x,y),exp(x),ln(x),sqrt(x)
 
     if (strcmp(func,"add") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      strarg2 = evaluate(arg2);
+	tree->type = ADD;
-      value1 = atof(strarg1);
+	tree->left = new Tree();
-      value2 = atof(strarg2);
+	tree->right = new Tree();
-      answer = value1 + value2;
+	value1 = evaluate(arg1,tree->left);
+	value2 = evaluate(arg2,tree->right);
+      } else answer = evaluate(arg1,NULL) + evaluate(arg2,NULL);
 
     } else if (strcmp(func,"sub") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      strarg2 = evaluate(arg2);
+	tree->type = SUB;
-      value1 = atof(strarg1);
+	tree->left = new Tree();
-      value2 = atof(strarg2);
+	tree->right = new Tree();
-      answer = value1 - value2;
+	value1 = evaluate(arg1,tree->left);
+	value2 = evaluate(arg2,tree->right);
+      } else answer = evaluate(arg1,NULL) - evaluate(arg2,NULL);
 
     } else if (strcmp(func,"mult") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      strarg2 = evaluate(arg2);
+	tree->type = MULT;
-      value1 = atof(strarg1);
+	tree->left = new Tree();
-      value2 = atof(strarg2);
+	tree->right = new Tree();
-      answer = value1 * value2;
+	value1 = evaluate(arg1,tree->left);
+	value2 = evaluate(arg2,tree->right);
+      } else answer = evaluate(arg1,NULL) * evaluate(arg2,NULL);
 
     } else if (strcmp(func,"div") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      strarg2 = evaluate(arg2);
+	tree->type = DIV;
-      value1 = atof(strarg1);
+	tree->left = new Tree();
-      value2 = atof(strarg2);
+	tree->right = new Tree();
-      if (value2 == 0.0)
+	value1 = evaluate(arg1,tree->left);
-	error->all("Cannot evaluate variable equal command");
+	value2 = evaluate(arg2,tree->right);
+      } else {
+	value1 = evaluate(arg1,NULL);
+	value2 = evaluate(arg2,NULL);
+	if (value2 == 0.0) error->all("Cannot evaluate variable");
 	answer = value1 / value2;
+      }
 
     } else if (strcmp(func,"neg") == 0) {
-      if (arg2) error->all("Cannot evaluate variable equal command");
+      if (arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      value1 = atof(strarg1);
+	tree->type = NEG;
-      answer = -value1;
+	tree->left = new Tree();
+	value1 = evaluate(arg1,tree->left);
+      } else answer = -evaluate(arg1,NULL);
 
     } else if (strcmp(func,"pow") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      strarg2 = evaluate(arg2);
+	tree->type = POW;
-      value1 = atof(strarg1);
+	tree->left = new Tree();
-      value2 = atof(strarg2);
+	tree->right = new Tree();
-      if (value2 == 0.0)
+	value1 = evaluate(arg1,tree->left);
-	error->all("Cannot evaluate variable equal command");
+	value2 = evaluate(arg2,tree->right);
+      } else {
+	value1 = evaluate(arg1,NULL);
+	value2 = evaluate(arg2,NULL);
+	if (value2 == 0.0) error->all("Cannot evaluate variable");
 	answer = pow(value1,value2);
+      }
 
     } else if (strcmp(func,"exp") == 0) {
-      if (arg2) error->all("Cannot evaluate variable equal command");
+      if (arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      value1 = atof(strarg1);
+	tree->type = EXP;
-      answer = exp(value1);
+	tree->left = new Tree();
+	value1 = evaluate(arg1,tree->left);
+      } else answer = exp(evaluate(arg1,NULL));
 
     } else if (strcmp(func,"ln") == 0) {
-      if (arg2) error->all("Cannot evaluate variable equal command");
+      if (arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      value1 = atof(strarg1);
+	tree->type = LN;
-      if (value1 == 0.0)
+	tree->left = new Tree();
-	error->all("Cannot evaluate variable equal command");
+	value1 = evaluate(arg1,tree->left);
+      } else {
+	value1 = evaluate(arg1,NULL);
+	if (value1 == 0.0) error->all("Cannot evaluate variable");
 	answer = log(value1);
+      }
 
     } else if (strcmp(func,"sqrt") == 0) {
-      if (arg2) error->all("Cannot evaluate variable equal command");
+      if (arg2) error->all("Cannot evaluate variable");
-      strarg1 = evaluate(arg1);
+      if (tree) {
-      value1 = atof(strarg1);
+	tree->type = SQRT;
-      if (value1 == 0.0)
+	tree->left = new Tree();
-	error->all("Cannot evaluate variable equal command");
+	value1 = evaluate(arg1,tree->left);
+      } else {
+	value1 = evaluate(arg1,NULL);
+	if (value1 == 0.0) error->all("Cannot evaluate variable");
 	answer = sqrt(value1);
+      }
 
     // customize by adding group function to this list and to if statement
     //   mass(group),charge(group),xcm(group,dim),vcm(group,dim),
     //   bound(group,xmin),gyration(group)
 
     } else if (strcmp(func,"mass") == 0) {
-      if (arg2) error->all("Cannot evaluate variable equal command");
+      if (arg2) error->all("Cannot evaluate variable");
       int igroup = group->find(arg1);
-      if (igroup == -1) error->all("Variable equal group ID does not exist");
+      if (igroup == -1) error->all("Variable group ID does not exist");
       atom->check_mass();
       answer = group->mass(igroup);
 
     } else if (strcmp(func,"charge") == 0) {
-      if (arg2) error->all("Cannot evaluate variable equal command");
+      if (arg2) error->all("Cannot evaluate variable");
       int igroup = group->find(arg1);
-      if (igroup == -1) error->all("Variable equal group ID does not exist");
+      if (igroup == -1) error->all("Variable group ID does not exist");
       answer = group->charge(igroup);
 
     } else if (strcmp(func,"xcm") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
       int igroup = group->find(arg1);
-      if (igroup == -1) error->all("Variable equal group ID does not exist");
+      if (igroup == -1) error->all("Variable group ID does not exist");
       atom->check_mass();
       double masstotal = group->mass(igroup);
       double xcm[3];
@@ -574,12 +595,12 @@ char *Variable::evaluate(char *str)
       if (strcmp(arg2,"x") == 0) answer = xcm[0];
       else if (strcmp(arg2,"y") == 0) answer = xcm[1];
       else if (strcmp(arg2,"z") == 0) answer = xcm[2];
-      else error->all("Cannot evaluate variable equal command");
+      else error->all("Cannot evaluate variable");
 
     } else if (strcmp(func,"vcm") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
       int igroup = group->find(arg1);
-      if (igroup == -1) error->all("Variable equal group ID does not exist");
+      if (igroup == -1) error->all("Variable group ID does not exist");
       atom->check_mass();
       double masstotal = group->mass(igroup);
       double vcm[3];
@@ -587,12 +608,12 @@ char *Variable::evaluate(char *str)
       if (strcmp(arg2,"x") == 0) answer = vcm[0];
       else if (strcmp(arg2,"y") == 0) answer = vcm[1];
       else if (strcmp(arg2,"z") == 0) answer = vcm[2];
-      else error->all("Cannot evaluate variable equal command");
+      else error->all("Cannot evaluate variable");
 
     } else if (strcmp(func,"bound") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
       int igroup = group->find(arg1);
-      if (igroup == -1) error->all("Variable equal group ID does not exist");
+      if (igroup == -1) error->all("Variable group ID does not exist");
       double minmax[6];
       group->bounds(igroup,minmax);
       if (strcmp(arg2,"xmin") == 0) answer = minmax[0];
@@ -601,28 +622,23 @@ char *Variable::evaluate(char *str)
       else if (strcmp(arg2,"ymax") == 0) answer = minmax[3];
       else if (strcmp(arg2,"zmin") == 0) answer = minmax[4];
       else if (strcmp(arg2,"zmax") == 0) answer = minmax[5];
-      else error->all("Cannot evaluate variable equal command");
+      else error->all("Cannot evaluate variable");
 
     } else if (strcmp(func,"gyration") == 0) {
-      if (!arg2) error->all("Cannot evaluate variable equal command");
+      if (!arg2) error->all("Cannot evaluate variable");
       int igroup = group->find(arg1);
-      if (igroup == -1) error->all("Variable equal group ID does not exist");
+      if (igroup == -1) error->all("Variable group ID does not exist");
       atom->check_mass();
       double masstotal = group->mass(igroup);
       double xcm[3];
       group->xcm(igroup,masstotal,xcm);
       answer = group->gyration(igroup,masstotal,xcm);
 
-    } else error->all("Invalid math/group function in variable equal command");
+    } else error->all("Invalid math/group function in variable");
 
     delete [] func;
     delete [] arg1;
-    delete [] strarg1;
-    if (arg2) {
     delete [] arg2;
-      delete [] strarg2;
-    }
-    sprintf(result,"%.10g",answer);
 
   // string is a "vector" since contains []
   // index = everything between [] evaluated as integer
@@ -633,8 +649,7 @@ char *Variable::evaluate(char *str)
   //   grab atom-based value with index as global atom ID
 
   } else if (strchr(str,'[')) {
-    if (str[strlen(str)-1] != ']')
+    if (str[strlen(str)-1] != ']') error->all("Cannot evaluate variable");
-      error->all("Cannot evaluate variable equal command");
 
     char *ptr = strchr(str,'[');
     int n = ptr - str;
@@ -659,7 +674,7 @@ char *Variable::evaluate(char *str)
       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 equal");
+	error->all("Invalid compute ID in variable");
       delete [] id;
       modify->compute[icompute]->init();
 
@@ -672,7 +687,7 @@ char *Variable::evaluate(char *str)
 	  error->all("Variable compute ID does not compute scalar info");
 	if (modify->compute[icompute]->id_pre) {
 	  int ipre = modify->find_compute(modify->compute[icompute]->id_pre);
-	  if (ipre < 0) error->all("Could not pre-compute in variable equal");
+	  if (ipre < 0) error->all("Could not pre-compute in variable");
 	  answer = modify->compute[ipre]->compute_scalar();
 	}
 	answer = modify->compute[icompute]->compute_scalar();
@@ -683,28 +698,62 @@ char *Variable::evaluate(char *str)
 	  error->all("Variable compute ID vector is not large enough");
 	if (modify->compute[icompute]->id_pre) {
 	  int ipre = modify->find_compute(modify->compute[icompute]->id_pre);
-	  if (ipre < 0) error->all("Could not pre-compute in variable equal");
+	  if (ipre < 0) error->all("Could not pre-compute in variable");
 	  modify->compute[ipre]->compute_vector();
 	}
 	modify->compute[icompute]->compute_vector();
 	answer = modify->compute[icompute]->vector[index-1];
-      } else error->all("Invalid compute ID index in variable equal");
+      } else error->all("Invalid compute ID index in variable");
+
+    } else if (tree) {
+
+      if (strlen(arg)) error->all("Invalid atom vector argument in variable");
+
+      // customize by adding atom vector to this list and to if statement
+      // mass,x,y,z,vx,vy,vz,fx,fy,fz
+
+      tree->type = ATOMARRAY;
+      tree->nstride = 3;
+
+      if (strcmp(vector,"mass") == 0) {
+	if (atom->mass) {
+	  tree->type = TYPEARRAY;
+	  tree->array = atom->mass;
+	} else {
+	  tree->nstride = 1;
+	  tree->array = atom->rmass;
+	}
+      }
+      else if (strcmp(vector,"x") == 0) tree->array = &atom->x[0][0];
+      else if (strcmp(vector,"y") == 0) tree->array = &atom->x[0][1];
+      else if (strcmp(vector,"z") == 0) tree->array = &atom->x[0][2];
+      else if (strcmp(vector,"vx") == 0) tree->array = &atom->v[0][0];
+      else if (strcmp(vector,"vy") == 0) tree->array = &atom->v[0][1];
+      else if (strcmp(vector,"vz") == 0) tree->array = &atom->v[0][2];
+      else if (strcmp(vector,"fx") == 0) tree->array = &atom->f[0][0];
+      else if (strcmp(vector,"fy") == 0) tree->array = &atom->f[0][1];
+      else if (strcmp(vector,"fz") == 0) tree->array = &atom->f[0][2];
+      
+      else error->all("Invalid atom vector in variable");
 
     } else {
 
       if (atom->map_style == 0)
-	error->all("Cannot use atom vectors in variable "
+	error->all("Cannot use atom vector in variable unless atom map exists");
-                   "unless atom map exists");
 
       int index = atom->map(atoi(arg));
 
       // customize by adding atom vector to this list and to if statement
-      // x,y,z,vx,vy,vz,fx,fy,fz
+      // mass,x,y,z,vx,vy,vz,fx,fy,fz
 
       double mine;
       if (index >= 0 && index < atom->nlocal) {
 
-	if (strcmp(vector,"x") == 0) mine = atom->x[index][0];
+	if (strcmp(vector,"mass") == 0) {
+	  if (atom->mass) mine = atom->mass[atom->type[index]];
+	  else mine = atom->rmass[index];
+	}
+        else if (strcmp(vector,"x") == 0) mine = atom->x[index][0];
 	else if (strcmp(vector,"y") == 0) mine = atom->x[index][1];
 	else if (strcmp(vector,"z") == 0) mine = atom->x[index][2];
 	else if (strcmp(vector,"vx") == 0) mine = atom->v[index][0];
@@ -714,7 +763,7 @@ char *Variable::evaluate(char *str)
 	else if (strcmp(vector,"fy") == 0) mine = atom->f[index][1];
 	else if (strcmp(vector,"fz") == 0) mine = atom->f[index][2];
 	
-	else error->one("Invalid atom vector in variable equal command");
+	else error->one("Invalid atom vector in variable");
 	
       } else mine = 0.0;
 
@@ -723,7 +772,6 @@ char *Variable::evaluate(char *str)
 
     delete [] vector;
     delete [] arg;
-    sprintf(result,"%.10g",answer);
 
   // string is "keyword" since starts with lowercase letter
   // if keyword starts with "v_", trailing chars are variable name
@@ -731,8 +779,6 @@ char *Variable::evaluate(char *str)
   // else is thermo keyword
   //   evaluate it via evaluate_keyword()
 
-  // compute appropriate value via thermo
-
   } else if (islower(str[0])) {
 
     if (strncmp(str,"v_",2) == 0) {
@@ -741,31 +787,98 @@ char *Variable::evaluate(char *str)
       strcpy(id,&str[2]);
 
       char *v = retrieve(id);
-      if (v == NULL)
+      if (v == NULL) error->all("Invalid variable name in variable");
-	error->all("Invalid variable name in variable equal command");
       delete [] id;
 
       answer = atof(v);
 
     } else {
       int flag = output->thermo->evaluate_keyword(str,&answer);
-      if (flag) error->all("Invalid thermo keyword in variable equal command");
+      if (flag) error->all("Invalid thermo keyword in variable");
     }
 
-    sprintf(result,"%.10g",answer);
-
   // string is a number since starts with digit or "." or "-"
   // just copy to result
 
   } else if (isdigit(str[0]) || str[0] == '.' || str[0] == '-') {
 
-    strcpy(result,str);
+    answer = atof(str);
 
   // string is an error
 
-  } else error->all("Cannot evaluate variable equal command");
+  } else error->all("Cannot evaluate variable");
 
-  // return newly allocated string
+  // store answer in tree and return it
 
-  return result;
+  if (tree) tree->value = answer;
+  return answer;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Variable::build_parse_tree(int ivar)
+{
+  if (style[ivar] != ATOM) error->all("");
+  ptree = new Tree();
+  double tmp = evaluate(data[ivar][0],ptree);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Variable::evaluate_parse_tree(int igroup, double *result)
+{
+  int groupbit = group->bitmask[igroup];
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] && groupbit) result[i] = eval_tree(ptree,i);
+    else result[i] = 0.0;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Variable::free_parse_tree()
+{
+  free_tree(ptree);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double Variable::eval_tree(Tree *tree, int i)
+{
+  if (tree->type == VALUE) return tree->value;
+  if (tree->type == ATOMARRAY) return tree->array[i*tree->nstride];
+  if (tree->type == TYPEARRAY) return tree->array[atom->type[i]];
+
+  if (tree->type == ADD)
+    return eval_tree(tree->left,i) + eval_tree(tree->right,i);
+  if (tree->type == SUB)
+    return eval_tree(tree->left,i) - eval_tree(tree->right,i);
+  if (tree->type == MULT)
+    return eval_tree(tree->left,i) * eval_tree(tree->right,i);
+  if (tree->type == DIV)
+    return eval_tree(tree->left,i) / eval_tree(tree->right,i);
+  if (tree->type == NEG)
+    return -eval_tree(tree->left,i);
+  if (tree->type == POW)
+    return pow(eval_tree(tree->left,i),eval_tree(tree->right,i));
+  if (tree->type == EXP)
+    return exp(eval_tree(tree->left,i));
+  if (tree->type == LN)
+    return log(eval_tree(tree->left,i));
+  if (tree->type == SQRT)
+    return sqrt(eval_tree(tree->left,i));
+
+  return 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Variable::free_tree(Tree *tree)
+{
+  if (tree->left) free_tree(tree->left);
+  if (tree->right) free_tree(tree->right);
+  delete tree;
 }

src/variable.h

@@ -28,6 +28,10 @@ class Variable : protected Pointers {
   int find(char *);
   char *retrieve(char *);
 
+  void build_parse_tree(int);
+  void evaluate_parse_tree(int, double *);
+  void free_parse_tree();
+
  private:
   int me;
   int nvar;                // # of defined variables
@@ -38,9 +42,21 @@ class Variable : protected Pointers {
   int *index;              // next available value for each variable
   char ***data;            // str value of each variable's values
 
+  struct Tree {
+    double value;
+    double *array;
+    int nstride;
+    int type;
+    Tree *left,*right;
+  };
+
+  Tree *ptree;             // parse tree for an ATOM variable
+
   void copy(int, char **, char **);
-  char *evaluate(char *);
+  double evaluate(char *, Tree *);
   void remove(int);
+  double eval_tree(Tree *, int);
+  void free_tree(Tree *);
 };
 
 }

src/verlet.h

@@ -21,19 +21,25 @@ namespace LAMMPS_NS {
 class Verlet : public Integrate {
  public:
   Verlet(class LAMMPS *, int, char **);
+  ~Verlet();
   void init();
   void setup();
   void iterate(int);
 
  private:
-  int virial_every;       // what vflag should be on every timestep (0,1,2)
+  int virial_style;                 // compute virial explicitly or implicitly
-  int virial_thermo;      // what vflag should be on thermo steps (1,2)
+  int virial_every;                 // 1 if virial computed every step
-  int pairflag,torqueflag,granflag;
+  int next_virial;                  // next timestep to compute virial
+  int nfix_virial;                  // # of fixes that need virial occasionally
+  int *fix_virial_every;            // frequency they require it
+  int *next_fix_virial;             // next timestep they need it
 
-  int maxpair;            // copies of Update quantities
+  int pairflag,torqueflag,granflag; // arrays to zero out every step
+  int maxpair;                      // local copies of Update quantities
   double **f_pair;
 
   void force_clear(int);
+  int fix_virial(int);
 };
 
 }