diff --git a/src/AMOEBA/amoeba_charge_transfer.cpp b/src/AMOEBA/amoeba_charge_transfer.cpp
index 3ab5171d93..a0928e1e63 100644
--- a/src/AMOEBA/amoeba_charge_transfer.cpp
+++ b/src/AMOEBA/amoeba_charge_transfer.cpp
@@ -148,25 +148,20 @@ void PairAmoeba::charge_transfer()
 
       // increment the internal virial tensor components
 
-      vxx = xr * frcx;
-      vxy = yr * frcx;
-      vxz = zr * frcx;
-      vyy = yr * frcy;
-      vyz = zr * frcy;
-      vzz = zr * frcz;
+      if (vflag_global) {
+        vxx = xr * frcx;
+        vxy = yr * frcx;
+        vxz = zr * frcx;
+        vyy = yr * frcy;
+        vyz = zr * frcy;
+        vzz = zr * frcz;
 
-      virqxfer[0] -= vxx;
-      virqxfer[1] -= vyy;
-      virqxfer[2] -= vzz;
-      virqxfer[3] -= vxy;
-      virqxfer[4] -= vxz;
-      virqxfer[5] -= vyz;
-
-      // energy = e
-      // virial = 6-vec vir
-      // NOTE: add tally function
-
-      if (evflag) {
+        virqxfer[0] -= vxx;
+        virqxfer[1] -= vyy;
+        virqxfer[2] -= vzz;
+        virqxfer[3] -= vxy;
+        virqxfer[4] -= vxz;
+        virqxfer[5] -= vyz;
       }
     }
   }
diff --git a/src/AMOEBA/amoeba_dispersion.cpp b/src/AMOEBA/amoeba_dispersion.cpp
index 5dd7aacc88..83405cfab9 100644
--- a/src/AMOEBA/amoeba_dispersion.cpp
+++ b/src/AMOEBA/amoeba_dispersion.cpp
@@ -212,25 +212,20 @@ void PairAmoeba::dispersion_real()
 
       // increment the internal virial tensor components
 
-      vxx = xr * dedx;
-      vyx = yr * dedx;
-      vzx = zr * dedx;
-      vyy = yr * dedy;
-      vzy = zr * dedy;
-      vzz = zr * dedz;
+      if (vflag_global) {
+        vxx = xr * dedx;
+        vyx = yr * dedx;
+        vzx = zr * dedx;
+        vyy = yr * dedy;
+        vzy = zr * dedy;
+        vzz = zr * dedz;
 
-      virdisp[0] -= vxx;
-      virdisp[1] -= vyy;
-      virdisp[2] -= vzz;
-      virdisp[3] -= vyx;
-      virdisp[4] -= vzx;
-      virdisp[5] -= vzy;
-
-      // energy = e
-      // virial = 6-vec vir
-      // NOTE: add tally function
-
-      if (evflag) {
+        virdisp[0] -= vxx;
+        virdisp[1] -= vyy;
+        virdisp[2] -= vzz;
+        virdisp[3] -= vyx;
+        virdisp[4] -= vzx;
+        virdisp[5] -= vzy;
       }
     }
   }
@@ -350,13 +345,15 @@ void PairAmoeba::dispersion_kspace()
           struc2 = gridfft[n]*gridfft[n] + gridfft[n+1]*gridfft[n+1];
           e = -(term1 / denom) * struc2;
           edisp += e;
-          vterm = 3.0 * (fac1*erfcterm*h + fac3*expterm) * struc2/denom;
-          virdisp[0] -= h1*h1*vterm - e;
-          virdisp[1] -= h2*h2*vterm - e;
-          virdisp[2] -= h3*h3*vterm - e;
-          virdisp[3] -= h1*h2*vterm;
-          virdisp[4] -= h1*h3*vterm;
-          virdisp[5] -= h2*h3*vterm;
+          if (vflag_global) {
+            vterm = 3.0 * (fac1*erfcterm*h + fac3*expterm) * struc2/denom;
+            virdisp[0] -= h1*h1*vterm - e;
+            virdisp[1] -= h2*h2*vterm - e;
+            virdisp[2] -= h3*h3*vterm - e;
+            virdisp[3] -= h1*h2*vterm;
+            virdisp[4] -= h1*h3*vterm;
+            virdisp[5] -= h2*h3*vterm;
+          }
         } else term1 = 0.0;
         // NOTE: pre-calc this division only once
         gridfft[n] *= -(term1/denom);
@@ -418,8 +415,10 @@ void PairAmoeba::dispersion_kspace()
 
   if (me == 0) {
     edisp -= term;
-    virdisp[0] -= term;
-    virdisp[1] -= term;
-    virdisp[2] -= term;
+    if (vflag_global) {
+      virdisp[0] -= term;
+      virdisp[1] -= term;
+      virdisp[2] -= term;
+    }
   }
 }
diff --git a/src/AMOEBA/amoeba_hal.cpp b/src/AMOEBA/amoeba_hal.cpp
index 561a2869a2..d73a24ab21 100644
--- a/src/AMOEBA/amoeba_hal.cpp
+++ b/src/AMOEBA/amoeba_hal.cpp
@@ -188,25 +188,20 @@ void PairAmoeba::hal()
 
       // increment the internal virial tensor components
 
-      vxx = xr * dedx;
-      vyx = yr * dedx;
-      vzx = zr * dedx;
-      vyy = yr * dedy;
-      vzy = zr * dedy;
-      vzz = zr * dedz;
+      if (vflag_global) {
+        vxx = xr * dedx;
+        vyx = yr * dedx;
+        vzx = zr * dedx;
+        vyy = yr * dedy;
+        vzy = zr * dedy;
+        vzz = zr * dedz;
 
-      virhal[0] -= vxx;
-      virhal[1] -= vyy;
-      virhal[2] -= vzz;
-      virhal[3] -= vyx;
-      virhal[4] -= vzx;
-      virhal[5] -= vzy;
-
-      // energy = e
-      // virial = 6-vec vir
-      // NOTE: add tally function
-
-      if (evflag) {
+        virhal[0] -= vxx;
+        virhal[1] -= vyy;
+        virhal[2] -= vzz;
+        virhal[3] -= vyx;
+        virhal[4] -= vzx;
+        virhal[5] -= vzy;
       }
     }
   }
diff --git a/src/AMOEBA/amoeba_multipole.cpp b/src/AMOEBA/amoeba_multipole.cpp
index 2c666419d1..faf54584c6 100644
--- a/src/AMOEBA/amoeba_multipole.cpp
+++ b/src/AMOEBA/amoeba_multipole.cpp
@@ -497,25 +497,20 @@ void PairAmoeba::multipole_real()
 
       // increment the virial due to pairwise Cartesian forces
 
-      vxx = -xr * frcx;
-      vxy = -0.5 * (yr*frcx+xr*frcy);
-      vxz = -0.5 * (zr*frcx+xr*frcz);
-      vyy = -yr * frcy;
-      vyz = -0.5 * (zr*frcy+yr*frcz);
-      vzz = -zr * frcz;
-
-      virmpole[0] -= vxx;
-      virmpole[1] -= vyy;
-      virmpole[2] -= vzz;
-      virmpole[3] -= vxy;
-      virmpole[4] -= vxz;
-      virmpole[5] -= vyz;
-
-      // energy = e
-      // virial = 6-vec vir
-      // NOTE: add tally function
-
-      if (evflag) {
+      if (vflag_global) {
+        vxx = -xr * frcx;
+        vxy = -0.5 * (yr*frcx+xr*frcy);
+        vxz = -0.5 * (zr*frcx+xr*frcz);
+        vyy = -yr * frcy;
+        vyz = -0.5 * (zr*frcy+yr*frcz);
+        vzz = -zr * frcz;
+      
+        virmpole[0] -= vxx;
+        virmpole[1] -= vyy;
+        virmpole[2] -= vzz;
+        virmpole[3] -= vxy;
+        virmpole[4] -= vxz;
+        virmpole[5] -= vyz;
       }
     }
   }
@@ -530,6 +525,8 @@ void PairAmoeba::multipole_real()
   for (i = 0; i < nlocal; i++) {
     torque2force(i,tq[i],fix,fiy,fiz,f);
 
+    if (!vflag_global) continue;
+
     iz = zaxis2local[i];
     ix = xaxis2local[i];
     iy = yaxis2local[i];
@@ -543,7 +540,7 @@ void PairAmoeba::multipole_real()
     xiy = x[iy][0] - x[i][0];
     yiy = x[iy][1] - x[i][1];
     ziy = x[iy][2] - x[i][2];
-
+    
     vxx = xix*fix[0] + xiy*fiy[0] + xiz*fiz[0];
     vxy = 0.5 * (yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
                  xix*fix[1] + xiy*fiy[1] + xiz*fiz[1]);
@@ -553,7 +550,7 @@ void PairAmoeba::multipole_real()
     vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
                  yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
     vzz = zix*fix[2] + ziy*fiy[2] + ziz*fiz[2];
-
+    
     virmpole[0] -= vxx;
     virmpole[1] -= vyy;
     virmpole[2] -= vzz;
@@ -768,22 +765,24 @@ void PairAmoeba::multipole_kspace()
 
   // augment the permanent multipole virial contributions
 
-  for (i = 0; i < nlocal; i++) {
-    vxx = vxx - cmp[i][1]*cphi[i][1] - 2.0*cmp[i][4]*cphi[i][4] -
-      cmp[i][7]*cphi[i][7] - cmp[i][8]*cphi[i][8];
-    vxy = vxy - 0.5*(cmp[i][2]*cphi[i][1]+cmp[i][1]*cphi[i][2]) -
-      (cmp[i][4]+cmp[i][5])*cphi[i][7] - 0.5*cmp[i][7]*(cphi[i][4]+cphi[i][5]) -
-      0.5*(cmp[i][8]*cphi[i][9]+cmp[i][9]*cphi[i][8]);
-    vxz = vxz - 0.5*(cmp[i][3]*cphi[i][1]+cmp[i][1]*cphi[i][3]) -
-      (cmp[i][4]+cmp[i][6])*cphi[i][8] - 0.5*cmp[i][8]*(cphi[i][4]+cphi[i][6]) -
-      0.5*(cmp[i][7]*cphi[i][9]+cmp[i][9]*cphi[i][7]);
-    vyy = vyy - cmp[i][2]*cphi[i][2] - 2.0*cmp[i][5]*cphi[i][5] -
-      cmp[i][7]*cphi[i][7] - cmp[i][9]*cphi[i][9];
-    vyz = vyz - 0.5*(cmp[i][3]*cphi[i][2]+cmp[i][2]*cphi[i][3]) -
-      (cmp[i][5]+cmp[i][6])*cphi[i][9] - 0.5*cmp[i][9]*(cphi[i][5]+cphi[i][6]) -
-      0.5*(cmp[i][7]*cphi[i][8]+cmp[i][8]*cphi[i][7]);
-    vzz = vzz - cmp[i][3]*cphi[i][3] - 2.0*cmp[i][6]*cphi[i][6] -
-      cmp[i][8]*cphi[i][8] - cmp[i][9]*cphi[i][9];
+  if (vflag_global) {
+    for (i = 0; i < nlocal; i++) {
+      vxx = vxx - cmp[i][1]*cphi[i][1] - 2.0*cmp[i][4]*cphi[i][4] -
+        cmp[i][7]*cphi[i][7] - cmp[i][8]*cphi[i][8];
+      vxy = vxy - 0.5*(cmp[i][2]*cphi[i][1]+cmp[i][1]*cphi[i][2]) -
+        (cmp[i][4]+cmp[i][5])*cphi[i][7] - 0.5*cmp[i][7]*(cphi[i][4]+cphi[i][5]) -
+        0.5*(cmp[i][8]*cphi[i][9]+cmp[i][9]*cphi[i][8]);
+      vxz = vxz - 0.5*(cmp[i][3]*cphi[i][1]+cmp[i][1]*cphi[i][3]) -
+        (cmp[i][4]+cmp[i][6])*cphi[i][8] - 0.5*cmp[i][8]*(cphi[i][4]+cphi[i][6]) -
+        0.5*(cmp[i][7]*cphi[i][9]+cmp[i][9]*cphi[i][7]);
+      vyy = vyy - cmp[i][2]*cphi[i][2] - 2.0*cmp[i][5]*cphi[i][5] -
+        cmp[i][7]*cphi[i][7] - cmp[i][9]*cphi[i][9];
+      vyz = vyz - 0.5*(cmp[i][3]*cphi[i][2]+cmp[i][2]*cphi[i][3]) -
+        (cmp[i][5]+cmp[i][6])*cphi[i][9] - 0.5*cmp[i][9]*(cphi[i][5]+cphi[i][6]) -
+        0.5*(cmp[i][7]*cphi[i][8]+cmp[i][8]*cphi[i][7]);
+      vzz = vzz - cmp[i][3]*cphi[i][3] - 2.0*cmp[i][6]*cphi[i][6] -
+        cmp[i][8]*cphi[i][8] - cmp[i][9]*cphi[i][9];
+    }
   }
 
   // resolve site torques then increment forces and virial
@@ -804,39 +803,43 @@ void PairAmoeba::multipole_kspace()
 
     torque2force(i,tem,fix,fiy,fiz,f);
 
-    iz = zaxis2local[i];
-    ix = xaxis2local[i];
-    iy = yaxis2local[i];
+    if (vflag_global) {
+      iz = zaxis2local[i];
+      ix = xaxis2local[i];
+      iy = yaxis2local[i];
 
-    xiz = x[iz][0] - x[i][0];
-    yiz = x[iz][1] - x[i][1];
-    ziz = x[iz][2] - x[i][2];
-    xix = x[ix][0] - x[i][0];
-    yix = x[ix][1] - x[i][1];
-    zix = x[ix][2] - x[i][2];
-    xiy = x[iy][0] - x[i][0];
-    yiy = x[iy][1] - x[i][1];
-    ziy = x[iy][2] - x[i][2];
+      xiz = x[iz][0] - x[i][0];
+      yiz = x[iz][1] - x[i][1];
+      ziz = x[iz][2] - x[i][2];
+      xix = x[ix][0] - x[i][0];
+      yix = x[ix][1] - x[i][1];
+      zix = x[ix][2] - x[i][2];
+      xiy = x[iy][0] - x[i][0];
+      yiy = x[iy][1] - x[i][1];
+      ziy = x[iy][2] - x[i][2];
 
-    vxx += xix*fix[0] + xiy*fiy[0] + xiz*fiz[0];
-    vxy += 0.5*(yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
-                xix*fix[1] + xiy*fiy[1] + xiz*fiz[1]);
-    vxz += 0.5*(zix*fix[0] + ziy*fiy[0] + ziz*fiz[0] +
-                xix*fix[2] + xiy*fiy[2] + xiz*fiz[2]);
-    vyy += yix*fix[1] + yiy*fiy[1] + yiz*fiz[1];
-    vyz += 0.5*(zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
-                yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
-    vzz += zix*fix[2] + ziy*fiy[2] + ziz*fiz[2];
+      vxx += xix*fix[0] + xiy*fiy[0] + xiz*fiz[0];
+      vxy += 0.5*(yix*fix[0] + yiy*fiy[0] + yiz*fiz[0] +
+                  xix*fix[1] + xiy*fiy[1] + xiz*fiz[1]);
+      vxz += 0.5*(zix*fix[0] + ziy*fiy[0] + ziz*fiz[0] +
+                  xix*fix[2] + xiy*fiy[2] + xiz*fiz[2]);
+      vyy += yix*fix[1] + yiy*fiy[1] + yiz*fiz[1];
+      vyz += 0.5*(zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
+                  yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
+      vzz += zix*fix[2] + ziy*fiy[2] + ziz*fiz[2];
+    }
   }
 
   // increment total internal virial tensor components
 
-  virmpole[0] -= vxx;
-  virmpole[1] -= vyy;
-  virmpole[2] -= vzz;
-  virmpole[3] -= vxy;
-  virmpole[4] -= vxz;
-  virmpole[5] -= vyz;
+  if (vflag_global) {
+    virmpole[0] -= vxx;
+    virmpole[1] -= vyy;
+    virmpole[2] -= vzz;
+    virmpole[3] -= vxy;
+    virmpole[4] -= vxz;
+    virmpole[5] -= vyz;
+  }
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/AMOEBA/amoeba_polar.cpp b/src/AMOEBA/amoeba_polar.cpp
index e1c90a3f34..75acd402f9 100644
--- a/src/AMOEBA/amoeba_polar.cpp
+++ b/src/AMOEBA/amoeba_polar.cpp
@@ -100,6 +100,8 @@ void PairAmoeba::polar()
 
     torque2force(i,tep,fix,fiy,fiz,f);
 
+    if (!vflag_global) continue;
+
     iz = zaxis2local[i];
     ix = xaxis2local[i];
     iy = yaxis2local[i];
@@ -113,7 +115,7 @@ void PairAmoeba::polar()
     xiy = x[iy][0] - x[i][0];
     yiy = x[iy][1] - x[i][1];
     ziy = x[iy][2] - x[i][2];
-
+    
     vxx = xix*fix[0] + xiy*fiy[0] + xiz*fiz[0];
     vyy = yix*fix[1] + yiy*fiy[1] + yiz*fiz[1];
     vzz = zix*fix[2] + ziy*fiy[2] + ziz*fiz[2];
@@ -124,7 +126,7 @@ void PairAmoeba::polar()
                  xix*fix[2] + xiy*fiy[2] + xiz*fiz[2]);
     vyz = 0.5 * (zix*fix[1] + ziy*fiy[1] + ziz*fiz[1] +
                  yix*fix[2] + yiy*fiy[2] + yiz*fiz[2]);
-
+    
     virpolar[0] -= vxx;
     virpolar[1] -= vyy;
     virpolar[2] -= vzz;
@@ -1139,25 +1141,20 @@ void PairAmoeba::polar_real()
 
       // increment the virial due to pairwise Cartesian forces
 
-      vxx = xr * frcx;
-      vxy = 0.5 * (yr*frcx+xr*frcy);
-      vxz = 0.5 * (zr*frcx+xr*frcz);
-      vyy = yr * frcy;
-      vyz = 0.5 * (zr*frcy+yr*frcz);
-      vzz = zr * frcz;
+      if (vflag_global) {
+        vxx = xr * frcx;
+        vxy = 0.5 * (yr*frcx+xr*frcy);
+        vxz = 0.5 * (zr*frcx+xr*frcz);
+        vyy = yr * frcy;
+        vyz = 0.5 * (zr*frcy+yr*frcz);
+        vzz = zr * frcz;
 
-      virpolar[0] -= vxx;
-      virpolar[1] -= vyy;
-      virpolar[2] -= vzz;
-      virpolar[3] -= vxy;
-      virpolar[4] -= vxz;
-      virpolar[5] -= vyz;
-
-      // energy = e
-      // virial = 6-vec vir
-      // NOTE: add tally function
-
-      if (evflag) {
+        virpolar[0] -= vxx;
+        virpolar[1] -= vyy;
+        virpolar[2] -= vzz;
+        virpolar[3] -= vxy;
+        virpolar[4] -= vxz;
+        virpolar[5] -= vyz;
       }
     }
   }
@@ -1191,6 +1188,8 @@ void PairAmoeba::polar_real()
 
     torque2force(i,tep,fix,fiy,fiz,f);
 
+    if (!vflag_global) continue;
+
     iz = zaxis2local[i];
     ix = xaxis2local[i];
     iy = yaxis2local[i];
@@ -2163,10 +2162,12 @@ void PairAmoeba::polar_kspace()
 
   // increment the total internal virial tensor components
 
-  virpolar[0] -= vxx;
-  virpolar[1] -= vyy;
-  virpolar[2] -= vzz;
-  virpolar[3] -= vxy;
-  virpolar[4] -= vxz;
-  virpolar[5] -= vyz;
+  if (vflag_global) {
+    virpolar[0] -= vxx;
+    virpolar[1] -= vyy;
+    virpolar[2] -= vzz;
+    virpolar[3] -= vxy;
+    virpolar[4] -= vxz;
+    virpolar[5] -= vyz;
+  }
 }
diff --git a/src/AMOEBA/amoeba_repulsion.cpp b/src/AMOEBA/amoeba_repulsion.cpp
index c2bd490715..632996c264 100644
--- a/src/AMOEBA/amoeba_repulsion.cpp
+++ b/src/AMOEBA/amoeba_repulsion.cpp
@@ -354,32 +354,24 @@ void PairAmoeba::repulsion()
 
       // increment the virial due to pairwise Cartesian forces
 
-      vxx = -xr * frcx;
-      vxy = -0.5 * (yr*frcx+xr*frcy);
-      vxz = -0.5 * (zr*frcx+xr*frcz);
-      vyy = -yr * frcy;
-      vyz = -0.5 * (zr*frcy+yr*frcz);
-      vzz = -zr * frcz;
+      if (vflag_global) {
+        vxx = -xr * frcx;
+        vxy = -0.5 * (yr*frcx+xr*frcy);
+        vxz = -0.5 * (zr*frcx+xr*frcz);
+        vyy = -yr * frcy;
+        vyz = -0.5 * (zr*frcy+yr*frcz);
+        vzz = -zr * frcz;
 
-      virrepulse[0] -= vxx;
-      virrepulse[1] -= vyy;
-      virrepulse[2] -= vzz;
-      virrepulse[3] -= vxy;
-      virrepulse[4] -= vxz;
-      virrepulse[5] -= vyz;
-
-      // energy = e
-      // virial = 6-vec vir
-      // NOTE: add tally function
-
-      if (evflag) {
+        virrepulse[0] -= vxx;
+        virrepulse[1] -= vyy;
+        virrepulse[2] -= vzz;
+        virrepulse[3] -= vxy;
+        virrepulse[4] -= vxz;
+        virrepulse[5] -= vyz;
       }
     }
   }
 
-  // DEBUG
-  //fclose(fp);
-
   // reverse comm to sum torque from ghost atoms to owned atoms
 
   crstyle = TORQUE;
@@ -390,6 +382,8 @@ void PairAmoeba::repulsion()
   for (i = 0; i < nlocal; i++) {
     torque2force(i,tq[i],fix,fiy,fiz,f);
 
+    if (!vflag_global) continue;
+
     iz = zaxis2local[i];
     ix = xaxis2local[i];
     iy = yaxis2local[i];
@@ -420,12 +414,6 @@ void PairAmoeba::repulsion()
     virrepulse[3] -= vxy;
     virrepulse[4] -= vxz;
     virrepulse[5] -= vyz;
-
-    // virial = 6-vec vir
-    // NOTE: add tally function
-
-    if (evflag) {
-    }
   }
 }
 
diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp
index 21feb94af8..74013fd5e0 100644
--- a/src/AMOEBA/pair_amoeba.cpp
+++ b/src/AMOEBA/pair_amoeba.cpp
@@ -49,8 +49,6 @@ enum{GEAR,ASPC,LSQR};
 
 #define DELTASTACK 16
 
-#define UIND_DEBUG 0       // also in amoeba_induce.cpp
-
 /* ---------------------------------------------------------------------- */
 
 PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
@@ -67,10 +65,16 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
   me = comm->me;
   nprocs = comm->nprocs;
 
-  // force field settings
+  // pair style settings
 
   one_coeff = 1;
   single_enable = 0;
+  no_virial_fdotr_compute = 1;
+
+  nextra = 6;
+  pvector = new double[nextra];
+
+  // force field settings
 
   amoeba = 1;
   hippo = 0;
@@ -165,6 +169,8 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
 
 PairAmoeba::~PairAmoeba()
 {
+  delete [] pvector;
+
   // check nfix in case all fixes have already been deleted
 
   if (modify->nfix) {
@@ -275,26 +281,12 @@ PairAmoeba::~PairAmoeba()
   }
 
   delete [] factors;
-
-  // DEBUG
-
-  if (me == 0 && UIND_DEBUG) fclose(fp_uind);
 }
 
 /* ---------------------------------------------------------------------- */
 
 void PairAmoeba::compute(int eflag, int vflag)
 {
-  // DEBUG timer init
-
-  if (update->ntimestep <= 1) {
-    time_init = time_hal = time_repulse = time_disp = time_mpole = 0.0;
-    time_induce = time_polar = time_qxfer = 0.0;
-  }
-
-  double evdwl;
-
-  evdwl = 0.0;
   ev_init(eflag,vflag);
 
   // zero energy/virial components
@@ -358,6 +350,13 @@ void PairAmoeba::compute(int eflag, int vflag)
   // end of one-time initializations
   // -------------------------------------------------------------------
 
+  // initialize timers on first compute() call after setup
+
+  if (update->ntimestep <= update->beginstep+1) {
+    time_init = time_hal = time_repulse = time_disp = time_mpole = 0.0;
+    time_induce = time_polar = time_qxfer = 0.0;
+  }
+
   double time0,time1,time2,time3,time4,time5,time6,time7,time8;
 
   MPI_Barrier(world);
@@ -462,29 +461,27 @@ void PairAmoeba::compute(int eflag, int vflag)
   // charge transfer, pairwise
 
   if (hippo && qxfer_flag) charge_transfer();
-
   time8 = MPI_Wtime();
 
-  // energy, force, virial summations
+  // store energy components for output by compute pair command
 
-  eng_vdwl = ehal + erepulse + edisp + empole + epolar + eqxfer;
+  pvector[0] = ehal;
+  pvector[2] = erepulse;
+  pvector[3] = edisp;
+  pvector[4] = empole;
+  pvector[5] = epolar;
+  pvector[6] = eqxfer;
 
-  double **f = atom->f;
-  nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
+  // energy & virial summations
+
+  eng_vdwl = ehal + edisp;
+  eng_coul = erepulse + empole + epolar + eqxfer;
 
   for (int i = 0; i < 6; i++)
     virial[i] = virhal[i] + virrepulse[i] + virdisp[i] +
       virpolar[i] + virmpole[i] + virqxfer[i];
 
-  // virial computation
-  // NOTE: how does this work for AMOEBA ?
-  // do all terms get summed this way, or only pairwise
-  // it is 2x what summed virial[6] above is
-
-  // if (vflag_fdotr) virial_fdotr_compute();
-
-  // timing information
+  // accumulate timing information
 
   time_init    += time1 - time0;
   time_hal     += time2 - time1;
@@ -494,15 +491,18 @@ void PairAmoeba::compute(int eflag, int vflag)
   time_induce  += time6 - time5;
   time_polar   += time7 - time6;
   time_qxfer   += time8 - time7;
+}
 
-  // timing output
-
-  if (update->ntimestep < update->laststep) return;
+/* ----------------------------------------------------------------------
+   print out AMOEBA/HIPPO timing info at end of run
+------------------------------------------------------------------------- */
 
+void PairAmoeba::finish()
+{
   double ave;
   MPI_Allreduce(&time_init,&ave,1,MPI_DOUBLE,MPI_SUM,world);
   time_init = ave/nprocs;
-
+  
   MPI_Allreduce(&time_hal,&ave,1,MPI_DOUBLE,MPI_SUM,world);
   time_hal = ave/nprocs;
 
@@ -524,32 +524,32 @@ void PairAmoeba::compute(int eflag, int vflag)
   MPI_Allreduce(&time_qxfer,&ave,1,MPI_DOUBLE,MPI_SUM,world);
   time_qxfer = ave/nprocs;
 
-  double time_amtot = time_init + time_hal + time_repulse + time_disp +
+  double time_total = time_init + time_hal + time_repulse + time_disp +
     time_mpole + time_induce + time_polar + time_qxfer;
 
   if (me == 0) {
-    utils::logmesg(lmp,"\nAMEOBA/HIPPO timing info:\n");
-    utils::logmesg(lmp,"  Init   time: {:.6g} {:.6g}\n",
-                   time_init,time_init/time_amtot);
+    utils::logmesg(lmp,"\nAMEOBA/HIPPO timing breakdown:\n");
+    utils::logmesg(lmp,"  Init    time: {:.6g} {:.6g}\n",
+                   time_init,time_init/time_total);
     if (amoeba)
-      utils::logmesg(lmp,"  Hal    time: {:.6g} {:.6g}\n",
-                     time_hal,time_hal/time_amtot*100);
+      utils::logmesg(lmp,"  Hal     time: {:.6g} {:.6g}\n",
+                     time_hal,time_hal/time_total*100);
     if (hippo)
-      utils::logmesg(lmp,"  Repls  time: {:.6g} {:.6g}\n",
-                     time_repulse,time_repulse/time_amtot*100);
+      utils::logmesg(lmp,"  Repulse time: {:.6g} {:.6g}\n",
+                     time_repulse,time_repulse/time_total*100);
     if (hippo)
-      utils::logmesg(lmp,"  Disp   time: {:.6g} {:.6g}\n",
-                     time_disp,time_disp/time_amtot*100);
-    utils::logmesg(lmp,"  Mpole  time: {:.6g} {:.6g}\n",
-                   time_mpole,time_mpole/time_amtot*100);
-    utils::logmesg(lmp,"  Induce time: {:.6g} {:.6g}\n",
-                   time_induce,time_induce/time_amtot*100);
-    utils::logmesg(lmp,"  Polar  time: {:.6g} {:.6g}\n",
-                   time_polar,time_polar/time_amtot*100);
+      utils::logmesg(lmp,"  Disp    time: {:.6g} {:.6g}\n",
+                     time_disp,time_disp/time_total*100);
+    utils::logmesg(lmp,"  Mpole   time: {:.6g} {:.6g}\n",
+                   time_mpole,time_mpole/time_total*100);
+    utils::logmesg(lmp,"  Induce  time: {:.6g} {:.6g}\n",
+                   time_induce,time_induce/time_total*100);
+    utils::logmesg(lmp,"  Polar   time: {:.6g} {:.6g}\n",
+                   time_polar,time_polar/time_total*100);
     if (hippo)
-      utils::logmesg(lmp,"  Qxfer  time: {:.6g} {:.6g}\n",
-                     time_qxfer,time_qxfer/time_amtot*100);
-    utils::logmesg(lmp,"  Total  time: {:.6g}\n\n",time_amtot);
+      utils::logmesg(lmp,"  Qxfer   time: {:.6g} {:.6g}\n",
+                     time_qxfer,time_qxfer/time_total*100);
+    utils::logmesg(lmp,"  Total   time: {:.6g}\n",time_total);
   }
 }
 
@@ -1016,15 +1016,6 @@ void PairAmoeba::init_style()
   // request neighbor lists
 
   int irequest = neighbor->request(this,instance_me);
-
-  // open debug output files
-  // names are hard-coded
-
-  if (me == 0) {
-    char fname[32];
-    sprintf(fname,"tmp.uind.kspace.%d",nprocs);
-    if (UIND_DEBUG) fp_uind = fopen(fname,"w");
-  }
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h
index f557f81772..21811d2610 100644
--- a/src/AMOEBA/pair_amoeba.h
+++ b/src/AMOEBA/pair_amoeba.h
@@ -46,6 +46,7 @@ class PairAmoeba : public Pair {
   void coeff(int, char **);
   void init_style();
   double init_one(int, int);
+  void finish();
 
   int pack_forward_comm(int, int *, double *, int, int *);
   void unpack_forward_comm(int, int, double *);
diff --git a/src/finish.cpp b/src/finish.cpp
index ea472682a2..b74bcf3a0e 100644
--- a/src/finish.cpp
+++ b/src/finish.cpp
@@ -28,6 +28,7 @@
 #include "neigh_request.h"
 #include "neighbor.h"           // IWYU pragma: keep
 #include "output.h"
+#include "pair.h"
 #include "thermo.h"
 #include "timer.h"              // IWYU pragma: keep
 #include "universe.h"
@@ -214,6 +215,10 @@ void Finish::end(int flag)
     }
   }
 
+  // pair_style timing stats if provided
+
+  if (force->pair) force->pair->finish();
+
   // PRD stats
 
   if (prdflag) {
diff --git a/src/pair.h b/src/pair.h
index 8c1c8340eb..cb36009af1 100644
--- a/src/pair.h
+++ b/src/pair.h
@@ -169,6 +169,7 @@ class Pair : protected Pointers {
     return 0.0;
   }
 
+  virtual void finish() {}
   virtual void settings(int, char **) = 0;
   virtual void coeff(int, char **) = 0;
 
diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp
index 82dce4d5f2..5977a8b101 100644
--- a/src/pair_hybrid.cpp
+++ b/src/pair_hybrid.cpp
@@ -192,6 +192,12 @@ void PairHybrid::compute(int eflag, int vflag)
   if (vflag_fdotr) virial_fdotr_compute();
 }
 
+/* ---------------------------------------------------------------------- */
+
+void PairHybrid::finish()
+{
+  for (int m = 0; m < nstyles; m++) styles[m]->finish();
+}
 
 /* ---------------------------------------------------------------------- */
 
diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h
index ee56783700..358e8d28d0 100644
--- a/src/pair_hybrid.h
+++ b/src/pair_hybrid.h
@@ -46,6 +46,7 @@ class PairHybrid : public Pair {
   void init_style() override;
   double init_one(int, int) override;
   void setup() override;
+  void finish();
   void write_restart(FILE *) override;
   void read_restart(FILE *) override;
   double single(int, int, int, int, double, double, double, double &) override;