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

src/ASPHERE/compute_temp_asphere.cpp

@@ -29,6 +29,9 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+
 /* ---------------------------------------------------------------------- */
 
 ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
@@ -41,7 +44,8 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 1;
   tempflag = 1;
 
   vector = new double[6];
@@ -110,6 +114,8 @@ void ComputeTempAsphere::recount()
 
 double ComputeTempAsphere::compute_scalar()
 {
+  invoked |= INVOKED_SCALAR;
+
   double **v = atom->v;
   double **quat = atom->quat;
   double **angmom = atom->angmom;
@@ -163,6 +169,8 @@ void ComputeTempAsphere::compute_vector()
 {
   int i;
 
+  invoked |= INVOKED_VECTOR;
+
   double **v = atom->v;
   double **quat = atom->quat;
   double **angmom = atom->angmom;

src/ASPHERE/fix_npt_asphere.cpp

@@ -172,18 +172,16 @@ void FixNPTASphere::final_integrate()
 
   t_current = temperature->compute_scalar();
   if (press_couple == 0) {
-    if (ptemperature) double ptmp = ptemperature->compute_scalar();
     double tmp = pressure->compute_scalar();
   } else {
     temperature->compute_vector();
-    if (ptemperature) ptemperature->compute_vector();
     pressure->compute_vector();
   }
   couple();
 
   // trigger virial computation on next timestep
 
-  pressure->add_step(update->ntimestep+1);
+  pressure->addstep(update->ntimestep+1);
 
   // update eta_dot
 

src/CLASS2/bond_class2.cpp

@@ -198,12 +198,12 @@ void BondClass2::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void BondClass2::single(int type, double rsq, int i, int j, double &eng)
+double BondClass2::single(int type, double rsq, int i, int j)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double dr2 = dr*dr;
   double dr3 = dr2*dr;
   double dr4 = dr3*dr;
-  eng = k2[type]*dr2 + k3[type]*dr3 + k4[type]*dr4;
+  return (k2[type]*dr2 + k3[type]*dr3 + k4[type]*dr4);
 }

src/CLASS2/bond_class2.h

@@ -28,7 +28,7 @@ class BondClass2 : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void single(int, double, int, int, double &);
+  double single(int, double, int, int);
 
  private:
   double *r0,*k2,*k3,*k4;

src/CLASS2/pair_lj_class2.cpp

@@ -354,9 +354,9 @@ void PairLJClass2::read_restart_settings(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJClass2::single(int i, int j, int itype, int jtype, double rsq,
+double PairLJClass2::single(int i, int j, int itype, int jtype, double rsq,
-		       double factor_coul, double factor_lj, int eflag,
+			    double factor_coul, double factor_lj,
-		       One &one)
+			    double &fforce)
 {
   double r2inv,rinv,r3inv,r6inv,forcelj,philj;
 
@@ -365,12 +365,9 @@ void PairLJClass2::single(int i, int j, int itype, int jtype, double rsq,
   r3inv = r2inv*rinv;
   r6inv = r3inv*r3inv;
   forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
-  one.fforce = factor_lj*forcelj*r2inv;
+  fforce = factor_lj*forcelj*r2inv;
 
-  if (eflag) {
+  philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
-    philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
+    offset[itype][jtype];
-      offset[itype][jtype];
+  return factor_lj*philj;
-    one.eng_vdwl = factor_lj*philj;
-    one.eng_coul = 0.0;
-  }
 }

src/CLASS2/pair_lj_class2.h

@@ -30,7 +30,7 @@ class PairLJClass2 : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
 
  private:
   double cut_global;

src/CLASS2/pair_lj_class2_coul_cut.cpp

@@ -412,9 +412,10 @@ void PairLJClass2CoulCut::read_restart_settings(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJClass2CoulCut::single(int i, int j, int itype, int jtype,
+double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype,
-				 double rsq, double factor_coul,
+				   double rsq,
-				 double factor_lj, int eflag, One &one)
+				   double factor_coul, double factor_lj,
+				   double &fforce)
 {
   double r2inv,rinv,r3inv,r6inv,forcecoul,forcelj,phicoul,philj;
 
@@ -428,17 +429,18 @@ void PairLJClass2CoulCut::single(int i, int j, int itype, int jtype,
     r6inv = r3inv*r3inv;
     forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
   } else forcelj = 0.0;
-  one.fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
+  fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
 
-  if (eflag) {
+  double eng = 0.0;
-    if (rsq < cut_coulsq[itype][jtype]) {
+  if (rsq < cut_coulsq[itype][jtype]) {
-      phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
+    phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
-      one.eng_coul = factor_coul*phicoul;
+    eng += factor_coul*phicoul;
-    } else one.eng_coul = 0.0;
-    if (rsq < cut_ljsq[itype][jtype]) {
-      philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
-	offset[itype][jtype];
-      one.eng_vdwl = factor_lj*philj;
-    } else one.eng_vdwl = 0.0;
   }
+  if (rsq < cut_ljsq[itype][jtype]) {
+    philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
+      offset[itype][jtype];
+    eng += factor_lj*philj;
+  }
+
+  return eng;
 }

src/CLASS2/pair_lj_class2_coul_cut.h

@@ -31,7 +31,7 @@ class PairLJClass2CoulCut : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
 
  private:
   double cut_lj_global,cut_coul_global;

src/CLASS2/pair_lj_class2_coul_long.cpp

@@ -423,9 +423,10 @@ void PairLJClass2CoulLong::read_restart_settings(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJClass2CoulLong::single(int i, int j, int itype, int jtype,
+double PairLJClass2CoulLong::single(int i, int j, int itype, int jtype,
-				  double rsq, double factor_coul,
+				    double rsq,
-				  double factor_lj, int eflag, One &one)
+				    double factor_coul, double factor_lj,
+				    double &fforce)
 {
   double r2inv,r,rinv,r3inv,r6inv,grij,expm2,t,erfc,prefactor;
   double forcecoul,forcelj,phicoul,philj;
@@ -447,20 +448,21 @@ void PairLJClass2CoulLong::single(int i, int j, int itype, int jtype,
     r6inv = r3inv*r3inv;
     forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
   } else forcelj = 0.0;
-  one.fforce = (forcecoul + factor_lj*forcelj) * r2inv;
+  fforce = (forcecoul + factor_lj*forcelj) * r2inv;
 
-  if (eflag) {
+  double eng = 0.0;
-    if (rsq < cut_coulsq) {
+  if (rsq < cut_coulsq) {
-      phicoul = prefactor*erfc;
+    phicoul = prefactor*erfc;
-      if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
+    if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
-      one.eng_coul = phicoul;
+    eng += phicoul;
-    } else one.eng_coul = 0.0;
-    if (rsq < cut_ljsq[itype][jtype]) {
-      philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
-	offset[itype][jtype];
-      one.eng_vdwl = factor_lj*philj;
-    } else one.eng_vdwl = 0.0;
   }
+  if (rsq < cut_ljsq[itype][jtype]) {
+    philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
+      offset[itype][jtype];
+    eng += factor_lj*philj;
+  }
+
+  return eng;
 }
 
 /* ---------------------------------------------------------------------- */

src/CLASS2/pair_lj_class2_coul_long.h

@@ -31,7 +31,7 @@ class PairLJClass2CoulLong : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
   void *extract(char *);
 
  private:

src/COLLOID/pair_colloid.cpp

@@ -351,9 +351,8 @@ double PairColloid::init_one(int i, int j)
 
   offset[j][i] = offset[i][j] = 0.0;
   if (offset_flag) {
-    One one;
+    double tmp;
-    single(0,0,i,j,cutsq[i][j],0.0,1.0,1,one);
+    offset[j][i] = offset[i][j] = single(0,0,i,j,cutsq[i][j],0.0,1.0,tmp);
-    offset[j][i] = offset[i][j] = one.eng_vdwl;
   }
 
   return cut[i][j];
@@ -443,9 +442,9 @@ void PairColloid::read_restart_settings(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void PairColloid::single(int i, int j, int itype, int jtype, double rsq,
+double PairColloid::single(int i, int j, int itype, int jtype, double rsq,
-		       double factor_coul, double factor_lj, int eflag,
+			   double factor_coul, double factor_lj,
-		       One &one)
+			   double &fforce)
 {
   double K[9],h[4],g[4];
   double r,r2inv,r6inv,forcelj,c1,c2,phi,fR,dUR,dUA;
@@ -455,9 +454,9 @@ void PairColloid::single(int i, int j, int itype, int jtype, double rsq,
     r2inv = 1.0/rsq;
     r6inv = r2inv*r2inv*r2inv;
     forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
-    one.fforce = factor_lj*forcelj*r2inv;
+    fforce = factor_lj*forcelj*r2inv;
-    if (eflag) phi = r6inv*(r6inv*lj3[itype][jtype]-lj4[itype][jtype]) -
+    phi = r6inv*(r6inv*lj3[itype][jtype]-lj4[itype][jtype]) -
-		     offset[itype][jtype];
+      offset[itype][jtype];
     break;
 
   case SMALL_LARGE:
@@ -470,13 +469,12 @@ void PairColloid::single(int i, int j, int itype, int jtype, double rsq,
     K[3] *= K[3]*K[3];
     K[6] = K[3]*K[3];
     fR = sigma3[itype][jtype]*a12[itype][jtype]*c2*K[1]/K[3];
-    one.fforce = 4.0/15.0*r*fR*factor_lj * 
+    fforce = 4.0/15.0*r*fR*factor_lj * 
       (2.0*(K[1]+K[2])*(K[1]*(5.0*K[1]+22.0*K[2])+5.0*K[4]) * 
        sigma6[itype][jtype]/K[6] - 5.0)/K[0];
-    if (eflag) 
+    phi = 2.0/9.0*fR * 
-      phi = 2.0/9.0*fR * 
+      (1.0-(K[1]*(K[1]*(K[1]/3.0+3.0*K[2])+4.2*K[4])+K[2]*K[4]) * 
-	(1.0-(K[1]*(K[1]*(K[1]/3.0+3.0*K[2])+4.2*K[4])+K[2]*K[4]) * 
+       sigma6[itype][jtype]/K[6]) - offset[itype][jtype];
-	 sigma6[itype][jtype]/K[6]) - offset[itype][jtype];
     break;
 
   case LARGE_LARGE:
@@ -510,16 +508,11 @@ void PairColloid::single(int i, int j, int itype, int jtype, double rsq,
     dUR = phi/r + 5.0*fR*(g[0]+g[1]-g[2]-g[3]);
     dUA = -a12[itype][jtype]/3.0*r*((2.0*K[0]*K[7]+1.0)*K[7] + 
 				    (2.0*K[0]*K[8]-1.0)*K[8]);
-    one.fforce = factor_lj*(dUR+dUA)/r;
+    fforce = factor_lj*(dUR+dUA)/r;
-    
+    phi += a12[itype][jtype]/6.0*(2.0*K[0]*(K[7]+K[8])-log(K[8]/K[7])) - 
-    if (eflag)
+      offset[itype][jtype];
-      phi += a12[itype][jtype]/6.0*(2.0*K[0]*(K[7]+K[8])-log(K[8]/K[7])) - 
-	offset[itype][jtype];
     break;
   }
 
-  if (eflag) {
+  return factor_lj*phi;
-    one.eng_vdwl = factor_lj*phi;
-    one.eng_coul = 0.0;
-  }
 }

src/COLLOID/pair_colloid.h

@@ -30,7 +30,7 @@ class PairColloid : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
 
  private:
   double cut_global;

src/DIPOLE/compute_temp_dipole.cpp

@@ -23,9 +23,9 @@
 
 using namespace LAMMPS_NS;
 
-// moment of inertia for a sphere
+#define INVOKED_SCALAR 1
-
+#define INVOKED_VECTOR 2
-#define INERTIA 0.4
+#define INERTIA 0.4             // moment of inertia for a sphere
 
 /* ---------------------------------------------------------------------- */
 
@@ -39,7 +39,8 @@ ComputeTempDipole::ComputeTempDipole(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 1;
   tempflag = 1;
 
   vector = new double[6];
@@ -87,6 +88,8 @@ void ComputeTempDipole::recount()
 
 double ComputeTempDipole::compute_scalar()
 {
+  invoked |= INVOKED_SCALAR;
+
   double **v = atom->v;
   double *mass = atom->mass;
   double **omega = atom->omega;
@@ -116,6 +119,8 @@ void ComputeTempDipole::compute_vector()
 {
   int i;
 
+  invoked |= INVOKED_VECTOR;
+
   double **v = atom->v;
   double *mass = atom->mass;
   double **omega = atom->omega;

src/DPD/pair_dpd.cpp

@@ -260,7 +260,7 @@ void PairDPD::init_style()
   // using different random numbers
 
   if (force->newton_pair == 0 && comm->me == 0) error->warning(
-      "DPD potential needs newton pair on for momentum conservation");
+      "Pair style dpd needs newton pair on for momentum conservation");
 
   int irequest = neighbor->request(this);
 }
@@ -370,27 +370,21 @@ void PairDPD::read_restart_settings(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void PairDPD::single(int i, int j, int itype, int jtype, double rsq,
+double PairDPD::single(int i, int j, int itype, int jtype, double rsq,
-		     double factor_coul, double factor_dpd, int eflag,
+		       double factor_coul, double factor_dpd, double &fforce)
-		     One &one)
 {
   double r,rinv,wd,phi;
 
   r = sqrt(rsq);
   if (r < EPSILON) {
-    one.fforce = 0.0;
+    fforce = 0.0;
-    if (eflag) one.eng_vdwl = one.eng_coul = 0.0;
+    return 0.0;
-    return;
   }
 
   rinv = 1.0/r;
   wd = 1.0 - r/cut[itype][jtype];
-
+  fforce = a0[itype][jtype]*wd * factor_dpd*rinv;
-  one.fforce = a0[itype][jtype]*wd * factor_dpd*rinv;
   
-  if (eflag) {
+  phi = a0[itype][jtype] * r * (1.0 - 0.5*r/cut[itype][jtype]);
-    phi = a0[itype][jtype] * r * (1.0 - 0.5*r/cut[itype][jtype]);
+  return factor_dpd*phi;
-    one.eng_vdwl = factor_dpd*phi;
-    one.eng_coul = 0.0;
-  }
 }

src/DPD/pair_dpd.h

@@ -31,7 +31,7 @@ class PairDPD : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
 
  private:
   double cut_global,temperature;

src/KSPACE/pair_buck_coul_long.cpp

@@ -394,9 +394,10 @@ void PairBuckCoulLong::read_restart_settings(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void PairBuckCoulLong::single(int i, int j, int itype, int jtype,
+double PairBuckCoulLong::single(int i, int j, int itype, int jtype,
-			     double rsq, double factor_coul, double factor_lj,
+				double rsq, 
-			     int eflag, One &one)
+				double factor_coul, double factor_lj,
+				double &fforce)
 {
   double r2inv,r6inv,r,rexp,grij,expm2,t,erfc,prefactor;
   double forcecoul,forcebuck,phicoul,phibuck;
@@ -418,20 +419,20 @@ void PairBuckCoulLong::single(int i, int j, int itype, int jtype,
     rexp = exp(-r*rhoinv[itype][jtype]);
     forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv;
   } else forcebuck = 0.0;
-  one.fforce = (forcecoul + factor_lj*forcebuck) * r2inv;
+  fforce = (forcecoul + factor_lj*forcebuck) * r2inv;
   
-  if (eflag) {
+  double eng = 0.0;
-    if (rsq < cut_coulsq) {
+  if (rsq < cut_coulsq) {
-      phicoul = prefactor*erfc;
+    phicoul = prefactor*erfc;
-      if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
+    if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
-      one.eng_coul = phicoul;
+    eng += phicoul;
-    } else one.eng_coul = 0.0;
-    if (rsq < cut_ljsq[itype][jtype]) {
-      phibuck = a[itype][jtype]*rexp - c[itype][jtype]*r6inv -
-	offset[itype][jtype];
-      one.eng_vdwl = factor_lj*phibuck;
-    } else one.eng_vdwl = 0.0;
   }
+  if (rsq < cut_ljsq[itype][jtype]) {
+    phibuck = a[itype][jtype]*rexp - c[itype][jtype]*r6inv -
+      offset[itype][jtype];
+    eng += factor_lj*phibuck;
+  }
+  return eng;
 }
 
 /* ---------------------------------------------------------------------- */

src/KSPACE/pair_buck_coul_long.h

@@ -31,7 +31,7 @@ class PairBuckCoulLong : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
   void *extract(char *);
 
  private:

src/KSPACE/pair_coul_long.cpp

@@ -512,10 +512,10 @@ void PairCoulLong::free_tables()
 
 /* ---------------------------------------------------------------------- */
 
-void PairCoulLong::single(int i, int j, int itype, int jtype,
+double PairCoulLong::single(int i, int j, int itype, int jtype,
-			       double rsq,
+			    double rsq,
-			       double factor_coul, double factor_lj,
+			    double factor_coul, double factor_lj,
-			       int eflag, One &one)
+			    double &fforce)
 {
   double r2inv,r,grij,expm2,t,erfc,prefactor;
   double fraction,table,forcecoul,phicoul;
@@ -545,19 +545,17 @@ void PairCoulLong::single(int i, int j, int itype, int jtype,
       forcecoul -= (1.0-factor_coul)*prefactor;
     }
   }
-  one.fforce = forcecoul * r2inv;
+  fforce = forcecoul * r2inv;
   
-  if (eflag) {
+  if (!ncoultablebits || rsq <= tabinnersq)
-    if (!ncoultablebits || rsq <= tabinnersq)
+    phicoul = prefactor*erfc;
-      phicoul = prefactor*erfc;
+  else {
-    else {
+    table = etable[itable] + fraction*detable[itable];
-      table = etable[itable] + fraction*detable[itable];
+    phicoul = atom->q[i]*atom->q[j] * table;
-      phicoul = atom->q[i]*atom->q[j] * table;
-    }
-    if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
-    one.eng_coul = phicoul;
-    one.eng_vdwl = 0.0;
   }
+  if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
+
+  return phicoul;
 }
 
 /* ---------------------------------------------------------------------- */

src/KSPACE/pair_coul_long.h

@@ -31,7 +31,7 @@ class PairCoulLong : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
   void *extract(char *);
 
  private:

src/KSPACE/pair_lj_charmm_coul_long.cpp

@@ -1127,10 +1127,10 @@ void PairLJCharmmCoulLong::free_tables()
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
+double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
-				  double rsq,
+				    double rsq,
-				  double factor_coul, double factor_lj,
+				    double factor_coul, double factor_lj,
-				  int eflag, One &one)
+				    double &fforce)
 {
   double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor;
   double switch1,switch2,fraction,table,forcecoul,forcelj,phicoul,philj;
@@ -1174,29 +1174,31 @@ void PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
       forcelj = forcelj*switch1 + philj*switch2;
     }
   } else forcelj = 0.0;
-  one.fforce = (forcecoul + factor_lj*forcelj) * r2inv;
+  fforce = (forcecoul + factor_lj*forcelj) * r2inv;
   
-  if (eflag) {
+  double eng = 0.0;
-    if (rsq < cut_coulsq) {
+  if (rsq < cut_coulsq) {
-      if (!ncoultablebits || rsq <= tabinnersq)
+    if (!ncoultablebits || rsq <= tabinnersq)
-	phicoul = prefactor*erfc;
+      phicoul = prefactor*erfc;
-      else {
+    else {
-	table = etable[itable] + fraction*detable[itable];
+      table = etable[itable] + fraction*detable[itable];
-	phicoul = atom->q[i]*atom->q[j] * table;
+      phicoul = atom->q[i]*atom->q[j] * table;
-      }
+    }
-      if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
+    if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
-      one.eng_coul = phicoul;
+    eng += phicoul;
-    } else one.eng_coul = 0.0;
-    if (rsq < cut_ljsq) {
-      philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
-      if (rsq > cut_lj_innersq) {
-	switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
-	  (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
-	philj *= switch1;
-      }
-      one.eng_vdwl = factor_lj*philj;
-    } else one.eng_vdwl = 0.0;
   }
+
+  if (rsq < cut_ljsq) {
+    philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+    if (rsq > cut_lj_innersq) {
+      switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+	(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+      philj *= switch1;
+    }
+    eng += factor_lj*philj;
+  }
+
+  return eng;
 }
 
 /* ---------------------------------------------------------------------- */

src/KSPACE/pair_lj_charmm_coul_long.h

@@ -32,7 +32,7 @@ class PairLJCharmmCoulLong : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
 
   void compute_inner();
   void compute_middle();

src/KSPACE/pair_lj_cut_coul_long.cpp

@@ -1080,10 +1080,10 @@ void PairLJCutCoulLong::free_tables()
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJCutCoulLong::single(int i, int j, int itype, int jtype,
+double PairLJCutCoulLong::single(int i, int j, int itype, int jtype,
-			       double rsq,
+				 double rsq,
-			       double factor_coul, double factor_lj,
+				 double factor_coul, double factor_lj,
-			       int eflag, One &one)
+				 double &fforce)
 {
   double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor;
   double fraction,table,forcecoul,forcelj,phicoul,philj;
@@ -1115,29 +1115,33 @@ void PairLJCutCoulLong::single(int i, int j, int itype, int jtype,
       }
     }
   } else forcecoul = 0.0;
+
   if (rsq < cut_ljsq[itype][jtype]) {
     r6inv = r2inv*r2inv*r2inv;
     forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
   } else forcelj = 0.0;
-  one.fforce = (forcecoul + factor_lj*forcelj) * r2inv;
+
-  
+  fforce = (forcecoul + factor_lj*forcelj) * r2inv;
-  if (eflag) {
+
-    if (rsq < cut_coulsq) {
+  double eng = 0.0;
-      if (!ncoultablebits || rsq <= tabinnersq)
+  if (rsq < cut_coulsq) {
-	phicoul = prefactor*erfc;
+    if (!ncoultablebits || rsq <= tabinnersq)
-      else {
+      phicoul = prefactor*erfc;
-	table = etable[itable] + fraction*detable[itable];
+    else {
-	phicoul = atom->q[i]*atom->q[j] * table;
+      table = etable[itable] + fraction*detable[itable];
-      }
+      phicoul = atom->q[i]*atom->q[j] * table;
-      if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
+    }
-      one.eng_coul = phicoul;
+    if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
-    } else one.eng_coul = 0.0;
+    eng += phicoul;
-    if (rsq < cut_ljsq[itype][jtype]) {
-      philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
-	offset[itype][jtype];
-      one.eng_vdwl = factor_lj*philj;
-    } else one.eng_vdwl = 0.0;
   }
+
+  if (rsq < cut_ljsq[itype][jtype]) {
+    philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
+      offset[itype][jtype];
+    eng += factor_lj*philj;
+  }
+
+  return eng;
 }
 
 /* ---------------------------------------------------------------------- */

src/KSPACE/pair_lj_cut_coul_long.h

@@ -32,7 +32,7 @@ class PairLJCutCoulLong : public Pair {
   void read_restart(FILE *);
   virtual void write_restart_settings(FILE *);
   virtual void read_restart_settings(FILE *);
-  virtual void single(int, int, int, int, double, double, double, int, One &);
+  virtual double single(int, int, int, int, double, double, double, double &);
 
   void compute_inner();
   void compute_middle();

src/MANYBODY/pair_airebo.h

@@ -21,7 +21,7 @@ namespace LAMMPS_NS {
 class PairAIREBO : public Pair {
  public:
   PairAIREBO(class LAMMPS *);
-  virtual ~PairAIREBO();
+  ~PairAIREBO();
   void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
@@ -31,7 +31,6 @@ class PairAIREBO : public Pair {
   void read_restart(FILE *) {}
   void write_restart_settings(FILE *) {}
   void read_restart_settings(FILE *) {}
-  void single(int, int, int, int, double, double, double, int, One &) {}
 
  private:
   int me;

src/MANYBODY/pair_eam.cpp

@@ -764,9 +764,9 @@ void PairEAM::grab(FILE *fp, int n, double *list)
 
 /* ---------------------------------------------------------------------- */
 
-void PairEAM::single(int i, int j, int itype, int jtype,
+double PairEAM::single(int i, int j, int itype, int jtype,
-		      double rsq, double factor_coul, double factor_lj,
+		       double rsq, double factor_coul, double factor_lj,
-		      int eflag, One &one)
+		       double &fforce)
 {
   int m;
   double r,p,rhoip,rhojp,z2,z2p,recip,phi,phip,psip;
@@ -791,12 +791,9 @@ void PairEAM::single(int i, int j, int itype, int jtype,
   phi = z2*recip;
   phip = z2p*recip - phi*recip;
   psip = fp[i]*rhojp + fp[j]*rhoip + phip;
-  one.fforce = -psip*recip;
+  fforce = -psip*recip;
 
-  if (eflag) {
+  return phi;
-    one.eng_vdwl = phi;
-    one.eng_coul = 0.0;
-  }
 }
 
 /* ---------------------------------------------------------------------- */

src/MANYBODY/pair_eam.h

@@ -27,7 +27,7 @@ class PairEAM : public Pair {
   virtual void coeff(int, char **);
   void init_style();
   double init_one(int, int);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
 
   int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);

src/MOLECULE/angle_hybrid.cpp

@@ -259,9 +259,8 @@ void AngleHybrid::read_restart(FILE *fp)
 
 double AngleHybrid::single(int type, int i1, int i2, int i3)
 {
-  if (styles[map[type]]) 
+  if (styles[map[type]]) return styles[map[type]]->single(type,i1,i2,i3);
-    return styles[map[type]]->single(type,i1,i2,i3);
+  return 0.0;
-  else return 0.0;
 }
 
 /* ----------------------------------------------------------------------

src/MOLECULE/bond_fene.cpp

@@ -217,7 +217,7 @@ void BondFENE::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void BondFENE::single(int type, double rsq, int i, int j, double &eng)
+double BondFENE::single(int type, double rsq, int i, int j)
 {
   double r0sq = r0[type] * r0[type];
   double rlogarg = 1.0 - rsq/r0sq;
@@ -234,11 +234,13 @@ void BondFENE::single(int type, double rsq, int i, int j, double &eng)
     rlogarg = 0.1;
   }
 
-  eng = -0.5 * k[type]*r0sq*log(rlogarg);
+  double eng = -0.5 * k[type]*r0sq*log(rlogarg);
   if (rsq < TWO_1_3*sigma[type]*sigma[type]) {
     double sr2,sr6;
     sr2 = sigma[type]*sigma[type]/rsq;
     sr6 = sr2*sr2*sr2;
     eng += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
   }
+
+  return eng;
 }

src/MOLECULE/bond_fene.h

@@ -28,7 +28,7 @@ class BondFENE : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void single(int, double, int, int, double &);
+  double single(int, double, int, int);
 
  private:
   double TWO_1_3;

src/MOLECULE/bond_fene_expand.cpp

@@ -228,7 +228,7 @@ void BondFENEExpand::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void BondFENEExpand::single(int type, double rsq, int i, int j, double &eng)
+double BondFENEExpand::single(int type, double rsq, int i, int j)
 {
   double r = sqrt(rsq);
   double rshift = r - shift[type];
@@ -248,11 +248,13 @@ void BondFENEExpand::single(int type, double rsq, int i, int j, double &eng)
     rlogarg = 0.1;
   }
 
-  eng = -0.5 * k[type]*r0sq*log(rlogarg);
+  double eng = -0.5 * k[type]*r0sq*log(rlogarg);
   if (rshiftsq < TWO_1_3*sigma[type]*sigma[type]) {
     double sr2,sr6;
     sr2 = sigma[type]*sigma[type]/rshiftsq;
     sr6 = sr2*sr2*sr2;
     eng += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
   }
+
+  return eng;
 }

src/MOLECULE/bond_fene_expand.h

@@ -28,7 +28,7 @@ class BondFENEExpand : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void single(int, double, int, int, double &);
+  double single(int, double, int, int);
 
  private:
   double TWO_1_3;

src/MOLECULE/bond_harmonic.cpp

@@ -177,10 +177,10 @@ void BondHarmonic::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void BondHarmonic::single(int type, double rsq, int i, int j, double &eng)
+double BondHarmonic::single(int type, double rsq, int i, int j)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double rk = k[type] * dr;
-  eng = rk*dr;
+  return rk*dr;
 }

src/MOLECULE/bond_harmonic.h

@@ -28,7 +28,7 @@ class BondHarmonic : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void single(int, double, int, int, double &);
+  double single(int, double, int, int);
 
  private:
   double *k,*r0;

src/MOLECULE/bond_hybrid.cpp

@@ -263,9 +263,10 @@ void BondHybrid::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void BondHybrid::single(int type, double rsq, int i, int j, double &eng)
+double BondHybrid::single(int type, double rsq, int i, int j)
 {
-  if (styles[map[type]]) styles[map[type]]->single(type,rsq,i,j,eng);
+  if (styles[map[type]]) return styles[map[type]]->single(type,rsq,i,j);
+  return 0.0;
 }
 
 /* ----------------------------------------------------------------------

src/MOLECULE/bond_morse.cpp

@@ -187,10 +187,10 @@ void BondMorse::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void BondMorse::single(int type, double rsq, int i, int j, double &eng)
+double BondMorse::single(int type, double rsq, int i, int j)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double ralpha = exp(-alpha[type]*dr);
-  eng = d0[type]*(1-ralpha)*(1-ralpha);
+  return d0[type]*(1-ralpha)*(1-ralpha);
 }

src/MOLECULE/bond_morse.h

@@ -28,7 +28,7 @@ class BondMorse : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void single(int, double, int, int, double &);
+  double single(int, double, int, int);
 
  private:
   double *d0,*alpha,*r0;

src/MOLECULE/bond_nonlinear.cpp

@@ -182,12 +182,12 @@ void BondNonlinear::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void BondNonlinear::single(int type, double rsq, int i, int j, double &eng)
+double BondNonlinear::single(int type, double rsq, int i, int j)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double drsq = dr*dr;
   double lamdasq = lamda[type]*lamda[type];
   double denom = lamdasq - drsq;
-  eng = epsilon[type] * drsq / denom;
+  return epsilon[type] * drsq / denom;
 }

src/MOLECULE/bond_nonlinear.h

@@ -28,7 +28,7 @@ class BondNonlinear : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void single(int, double, int, int, double &);
+  double single(int, double, int, int);
 
  private:
   double *epsilon,*r0,*lamda;

src/MOLECULE/bond_quartic.cpp

@@ -58,7 +58,6 @@ void BondQuartic::compute(int eflag, int vflag)
   int i1,i2,n,m,type,itype,jtype;
   double delx,dely,delz,ebond,fbond,evdwl,fpair;
   double r,rsq,dr,r2,ra,rb,sr2,sr6;
-  Pair::One one;
 
   ebond = evdwl = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
@@ -158,9 +157,8 @@ void BondQuartic::compute(int eflag, int vflag)
     jtype = atom->type[i2];
 
     if (rsq < cutsq[itype][jtype]) {
-      force->pair->single(i1,i2,itype,jtype,rsq,1.0,1.0,eflag,one);
+      evdwl = -force->pair->single(i1,i2,itype,jtype,rsq,1.0,1.0,fpair);
-      fpair = -one.fforce;
+      fpair = -fpair;
-      if (eflag) evdwl = -(one.eng_vdwl + one.eng_coul);
 
       if (newton_bond || i1 < nlocal) {
 	f[i1][0] += delx*fpair;
@@ -299,12 +297,12 @@ void BondQuartic::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void BondQuartic::single(int type, double rsq, int i, int j, double &eng)
+double BondQuartic::single(int type, double rsq, int i, int j)
 {
   double r,dr,r2,ra,rb,sr2,sr6;
 
-  eng = 0.0;
+  if (type <= 0) return 0.0;
-  if (type <= 0) return;
+  double eng = 0.0;
 
   // subtract out pairwise contribution from 2 atoms via pair->single()
   // required since special_bond = 1,1,1
@@ -313,9 +311,8 @@ void BondQuartic::single(int type, double rsq, int i, int j, double &eng)
   int jtype = atom->type[j];
   
   if (rsq < force->pair->cutsq[itype][jtype]) {
-    Pair::One one;
+    double tmp;
-    force->pair->single(i,j,itype,jtype,rsq,1.0,1.0,1,one);
+    eng = -force->pair->single(i,j,itype,jtype,rsq,1.0,1.0,tmp);
-    eng = -one.eng_coul - one.eng_vdwl;
   }
 
   // quartic bond
@@ -334,4 +331,6 @@ void BondQuartic::single(int type, double rsq, int i, int j, double &eng)
     sr6 = sr2*sr2*sr2;
     eng += 4.0*sr6*(sr6-1.0) + 1.0;
   }
+
+  return eng;
 }

src/MOLECULE/bond_quartic.h

@@ -29,7 +29,7 @@ class BondQuartic : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void single(int, double, int, int, double &);
+  double single(int, double, int, int);
 
  private:
   double TWO_1_3;

src/MOLECULE/improper_harmonic.cpp

@@ -162,9 +162,9 @@ void ImproperHarmonic::compute(int eflag, int vflag)
     sy2  = a22*vb2y + a23*vb3y + a12*vb1y;
     sz2  = a22*vb2z + a23*vb3z + a12*vb1z;
 
-    f1[0]  = a12*vb2x + a13*vb3x + a11*vb1x;
+    f1[0] = a12*vb2x + a13*vb3x + a11*vb1x;
-    f1[1]  = a12*vb2y + a13*vb3y + a11*vb1y;
+    f1[1] = a12*vb2y + a13*vb3y + a11*vb1y;
-    f1[2]  = a12*vb2z + a13*vb3z + a11*vb1z;
+    f1[2] = a12*vb2z + a13*vb3z + a11*vb1z;
 
     f2[0] = -sx2 - f1[0];
     f2[1] = -sy2 - f1[1];

src/MOLECULE/pair_lj_charmm_coul_charmm.cpp

@@ -436,10 +436,10 @@ void PairLJCharmmCoulCharmm::read_restart_settings(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJCharmmCoulCharmm::single(int i, int j, int itype, int jtype,
+double PairLJCharmmCoulCharmm::single(int i, int j, int itype, int jtype,
-				    double rsq, double factor_coul,
+				      double rsq,
-				    double factor_lj,
+				      double factor_coul, double factor_lj,
-				    int eflag, One &one)
+				      double &fforce)
 {
   double r2inv,r6inv,switch1,switch2,forcecoul,forcelj,phicoul,philj;
 
@@ -466,29 +466,30 @@ void PairLJCharmmCoulCharmm::single(int i, int j, int itype, int jtype,
       forcelj = forcelj*switch1 + philj*switch2;
     }
   } else forcelj = 0.0;
-  one.fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
+  fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
-  
+
-  if (eflag) {
+  double eng = 0.0;
-    if (rsq < cut_coulsq) {
+  if (rsq < cut_coulsq) {
-      phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
+    phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
-      if (rsq > cut_coul_innersq) {
+    if (rsq > cut_coul_innersq) {
-	switch1 = (cut_coulsq-rsq) * (cut_coulsq-rsq) *
+      switch1 = (cut_coulsq-rsq) * (cut_coulsq-rsq) *
-	  (cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) / 
+	(cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) / 
-	  denom_coul;
+	denom_coul;
-	phicoul *= switch1;
+      phicoul *= switch1;
-      }
+    }
-      one.eng_coul = factor_coul*phicoul;
+    eng += factor_coul*phicoul;
-    } else one.eng_coul = 0.0;
-    if (rsq < cut_ljsq) {
-      philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
-      if (rsq > cut_lj_innersq) {
-	switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
-	  (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
-	philj *= switch1;
-      }
-      one.eng_vdwl = factor_lj*philj;
-    } else one.eng_vdwl = 0.0;
   }
+  if (rsq < cut_ljsq) {
+    philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+    if (rsq > cut_lj_innersq) {
+      switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+	(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+      philj *= switch1;
+    }
+    eng += factor_lj*philj;
+  }
+
+  return eng;
 }
 
 /* ---------------------------------------------------------------------- */

src/MOLECULE/pair_lj_charmm_coul_charmm.h

@@ -31,7 +31,7 @@ class PairLJCharmmCoulCharmm : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  virtual void single(int, int, int, int, double, double, double, int, One &);
+  virtual double single(int, int, int, int, double, double, double, double &);
   virtual void *extract(char *);
 
  protected:

src/MOLECULE/pair_lj_charmm_coul_charmm_implicit.cpp

@@ -157,10 +157,12 @@ void PairLJCharmmCoulCharmmImplicit::compute(int eflag, int vflag)
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJCharmmCoulCharmmImplicit::single(int i, int j, int itype, int jtype,
+double PairLJCharmmCoulCharmmImplicit::single(int i, int j,
-					    double rsq, double factor_coul,
+					      int itype, int jtype,
-					    double factor_lj,
+					      double rsq,
-					    int eflag, One &one)
+					      double factor_coul,
+					      double factor_lj,
+					      double &fforce)
 {
   double r2inv,r6inv,switch1,switch2,forcecoul,forcelj,phicoul,philj;
 
@@ -187,29 +189,30 @@ void PairLJCharmmCoulCharmmImplicit::single(int i, int j, int itype, int jtype,
       forcelj = forcelj*switch1 + philj*switch2;
     }
   } else forcelj = 0.0;
-  one.fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
+  fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
-  
+
-  if (eflag) {
+  double eng = 0.0;
-    if (rsq < cut_coulsq) {
+  if (rsq < cut_coulsq) {
-      phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*r2inv;
+    phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*r2inv;
-      if (rsq > cut_coul_innersq) {
+    if (rsq > cut_coul_innersq) {
-	switch1 = (cut_coulsq-rsq) * (cut_coulsq-rsq) *
+      switch1 = (cut_coulsq-rsq) * (cut_coulsq-rsq) *
-	  (cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) / 
+	(cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) / 
-	  denom_coul;
+	denom_coul;
-	phicoul *= switch1;
+      phicoul *= switch1;
-      }
+    }
-      one.eng_coul = factor_coul*phicoul;
+    eng += factor_coul*phicoul;
-    } else one.eng_coul = 0.0;
-    if (rsq < cut_ljsq) {
-      philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
-      if (rsq > cut_lj_innersq) {
-	switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
-	  (cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
-	philj *= switch1;
-      }
-      one.eng_vdwl = factor_lj*philj;
-    } else one.eng_vdwl = 0.0;
   }
+  if (rsq < cut_ljsq) {
+    philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+    if (rsq > cut_lj_innersq) {
+      switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
+	(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
+      philj *= switch1;
+    }
+    eng += factor_lj*philj;
+  }
+
+  return eng;
 }
 
 /* ---------------------------------------------------------------------- */

src/MOLECULE/pair_lj_charmm_coul_charmm_implicit.h

@@ -22,7 +22,7 @@ class PairLJCharmmCoulCharmmImplicit : public PairLJCharmmCoulCharmm {
  public:
   PairLJCharmmCoulCharmmImplicit(class LAMMPS *);
   void compute(int, int);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
   void *extract(char *);
 };
 

src/USER-ACKLAND/compute_ackland_atom.cpp

@@ -34,6 +34,7 @@
 using namespace LAMMPS_NS;
 
 enum{UNKNOWN,BCC,FCC,HCP,ICO};
+#define INVOKED_PERATOM 4
 
 /* ---------------------------------------------------------------------- */
 
@@ -93,6 +94,7 @@ void ComputeAcklandAtom::init_list(int id, NeighList *ptr)
 }
 
 /* ---------------------------------------------------------------------- */
+
 void ComputeAcklandAtom::compute_peratom()
 {
   int i,j,ii,jj,k,n,inum,jnum;
@@ -100,6 +102,8 @@ void ComputeAcklandAtom::compute_peratom()
   int *ilist,*jlist,*numneigh,**firstneigh;
   int chi[8];
 
+  invoked |= INVOKED_PERATOM;
+
   // grow structure array if necessary
 
   if (atom->nlocal > nmax) {

src/USER-EWALDN/pair_buck_coul.cpp

@@ -1116,24 +1116,24 @@ void PairBuckCoul::free_tables()
 
 /* ---------------------------------------------------------------------- */
 
-void PairBuckCoul::single(int i, int j, int itype, int jtype,
+double PairBuckCoul::single(int i, int j, int itype, int jtype,
-			       double rsq,
+			    double rsq, double factor_coul, double factor_buck,
-			       double factor_coul, double factor_buck,
+			    double &fforce)
-			       int eflag, One &one)
 {
   double f, r, r2inv, r6inv, force_coul, force_buck;
   double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2, *q = atom->q;
 
   r = sqrt(rsq);
   r2inv = 1.0/rsq;
-  one.eng_coul = 0.0;
+  double eng = 0.0;
+
   if ((ewald_order&2) && (rsq < cut_coulsq)) {		// coulombic
     if (!ncoultablebits || rsq <= tabinnersq) {		// series real space
       register double x = g_ewald*r;
       register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
       f = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x);
       force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f;
-      if (eflag) one.eng_coul = t-f;
+      eng += t-f;
     }
     else {						// table real space
       register float t = rsq;
@@ -1141,12 +1141,10 @@ void PairBuckCoul::single(int i, int j, int itype, int jtype,
       register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
       t = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
       force_coul = qiqj*(ftable[k]+f*dftable[k]-t);
-      if (eflag) one.eng_coul = qiqj*(etable[k]+f*detable[k]-t);
+      eng += qiqj*(etable[k]+f*detable[k]-t);
     }
-  }
+  } else force_coul = 0.0;
-  else force_coul = 0.0;
   
-  one.eng_vdwl = 0.0;
   if (rsq < cut_bucksq[itype][jtype]) {			// buckingham
     register double expr = factor_buck*exp(-sqrt(rsq)*rhoinv[itype][jtype]);
     r6inv = r2inv*r2inv*r2inv;
@@ -1155,18 +1153,18 @@ void PairBuckCoul::single(int i, int j, int itype, int jtype,
       x2 = a2*exp(-x2)*buck_c[itype][jtype];
       force_buck = buck1[itype][jtype]*r*expr-
        	g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*buck2[itype][jtype];
-      if (eflag) one.eng_vdwl = buck_a[itype][jtype]*expr-
+      eng += buck_a[itype][jtype]*expr-
 	g6*((a2+1.0)*a2+0.5)*x2+t*buck_c[itype][jtype];
     }
     else {						// cut
       force_buck = 
 	buck1[itype][jtype]*r*expr-factor_buck*buck_c[itype][jtype]*r6inv;
-      if (eflag) one.eng_vdwl = buck_a[itype][jtype]*expr-
+      eng += buck_a[itype][jtype]*expr-
-	    factor_buck*(buck_c[itype][jtype]*r6inv-offset[itype][jtype]);
+	factor_buck*(buck_c[itype][jtype]*r6inv-offset[itype][jtype]);
     }
-  } 
+  } else force_buck = 0.0;
-  else force_buck = 0.0;
 
-  one.fforce = (force_coul+force_buck)*r2inv;
+  fforce = (force_coul+force_buck)*r2inv;
+  return eng;
 }
 

src/USER-EWALDN/pair_buck_coul.h

@@ -36,7 +36,7 @@ class PairBuckCoul : public Pair {
   
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
   void *extract(char *);
 
   void compute_inner();

src/USER-EWALDN/pair_lj_coul.cpp

@@ -1107,23 +1107,23 @@ void PairLJCoul::free_tables()
 
 /* ---------------------------------------------------------------------- */
 
-void PairLJCoul::single(int i, int j, int itype, int jtype,
+double PairLJCoul::single(int i, int j, int itype, int jtype,
-			       double rsq,
+			  double rsq, double factor_coul, double factor_lj,
-			       double factor_coul, double factor_lj,
+			  double &fforce)
-			       int eflag, One &one)
 {
   double r2inv, r6inv, force_coul, force_lj;
   double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2, *q = atom->q;
 
+  double eng = 0.0;
+
   r2inv = 1.0/rsq;
-  one.eng_coul = 0.0;
   if ((ewald_order&2) && (rsq < cut_coulsq)) {		// coulombic
     if (!ncoultablebits || rsq <= tabinnersq) {		// series real space
       register double r = sqrt(rsq), x = g_ewald*r;
       register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
       r = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x);
       force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r;
-      if (eflag) one.eng_coul = t-r;
+      eng += t-r;
     }
     else {						// table real space
       register float t = rsq;
@@ -1131,12 +1131,10 @@ void PairLJCoul::single(int i, int j, int itype, int jtype,
       register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
       t = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
       force_coul = qiqj*(ftable[k]+f*dftable[k]-t);
-      if (eflag) one.eng_coul = qiqj*(etable[k]+f*detable[k]-t);
+      eng += qiqj*(etable[k]+f*detable[k]-t);
     }
-  }
+  } else force_coul = 0.0;
-  else force_coul = 0.0;
   
-  one.eng_vdwl = 0.0;
   if (rsq < cut_ljsq[itype][jtype]) {			// lennard-jones
     r6inv = r2inv*r2inv*r2inv;
     if (ewald_order&64) {				// long-range
@@ -1144,17 +1142,16 @@ void PairLJCoul::single(int i, int j, int itype, int jtype,
       x2 = a2*exp(-x2)*lj4[itype][jtype];
       force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]-
        	g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype];
-      if (eflag) one.eng_vdwl = factor_lj*r6inv*lj3[itype][jtype]-
+      eng += factor_lj*r6inv*lj3[itype][jtype]-
 	g6*((a2+1.0)*a2+0.5)*x2+t*lj4[itype][jtype];
     }
     else {						// cut
       force_lj = factor_lj*r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
-      if (eflag) one.eng_vdwl = factor_lj*(r6inv*(r6inv*lj3[itype][jtype]-
+      eng += factor_lj*(r6inv*(r6inv*lj3[itype][jtype]-
-	    lj4[itype][jtype])-offset[itype][jtype]);
+			       lj4[itype][jtype])-offset[itype][jtype]);
     }
-  } 
+  } else force_lj = 0.0;
-  else force_lj = 0.0;
 
-  one.fforce = (force_coul+force_lj)*r2inv;
+  fforce = (force_coul+force_lj)*r2inv;
+  return eng;
 }
-

src/USER-EWALDN/pair_lj_coul.h

@@ -35,7 +35,7 @@ class PairLJCoul : public Pair {
   
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
-  void single(int, int, int, int, double, double, double, int, One &);
+  double single(int, int, int, int, double, double, double, double &);
   void *extract(char *);
 
   void compute_inner();

src/bond.h

@@ -38,7 +38,7 @@ class Bond : protected Pointers {
   virtual double equilibrium_distance(int) = 0;
   virtual void write_restart(FILE *) = 0;
   virtual void read_restart(FILE *) = 0;
-  virtual void single(int, double, int, int, double &) = 0;
+  virtual double single(int, double, int, int) = 0;
   virtual double memory_usage();
 
  protected:

src/bond_hybrid.h

@@ -32,7 +32,7 @@ class BondHybrid : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  void single(int, double, int, int, double &);
+  double single(int, double, int, int);
   double memory_usage();
 
  private:

src/compute.cpp

@@ -14,6 +14,7 @@
 #include "mpi.h"
 #include "stdlib.h"
 #include "string.h"
+#include "ctype.h"
 #include "compute.h"
 #include "group.h"
 #include "domain.h"
@@ -32,12 +33,18 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   if (narg < 3) error->all("Illegal compute command");
 
   // compute ID, group, and style
+  // ID must be all alphanumeric chars or underscores
 
   int n = strlen(arg[0]) + 1;
   id = new char[n];
   strcpy(id,arg[0]);
 
+  for (int i = 0; i < n-1; i++)
+    if (!isalnum(id[i]) && id[i] != '_')
+      error->all("Compute ID must be alphanumeric or underscore characters");
+
   igroup = group->find(arg[1]);
+  if (igroup == -1) error->all("Could not find compute group ID");
   groupbit = group->bitmask[igroup];
 
   n = strlen(arg[2]) + 1;
@@ -53,10 +60,9 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   scalar_flag = vector_flag = peratom_flag = 0;
   tempflag = pressflag = peflag = 0;
   pressatomflag = peatomflag = 0;
+  id_pre = NULL;
   timeflag = 0;
   invoked = 0;
-  npre = 0;
-  id_pre = NULL;
   comm_forward = comm_reverse = 0;
 
   // set modify defaults
@@ -76,8 +82,6 @@ Compute::~Compute()
 {
   delete [] id;
   delete [] style;
-
-  for (int i = 0; i < npre; i++) delete [] id_pre[i];
   delete [] id_pre;
 
   memory->sfree(tlist);
@@ -117,7 +121,7 @@ void Compute::modify_params(int narg, char **arg)
    search from top downward, since list of times is in decreasing order
 ------------------------------------------------------------------------- */
 
-void Compute::add_step(int ntimestep)
+void Compute::addstep(int ntimestep)
 {
   // i = location in list to insert ntimestep
 
@@ -149,7 +153,7 @@ void Compute::add_step(int ntimestep)
    search from top downward, since list of times is in decreasing order
 ------------------------------------------------------------------------- */
 
-int Compute::match_step(int ntimestep)
+int Compute::matchstep(int ntimestep)
 {
   for (int i = ntime-1; i >= 0; i--) {
     if (ntimestep < tlist[i]) return 0;

src/compute.h

@@ -33,8 +33,10 @@ class Compute : protected Pointers {
   int size_vector;          // N = size of global vector
   int peratom_flag;         // 0/1 if compute_peratom() function exists
   int size_peratom;         // 0 = scalar_atom, N = size of vector_atom
+  int extscalar;            // 0/1 if scalar is intensive/extensive
+  int extvector;            // 0/1/-1 if vector is all int/ext/extlist
+  int *extlist;             // list of 0/1 int/ext for each vec component
 
-  int extensive;      // 0/1 if scalar,vector are intensive/extensive values
   int tempflag;       // 1 if Compute can be used as temperature
                       // must have both compute_scalar, compute_vector
   int pressflag;      // 1 if Compute can be used as pressure (uses virial)
@@ -43,17 +45,17 @@ class Compute : protected Pointers {
   int peflag;         // 1 if Compute calculates PE (uses Force energies)
   int peatomflag;     // 1 if Compute calculates per-atom PE
 
+  char *id_pre;       // ID of pre-compute the Compute may store
+
   int timeflag;       // 1 if Compute stores list of timesteps it's called on
   int ntime;          // # of entries in time list
   int maxtime;        // max # of entries time list can hold
   int *tlist;         // time list of steps the Compute is called on
-  int invoked;        // 1 if Compute was invoked (e.g. by a variable)
+
+  int invoked;        // set when Compute is invoked, to avoid re-invoking
 
   double dof;         // degrees-of-freedom for temperature
 
-  int npre;           // # of computes to compute before this one
-  char **id_pre;      // IDs of Computes to compute before this one
-
   int comm_forward;     // size of forward communication (0 if none)
   int comm_reverse;     // size of reverse communication (0 if none)
 
@@ -71,8 +73,8 @@ class Compute : protected Pointers {
   virtual int pack_reverse_comm(int, int, double *) {return 0;}
   virtual void unpack_reverse_comm(int, int *, double *) {}
 
-  void add_step(int);
+  void addstep(int);
-  int match_step(int);
+  int matchstep(int);
 
   virtual double memory_usage() {return 0.0;}
 

src/compute_attribute_atom.cpp

@@ -1,211 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 "string.h"
-#include "compute_attribute_atom.h"
-#include "atom.h"
-#include "domain.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-enum{X,Y,Z,XU,YU,ZU,VX,VY,VZ,FX,FY,FZ,XYZ,V,F};
-
-/* --------------------------------------------------------------------- */
-
-ComputeAttributeAtom::ComputeAttributeAtom(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg)
-{
-  if (narg != 4) error->all("Illegal compute ke/atom command");
-
-  peratom_flag = 1;
-  size_peratom = 0;
-
-  if (strcmp(arg[3],"x") == 0) which = X;
-  else if (strcmp(arg[3],"y") == 0) which = Y;
-  else if (strcmp(arg[3],"z") == 0) which = Z;
-  else if (strcmp(arg[3],"xu") == 0) which = XU;
-  else if (strcmp(arg[3],"yu") == 0) which = YU;
-  else if (strcmp(arg[3],"zu") == 0) which = ZU;
-  else if (strcmp(arg[3],"vx") == 0) which = VX;
-  else if (strcmp(arg[3],"vy") == 0) which = VY;
-  else if (strcmp(arg[3],"vz") == 0) which = VZ;
-  else if (strcmp(arg[3],"fx") == 0) which = FX;
-  else if (strcmp(arg[3],"fy") == 0) which = FY;
-  else if (strcmp(arg[3],"fz") == 0) which = FZ;
-
-  else if (strcmp(arg[3],"xyz") == 0) {
-    which = XYZ;
-    size_peratom = 3;
-  } else if (strcmp(arg[3],"v") == 0) {
-    which = V;
-    size_peratom = 3;
-  } else if (strcmp(arg[3],"f") == 0) {
-    which = F;
-    size_peratom = 3;
-  } else error->all("Illegal compute attribute/atom command");
-
-  nmax = 0;
-  s_attribute = NULL;
-  v_attribute = NULL;
-}
-
-/* ---------------------------------------------------------------------- */
-
-ComputeAttributeAtom::~ComputeAttributeAtom()
-{
-  memory->sfree(s_attribute);
-  memory->destroy_2d_double_array(v_attribute);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeAttributeAtom::compute_peratom()
-{
-  // grow attribute array if necessary
-
-  if (atom->nlocal > nmax) {
-    if (size_peratom == 0) {
-      memory->sfree(s_attribute);
-      nmax = atom->nmax;
-      s_attribute = (double *) 
-	memory->smalloc(nmax*sizeof(double),
-			"compute/attribute/atom:s_attribute");
-      scalar_atom = s_attribute;
-    } else {
-      memory->destroy_2d_double_array(v_attribute);
-      nmax = atom->nmax;
-      v_attribute =  
-	memory->create_2d_double_array(nmax,size_peratom,
-				       "compute/attribute/atom:v_attribute");
-      vector_atom = v_attribute;
-    }
-  }
-
-  // fill attribute vector with appropriate atom value
-  // or simply set pointer to exisitng atom vector
-
-  double xprd = domain->xprd;
-  double yprd = domain->yprd;
-  double zprd = domain->zprd;
-
-  double **x = atom->x;
-  double **v = atom->v;
-  double **f = atom->f;
-  int *mask = atom->mask;
-  int *image = atom->image;
-  int nlocal = atom->nlocal;
-
-  if (which == X) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = x[i][0];
-      else s_attribute[i] = 0.0;
-  } else if (which == Y) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = x[i][1];
-      else s_attribute[i] = 0.0;
-  } else if (which == Z) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = x[i][2];
-      else s_attribute[i] = 0.0;
-
-  } else if (which == XU) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit)
-	s_attribute[i] = x[i][0] + ((image[i] & 1023) - 512) * xprd;
-      else s_attribute[i] = 0.0;
-  } else if (which == YU) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit)
-	s_attribute[i] = x[i][1] + ((image[i] >> 10 & 1023) - 512) * yprd;
-      else s_attribute[i] = 0.0;
-  } else if (which == ZU) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit)
-	s_attribute[i] = x[i][2] + ((image[i] >> 20) - 512) * zprd;
-      else s_attribute[i] = 0.0;
-
-  } else if (which == VX) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = v[i][0];
-      else s_attribute[i] = 0.0;
-  } else if (which == VY) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = v[i][1];
-      else s_attribute[i] = 0.0;
-  } else if (which == VZ) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = v[i][2];
-      else s_attribute[i] = 0.0;
-
-  } else if (which == FX) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = f[i][0];
-      else s_attribute[i] = 0.0;
-  } else if (which == FY) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = f[i][1];
-      else s_attribute[i] = 0.0;
-  } else if (which == FZ) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_attribute[i] = f[i][2];
-      else s_attribute[i] = 0.0;
-
-  } else if (which == XYZ) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-	v_attribute[i][0] = x[i][0];
-	v_attribute[i][1] = x[i][1];
-	v_attribute[i][2] = x[i][2];
-      } else {
-	v_attribute[i][0] = 0.0;
-	v_attribute[i][1] = 0.0;
-	v_attribute[i][2] = 0.0;
-      }
-  } else if (which == V) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-	v_attribute[i][0] = v[i][0];
-	v_attribute[i][1] = v[i][1];
-	v_attribute[i][2] = v[i][2];
-      } else {
-	v_attribute[i][0] = 0.0;
-	v_attribute[i][1] = 0.0;
-	v_attribute[i][2] = 0.0;
-      }
-  } else if (which == F) {
-    for (int i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) {
-	v_attribute[i][0] = f[i][0];
-	v_attribute[i][1] = f[i][1];
-	v_attribute[i][2] = f[i][2];
-      } else {
-	v_attribute[i][0] = 0.0;
-	v_attribute[i][1] = 0.0;
-	v_attribute[i][2] = 0.0;
-      }
-  }
-}
-
-/* ----------------------------------------------------------------------
-   memory usage of local atom-based array
-------------------------------------------------------------------------- */
-
-double ComputeAttributeAtom::memory_usage()
-{
-  double bytes = 0.0;
-  if (size_peratom == 0) bytes = nmax * sizeof(double);
-  else bytes = size_peratom * nmax * sizeof(double);
-  return bytes;
-}

src/compute_attribute_atom.h

@@ -1,36 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifndef COMPUTE_ATTRIBUTE_ATOM_H
-#define COMPUTE_ATTRIBUTE_ATOM_H
-
-#include "compute.h"
-
-namespace LAMMPS_NS {
-
-class ComputeAttributeAtom : public Compute {
- public:
-  ComputeAttributeAtom(class LAMMPS *, int, char **);
-  ~ComputeAttributeAtom();
-  void init() {}
-  void compute_peratom();
-  double memory_usage();
-
- private:
-  int which,allocate,nmax;
-  double *s_attribute,**v_attribute;
-};
-
-}
-
-#endif

src/compute_centro_atom.cpp

@@ -27,6 +27,8 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_PERATOM 4
+
 /* ---------------------------------------------------------------------- */
 
 ComputeCentroAtom::ComputeCentroAtom(LAMMPS *lmp, int narg, char **arg) :
@@ -89,6 +91,8 @@ void ComputeCentroAtom::compute_peratom()
   int *ilist,*jlist,*numneigh,**firstneigh;
   double pairs[66];
 
+  invoked |= INVOKED_PERATOM;
+
   // grow centro array if necessary
 
   if (atom->nlocal > nmax) {

src/compute_coord_atom.cpp

@@ -27,6 +27,8 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_PERATOM 4
+
 /* ---------------------------------------------------------------------- */
 
 ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
@@ -88,6 +90,8 @@ void ComputeCoordAtom::compute_peratom()
   double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
+  invoked |= INVOKED_PERATOM;
+
   // grow coordination array if necessary
 
   if (atom->nlocal > nmax) {

src/compute_ke_atom.cpp

@@ -22,6 +22,8 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_PERATOM 4
+
 /* ---------------------------------------------------------------------- */
 
 ComputeKEAtom::ComputeKEAtom(LAMMPS *lmp, int narg, char **arg) :
@@ -58,6 +60,8 @@ void ComputeKEAtom::init()
 
 void ComputeKEAtom::compute_peratom()
 {
+  invoked |= INVOKED_PERATOM;
+
   // grow ke array if necessary
 
   if (atom->nlocal > nmax) {

src/compute_pe.cpp

@@ -28,6 +28,8 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
+
 /* ---------------------------------------------------------------------- */
 
 ComputePE::ComputePE(LAMMPS *lmp, int narg, char **arg) : 
@@ -62,7 +64,7 @@ ComputePE::ComputePE(LAMMPS *lmp, int narg, char **arg) :
   // settings
 
   scalar_flag = 1;
-  extensive = 1;
+  extscalar = 1;
   peflag = 1;
   timeflag = 1;
 }
@@ -71,9 +73,9 @@ ComputePE::ComputePE(LAMMPS *lmp, int narg, char **arg) :
 
 double ComputePE::compute_scalar()
 {
-  invoked = 1;
+  invoked |= INVOKED_SCALAR;
-  double one = 0.0;
 
+  double one = 0.0;
   if (pairflag && force->pair)
     one += force->pair->eng_vdwl + force->pair->eng_coul;
 

src/compute_pe_atom.cpp

@@ -26,6 +26,8 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_PERATOM 4
+
 /* ---------------------------------------------------------------------- */
 
 ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) : 
@@ -74,7 +76,7 @@ void ComputePEAtom::compute_peratom()
 {
   int i;
 
-  invoked = 1;
+  invoked |= INVOKED_PERATOM;
 
   // grow local energy array if necessary
 

src/compute_pressure.cpp

@@ -30,6 +30,9 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+
 /* ---------------------------------------------------------------------- */
 
 ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
@@ -40,21 +43,19 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 1;
   pressflag = 1;
   timeflag = 1;
 
   // store temperature ID used by pressure computation
   // insure it is valid for temperature computation
 
-  npre = 1;
-  id_pre = new char*[1];
   int n = strlen(arg[3]) + 1;
-  id_pre[0] = new char[n];
+  id_pre = new char[n];
-  strcpy(id_pre[0],arg[3]);
+  strcpy(id_pre,arg[3]);
-
-  int icompute = modify->find_compute(id_pre[0]);
 
+  int icompute = modify->find_compute(id_pre);
   if (icompute < 0) error->all("Could not find compute pressure temp ID");
   if (modify->compute[icompute]->tempflag == 0)
     error->all("Compute pressure temp ID does not compute temperature");
@@ -109,7 +110,7 @@ void ComputePressure::init()
 
   // set temperature used by pressure
 
-  int icompute = modify->find_compute(id_pre[0]);
+  int icompute = modify->find_compute(id_pre);
   if (icompute < 0) error->all("Could not find compute pressure temp ID");
   temperature = modify->compute[icompute];
 
@@ -153,16 +154,25 @@ void ComputePressure::init()
 
 /* ----------------------------------------------------------------------
    compute total pressure, averaged over Pxx, Pyy, Pzz
-   assume temperature has already been computed
 ------------------------------------------------------------------------- */
 
 double ComputePressure::compute_scalar()
 {
+  invoked |= INVOKED_SCALAR;
+
+  // invoke temperature it it hasn't been already
+
+  double t;
+  if (keflag) {
+    if (temperature->invoked & INVOKED_SCALAR) t = temperature->scalar;
+    else t = temperature->compute_scalar();
+  }
+
   if (dimension == 3) {
     inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
     virial_compute(3,3);
     if (keflag)
-      scalar = (temperature->dof * boltz * temperature->scalar + 
+      scalar = (temperature->dof * boltz * t + 
 		virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;
     else
       scalar = (virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p;
@@ -170,7 +180,7 @@ double ComputePressure::compute_scalar()
     inv_volume = 1.0 / (domain->xprd * domain->yprd);
     virial_compute(2,2);
     if (keflag)
-      scalar = (temperature->dof * boltz * temperature->scalar + 
+      scalar = (temperature->dof * boltz * t + 
 		virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;
     else
       scalar = (virial[0] + virial[1]) / 2.0 * inv_volume * nktv2p;
@@ -186,11 +196,21 @@ double ComputePressure::compute_scalar()
 
 void ComputePressure::compute_vector()
 {
+  invoked |= INVOKED_VECTOR;
+
+  // invoke temperature it it hasn't been already
+
+  double *ke_tensor;
+  if (keflag) {
+    if (!(temperature->invoked & INVOKED_VECTOR))
+      temperature->compute_vector();
+    ke_tensor = temperature->vector;
+  }
+
   if (dimension == 3) {
     inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
     virial_compute(6,3);
     if (keflag) {
-      double *ke_tensor = temperature->vector;
       for (int i = 0; i < 6; i++)
 	vector[i] = (ke_tensor[i] + virial[i]) * inv_volume * nktv2p;
     } else
@@ -200,7 +220,6 @@ void ComputePressure::compute_vector()
     inv_volume = 1.0 / (domain->xprd * domain->yprd);
     virial_compute(4,2);
     if (keflag) {
-      double *ke_tensor = temperature->vector;
       vector[0] = (ke_tensor[0] + virial[0]) * inv_volume * nktv2p;
       vector[1] = (ke_tensor[1] + virial[1]) * inv_volume * nktv2p;
       vector[3] = (ke_tensor[3] + virial[3]) * inv_volume * nktv2p;
@@ -219,7 +238,6 @@ void ComputePressure::virial_compute(int n, int ndiag)
   int i,j;
   double v[6],*vcomponent;
 
-  invoked = 1;
   for (i = 0; i < n; i++) v[i] = 0.0;
 
   // sum contributions to virial from forces and fixes

src/compute_rotate_dipole.cpp

@@ -20,6 +20,7 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
 #define INERTIA3D 0.4          // moments of inertia for sphere and disk
 #define INERTIA2D 0.5
 
@@ -34,7 +35,7 @@ ComputeRotateDipole::ComputeRotateDipole(LAMMPS *lmp, int narg, char **arg) :
     error->all("Compute rotate/dipole requires atom attributes dipole, omega");
 
   scalar_flag = 1;
-  extensive = 1;
+  extscalar = 1;
 
   inertia = NULL;
 }
@@ -68,6 +69,8 @@ void ComputeRotateDipole::init()
 
 double ComputeRotateDipole::compute_scalar()
 {
+  invoked |= INVOKED_SCALAR;
+
   double *dipole = atom->dipole;
   double **omega = atom->omega;
   int *type = atom->type;

src/compute_rotate_gran.cpp

@@ -21,6 +21,7 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
 #define INERTIA3D 0.4          // moments of inertia for sphere and disk
 #define INERTIA2D 0.5
 
@@ -36,7 +37,7 @@ ComputeRotateGran::ComputeRotateGran(LAMMPS *lmp, int narg, char **arg) :
 	       "radius, rmass, omega");
 
   scalar_flag = 1;
-  extensive = 1;
+  extscalar = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -51,6 +52,8 @@ void ComputeRotateGran::init()
 
 double ComputeRotateGran::compute_scalar()
 {
+  invoked |= INVOKED_SCALAR;
+
   double **omega = atom->omega;
   double *radius = atom->radius;
   double *rmass = atom->rmass;

src/compute_stress_atom.cpp

@@ -27,6 +27,8 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_PERATOM 4
+
 /* ---------------------------------------------------------------------- */
 
 ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) : 
@@ -78,7 +80,7 @@ void ComputeStressAtom::compute_peratom()
 {
   int i,j;
 
-  invoked = 1;
+  invoked |= INVOKED_PERATOM;
 
   // grow local stress array if necessary
 

src/compute_sum.cpp

@@ -12,16 +12,27 @@
 ------------------------------------------------------------------------- */
 
 #include "string.h"
+#include "stdlib.h"
 #include "compute_sum.h"
 #include "atom.h"
+#include "update.h"
 #include "modify.h"
+#include "fix.h"
 #include "force.h"
 #include "comm.h"
+#include "input.h"
+#include "variable.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+#define INVOKED_PERATOM 4
+
+enum{X,V,F,COMPUTE,FIX,VARIABLE};
+
 /* ---------------------------------------------------------------------- */
 
 ComputeSum::ComputeSum(LAMMPS *lmp, int narg, char **arg) :
@@ -29,74 +40,171 @@ ComputeSum::ComputeSum(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg < 4) error->all("Illegal compute sum command");
 
-  // store pre-compute IDs
+  // parse remaining values
 
-  npre = narg - 3;
+  which = new int[narg-3];
-  id_pre = new char*[npre];
+  argindex = new int[narg-3];
-  for (int i = 0; i < npre; i++) {
+  ids = new char*[narg-3];
-    int iarg = i + 3;
+  value2index = new int[narg-3];
-    int n = strlen(arg[iarg]) + 1;
+  nvalues = 0;
-    id_pre[i] = new char[n];
+
-    strcpy(id_pre[i],arg[iarg]);
+  int iarg = 3;
+  while (iarg < narg) {
+    ids[nvalues] = NULL;
+
+    if (strcmp(arg[iarg],"x") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"y") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"z") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 2;
+
+    } else if (strcmp(arg[iarg],"vx") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"vy") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"vz") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 2;
+
+    } else if (strcmp(arg[iarg],"fx") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"fy") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"fz") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 2;
+
+    } else if ((strncmp(arg[iarg],"c_",2) == 0) || 
+	       (strncmp(arg[iarg],"f_",2) == 0) || 
+	       (strncmp(arg[iarg],"v_",2) == 0)) {
+      if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
+      else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
+      else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
+
+      int n = strlen(arg[iarg]);
+      char *suffix = new char[n];
+      strcpy(suffix,&arg[iarg][2]);
+
+      char *ptr = strchr(suffix,'[');
+      if (ptr) {
+	if (suffix[strlen(suffix)-1] != ']')
+	  error->all("Illegal compute sum command");
+	argindex[nvalues] = atoi(ptr+1);
+	*ptr = '\0';
+      } else argindex[nvalues] = 0;
+
+      n = strlen(suffix) + 1;
+      ids[nvalues] = new char[n];
+      strcpy(ids[nvalues],suffix);
+      nvalues++;
+      delete [] suffix;
+
+    } else error->all("Illegal compute sum command");
+
+    iarg++;
   }
 
-  compute = new Compute*[npre];
+  // setup and error check
 
-  // all sub-computes must be peratom
+  for (int i = 0; i < nvalues; i++) {
-  // check consistency of sub-computes for scalar & vector output
+    if (which[i] == COMPUTE) {
-
+      int icompute = modify->find_compute(ids[i]);
-  int icompute;
+      if (icompute < 0)
-  for (int i = 0; i < npre; i++) {
+	error->all("Compute ID for compute sum does not exist");
-    icompute = modify->find_compute(id_pre[i]);
+      if (modify->compute[icompute]->peratom_flag == 0)
-    if (icompute < 0)
+	error->all("Compute sum compute does not calculate per-atom values");
-      error->all("Could not find compute sum/atom pre-compute ID");
+      if (argindex[i] == 0 && modify->compute[icompute]->size_peratom != 0)
-    if (modify->compute[icompute]->peratom_flag == 0)
+	error->all("Compute sum compute does not calculate a per-atom scalar");
-      error->all("Compute sum compute is not a per-atom compute");
+      if (argindex[i] && modify->compute[icompute]->size_peratom == 0)
+	error->all("Compute sum compute does not calculate a per-atom vector");
+    } else if (which[i] == FIX) {
+      int ifix = modify->find_fix(ids[i]);
+      if (ifix < 0)
+	error->all("Fix ID for compute sum does not exist");
+      if (modify->fix[ifix]->peratom_flag == 0)
+	error->all("Compute sum fix does not calculate per-atom values");
+      if (argindex[i] == 0 && modify->fix[ifix]->size_peratom != 0)
+	error->all("Compute sum fix does not calculate a per-atom scalar");
+      if (argindex[i] && modify->fix[ifix]->size_peratom == 0)
+	error->all("Compute sum fix does not calculate a per-atom vector");
+    } else if (which[i] == VARIABLE) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0)
+	error->all("Variable name for compute sum does not exist");
+      if (input->variable->atomstyle(ivariable) == 0)
+	error->all("Compute sum variable is not atom-style variable");
+    }
   }
 
-  peratom_flag = 0;
+  // this compute produces either a scalar or vector
-  extensive = 0;
 
-  icompute = modify->find_compute(id_pre[0]);
+  if (nvalues == 1) {
-  int size = modify->compute[icompute]->size_peratom;
-  if (size == 0) {
     scalar_flag = 1;
+    extscalar = 1;
     vector = NULL;
     onevec = NULL;
   } else {
     vector_flag = 1;
-    size_vector = size;
+    size_vector = nvalues;
+    extvector = 1;
     vector = new double[size_vector];
     onevec = new double[size_vector];
   }
 
-  for (int i = 1; i < npre; i++) {
+  maxatom = 0;
-    icompute = modify->find_compute(id_pre[i]);
+  varatom = NULL;
-    if (modify->compute[icompute]->size_peratom != size)
-      error->all("Inconsistent sizes of compute sum compute quantities");
-  }
 }
 
 /* ---------------------------------------------------------------------- */
 
 ComputeSum::~ComputeSum()
 {
-  delete [] compute;
+  delete [] which;
+  delete [] argindex;
+  for (int m = 0; m < nvalues; m++) delete [] ids[m];
+  delete [] ids;
+  delete [] value2index;
+
   delete [] vector;
   delete [] onevec;
+
+  memory->sfree(varatom);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void ComputeSum::init()
 {
-  // set ptrs to Computes used as pre-computes by this compute
+  // set indices and check validity of all computes,fixes,variables
 
-  for (int i = 0; i < npre; i++) {
+  for (int m = 0; m < nvalues; m++) {
-    int icompute = modify->find_compute(id_pre[i]);
+    if (which[m] == COMPUTE) {
-    if (icompute < 0)
+      int icompute = modify->find_compute(ids[m]);
-      error->all("Could not find compute sum/atom pre-compute ID");
+      if (icompute < 0)
-    compute[i] = modify->compute[icompute];
+	error->all("Compute ID for compute sum does not exist");
+      value2index[m] = icompute;
+      
+    } else if (which[m] == FIX) {
+      int ifix = modify->find_fix(ids[m]);
+      if (ifix < 0) 
+	error->all("Fix ID for compute sum does not exist");
+      value2index[m] = ifix;
+
+    } else if (which[m] == VARIABLE) {
+      int ivariable = input->variable->find(ids[m]);
+      if (ivariable < 0) 
+	error->all("Variable name for compute sum does not exist");
+      value2index[m] = ivariable;
+
+    } else value2index[m] = -1;
   }
 }
 
@@ -104,31 +212,10 @@ void ComputeSum::init()
 
 double ComputeSum::compute_scalar()
 {
-  int i;
+  invoked |= INVOKED_SCALAR;
-
-  // invoke all the pre-computes
-  // this is the only compute that does this
-  // done b/c pre-computes are per-atom and this compute is not
-
-  for (int icompute = 0; icompute < npre; icompute++)
-    compute[icompute]->compute_peratom();
-
-  // compute scalar quantity by summing over atom scalars
-  // only include atoms in group
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  double one = 0.0;
-
-  for (int icompute = 0; icompute < npre; icompute++) {
-    double *scalar_atom = compute[icompute]->scalar_atom;
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) one += scalar_atom[i];
-  }
 
+  double one = compute_one(0);
   MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
-
   return scalar;
 }
 
@@ -136,30 +223,100 @@ double ComputeSum::compute_scalar()
 
 void ComputeSum::compute_vector()
 {
-  int i,j;
+  invoked |= INVOKED_VECTOR;
 
-  // invoke all the pre-computes
+  for (int m = 0; m < nvalues; m++) onevec[m] = compute_one(m);
-  // this is the only compute that does this
+  MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
-  // done b/c pre-computes are per-atom and this compute is not
+}
 
-  for (int icompute = 0; icompute < npre; icompute++)
+/* ---------------------------------------------------------------------- */
-    compute[icompute]->compute_peratom();
 
-  // compute vector quantity by summing over atom vectors
+double ComputeSum::compute_one(int m)
+{
+  int i;
+
+  // invoke the appropriate attribute,compute,fix,variable
+  // compute scalar quantity by summing over atom scalars
   // only include atoms in group
 
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  for (j = 0; j < size_vector; j++) onevec[j] = 0.0;
+  int n = value2index[m];
+  int j = argindex[m];
+  double one = 0.0;
 
-  for (int icompute = 0; icompute < npre; icompute++) {
+  if (which[m] == X) {
-    double **vector_atom = compute[icompute]->vector_atom;
+    double **x = atom->x;
     for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) 
+      if (mask[i] & groupbit) one += x[i][j];
-	for (j = 0; j < size_vector; j++)
+  } else if (which[m] == V) {
-	  onevec[j] += vector_atom[i][j];
+    double **v = atom->v;
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) one += v[i][j];
+  } else if (which[m] == F) {
+    double **f = atom->f;
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) one += f[i][j];
+    
+  // invoke compute if not previously invoked
+
+  } else if (which[m] == COMPUTE) {
+    if (!(modify->compute[n]->invoked & INVOKED_PERATOM))
+      modify->compute[n]->compute_peratom();
+
+    if (j == 0) {
+      double *compute_scalar = modify->compute[n]->scalar_atom;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) one += compute_scalar[i];
+    } else {
+      double **compute_vector = modify->compute[n]->vector_atom;
+      int jm1 = j - 1;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) one += compute_vector[i][jm1];
+    }
+
+  // access fix fields, check if frequency is a match
+
+  } else if (which[m] == FIX) {
+    if (update->ntimestep % modify->fix[n]->peratom_freq)
+      error->all("Fix used in compute sum not computed at compatible time");
+
+    if (j == 0) {
+      double *fix_scalar = modify->fix[n]->scalar_atom;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) one += fix_scalar[i];
+    } else {
+      double **fix_vector = modify->fix[n]->vector_atom;
+      int jm1 = j - 1;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) one += fix_vector[i][jm1];
+    }
+    
+  // evaluate atom-style variable
+
+  } else if (which[m] == VARIABLE) {
+    if (nlocal > maxatom) {
+      maxatom = atom->nmax;
+      memory->sfree(varatom);
+      varatom =	(double *) 
+	memory->smalloc(maxatom*sizeof(double),"compute/sum:varatom");
+    }
+
+    input->variable->compute_atom(n,igroup,varatom,1,0);
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) one += varatom[i];
   }
 
-  MPI_Allreduce(onevec,vector,size_vector,MPI_DOUBLE,MPI_SUM,world);
+  return one;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of varatom
+------------------------------------------------------------------------- */
+
+double ComputeSum::memory_usage()
+{
+  double bytes = maxatom * sizeof(double);
+  return bytes;
 }

src/compute_sum.h

@@ -25,10 +25,18 @@ class ComputeSum : public Compute {
   void init();
   double compute_scalar();
   void compute_vector();
+  double memory_usage();
 
  private:
+  int nvalues;
+  int *which,*argindex,*value2index;
+  char **ids;
   double *onevec;
-  class Compute **compute;
+
+  int maxatom;
+  double *varatom;
+
+  double compute_one(int);
 };
 
 }

src/compute_sum_atom.cpp

@@ -1,166 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 "string.h"
-#include "compute_sum_atom.h"
-#include "atom.h"
-#include "modify.h"
-#include "force.h"
-#include "comm.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-/* ---------------------------------------------------------------------- */
-
-ComputeSumAtom::ComputeSumAtom(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg)
-{
-  if (narg < 5) error->all("Illegal compute sum/atom command");
-
-  // store pre-compute IDs
-
-  npre = narg - 3;
-  id_pre = new char*[npre];
-  for (int i = 0; i < npre; i++) {
-    int iarg = i + 3;
-    int n = strlen(arg[iarg]) + 1;
-    id_pre[i] = new char[n];
-    strcpy(id_pre[i],arg[iarg]);
-  }
-
-  compute = new Compute*[npre];
-
-  // all sub-computes must be peratom
-  // check consistency of sub-computes for scalar & vector output
-
-  int icompute;
-  for (int i = 0; i < npre; i++) {
-    icompute = modify->find_compute(id_pre[i]);
-    if (icompute < 0)
-      error->all("Could not find compute sum/atom pre-compute ID");
-    if (modify->compute[icompute]->peratom_flag == 0)
-      error->all("Compute sum/atom compute is not a per-atom compute");
-  }
-
-  peratom_flag = 1;
-  icompute = modify->find_compute(id_pre[0]);
-  size_peratom = modify->compute[icompute]->size_peratom;
-
-  for (int i = 1; i < npre; i++) {
-    icompute = modify->find_compute(id_pre[i]);
-    if (modify->compute[icompute]->size_peratom != size_peratom)
-      error->all("Inconsistent sizes of compute sum/atom compute quantities");
-  }
-
-  // values generated by this compute
-
-  nmax = 0;
-  s_value = NULL;
-  v_value = NULL;
-}
-
-/* ---------------------------------------------------------------------- */
-
-ComputeSumAtom::~ComputeSumAtom()
-{
-  delete [] compute;
-  memory->sfree(s_value);
-  memory->destroy_2d_double_array(v_value);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeSumAtom::init()
-{
-  // set ptrs to Computes used as pre-computes by this compute
-
-  for (int i = 0; i < npre; i++) {
-    int icompute = modify->find_compute(id_pre[i]);
-    if (icompute < 0)
-      error->all("Could not find compute sum/atom pre-compute ID");
-    compute[i] = modify->compute[icompute];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeSumAtom::compute_peratom()
-{
-  int i,j,m;
-
-  // grow sum array if necessary
-
-  if (atom->nlocal > nmax) {
-    nmax = atom->nmax;
-    if (size_peratom == 0) {
-      memory->sfree(s_value);
-      s_value = (double *) 
-	memory->smalloc(nmax*sizeof(double),"compute/sum/atom:s_value");
-      scalar_atom = s_value;
-    } else {
-      memory->destroy_2d_double_array(v_value);
-      v_value = memory->create_2d_double_array(nmax,size_peratom,
-					       "compute/sum/atom:v_value");
-      vector_atom = v_value;
-    }
-  }
-
-  // sum over pre-computes
-  // pre-computes of the pre-computes are not invoked
-
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  if (size_peratom == 0) {
-    double *scalar = compute[0]->scalar_atom;
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) s_value[i] = scalar[i];
-      else s_value[i] = 0.0;
-    
-    for (m = 1; m < npre; m++) {
-      scalar = compute[m]->scalar_atom;
-      for (i = 0; i < nlocal; i++)
-	if (mask[i] & groupbit) s_value[i] += scalar[i];
-    }
-
-  } else {
-    double **vector = compute[0]->vector_atom;
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit) 
-	for (j = 0; j < size_peratom; j++)
-	  v_value[i][j] = vector[i][j];
-      else 
-	for (j = 0; j < size_peratom; j++)
-	  v_value[i][j] = 0.0;
-
-    for (m = 1; m < npre; m++) {
-      vector = compute[m]->vector_atom;
-      for (j = 0; j < size_peratom; j++)
-	if (mask[i] & groupbit) v_value[i][j] += vector[i][j];
-    }
-  }
-}
-
-/* ----------------------------------------------------------------------
-   memory usage of local atom-based array
-------------------------------------------------------------------------- */
-
-double ComputeSumAtom::memory_usage()
-{
-  double bytes = 0.0;
-  if (size_peratom == 0) bytes = nmax * sizeof(double);
-  else bytes = nmax*size_peratom * sizeof(double);
-  return bytes;
-}

src/compute_sum_atom.h

@@ -1,37 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifndef COMPUTE_SUM_ATOM_H
-#define COMPUTE_SUM_ATOM_H
-
-#include "compute.h"
-
-namespace LAMMPS_NS {
-
-class ComputeSumAtom : public Compute {
- public:
-  ComputeSumAtom(class LAMMPS *, int, char **);
-  ~ComputeSumAtom();
-  void init();
-  void compute_peratom();
-  double memory_usage();
-
- private:
-  int nmax;
-  class Compute **compute;
-  double *s_value,**v_value;
-};
-
-}
-
-#endif

src/compute_temp.cpp

@@ -23,6 +23,9 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+
 /* ---------------------------------------------------------------------- */
 
 ComputeTemp::ComputeTemp(LAMMPS *lmp, int narg, char **arg) : 
@@ -32,7 +35,8 @@ ComputeTemp::ComputeTemp(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 1;
   tempflag = 1;
 
   vector = new double[6];
@@ -70,6 +74,8 @@ void ComputeTemp::recount()
 
 double ComputeTemp::compute_scalar()
 {
+  invoked |= INVOKED_SCALAR;
+
   double **v = atom->v;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
@@ -102,6 +108,8 @@ void ComputeTemp::compute_vector()
 {
   int i;
 
+  invoked |= INVOKED_VECTOR;
+
   double **v = atom->v;
   double *mass = atom->mass;
   double *rmass = atom->rmass;

src/compute_temp_deform.cpp

@@ -30,6 +30,8 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
 enum{NO_REMAP,X_REMAP,V_REMAP};                   // same as fix_deform.cpp
 
 /* ---------------------------------------------------------------------- */
@@ -41,7 +43,8 @@ ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 1;
   tempflag = 1;
 
   vector = new double[6];
@@ -95,6 +98,8 @@ double ComputeTempDeform::compute_scalar()
 {
   double lamda[3],vstream[3],vthermal[3];
 
+  invoked |= INVOKED_SCALAR;
+
   double **x = atom->x;
   double **v = atom->v;
   double *mass = atom->mass;
@@ -153,6 +158,8 @@ void ComputeTempDeform::compute_vector()
 {
   double lamda[3],vstream[3],vthermal[3];
 
+  invoked |= INVOKED_VECTOR;
+
   double **x = atom->x;
   double **v = atom->v;
   double *mass = atom->mass;

src/compute_temp_partial.cpp

@@ -23,6 +23,9 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+
 /* ---------------------------------------------------------------------- */
 
 ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) : 
@@ -36,7 +39,8 @@ ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 1;
   tempflag = 1;
 
   vector = new double[6];
@@ -74,6 +78,8 @@ void ComputeTempPartial::recount()
 
 double ComputeTempPartial::compute_scalar()
 {
+  invoked |= INVOKED_SCALAR;
+
   double **v = atom->v;
   double *mass = atom->mass;
   double *rmass = atom->rmass;
@@ -107,6 +113,8 @@ void ComputeTempPartial::compute_vector()
 {
   int i;
 
+  invoked |= INVOKED_VECTOR;
+
   double **v = atom->v;
   double *mass = atom->mass;
   double *rmass = atom->rmass;

src/compute_temp_ramp.cpp

@@ -29,6 +29,9 @@ using namespace LAMMPS_NS;
 #define MIN(A,B) ((A) < (B)) ? (A) : (B)
 #define MAX(A,B) ((A) > (B)) ? (A) : (B)
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+
 /* ---------------------------------------------------------------------- */
 
 ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) : 
@@ -101,7 +104,8 @@ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 1;
   tempflag = 1;
 
   vector = new double[6];
@@ -141,6 +145,8 @@ double ComputeTempRamp::compute_scalar()
 {
   double fraction,vramp,vthermal[3];
 
+  invoked |= INVOKED_SCALAR;
+
   double **x = atom->x;
   double **v = atom->v;
   double *mass = atom->mass;
@@ -181,6 +187,8 @@ void ComputeTempRamp::compute_vector()
   int i;
   double fraction,vramp,vthermal[3];
 
+  invoked |= INVOKED_VECTOR;
+
   double **x = atom->x;
   double **v = atom->v;
   double *mass = atom->mass;

src/compute_temp_region.cpp

@@ -22,6 +22,9 @@
 
 using namespace LAMMPS_NS;
 
+#define INVOKED_SCALAR 1
+#define INVOKED_VECTOR 2
+
 /* ---------------------------------------------------------------------- */
 
 ComputeTempRegion::ComputeTempRegion(LAMMPS *lmp, int narg, char **arg) : 
@@ -34,7 +37,8 @@ ComputeTempRegion::ComputeTempRegion(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = vector_flag = 1;
   size_vector = 6;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 1;
   tempflag = 1;
 
   vector = new double[6];
@@ -58,6 +62,8 @@ void ComputeTempRegion::init()
 
 double ComputeTempRegion::compute_scalar()
 {
+  invoked |= INVOKED_SCALAR;
+
   double **x = atom->x;
   double **v = atom->v;
   double *mass = atom->mass;
@@ -102,6 +108,8 @@ void ComputeTempRegion::compute_vector()
 {
   int i;
 
+  invoked |= INVOKED_VECTOR;
+
   double **x = atom->x;
   double **v = atom->v;
   double *mass = atom->mass;

src/compute_variable.cpp

@@ -1,65 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 "mpi.h"
-#include "string.h"
-#include "stdlib.h"
-#include "compute_variable.h"
-#include "input.h"
-#include "variable.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-/* ---------------------------------------------------------------------- */
-
-ComputeVariable::ComputeVariable(LAMMPS *lmp, int narg, char **arg) : 
-  Compute(lmp, narg, arg)
-{
-  if (narg != 4) error->all("Illegal compute variable command");
-
-  // store variable name
-
-  int n = strlen(arg[3]) + 1;
-  varname = new char[n];
-  strcpy(varname,arg[3]);
-
-  scalar_flag = 1;
-  extensive = 0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-ComputeVariable::~ComputeVariable()
-{
-  delete [] varname;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeVariable::init()
-{
-  // check if variable exists
-
-  int ivariable = input->variable->find(varname);
-  if (ivariable < 0)
-    error->all("Could not find compute variable name");
-}
-
-/* ---------------------------------------------------------------------- */
-
-double ComputeVariable::compute_scalar()
-{
-  scalar = atof(input->variable->retrieve(varname));
-  return scalar;
-}

src/compute_variable.h

@@ -1,34 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifndef COMPUTE_VARIABLE_H
-#define COMPUTE_VARIABLE_H
-
-#include "compute.h"
-
-namespace LAMMPS_NS {
-
-class ComputeVariable : public Compute {
- public:
-  ComputeVariable(class LAMMPS *, int, char **);
-  ~ComputeVariable();
-  void init();
-  double compute_scalar();
-
- private:
-  char *varname;
-};
-
-}
-
-#endif

src/compute_variable_atom.cpp

@@ -1,94 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 "string.h"
-#include "compute_variable_atom.h"
-#include "atom.h"
-#include "input.h"
-#include "variable.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-/* ---------------------------------------------------------------------- */
-
-ComputeVariableAtom::ComputeVariableAtom(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg)
-{
-  if (narg != 4) error->all("Illegal compute variable/atom command");
-
-  // store variable name
-
-  int n = strlen(arg[3]) + 1;
-  varname = new char[n];
-  strcpy(varname,arg[3]);
-
-  peratom_flag = 1;
-  size_peratom = 0;
-
-  nmax = 0;
-  result = NULL;
-}
-
-/* ---------------------------------------------------------------------- */
-
-ComputeVariableAtom::~ComputeVariableAtom()
-{
-  delete [] varname;
-  memory->sfree(result);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeVariableAtom::init()
-{
-  // set ivariable used by this compute and check if it exists
-
-  ivariable = input->variable->find(varname);
-  if (ivariable < 0)
-    error->all("Could not find compute variable name");
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeVariableAtom::compute_peratom()
-{
-  // grow result array if necessary
-
-  if (atom->nlocal > nmax) {
-    memory->sfree(result);
-    nmax = atom->nmax;
-    result = (double *) memory->smalloc(nmax*sizeof(double),
-					"compute/variable/atom:result");
-    scalar_atom = result;
-  }
-
-  // parse variable once to create parse tree
-  // evaluate tree for all atoms, will be zero for atoms not in group
-  // free parse tree memory stored by Variable
-
-  input->variable->build_parse_tree(ivariable);
-  input->variable->evaluate_parse_tree(igroup,result);
-  input->variable->free_parse_tree();
-}
-
-/* ----------------------------------------------------------------------
-   memory usage of local atom-based array
-------------------------------------------------------------------------- */
-
-double ComputeVariableAtom::memory_usage()
-{
-  double bytes = nmax * sizeof(double);
-  return bytes;
-}

src/compute_variable_atom.h

@@ -1,37 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifndef COMPUTE_VARIABLE_ATOM_H
-#define COMPUTE_VARIABLE_ATOM_H
-
-#include "compute.h"
-
-namespace LAMMPS_NS {
-
-class ComputeVariableAtom : public Compute {
- public:
-  ComputeVariableAtom(class LAMMPS *, int, char **);
-  ~ComputeVariableAtom();
-  void init();
-  void compute_peratom();
-  double memory_usage();
-
- private:
-  int nmax,ivariable;
-  char *varname;
-  double *result;
-};
-
-}
-
-#endif

src/dump_custom.cpp

@@ -40,6 +40,8 @@ enum{TAG,MOL,TYPE,X,Y,Z,XS,YS,ZS,XU,YU,ZU,IX,IY,IZ,
 enum{LT,LE,GT,GE,EQ,NEQ};
 enum{INT,DOUBLE};
 
+#define INVOKED_PERATOM 4          // same as in computes
+
 /* ---------------------------------------------------------------------- */
 
 DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
@@ -247,7 +249,7 @@ int DumpCustom::count()
 {
   int i;
 
-  // grow choose and variable storage arrays if needed
+  // grow choose and variable vbuf arrays if needed
 
   int nlocal = atom->nlocal;
   if (nlocal > maxlocal) {
@@ -267,22 +269,19 @@ int DumpCustom::count()
   }
 
   // invoke Computes for per-atom dump quantities
+  // only if not already invoked
 
   if (ncompute)
-    for (i = 0; i < ncompute; i++) compute[i]->compute_peratom();
+    for (i = 0; i < ncompute; i++)
+      if (!(compute[i]->invoked & INVOKED_PERATOM))
+	compute[i]->compute_peratom();
 
-  // invoke Variables for per-atom dump quantities
+  // evaluate atom-style Variables for per-atom dump quantities
-  // parse variable once to create parse tree
-  // evaluate tree for all atoms, will be zero for atoms not in group
-  // free parse tree memory stored by Variable
 
   if (nvariable)
-    for (i = 0; i < nvariable; i++) {
+    for (i = 0; i < nvariable; i++)
-      input->variable->build_parse_tree(variable[i]);
+      input->variable->compute_atom(variable[i],igroup,vbuf[i],1,0);
-      input->variable->evaluate_parse_tree(igroup,vbuf[i]);
+  
-      input->variable->free_parse_tree();
-    }
-
   // choose all local atoms for output
 
   for (i = 0; i < nlocal; i++) choose[i] = 1;
@@ -722,7 +721,6 @@ void DumpCustom::parse_fields(int narg, char **arg)
 
     // compute value = c_ID
     // if no trailing [], then arg is set to 0, else arg is between []
-    // if Compute has pre-computes, first add them to list
 
     } else if (strncmp(arg[iarg],"c_",2) == 0) {
       pack_choice[i] = &DumpCustom::pack_compute;
@@ -752,9 +750,6 @@ void DumpCustom::parse_fields(int narg, char **arg)
 	  argindex[i] > modify->compute[n]->size_peratom)
 	error->all("Dump custom compute ID vector is not large enough");
 
-      if (modify->compute[n]->npre)
-	for (int ic = 0; ic < modify->compute[n]->npre; ic++)
-	  int tmp = add_compute(modify->compute[n]->id_pre[ic]);
       field2index[i] = add_compute(suffix);
       delete [] suffix;
       
@@ -806,8 +801,8 @@ void DumpCustom::parse_fields(int narg, char **arg)
 
       n = input->variable->find(suffix);
       if (n < 0) error->all("Could not find dump custom variable name");
-      if (input->variable->peratom(n) == 0)
+      if (input->variable->atomstyle(n) == 0)
-	error->all("Dump custom variable does not compute peratom info");
+	error->all("Dump custom variable is not atom-style variable");
 
       field2index[i] = add_variable(suffix);
       delete [] suffix;
@@ -983,7 +978,6 @@ int DumpCustom::modify_param(int narg, char **arg)
     // compute value = c_ID
     // if no trailing [], then arg is set to 0, else arg is between []
     // must grow field2index and argindex arrays, since access is beyond nfield
-    // if Compute has pre-computes, first add them to list
 
     else if (strncmp(arg[1],"c_",2) == 0) {
       thresh_array[nthresh] = COMPUTE;
@@ -1020,9 +1014,6 @@ int DumpCustom::modify_param(int narg, char **arg)
 	  argindex[nfield+nthresh] > modify->compute[n]->size_peratom)
 	error->all("Dump custom compute ID vector is not large enough");
 
-      if (modify->compute[n]->npre)
-	for (int ic = 0; ic < modify->compute[n]->npre; ic++)
-	  int tmp = add_compute(modify->compute[n]->id_pre[ic]);
       field2index[nfield+nthresh] = add_compute(suffix);
       delete [] suffix;
 
@@ -1087,8 +1078,8 @@ int DumpCustom::modify_param(int narg, char **arg)
       
       n = input->variable->find(suffix);
       if (n < 0) error->all("Could not find dump custom variable name");
-      if (input->variable->peratom(n) == 0)
+      if (input->variable->atomstyle(n) == 0)
-	error->all("Dump custom variable does not compute peratom info");
+	error->all("Dump custom variable is not atom-style variable");
 
       field2index[nfield+nthresh] = add_variable(suffix);
       delete [] suffix;
@@ -1129,11 +1120,80 @@ double DumpCustom::memory_usage()
   return bytes;
 }
 
-// ----------------------------------------------------------------------
+/* ----------------------------------------------------------------------
-// one method for every keyword dump custom can output
+   extraction of Compute, Fix, Variable results
-// the atom quantity is packed into buf starting at n with stride size_one
+------------------------------------------------------------------------- */
-// customize by adding a method
+
-// ----------------------------------------------------------------------
+void DumpCustom::pack_compute(int n)
+{
+  double *vector = compute[field2index[n]]->scalar_atom;
+  double **array = compute[field2index[n]]->vector_atom;
+  int index = argindex[n];
+
+  int nlocal = atom->nlocal;
+
+  if (index == 0) {
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) {
+	buf[n] = vector[i];
+	n += size_one;
+      }
+  } else {
+    index--;
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) {
+	buf[n] = array[i][index];
+	n += size_one;
+      }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_fix(int n)
+{
+  double *vector = fix[field2index[n]]->scalar_atom;
+  double **array = fix[field2index[n]]->vector_atom;
+  int index = argindex[n];
+
+  int nlocal = atom->nlocal;
+
+  if (index == 0) {
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) {
+	buf[n] = vector[i];
+	n += size_one;
+      }
+  } else {
+    index--;
+    for (int i = 0; i < nlocal; i++)
+      if (choose[i]) {
+	buf[n] = array[i][index];
+	n += size_one;
+      }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DumpCustom::pack_variable(int n)
+{
+  double *vector = vbuf[field2index[n]];
+
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++)
+    if (choose[i]) {
+      buf[n] = vector[i];
+      n += size_one;
+    }
+}
+
+/* ----------------------------------------------------------------------
+   one method for every keyword dump custom can output
+   the atom quantity is packed into buf starting at n with stride size_one
+   customize a new keyword by adding a method
+------------------------------------------------------------------------- */
 
 /* ---------------------------------------------------------------------- */
 
@@ -1597,70 +1657,3 @@ void DumpCustom::pack_tqz(int n)
       n += size_one;
     }
 }
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_compute(int n)
-{
-  double *vector = compute[field2index[n]]->scalar_atom;
-  double **array = compute[field2index[n]]->vector_atom;
-  int index = argindex[n];
-
-  int nlocal = atom->nlocal;
-
-  if (index == 0) {
-    for (int i = 0; i < nlocal; i++)
-      if (choose[i]) {
-	buf[n] = vector[i];
-	n += size_one;
-      }
-  } else {
-    index--;
-    for (int i = 0; i < nlocal; i++)
-      if (choose[i]) {
-	buf[n] = array[i][index];
-	n += size_one;
-      }
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_fix(int n)
-{
-  double *vector = fix[field2index[n]]->scalar_atom;
-  double **array = fix[field2index[n]]->vector_atom;
-  int index = argindex[n];
-
-  int nlocal = atom->nlocal;
-
-  if (index == 0) {
-    for (int i = 0; i < nlocal; i++)
-      if (choose[i]) {
-	buf[n] = vector[i];
-	n += size_one;
-      }
-  } else {
-    index--;
-    for (int i = 0; i < nlocal; i++)
-      if (choose[i]) {
-	buf[n] = array[i][index];
-	n += size_one;
-      }
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void DumpCustom::pack_variable(int n)
-{
-  double *vector = vbuf[field2index[n]];
-
-  int nlocal = atom->nlocal;
-
-  for (int i = 0; i < nlocal; i++)
-    if (choose[i]) {
-      buf[n] = vector[i];
-      n += size_one;
-    }
-}

src/fix.cpp

@@ -12,6 +12,7 @@
 ------------------------------------------------------------------------- */
 
 #include "string.h"
+#include "ctype.h"
 #include "fix.h"
 #include "group.h"
 #include "error.h"
@@ -22,11 +23,19 @@ using namespace LAMMPS_NS;
 
 Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 {
+  // fix ID, group, and style
+  // ID must be all alphanumeric chars or underscores
+
   int n = strlen(arg[0]) + 1;
   id = new char[n];
   strcpy(id,arg[0]);
 
+  for (int i = 0; i < n-1; i++)
+    if (!isalnum(id[i]) && id[i] != '_')
+      error->all("Fix ID must be alphanumeric or underscore characters");
+
   igroup = group->find(arg[1]);
+  if (igroup == -1) error->all("Could not find fix group ID");
   groupbit = group->bitmask[igroup];
 
   n = strlen(arg[2]) + 1;

src/fix.h

@@ -38,7 +38,9 @@ class Fix : protected Pointers {
   int vector_flag;               // 0/1 if compute_vector() function exists
   int size_vector;               // N = size of global vector
   int scalar_vector_freq;        // frequency compute s/v data is available at
-  int extensive;                 // 0/1 if s/v are intensive/extensive values
+  int extscalar;                 // 0/1 if scalar is intensive/extensive
+  int extvector;                 // 0/1/-1 if vector is all int/ext/extlist
+  int *extlist;                  // list of 0/1 int/ext for each vec component
 
   int peratom_flag;              // 0/1 if per-atom data is stored
   int size_peratom;              // 0 = scalar_atom, N = size of vector_atom

src/fix_ave_atom.cpp

@@ -18,25 +18,99 @@
 #include "update.h"
 #include "modify.h"
 #include "compute.h"
+#include "input.h"
+#include "variable.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
+enum{X,V,F,COMPUTE,FIX,VARIABLE};
+
+#define INVOKED_PERATOM 4      // same as in computes
+
 /* ---------------------------------------------------------------------- */
 
 FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg != 7) error->all("Illegal fix ave/atom command");
+  if (narg < 7) error->all("Illegal fix ave/atom command");
 
   nevery = atoi(arg[3]);
   nrepeat = atoi(arg[4]);
   peratom_freq = atoi(arg[5]);
 
-  int n = strlen(arg[6]) + 1;
+  // parse remaining values
-  id_compute = new char[n];
+
-  strcpy(id_compute,arg[6]);
+  which = new int[narg-6];
+  argindex = new int[narg-6];
+  ids = new char*[narg-6];
+  value2index = new int[narg-6];
+  nvalues = 0;
+
+  int iarg = 6;
+  while (iarg < narg) {
+    ids[nvalues] = NULL;
+
+    if (strcmp(arg[iarg],"x") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"y") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"z") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 2;
+
+    } else if (strcmp(arg[iarg],"vx") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"vy") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"vz") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 2;
+
+    } else if (strcmp(arg[iarg],"fx") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"fy") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"fz") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 2;
+
+    } else if ((strncmp(arg[iarg],"c_",2) == 0) || 
+	       (strncmp(arg[iarg],"f_",2) == 0) || 
+	       (strncmp(arg[iarg],"v_",2) == 0)) {
+      if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
+      else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
+      else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
+
+      int n = strlen(arg[iarg]);
+      char *suffix = new char[n];
+      strcpy(suffix,&arg[iarg][2]);
+
+      char *ptr = strchr(suffix,'[');
+      if (ptr) {
+	if (suffix[strlen(suffix)-1] != ']')
+	  error->all("Illegal fix ave/atom command");
+	argindex[nvalues] = atoi(ptr+1);
+	*ptr = '\0';
+      } else argindex[nvalues] = 0;
+
+      n = strlen(suffix) + 1;
+      ids[nvalues] = new char[n];
+      strcpy(ids[nvalues],suffix);
+      nvalues++;
+      delete [] suffix;
+
+    } else error->all("Illegal fix ave/atom command");
+
+    iarg++;
+  }
 
   // setup and error check
 
@@ -45,24 +119,49 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
       (nrepeat-1)*nevery >= peratom_freq)
     error->all("Illegal fix ave/atom command");
 
-  int icompute = modify->find_compute(id_compute);
+  for (int i = 0; i < nvalues; i++) {
-  if (icompute < 0) error->all("Compute ID for fix ave/atom does not exist");
+    if (which[i] == COMPUTE) {
+      int icompute = modify->find_compute(ids[i]);
+      if (icompute < 0)
+	error->all("Compute ID for fix ave/atom does not exist");
+      if (modify->compute[icompute]->peratom_flag == 0)
+	error->all("Fix ave/atom compute does not calculate per-atom values");
+      if (argindex[i] == 0 && modify->compute[icompute]->size_peratom != 0)
+	error->all("Fix ave/atom compute does not calculate a per-atom scalar");
+      if (argindex[i] && modify->compute[icompute]->size_peratom == 0)
+	error->all("Fix ave/atom compute does not calculate a per-atom vector");
+      if (argindex[i] && argindex[i] > modify->compute[icompute]->size_peratom)
+	error->all("Fix ave/atom compute vector is accessed out-of-range");
+    } else if (which[i] == FIX) {
+      int ifix = modify->find_fix(ids[i]);
+      if (ifix < 0)
+	error->all("Fix ID for fix ave/atom does not exist");
+      if (modify->fix[ifix]->peratom_flag == 0)
+	error->all("Fix ave/atom fix does not calculate per-atom values");
+      if (argindex[i] && modify->fix[ifix]->size_peratom != 0)
+	error->all("Fix ave/atom fix does not calculate a per-atom scalar");
+      if (argindex[i] && modify->fix[ifix]->size_peratom == 0)
+	error->all("Fix ave/atom fix does not calculate a per-atom vector");
+      if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom)
+	error->all("Fix ave/atom fix vector is accessed out-of-range");
+    } else if (which[i] == VARIABLE) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0)
+	error->all("Variable name for fix ave/atom does not exist");
+      if (input->variable->atomstyle(ivariable) == 0)
+	error->all("Fix ave/atom variable is not atom-style variable");
+    }
+  }
 
-  if (modify->compute[icompute]->peratom_flag == 0)
+  // this fix produces either a per-atom scalar or vector
-    error->all("Fix ave/atom compute does not calculate per-atom info");
 
   peratom_flag = 1;
-
+  if (nvalues == 1) size_peratom = 0;
-  // setup list of computes to call, including pre-computes
+  else size_peratom = nvalues;
-
-  ncompute = 1 + modify->compute[icompute]->npre;
-  compute = new Compute*[ncompute];
 
   // perform initial allocation of atom-based array
   // register with Atom class
 
-  size_peratom = modify->compute[icompute]->size_peratom;
-  scalar = NULL;
   vector = NULL;
   grow_arrays(atom->nmax);
   atom->add_callback(0);
@@ -70,13 +169,10 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
   // zero the array since dump may access it on timestep 0
 
   int nlocal = atom->nlocal;
-  if (size_peratom == 0)
+  for (int i = 0; i < nlocal; i++)
-    for (int i = 0; i < nlocal; i++) scalar[i] = 0.0;
+    for (int m = 0; m < nvalues; m++)
-  else
+      vector[i][m] = 0.0;
-    for (int i = 0; i < nlocal; i++)
+  
-      for (int m = 0; m < size_peratom; m++)
-	vector[i][m] = 0.0;
-
   // nvalid = next step on which end_of_step does something
   // can be this timestep if multiple of peratom_freq and nrepeat = 1
   // else backup from next multiple of peratom_freq
@@ -89,12 +185,11 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
     nvalid -= (nrepeat-1)*nevery;
   if (nvalid < update->ntimestep) nvalid += peratom_freq;
 
-  // set timestep for all computes that store invocation times
+  // add nvalid to ALL computes that store invocation times
-  //   since don't know a priori which are invoked by this fix
+  // since don't know a priori which are invoked by this fix
   // once in end_of_step() can just set timestep for ones actually invoked
 
-  for (int i = 0; i < modify->ncompute; i++)
+  modify->addstep_compute_all(nvalid);
-    if (modify->compute[i]->timeflag) modify->compute[i]->add_step(nvalid);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -105,9 +200,12 @@ FixAveAtom::~FixAveAtom()
  
   atom->delete_callback(id,0);
 
-  delete [] id_compute;
+  delete [] which;
-  delete [] compute;
+  delete [] argindex;
-  memory->sfree(scalar);
+  for (int m = 0; m < nvalues; m++) delete [] ids[m];
+  delete [] ids;
+  delete [] value2index;
+
   memory->destroy_2d_double_array(vector);
 }
 
@@ -124,29 +222,49 @@ int FixAveAtom::setmask()
 
 void FixAveAtom::init()
 {
-  // set ptrs to compute and its pre-computes called each end-of-step
+  // set indices and check validity of all computes,fixes,variables
-  // put pre-computes in list before compute
+  // check that fix frequency is acceptable
 
-  int icompute = modify->find_compute(id_compute);
+  for (int m = 0; m < nvalues; m++) {
-  if (icompute < 0) error->all("Compute ID for fix ave/atom does not exist");
+    if (which[m] == COMPUTE) {
-  
+      int icompute = modify->find_compute(ids[m]);
-  ncompute = 0;
+      if (icompute < 0)
-  if (modify->compute[icompute]->npre)
+	error->all("Compute ID for fix ave/atom does not exist");
-    for (int i = 0; i < modify->compute[icompute]->npre; i++) {
+      value2index[m] = icompute;
-      int ic = modify->find_compute(modify->compute[icompute]->id_pre[i]);
+      
-      if (ic < 0)
+    } else if (which[m] == FIX) {
-	error->all("Precompute ID of compute for fix ave/atom does not exist");
+      int ifix = modify->find_fix(ids[m]);
-      compute[ncompute++] = modify->compute[ic];
+      if (ifix < 0) 
-    }
+	error->all("Fix ID for fix ave/atom does not exist");
+      value2index[m] = ifix;
 
-  compute[ncompute++] = modify->compute[icompute];
+      if (nevery % modify->fix[ifix]->peratom_freq)
+	error->all("Fix for fix ave/atom not computed at compatible time");
+
+    } else if (which[m] == VARIABLE) {
+      int ivariable = input->variable->find(ids[m]);
+      if (ivariable < 0) 
+	error->all("Variable name for fix ave/atom does not exist");
+      value2index[m] = ivariable;
+
+    } else value2index[m] = -1;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   only does something if nvalid = current timestep
+------------------------------------------------------------------------- */
+
+void FixAveAtom::setup()
+{
+  end_of_step();
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixAveAtom::end_of_step()
 {
-  int i,m;
+  int i,j,m,n;
 
   // skip if not step which requires doing something
 
@@ -156,32 +274,70 @@ void FixAveAtom::end_of_step()
 
   int nlocal = atom->nlocal;
 
-  if (irepeat == 0) {
+  if (irepeat == 0)
-    if (size_peratom == 0)
+    for (i = 0; i < nlocal; i++)
-      for (i = 0; i < nlocal; i++) scalar[i] = 0.0;
+      for (m = 0; m < nvalues; m++)
-    else
+	vector[i][m] = 0.0;
-      for (i = 0; i < nlocal; i++)
-	for (m = 0; m < size_peratom; m++)
-	  vector[i][m] = 0.0;
-  }
-  
-  // accumulate results of compute to local copy
   
+  // accumulate results of attributes,computes,fixes,variables to local copy
+  // compute/fix/variable may invoke computes so wrap with clear/add
+
   modify->clearstep_compute();
-  for (i = 0; i < ncompute; i++) compute[i]->compute_peratom();
+
-  
   int *mask = atom->mask;
 
-  if (size_peratom == 0) {
+  for (m = 0; m < nvalues; m++) {
-    double *compute_scalar = compute[ncompute-1]->scalar_atom;
+    n = value2index[m];
-    for (i = 0; i < nlocal; i++)
+    j = argindex[m];
-      if (mask[i] & groupbit) scalar[i] += compute_scalar[i];
+
-  } else {
+    if (which[m] == X) {
-    double **compute_vector = compute[ncompute-1]->vector_atom;
+      double **x = atom->x;
-    for (i = 0; i < nlocal; i++)
+      for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit)
+	if (mask[i] & groupbit) vector[i][m] += x[i][j];
-	for (m = 0; m < size_peratom; m++)
+    } else if (which[m] == V) {
-	  vector[i][m] += compute_vector[i][m];
+      double **v = atom->v;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) vector[i][m] += v[i][j];
+    } else if (which[m] == F) {
+      double **f = atom->f;
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) vector[i][m] += f[i][j];
+
+    // invoke compute if not previously invoked
+
+    } else if (which[m] == COMPUTE) {
+      Compute *compute = modify->compute[n];
+      if (!(compute->invoked & INVOKED_PERATOM)) compute->compute_peratom();
+
+      if (j == 0) {
+	double *compute_scalar = compute->scalar_atom;
+	for (i = 0; i < nlocal; i++)
+	  if (mask[i] & groupbit) vector[i][m] += compute_scalar[i];
+      } else {
+	int jm1 = j - 1;
+	double **compute_vector = compute->vector_atom;
+	for (i = 0; i < nlocal; i++)
+	  if (mask[i] & groupbit) vector[i][m] += compute_vector[i][jm1];
+      }
+
+    // access fix fields, guaranteed to be ready
+
+    } else if (which[m] == FIX) {
+      if (j == 0) {
+	double *fix_scalar = modify->fix[n]->scalar_atom;
+	for (i = 0; i < nlocal; i++)
+	  if (mask[i] & groupbit) vector[i][m] += fix_scalar[i];
+      } else {
+	int jm1 = j - 1;
+	double **fix_vector = modify->fix[n]->vector_atom;
+	for (i = 0; i < nlocal; i++)
+	  if (mask[i] & groupbit) vector[i][m] += fix_vector[i][jm1];
+      }
+
+    // evaluate atom-style variable
+
+    } else if (which[m] == VARIABLE)
+      input->variable->compute_atom(n,igroup,&vector[0][m],nvalues,1);
   }
 
   // done if irepeat < nrepeat
@@ -201,13 +357,9 @@ void FixAveAtom::end_of_step()
   // average the final result for the Nfreq timestep
 
   double repeat = nrepeat;
-  if (size_peratom == 0)
+  for (i = 0; i < nlocal; i++)
-    for (i = 0; i < nlocal; i++)
+    for (m = 0; m < nvalues; m++)
-      scalar[i] /= repeat;
+      vector[i][m] /= repeat;
-  else
-    for (i = 0; i < nlocal; i++)
-      for (m = 0; m < size_peratom; m++)
-	vector[i][m] /= repeat;
 }
 
 /* ----------------------------------------------------------------------
@@ -217,8 +369,7 @@ void FixAveAtom::end_of_step()
 double FixAveAtom::memory_usage()
 {
   double bytes;
-  if (size_peratom == 0) bytes = atom->nmax * sizeof(double);
+  bytes = atom->nmax*nvalues * sizeof(double);
-  else bytes = atom->nmax*size_peratom * sizeof(double);
   return bytes;
 }
 
@@ -228,15 +379,10 @@ double FixAveAtom::memory_usage()
 
 void FixAveAtom::grow_arrays(int nmax)
 {
-  if (size_peratom == 0) {
+  vector = memory->grow_2d_double_array(vector,nmax,nvalues,
-    scalar = (double *) memory->srealloc(scalar,nmax*sizeof(double),
+					"fix_ave/atom:vector");
-					 "fix_ave/atom:scalar");
+  vector_atom = vector;
-    scalar_atom = scalar;
+  scalar_atom = vector[0];
-  } else {
-    vector = memory->grow_2d_double_array(vector,nmax,size_peratom,
-					  "fix_ave/atom:vector");
-    vector_atom = vector;
-  }
 }
 
 /* ----------------------------------------------------------------------
@@ -245,11 +391,8 @@ void FixAveAtom::grow_arrays(int nmax)
 
 void FixAveAtom::copy_arrays(int i, int j)
 {
-  if (size_peratom == 0)
+  for (int m = 0; m < nvalues; m++)
-    scalar[j] = scalar[i];
+    vector[j][m] = vector[i][m];
-  else
-    for (int m = 0; m <= size_peratom; m++)
-      vector[j][m] = vector[i][m];
 }
 
 /* ----------------------------------------------------------------------
@@ -258,13 +401,8 @@ void FixAveAtom::copy_arrays(int i, int j)
 
 int FixAveAtom::pack_exchange(int i, double *buf)
 {
-  if (size_peratom == 0) {
+  for (int m = 0; m < nvalues; m++) buf[m] = vector[i][m];
-    buf[0] = scalar[i];
+  return nvalues;
-    return 1;
-  }
-
-  for (int m = 0; m <= size_peratom; m++) buf[m] = vector[i][m];
-  return size_peratom;
 }
 
 /* ----------------------------------------------------------------------
@@ -273,11 +411,6 @@ int FixAveAtom::pack_exchange(int i, double *buf)
 
 int FixAveAtom::unpack_exchange(int nlocal, double *buf)
 {
-  if (size_peratom == 0) {
+  for (int m = 0; m < nvalues; m++) vector[nlocal][m] = buf[m];
-    scalar[nlocal] = buf[0];
+  return nvalues;
-    return 1;
-  }
-
-  for (int m = 0; m <= size_peratom; m++) vector[nlocal][m] = buf[m];
-  return size_peratom;
 }

src/fix_ave_atom.h

@@ -25,6 +25,7 @@ class FixAveAtom : public Fix {
   ~FixAveAtom();
   int setmask();
   void init();
+  void setup();
   void end_of_step();
 
   double memory_usage();
@@ -34,12 +35,11 @@ class FixAveAtom : public Fix {
   int unpack_exchange(int, double *);
 
  private:
-  int nrepeat,irepeat,nvalid;
+  int nvalues;
-  char *id_compute;
+  int nrepeat,nvalid,irepeat;
-  int ncompute;
+  int *which,*argindex,*value2index;
-  class Compute **compute;
+  char **ids;
 
-  double *scalar;
   double **vector;
 };
 

src/fix_ave_spatial.cpp

@@ -24,26 +24,28 @@
 #include "lattice.h"
 #include "modify.h"
 #include "compute.h"
-#include "group.h"
+#include "input.h"
+#include "variable.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 enum{LOWER,CENTER,UPPER,COORD};
-enum{DENSITY_MASS,DENSITY_NUM,COMPUTE,FIX};
+enum{X,V,F,DENSITY_NUMBER,DENSITY_MASS,COMPUTE,FIX,VARIABLE};
 enum{SAMPLE,ALL};
 enum{BOX,LATTICE,REDUCED};
 enum{ONE,RUNNING,WINDOW};
 
 #define BIG 1000000000
+#define INVOKED_PERATOM 4      // same as in computes
 
 /* ---------------------------------------------------------------------- */
 
 FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg < 11) error->all("Illegal fix ave/spatial command");
+  if (narg < 10) error->all("Illegal fix ave/spatial command");
 
   MPI_Comm_rank(world,&me);
 
@@ -62,34 +64,96 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
   if (strcmp(arg[7],"center") == 0) originflag = CENTER;
   if (strcmp(arg[7],"upper") == 0) originflag = UPPER;
   else originflag = COORD;
-  if (originflag == COORD) origin = atof(arg[6]);
+  if (originflag == COORD) origin = atof(arg[7]);
 
   delta = atof(arg[8]);
 
-  if (strcmp(arg[9],"density") == 0) {
+  // parse values until one isn't recognized
-    if (strcmp(arg[10],"mass") == 0) which = DENSITY_MASS;
-    else if (strcmp(arg[10],"number") == 0) which = DENSITY_NUM;
-    else error->all("Illegal fix ave/spatial command");
-  } else if (strcmp(arg[9],"compute") == 0) {
-    which = COMPUTE;
-    int n = strlen(arg[10]) + 1;
-    id_compute = new char[n];
-    strcpy(id_compute,arg[10]);
-  } else if (strcmp(arg[9],"fix") == 0) {
-    which = FIX;
-    int n = strlen(arg[10]) + 1;
-    id_fix = new char[n];
-    strcpy(id_fix,arg[10]);
-  } else error->all("Illegal fix ave/spatial command");
 
-  // parse optional args
+  which = new int[narg-9];
+  argindex = new int[narg-9];
+  ids = new char*[narg-9];
+  value2index = new int[narg-9];
+  nvalues = 0;
+
+  int iarg = 9;
+  while (iarg < narg) {
+    ids[nvalues] = NULL;
+
+    if (strcmp(arg[iarg],"x") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"y") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"z") == 0) {
+      which[nvalues] = X;
+      argindex[nvalues++] = 2;
+
+    } else if (strcmp(arg[iarg],"vx") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"vy") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"vz") == 0) {
+      which[nvalues] = V;
+      argindex[nvalues++] = 2;
+
+    } else if (strcmp(arg[iarg],"fx") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"fy") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 1;
+    } else if (strcmp(arg[iarg],"fz") == 0) {
+      which[nvalues] = F;
+      argindex[nvalues++] = 2;
+
+    } else if (strcmp(arg[iarg],"density/number") == 0) {
+      which[nvalues] = DENSITY_NUMBER;
+      argindex[nvalues++] = 0;
+    } else if (strcmp(arg[iarg],"density/mass") == 0) {
+      which[nvalues] = DENSITY_MASS;
+      argindex[nvalues++] = 0;
+
+    } else if ((strncmp(arg[iarg],"c_",2) == 0) || 
+	       (strncmp(arg[iarg],"f_",2) == 0) || 
+	       (strncmp(arg[iarg],"v_",2) == 0)) {
+      if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
+      else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
+      else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
+
+      int n = strlen(arg[iarg]);
+      char *suffix = new char[n];
+      strcpy(suffix,&arg[iarg][2]);
+
+      char *ptr = strchr(suffix,'[');
+      if (ptr) {
+	if (suffix[strlen(suffix)-1] != ']')
+	  error->all("Illegal fix ave/spatial command");
+	argindex[nvalues] = atoi(ptr+1);
+	*ptr = '\0';
+      } else argindex[nvalues] = 0;
+
+      n = strlen(suffix) + 1;
+      ids[nvalues] = new char[n];
+      strcpy(ids[nvalues],suffix);
+      nvalues++;
+      delete [] suffix;
+
+    } else break;
+
+    iarg++;
+  }
+
+  // optional args
 
   normflag = ALL;
   scaleflag = BOX;
   fp = NULL;
   ave = ONE;
 
-  int iarg = 11;
   while (iarg < narg) {
     if (strcmp(arg[iarg],"norm") == 0) {
       if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
@@ -131,11 +195,71 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
     } else error->all("Illegal fix ave/spatial command");
   }
 
-  // if density, no normalization by atom count should be done
+  // setup and error check
-  // thus ALL and SAMPLE should give same answer, but code does normalize
-  // thus only ALL is computed correctly, so force norm to be ALL
 
-  if (which == DENSITY_MASS || which == DENSITY_NUM) normflag = ALL;
+  if (nevery <= 0) error->all("Illegal fix ave/spatial command");
+  if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
+    error->all("Illegal fix ave/spatial command");
+
+  if (delta <= 0.0) error->all("Illegal fix ave/spatial command");
+  invdelta = 1.0/delta;
+
+  for (int i = 0; i < nvalues; i++) {
+    if (which[i] == COMPUTE) {
+      int icompute = modify->find_compute(ids[i]);
+      if (icompute < 0)
+	error->all("Compute ID for fix ave/spatial does not exist");
+      if (modify->compute[icompute]->peratom_flag == 0)
+	error->all("Fix ave/spatial compute does not calculate per-atom values");
+      if (argindex[i] == 0 && modify->compute[icompute]->size_peratom != 0)
+	error->all("Fix ave/spatial compute does not calculate a per-atom scalar");
+      if (argindex[i] && modify->compute[icompute]->size_peratom == 0)
+	error->all("Fix ave/spatial compute does not calculate a per-atom vector");
+      if (argindex[i] && argindex[i] > modify->compute[icompute]->size_peratom)
+	error->all("Fix ave/spatial compute vector is accessed out-of-range");
+
+    } else if (which[i] == FIX) {
+      int ifix = modify->find_fix(ids[i]);
+      if (ifix < 0)
+	error->all("Fix ID for fix ave/spatial does not exist");
+      if (modify->fix[ifix]->peratom_flag == 0)
+	error->all("Fix ave/spatial fix does not calculate per-atom values");
+      if (argindex[i] && modify->fix[ifix]->size_peratom != 0)
+	error->all("Fix ave/spatial fix does not calculate a per-atom scalar");
+      if (argindex[i] && modify->fix[ifix]->size_peratom == 0)
+	error->all("Fix ave/spatial fix does not calculate a per-atom vector");
+      if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom)
+	error->all("Fix ave/spatial fix vector is accessed out-of-range");
+    } else if (which[i] == VARIABLE) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0)
+	error->all("Variable name for fix ave/spatial does not exist");
+      if (input->variable->atomstyle(ivariable) == 0)
+	error->all("Fix ave/spatial variable is not atom-style variable");
+    }
+  }
+
+  // print header into file
+
+  if (fp && me == 0) {
+    fprintf(fp,"Spatial-averaged data for fix %s and group %s\n",id,arg[1]);
+    fprintf(fp,"TimeStep Number-of-layers\n");
+    fprintf(fp,"Layer Coordinate Natoms");
+    for (int i = 0; i < nvalues; i++)
+      if (which[i] == COMPUTE) fprintf(fp," c_%s",ids[i]);
+      else if (which[i] == FIX) fprintf(fp," f_%s",ids[i]);
+      else if (which[i] == VARIABLE) fprintf(fp," v_%s",ids[i]);
+      else fprintf(fp," %s",arg[9+i]);
+    fprintf(fp,"\n");
+  }
+  
+  // this fix produces a global vector
+  // set size_vector to BIG since compute_vector() checks bounds on-the-fly
+
+  vector_flag = 1;
+  size_vector = BIG;
+  scalar_vector_freq = nfreq;
+  extvector = 0;
 
   // setup scaling
 
@@ -162,60 +286,6 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
   delta *= scale;
   if (originflag == COORD) origin *= scale;
 
-  // setup and error check
-
-  if (nevery <= 0) error->all("Illegal fix ave/spatial command");
-  if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
-    error->all("Illegal fix ave/spatial command");
-
-  if (delta <= 0.0) error->all("Illegal fix ave/spatial command");
-  invdelta = 1.0/delta;
-
-  // nvalues = # of quantites per line of output file
-  // for COMPUTE, setup list of computes to call, including pre-computes
-
-  nvalues = 1;
-  compute = NULL;
-
-  if (which == COMPUTE) {
-    int icompute = modify->find_compute(id_compute);
-    if (icompute < 0)
-      error->all("Compute ID for fix ave/spatial does not exist");
-    if (modify->compute[icompute]->peratom_flag == 0)
-      error->all("Fix ave/spatial compute does not calculate per-atom info");
-    nvalues = size_peratom = modify->compute[icompute]->size_peratom;
-    if (nvalues == 0) nvalues = 1;
-    ncompute = 1 + modify->compute[icompute]->npre;
-    compute = new Compute*[ncompute];
-  }
-
-  if (which == FIX) {
-    int ifix = modify->find_fix(id_fix);
-    if (ifix < 0)
-      error->all("Fix ID for fix ave/spatial does not exist");
-    if (modify->fix[ifix]->peratom_flag == 0)
-      error->all("Fix ave/spatial fix does not calculate per-atom info");
-    nvalues = size_peratom = modify->fix[ifix]->size_peratom;
-    if (nvalues == 0) nvalues = 1;
-  }
-
-  // print header into file
-
-  if (fp && me == 0) {
-    fprintf(fp,"Spatial-averaged data for fix %s, group %s, and %s %s\n",
-	    id,group->names[igroup],arg[10],arg[11]);
-    fprintf(fp,"TimeStep Number-of-layers (one per snapshot)\n");
-    fprintf(fp,"Layer Coord Atoms Value(s) (one per layer)\n");
-  }
-  
-  // enable this fix to produce a global vector
-  // set size_vector to BIG since compute_vector() will check bounds
-
-  vector_flag = 1;
-  size_vector = BIG;
-  scalar_vector_freq = nfreq;
-  extensive = 0;
-
   // initializations
 
   irepeat = 0;
@@ -228,6 +298,12 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
   values_one = values_many = values_sum = values_total = NULL;
   values_list = NULL;
 
+  maxatomvar = 0;
+  varatom = NULL;
+
+  maxatomlayer = 0;
+  layer = NULL;
+
   // nvalid = next step on which end_of_step does something
   // can be this timestep if multiple of nfreq and nrepeat = 1
   // else backup from next multiple of nfreq
@@ -239,24 +315,24 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
     nvalid -= (nrepeat-1)*nevery;
   if (nvalid < update->ntimestep) nvalid += nfreq;
 
-  // set timestep for all computes that store invocation times
+  // add nvalid to ALL computes that store invocation times
-  //   since don't know a priori which are invoked by this fix
+  // since don't know a priori which are invoked by this fix
   // once in end_of_step() can just set timestep for ones actually invoked
 
-  if (which == COMPUTE)
+  modify->addstep_compute_all(nvalid);
-    for (int i = 0; i < modify->ncompute; i++)
-      if (modify->compute[i]->timeflag) modify->compute[i]->add_step(nvalid);
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixAveSpatial::~FixAveSpatial()
 {
-  if (which == COMPUTE) delete [] id_compute;
+  delete [] which;
-  if (which == FIX) delete [] id_fix;
+  delete [] argindex;
-  if (fp && me == 0) fclose(fp);
+  for (int i = 0; i < nvalues; i++) delete [] ids[i];
+  delete [] ids;
+  delete [] value2index;
 
-  delete [] compute;
+  if (fp && me == 0) fclose(fp);
 
   memory->sfree(coord);
   memory->sfree(count_one);
@@ -269,6 +345,9 @@ FixAveSpatial::~FixAveSpatial()
   memory->destroy_2d_double_array(values_sum);
   memory->destroy_2d_double_array(values_total);
   memory->destroy_3d_double_array(values_list);
+
+  memory->sfree(varatom);
+  memory->sfree(layer);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -291,36 +370,32 @@ void FixAveSpatial::init()
       error->all("Fix ave/spatial settings invalid with changing box");
   }
 
-  // set ptrs to compute and its pre-computes called each end-of-step
+  // set indices and check validity of all computes,fixes,variables
-  // put pre-computes in list before compute
-
-  if (which == COMPUTE) {
-    int icompute = modify->find_compute(id_compute);
-    if (icompute < 0)
-      error->all("Compute ID for fix ave/spatial does not exist");
-    
-    ncompute = 0;
-    if (modify->compute[icompute]->npre)
-      for (int i = 0; i < modify->compute[icompute]->npre; i++) {
-	int ic = modify->find_compute(modify->compute[icompute]->id_pre[i]);
-	if (ic < 0)
-	  error->all("Precompute ID for fix ave/spatial does not exist");
-	compute[ncompute++] = modify->compute[ic];
-      }
-    
-    compute[ncompute++] = modify->compute[icompute];
-  }
-
-  // set ptr to fix ID
   // check that fix frequency is acceptable
 
-  if (which == FIX) {
+  for (int m = 0; m < nvalues; m++) {
-    int ifix = modify->find_fix(id_fix);
+    if (which[m] == COMPUTE) {
-    if (ifix < 0) 
+      int icompute = modify->find_compute(ids[m]);
-      error->all("Fix ID for fix ave/spatial does not exist");
+      if (icompute < 0)
-    fix = modify->fix[ifix];
+	error->all("Compute ID for fix ave/spatial does not exist");
-    if (nevery % fix->peratom_freq)
+      value2index[m] = icompute;
-      error->all("Fix ave/spatial and fix not computed at compatible times");
+      
+    } else if (which[m] == FIX) {
+      int ifix = modify->find_fix(ids[m]);
+      if (ifix < 0) 
+	error->all("Fix ID for fix ave/spatial does not exist");
+      value2index[m] = ifix;
+
+      if (nevery % modify->fix[ifix]->peratom_freq)
+	error->all("Fix for fix ave/spatial not computed at compatible time");
+
+    } else if (which[m] == VARIABLE) {
+      int ivariable = input->variable->find(ids[m]);
+      if (ivariable < 0) 
+	error->all("Variable name for fix ave/spatial does not exist");
+      value2index[m] = ivariable;
+
+    } else value2index[m] = -1;
   }
 }
 
@@ -337,7 +412,7 @@ void FixAveSpatial::setup()
 
 void FixAveSpatial::end_of_step()
 {
-  int i,j,m,ilayer;
+  int i,j,m,n,ilayer;
   double lo,hi;
 
   // skip if not step which requires doing something
@@ -448,107 +523,129 @@ void FixAveSpatial::end_of_step()
     for (i = 0; i < nvalues; i++) values_one[m][i] = 0.0;
   }
 
-  // perform the computation for one sample
+  // assign each atom to a layer
-  // sum within each layer, only include atoms in fix group
   // insure array index is within bounds (since atoms can be outside box)
-  // if scaleflag = REDUCED, box coords -> lamda coords before computing layer
+  // if scaleflag = REDUCED, box coords -> lamda coords
 
   double **x = atom->x;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  // DENSITY_MASS adds mass to values
+  if (nlocal > maxatomlayer) {
-
+    maxatomlayer = atom->nmax;
-  if (which == DENSITY_MASS) {
+    memory->sfree(layer);
-    int *type = atom->type;
+    layer = (int *) 
-    double *mass = atom->mass;
+      memory->smalloc(maxatomlayer*sizeof(int),"ave/spatial:layer");
-    double *rmass = atom->rmass;
-
-    if (scaleflag == REDUCED) domain->x2lamda(nlocal);
-
-    for (i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-	ilayer = static_cast<int> ((x[i][dim] - offset) * invdelta);
-	if (ilayer < 0) ilayer = 0;
-	if (ilayer >= nlayers) ilayer = nlayers-1;
-	count_one[ilayer] += 1.0;
-	if (mass) values_one[ilayer][0] += mass[type[i]];
-	else values_one[ilayer][0] += rmass[i];
-      }
-    }
-
-    if (scaleflag == REDUCED) domain->lamda2x(nlocal);
-
-  // DENSITY_NUM adds 1 to values
-
-  } else if (which == DENSITY_NUM) {
-
-    if (scaleflag == REDUCED) domain->x2lamda(nlocal);
-
-    for (i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-	ilayer = static_cast<int> ((x[i][dim] - offset) * invdelta);
-	if (ilayer < 0) ilayer = 0;
-	if (ilayer >= nlayers) ilayer = nlayers-1;
-	count_one[ilayer] += 1.0;
-	values_one[ilayer][0] += 1.0;
-      }
-    }
-
-    if (scaleflag == REDUCED) domain->lamda2x(nlocal);
-
-  // COMPUTE adds its scalar or vector quantity to values
-
-  } else if (which == COMPUTE) {
-    modify->clearstep_compute();
-    for (i = 0; i < ncompute; i++) compute[i]->compute_peratom();
-    double *scalar = compute[ncompute-1]->scalar_atom;
-    double **vector = compute[ncompute-1]->vector_atom;
-
-    if (scaleflag == REDUCED) domain->x2lamda(nlocal);
-
-    m = 0;
-    for (i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-	ilayer = static_cast<int> ((x[i][dim] - offset) * invdelta);
-	if (ilayer < 0) ilayer = 0;
-	if (ilayer >= nlayers) ilayer = nlayers-1;
-	count_one[ilayer] += 1.0;
-	if (size_peratom == 0) values_one[ilayer][0] += scalar[i];
-	else
-	  for (j = 0; j < nvalues; j++)
-	    values_one[ilayer][j] += vector[i][j];
-      }
-    }
-
-    if (scaleflag == REDUCED) domain->lamda2x(nlocal);
-
-  // FIX adds its scalar or vector quantity to values
-
-  } else if (which == FIX) {
-    double *scalar = fix->scalar_atom;
-    double **vector = fix->vector_atom;
-
-    if (scaleflag == REDUCED) domain->x2lamda(nlocal);
-
-    m = 0;
-    for (i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-	ilayer = static_cast<int> ((x[i][dim] - offset) * invdelta);
-	if (ilayer < 0) ilayer = 0;
-	if (ilayer >= nlayers) ilayer = nlayers-1;
-	count_one[ilayer] += 1.0;
-	if (size_peratom == 0) values_one[ilayer][0] += scalar[i];
-	else
-	  for (j = 0; j < nvalues; j++)
-	    values_one[ilayer][j] += vector[i][j];
-      }
-    }
-
-    if (scaleflag == REDUCED) domain->lamda2x(nlocal);
   }
 
-  // average a single sample
+  if (scaleflag == REDUCED) domain->x2lamda(nlocal);
+  
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      ilayer = static_cast<int> ((x[i][dim] - offset) * invdelta);
+      if (ilayer < 0) ilayer = 0;
+      if (ilayer >= nlayers) ilayer = nlayers-1;
+      layer[i] = ilayer;
+      count_one[ilayer] += 1.0;
+    }
+
+  if (scaleflag == REDUCED) domain->lamda2x(nlocal);
+
+  // perform the computation for one sample
+  // accumulate results of attributes,computes,fixes,variables to local copy
+  // sum within each layer, only include atoms in fix group
+  // compute/fix/variable may invoke computes so wrap with clear/add
+
+  modify->clearstep_compute();
+
+  for (m = 0; m < nvalues; m++) {
+    n = value2index[m];
+    j = argindex[m];
+
+    // X,V,F adds coords,velocities,forces to values
+
+    if (which[m] == X || which[m] == V || which[m] == F) {
+      double **attribute;
+      if (which[m] == X) attribute = x;
+      else if (which[m] == V) attribute = atom->v;
+      else attribute = atom->f;
+
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit)
+	  values_one[layer[i]][m] += attribute[i][j];
+
+    // DENSITY_NUMBER adds 1 to values
+
+    } else if (which[m] == DENSITY_NUMBER) {
+
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit)
+	  values_one[layer[i]][m] += 1.0;
+
+    // DENSITY_MASS adds mass to values
+
+    } else if (which[m] == DENSITY_MASS) {
+      int *type = atom->type;
+      double *mass = atom->mass;
+      double *rmass = atom->rmass;
+
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit)
+	  if (mass) values_one[layer[i]][m] += mass[type[i]];
+	  else values_one[layer[i]][m] += rmass[i];
+
+    // COMPUTE adds its scalar or vector component to values
+    // invoke compute if not previously invoked
+
+    } else if (which[m] == COMPUTE) {
+      Compute *compute = modify->compute[n];
+      if (!(compute->invoked & INVOKED_PERATOM)) compute->compute_peratom();
+      double *scalar = compute->scalar_atom;
+      double **vector = compute->vector_atom;
+      int jm1 = j - 1;
+
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit)
+	  if (j == 0) values_one[layer[i]][m] += scalar[i];
+	  else values_one[layer[i]][m] += vector[i][jm1];
+      
+    // FIX adds its scalar or vector component to values
+    // access fix fields, guaranteed to be ready
+
+    } else if (which[m] == FIX) {
+      double *scalar = modify->fix[n]->scalar_atom;
+      double **vector = modify->fix[n]->vector_atom;
+      int jm1 = j - 1;
+
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit) {
+	  if (j == 0) values_one[layer[i]][m] += scalar[i];
+	  else values_one[layer[i]][m] += vector[i][jm1];
+	}
+
+    // VARIABLE adds its per-atom quantities to values
+    // evaluate atom-style variable
+
+    } else if (which[m] == VARIABLE) {
+      if (nlocal > maxatomvar) {
+	maxatomvar = atom->nmax;
+	memory->sfree(varatom);
+	varatom = (double *) 
+	  memory->smalloc(maxatomvar*sizeof(double),"ave/spatial:varatom");
+      }
+
+      input->variable->compute_atom(n,igroup,varatom,1,0);
+
+      for (i = 0; i < nlocal; i++)
+	if (mask[i] & groupbit)
+	  values_one[layer[i]][m] += varatom[i];
+    }
+  }
+
+  // process a single sample
+  // if normflag = ALL, accumulate values,count separately to many
+  // if normflag = SAMPLE, one = value/count, accumulate one to many
+  // exception is SAMPLE density: no normalization by atom count
 
   if (normflag == ALL) {
     for (m = 0; m < nlayers; m++) {
@@ -560,8 +657,11 @@ void FixAveSpatial::end_of_step()
     MPI_Allreduce(count_one,count_many,nlayers,MPI_DOUBLE,MPI_SUM,world);
     for (m = 0; m < nlayers; m++) {
       if (count_many[m] > 0.0)
-	for (j = 0; j < nvalues; j++)
+	for (j = 0; j < nvalues; j++) {
-	  values_many[m][j] += values_one[m][j]/count_many[m];
+	  if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS)
+	    values_many[m][j] += values_one[m][j];
+	  else values_many[m][j] += values_one[m][j]/count_many[m];
+	}
       count_sum[m] += count_many[m];
     }
   }
@@ -572,16 +672,18 @@ void FixAveSpatial::end_of_step()
   irepeat++;
   if (irepeat < nrepeat) {
     nvalid += nevery;
-    if (which == COMPUTE) modify->addstep_compute(nvalid);
+    modify->addstep_compute(nvalid);
     return;
   }
 
   irepeat = 0;
   nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery;
-  if (which == COMPUTE) modify->addstep_compute(nvalid);
+  modify->addstep_compute(nvalid);
 
   // time average across samples
-  // if density, also normalize by volume
+  // if normflag = ALL, final is total value / total count
+  // if normflag = SAMPLE, final is sum of ave / repeat
+  // exception is ALL density: normalized by repeat, not total count
 
   double repeat = nrepeat;
 
@@ -592,7 +694,9 @@ void FixAveSpatial::end_of_step()
     for (m = 0; m < nlayers; m++) {
       if (count_sum[m] > 0.0)
 	for (j = 0; j < nvalues; j++)
-	  values_sum[m][j] /= count_sum[m];
+	  if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS)
+	    values_sum[m][j] /= repeat;
+	  else values_sum[m][j] /= count_sum[m];
       count_sum[m] /= repeat;
     }
   } else {
@@ -605,10 +709,12 @@ void FixAveSpatial::end_of_step()
     }
   }
 
-  if (which == DENSITY_MASS || which == DENSITY_NUM) {
+  // density is additionally normalized by layer volume
-    for (m = 0; m < nlayers; m++)
+
-      values_sum[m][0] *= count_sum[m] / layer_volume;
+  for (j = 0; j < nvalues; j++)
-  }
+    if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS)
+      for (m = 0; m < nlayers; m++)
+	values_sum[m][j] /= layer_volume;
 
   // if ave = ONE, only single Nfreq timestep value is needed
   // if ave = RUNNING, combine with all previous Nfreq timestep values
@@ -667,13 +773,24 @@ void FixAveSpatial::end_of_step()
 /* ----------------------------------------------------------------------
    return Nth vector value
    since values_sum is 2d array, map N into ilayer and ivalue
-   if ilayer >= nlayers, just return 0, since nlayers can vary with time
+   if ilayer exceeds current layers, return 0.0 instead of generate an error
 ------------------------------------------------------------------------- */
 
 double FixAveSpatial::compute_vector(int n)
 {
   int ivalue = n % nvalues;
   int ilayer = n / nvalues;
-  if (ilayer < nlayers && norm) return values_total[ilayer][ivalue]/norm;
+  if (ilayer >= nlayers) return 0.0;
-  return 0.0;
+  return values_total[ilayer][ivalue]/norm;
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of varatom and layer
+------------------------------------------------------------------------- */
+
+double FixAveSpatial::memory_usage()
+{
+  double bytes = maxatomvar * sizeof(double);
+  bytes += maxatomlayer * sizeof(int);
+  return bytes;
 }

src/fix_ave_spatial.h

@@ -28,26 +28,31 @@ class FixAveSpatial : public Fix {
   void setup();
   void end_of_step();
   double compute_vector(int);
+  double memory_usage();
 
  private:
-  int me;
+  int me,nvalues;
-  int nrepeat,nfreq,irepeat,nvalid;
+  int nrepeat,nfreq,nvalid,irepeat;
-  int dim,originflag,which,normflag;
+  int dim,originflag,normflag;
   double origin,delta;
-  char *id_compute,*id_fix;
+  int *which,*argindex,*value2index;
+  char **ids;
   FILE *fp;
 
-  int nlayers,nvalues,ave,nwindow;
+  int nlayers,ave,nwindow;
-  int maxlayer,scaleflag,size_peratom;
+  int maxlayer,scaleflag;
   double xscale,yscale,zscale;
   double layer_volume;
   double *coord;
   double *count_one,*count_many,*count_sum;
   double **values_one,**values_many,**values_sum;
   double offset,invdelta; 
-  int ncompute;
+
-  class Compute **compute;
+  int maxatomvar;
-  class Fix *fix;
+  double *varatom;
+
+  int maxatomlayer;
+  int *layer;
 
   int norm,iwindow,window_limit;
   double *count_total;

src/fix_ave_time.cpp

@@ -22,21 +22,25 @@
 #include "modify.h"
 #include "compute.h"
 #include "group.h"
+#include "input.h"
+#include "variable.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
-enum{COMPUTE,FIX};
+enum{COMPUTE,FIX,VARIABLE};
-enum{SCALAR,VECTOR,BOTH};
 enum{ONE,RUNNING,WINDOW};
 
+#define INVOKED_SCALAR 1      // same as in computes
+#define INVOKED_VECTOR 2
+
 /* ---------------------------------------------------------------------- */
 
 FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg < 8) error->all("Illegal fix ave/time command");
+  if (narg < 7) error->all("Illegal fix ave/time command");
 
   MPI_Comm_rank(world,&me);
 
@@ -44,37 +48,55 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   nrepeat = atoi(arg[4]);
   nfreq = atoi(arg[5]);
 
-  if (strcmp(arg[6],"compute") == 0) which = COMPUTE;
+  // parse values until one isn't recognized
-  else if (strcmp(arg[6],"fix") == 0) which = FIX;
-  else error->all("Illegal fix ave/time command");
 
-  int n = strlen(arg[7]) + 1;
+  which = new int[narg-6];
-  id = new char[n];
+  argindex = new int[narg-6];
-  strcpy(id,arg[7]);
+  ids = new char*[narg-6];
+  value2index = new int[narg-6];
+  nvalues = 0;
+
+  int iarg = 6;
+  while (iarg < narg) {
+    if ((strncmp(arg[iarg],"c_",2) == 0) || 
+	(strncmp(arg[iarg],"f_",2) == 0) || 
+	(strncmp(arg[iarg],"v_",2) == 0)) {
+      if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
+      else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
+      else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
+
+      int n = strlen(arg[iarg]);
+      char *suffix = new char[n];
+      strcpy(suffix,&arg[iarg][2]);
+
+      char *ptr = strchr(suffix,'[');
+      if (ptr) {
+	if (suffix[strlen(suffix)-1] != ']')
+	  error->all("Illegal fix ave/time command");
+	argindex[nvalues] = atoi(ptr+1);
+	*ptr = '\0';
+      } else argindex[nvalues] = 0;
+
+      n = strlen(suffix) + 1;
+      ids[nvalues] = new char[n];
+      strcpy(ids[nvalues],suffix);
+      nvalues++;
+      delete [] suffix;
+
+    } else break;
+
+    iarg++;
+  }
 
   // option defaults
 
-  sflag = 1;
-  vflag = 0;
   fp = NULL;
   ave = ONE;
 
   // optional args
 
-  int iarg = 8;
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"type") == 0) {
+    if (strcmp(arg[iarg],"file") == 0) {
-      if (iarg+2 > narg) error->all("Illegal fix ave/time command");
-      if (strcmp(arg[iarg+1],"scalar") == 0) {
-	sflag = 1;
-	vflag = 0;
-      } else if (strcmp(arg[iarg+1],"vector") == 0) {
-	sflag = 0;
-	vflag = 1;
-      } else if (strcmp(arg[iarg+1],"both") == 0) sflag = vflag = 1;
-      else error->all("Illegal fix ave/time command");
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"file") == 0) {
       if (iarg+2 > narg) error->all("Illegal fix ave/time command");
       if (me == 0) {
 	fp = fopen(arg[iarg+1],"w");
@@ -107,85 +129,108 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
     error->all("Illegal fix ave/time command");
 
-  int icompute,ifix;
+  for (int i = 0; i < nvalues; i++) {
-  if (which == COMPUTE) {
+    if (which[i] == COMPUTE) {
-    icompute = modify->find_compute(id);
+      int icompute = modify->find_compute(ids[i]);
-    if (icompute < 0) error->all("Compute ID for fix ave/time does not exist");
+      if (icompute < 0)
-  } else {
+	error->all("Compute ID for fix ave/time does not exist");
-    ifix = modify->find_fix(id);
+      if (argindex[i] == 0 && modify->compute[icompute]->scalar_flag == 0)
-    if (ifix < 0) error->all("Fix ID for fix ave/time does not exist");
+	error->all("Fix ave/time compute does not calculate a scalar");
+      if (argindex[i] && modify->compute[icompute]->vector_flag == 0)
+	error->all("Fix ave/time compute does not calculate a vector");
+      if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector)
+	error->all("Fix ave/time compute vector is accessed out-of-range");
+    } else if (which[i] == FIX) {
+      int ifix = modify->find_fix(ids[i]);
+      if (ifix < 0)
+	error->all("Fix ID for fix ave/time does not exist");
+      if (argindex[i] && modify->fix[ifix]->scalar_flag == 0)
+	error->all("Fix ave/time fix does not calculate a scalar");
+      if (argindex[i] && modify->fix[ifix]->vector_flag == 0)
+	error->all("Fix ave/time fix does not calculate a vector");
+      if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector)
+	error->all("Fix ave/time fix vector is accessed out-of-range");
+    } else if (which[i] == VARIABLE) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0)
+	error->all("Variable name for fix ave/time does not exist");
+      if (input->variable->equalstyle(ivariable) == 0)
+	error->all("Fix ave/time variable is not equal-style variable");
+    }
   }
 
-  if (which == COMPUTE) {
-    if (sflag && modify->compute[icompute]->scalar_flag == 0)
-      error->all("Fix ave/time compute does not calculate a scalar");
-    if (vflag && modify->compute[icompute]->vector_flag == 0)
-      error->all("Fix ave/time compute does not calculate a vector");
-  } else {
-    if (sflag && modify->fix[ifix]->scalar_flag == 0)
-      error->all("Fix ave/time fix does not calculate a scalar");
-    if (vflag && modify->fix[ifix]->vector_flag == 0)
-      error->all("Fix ave/time fix does not calculate a vector");
-  }
-
-  // setup list of computes to call, including pre-computes
-
-  compute = NULL;
-  if (which == COMPUTE) {
-    ncompute = 1 + modify->compute[icompute]->npre;
-    compute = new Compute*[ncompute];
-  } else ncompute = 0;
-
   // print header into file
 
   if (fp && me == 0) {
-    if (which == COMPUTE)
+    fprintf(fp,"Time-averaged data for fix %s\n",id);
-      fprintf(fp,"Time-averaged data for fix %s, group %s, and compute %s\n",
+    fprintf(fp,"TimeStep");
-	      id,group->names[modify->compute[icompute]->igroup],id);
+    for (int i = 0; i < nvalues; i++)
-    else
+      if (which[i] == COMPUTE) fprintf(fp," c_%s",ids[i]);
-      fprintf(fp,"Time-averaged data for fix %s, group %s, and fix %s\n",
+      else if (which[i] == FIX) fprintf(fp," f_%s",ids[i]);
-	      id,group->names[modify->fix[ifix]->igroup],id);
+      else if (which[i] == VARIABLE) fprintf(fp," v_%s",ids[i]);
-    if (sflag and !vflag)
+    fprintf(fp,"\n");
-      fprintf(fp,"TimeStep Value\n");
-    else if (!sflag and vflag)
-      fprintf(fp,"TimeStep Vector-values\n");
-    else if (!sflag and vflag)
-      fprintf(fp,"TimeStep Scalar-value Vector-values\n");
   }
 
-  // allocate memory for averaging
+  // allocate and initialize memory for averaging
 
-  vector = vector_total = NULL;
+  vector = new double[nvalues];
-  if (vflag) {
+  vector_total = new double[nvalues];
-    if (which == COMPUTE) size_vector = modify->compute[icompute]->size_vector;
+  for (int i = 0; i < nvalues; i++) vector_total[i] = 0.0;
-    else size_vector = modify->fix[ifix]->size_vector;
-    vector = new double[size_vector];
-    vector_total = new double[size_vector];
-  }
 
-  scalar_list = NULL;
   vector_list = NULL;
-  if (sflag && ave == WINDOW) scalar_list = new double[nwindow];
+  if (ave == WINDOW)
-  if (vflag && ave == WINDOW)
+    vector_list = memory->create_2d_double_array(nwindow,nvalues,
-    vector_list = memory->create_2d_double_array(nwindow,size_vector,
 						 "ave/time:vector_list");
 
-  // enable this fix to produce a global scalar and/or vector
+  // this fix produces either a global scalar or vector
+  // intensive/extensive flags set by compute,fix,variable that produces value
 
-  if (sflag) scalar_flag = 1;
-  if (vflag) vector_flag = 1;
   scalar_vector_freq = nfreq;
-  if (which == COMPUTE) extensive = modify->compute[icompute]->extensive;
+  extlist = NULL;
-  else extensive = modify->fix[ifix]->extensive;
+
+  if (nvalues == 1) {
+    scalar_flag = 1;
+    if (which[0] == COMPUTE) {
+      Compute *compute = modify->compute[modify->find_compute(ids[0])];
+      if (argindex[0] == 0) extscalar = compute->extscalar;
+      else if (compute->extvector >= 0) extscalar = compute->extvector;
+      else extscalar = compute->extlist[argindex[0]-1];
+    } else if (which[0] == FIX) {
+      Fix *fix = modify->fix[modify->find_fix(ids[0])];
+      if (argindex[0] == 0) extscalar = fix->extscalar;
+      else if (fix->extvector >= 0) extscalar = fix->extvector;
+      else extscalar = fix->extlist[argindex[0]-1];
+    } else if (which[0] == VARIABLE)
+      extscalar = 0;
+
+  } else {
+    vector_flag = 1;
+    size_vector = nvalues;
+    extvector = -1;
+    extlist = new int[nvalues];
+    for (int i = 0; i < nvalues; i++) {
+      if (which[i] == COMPUTE) {
+	Compute *compute = modify->compute[modify->find_compute(ids[i])];
+	if (argindex[i] == 0) extlist[i] = compute->extscalar;
+	else if (compute->extvector >= 0) extlist[i] = compute->extvector;
+	else extlist[i] = compute->extlist[argindex[i]-1];
+      } else if (which[i] == FIX) {
+	Fix *fix = modify->fix[modify->find_fix(ids[i])];
+	if (argindex[i] == 0) extlist[i] = fix->extscalar;
+	else if (fix->extvector >= 0) extlist[i] = fix->extvector;
+	else extlist[i] = fix->extlist[argindex[i]-1];
+      } else if (which[i] == VARIABLE)
+	extlist[i] = 0;
+    }
+  }
 
   // initializations
-  // set scalar and vector total to zero since they accumulate
+  // set vector total to zero since it accumulates
 
   irepeat = 0;
   iwindow = window_limit = 0;
   norm = 0;
-  scalar_total = 0.0;
+  for (int i = 0; i < nvalues; i++) vector_total[i] = 0.0;
-  if (vflag) for (int i = 0; i < size_vector; i++) vector_total[i] = 0.0;
 
   // nvalid = next step on which end_of_step does something
   // can be this timestep if multiple of nfreq and nrepeat = 1
@@ -198,24 +243,27 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
     nvalid -= (nrepeat-1)*nevery;
   if (nvalid < update->ntimestep) nvalid += nfreq;
 
-  // set timestep for all computes that store invocation times
+  // add nvalid to ALL computes that store invocation times
-  //   since don't know a priori which are invoked by this fix
+  // since don't know a priori which are invoked by this fix
   // once in end_of_step() can just set timestep for ones actually invoked
 
-  for (int i = 0; i < modify->ncompute; i++)
+  modify->addstep_compute_all(nvalid);
-    if (modify->compute[i]->timeflag) modify->compute[i]->add_step(nvalid);
 }
 
 /* ---------------------------------------------------------------------- */
 
 FixAveTime::~FixAveTime()
 {
-  delete [] id;
+  delete [] which;
+  delete [] argindex;
+  for (int i = 0; i < nvalues; i++) delete [] ids[i];
+  delete [] ids;
+  delete [] value2index;
+  delete [] extlist;
+
   if (fp && me == 0) fclose(fp);
-  delete [] compute;
   delete [] vector;
   delete [] vector_total;
-  delete [] scalar_list;
   memory->destroy_2d_double_array(vector_list);
 }
 
@@ -232,36 +280,33 @@ int FixAveTime::setmask()
 
 void FixAveTime::init()
 {
-  // set ptrs to compute and its pre-computes called each end-of-step
+  // set indices and check validity of all computes,fixes,variables
-  // put pre-computes in list before compute
-
-  if (which == COMPUTE) {
-    int icompute = modify->find_compute(id);
-    if (icompute < 0)
-      error->all("Compute ID for fix ave/time does not exist");
-  
-    ncompute = 0;
-    if (modify->compute[icompute]->npre)
-      for (int i = 0; i < modify->compute[icompute]->npre; i++) {
-	int ic = modify->find_compute(modify->compute[icompute]->id_pre[i]);
-	if (ic < 0)
-	  error->all("Precompute ID for fix ave/time does not exist");
-	compute[ncompute++] = modify->compute[ic];
-      }
-    
-    compute[ncompute++] = modify->compute[icompute];
-  }
-
-  // set ptr to fix ID
   // check that fix frequency is acceptable
 
-  if (which == FIX) {
+  int ilist = 0;
-    int ifix = modify->find_fix(id);
+
-    if (ifix < 0) 
+  for (int i = 0; i < nvalues; i++) {
-      error->all("Fix ID for fix ave/time does not exist");
+    if (which[i] == COMPUTE) {
-    fix = modify->fix[ifix];
+      int icompute = modify->find_compute(ids[i]);
-    if (nevery % fix->scalar_vector_freq)
+      if (icompute < 0)
-      error->all("Fix ave/time and fix not computed at compatible times");
+	error->all("Compute ID for fix ave/time does not exist");
+      value2index[i] = icompute;
+
+    } else if (which[i] == FIX) {
+      int ifix = modify->find_fix(ids[i]);
+      if (ifix < 0) 
+	error->all("Fix ID for fix ave/time does not exist");
+      value2index[i] = ifix;
+
+      if (nevery % modify->fix[ifix]->scalar_vector_freq)
+	error->all("Fix for fix ave/time not computed at compatible time");
+
+    } else if (which[i] == VARIABLE) {
+      int ivariable = input->variable->find(ids[i]);
+      if (ivariable < 0) 
+	error->all("Variable name for fix ave/time does not exist");
+      value2index[i] = ivariable;
+    }
   }
 }
 
@@ -278,7 +323,8 @@ void FixAveTime::setup()
 
 void FixAveTime::end_of_step()
 {
-  int i;
+  int i,j,k,m;
+  double tmp;
 
   // skip if not step which requires doing something
 
@@ -286,33 +332,44 @@ void FixAveTime::end_of_step()
 
   // zero if first step
 
-  if (irepeat == 0) {
+  if (irepeat == 0)
-    scalar = 0.0;
+    for (i = 0; i < nvalues; i++) vector[i] = 0.0;
-    if (vflag) 
-      for (i = 0; i < size_vector; i++) vector[i] = 0.0;
-  }
 
-  // accumulate results of compute or fix to local copy
+  // accumulate results of computes,fixes,variables to local copy
-  
+  // compute/fix/variable may invoke computes so wrap with clear/add
-  if (which == COMPUTE) {
-    modify->clearstep_compute();
 
-    if (sflag) {
+  modify->clearstep_compute();
-      double value;
-      for (i = 0; i < ncompute; i++) value = compute[i]->compute_scalar();
-      scalar += value;
-    }
-    if (vflag) {
-      for (i = 0; i < ncompute; i++) compute[i]->compute_vector();
-      double *cvector = compute[ncompute-1]->vector;
-      for (i = 0; i < size_vector; i++) vector[i] += cvector[i];
-    }
 
-  } else {
+  int ilist = 0;
-    if (sflag) scalar += fix->compute_scalar();
+
-    if (vflag)
+  for (i = 0; i < nvalues; i++) {
-      for (i = 0; i < size_vector; i++)
+    m = value2index[i];
-	vector[i] += fix->compute_vector(i);
+
+    // invoke compute if not previously invoked
+
+    if (which[i] == COMPUTE) {
+      Compute *compute = modify->compute[m];
+
+      if (argindex[i] == 0) {
+	if (compute->invoked & INVOKED_SCALAR) vector[i] += compute->scalar;
+	else vector[i] += compute->compute_scalar();
+      } else {
+	if (!(compute->invoked & INVOKED_VECTOR)) compute->compute_vector();
+	vector[i] += compute->vector[argindex[i]-1];
+      }
+
+    // access fix fields, guaranteed to be ready
+
+    } else if (which[i] == FIX) {
+      if (argindex[i] == 0) 
+	vector[i] += modify->fix[m]->compute_scalar();
+      else
+	vector[i] += modify->fix[m]->compute_vector(argindex[i]-1);
+
+    // evaluate equal-style variable
+
+    } else if (which[i] == VARIABLE)
+      vector[i] += input->variable->compute_equal(m);
   }
 
   // done if irepeat < nrepeat
@@ -321,46 +378,36 @@ void FixAveTime::end_of_step()
   irepeat++;
   if (irepeat < nrepeat) {
     nvalid += nevery;
-    if (which == COMPUTE) modify->addstep_compute(nvalid);
+    modify->addstep_compute(nvalid);
     return;
   }
 
   irepeat = 0;
   nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery;
-  if (which == COMPUTE) modify->addstep_compute(nvalid);
+  modify->addstep_compute(nvalid);
 
   // average the final result for the Nfreq timestep
 
   double repeat = nrepeat;
-  if (sflag) scalar /= repeat;
+  for (i = 0; i < nvalues; i++) vector[i] /= repeat;
-  if (vflag) for (i = 0; i < size_vector; i++) vector[i] /= repeat;
 
   // if ave = ONE, only single Nfreq timestep value is needed
   // if ave = RUNNING, combine with all previous Nfreq timestep values
   // if ave = WINDOW, comine with nwindow most recent Nfreq timestep values
 
   if (ave == ONE) {
-    if (sflag) scalar_total = scalar;
+    for (i = 0; i < nvalues; i++) vector_total[i] = vector[i];
-    if (vflag) for (i = 0; i < size_vector; i++) vector_total[i] = vector[i];
     norm = 1;
 
   } else if (ave == RUNNING) {
-    if (sflag) scalar_total += scalar;
+    for (i = 0; i < nvalues; i++) vector_total[i] += vector[i];
-    if (vflag) for (i = 0; i < size_vector; i++) vector_total[i] += vector[i];
     norm++;
 
   } else if (ave == WINDOW) {
-    if (sflag) {
+    for (i = 0; i < nvalues; i++) {
-      scalar_total += scalar;
+      vector_total[i] += vector[i];
-      if (window_limit) scalar_total -= scalar_list[iwindow];
+      if (window_limit) vector_total[i] -= vector_list[iwindow][i];
-      scalar_list[iwindow] = scalar;
+      vector_list[iwindow][i] = vector[i];
-    }
-    if (vflag) {
-      for (i = 0; i < size_vector; i++) {
-	vector_total[i] += vector[i];
-	if (window_limit) vector_total[i] -= vector_list[iwindow][i];
-	vector_list[iwindow][i] = vector[i];
-      }
     }
 
     iwindow++;
@@ -376,9 +423,7 @@ void FixAveTime::end_of_step()
 
   if (fp && me == 0) {
     fprintf(fp,"%d",update->ntimestep);
-    if (sflag) fprintf(fp," %g",scalar_total/norm);
+    for (i = 0; i < nvalues; i++) fprintf(fp," %g",vector_total[i]/norm);
-    if (vflag)
-      for (i = 0; i < size_vector; i++) fprintf(fp," %g",vector_total[i]/norm);
     fprintf(fp,"\n");
     fflush(fp);
   }
@@ -391,7 +436,7 @@ void FixAveTime::end_of_step()
 
 double FixAveTime::compute_scalar()
 {
-  if (norm) return scalar_total/norm;
+  if (norm) return vector_total[0]/norm;
   return 0.0;
 }
 

src/fix_ave_time.h

@@ -31,20 +31,16 @@ class FixAveTime : public Fix {
   double compute_vector(int);
 
  private:
-  int me;
+  int me,nvalues;
-  int nrepeat,nfreq,nvalid,irepeat,which,ifix;
+  int nrepeat,nfreq,nvalid,irepeat;
-  char *id;
+  int *which,*argindex,*value2index;
+  char **ids;
   FILE *fp;
 
-  int sflag,vflag,ave,nwindow,nsum;
+  int ave,nwindow,nsum;
-  double scalar,*vector;
+  double *vector;
-  int ncompute;
-  class Compute **compute;
-  class Fix *fix;
 
   int norm,iwindow,window_limit;
-  double scalar_total;
-  double *scalar_list;
   double *vector_total;
   double **vector_list;
 };

src/fix_drag.cpp

@@ -33,7 +33,7 @@ FixDrag::FixDrag(LAMMPS *lmp, int narg, char **arg) :
   vector_flag = 1;
   size_vector = 3;
   scalar_vector_freq = 1;
-  extensive = 1;
+  extvector = 1;
 
   xflag = yflag = zflag = 1;
 
@@ -46,6 +46,9 @@ FixDrag::FixDrag(LAMMPS *lmp, int narg, char **arg) :
 
   f_mag = atof(arg[6]);
   delta = atof(arg[7]);
+
+  force_flag = 0;
+  ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_dt_reset.cpp

@@ -47,7 +47,8 @@ FixDtReset::FixDtReset(LAMMPS *lmp, int narg, char **arg) :
   vector_flag = 1;
   size_vector = 1;
   scalar_vector_freq = 1;
-  extensive = 0;
+  extscalar = 0;
+  extvector = 0;
 
   nevery = atoi(arg[3]);
   if (nevery <= 0) error->all("Illegal fix dt/reset command");

src/fix_indent.cpp

@@ -42,7 +42,8 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
   vector_flag = 1;
   size_vector = 3;
   scalar_vector_freq = 1;
-  extensive = 1;
+  extscalar = 1;
+  extvector = 1;
 
   k = atof(arg[3]);
 
@@ -103,6 +104,9 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
       r0_start *= zscale;
     }
   } else error->all("Illegal fix indent command");
+
+  indenter_flag = 0;
+  indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_nph.cpp

@@ -48,7 +48,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
   box_change = 1;
   scalar_flag = 1;
   scalar_vector_freq = 1;
-  extensive = 1;
+  extscalar = 1;
 
   double p_period[3];
   if (strcmp(arg[3],"xyz") == 0) {
@@ -242,7 +242,6 @@ void FixNPH::init()
     }
 
   // set temperature and pressure ptrs
-  // set ptemperature only if pressure's id_pre[0] is not id_temp
 
   int icompute = modify->find_compute(id_temp);
   if (icompute < 0) error->all("Temp ID for fix nph does not exist");
@@ -252,14 +251,6 @@ void FixNPH::init()
   if (icompute < 0) error->all("Press ID for fix nph does not exist");
   pressure = modify->compute[icompute];
 
-  if (strcmp(id_temp,pressure->id_pre[0]) == 0) ptemperature = NULL;
-  else {
-    icompute = modify->find_compute(pressure->id_pre[0]);
-    if (icompute < 0)
-      error->all("Temp ID of press ID for fix nph does not exist");
-    ptemperature = modify->compute[icompute];
-  }
-
   // set timesteps and frequencies
   // Nkt = initial value for piston mass and energy conservation
   // guesstimate a unit-dependent t_initial if actual T = 0.0
@@ -323,18 +314,16 @@ void FixNPH::setup()
 
   double tmp = temperature->compute_scalar();
   if (press_couple == 0) {
-    if (ptemperature) tmp = ptemperature->compute_scalar();
     tmp = pressure->compute_scalar();
   } else {
     temperature->compute_vector();
-    if (ptemperature) ptemperature->compute_vector();
     pressure->compute_vector();
   }
   couple();
 
   // trigger virial computation on next timestep
 
-  pressure->add_step(update->ntimestep+1);
+  pressure->addstep(update->ntimestep+1);
 }
 
 /* ----------------------------------------------------------------------
@@ -438,18 +427,16 @@ void FixNPH::final_integrate()
 
   double tmp = temperature->compute_scalar();
   if (press_couple == 0) {
-    if (ptemperature) tmp = ptemperature->compute_scalar();
     tmp = pressure->compute_scalar();
   } else {
     temperature->compute_vector();
-    if (ptemperature) ptemperature->compute_vector();
     pressure->compute_vector();
   }
   couple();
 
   // trigger virial computation on next timestep
 
-  pressure->add_step(update->ntimestep+1);
+  pressure->addstep(update->ntimestep+1);
 
   // update omega_dot
   // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
@@ -752,13 +739,13 @@ int FixNPH::modify_param(int narg, char **arg)
     if (temperature->igroup != 0 && comm->me == 0)
       error->warning("Temperature for NPH is not for group all");
 
-    // reset id_pre[0] of pressure to new temp ID
+    // reset id_pre of pressure to new temp ID
     
     icompute = modify->find_compute(id_press);
     if (icompute < 0) error->all("Press ID for fix npt does not exist");
-    delete [] modify->compute[icompute]->id_pre[0];
+    delete [] modify->compute[icompute]->id_pre;
-    modify->compute[icompute]->id_pre[0] = new char[n];
+    modify->compute[icompute]->id_pre = new char[n];
-    strcpy(modify->compute[icompute]->id_pre[0],id_temp);
+    strcpy(modify->compute[icompute]->id_pre,id_temp);
 
     return 2;
 

src/fix_nph.h

@@ -56,7 +56,7 @@ class FixNPH : public Fix {
   double *step_respa;
 
   char *id_temp,*id_press;
-  class Compute *temperature,*pressure,*ptemperature;
+  class Compute *temperature,*pressure;
   int tflag,pflag;
 
   void couple();

src/fix_npt.cpp

@@ -47,7 +47,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
   box_change = 1;
   scalar_flag = 1;
   scalar_vector_freq = 1;
-  extensive = 1;
+  extscalar = 1;
 
   t_start = atof(arg[3]);
   t_stop = atof(arg[4]);
@@ -252,7 +252,6 @@ void FixNPT::init()
     }
 
   // set temperature and pressure ptrs
-  // set ptemperature only if pressure's id_pre[0] is not id_temp
 
   int icompute = modify->find_compute(id_temp);
   if (icompute < 0) error->all("Temp ID for fix npt does not exist");
@@ -262,14 +261,6 @@ void FixNPT::init()
   if (icompute < 0) error->all("Press ID for fix npt does not exist");
   pressure = modify->compute[icompute];
 
-  if (strcmp(id_temp,pressure->id_pre[0]) == 0) ptemperature = NULL;
-  else {
-    icompute = modify->find_compute(pressure->id_pre[0]);
-    if (icompute < 0)
-      error->all("Temp ID of press ID for fix npt does not exist");
-    ptemperature = modify->compute[icompute];
-  }
-
   // set timesteps and frequencies
 
   dtv = update->dt;
@@ -325,18 +316,16 @@ void FixNPT::setup()
 
   t_current = temperature->compute_scalar();
   if (press_couple == 0) {
-    if (ptemperature) double ptmp = ptemperature->compute_scalar();
     double tmp = pressure->compute_scalar();
   } else {
     temperature->compute_vector();
-    if (ptemperature) ptemperature->compute_vector();
     pressure->compute_vector();
   }
   couple();
 
   // trigger virial computation on next timestep
 
-  pressure->add_step(update->ntimestep+1);
+  pressure->addstep(update->ntimestep+1);
 }
 
 /* ----------------------------------------------------------------------
@@ -448,18 +437,16 @@ void FixNPT::final_integrate()
 
   t_current = temperature->compute_scalar();
   if (press_couple == 0) {
-    if (ptemperature) double ptmp = ptemperature->compute_scalar();
     double tmp = pressure->compute_scalar();
   } else {
     temperature->compute_vector();
-    if (ptemperature) ptemperature->compute_vector();
     pressure->compute_vector();
   }
   couple();
 
   // trigger virial computation on next timestep
 
-  pressure->add_step(update->ntimestep+1);
+  pressure->addstep(update->ntimestep+1);
 
   // update eta_dot
 
@@ -781,13 +768,13 @@ int FixNPT::modify_param(int narg, char **arg)
     if (temperature->igroup != 0 && comm->me == 0)
       error->warning("Temperature for NPT is not for group all");
 
-    // reset id_pre[0] of pressure to new temp ID
+    // reset id_pre of pressure to new temp ID
     
     icompute = modify->find_compute(id_press);
     if (icompute < 0) error->all("Press ID for fix npt does not exist");
-    delete [] modify->compute[icompute]->id_pre[0];
+    delete [] modify->compute[icompute]->id_pre;
-    modify->compute[icompute]->id_pre[0] = new char[n];
+    modify->compute[icompute]->id_pre = new char[n];
-    strcpy(modify->compute[icompute]->id_pre[0],id_temp);
+    strcpy(modify->compute[icompute]->id_pre,id_temp);
 
     return 2;
 

src/fix_npt.h

@@ -62,7 +62,7 @@ class FixNPT : public Fix {
   double *step_respa;
 
   char *id_temp,*id_press;
-  class Compute *temperature,*pressure,*ptemperature;
+  class Compute *temperature,*pressure;
   int tflag,pflag;
 
   void couple();

src/fix_nvt.cpp

@@ -41,7 +41,7 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) :
   restart_global = 1;
   scalar_flag = 1;
   scalar_vector_freq = 1;
-  extensive = 1;
+  extscalar = 1;
 
   t_start = atof(arg[3]);
   t_stop = atof(arg[4]);

src/fix_orient_fcc.cpp

@@ -48,7 +48,7 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = 1;
   scalar_vector_freq = 1;
-  extensive = 1;
+  extscalar = 1;
 
   nstats = atoi(arg[3]);
   direction_of_motion = atoi(arg[4]);
@@ -149,6 +149,8 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
 
   if (use_xismooth) comm_forward = 62;
   else comm_forward = 50;
+
+  added_energy = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_print.cpp

@@ -101,7 +101,7 @@ void FixPrint::end_of_step()
   // make a copy of string to work on
   // substitute for $ variables (no printing)
   // append a newline and print final copy
-  // variable evaluation may invoke a compute that affects Verlet::eflag,vflag
+  // variable evaluation may invoke computes so wrap with clear/add
 
   modify->clearstep_compute();
 

src/fix_set_force.cpp

@@ -31,7 +31,7 @@ FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) :
   vector_flag = 1;
   size_vector = 3;
   scalar_vector_freq = 1;
-  extensive = 1;
+  extvector = 1;
 
   flagx = flagy = flagz = 1;
   if (strcmp(arg[3],"NULL") == 0) flagx = 0;
@@ -40,6 +40,9 @@ FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) :
   else yvalue = atof(arg[4]);
   if (strcmp(arg[5],"NULL") == 0) flagz = 0;
   else zvalue = atof(arg[5]);
+
+  force_flag = 0;
+  foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_spring.cpp

@@ -41,7 +41,7 @@ FixSpring::FixSpring(LAMMPS *lmp, int narg, char **arg) :
   vector_flag = 1;
   size_vector = 3;
   scalar_vector_freq = 1;
-  extensive = 1;
+  extvector = 1;
 
   if (strcmp(arg[3],"tether") == 0) {
     if (narg != 9) error->all("Illegal fix spring command");
@@ -78,6 +78,9 @@ FixSpring::FixSpring(LAMMPS *lmp, int narg, char **arg) :
     if (r0 < 0) error->all("R0 < 0 for fix spring command");
 
   } else error->all("Illegal fix spring command");
+
+  force_flag = 0;
+  ftotal[0] = ftotal[1] = ftotal[2] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_temp_rescale.cpp

@@ -42,7 +42,7 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
 
   scalar_flag = 1;
   scalar_vector_freq = nevery;
-  extensive = 1;
+  extscalar = 1;
 
   t_start = atof(arg[4]);
   t_end = atof(arg[5]);
@@ -109,6 +109,8 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
   }
   delete [] newarg;
   tflag = 1;
+
+  energy = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */