diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp
index 2d5cb0652b..acc162ecbf 100644
--- a/src/BPM/bond_bpm_rotational.cpp
+++ b/src/BPM/bond_bpm_rotational.cpp
@@ -454,14 +454,12 @@ void BondBPMRotational::damping_forces(int i1, int i2, int type, double *rhat, d
 
 void BondBPMRotational::compute(int eflag, int vflag)
 {
-
   if (!fix_bond_history->stored_flag) {
     fix_bond_history->stored_flag = true;
     store_data();
   }
 
-  if (hybrid_flag)
-    fix_bond_history->compress_history();
+  if (hybrid_flag) fix_bond_history->compress_history();
 
   int i1, i2, itmp, n, type;
   double r[3], r0[3], rhat[3];
@@ -544,33 +542,36 @@ void BondBPMRotational::compute(int eflag, int vflag)
     //  Apply forces and torques to particles
     // ------------------------------------------------------//
 
-    if (newton_bond || i1 < nlocal) {
-      f[i1][0] -= force1on2[0] * smooth;
-      f[i1][1] -= force1on2[1] * smooth;
-      f[i1][2] -= force1on2[2] * smooth;
+    MathExtra::scale3(smooth, force1on2);
+    MathExtra::scale3(smooth, torque2on1);
 
-      torque[i1][0] += torque2on1[0] * smooth;
-      torque[i1][1] += torque2on1[1] * smooth;
-      torque[i1][2] += torque2on1[2] * smooth;
+    if (newton_bond || i1 < nlocal) {
+      f[i1][0] -= force1on2[0];
+      f[i1][1] -= force1on2[1];
+      f[i1][2] -= force1on2[2];
+
+      torque[i1][0] += torque2on1[0];
+      torque[i1][1] += torque2on1[1];
+      torque[i1][2] += torque2on1[2];
     }
 
     if (newton_bond || i2 < nlocal) {
-      f[i2][0] += force1on2[0] * smooth;
-      f[i2][1] += force1on2[1] * smooth;
-      f[i2][2] += force1on2[2] * smooth;
+      f[i2][0] += force1on2[0];
+      f[i2][1] += force1on2[1];
+      f[i2][2] += force1on2[2];
 
-      torque[i2][0] += torque1on2[0] * smooth;
-      torque[i2][1] += torque1on2[1] * smooth;
-      torque[i2][2] += torque1on2[2] * smooth;
+      MathExtra::scale3(smooth, torque1on2);
+      torque[i2][0] += torque1on2[0];
+      torque[i2][1] += torque1on2[1];
+      torque[i2][2] += torque1on2[2];
     }
 
     if (evflag)
-      ev_tally_xyz(i1, i2, nlocal, newton_bond, 0.0, -force1on2[0] * smooth, -force1on2[1] * smooth,
-                   -force1on2[2] * smooth, r[0], r[1], r[2]);
+      ev_tally_xyz(i1, i2, nlocal, newton_bond, 0.0, -force1on2[0], -force1on2[1],
+                   -force1on2[2], r[0], r[1], r[2]);
   }
 
-  if (hybrid_flag)
-    fix_bond_history->uncompress_history();
+  if (hybrid_flag) fix_bond_history->uncompress_history();
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/BPM/bond_bpm_spring.cpp b/src/BPM/bond_bpm_spring.cpp
index 37437d5b60..0c61c20bf6 100644
--- a/src/BPM/bond_bpm_spring.cpp
+++ b/src/BPM/bond_bpm_spring.cpp
@@ -149,8 +149,6 @@ void BondBPMSpring::store_data()
 
 void BondBPMSpring::compute(int eflag, int vflag)
 {
-  if (hybrid_flag) fix_bond_history->compress_history();
-
   int i, bond_change_flag;
   double *vol0, *vol;
 
@@ -195,6 +193,8 @@ void BondBPMSpring::compute(int eflag, int vflag)
     }
   }
 
+  if (hybrid_flag) fix_bond_history->compress_history();
+
   int i1, i2, itmp, n, type;
   double delx, dely, delz, delvx, delvy, delvz;
   double e, rsq, r, r0, rinv, smooth, fbond, dot;