diff --git a/src/USER-REAXC/reaxc_bond_orders.cpp b/src/USER-REAXC/reaxc_bond_orders.cpp
index 9d12618373..6fdea3601e 100644
--- a/src/USER-REAXC/reaxc_bond_orders.cpp
+++ b/src/USER-REAXC/reaxc_bond_orders.cpp
@@ -557,6 +557,7 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
   /* Virial Tallying variables */
   int ii;
   real f_scaler, fi_tmp[3], fj_tmp[3], fk_tmp[3];
+  real delij[3], delki[3], delkj[3];
 
   /* Initializations */ 
   nbr_j = &(bonds->select.bond_list[pj]);
@@ -599,7 +600,10 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
     if( system->pair_ptr->vflag_atom ) {
       f_scaler = -(coef.C2dbo + coef.C2dDelta + coef.C3dbopi + coef.C3dbopi2);
       rvec_Scale( fk_tmp, -f_scaler, nbr_k->bo_data.dBOp);
-      system->pair_ptr->v_tally(k, fk_tmp);
+
+      rvec_ScaledSum( delki, 1., system->my_atoms[k].x,-1., system->my_atoms[i].x );
+      system->pair_ptr->ev_tally_xyz(k,i,system->n,1,0.0,0.0,
+	fk_tmp[0],fk_tmp[1],fk_tmp[2],delki[0],delki[1],delki[2]);
     }
   }
  
@@ -645,7 +649,10 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
     if( system->pair_ptr->vflag_atom ) {
       f_scaler = -(coef.C3dbo + coef.C3dDelta + coef.C4dbopi + coef.C4dbopi2);
       rvec_Scale( fk_tmp, -f_scaler, nbr_k->bo_data.dBOp);
-      system->pair_ptr->v_tally(k, fk_tmp);
+
+      rvec_ScaledSum( delkj, 1., system->my_atoms[k].x,-1., system->my_atoms[j].x );
+      system->pair_ptr->ev_tally_xyz(k,j,system->n,1,0.0,0.0,
+	fk_tmp[0],fk_tmp[1],fk_tmp[2],delkj[0],delkj[1],delkj[2]);
     }	
   }
   
@@ -684,6 +691,11 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
     rvec_ScaledAdd( fi_tmp, -f_scaler, bo_ij->dln_BOp_pi);
     f_scaler = coef.C1dbopi2 ;
     rvec_ScaledAdd( fi_tmp, -f_scaler, bo_ij->dln_BOp_pi2);
+    rvec_Scale( fi_tmp, 0.5, fi_tmp);
+
+    rvec_ScaledSum( delij, 1., system->my_atoms[i].x,-1., system->my_atoms[j].x );
+    system->pair_ptr->ev_tally_xyz(i,j,system->n,1,0.0,0.0,
+	fi_tmp[0],fi_tmp[1],fi_tmp[2],delij[0],delij[1],delij[2]);
 
     // forces on j
     f_scaler = -(coef.C1dbo + coef.C1dDelta + coef.C2dbopi + coef.C2dbopi2);
@@ -694,9 +706,11 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj,
     rvec_ScaledAdd( fj_tmp, -f_scaler, bo_ij->dln_BOp_pi);
     f_scaler = -coef.C1dbopi2 ;
     rvec_ScaledAdd( fj_tmp, -f_scaler, bo_ij->dln_BOp_pi2);
+    rvec_Scale( fj_tmp, 0.5, fj_tmp);
 
-    system->pair_ptr->v_tally(i, fi_tmp);
-    system->pair_ptr->v_tally(j, fj_tmp);
+    rvec_ScaledSum( delij, 1., system->my_atoms[j].x,-1., system->my_atoms[i].x );
+    system->pair_ptr->ev_tally_xyz(j,i,system->n,1,0.0,0.0,
+	fj_tmp[0],fj_tmp[1],fj_tmp[2],delij[0],delij[1],delij[2]);
 
   }
 }
diff --git a/src/compute_group_group.cpp b/src/compute_group_group.cpp
index 5393fae0be..9ee8075fad 100644
--- a/src/compute_group_group.cpp
+++ b/src/compute_group_group.cpp
@@ -47,7 +47,8 @@ ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) :
   strcpy(group2,arg[3]);
 
   jgroup = group->find(group2);
-  if (jgroup == -1) error->all(FLERR,"Compute group/group group ID does not exist");
+  if (jgroup == -1) 
+    error->all(FLERR,"Compute group/group group ID does not exist");
   jgroupbit = group->bitmask[jgroup];
 
   vector = new double[3];
@@ -80,7 +81,8 @@ void ComputeGroupGroup::init()
   // recheck that group 2 has not been deleted
 
   jgroup = group->find(group2);
-  if (jgroup == -1) error->all(FLERR,"Compute group/group group ID does not exist");
+  if (jgroup == -1) 
+    error->all(FLERR,"Compute group/group group ID does not exist");
   jgroupbit = group->bitmask[jgroup];
 
   // need an occasional half neighbor list
diff --git a/src/pair.cpp b/src/pair.cpp
index 8a8b29274a..234235211f 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -764,28 +764,6 @@ void Pair::ev_tally_list(int n, int *list, double ecoul, double *v)
   }
 }
 
-/* ----------------------------------------------------------------------
-   tally virial into per-atom accumulators
-   called by REAX/C potential, newton_pair is always on
-   fi is magnitude of force on atom i
-------------------------------------------------------------------------- */
-
-void Pair::v_tally(int i, double *fi)
-{
-  double v[6];
-  double **x = atom->x;
-  
-  v[0] = x[i][0]*fi[0];
-  v[1] = x[i][1]*fi[1];
-  v[2] = x[i][2]*fi[2];
-  v[3] = x[i][0]*fi[1];
-  v[4] = x[i][0]*fi[2];
-  v[5] = x[i][1]*fi[2];
-
-  vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
-  vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
-}
-
 /* ----------------------------------------------------------------------
    tally virial into per-atom accumulators
    called by AIREBO potential, newton_pair is always on
diff --git a/src/pair.h b/src/pair.h
index 9f531958a9..53fa3e9c43 100644
--- a/src/pair.h
+++ b/src/pair.h
@@ -87,7 +87,6 @@ class Pair : protected Pointers {
 
   // need to be public, so can be called by pair_style reaxc
 
-  void v_tally(int, double *);
   void ev_tally(int, int, int, int, double, double, double,
 		double, double, double);
   void ev_tally3(int, int, int, double, double,
@@ -95,6 +94,8 @@ class Pair : protected Pointers {
   void v_tally3(int, int, int, double *, double *, double *, double *);
   void v_tally4(int, int, int, int, double *, double *, double *,
 		double *, double *, double *);
+  void ev_tally_xyz(int, int, int, int, double, double,
+		    double, double, double, double, double, double);
 
   // general child-class methods
 
@@ -153,8 +154,6 @@ class Pair : protected Pointers {
 
   virtual void ev_setup(int, int);
   void ev_tally_full(int, double, double, double, double, double, double);
-  void ev_tally_xyz(int, int, int, int, double, double,
-		    double, double, double, double, double, double);
   void ev_tally_xyz_full(int, double, double,
 			 double, double, double, double, double, double);
   void ev_tally4(int, int, int, int, double,