add variable option to compute bond/angle/dihedral local

doc/src/compute_angle_local.txt

@@ -10,20 +10,27 @@ compute angle/local command :h3
 
 [Syntax:]
 
-compute ID group-ID angle/local value1 value2 ... :pre
+compute ID group-ID angle/local value1 value2 ... keyword args ... :pre
 
 ID, group-ID are documented in "compute"_compute.html command :ulb,l
 angle/local = style name of this compute command :l
 one or more values may be appended :l
-value = {theta} or {eng} :l
+value = {theta} or {eng} or {v_name} :l
   {theta} = tabulate angles
-  {eng} = tabulate angle energies :pre
+  {eng} = tabulate angle energies
+  {v_name} = equal-style variable with name (see below) :pre
+zero or more keyword/args pairs may be appended :l
+keyword = {set} :l
+  {set} args = theta name
+    theta = only currently allowed arg
+    name = name of variable to set with theta :pre
 :ule
 
 [Examples:]
 
 compute 1 all angle/local theta
-compute 1 all angle/local eng theta :pre
+compute 1 all angle/local eng theta 
+compute 1 all angle/local theta v_cos set theta t :pre
 
 [Description:]
 
@@ -36,6 +43,47 @@ The value {theta} is the angle for the 3 atoms in the interaction.
 
 The value {eng} is the interaction energy for the angle.
 
+The value {v_name} can be used together with the {set} keyword to
+compute a user-specified function of the angle theta.  The {name}
+specified for the {v_name} value is the name of an "equal-style
+variable"_variable.html which should evaluate a formula based on a
+variable which will store the angle theta.  This other variable must
+be an "internal-style variable"_variable.html defined in the input
+script; its initial numeric value can be anything.  It must be an
+internal-style variable, because this command resets its value
+directly.  The {set} keyword is used to identify the name of this
+other variable associated with theta.
+
+Note that the value of theta for each angle which stored in the
+internal variable is in radians, not degrees.
+
+As an example, these commands can be added to the bench/in.rhodo
+script to compute the cosine and cosine^2 of every angle in the system
+and output the statistics in various ways:
+
+variable t internal 0.0
+variable cos equal cos(v_t)
+variable cossq equal cos(v_t)*cos(v_t) :pre
+
+compute 1 all property/local aatom1 aatom2 aatom3 atype
+compute 2 all angle/local eng theta v_cos v_cossq set theta t
+dump 1 all local 100 tmp.dump c_1[*] c_2[*] :pre
+
+compute 3 all reduce ave c_2[*]
+thermo_style custom step temp press c_3[*] :pre
+
+fix 10 all ave/histo 10 10 100 -1 1 20 c_2[3] mode vector file tmp.histo :pre
+
+The "dump local"_dump.html command will output the energy, angle,
+cosine(angle), cosine^2(angle) for every angle in the system.  The
+"thermo_style"_thermo_style.html command will print the average of
+those quantities via the "compute reduce"_compute_reduce.html command
+with thermo output.  And the "fix ave/histo"_fix_ave_histo.html
+command will histogram the cosine(angle) values and write them to a
+file.
+
+:line
+
 The local data stored by this command is generated by looping over all
 the atoms owned on a processor and their angles.  An angle will only
 be included if all 3 atoms in the angle are in the specified compute
@@ -65,12 +113,12 @@ dump 1 all local 1000 tmp.dump index c_1\[1\] c_1\[2\] c_1\[3\] c_1\[4\] c_2\[1\
 [Output info:]
 
 This compute calculates a local vector or local array depending on the
-number of keywords.  The length of the vector or number of rows in the
+number of values.  The length of the vector or number of rows in the
-array is the number of angles.  If a single keyword is specified, a
+array is the number of angles.  If a single value is specified, a
-local vector is produced.  If two or more keywords are specified, a
+local vector is produced.  If two or more values are specified, a
 local array is produced where the number of columns = the number of
-keywords.  The vector or array can be accessed by any command that
+values.  The vector or array can be accessed by any command that uses
-uses local values from a compute as input.  See the "Howto
+local values from a compute as input.  See the "Howto
 output"_Howto_output.html doc page for an overview of LAMMPS output
 options.
 

doc/src/compute_bond_local.txt

@@ -10,12 +10,12 @@ compute bond/local command :h3
 
 [Syntax:]
 
-compute ID group-ID bond/local value1 value2 ... :pre
+compute ID group-ID bond/local value1 value2 ... keyword args ... :pre
 
 ID, group-ID are documented in "compute"_compute.html command :ulb,l
 bond/local = style name of this compute command :l
 one or more values may be appended :l
-value = {dist} or {engpot} or {force} or {engvib} or {engrot} or {engtrans} or {omega} or {velvib} :l
+value = {dist} or {engpot} or {force} or {engvib} or {engrot} or {engtrans} or {omega} or {velvib} or {v_name} :l
   {dist} = bond distance
   {engpot} = bond potential energy
   {force} = bond force :pre
@@ -23,13 +23,22 @@ value = {dist} or {engpot} or {force} or {engvib} or {engrot} or {engtrans} or {
   {engrot} = bond kinetic energy of rotation
   {engtrans} = bond kinetic energy of translation
   {omega} = magnitude of bond angular velocity
-  {velvib} = vibrational velocity along the bond length :pre
+  {velvib} = vibrational velocity along the bond length
+  {v_name} = equal-style variable with name (see below) :pre
+zero or more keyword/args pairs may be appended :l
+keyword = {set} :l
+  {set} args = dist name
+    dist = only currently allowed arg
+    name = name of variable to set with distance (dist) :pre
+:ule
+
 :ule
 
 [Examples:]
 
 compute 1 all bond/local engpot
 compute 1 all bond/local dist engpot force :pre
+compute 1 all angle/local dist v_distsq set dist d :pre
 
 [Description:]
 
@@ -38,6 +47,10 @@ interactions.  The number of datums generated, aggregated across all
 processors, equals the number of bonds in the system, modified by the
 group parameter as explained below.
 
+All these properties are computed for the pair of atoms in a bond,
+whether the 2 atoms represent a simple diatomic molecule, or are part
+of some larger molecule.
+
 The value {dist} is the current length of the bond.
 
 The value {engpot} is the potential energy for the bond,
@@ -79,9 +92,41 @@ two atoms in the bond towards each other.  A negative value means the
 2 atoms are moving toward each other; a positive value means they are
 moving apart.
 
-Note that all these properties are computed for the pair of atoms in a
+The value {v_name} can be used together with the {set} keyword to
-bond, whether the 2 atoms represent a simple diatomic molecule, or are
+compute a user-specified function of the bond distance.  The {name}
-part of some larger molecule.
+specified for the {v_name} value is the name of an "equal-style
+variable"_variable.html which should evaluate a formula based on a
+variable which will store the bond distance.  This other variable must
+be an "internal-style variable"_variable.html defined in the input
+script; its initial numeric value can be anything.  It must be an
+internal-style variable, because this command resets its value
+directly.  The {set} keyword is used to identify the name of this
+other variable associated with theta.
+
+As an example, these commands can be added to the bench/in.rhodo
+script to compute the distance^2 of every bond in the system and
+output the statistics in various ways:
+
+variable d internal 0.0
+variable dsq equal v_d*v_d :pre
+
+compute 1 all property/local batom1 batom2 btype
+compute 2 all bond/local engpot dist v_dsq set dist d
+dump 1 all local 100 tmp.dump c_1[*] c_2[*] :pre
+
+compute 3 all reduce ave c_2[*]
+thermo_style custom step temp press c_3[*] :pre
+
+fix 10 all ave/histo 10 10 100 0 6 20 c_2[3] mode vector file tmp.histo :pre
+
+The "dump local"_dump.html command will output the energy, distance,
+distance^2 for every bond in the system.  The
+"thermo_style"_thermo_style.html command will print the average of
+those quantities via the "compute reduce"_compute_reduce.html command
+with thermo output.  And the "fix ave/histo"_fix_ave_histo.html
+command will histogram the distance^2 values and write them to a file.
+
+:line
 
 The local data stored by this command is generated by looping over all
 the atoms owned on a processor and their bonds.  A bond will only be
@@ -111,12 +156,12 @@ dump 1 all local 1000 tmp.dump index c_1\[*\] c_2\[*\] :pre
 [Output info:]
 
 This compute calculates a local vector or local array depending on the
-number of keywords.  The length of the vector or number of rows in the
+number of values.  The length of the vector or number of rows in the
-array is the number of bonds.  If a single keyword is specified, a
+array is the number of bonds.  If a single value is specified, a local
-local vector is produced.  If two or more keywords are specified, a
+vector is produced.  If two or more values are specified, a local
-local array is produced where the number of columns = the number of
+array is produced where the number of columns = the number of values.
-keywords.  The vector or array can be accessed by any command that
+The vector or array can be accessed by any command that uses local
-uses local values from a compute as input.  See the "Howto
+values from a compute as input.  See the "Howto
 output"_Howto_output.html doc page for an overview of LAMMPS output
 options.
 

doc/src/compute_dihedral_local.txt

@@ -10,18 +10,25 @@ compute dihedral/local command :h3
 
 [Syntax:]
 
-compute ID group-ID dihedral/local value1 value2 ... :pre
+compute ID group-ID dihedral/local value1 value2 ... keyword args ... :pre
 
 ID, group-ID are documented in "compute"_compute.html command :ulb,l
 dihedral/local = style name of this compute command :l
 one or more values may be appended :l
-value = {phi} :l
+value = {phi} or {v_name} :l
-  {phi} = tabulate dihedral angles :pre
+  {phi} = tabulate dihedral angles
+  {v_name} = equal-style variable with name (see below) :pre
+zero or more keyword/args pairs may be appended :l
+keyword = {set} :l
+  {set} args = phi name
+    phi = only currently allowed arg
+    name = name of variable to set with phi :pre
 :ule
 
 [Examples:]
 
 compute 1 all dihedral/local phi :pre
+compute 1 all dihedral/local phi v_cos set phi p :pre
 
 [Description:]
 
@@ -33,6 +40,47 @@ by the group parameter as explained below.
 The value {phi} is the dihedral angle, as defined in the diagram on
 the "dihedral_style"_dihedral_style.html doc page.
 
+The value {v_name} can be used together with the {set} keyword to
+compute a user-specified function of the dihedral angle phi.  The
+{name} specified for the {v_name} value is the name of an "equal-style
+variable"_variable.html which should evaluate a formula based on a
+variable which will store the angle phi.  This other variable must
+be an "internal-style variable"_variable.html defined in the input
+script; its initial numeric value can be anything.  It must be an
+internal-style variable, because this command resets its value
+directly.  The {set} keyword is used to identify the name of this
+other variable associated with phi.
+
+Note that the value of phi for each angle which stored in the internal
+variable is in radians, not degrees.
+
+As an example, these commands can be added to the bench/in.rhodo
+script to compute the cosine and cosine^2 of every dihedral angle in
+the system and output the statistics in various ways:
+
+variable p internal 0.0
+variable cos equal cos(v_p)
+variable cossq equal cos(v_p)*cos(v_p) :pre
+
+compute 1 all property/local datom1 datom2 datom3 datom4 dtype
+compute 2 all dihedral/local phi v_cos v_cossq set phi p
+dump 1 all local 100 tmp.dump c_1[*] c_2[*] :pre
+
+compute 3 all reduce ave c_2[*]
+thermo_style custom step temp press c_3[*] :pre
+
+fix 10 all ave/histo 10 10 100 -1 1 20 c_2[2] mode vector file tmp.histo :pre
+
+The "dump local"_dump.html command will output the angle,
+cosine(angle), cosine^2(angle) for every dihedral in the system.  The
+"thermo_style"_thermo_style.html command will print the average of
+those quantities via the "compute reduce"_compute_reduce.html command
+with thermo output.  And the "fix ave/histo"_fix_ave_histo.html
+command will histogram the cosine(angle) values and write them to a
+file.
+
+:line
+
 The local data stored by this command is generated by looping over all
 the atoms owned on a processor and their dihedrals.  A dihedral will
 only be included if all 4 atoms in the dihedral are in the specified
@@ -57,12 +105,12 @@ dump 1 all local 1000 tmp.dump index c_1\[1\] c_1\[2\] c_1\[3\] c_1\[4\] c_1\[5\
 [Output info:]
 
 This compute calculates a local vector or local array depending on the
-number of keywords.  The length of the vector or number of rows in the
+number of values.  The length of the vector or number of rows in the
-array is the number of dihedrals.  If a single keyword is specified, a
+array is the number of dihedrals.  If a single value is specified, a
-local vector is produced.  If two or more keywords are specified, a
+local vector is produced.  If two or more values are specified, a
 local array is produced where the number of columns = the number of
-keywords.  The vector or array can be accessed by any command that
+values.  The vector or array can be accessed by any command that uses
-uses local values from a compute as input.  See the "Howto
+local values from a compute as input.  See the "Howto
 output"_Howto_output.html doc page for an overview of LAMMPS output
 options.
 

src/compute_angle_local.cpp

@@ -21,6 +21,8 @@
 #include "domain.h"
 #include "force.h"
 #include "angle.h"
+#include "input.h"
+#include "variable.h"
 #include "math_const.h"
 #include "memory.h"
 #include "error.h"
@@ -30,11 +32,13 @@ using namespace MathConst;
 
 #define DELTA 10000
 
+enum{THETA,ENG,VARIABLE};
+
 /* ---------------------------------------------------------------------- */
 
 ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
-  vlocal(NULL), alocal(NULL)
+  bstyle(NULL), vstr(NULL), vvar(NULL), tstr(NULL), vlocal(NULL), alocal(NULL)
 {
   if (narg < 4) error->all(FLERR,"Illegal compute angle/local command");
 
@@ -42,19 +46,82 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
     error->all(FLERR,"Compute angle/local used when angles are not allowed");
 
   local_flag = 1;
+
+  // style args
+
   nvalues = narg - 3;
+  bstyle = new int[nvalues];
+  vstr = new char*[nvalues];
+  vvar = new int[nvalues];
+  
+  nvalues = 0;
+  tflag = 0;
+  nvar = 0;
+
+  int iarg;
+  for (iarg = 3; iarg < narg; iarg++) {
+    if (strcmp(arg[iarg],"theta") == 0) {
+      bstyle[nvalues++] = THETA;
+      tflag = 1;
+    } else if (strcmp(arg[iarg],"eng") == 0) {
+      bstyle[nvalues++] = ENG;
+    } else if (strncmp(arg[iarg],"v_",2) == 0) {
+      bstyle[nvalues++] = VARIABLE;
+      int n = strlen(arg[iarg]);
+      vstr[nvar] = new char[n];
+      strcpy(vstr[nvar],&arg[iarg][2]);
+      nvar++;
+    } else break;
+  }
+
+  // optional args
+
+  setflag = 0;
+  tstr = NULL;
+  
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"set") == 0) {
+      setflag = 1;
+      if (iarg+3 > narg) error->all(FLERR,"Illegal compute angle/local command");
+      if (strcmp(arg[iarg+1],"theta") == 0) {
+        delete [] tstr;
+        int n = strlen(arg[iarg+2]) + 1;
+        tstr = new char[n];
+        strcpy(tstr,arg[iarg+2]);
+	tflag = 1;
+      } else error->all(FLERR,"Illegal compute angle/local command");
+      iarg += 3;
+    } else error->all(FLERR,"Illegal compute angle/local command");
+  }
+
+  // error check
+
+  if (nvar) {
+    if (!setflag)
+      error->all(FLERR,"Compute angle/local variable requires a set variable");
+    for (int i = 0; i < nvar; i++) {
+      vvar[i] = input->variable->find(vstr[i]);
+      if (vvar[i] < 0)
+	error->all(FLERR,"Variable name for copute angle/local does not exist");
+      if (!input->variable->equalstyle(vvar[i]))
+	error->all(FLERR,"Variable for compute angle/local is invalid style");
+    }
+
+    if (tstr) {
+      tvar = input->variable->find(tstr);
+      if (tvar < 0)
+        error->all(FLERR,"Variable name for compute angle/local does not exist");
+      if (!input->variable->internalstyle(tvar))
+        error->all(FLERR,"Variable for compute angle/local is invalid style");
+    }
+  } else if (setflag) 
+    error->all(FLERR,"Compute angle/local set with no variable");
+
+  // initialize output
+
   if (nvalues == 1) size_local_cols = 0;
   else size_local_cols = nvalues;
 
-  tflag = eflag = -1;
-  nvalues = 0;
-
-  for (int iarg = 3; iarg < narg; iarg++) {
-    if (strcmp(arg[iarg],"theta") == 0) tflag = nvalues++;
-    else if (strcmp(arg[iarg],"eng") == 0) eflag = nvalues++;
-    else error->all(FLERR,"Invalid keyword in compute angle/local command");
-  }
-
   nmax = 0;
   vlocal = NULL;
   alocal = NULL;
@@ -64,6 +131,13 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
 
 ComputeAngleLocal::~ComputeAngleLocal()
 {
+  delete [] bstyle;
+  for (int i = 0; i < nvar; i++) delete [] vstr[i];
+  delete [] vstr;
+  delete [] vvar;
+
+  delete [] tstr;
+
   memory->destroy(vlocal);
   memory->destroy(alocal);
 }
@@ -75,6 +149,20 @@ void ComputeAngleLocal::init()
   if (force->angle == NULL)
     error->all(FLERR,"No angle style is defined for compute angle/local");
 
+  if (nvar) {
+    for (int i = 0; i < nvar; i++) {
+      vvar[i] = input->variable->find(vstr[i]);
+      if (vvar[i] < 0)
+	error->all(FLERR,"Variable name for compute angle/local does not exist");
+    }
+
+    if (tstr) {
+      tvar = input->variable->find(tstr);
+      if (tvar < 0)
+        error->all(FLERR,"Variable name for compute angle/local does not exist");
+    }
+  }
+
   // do initial memory allocation so that memory_usage() is correct
 
   ncount = compute_angles(0);
@@ -109,11 +197,11 @@ void ComputeAngleLocal::compute_local()
 
 int ComputeAngleLocal::compute_angles(int flag)
 {
-  int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype;
+  int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype,ivar;
   tagint tagprev;
   double delx1,dely1,delz1,delx2,dely2,delz2;
-  double rsq1,rsq2,r1,r2,c;
+  double rsq1,rsq2,r1,r2,c,theta;
-  double *tbuf,*ebuf;
+  double *ptr;
 
   double **x = atom->x;
   tagint *tag = atom->tag;
@@ -131,17 +219,7 @@ int ComputeAngleLocal::compute_angles(int flag)
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
 
-  if (flag) {
+  // loop over all atoms and their angles
-    if (nvalues == 1) {
-      if (tflag >= 0) tbuf = vlocal;
-      if (eflag >= 0) ebuf = vlocal;
-    } else {
-      if (tflag >= 0 && alocal) tbuf = &alocal[0][tflag];
-      else tbuf = NULL;
-      if (eflag >= 0 && alocal) ebuf = &alocal[0][eflag];
-      else ebuf = NULL;
-    }
-  }
 
   Angle *angle = force->angle;
 
@@ -175,39 +253,62 @@ int ComputeAngleLocal::compute_angles(int flag)
       if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
       if (atype == 0) continue;
 
-      if (flag) {
+      if (!flag) {
-        if (tflag >= 0) {
+        m++;
-          delx1 = x[atom1][0] - x[atom2][0];
+        continue;
-          dely1 = x[atom1][1] - x[atom2][1];
+      }
-          delz1 = x[atom1][2] - x[atom2][2];
-          domain->minimum_image(delx1,dely1,delz1);
 
-          rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
+      // theta needed by one or more outputs
-          r1 = sqrt(rsq1);
 
-          delx2 = x[atom3][0] - x[atom2][0];
+      if (tflag) {
-          dely2 = x[atom3][1] - x[atom2][1];
+	delx1 = x[atom1][0] - x[atom2][0];
-          delz2 = x[atom3][2] - x[atom2][2];
+	dely1 = x[atom1][1] - x[atom2][1];
-          domain->minimum_image(delx2,dely2,delz2);
+	delz1 = x[atom1][2] - x[atom2][2];
+	domain->minimum_image(delx1,dely1,delz1);
 
-          rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+	rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
-          r2 = sqrt(rsq2);
+	r1 = sqrt(rsq1);
+	
+	delx2 = x[atom3][0] - x[atom2][0];
+	dely2 = x[atom3][1] - x[atom2][1];
+	delz2 = x[atom3][2] - x[atom2][2];
+	domain->minimum_image(delx2,dely2,delz2);
 
-          // c = cosine of angle
+	rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
+	r2 = sqrt(rsq2);
 
-          c = delx1*delx2 + dely1*dely2 + delz1*delz2;
+	// c = cosine of angle
-          c /= r1*r2;
+	// theta = angle in radians
-          if (c > 1.0) c = 1.0;
+
-          if (c < -1.0) c = -1.0;
+	c = delx1*delx2 + dely1*dely2 + delz1*delz2;
-          tbuf[n] = 180.0*acos(c)/MY_PI;
+	c /= r1*r2;
+	if (c > 1.0) c = 1.0;
+	if (c < -1.0) c = -1.0;
+	theta = acos(c);
+      }
+
+      if (nvalues == 1) ptr = &vlocal[m];
+      else ptr = alocal[m];
+      
+      if (nvar) {
+	ivar = 0;
+	if (tstr) input->variable->internal_set(tvar,theta);
+      }
+
+      for (n = 0; n < nvalues; n++) {
+	switch (bstyle[n]) {
+	case THETA:
+	  ptr[n] = 180.0*theta/MY_PI;
+	  break;
+	case ENG:
+	  if (atype > 0) ptr[n] = angle->single(atype,atom1,atom2,atom3);
+	  else ptr[n] = 0.0;
+	  break;
+	case VARIABLE:
+	  ptr[n] = input->variable->compute_equal(vvar[ivar]);
+	  ivar++;
+	  break;
         }
-
-        if (eflag >= 0) {
-          if (atype > 0)
-            ebuf[n] = angle->single(atype,atom1,atom2,atom3);
-          else ebuf[n] = 0.0;
-        }
-        n += nvalues;
       }
 
       m++;

src/compute_angle_local.h

@@ -33,8 +33,12 @@ class ComputeAngleLocal : public Compute {
   double memory_usage();
 
  private:
-  int nvalues,tflag,eflag;
+  int nvalues,nvar,ncount,setflag,tflag;
-  int ncount;
+
+  int tvar;
+  int *bstyle,*vvar;
+  char *tstr;
+  char **vstr;
 
   int nmax;
   double *vlocal;

src/compute_bond_local.cpp

@@ -21,8 +21,10 @@
 #include "domain.h"
 #include "force.h"
 #include "bond.h"
-#include "math_extra.h"
 #include "comm.h"
+#include "input.h"
+#include "variable.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "error.h"
 
@@ -31,13 +33,13 @@ using namespace LAMMPS_NS;
 #define DELTA 10000
 #define EPSILON 1.0e-12
 
-enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE};
+enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,VARIABLE};
 
 /* ---------------------------------------------------------------------- */
 
 ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
-  bstyle(NULL), vlocal(NULL), alocal(NULL)
+  bstyle(NULL), vstr(NULL), vvar(NULL), dstr(NULL), vlocal(NULL), alocal(NULL)
 {
   if (narg < 4) error->all(FLERR,"Illegal compute bond/local command");
 
@@ -47,14 +49,18 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
   local_flag = 1;
   comm_forward = 3;
 
+  // style args
+
   nvalues = narg - 3;
-  if (nvalues == 1) size_local_cols = 0;
-  else size_local_cols = nvalues;
-
   bstyle = new int[nvalues];
-
+  vstr = new char*[nvalues];
+  vvar = new int[nvalues];
+  
   nvalues = 0;
-  for (int iarg = 3; iarg < narg; iarg++) {
+  nvar = 0;
+
+  int iarg;
+  for (iarg = 3; iarg < narg; iarg++) {
     if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST;
     else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT;
     else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE;
@@ -63,9 +69,58 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
     else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS;
     else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA;
     else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB;
-    else error->all(FLERR,"Invalid keyword in compute bond/local command");
+    else if (strncmp(arg[iarg],"v_",2) == 0) {
+      bstyle[nvalues++] = VARIABLE;
+      int n = strlen(arg[iarg]);
+      vstr[nvar] = new char[n];
+      strcpy(vstr[nvar],&arg[iarg][2]);
+      nvar++;
+    } else break;
   }
 
+  // optional args
+
+  setflag = 0;
+  dstr = NULL;
+  
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"set") == 0) {
+      setflag = 1;
+      if (iarg+3 > narg) error->all(FLERR,"Illegal compute bond/local command");
+      if (strcmp(arg[iarg+1],"dist") == 0) {
+        delete [] dstr;
+        int n = strlen(arg[iarg+2]) + 1;
+        dstr = new char[n];
+        strcpy(dstr,arg[iarg+2]);
+      } else error->all(FLERR,"Illegal compute bond/local command");
+      iarg += 3;
+    } else error->all(FLERR,"Illegal compute bond/local command");
+  }
+
+  // error check
+
+  if (nvar) {
+    if (!setflag)
+      error->all(FLERR,"Compute bond/local variable requires a set variable");
+    for (int i = 0; i < nvar; i++) {
+      vvar[i] = input->variable->find(vstr[i]);
+      if (vvar[i] < 0)
+	error->all(FLERR,"Variable name for copute bond/local does not exist");
+      if (!input->variable->equalstyle(vvar[i]))
+	error->all(FLERR,"Variable for compute bond/local is invalid style");
+    }
+
+    if (dstr) {
+      dvar = input->variable->find(dstr);
+      if (dvar < 0)
+        error->all(FLERR,"Variable name for compute bond/local does not exist");
+      if (!input->variable->internalstyle(dvar))
+        error->all(FLERR,"Variable for compute bond/local is invalid style");
+    }
+  } else if (setflag) 
+    error->all(FLERR,"Compute bond/local set with no variable");
+
+  
   // set singleflag if need to call bond->single()
   // set velflag if compute any quantities based on velocities
 
@@ -77,6 +132,11 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
         bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1;
   }
 
+  // initialize output
+
+  if (nvalues == 1) size_local_cols = 0;
+  else size_local_cols = nvalues;
+
   nmax = 0;
   vlocal = NULL;
   alocal = NULL;
@@ -86,9 +146,15 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
 
 ComputeBondLocal::~ComputeBondLocal()
 {
+  delete [] bstyle;
+  for (int i = 0; i < nvar; i++) delete [] vstr[i];
+  delete [] vstr;
+  delete [] vvar;
+
+  delete [] dstr;
+
   memory->destroy(vlocal);
   memory->destroy(alocal);
-  delete [] bstyle;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -98,6 +164,20 @@ void ComputeBondLocal::init()
   if (force->bond == NULL)
     error->all(FLERR,"No bond style is defined for compute bond/local");
 
+  if (nvar) {
+    for (int i = 0; i < nvar; i++) {
+      vvar[i] = input->variable->find(vstr[i]);
+      if (vvar[i] < 0)
+	error->all(FLERR,"Variable name for compute bond/local does not exist");
+    }
+
+    if (dstr) {
+      dvar = input->variable->find(dstr);
+      if (dvar < 0)
+        error->all(FLERR,"Variable name for compute bond/local does not exist");
+    }
+  }
+
   // set ghostvelflag if need to acquire ghost atom velocities
 
   if (velflag && !comm->ghost_velocity) ghostvelflag = 1;
@@ -140,7 +220,7 @@ void ComputeBondLocal::compute_local()
 
 int ComputeBondLocal::compute_bonds(int flag)
 {
-  int i,m,n,nb,atom1,atom2,imol,iatom,btype;
+  int i,m,n,nb,atom1,atom2,imol,iatom,btype,ivar;
   tagint tagprev;
   double dx,dy,dz,rsq;
   double mass1,mass2,masstotal,invmasstotal;
@@ -297,6 +377,11 @@ int ComputeBondLocal::compute_bonds(int flag)
         if (nvalues == 1) ptr = &vlocal[m];
         else ptr = alocal[m];
 
+	if (nvar) {
+	  ivar = 0;
+	  if (dstr) input->variable->internal_set(dvar,sqrt(rsq));
+	}
+
         for (n = 0; n < nvalues; n++) {
           switch (bstyle[n]) {
           case DIST:
@@ -323,6 +408,10 @@ int ComputeBondLocal::compute_bonds(int flag)
           case VELVIB:
             ptr[n] = vvib;
             break;
+	  case VARIABLE:
+	    ptr[n] = input->variable->compute_equal(vvar[ivar]);
+	    ivar++;
+	    break;
           }
         }
       }

src/compute_bond_local.h

@@ -35,10 +35,13 @@ class ComputeBondLocal : public Compute {
   double memory_usage();
 
  private:
-  int nvalues;
+  int nvalues,nvar,ncount,setflag;
-  int ncount;
+
-  int *bstyle;
   int singleflag,velflag,ghostvelflag,initflag;
+  int dvar;
+  int *bstyle,*vvar;
+  char *dstr;
+  char **vstr;
 
   int nmax;
   double *vlocal;

src/compute_dihedral_local.cpp

@@ -21,6 +21,9 @@
 #include "domain.h"
 #include "force.h"
 #include "dihedral.h"
+#include "input.h"
+#include "variable.h"
+
 #include "math_const.h"
 #include "memory.h"
 #include "error.h"
@@ -31,11 +34,13 @@ using namespace MathConst;
 #define DELTA 10000
 #define SMALL 0.001
 
+enum{PHI,VARIABLE};
+
 /* ---------------------------------------------------------------------- */
 
 ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
   Compute(lmp, narg, arg),
-  vlocal(NULL), alocal(NULL)
+  bstyle(NULL), vstr(NULL), vvar(NULL), pstr(NULL), vlocal(NULL), alocal(NULL)
 {
   if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command");
 
@@ -44,18 +49,80 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
                "Compute dihedral/local used when dihedrals are not allowed");
 
   local_flag = 1;
+
+  // style args
+
   nvalues = narg - 3;
+  bstyle = new int[nvalues];
+  vstr = new char*[nvalues];
+  vvar = new int[nvalues];
+  
+  nvalues = 0;
+  nvar = 0;
+
+  int iarg;
+  for (iarg = 3; iarg < narg; iarg++) {
+    if (strcmp(arg[iarg],"phi") == 0) {
+      bstyle[nvalues++] = PHI;
+    } else if (strncmp(arg[iarg],"v_",2) == 0) {
+      bstyle[nvalues++] = VARIABLE;
+      int n = strlen(arg[iarg]);
+      vstr[nvar] = new char[n];
+      strcpy(vstr[nvar],&arg[iarg][2]);
+      nvar++;
+    } else break;
+  }
+
+  // optional args
+
+  setflag = 0;
+  pstr = NULL;
+  
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"set") == 0) {
+      setflag = 1;
+      if (iarg+3 > narg) 
+        error->all(FLERR,"Illegal compute dihedral/local command");
+      if (strcmp(arg[iarg+1],"phi") == 0) {
+        delete [] pstr;
+        int n = strlen(arg[iarg+2]) + 1;
+        pstr = new char[n];
+        strcpy(pstr,arg[iarg+2]);
+      } else error->all(FLERR,"Illegal compute dihedral/local command");
+      iarg += 3;
+    } else error->all(FLERR,"Illegal compute dihedral/local command");
+  }
+
+  // error check
+
+  if (nvar) {
+    if (!setflag)
+      error->all(FLERR,"Compute dihedral/local variable requires a set variable");
+    for (int i = 0; i < nvar; i++) {
+      vvar[i] = input->variable->find(vstr[i]);
+      if (vvar[i] < 0)
+	error->all(FLERR,
+                   "Variable name for copute dihedral/local does not exist");
+      if (!input->variable->equalstyle(vvar[i]))
+	error->all(FLERR,"Variable for compute dihedral/local is invalid style");
+    }
+
+    if (pstr) {
+      pvar = input->variable->find(pstr);
+      if (pvar < 0)
+        error->all(FLERR,
+                   "Variable name for compute dihedral/local does not exist");
+      if (!input->variable->internalstyle(pvar))
+        error->all(FLERR,"Variable for compute dihedral/local is invalid style");
+    }
+  } else if (setflag) 
+    error->all(FLERR,"Compute dihedral/local set with no variable");
+
+  // initialize output
+
   if (nvalues == 1) size_local_cols = 0;
   else size_local_cols = nvalues;
 
-  pflag = -1;
-  nvalues = 0;
-
-  for (int iarg = 3; iarg < narg; iarg++) {
-    if (strcmp(arg[iarg],"phi") == 0) pflag = nvalues++;
-    else error->all(FLERR,"Invalid keyword in compute dihedral/local command");
-  }
-
   nmax = 0;
   vlocal = NULL;
   alocal = NULL;
@@ -65,6 +132,13 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
 
 ComputeDihedralLocal::~ComputeDihedralLocal()
 {
+  delete [] bstyle;
+  for (int i = 0; i < nvar; i++) delete [] vstr[i];
+  delete [] vstr;
+  delete [] vvar;
+
+  delete [] pstr;
+
   memory->destroy(vlocal);
   memory->destroy(alocal);
 }
@@ -76,6 +150,22 @@ void ComputeDihedralLocal::init()
   if (force->dihedral == NULL)
     error->all(FLERR,"No dihedral style is defined for compute dihedral/local");
 
+  if (nvar) {
+    for (int i = 0; i < nvar; i++) {
+      vvar[i] = input->variable->find(vstr[i]);
+      if (vvar[i] < 0)
+	error->all(FLERR,
+                   "Variable name for compute dihedral/local does not exist");
+    }
+
+    if (pstr) {
+      pvar = input->variable->find(pstr);
+      if (pvar < 0)
+        error->all(FLERR,
+                   "Variable name for compute dihedral/local does not exist");
+    }
+  }
+
   // do initial memory allocation so that memory_usage() is correct
 
   ncount = compute_dihedrals(0);
@@ -107,12 +197,12 @@ void ComputeDihedralLocal::compute_local()
 
 int ComputeDihedralLocal::compute_dihedrals(int flag)
 {
-  int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom;
+  int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom,ivar;
   tagint tagprev;
   double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
   double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv;
-  double s,c;
+  double s,c,phi;
-  double *pbuf;
+  double *ptr;
 
   double **x = atom->x;
   tagint *tag = atom->tag;
@@ -130,14 +220,7 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
   int nlocal = atom->nlocal;
   int molecular = atom->molecular;
 
-  if (flag) {
+  // loop over all atoms and their dihedrals
-    if (nvalues == 1) {
-      if (pflag >= 0) pbuf = vlocal;
-    } else {
-      if (pflag >= 0 && alocal) pbuf = &alocal[0][pflag];
-      else pbuf = NULL;
-    }
-  }
 
   m = n = 0;
   for (atom2 = 0; atom2 < nlocal; atom2++) {
@@ -169,56 +252,75 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
       if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
       if (atom4 < 0 || !(mask[atom4] & groupbit)) continue;
 
-      if (flag) {
+      if (!flag) {
+        m++;
+        continue;
+      }
 
-        // phi calculation from dihedral style harmonic
+      // phi calculation from dihedral style harmonic
 
-        if (pflag >= 0) {
+      vb1x = x[atom1][0] - x[atom2][0];
-          vb1x = x[atom1][0] - x[atom2][0];
+      vb1y = x[atom1][1] - x[atom2][1];
-          vb1y = x[atom1][1] - x[atom2][1];
+      vb1z = x[atom1][2] - x[atom2][2];
-          vb1z = x[atom1][2] - x[atom2][2];
+      domain->minimum_image(vb1x,vb1y,vb1z);
-          domain->minimum_image(vb1x,vb1y,vb1z);
 
-          vb2x = x[atom3][0] - x[atom2][0];
+      vb2x = x[atom3][0] - x[atom2][0];
-          vb2y = x[atom3][1] - x[atom2][1];
+      vb2y = x[atom3][1] - x[atom2][1];
-          vb2z = x[atom3][2] - x[atom2][2];
+      vb2z = x[atom3][2] - x[atom2][2];
-          domain->minimum_image(vb2x,vb2y,vb2z);
+      domain->minimum_image(vb2x,vb2y,vb2z);
 
-          vb2xm = -vb2x;
+      vb2xm = -vb2x;
-          vb2ym = -vb2y;
+      vb2ym = -vb2y;
-          vb2zm = -vb2z;
+      vb2zm = -vb2z;
-          domain->minimum_image(vb2xm,vb2ym,vb2zm);
+      domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
-          vb3x = x[atom4][0] - x[atom3][0];
+      vb3x = x[atom4][0] - x[atom3][0];
-          vb3y = x[atom4][1] - x[atom3][1];
+      vb3y = x[atom4][1] - x[atom3][1];
-          vb3z = x[atom4][2] - x[atom3][2];
+      vb3z = x[atom4][2] - x[atom3][2];
-          domain->minimum_image(vb3x,vb3y,vb3z);
+      domain->minimum_image(vb3x,vb3y,vb3z);
 
-          ax = vb1y*vb2zm - vb1z*vb2ym;
+      ax = vb1y*vb2zm - vb1z*vb2ym;
-          ay = vb1z*vb2xm - vb1x*vb2zm;
+      ay = vb1z*vb2xm - vb1x*vb2zm;
-          az = vb1x*vb2ym - vb1y*vb2xm;
+      az = vb1x*vb2ym - vb1y*vb2xm;
-          bx = vb3y*vb2zm - vb3z*vb2ym;
+      bx = vb3y*vb2zm - vb3z*vb2ym;
-          by = vb3z*vb2xm - vb3x*vb2zm;
+      by = vb3z*vb2xm - vb3x*vb2zm;
-          bz = vb3x*vb2ym - vb3y*vb2xm;
+      bz = vb3x*vb2ym - vb3y*vb2xm;
 
-          rasq = ax*ax + ay*ay + az*az;
+      rasq = ax*ax + ay*ay + az*az;
-          rbsq = bx*bx + by*by + bz*bz;
+      rbsq = bx*bx + by*by + bz*bz;
-          rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
+      rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm;
-          rg = sqrt(rgsq);
+      rg = sqrt(rgsq);
 
-          ra2inv = rb2inv = 0.0;
+      ra2inv = rb2inv = 0.0;
-          if (rasq > 0) ra2inv = 1.0/rasq;
+      if (rasq > 0) ra2inv = 1.0/rasq;
-          if (rbsq > 0) rb2inv = 1.0/rbsq;
+      if (rbsq > 0) rb2inv = 1.0/rbsq;
-          rabinv = sqrt(ra2inv*rb2inv);
+      rabinv = sqrt(ra2inv*rb2inv);
+      
+      c = (ax*bx + ay*by + az*bz)*rabinv;
+      s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
+      
+      if (c > 1.0) c = 1.0;
+      if (c < -1.0) c = -1.0;
+      phi = atan2(s,c);
 
-          c = (ax*bx + ay*by + az*bz)*rabinv;
+      if (nvalues == 1) ptr = &vlocal[m];
-          s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z);
+      else ptr = alocal[m];
+      
+      if (nvar) {
+	ivar = 0;
+	if (pstr) input->variable->internal_set(pvar,phi);
+      }
 
-          if (c > 1.0) c = 1.0;
+      for (n = 0; n < nvalues; n++) {
-          if (c < -1.0) c = -1.0;
+	switch (bstyle[n]) {
-          pbuf[n] = 180.0*atan2(s,c)/MY_PI;
+	case PHI:
+	  ptr[n] = 180.0*phi/MY_PI;
+	  break;
+	case VARIABLE:
+	  ptr[n] = input->variable->compute_equal(vvar[ivar]);
+	  ivar++;
+	  break;
         }
-        n += nvalues;
       }
 
       m++;

src/compute_dihedral_local.h

@@ -33,8 +33,12 @@ class ComputeDihedralLocal : public Compute {
   double memory_usage();
 
  private:
-  int nvalues,pflag;
+  int nvalues,nvar,ncount,setflag,tflag;
-  int ncount;
+
+  int pvar;
+  int *bstyle,*vvar;
+  char *pstr;
+  char **vstr;
 
   int nmax;
   double *vlocal;