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

src/ASPHERE/fix_npt_asphere.cpp

@@ -43,7 +43,7 @@ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) :
    1st half of Verlet update 
 ------------------------------------------------------------------------- */
 
-void FixNPTASphere::initial_integrate()
+void FixNPTASphere::initial_integrate(int vflag)
 {
   int i;
 

src/ASPHERE/fix_npt_asphere.h

@@ -22,7 +22,7 @@ class FixNPTASphere : public FixNPT {
  public:
   FixNPTASphere(class LAMMPS *, int, char **);
   ~FixNPTASphere() {}
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
 
  private:

src/ASPHERE/fix_nve_asphere.cpp

@@ -60,7 +60,7 @@ void FixNVEASphere::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVEASphere::initial_integrate()
+void FixNVEASphere::initial_integrate(int vflag)
 {
   double dtfm;
 

src/ASPHERE/fix_nve_asphere.h

@@ -23,7 +23,7 @@ class FixNVEASphere : public FixNVE {
   FixNVEASphere(class LAMMPS *, int, char **);
   ~FixNVEASphere();
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
 
  private:

src/ASPHERE/fix_nvt_asphere.cpp

@@ -43,7 +43,7 @@ FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) :
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVTASphere::initial_integrate()
+void FixNVTASphere::initial_integrate(int vflag)
 {
   double dtfm;
 

src/ASPHERE/fix_nvt_asphere.h

@@ -22,7 +22,7 @@ class FixNVTASphere : public FixNVT {
  public:
   FixNVTASphere(class LAMMPS *, int, char **);
   ~FixNVTASphere() {}
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
 
  private:

src/DIPOLE/fix_nve_dipole.cpp

@@ -79,7 +79,7 @@ void FixNVEDipole::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVEDipole::initial_integrate()
+void FixNVEDipole::initial_integrate(int vflag)
 {
   double dtfm,msq,scale;
 
@@ -169,14 +169,14 @@ void FixNVEDipole::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVEDipole::initial_integrate_respa(int ilevel, int flag)
+void FixNVEDipole::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   if (flag) return;             // only used by NPT,NPH
 
   dtv = step_respa[ilevel];
   dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
 
-  if (ilevel == 0) initial_integrate();
+  if (ilevel == 0) initial_integrate(vflag);
   else final_integrate();
 }
 

src/DIPOLE/fix_nve_dipole.h

@@ -24,9 +24,9 @@ class FixNVEDipole : public Fix {
   ~FixNVEDipole();
   int setmask();
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
-  void initial_integrate_respa(int,int);
+  void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int);
   void reset_dt();
 

src/GRANULAR/fix_freeze.cpp

@@ -55,9 +55,9 @@ void FixFreeze::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixFreeze::setup()
+void FixFreeze::setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/GRANULAR/fix_freeze.h

@@ -23,7 +23,7 @@ class FixFreeze : public Fix {
   FixFreeze(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
 };
 

src/GRANULAR/fix_nve_gran.cpp

@@ -61,7 +61,7 @@ void FixNVEGran::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVEGran::initial_integrate()
+void FixNVEGran::initial_integrate(int vflag)
 {
   double dtfm;
 

src/GRANULAR/fix_nve_gran.h

@@ -23,7 +23,7 @@ class FixNVEGran : public Fix {
   FixNVEGran(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
   void reset_dt();
 

src/GRANULAR/fix_wall_gran.cpp

@@ -196,9 +196,9 @@ void FixWallGran::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallGran::setup()
+void FixWallGran::setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/GRANULAR/fix_wall_gran.h

@@ -24,7 +24,7 @@ class FixWallGran : public Fix {
   ~FixWallGran();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
 
   double memory_usage();

src/compute_stress_atom.cpp

@@ -22,6 +22,8 @@
 #include "angle.h"
 #include "dihedral.h"
 #include "improper.h"
+#include "modify.h"
+#include "fix.h"
 #include "memory.h"
 #include "error.h"
 
@@ -46,10 +48,12 @@ ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) :
     keflag = 1;
     pairflag = 1;
     bondflag = angleflag = dihedralflag = improperflag = 1;
+    fixflag = 1;
   } else {
     keflag = 0;
     pairflag = 0;
     bondflag = angleflag = dihedralflag = improperflag = 0;
+    fixflag = 0;
     int iarg = 3;
     while (iarg < narg) {
       if (strcmp(arg[iarg],"ke") == 0) keflag = 1;
@@ -58,6 +62,7 @@ ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) :
       else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
       else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;
       else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;
+      else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1;
       else error->all("Illegal compute stress/atom command");
       iarg++;
     }
@@ -148,7 +153,19 @@ void ComputeStressAtom::compute_peratom()
 	stress[i][j] += vatom[i][j];
   }
 
-  // communicate ghost energy between neighbor procs
+  // add in per-atom contributions from relevant fixes
+
+  if (fixflag) {
+    for (int i = 0; i < modify->nfix; i++)
+      if (modify->fix[i]->virial_flag) {
+	double **vatom = modify->fix[i]->vatom;
+	for (i = 0; i < nlocal; i++)
+	  for (j = 0; j < 6; j++)
+	    stress[i][j] += vatom[i][j];
+      }
+  }
+
+  // communicate ghost atom virials between neighbor procs
 
   if (force->newton) comm->reverse_comm_compute(this);
 

src/compute_stress_atom.h

@@ -29,7 +29,7 @@ class ComputeStressAtom : public Compute {
   double memory_usage();
 
  private:
-  int keflag,pairflag,bondflag,angleflag,dihedralflag,improperflag;
+  int keflag,pairflag,bondflag,angleflag,dihedralflag,improperflag,fixflag;
   int nmax;
   double **stress;
 };

src/fix.cpp

@@ -14,7 +14,9 @@
 #include "string.h"
 #include "ctype.h"
 #include "fix.h"
+#include "atom.h"
 #include "group.h"
+#include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
@@ -55,6 +57,9 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
 
   comm_forward = comm_reverse = 0;
 
+  maxvatom = 0;
+  vatom = NULL;
+
   // mask settings - same as in modify.cpp
 
   INITIAL_INTEGRATE = 1;
@@ -77,6 +82,7 @@ Fix::~Fix()
 {
   delete [] id;
   delete [] style;
+  memory->destroy_2d_double_array(vatom);
 }
 
 /* ----------------------------------------------------------------------
@@ -103,3 +109,78 @@ void Fix::modify_params(int narg, char **arg)
     }
   }
 }
+
+/* ----------------------------------------------------------------------
+   setup for virial computation
+   see integrate::ev_set() for values of vflag (0-6)
+------------------------------------------------------------------------- */
+
+void Fix::v_setup(int vflag)
+{
+  int i,n;
+
+  evflag = 1;
+
+  vflag_global = vflag % 4;
+  vflag_atom = vflag / 4;
+
+  // reallocate per-atom array if necessary
+  
+  if (vflag_atom && atom->nlocal > maxvatom) {
+    maxvatom = atom->nmax;
+    memory->destroy_2d_double_array(vatom);
+    vatom = memory->create_2d_double_array(maxvatom,6,"bond:vatom");
+  }
+
+  // zero accumulators
+
+  if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0;
+  if (vflag_atom) {
+    n = atom->nlocal;
+    for (i = 0; i < n; i++) {
+      vatom[i][0] = 0.0;
+      vatom[i][1] = 0.0;
+      vatom[i][2] = 0.0;
+      vatom[i][3] = 0.0;
+      vatom[i][4] = 0.0;
+      vatom[i][5] = 0.0;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   tally virial into global and per-atom accumulators
+   v = total virial for the interaction involving total atoms
+   n = # of local atoms involved, with local indices in list
+   increment global virial by n/total fraction
+   increment per-atom virial of each atom in list by 1/total fraction
+   assumes other procs will tally left-over fractions
+------------------------------------------------------------------------- */
+
+void Fix::v_tally(int n, int *list, double total, double *v)
+{
+  int m;
+
+  if (vflag_global) {
+    double fraction = n/total;
+    virial[0] += fraction*v[0];
+    virial[1] += fraction*v[1];
+    virial[2] += fraction*v[2];
+    virial[3] += fraction*v[3];
+    virial[4] += fraction*v[4];
+    virial[5] += fraction*v[5];
+  }
+
+  if (vflag_atom) {
+    double fraction = 1.0/total;
+    for (int i = 0; i < n; i++) {
+      m = list[i];
+      vatom[m][0] += fraction*v[0];
+      vatom[m][1] += fraction*v[1];
+      vatom[m][2] += fraction*v[2];
+      vatom[m][3] += fraction*v[3];
+      vatom[m][4] += fraction*v[4];
+      vatom[m][5] += fraction*v[5];
+    }
+  }
+}

src/fix.h

@@ -51,7 +51,8 @@ class Fix : protected Pointers {
   int comm_forward;              // size of forward communication (0 if none)
   int comm_reverse;              // size of reverse communication (0 if none)
 
-  double virial[6];              // fix contribution to pressure virial
+  double virial[6];              // accumlated virial
+  double **vatom;                // accumulated per-atom virial
 
   int INITIAL_INTEGRATE,PRE_EXCHANGE,PRE_NEIGHBOR;    // mask settings
   int POST_FORCE,FINAL_INTEGRATE,END_OF_STEP,THERMO_ENERGY;
@@ -66,9 +67,9 @@ class Fix : protected Pointers {
 
   virtual void init() {}
   virtual void init_list(int, class NeighList *) {}
-  virtual void setup() {}
+  virtual void setup(int) {}
-  virtual void min_setup() {}
+  virtual void min_setup(int) {}
-  virtual void initial_integrate() {}
+  virtual void initial_integrate(int) {}
   virtual void pre_exchange() {}
   virtual void pre_neighbor() {}
   virtual void post_force(int) {}
@@ -86,7 +87,7 @@ class Fix : protected Pointers {
   virtual int size_restart(int) {return 0;}
   virtual int maxsize_restart() {return 0;}
 
-  virtual void initial_integrate_respa(int, int) {}
+  virtual void initial_integrate_respa(int, int, int) {}
   virtual void post_force_respa(int, int, int) {}
   virtual void final_integrate_respa(int) {}
 
@@ -110,6 +111,14 @@ class Fix : protected Pointers {
   virtual int modify_param(int, char **) {return 0;}
 
   virtual double memory_usage() {return 0.0;}
+
+ protected:
+  int evflag;
+  int vflag_global,vflag_atom;
+  int maxvatom;
+
+  void v_setup(int);
+  void v_tally(int, int *, double, double *);
 };
 
 }

src/fix_add_force.cpp

@@ -53,22 +53,22 @@ void FixAddForce::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixAddForce::setup()
+void FixAddForce::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixAddForce::min_setup()
+void FixAddForce::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_add_force.h

@@ -23,8 +23,8 @@ class FixAddForce : public Fix {
   FixAddForce(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_ave_atom.cpp

@@ -255,7 +255,7 @@ void FixAveAtom::init()
    only does something if nvalid = current timestep
 ------------------------------------------------------------------------- */
 
-void FixAveAtom::setup()
+void FixAveAtom::setup(int vflag)
 {
   end_of_step();
 }

src/fix_ave_atom.h

@@ -25,7 +25,7 @@ class FixAveAtom : public Fix {
   ~FixAveAtom();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void end_of_step();
 
   double memory_usage();

src/fix_ave_force.cpp

@@ -69,23 +69,23 @@ void FixAveForce::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixAveForce::setup()
+void FixAveForce::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else
     for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
       ((Respa *) update->integrate)->copy_flevel_f(ilevel);
-      post_force_respa(1,ilevel,0);
+      post_force_respa(vflag,ilevel,0);
       ((Respa *) update->integrate)->copy_f_flevel(ilevel);
     }
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixAveForce::min_setup()
+void FixAveForce::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_ave_force.h

@@ -23,8 +23,8 @@ class FixAveForce : public Fix {
   FixAveForce(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_ave_spatial.cpp

@@ -403,7 +403,7 @@ void FixAveSpatial::init()
    only does something if nvalid = current timestep
 ------------------------------------------------------------------------- */
 
-void FixAveSpatial::setup()
+void FixAveSpatial::setup(int vflag)
 {
   end_of_step();
 }

src/fix_ave_spatial.h

@@ -25,7 +25,7 @@ class FixAveSpatial : public Fix {
   ~FixAveSpatial();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void end_of_step();
   double compute_vector(int);
   double memory_usage();

src/fix_ave_time.cpp

@@ -314,7 +314,7 @@ void FixAveTime::init()
    only does something if nvalid = current timestep
 ------------------------------------------------------------------------- */
 
-void FixAveTime::setup()
+void FixAveTime::setup(int vflag)
 {
   end_of_step();
 }

src/fix_ave_time.h

@@ -25,7 +25,7 @@ class FixAveTime : public Fix {
   ~FixAveTime();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void end_of_step();
   double compute_scalar();
   double compute_vector(int);

src/fix_com.cpp

@@ -72,7 +72,7 @@ void FixCOM::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixCOM::setup()
+void FixCOM::setup(int vflag)
 {
   if (first) end_of_step();
   first = 0;

src/fix_com.h

@@ -25,7 +25,7 @@ class FixCOM : public Fix {
   ~FixCOM();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void end_of_step();
 
  private:

src/fix_drag.cpp

@@ -71,13 +71,13 @@ void FixDrag::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixDrag::setup()
+void FixDrag::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }

src/fix_drag.h

@@ -23,7 +23,7 @@ class FixDrag : public Fix {
   FixDrag(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   double compute_vector(int);

src/fix_dt_reset.cpp

@@ -120,7 +120,7 @@ void FixDtReset::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixDtReset::setup()
+void FixDtReset::setup(int vflag)
 {
   end_of_step();
 }

src/fix_dt_reset.h

@@ -24,7 +24,7 @@ class FixDtReset : public Fix {
   ~FixDtReset() {}
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void end_of_step();
   double compute_scalar();
   double compute_vector(int);

src/fix_efield.cpp

@@ -65,13 +65,13 @@ void FixEfield::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixEfield::setup()
+void FixEfield::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }

src/fix_efield.h

@@ -23,7 +23,7 @@ class FixEfield : public Fix {
   FixEfield(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
 

src/fix_enforce2d.cpp

@@ -41,15 +41,15 @@ int FixEnforce2D::setmask()
 
 /* ---------------------------------------------------------------------- */
 
-void FixEnforce2D::setup()
+void FixEnforce2D::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     int nlevels_respa = ((Respa *) update->integrate)->nlevels;
     for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
       ((Respa *) update->integrate)->copy_flevel_f(ilevel);
-      post_force_respa(1,ilevel,0);
+      post_force_respa(vflag,ilevel,0);
       ((Respa *) update->integrate)->copy_f_flevel(ilevel);
     }
   }
@@ -57,9 +57,9 @@ void FixEnforce2D::setup()
 
 /* ---------------------------------------------------------------------- */
 
-void FixEnforce2D::min_setup()
+void FixEnforce2D::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_enforce2d.h

@@ -22,8 +22,8 @@ class FixEnforce2D : public Fix {
  public:
   FixEnforce2D(class LAMMPS *, int, char **);
   int setmask();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_gravity.cpp

@@ -111,9 +111,9 @@ void FixGravity::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixGravity::setup()
+void FixGravity::setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_gravity.h

@@ -25,7 +25,7 @@ class FixGravity : public Fix {
   FixGravity(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
 
  private:

src/fix_gyration.cpp

@@ -72,7 +72,7 @@ void FixGyration::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixGyration::setup()
+void FixGyration::setup(int vflag)
 {
   if (first) end_of_step();
   first = 0;

src/fix_gyration.h

@@ -25,7 +25,7 @@ class FixGyration : public Fix {
   ~FixGyration();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void end_of_step();
 
  private:

src/fix_indent.cpp

@@ -131,22 +131,22 @@ void FixIndent::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixIndent::setup()
+void FixIndent::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixIndent::min_setup()
+void FixIndent::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_indent.h

@@ -23,8 +23,8 @@ class FixIndent : public Fix {
   FixIndent(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   virtual void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_langevin.cpp

@@ -126,13 +126,13 @@ void FixLangevin::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixLangevin::setup()
+void FixLangevin::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }

src/fix_langevin.h

@@ -24,7 +24,7 @@ class FixLangevin : public Fix {
   ~FixLangevin();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void reset_target(double);

src/fix_line_force.cpp

@@ -45,15 +45,15 @@ int FixLineForce::setmask()
 
 /* ---------------------------------------------------------------------- */
 
-void FixLineForce::setup()
+void FixLineForce::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     int nlevels_respa = ((Respa *) update->integrate)->nlevels;
     for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
       ((Respa *) update->integrate)->copy_flevel_f(ilevel);
-      post_force_respa(1,ilevel,0);
+      post_force_respa(vflag,ilevel,0);
       ((Respa *) update->integrate)->copy_f_flevel(ilevel);
     }
   }
@@ -61,9 +61,9 @@ void FixLineForce::setup()
 
 /* ---------------------------------------------------------------------- */
 
-void FixLineForce::min_setup()
+void FixLineForce::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_line_force.h

@@ -22,8 +22,8 @@ class FixLineForce : public Fix {
  public:
   FixLineForce(class LAMMPS *, int, char **);
   int setmask();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_msd.cpp

@@ -123,7 +123,7 @@ void FixMSD::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixMSD::setup()
+void FixMSD::setup(int vflag)
 {
   if (first) end_of_step();
   first = 0;

src/fix_msd.h

@@ -25,7 +25,7 @@ class FixMSD : public Fix {
   ~FixMSD();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void end_of_step();
 
   double memory_usage();

src/fix_nph.cpp

@@ -306,7 +306,7 @@ void FixNPH::init()
    compute T,P before integrator starts 
 ------------------------------------------------------------------------- */
 
-void FixNPH::setup()
+void FixNPH::setup(int vflag)
 {
   p_target[0] = p_start[0];                 // used by compute_scalar()
   p_target[1] = p_start[1];
@@ -330,7 +330,7 @@ void FixNPH::setup()
    1st half of Verlet update 
 ------------------------------------------------------------------------- */
 
-void FixNPH::initial_integrate()
+void FixNPH::initial_integrate(int vflag)
 {
   int i;
 
@@ -455,7 +455,7 @@ void FixNPH::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNPH::initial_integrate_respa(int ilevel, int flag)
+void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   // if flag = 1, then is 2nd call at outermost level from rRESPA
   // perform 2nd half of box remap on own + ghost atoms and return

src/fix_nph.h

@@ -24,10 +24,10 @@ class FixNPH : public Fix {
   ~FixNPH();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
-  void initial_integrate_respa(int,int);
+  void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int);
   double compute_scalar();
   void write_restart(FILE *);

src/fix_npt.cpp

@@ -307,7 +307,7 @@ void FixNPT::init()
    compute T,P before integrator starts 
 ------------------------------------------------------------------------- */
 
-void FixNPT::setup()
+void FixNPT::setup(int vflag)
 {
   t_target = t_start;                      // used by compute_scalar()
   p_target[0] = p_start[0];
@@ -332,7 +332,7 @@ void FixNPT::setup()
    1st half of Verlet update 
 ------------------------------------------------------------------------- */
 
-void FixNPT::initial_integrate()
+void FixNPT::initial_integrate(int vflag)
 {
   int i;
 
@@ -471,7 +471,7 @@ void FixNPT::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNPT::initial_integrate_respa(int ilevel, int flag)
+void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   // if flag = 1, then is 2nd call at outermost level from rRESPA
   // perform 2nd half of box remap on own + ghost atoms and return

src/fix_npt.h

@@ -24,10 +24,10 @@ class FixNPT : public Fix {
   virtual ~FixNPT();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  virtual void initial_integrate();
+  void initial_integrate(int);
   virtual void final_integrate();
-  void initial_integrate_respa(int, int);
+  void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int);
   double compute_scalar();
   void write_restart(FILE *);

src/fix_nve.cpp

@@ -58,7 +58,7 @@ void FixNVE::init()
    allow for both per-type and per-atom mass
 ------------------------------------------------------------------------- */
 
-void FixNVE::initial_integrate()
+void FixNVE::initial_integrate(int vflag)
 {
   double dtfm;
 
@@ -137,14 +137,14 @@ void FixNVE::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVE::initial_integrate_respa(int ilevel, int flag)
+void FixNVE::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   if (flag) return;             // only used by NPT,NPH
 
   dtv = step_respa[ilevel];
   dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
 
-  if (ilevel == 0) initial_integrate();
+  if (ilevel == 0) initial_integrate(vflag);
   else final_integrate();
 }
 

src/fix_nve.h

@@ -23,9 +23,9 @@ class FixNVE : public Fix {
   FixNVE(class LAMMPS *, int, char **);
   int setmask();
   virtual void init();
-  virtual void initial_integrate();
+  virtual void initial_integrate(int);
   virtual void final_integrate();
-  void initial_integrate_respa(int, int);
+  void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int);
   void reset_dt();
 

src/fix_nve_limit.cpp

@@ -67,7 +67,7 @@ void FixNVELimit::init()
    allow for both per-type and per-atom mass
 ------------------------------------------------------------------------- */
 
-void FixNVELimit::initial_integrate()
+void FixNVELimit::initial_integrate(int vflag)
 {
   double dtfm,vsq,scale;
 
@@ -184,14 +184,14 @@ void FixNVELimit::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVELimit::initial_integrate_respa(int ilevel, int flag)
+void FixNVELimit::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   if (flag) return;             // only used by NPT,NPH
 
   dtv = step_respa[ilevel];
   dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
 
-  if (ilevel == 0) initial_integrate();
+  if (ilevel == 0) initial_integrate(vflag);
   else final_integrate();
 }
 

src/fix_nve_limit.h

@@ -23,9 +23,9 @@ class FixNVELimit : public Fix {
   FixNVELimit(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
-  void initial_integrate_respa(int, int);
+  void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int);
   void reset_dt();
   double compute_scalar();

src/fix_nve_noforce.cpp

@@ -51,7 +51,7 @@ void FixNVENoforce::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVENoforce::initial_integrate()
+void FixNVENoforce::initial_integrate(int vflag)
 {
   double **x = atom->x;
   double **v = atom->v;
@@ -69,13 +69,13 @@ void FixNVENoforce::initial_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVENoforce::initial_integrate_respa(int ilevel, int flag)
+void FixNVENoforce::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   if (flag) return;             // only used by NPT,NPH
 
   dtv = step_respa[ilevel];
 
-  if (ilevel == 0) initial_integrate();
+  if (ilevel == 0) initial_integrate(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_nve_noforce.h

@@ -23,8 +23,8 @@ class FixNVENoforce : public Fix {
   FixNVENoforce(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
-  void initial_integrate_respa(int, int);
+  void initial_integrate_respa(int, int, int);
   void reset_dt();
 
  private:

src/fix_nvt.cpp

@@ -134,7 +134,7 @@ void FixNVT::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVT::setup()
+void FixNVT::setup(int vflag)
 {
   t_target = t_start;                         // used by compute_scalar()
   t_current = temperature->compute_scalar();
@@ -142,7 +142,7 @@ void FixNVT::setup()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVT::initial_integrate()
+void FixNVT::initial_integrate(int vflag)
 {
   double dtfm;
 
@@ -218,7 +218,7 @@ void FixNVT::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVT::initial_integrate_respa(int ilevel, int flag)
+void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   if (flag) return;             // only used by NPT,NPH
 

src/fix_nvt.h

@@ -24,10 +24,10 @@ class FixNVT : public Fix {
   virtual ~FixNVT();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  virtual void initial_integrate();
+  virtual void initial_integrate(int);
   virtual void final_integrate();
-  virtual void initial_integrate_respa(int,int);
+  virtual void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int);
   double compute_scalar();
   void write_restart(FILE *);

src/fix_nvt_sllod.cpp

@@ -62,7 +62,7 @@ void FixNVTSlodd::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVTSlodd::initial_integrate()
+void FixNVTSlodd::initial_integrate(int vflag)
 {
   double dtfm;
 
@@ -191,7 +191,7 @@ void FixNVTSlodd::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixNVTSlodd::initial_integrate_respa(int ilevel, int flag)
+void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   if (flag) return;             // only used by NPT,NPH
 

src/fix_nvt_sllod.h

@@ -23,9 +23,9 @@ class FixNVTSlodd : public FixNVT {
   FixNVTSlodd(class LAMMPS *, int, char **);
   ~FixNVTSlodd() {}
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
-  void initial_integrate_respa(int,int);
+  void initial_integrate_respa(int, int, int);
 };
 
 }

src/fix_orient_fcc.cpp

@@ -198,13 +198,13 @@ void FixOrientFCC::init_list(int id, NeighList *ptr)
 
 /* ---------------------------------------------------------------------- */
 
-void FixOrientFCC::setup()
+void FixOrientFCC::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }

src/fix_orient_fcc.h

@@ -42,7 +42,7 @@ class FixOrientFCC : public Fix {
   int setmask();
   void init();
   void init_list(int, class NeighList *);
-  void setup();
+  void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   double compute_scalar();

src/fix_plane_force.cpp

@@ -45,15 +45,15 @@ int FixPlaneForce::setmask()
 
 /* ---------------------------------------------------------------------- */
 
-void FixPlaneForce::setup()
+void FixPlaneForce::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     int nlevels_respa = ((Respa *) update->integrate)->nlevels;
     for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
       ((Respa *) update->integrate)->copy_flevel_f(ilevel);
-      post_force_respa(1,ilevel,0);
+      post_force_respa(vflag,ilevel,0);
       ((Respa *) update->integrate)->copy_f_flevel(ilevel);
     }
   }
@@ -61,9 +61,9 @@ void FixPlaneForce::setup()
 
 /* ---------------------------------------------------------------------- */
 
-void FixPlaneForce::min_setup()
+void FixPlaneForce::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_plane_force.h

@@ -22,8 +22,8 @@ class FixPlaneForce : public Fix {
  public:
   FixPlaneForce(class LAMMPS *, int, char **);
   int setmask();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_rdf.cpp

@@ -161,7 +161,7 @@ void FixRDF::init_list(int id, NeighList *ptr)
 
 /* ---------------------------------------------------------------------- */
 
-void FixRDF::setup()
+void FixRDF::setup(int vflag)
 {
   if (first) end_of_step();
   first = 0;

src/fix_rdf.h

@@ -26,7 +26,7 @@ class FixRDF : public Fix {
   int setmask();
   void init();
   void init_list(int, class NeighList *);
-  void setup();
+  void setup(int);
   void end_of_step();
 
  private:

src/fix_recenter.cpp

@@ -132,7 +132,7 @@ void FixRecenter::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixRecenter::initial_integrate()
+void FixRecenter::initial_integrate(int vflag)
 {
   // target COM
   // bounding box around domain works for both orthogonal and triclinic

src/fix_recenter.h

@@ -23,7 +23,7 @@ class FixRecenter : public Fix {
   FixRecenter(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
 
  private:
   int group2bit,scaleflag;

src/fix_rigid.cpp

@@ -219,10 +219,6 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
     if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum);
     if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum);
   }
-
-  // zero fix_rigid virial in case pressure uses it before 1st fix_rigid call
-
-  for (int n = 0; n < 6; n++) virial[n] = 0.0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -560,9 +556,9 @@ void FixRigid::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixRigid::setup()
+void FixRigid::setup(int vflag)
 {
-  int i,ibody;
+  int i,n,ibody;
   
   // vcm = velocity of center-of-mass of each rigid body
   // fcm = force on center-of-mass of each rigid body
@@ -660,21 +656,32 @@ void FixRigid::setup()
     torque[ibody][2] = all[ibody][5];
   }
 
+  // virial setup before call to set_v
+
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   // set velocities from angmom & omega
-  // guestimate virial as 2x the set_v contribution
 
   for (ibody = 0; ibody < nbody; ibody++)
     omega_from_mq(angmom[ibody],ex_space[ibody],ey_space[ibody],
 		  ez_space[ibody],inertia[ibody],omega[ibody]);
+  set_v();
 
-  for (int n = 0; n < 6; n++) virial[n] = 0.0;
+  // guestimate virial as 2x the set_v contribution
-  set_v(1);
+
-  for (int n = 0; n < 6; n++) virial[n] *= 2.0;
+  if (vflag_global)
+    for (n = 0; n < 6; n++) virial[n] *= 2.0;
+  if (vflag_atom) {
+    for (i = 0; i < nlocal; i++)
+      for (n = 0; n < 6; n++)
+	vatom[i][n] *= 2.0;
+  }
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixRigid::initial_integrate()
+void FixRigid::initial_integrate(int vflag)
 {
   double dtfm;
 
@@ -708,12 +715,15 @@ void FixRigid::initial_integrate()
 	       ex_space[ibody],ey_space[ibody],ez_space[ibody]);
   }
 
+  // virial setup before call to set_xv
+
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   // set coords and velocities if atoms in rigid bodies
   // from quarternion and omega
   
-  int vflag = 0;
+  set_xv();
-  if (pressure_flag || output->next_thermo == update->ntimestep) vflag = 1;
-  set_xv(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -862,19 +872,18 @@ void FixRigid::final_integrate()
   }
 
   // set velocities from angmom & omega
-  
+  // virial is already setup from initial_integrate
+
   for (ibody = 0; ibody < nbody; ibody++) 
     omega_from_mq(angmom[ibody],ex_space[ibody],ey_space[ibody],
 		  ez_space[ibody],inertia[ibody],omega[ibody]);
 
-  int vflag = 0;
+  set_v();
-  if (pressure_flag || output->next_thermo == update->ntimestep) vflag = 1;
-  set_v(vflag);
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixRigid::initial_integrate_respa(int ilevel, int flag)
+void FixRigid::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   if (flag) return;             // only used by NPT,NPH
 
@@ -882,7 +891,7 @@ void FixRigid::initial_integrate_respa(int ilevel, int flag)
   dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
   dtq = 0.5 * step_respa[ilevel];
 
-  if (ilevel == 0) initial_integrate();
+  if (ilevel == 0) initial_integrate(vflag);
   else final_integrate();
 }
 
@@ -1220,12 +1229,13 @@ void FixRigid::exyz_from_q(double *q, double *ex, double *ey, double *ez)
    v = Vcm + (W cross (x - Xcm))
 ------------------------------------------------------------------------- */
 
-void FixRigid::set_xv(int vflag)
+void FixRigid::set_xv()
 {
   int ibody;
   int xbox,ybox,zbox;
-  double vold0,vold1,vold2,fc0,fc1,fc2,massone,x0,x1,x2;
+  double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
   double xy,xz,yz;
+  double vr[6];
 
   int *image = atom->image;
   double **x = atom->x;
@@ -1246,10 +1256,6 @@ void FixRigid::set_xv(int vflag)
   double **f = atom->f;
   int *type = atom->type;
   
-  // zero out fix_rigid virial
-
-  if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0;
-
   for (int i = 0; i < nlocal; i++) {
     if (body[i] < 0) continue;
     ibody = body[i];
@@ -1258,16 +1264,21 @@ void FixRigid::set_xv(int vflag)
     ybox = (image[i] >> 10 & 1023) - 512;
     zbox = (image[i] >> 20) - 512;
 
-    // save old positions and velocities for virial contribution
+    // save old positions and velocities for virial
 
-    if (vflag) {
+    if (evflag) {
-      x0 = x[i][0] + xbox*xprd;
+      if (triclinic == 0) {
-      x1 = x[i][1] + ybox*yprd;
+	x0 = x[i][0] + xbox*xprd;
-      x2 = x[i][2] + zbox*zprd;
+	x1 = x[i][1] + ybox*yprd;
-
+	x2 = x[i][2] + zbox*zprd;
-      vold0 = v[i][0];
+      } else {
-      vold1 = v[i][1];
+	x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
-      vold2 = v[i][2];
+	x1 = x[i][1] + ybox*yprd + zbox*yz;
+	x2 = x[i][2] + zbox*zprd;
+      }
+      v0 = v[i][0];
+      v1 = v[i][1];
+      v2 = v[i][2];
     }
 
     // x = displacement from center-of-mass, based on body orientation
@@ -1304,20 +1315,26 @@ void FixRigid::set_xv(int vflag)
       x[i][2] += xcm[ibody][2] - zbox*zprd;
     }
 
-    // compute body constraint forces for virial
+    // virial = unwrapped coords dotted into body constraint force
+    // body constraint force = implied force due to v change minus f external
+    // assume f does not include forces internal to body
+    // 1/2 factor b/c final_integrate contributes other half
+    // assume per-atom contribution is due to constraint force on that atom
 
-    if (vflag) {
+    if (evflag) {
       massone = mass[type[i]];
-      fc0 = massone*(v[i][0] - vold0)/dtf - f[i][0];
+      fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
-      fc1 = massone*(v[i][1] - vold1)/dtf - f[i][1];
+      fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
-      fc2 = massone*(v[i][2] - vold2)/dtf - f[i][2]; 
+      fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; 
 
-      virial[0] += 0.5*fc0*x0;
+      vr[0] = 0.5*x0*fc0;
-      virial[1] += 0.5*fc1*x1;
+      vr[1] = 0.5*x1*fc1;
-      virial[2] += 0.5*fc2*x2;
+      vr[2] = 0.5*x2*fc2;
-      virial[3] += 0.5*fc1*x0;
+      vr[3] = 0.5*x0*fc1;
-      virial[4] += 0.5*fc2*x0;
+      vr[4] = 0.5*x0*fc2;
-      virial[5] += 0.5*fc2*x1;
+      vr[5] = 0.5*x1*fc2;
+
+      v_tally(1,&i,1.0,vr);
     }
   }
 }
@@ -1327,22 +1344,22 @@ void FixRigid::set_xv(int vflag)
    v = Vcm + (W cross (x - Xcm))
 ------------------------------------------------------------------------- */
 
-void FixRigid::set_v(int vflag)
+void FixRigid::set_v()
 {
-  double **v = atom->v;
-  int nlocal = atom->nlocal;
-
   int ibody;
+  int xbox,ybox,zbox;
   double xunwrap,yunwrap,zunwrap,dx,dy,dz;
-  double vold0,vold1,vold2,fc0,fc1,fc2,massone;
+  double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
   double xy,xz,yz;
+  double vr[6];
 
   double *mass = atom->mass; 
   double **f = atom->f;
+  double **v = atom->v;
   double **x = atom->x;
   int *type = atom->type;
   int *image = atom->image;
-  int xbox,ybox,zbox;
+  int nlocal = atom->nlocal;
 
   double xprd = domain->xprd;
   double yprd = domain->yprd;
@@ -1369,45 +1386,50 @@ void FixRigid::set_v(int vflag)
 
     // save old velocities for virial
 
-    if (vflag) {
+    if (evflag) {
-      vold0 = v[i][0];
+      v0 = v[i][0];
-      vold1 = v[i][1];
+      v1 = v[i][1];
-      vold2 = v[i][2];
+      v2 = v[i][2];
     }
 
     v[i][0] = omega[ibody][1]*dz - omega[ibody][2]*dy + vcm[ibody][0];
     v[i][1] = omega[ibody][2]*dx - omega[ibody][0]*dz + vcm[ibody][1];
     v[i][2] = omega[ibody][0]*dy - omega[ibody][1]*dx + vcm[ibody][2];
 
-    // compute body constraint forces for virial
+    // virial = unwrapped coords dotted into body constraint force
-    // use unwrapped atom positions
+    // body constraint force = implied force due to v change minus f external
+    // assume f does not include forces internal to body
+    // 1/2 factor b/c initial_integrate contributes other half
+    // assume per-atom contribution is due to constraint force on that atom
 
-    if (vflag) {
+    if (evflag) {
       massone = mass[type[i]];
-      fc0 = massone*(v[i][0] - vold0)/dtf - f[i][0];
+      fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
-      fc1 = massone*(v[i][1] - vold1)/dtf - f[i][1];
+      fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
-      fc2 = massone*(v[i][2] - vold2)/dtf - f[i][2]; 
+      fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; 
 
       xbox = (image[i] & 1023) - 512;
       ybox = (image[i] >> 10 & 1023) - 512;
       zbox = (image[i] >> 20) - 512;
 
       if (triclinic == 0) {
-	xunwrap = x[i][0] + xbox*xprd;
+	x0 = x[i][0] + xbox*xprd;
-	yunwrap = x[i][1] + ybox*yprd;
+	x1 = x[i][1] + ybox*yprd;
-	zunwrap = x[i][2] + zbox*zprd;
+	x2 = x[i][2] + zbox*zprd;
       } else {
-	xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
+	x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
-	yunwrap = x[i][1] + ybox*yprd + zbox*yz;
+	x1 = x[i][1] + ybox*yprd + zbox*yz;
-	zunwrap = x[i][2] + zbox*zprd;
+	x2 = x[i][2] + zbox*zprd;
       }
 
-      virial[0] += 0.5*fc0*xunwrap;
+      vr[0] = 0.5*x0*fc0;
-      virial[1] += 0.5*fc1*yunwrap;
+      vr[1] = 0.5*x1*fc1;
-      virial[2] += 0.5*fc2*zunwrap;
+      vr[2] = 0.5*x2*fc2;
-      virial[3] += 0.5*fc1*xunwrap;
+      vr[3] = 0.5*x0*fc1;
-      virial[4] += 0.5*fc2*xunwrap;
+      vr[4] = 0.5*x0*fc2;
-      virial[5] += 0.5*fc2*yunwrap;
+      vr[5] = 0.5*x1*fc2;
+
+      v_tally(1,&i,1.0,vr);
     }
   }
 }
@@ -1421,6 +1443,7 @@ double FixRigid::memory_usage()
   int nmax = atom->nmax;
   double bytes = nmax * sizeof(int);
   bytes += nmax*3 * sizeof(double);
+  bytes += maxvatom*6 * sizeof(double);
   return bytes;
 } 
 

src/fix_rigid.h

@@ -24,10 +24,10 @@ class FixRigid : public Fix {
   ~FixRigid();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  void initial_integrate();
+  void initial_integrate(int);
   void final_integrate();
-  void initial_integrate_respa(int, int);
+  void initial_integrate_respa(int, int, int);
   void final_integrate_respa(int);
 
   double memory_usage();
@@ -77,8 +77,8 @@ class FixRigid : public Fix {
   void omega_from_mq(double *, double *, double *,
 		     double *, double *, double *);
   void exyz_from_q(double *, double *, double *, double *);
-  void set_xv(int);
+  void set_xv();
-  void set_v(int);
+  void set_v();
 };
 
 }

src/fix_set_force.cpp

@@ -66,23 +66,23 @@ void FixSetForce::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixSetForce::setup()
+void FixSetForce::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else
     for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) {
       ((Respa *) update->integrate)->copy_flevel_f(ilevel);
-      post_force_respa(1,ilevel,0);
+      post_force_respa(vflag,ilevel,0);
       ((Respa *) update->integrate)->copy_f_flevel(ilevel);
     }
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixSetForce::min_setup()
+void FixSetForce::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_set_force.h

@@ -23,8 +23,8 @@ class FixSetForce : public Fix {
   FixSetForce(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_shake.cpp

@@ -379,7 +379,7 @@ void FixShake::init()
    SHAKE as pre-integrator constraint 
 ------------------------------------------------------------------------- */
 
-void FixShake::setup()
+void FixShake::setup(int vflag)
 {
   pre_neighbor();
 
@@ -398,14 +398,14 @@ void FixShake::setup()
   if (strcmp(update->integrate_style,"verlet") == 0) {
     dtv = update->dt;
     dtfsq = 0.5 * update->dt * update->dt * force->ftm2v;
-    post_force(1);
+    post_force(vflag);
     dtfsq = update->dt * update->dt * force->ftm2v;
   } else {
     dtv = step_respa[0];
     dtf_innerhalf = 0.5 * step_respa[0] * force->ftm2v;
     dtf_inner = dtf_innerhalf;
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
     dtf_inner = step_respa[0] * force->ftm2v;
   }
@@ -490,7 +490,7 @@ void FixShake::pre_neighbor()
    compute the force adjustment for SHAKE constraint 
 ------------------------------------------------------------------------- */
 
-void FixShake::post_force(int vflag_in)
+void FixShake::post_force(int vflag)
 {
   if (update->ntimestep == next_output) stats();
 
@@ -502,10 +502,10 @@ void FixShake::post_force(int vflag_in)
 
   if (nprocs > 1) comm->comm_fix(this);
 
-  // zero out SHAKE contribution to virial
+  // virial setup
 
-  vflag = vflag_in;
+  if (vflag) v_setup(vflag);
-  if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0;
+  else evflag = 0;
 
   // loop over clusters
 
@@ -1148,6 +1148,9 @@ void FixShake::unconstrained_update()
 
 void FixShake::shake2(int m)
 {
+  int nlist,list[2];
+  double v[6];
+
   // local atom IDs and constraint distances
 
   int i0 = atom->map(shake_atom[m][0]);
@@ -1218,18 +1221,19 @@ void FixShake::shake2(int m)
     f[i1][2] -= lamda*r01[2];
   }
 
-  if (vflag) {
+  if (evflag) {
-    int factor = 0;
+    nlist = 0;
-    if (i0 < nlocal) factor++;
+    if (i0 < nlocal) list[nlist++] = i0;
-    if (i1 < nlocal) factor++;
+    if (i1 < nlocal) list[nlist++] = i1;
-    double rfactor = 0.5 * factor;
 
-    virial[0] += rfactor*lamda*r01[0]*r01[0];
+    v[0] = lamda*r01[0]*r01[0];
-    virial[1] += rfactor*lamda*r01[1]*r01[1];
+    v[1] = lamda*r01[1]*r01[1];
-    virial[2] += rfactor*lamda*r01[2]*r01[2];
+    v[2] = lamda*r01[2]*r01[2];
-    virial[3] += rfactor*lamda*r01[0]*r01[1];
+    v[3] = lamda*r01[0]*r01[1];
-    virial[4] += rfactor*lamda*r01[0]*r01[2];
+    v[4] = lamda*r01[0]*r01[2];
-    virial[5] += rfactor*lamda*r01[1]*r01[2];
+    v[5] = lamda*r01[1]*r01[2];
+
+    v_tally(nlist,list,2.0,v);
   }
 }
 
@@ -1237,6 +1241,9 @@ void FixShake::shake2(int m)
 
 void FixShake::shake3(int m)
 {
+  int nlist,list[3];
+  double v[6];
+
   // local atom IDs and constraint distances
 
   int i0 = atom->map(shake_atom[m][0]);
@@ -1371,26 +1378,20 @@ void FixShake::shake3(int m)
     f[i2][2] -= lamda02*r02[2];
   }
 
-  if (vflag) {
+  if (evflag) {
-    int factor = 0;
+    nlist = 0;
-    if (i0 < nlocal) factor++;
+    if (i0 < nlocal) list[nlist++] = i0;
-    if (i1 < nlocal) factor++;
+    if (i1 < nlocal) list[nlist++] = i1;
-    if (i2 < nlocal) factor++;
+    if (i2 < nlocal) list[nlist++] = i2;
-    double rfactor = factor/3.0;
 
-    virial[0] += rfactor*lamda01*r01[0]*r01[0];
+    v[0] = lamda01*r01[0]*r01[0] + lamda02*r02[0]*r02[0];
-    virial[1] += rfactor*lamda01*r01[1]*r01[1];
+    v[1] = lamda01*r01[1]*r01[1] + lamda02*r02[1]*r02[1];
-    virial[2] += rfactor*lamda01*r01[2]*r01[2];
+    v[2] = lamda01*r01[2]*r01[2] + lamda02*r02[2]*r02[2];
-    virial[3] += rfactor*lamda01*r01[0]*r01[1];
+    v[3] = lamda01*r01[0]*r01[1] + lamda02*r02[0]*r02[1];
-    virial[4] += rfactor*lamda01*r01[0]*r01[2];
+    v[4] = lamda01*r01[0]*r01[2] + lamda02*r02[0]*r02[2];
-    virial[5] += rfactor*lamda01*r01[1]*r01[2];
+    v[5] = lamda01*r01[1]*r01[2] + lamda02*r02[1]*r02[2];
 
-    virial[0] += rfactor*lamda02*r02[0]*r02[0];
+    v_tally(nlist,list,3.0,v);
-    virial[1] += rfactor*lamda02*r02[1]*r02[1];
-    virial[2] += rfactor*lamda02*r02[2]*r02[2];
-    virial[3] += rfactor*lamda02*r02[0]*r02[1];
-    virial[4] += rfactor*lamda02*r02[0]*r02[2];
-    virial[5] += rfactor*lamda02*r02[1]*r02[2];
   }
 }
 
@@ -1398,6 +1399,9 @@ void FixShake::shake3(int m)
 
 void FixShake::shake4(int m)
 {
+  int nlist,list[4];
+  double v[6];
+
   // local atom IDs and constraint distances
 
   int i0 = atom->map(shake_atom[m][0]);
@@ -1608,34 +1612,21 @@ void FixShake::shake4(int m)
     f[i3][2] -= lamda03*r03[2];
   }
 
-  if (vflag) {
+  if (evflag) {
-    int factor = 0;
+    nlist = 0;
-    if (i0 < nlocal) factor++;
+    if (i0 < nlocal) list[nlist++] = i0;
-    if (i1 < nlocal) factor++;
+    if (i1 < nlocal) list[nlist++] = i1;
-    if (i2 < nlocal) factor++;
+    if (i2 < nlocal) list[nlist++] = i2;
-    if (i3 < nlocal) factor++;
+    if (i3 < nlocal) list[nlist++] = i3;
-    double rfactor = 0.25*factor;
 
-    virial[0] += rfactor*lamda01*r01[0]*r01[0];
+    v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda03*r03[0]*r03[0];
-    virial[1] += rfactor*lamda01*r01[1]*r01[1];
+    v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda03*r03[1]*r03[1];
-    virial[2] += rfactor*lamda01*r01[2]*r01[2];
+    v[2] = lamda01*r01[2]*r01[2]+lamda02*r02[2]*r02[2]+lamda03*r03[2]*r03[2];
-    virial[3] += rfactor*lamda01*r01[0]*r01[1];
+    v[3] = lamda01*r01[0]*r01[1]+lamda02*r02[0]*r02[1]+lamda03*r03[0]*r03[1];
-    virial[4] += rfactor*lamda01*r01[0]*r01[2];
+    v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda03*r03[0]*r03[2];
-    virial[5] += rfactor*lamda01*r01[1]*r01[2];
+    v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda03*r03[1]*r03[2];
 
-    virial[0] += rfactor*lamda02*r02[0]*r02[0];
+    v_tally(nlist,list,4.0,v);
-    virial[1] += rfactor*lamda02*r02[1]*r02[1];
-    virial[2] += rfactor*lamda02*r02[2]*r02[2];
-    virial[3] += rfactor*lamda02*r02[0]*r02[1];
-    virial[4] += rfactor*lamda02*r02[0]*r02[2];
-    virial[5] += rfactor*lamda02*r02[1]*r02[2];
-
-    virial[0] += rfactor*lamda03*r03[0]*r03[0];
-    virial[1] += rfactor*lamda03*r03[1]*r03[1];
-    virial[2] += rfactor*lamda03*r03[2]*r03[2];
-    virial[3] += rfactor*lamda03*r03[0]*r03[1];
-    virial[4] += rfactor*lamda03*r03[0]*r03[2];
-    virial[5] += rfactor*lamda03*r03[1]*r03[2];
   }
 }
 
@@ -1643,6 +1634,9 @@ void FixShake::shake4(int m)
 
 void FixShake::shake3angle(int m)
 {
+  int nlist,list[3];
+  double v[6];
+
   // local atom IDs and constraint distances
 
   int i0 = atom->map(shake_atom[m][0]);
@@ -1845,33 +1839,20 @@ void FixShake::shake3angle(int m)
     f[i2][2] -= lamda02*r02[2] + lamda12*r12[2];
   }
 
-  if (vflag) {
+  if (evflag) {
-    int factor = 0;
+    nlist = 0;
-    if (i0 < nlocal) factor++;
+    if (i0 < nlocal) list[nlist++] = i0;
-    if (i1 < nlocal) factor++;
+    if (i1 < nlocal) list[nlist++] = i1;
-    if (i2 < nlocal) factor++;
+    if (i2 < nlocal) list[nlist++] = i2;
-    double rfactor = factor/3.0;
 
-    virial[0] += rfactor*lamda01*r01[0]*r01[0];
+    v[0] = lamda01*r01[0]*r01[0]+lamda02*r02[0]*r02[0]+lamda12*r12[0]*r12[0];
-    virial[1] += rfactor*lamda01*r01[1]*r01[1];
+    v[1] = lamda01*r01[1]*r01[1]+lamda02*r02[1]*r02[1]+lamda12*r12[1]*r12[1];
-    virial[2] += rfactor*lamda01*r01[2]*r01[2];
+    v[2] = lamda01*r01[2]*r01[2]+lamda02*r02[2]*r02[2]+lamda12*r12[2]*r12[2];
-    virial[3] += rfactor*lamda01*r01[0]*r01[1];
+    v[3] = lamda01*r01[0]*r01[1]+lamda02*r02[0]*r02[1]+lamda12*r12[0]*r12[1];
-    virial[4] += rfactor*lamda01*r01[0]*r01[2];
+    v[4] = lamda01*r01[0]*r01[2]+lamda02*r02[0]*r02[2]+lamda12*r12[0]*r12[2];
-    virial[5] += rfactor*lamda01*r01[1]*r01[2];
+    v[5] = lamda01*r01[1]*r01[2]+lamda02*r02[1]*r02[2]+lamda12*r12[1]*r12[2];
 
-    virial[0] += rfactor*lamda02*r02[0]*r02[0];
+    v_tally(nlist,list,3.0,v);
-    virial[1] += rfactor*lamda02*r02[1]*r02[1];
-    virial[2] += rfactor*lamda02*r02[2]*r02[2];
-    virial[3] += rfactor*lamda02*r02[0]*r02[1];
-    virial[4] += rfactor*lamda02*r02[0]*r02[2];
-    virial[5] += rfactor*lamda02*r02[1]*r02[2];
-
-    virial[0] += rfactor*lamda12*r12[0]*r12[0];
-    virial[1] += rfactor*lamda12*r12[1]*r12[1];
-    virial[2] += rfactor*lamda12*r12[2]*r12[2];
-    virial[3] += rfactor*lamda12*r12[0]*r12[1];
-    virial[4] += rfactor*lamda12*r12[0]*r12[2];
-    virial[5] += rfactor*lamda12*r12[1]*r12[2];
   }
 }
 
@@ -2064,6 +2045,7 @@ double FixShake::memory_usage()
   bytes += nmax*4 * sizeof(int);
   bytes += nmax*3 * sizeof(int);
   bytes += nmax*3 * sizeof(double);
+  bytes += maxvatom*6 * sizeof(double);
   return bytes;
 }
 
@@ -2199,7 +2181,7 @@ int FixShake::unpack_exchange(int nlocal, double *buf)
    rRESPA updating of atom coords is done with full v, but only portions of f
 ------------------------------------------------------------------------- */
 
-void FixShake::post_force_respa(int vflag_in, int ilevel, int iloop)
+void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
 {
   // call stats only on outermost level
 
@@ -2242,10 +2224,10 @@ void FixShake::post_force_respa(int vflag_in, int ilevel, int iloop)
 
   if (nprocs > 1) comm->comm_fix(this);
 
-  // zero out SHAKE contribution to virial
+  // virial setup
 
-  vflag = vflag_in;
+  if (vflag) v_setup(vflag);
-  if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0;
+  else evflag = 0;
 
   // loop over clusters
 

src/fix_shake.h

@@ -24,7 +24,7 @@ class FixShake : public Fix {
   ~FixShake();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void pre_neighbor();
   void post_force(int);
   void post_force_respa(int, int, int);

src/fix_spring.cpp

@@ -106,13 +106,13 @@ void FixSpring::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixSpring::setup()
+void FixSpring::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }

src/fix_spring.h

@@ -23,7 +23,7 @@ class FixSpring : public Fix {
   FixSpring(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   double compute_vector(int);

src/fix_spring_rg.cpp

@@ -74,13 +74,13 @@ void FixSpringRG::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixSpringRG::setup()
+void FixSpringRG::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }

src/fix_spring_rg.h

@@ -23,7 +23,7 @@ class FixSpringRG : public Fix {
   FixSpringRG(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
 

src/fix_spring_self.cpp

@@ -105,13 +105,13 @@ void FixSpringSelf::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixSpringSelf::setup()
+void FixSpringSelf::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }

src/fix_spring_self.h

@@ -24,7 +24,7 @@ class FixSpringSelf : public Fix {
   ~FixSpringSelf();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
 

src/fix_tmd.cpp

@@ -180,7 +180,7 @@ void FixTMD::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixTMD::initial_integrate()
+void FixTMD::initial_integrate(int vflag)
 {
   double a,b,c,d,e;
   double dx,dy,dz,dxkt,dykt,dzkt;
@@ -321,14 +321,14 @@ void FixTMD::initial_integrate()
 
 /* ---------------------------------------------------------------------- */
 
-void FixTMD::initial_integrate_respa(int ilevel, int flag)
+void FixTMD::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   if (flag) return;             // only used by NPT,NPH
 
   dtv = step_respa[ilevel];
   dtf = step_respa[ilevel] * force->ftm2v;
 
-  if (ilevel == 0) initial_integrate();
+  if (ilevel == 0) initial_integrate(vflag);
 }
 
 /* ----------------------------------------------------------------------

src/fix_tmd.h

@@ -25,8 +25,8 @@ class FixTMD : public Fix {
   ~FixTMD();
   int setmask();
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
-  void initial_integrate_respa(int,int);
+  void initial_integrate_respa(int, int, int);
 
   double memory_usage();
   void grow_arrays(int);

src/fix_viscous.cpp

@@ -76,13 +76,13 @@ void FixViscous::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixViscous::setup()
+void FixViscous::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }

src/fix_viscous.h

@@ -24,7 +24,7 @@ class FixViscous : public Fix {
   ~FixViscous();
   int setmask();
   void init();
-  void setup();
+  void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
 

src/fix_wall_lj126.cpp

@@ -109,22 +109,22 @@ void FixWallLJ126::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallLJ126::setup()
+void FixWallLJ126::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallLJ126::min_setup()
+void FixWallLJ126::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_wall_lj126.h

@@ -23,8 +23,8 @@ class FixWallLJ126 : public Fix {
   FixWallLJ126(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_wall_lj93.cpp

@@ -106,22 +106,22 @@ void FixWallLJ93::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallLJ93::setup()
+void FixWallLJ93::setup(int vflag)
 {
   if (strcmp(update->integrate_style,"verlet") == 0)
-    post_force(1);
+    post_force(vflag);
   else {
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
-    post_force_respa(1,nlevels_respa-1,0);
+    post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallLJ93::min_setup()
+void FixWallLJ93::min_setup(int vflag)
 {
-  post_force(1);
+  post_force(vflag);
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_wall_lj93.h

@@ -24,8 +24,8 @@ class FixWallLJ93 : public Fix {
   ~FixWallLJ93() {}
   int setmask();
   void init();
-  void setup();
+  void setup(int);
-  void min_setup();
+  void min_setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);

src/fix_wall_reflect.cpp

@@ -74,7 +74,7 @@ void FixWallReflect::init()
 
 /* ---------------------------------------------------------------------- */
 
-void FixWallReflect::initial_integrate()
+void FixWallReflect::initial_integrate(int vflag)
 {
   double xlo = domain->boxlo[0];
   double xhi = domain->boxhi[0];

src/fix_wall_reflect.h

@@ -23,7 +23,7 @@ class FixWallReflect : public Fix {
   FixWallReflect(class LAMMPS *, int, char **);
   int setmask();
   void init();
-  void initial_integrate();
+  void initial_integrate(int);
 
  private:
   int xloflag,xhiflag,yloflag,yhiflag,zloflag,zhiflag;

src/min_cg.cpp

@@ -259,7 +259,7 @@ void MinCG::setup()
 
   if (force->newton) comm->reverse_communicate();
 
-  modify->setup();
+  modify->setup(vflag);
   output->setup(1);
 }
 

src/modify.cpp

@@ -171,22 +171,22 @@ void Modify::init()
    setup for run, calls setup() of all fixes
 ------------------------------------------------------------------------- */
 
-void Modify::setup()
+void Modify::setup(int vflag)
 {
   if (update->whichflag == 0)
-    for (int i = 0; i < nfix; i++) fix[i]->setup();
+    for (int i = 0; i < nfix; i++) fix[i]->setup(vflag);
   else
-    for (int i = 0; i < nfix; i++) fix[i]->min_setup();
+    for (int i = 0; i < nfix; i++) fix[i]->min_setup(vflag);
 }
 
 /* ----------------------------------------------------------------------
    1st half of integrate call only for relevant fixes
 ------------------------------------------------------------------------- */
 
-void Modify::initial_integrate()
+void Modify::initial_integrate(int vflag)
 {
   for (int i = 0; i < n_initial_integrate; i++)
-    fix[list_initial_integrate[i]]->initial_integrate();
+    fix[list_initial_integrate[i]]->initial_integrate(vflag);
 }
 
 /* ----------------------------------------------------------------------
@@ -259,10 +259,11 @@ double Modify::thermo_energy()
    1st half of rRESPA integrate call only for relevant fixes
 ------------------------------------------------------------------------- */
 
-void Modify::initial_integrate_respa(int ilevel, int flag)
+void Modify::initial_integrate_respa(int vflag, int ilevel, int flag)
 {
   for (int i = 0; i < n_initial_integrate_respa; i++)
-    fix[list_initial_integrate_respa[i]]->initial_integrate_respa(ilevel,flag);
+    fix[list_initial_integrate_respa[i]]->
+      initial_integrate_respa(vflag,ilevel,flag);
 }
 
 /* ----------------------------------------------------------------------

src/modify.h

@@ -37,8 +37,8 @@ class Modify : protected Pointers {
   Modify(class LAMMPS *);
   ~Modify();
   void init();
-  void setup();
+  void setup(int);
-  void initial_integrate();
+  void initial_integrate(int);
   void pre_decide();
   void pre_exchange();
   void pre_neighbor();
@@ -47,7 +47,7 @@ class Modify : protected Pointers {
   void end_of_step();
   double thermo_energy();
 
-  void initial_integrate_respa(int,int);
+  void initial_integrate_respa(int,int,int);
   void post_force_respa(int,int,int);
   void final_integrate_respa(int);
 

src/respa.cpp

@@ -357,7 +357,7 @@ void Respa::setup()
     copy_f_flevel(ilevel);
   }
   
-  modify->setup();
+  modify->setup(vflag);
   sum_flevel_f();
   output->setup(1);
 }
@@ -373,7 +373,6 @@ void Respa::iterate(int n)
   for (int i = 0; i < n; i++) {
 
     ntimestep = ++update->ntimestep;
-
     ev_set(ntimestep);
 
     recurse(nlevels-1);
@@ -415,7 +414,7 @@ void Respa::recurse(int ilevel)
 
   for (int iloop = 0; iloop < loop[ilevel]; iloop++) {
 
-    modify->initial_integrate_respa(ilevel,0);
+    modify->initial_integrate_respa(vflag,ilevel,0);
 
     if (ilevel) recurse(ilevel-1);
 
@@ -427,7 +426,7 @@ void Respa::recurse(int ilevel)
     // at middle levels, do nothing
 
     if (ilevel == nlevels-1) {
-      modify->initial_integrate_respa(ilevel,1);
+      modify->initial_integrate_respa(vflag,ilevel,1);
       int nflag = neighbor->decide();
       if (nflag) {
 	if (modify->n_pre_exchange) modify->pre_exchange();

src/verlet.cpp

@@ -117,7 +117,7 @@ void Verlet::setup()
 
   if (force->newton) comm->reverse_communicate();
 
-  modify->setup();
+  modify->setup(vflag);
   output->setup(1);
 }
 
@@ -132,10 +132,11 @@ void Verlet::iterate(int n)
   for (int i = 0; i < n; i++) {
 
     ntimestep = ++update->ntimestep;
+    ev_set(ntimestep);
 
     // initial time integration
 
-    modify->initial_integrate();
+    modify->initial_integrate(vflag);
 
     // regular communication vs neighbor list rebuild
 
@@ -166,7 +167,6 @@ void Verlet::iterate(int n)
 
     // force computations
 
-    ev_set(ntimestep);
     force_clear();
 
     timer->stamp();