diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp
index b18dbfef90..67ff81d4a3 100644
--- a/src/MANYBODY/pair_airebo.cpp
+++ b/src/MANYBODY/pair_airebo.cpp
@@ -85,10 +85,10 @@ PairAIREBO::~PairAIREBO()
 {
   memory->destroy(REBO_numneigh);
   memory->sfree(REBO_firstneigh);
-  delete [] ipage;
+  delete[] ipage;
   memory->destroy(nC);
   memory->destroy(nH);
-  delete [] pvector;
+  delete[] pvector;
 
   if (allocated) {
     memory->destroy(setflag);
@@ -100,7 +100,7 @@ PairAIREBO::~PairAIREBO()
     memory->destroy(lj2);
     memory->destroy(lj3);
     memory->destroy(lj4);
-    delete [] map;
+    delete[] map;
   }
 }
 
@@ -112,9 +112,9 @@ void PairAIREBO::compute(int eflag, int vflag)
   pvector[0] = pvector[1] = pvector[2] = 0.0;
 
   REBO_neigh();
-  FREBO(eflag,vflag);
-  if (ljflag) FLJ(eflag,vflag);
-  if (torflag) TORSION(eflag,vflag);
+  FREBO(eflag);
+  if (ljflag) FLJ(eflag);
+  if (torflag) TORSION(eflag);
 
   if (vflag_fdotr) virial_fdotr_compute();
 }
@@ -252,7 +252,7 @@ void PairAIREBO::init_style()
   if (oneatom != neighbor->oneatom) create = 1;
 
   if (create) {
-    delete [] ipage;
+    delete[] ipage;
     pgsize = neighbor->pgsize;
     oneatom = neighbor->oneatom;
 
@@ -421,7 +421,7 @@ void PairAIREBO::REBO_neigh()
    REBO forces and energy
 ------------------------------------------------------------------------- */
 
-void PairAIREBO::FREBO(int eflag, int /*vflag*/)
+void PairAIREBO::FREBO(int eflag)
 {
   int i,j,k,m,ii,inum,itype,jtype;
   tagint itag,jtag;
@@ -497,7 +497,7 @@ void PairAIREBO::FREBO(int eflag, int /*vflag*/)
       del[0] = delx;
       del[1] = dely;
       del[2] = delz;
-      bij = bondorder(i,j,del,rij,VA,f,vflag_atom);
+      bij = bondorder(i,j,del,rij,VA,f);
       dVAdi = bij*dVA;
 
       fpair = -(dVRdi+dVAdi) / rij;
@@ -520,7 +520,7 @@ void PairAIREBO::FREBO(int eflag, int /*vflag*/)
    find 3- and 4-step paths between atoms I,J via REBO neighbor lists
 ------------------------------------------------------------------------- */
 
-void PairAIREBO::FLJ(int eflag, int /*vflag*/)
+void PairAIREBO::FLJ(int eflag)
 {
   int i,j,k,m,ii,jj,kk,mm,inum,jnum,itype,jtype,ktype,mtype;
   int atomi,atomj,atomk,atomm;
@@ -786,8 +786,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
         delscale[0] = scale * delij[0];
         delscale[1] = scale * delij[1];
         delscale[2] = scale * delij[2];
-        Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA,
-                          delij,rij,f,vflag_atom);
+        Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA,delij,rij,f);
       } else Stb = 0.0;
 
       fpair = -(dStr * (Stb*cij*VLJ - cij*VLJ) +
@@ -815,7 +814,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
           f[atomj][1] -= delij[1]*fpair;
           f[atomj][2] -= delij[2]*fpair;
 
-          if (vflag_atom) v_tally2(atomi,atomj,fpair,delij);
+          if (vflag_either) v_tally2(atomi,atomj,fpair,delij);
 
         } else if (npath == 3) {
           fpair1 = dC*dwikS*wkjS / rikS;
@@ -837,7 +836,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
           f[atomk][1] -= fi[1] + fj[1];
           f[atomk][2] -= fi[2] + fj[2];
 
-          if (vflag_atom)
+          if (vflag_either)
             v_tally3(atomi,atomj,atomk,fi,fj,delikS,deljkS);
 
         } else if (npath == 4) {
@@ -873,7 +872,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
           f[atomm][1] += fm[1];
           f[atomm][2] += fm[2];
 
-          if (vflag_atom) {
+          if (vflag_either) {
             delimS[0] = delikS[0] + delkmS[0];
             delimS[1] = delikS[1] + delkmS[1];
             delimS[2] = delikS[2] + delkmS[2];
@@ -889,7 +888,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
    torsional forces and energy
 ------------------------------------------------------------------------- */
 
-void PairAIREBO::TORSION(int eflag, int /*vflag*/)
+void PairAIREBO::TORSION(int eflag)
 {
   int i,j,k,l,ii,inum;
   tagint itag,jtag;
@@ -1249,9 +1248,7 @@ void PairAIREBO::TORSION(int eflag, int /*vflag*/)
    Bij function
 ------------------------------------------------------------------------- */
 
-double PairAIREBO::bondorder(int i, int j, double rij[3],
-                             double rijmag, double VA,
-                             double **f, int vflag_atom)
+double PairAIREBO::bondorder(int i, int j, double rij[3], double rijmag, double VA, double **f)
 {
   int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
   int itype,jtype,ktype,ltype,ntype;
@@ -1435,7 +1432,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
       f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
       f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2];
 
-      if (vflag_atom) {
+      if (vflag_either) {
         rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2];
         rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2];
         v_tally3(atomi,atomj,atomk,fj,fk,rji,rki);
@@ -1577,7 +1574,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
       f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
       f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2];
 
-      if (vflag_atom) {
+      if (vflag_either) {
         rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2];
         v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj);
       }
@@ -1613,7 +1610,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
       f[atomk][1] += tmp2*rik[1];
       f[atomk][2] += tmp2*rik[2];
 
-      if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
+      if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
 
       // due to kronecker(ktype, 0) term in contribution
       // to NconjtmpI and later Nijconj
@@ -1627,7 +1624,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
       f[atomk][1] += tmp2*rik[1];
       f[atomk][2] += tmp2*rik[2];
 
-      if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
+      if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
 
       if (fabs(dNki) > TOL) {
         REBO_neighs_k = REBO_firstneigh[atomk];
@@ -1649,7 +1646,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
             f[atomn][1] += tmp2*rkn[1];
             f[atomn][2] += tmp2*rkn[2];
 
-            if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn);
+            if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn);
           }
         }
       }
@@ -1680,7 +1677,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
       f[atoml][1] += tmp2*rjl[1];
       f[atoml][2] += tmp2*rjl[2];
 
-      if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
+      if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
 
       // due to kronecker(ltype, 0) term in contribution
       // to NconjtmpJ and later Nijconj
@@ -1694,7 +1691,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
       f[atoml][1] += tmp2*rjl[1];
       f[atoml][2] += tmp2*rjl[2];
 
-      if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
+      if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
 
       if (fabs(dNlj) > TOL) {
         REBO_neighs_l = REBO_firstneigh[atoml];
@@ -1716,7 +1713,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
             f[atomn][1] += tmp2*rln[1];
             f[atomn][2] += tmp2*rln[2];
 
-            if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln);
+            if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln);
           }
         }
       }
@@ -1928,7 +1925,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
                 f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
                 f[atom4][2] += f4[2];
 
-                if (vflag_atom) {
+                if (vflag_either) {
                   r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2];
                   r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2];
                   v_tally4(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43);
@@ -1964,7 +1961,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
         f[atomk][1] += tmp2*rik[1];
         f[atomk][2] += tmp2*rik[2];
 
-        if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
+        if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
 
         // due to kronecker(ktype, 0) term in contribution
         // to NconjtmpI and later Nijconj
@@ -1978,7 +1975,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
         f[atomk][1] += tmp2*rik[1];
         f[atomk][2] += tmp2*rik[2];
 
-        if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
+        if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
 
         if (fabs(dNki) > TOL) {
           REBO_neighs_k = REBO_firstneigh[atomk];
@@ -2000,7 +1997,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
               f[atomn][1] += tmp2*rkn[1];
               f[atomn][2] += tmp2*rkn[2];
 
-              if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn);
+              if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn);
             }
           }
         }
@@ -2031,7 +2028,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
         f[atoml][1] += tmp2*rjl[1];
         f[atoml][2] += tmp2*rjl[2];
 
-        if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
+        if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
 
         // due to kronecker(ltype, 0) term in contribution
         // to NconjtmpJ and later Nijconj
@@ -2045,7 +2042,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
         f[atoml][1] += tmp2*rjl[1];
         f[atoml][2] += tmp2*rjl[2];
 
-        if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
+        if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
 
         if (fabs(dNlj) > TOL) {
           REBO_neighs_l = REBO_firstneigh[atoml];
@@ -2067,7 +2064,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
               f[atomn][1] += tmp2*rln[1];
               f[atomn][2] += tmp2*rln[2];
 
-              if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln);
+              if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln);
             }
           }
         }
@@ -2113,8 +2110,7 @@ but of the vector r_ij.
 */
 
 double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rijmag_mod,
-                               double VA, double rij[3], double rijmag,
-                               double **f, int vflag_atom)
+                               double VA, double rij[3], double rijmag, double **f)
 {
   int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
   int itype,jtype,ktype,ltype,ntype;
@@ -2438,7 +2434,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
         f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
         f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2];
 
-        if (vflag_atom) {
+        if (vflag_either) {
           rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2];
           rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2];
           v_tally3(atomi,atomj,atomk,fj,fk,rji,rki);
@@ -2542,7 +2538,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
         f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
         f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2];
 
-        if (vflag_atom) {
+        if (vflag_either) {
           rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2];
           v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj);
         }
@@ -2577,7 +2573,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
         f[atomk][1] += tmp2*rik[1];
         f[atomk][2] += tmp2*rik[2];
 
-        if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
+        if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
 
         // due to kronecker(ktype, 0) term in contribution
         // to NconjtmpI and later Nijconj
@@ -2591,7 +2587,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
         f[atomk][1] += tmp2*rik[1];
         f[atomk][2] += tmp2*rik[2];
 
-        if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
+        if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
 
         if (fabs(dNki) > TOL) {
           REBO_neighs_k = REBO_firstneigh[atomk];
@@ -2613,7 +2609,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
               f[atomn][1] += tmp2*rkn[1];
               f[atomn][2] += tmp2*rkn[2];
 
-              if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn);
+              if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn);
             }
           }
         }
@@ -2644,7 +2640,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
         f[atoml][1] += tmp2*rjl[1];
         f[atoml][2] += tmp2*rjl[2];
 
-        if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
+        if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
 
         // due to kronecker(ltype, 0) term in contribution
         // to NconjtmpJ and later Nijconj
@@ -2658,7 +2654,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
         f[atoml][1] += tmp2*rjl[1];
         f[atoml][2] += tmp2*rjl[2];
 
-        if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
+        if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
 
         if (fabs(dNlj) > TOL) {
           REBO_neighs_l = REBO_firstneigh[atoml];
@@ -2680,7 +2676,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
               f[atomn][1] += tmp2*rln[1];
               f[atomn][2] += tmp2*rln[2];
 
-              if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln);
+              if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln);
             }
           }
         }
@@ -2885,7 +2881,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
                   f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
                   f[atom4][2] += f4[2];
 
-                  if (vflag_atom) {
+                  if (vflag_either) {
                     r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2];
                     r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2];
                     v_tally4(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43);
@@ -2919,7 +2915,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
           f[atomk][1] += tmp2*rik[1];
           f[atomk][2] += tmp2*rik[2];
 
-          if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
+          if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
 
           // due to kronecker(ktype, 0) term in contribution
           // to NconjtmpI and later Nijconj
@@ -2933,7 +2929,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
           f[atomk][1] += tmp2*rik[1];
           f[atomk][2] += tmp2*rik[2];
 
-          if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
+          if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
 
           if (fabs(dNki) > TOL) {
             REBO_neighs_k = REBO_firstneigh[atomk];
@@ -2955,7 +2951,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
                 f[atomn][1] += tmp2*rkn[1];
                 f[atomn][2] += tmp2*rkn[2];
 
-                if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn);
+                if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn);
               }
             }
           }
@@ -2986,7 +2982,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
           f[atoml][1] += tmp2*rjl[1];
           f[atoml][2] += tmp2*rjl[2];
 
-          if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
+          if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
 
           // due to kronecker(ltype, 0) term in contribution
           // to NconjtmpJ and later Nijconj
@@ -3000,7 +2996,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
           f[atoml][1] += tmp2*rjl[1];
           f[atoml][2] += tmp2*rjl[2];
 
-          if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
+          if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
 
           if (fabs(dNlj) > TOL) {
             REBO_neighs_l = REBO_firstneigh[atoml];
@@ -3022,7 +3018,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
                 f[atomn][1] += tmp2*rln[1];
                 f[atomn][2] += tmp2*rln[2];
 
-                if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln);
+                if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln);
               }
             }
           }
diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h
index fa1572b1e4..fc511b7bb0 100644
--- a/src/MANYBODY/pair_airebo.h
+++ b/src/MANYBODY/pair_airebo.h
@@ -96,12 +96,12 @@ class PairAIREBO : public Pair {
   double Tf[5][5][10], Tdfdx[5][5][10], Tdfdy[5][5][10], Tdfdz[5][5][10];
 
   void REBO_neigh();
-  void FREBO(int, int);
-  void FLJ(int, int);
-  void TORSION(int, int);
+  void FREBO(int);
+  void FLJ(int);
+  void TORSION(int);
 
-  double bondorder(int, int, double *, double, double, double **, int);
-  double bondorderLJ(int, int, double *, double, double, double *, double, double **, int);
+  double bondorder(int, int, double *, double, double, double **);
+  double bondorderLJ(int, int, double *, double, double, double *, double, double **);
 
   double gSpline(double, double, int, double *, double *);
   double PijSpline(double, double, int, int, double *);
diff --git a/src/OPENMP/pair_airebo_omp.cpp b/src/OPENMP/pair_airebo_omp.cpp
index 2fff486343..6e7491787a 100644
--- a/src/OPENMP/pair_airebo_omp.cpp
+++ b/src/OPENMP/pair_airebo_omp.cpp
@@ -70,9 +70,9 @@ void PairAIREBOOMP::compute(int eflag, int vflag)
     thr->timer(Timer::START);
     ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
 
-    FREBO_thr(ifrom,ito,evflag,eflag,vflag_atom,&pv0,thr);
-    if (ljflag) FLJ_thr(ifrom,ito,evflag,eflag,vflag_atom,&pv1,thr);
-    if (torflag) TORSION_thr(ifrom,ito,evflag,eflag,&pv2,thr);
+    FREBO_thr(ifrom,ito,eflag,&pv0,thr);
+    if (ljflag) FLJ_thr(ifrom,ito,eflag,&pv1,thr);
+    if (torflag) TORSION_thr(ifrom,ito,eflag,&pv2,thr);
 
     thr->timer(Timer::PAIR);
     reduce_thr(this, eflag, vflag, thr);
@@ -184,8 +184,7 @@ void PairAIREBOOMP::REBO_neigh_thr()
    REBO forces and energy
 ------------------------------------------------------------------------- */
 
-void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
-                              int vflag_atom, double *pv0, ThrData * const thr)
+void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int eflag, double *pv0, ThrData * const thr)
 {
   int i,j,k,m,ii,itype,jtype;
   tagint itag,jtag;
@@ -259,7 +258,7 @@ void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
       del[0] = delx;
       del[1] = dely;
       del[2] = delz;
-      bij = bondorder_thr(i,j,del,rij,VA,vflag_atom,thr);
+      bij = bondorder_thr(i,j,del,rij,VA,thr);
       dVAdi = bij*dVA;
 
       fpair = -(dVRdi+dVAdi) / rij;
@@ -282,8 +281,7 @@ void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
    find 3- and 4-step paths between atoms I,J via REBO neighbor lists
 ------------------------------------------------------------------------- */
 
-void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
-                            int vflag_atom, double *pv1, ThrData * const thr)
+void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int eflag, double *pv1, ThrData * const thr)
 {
   int i,j,k,m,ii,jj,kk,mm,jnum,itype,jtype,ktype,mtype;
   tagint itag,jtag;
@@ -547,8 +545,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
         delscale[0] = scale * delij[0];
         delscale[1] = scale * delij[1];
         delscale[2] = scale * delij[2];
-        Stb = bondorderLJ_thr(i,j,delscale,rcmin[itype][jtype],VA,
-                              delij,rij,vflag_atom,thr);
+        Stb = bondorderLJ_thr(i,j,delscale,rcmin[itype][jtype],VA,delij,rij,thr);
       } else Stb = 0.0;
 
       fpair = -(dStr * (Stb*cij*VLJ - cij*VLJ) +
@@ -576,7 +573,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
           f[atomj][1] -= delij[1]*fpair;
           f[atomj][2] -= delij[2]*fpair;
 
-          if (vflag_atom) v_tally2_thr(atomi,atomj,fpair,delij,thr);
+          if (vflag_either) v_tally2_thr(this,atomi,atomj,fpair,delij,thr);
 
         } else if (npath == 3) {
           fpair1 = dC*dwikS*wkjS / rikS;
@@ -634,11 +631,11 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
           f[atomm][1] += fm[1];
           f[atomm][2] += fm[2];
 
-          if (vflag_atom) {
+          if (vflag_either) {
             delimS[0] = delikS[0] + delkmS[0];
             delimS[1] = delikS[1] + delkmS[1];
             delimS[2] = delikS[2] + delkmS[2];
-            v_tally4_thr(atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS,thr);
+            v_tally4_thr(this,atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS,thr);
           }
         }
       }
@@ -650,8 +647,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
    torsional forces and energy
 ------------------------------------------------------------------------- */
 
-void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int evflag, int eflag,
-                                double *pv2, ThrData * const thr)
+void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int eflag, double *pv2, ThrData * const thr)
 {
   int i,j,k,l,ii;
   tagint itag,jtag;
@@ -1011,7 +1007,7 @@ void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int evflag, int eflag,
 ------------------------------------------------------------------------- */
 
 double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
-                                    double VA, int vflag_atom, ThrData * const thr)
+                                    double VA, ThrData * const thr)
 {
   int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
   int itype,jtype,ktype,ltype,ntype;
@@ -1378,7 +1374,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
       f[atomk][1] += tmp2*rik[1];
       f[atomk][2] += tmp2*rik[2];
 
-      if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
+      if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
 
       // due to kronecker(ktype, 0) term in contribution
       // to NconjtmpI and later Nijconj
@@ -1392,7 +1388,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
       f[atomk][1] += tmp2*rik[1];
       f[atomk][2] += tmp2*rik[2];
 
-      if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
+      if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
 
       if (fabs(dNki) > TOL) {
         REBO_neighs_k = REBO_firstneigh[atomk];
@@ -1414,7 +1410,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
             f[atomn][1] += tmp2*rkn[1];
             f[atomn][2] += tmp2*rkn[2];
 
-            if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
+            if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr);
           }
         }
       }
@@ -1445,7 +1441,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
       f[atoml][1] += tmp2*rjl[1];
       f[atoml][2] += tmp2*rjl[2];
 
-      if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
+      if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
 
       // due to kronecker(ltype, 0) term in contribution
       // to NconjtmpJ and later Nijconj
@@ -1459,7 +1455,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
       f[atoml][1] += tmp2*rjl[1];
       f[atoml][2] += tmp2*rjl[2];
 
-      if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
+      if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
 
       if (fabs(dNlj) > TOL) {
         REBO_neighs_l = REBO_firstneigh[atoml];
@@ -1481,7 +1477,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
             f[atomn][1] += tmp2*rln[1];
             f[atomn][2] += tmp2*rln[2];
 
-            if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
+            if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr);
           }
         }
       }
@@ -1693,10 +1689,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
                 f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
                 f[atom4][2] += f4[2];
 
-                if (vflag_atom) {
+                if (vflag_either) {
                   r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2];
                   r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2];
-                  v_tally4_thr(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
+                  v_tally4_thr(this,atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
                 }
               }
             }
@@ -1729,7 +1725,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
         f[atomk][1] += tmp2*rik[1];
         f[atomk][2] += tmp2*rik[2];
 
-        if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
+        if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
 
         // due to kronecker(ktype, 0) term in contribution
         // to NconjtmpI and later Nijconj
@@ -1743,7 +1739,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
         f[atomk][1] += tmp2*rik[1];
         f[atomk][2] += tmp2*rik[2];
 
-        if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
+        if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
 
         if (fabs(dNki) > TOL) {
           REBO_neighs_k = REBO_firstneigh[atomk];
@@ -1765,7 +1761,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
               f[atomn][1] += tmp2*rkn[1];
               f[atomn][2] += tmp2*rkn[2];
 
-              if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
+              if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr);
             }
           }
         }
@@ -1796,7 +1792,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
         f[atoml][1] += tmp2*rjl[1];
         f[atoml][2] += tmp2*rjl[2];
 
-        if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
+        if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
 
         // due to kronecker(ltype, 0) term in contribution
         // to NconjtmpJ and later Nijconj
@@ -1810,7 +1806,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
         f[atoml][1] += tmp2*rjl[1];
         f[atoml][2] += tmp2*rjl[2];
 
-        if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
+        if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
 
         if (fabs(dNlj) > TOL) {
           REBO_neighs_l = REBO_firstneigh[atoml];
@@ -1832,7 +1828,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
               f[atomn][1] += tmp2*rln[1];
               f[atomn][2] += tmp2*rln[2];
 
-              if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
+              if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr);
             }
           }
         }
@@ -1870,8 +1866,7 @@ there probably also need to be performed here.
 */
 
 double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], double rijmag_mod,
-                                      double VA, double rij[3], double rijmag,
-                                      int vflag_atom, ThrData * const thr)
+                                      double VA, double rij[3], double rijmag, ThrData * const thr)
 {
   int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
   int itype,jtype,ktype,ltype,ntype;
@@ -2335,7 +2330,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
         f[atomk][1] += tmp2*rik[1];
         f[atomk][2] += tmp2*rik[2];
 
-        if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
+        if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
 
         // due to kronecker(ktype, 0) term in contribution
         // to NconjtmpI and later Nijconj
@@ -2349,7 +2344,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
         f[atomk][1] += tmp2*rik[1];
         f[atomk][2] += tmp2*rik[2];
 
-        if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
+        if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
 
         if (fabs(dNki) > TOL) {
           REBO_neighs_k = REBO_firstneigh[atomk];
@@ -2371,7 +2366,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
               f[atomn][1] += tmp2*rkn[1];
               f[atomn][2] += tmp2*rkn[2];
 
-              if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
+              if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr);
             }
           }
         }
@@ -2402,7 +2397,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
         f[atoml][1] += tmp2*rjl[1];
         f[atoml][2] += tmp2*rjl[2];
 
-        if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
+        if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
 
         // due to kronecker(ltype, 0) term in contribution
         // to NconjtmpJ and later Nijconj
@@ -2416,7 +2411,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
         f[atoml][1] += tmp2*rjl[1];
         f[atoml][2] += tmp2*rjl[2];
 
-        if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
+        if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
 
         if (fabs(dNlj) > TOL) {
           REBO_neighs_l = REBO_firstneigh[atoml];
@@ -2438,7 +2433,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
               f[atomn][1] += tmp2*rln[1];
               f[atomn][2] += tmp2*rln[2];
 
-              if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
+              if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr);
             }
           }
         }
@@ -2643,10 +2638,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
                   f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
                   f[atom4][2] += f4[2];
 
-                  if (vflag_atom) {
+                  if (vflag_either) {
                     r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2];
                     r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2];
-                    v_tally4_thr(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
+                    v_tally4_thr(this,atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
                   }
                 }
               }
@@ -2677,7 +2672,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
           f[atomk][1] += tmp2*rik[1];
           f[atomk][2] += tmp2*rik[2];
 
-          if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
+          if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
 
           // due to kronecker(ktype, 0) term in contribution
           // to NconjtmpI and later Nijconj
@@ -2691,7 +2686,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
           f[atomk][1] += tmp2*rik[1];
           f[atomk][2] += tmp2*rik[2];
 
-          if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
+          if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
 
           if (fabs(dNki) > TOL) {
             REBO_neighs_k = REBO_firstneigh[atomk];
@@ -2713,7 +2708,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
                 f[atomn][1] += tmp2*rkn[1];
                 f[atomn][2] += tmp2*rkn[2];
 
-                if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
+                if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr);
               }
             }
           }
@@ -2744,7 +2739,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
           f[atoml][1] += tmp2*rjl[1];
           f[atoml][2] += tmp2*rjl[2];
 
-          if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
+          if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
 
           // due to kronecker(ltype, 0) term in contribution
           // to NconjtmpJ and later Nijconj
@@ -2758,7 +2753,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
           f[atoml][1] += tmp2*rjl[1];
           f[atoml][2] += tmp2*rjl[2];
 
-          if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
+          if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
 
           if (fabs(dNlj) > TOL) {
             REBO_neighs_l = REBO_firstneigh[atoml];
@@ -2780,7 +2775,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
                 f[atomn][1] += tmp2*rln[1];
                 f[atomn][2] += tmp2*rln[2];
 
-                if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
+                if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr);
               }
             }
           }
diff --git a/src/OPENMP/pair_airebo_omp.h b/src/OPENMP/pair_airebo_omp.h
index 9c5b29545f..355614552e 100644
--- a/src/OPENMP/pair_airebo_omp.h
+++ b/src/OPENMP/pair_airebo_omp.h
@@ -33,16 +33,13 @@ class PairAIREBOOMP : public PairAIREBO, public ThrOMP {
   virtual double memory_usage();
 
  protected:
-  double bondorder_thr(int i, int j, double rij[3], double rijmag, double VA, int vflag_atom,
-                       ThrData *const thr);
+  double bondorder_thr(int i, int j, double rij[3], double rijmag, double VA, ThrData *const thr);
   double bondorderLJ_thr(int i, int j, double rij[3], double rijmag, double VA, double rij0[3],
-                         double rijmag0, int vflag_atom, ThrData *const thr);
+                         double rijmag0, ThrData *const thr);
 
-  void FREBO_thr(int ifrom, int ito, int evflag, int eflag, int vflag_atom, double *pv0,
-                 ThrData *const thr);
-  void FLJ_thr(int ifrom, int ito, int evflag, int eflag, int vflag_atom, double *pv1,
-               ThrData *const thr);
-  void TORSION_thr(int ifrom, int ito, int evflag, int eflag, double *pv2, ThrData *const thr);
+  void FREBO_thr(int ifrom, int ito, int eflag, double *pv0, ThrData *const thr);
+  void FLJ_thr(int ifrom, int ito, int eflag, double *pv1, ThrData *const thr);
+  void TORSION_thr(int ifrom, int ito, int eflag, double *pv2, ThrData *const thr);
   void REBO_neigh_thr();
 };
 
diff --git a/src/OPENMP/thr_omp.cpp b/src/OPENMP/thr_omp.cpp
index 867623361c..319a1f45df 100644
--- a/src/OPENMP/thr_omp.cpp
+++ b/src/OPENMP/thr_omp.cpp
@@ -1465,20 +1465,31 @@ void ThrOMP::ev_tally_thr(Improper * const imprp, const int i1, const int i2,
    fpair is magnitude of force on atom I
 ------------------------------------------------------------------------- */
 
-void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair,
+void ThrOMP::v_tally2_thr(Pair *const pair, const int i, const int j, const double fpair,
                           const double * const drij, ThrData * const thr)
 {
   double v[6];
 
-  v[0] = 0.5 * drij[0]*drij[0]*fpair;
-  v[1] = 0.5 * drij[1]*drij[1]*fpair;
-  v[2] = 0.5 * drij[2]*drij[2]*fpair;
-  v[3] = 0.5 * drij[0]*drij[1]*fpair;
-  v[4] = 0.5 * drij[0]*drij[2]*fpair;
-  v[5] = 0.5 * drij[1]*drij[2]*fpair;
+  v[0] = drij[0]*drij[0]*fpair;
+  v[1] = drij[1]*drij[1]*fpair;
+  v[2] = drij[2]*drij[2]*fpair;
+  v[3] = drij[0]*drij[1]*fpair;
+  v[4] = drij[0]*drij[2]*fpair;
+  v[5] = drij[1]*drij[2]*fpair;
 
-  v_tally(thr->vatom_pair[i],v);
-  v_tally(thr->vatom_pair[j],v);
+  if (pair->vflag_global) v_tally(thr->virial_pair,v);
+
+  if (pair->vflag_atom) {
+    v[0] *= 0.5;
+    v[1] *= 0.5;
+    v[2] *= 0.5;
+    v[3] *= 0.5;
+    v[4] *= 0.5;
+    v[5] *= 0.5;
+
+    v_tally(thr->vatom_pair[i],v);
+    v_tally(thr->vatom_pair[j],v);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1486,7 +1497,7 @@ void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair,
    called by AIREBO and Tersoff potential, newton_pair is always on
 ------------------------------------------------------------------------- */
 
-void ThrOMP::v_tally3_thr(Pair *pair, const int i, const int j, const int k,
+void ThrOMP::v_tally3_thr(Pair *const pair, const int i, const int j, const int k,
                           const double * const fi, const double * const fj,
                           const double * const drik, const double * const drjk,
                           ThrData * const thr)
@@ -1521,7 +1532,7 @@ void ThrOMP::v_tally3_thr(Pair *pair, const int i, const int j, const int k,
    called by AIREBO potential, newton_pair is always on
 ------------------------------------------------------------------------- */
 
-void ThrOMP::v_tally4_thr(const int i, const int j, const int k, const int m,
+void ThrOMP::v_tally4_thr(Pair *const pair, const int i, const int j, const int k, const int m,
                           const double * const fi, const double * const fj,
                           const double * const fk, const double * const drim,
                           const double * const drjm, const double * const drkm,
@@ -1529,17 +1540,26 @@ void ThrOMP::v_tally4_thr(const int i, const int j, const int k, const int m,
 {
   double v[6];
 
-  v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
-  v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
-  v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
-  v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
-  v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
-  v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
+  v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
+  v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
+  v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
+  v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
+  v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
+  v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
+  if (pair->vflag_global) v_tally(thr->virial_pair,v);
 
-  v_tally(thr->vatom_pair[i],v);
-  v_tally(thr->vatom_pair[j],v);
-  v_tally(thr->vatom_pair[k],v);
-  v_tally(thr->vatom_pair[m],v);
+  if (pair->vflag_atom) {
+    v[0] *= 0.25;
+    v[1] *= 0.25;
+    v[2] *= 0.25;
+    v[3] *= 0.25;
+    v[4] *= 0.25;
+    v[5] *= 0.25;
+    v_tally(thr->vatom_pair[i],v);
+    v_tally(thr->vatom_pair[j],v);
+    v_tally(thr->vatom_pair[k],v);
+    v_tally(thr->vatom_pair[m],v);
+  }
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/OPENMP/thr_omp.h b/src/OPENMP/thr_omp.h
index 22fdb36f38..35cd43be5d 100644
--- a/src/OPENMP/thr_omp.h
+++ b/src/OPENMP/thr_omp.h
@@ -135,6 +135,8 @@ class ThrOMP {
   void ev_tally_xyz_full_thr(Pair *const, const int, const double, const double, const double,
                              const double, const double, const double, const double, const double,
                              ThrData *const);
+  void v_tally2_thr(Pair *const, const int, const int, const double, const double *const,
+                    ThrData *const);
   void ev_tally3_thr(Pair *const, const int, const int, const int, const double, const double,
                      const double *const, const double *const, const double *const,
                      const double *const, ThrData *const);
@@ -143,6 +145,9 @@ class ThrOMP {
   void ev_tally4_thr(Pair *const, const int, const int, const int, const int, const double,
                      const double *const, const double *const, const double *const,
                      const double *const, const double *const, const double *const, ThrData *const);
+  void v_tally4_thr(Pair *const, const int, const int, const int, const int, const double *const,
+                    const double *const, const double *const, const double *const,
+                    const double *const, const double *const, ThrData *const);
 
   // Bond
   void ev_tally_thr(Bond *const, const int, const int, const int, const int, const double,
@@ -171,10 +176,6 @@ class ThrOMP {
                     ThrData *const);
 
   // style independent versions
-  void v_tally2_thr(const int, const int, const double, const double *const, ThrData *const);
-  void v_tally4_thr(const int, const int, const int, const int, const double *const,
-                    const double *const, const double *const, const double *const,
-                    const double *const, const double *const, ThrData *const);
   void ev_tally_list_thr(Pair *const, const int, const int *const, const double *const,
                          const double, const double, ThrData *const);
 };
diff --git a/src/pair.cpp b/src/pair.cpp
index f27d2ac370..1a6a476b4f 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -1499,17 +1499,35 @@ void Pair::v_tally2(int i, int j, double fpair, double *drij)
 {
   double v[6];
 
-  v[0] = 0.5 * drij[0]*drij[0]*fpair;
-  v[1] = 0.5 * drij[1]*drij[1]*fpair;
-  v[2] = 0.5 * drij[2]*drij[2]*fpair;
-  v[3] = 0.5 * drij[0]*drij[1]*fpair;
-  v[4] = 0.5 * drij[0]*drij[2]*fpair;
-  v[5] = 0.5 * drij[1]*drij[2]*fpair;
+  v[0] = drij[0]*drij[0]*fpair;
+  v[1] = drij[1]*drij[1]*fpair;
+  v[2] = drij[2]*drij[2]*fpair;
+  v[3] = drij[0]*drij[1]*fpair;
+  v[4] = drij[0]*drij[2]*fpair;
+  v[5] = drij[1]*drij[2]*fpair;
 
-  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];
-  vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
-  vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
+  if (vflag_global) {
+    virial[0] += v[0];
+    virial[1] += v[1];
+    virial[2] += v[2];
+    virial[3] += v[3];
+    virial[4] += v[4];
+    virial[5] += v[5];
+  }
+  
+  if (vflag_atom) {
+    v[0] *= 0.5;
+    v[1] *= 0.5;
+    v[2] *= 0.5;
+    v[3] *= 0.5;
+    v[4] *= 0.5;
+    v[5] *= 0.5;
+
+    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];
+    vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
+    vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -1522,23 +1540,29 @@ void Pair::v_tally3(int i, int j, int k,
 {
   double v[6];
 
-  v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]);
-  v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]);
-  v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]);
-  v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]);
-  v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]);
-  v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]);
+  v[0] = (drik[0]*fi[0] + drjk[0]*fj[0]);
+  v[1] = (drik[1]*fi[1] + drjk[1]*fj[1]);
+  v[2] = (drik[2]*fi[2] + drjk[2]*fj[2]);
+  v[3] = (drik[0]*fi[1] + drjk[0]*fj[1]);
+  v[4] = (drik[0]*fi[2] + drjk[0]*fj[2]);
+  v[5] = (drik[1]*fi[2] + drjk[1]*fj[2]);
 
   if (vflag_global) {
-      virial[0] += 3.0*v[0];
-      virial[1] += 3.0*v[1];
-      virial[2] += 3.0*v[2];
-      virial[3] += 3.0*v[3];
-      virial[4] += 3.0*v[4];
-      virial[5] += 3.0*v[5];
+      virial[0] += v[0];
+      virial[1] += v[1];
+      virial[2] += v[2];
+      virial[3] += v[3];
+      virial[4] += v[4];
+      virial[5] += v[5];
   }
 
   if (vflag_atom) {
+    v[0] *= THIRD;
+    v[1] *= THIRD;
+    v[2] *= THIRD;
+    v[3] *= THIRD;
+    v[4] *= THIRD;
+    v[5] *= THIRD;
     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];
     vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
@@ -1559,21 +1583,39 @@ void Pair::v_tally4(int i, int j, int k, int m,
 {
   double v[6];
 
-  v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
-  v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
-  v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
-  v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
-  v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
-  v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
+  v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
+  v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
+  v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
+  v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
+  v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
+  v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[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];
-  vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
-  vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
-  vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2];
-  vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5];
-  vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2];
-  vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5];
+  if (vflag_global) {
+      virial[0] += v[0];
+      virial[1] += v[1];
+      virial[2] += v[2];
+      virial[3] += v[3];
+      virial[4] += v[4];
+      virial[5] += v[5];
+  }
+
+  if (vflag_atom) {
+    v[0] *= 0.25;
+    v[1] *= 0.25;
+    v[2] *= 0.25;
+    v[3] *= 0.25;
+    v[4] *= 0.25;
+    v[5] *= 0.25;
+
+    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];
+    vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
+    vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
+    vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2];
+    vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5];
+    vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2];
+    vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5];
+  }
 }
 
 /* ----------------------------------------------------------------------
diff --git a/unittest/force-styles/test_pair_style.cpp b/unittest/force-styles/test_pair_style.cpp
index 9f70a22daf..d256685659 100644
--- a/unittest/force-styles/test_pair_style.cpp
+++ b/unittest/force-styles/test_pair_style.cpp
@@ -1089,12 +1089,14 @@ TEST(PairStyle, intel)
     EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
     if (print_stats) std::cerr << "run_energy  stats:" << stats << std::endl;
 
-    // pair styles sw and tersoff require newton on, but that also requires fdotr for /intel
-    std::cerr << "pair style : " << test_config.pair_style << "\n";
-    if ((test_config.pair_style != "sw") && (test_config.pair_style != "tersoff")) {
+    // pair styles sw and tersoff and airebo INTEL package variants require newton on,
+    // but that also requires fdotr for /intel
+    if ((test_config.pair_style != "sw") && (test_config.pair_style != "tersoff")
+        && (test_config.pair_style != "rebo") && (test_config.pair_style != "airebo")
+        && (test_config.pair_style != "airebo/morse")) {
 
         if (!verbose) ::testing::internal::CaptureStdout();
-        restart_lammps(lmp, test_config, true, false);
+        restart_lammps(lmp, test_config, true, true);
         if (!verbose) ::testing::internal::GetCapturedStdout();
         f   = lmp->atom->f;
         tag = lmp->atom->tag;