diff --git a/src/MOLECULE/improper_umbrella.cpp b/src/MOLECULE/improper_umbrella.cpp
index 31d7cba54dca313ae2016fe477b9882e0787f515..8ed3d3a84cd6437a6c0894f9d7140c231234aa75 100644
--- a/src/MOLECULE/improper_umbrella.cpp
+++ b/src/MOLECULE/improper_umbrella.cpp
@@ -189,17 +189,17 @@ void ImproperUmbrella::compute(int eflag, int vflag)
     dahy = ary-c*hry;
     dahz = arz-c*hrz;
 
-    f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
-    f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
-    f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
+    f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
+    f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
+    f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
 
-    f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
-    f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
-    f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
+    f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
+    f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
+    f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
 
-    f4[0] = dahx*rhr;
-    f4[1] = dahy*rhr;
-    f4[2] = dahz*rhr;
+    f4[0] = dahx*rhr*a;
+    f4[1] = dahy*rhr*a;
+    f4[2] = dahz*rhr*a;
 
     f1[0] = -(f2[0] + f3[0] + f4[0]);
     f1[1] = -(f2[1] + f3[1] + f4[1]);
@@ -208,27 +208,27 @@ void ImproperUmbrella::compute(int eflag, int vflag)
     // apply force to each of 4 atoms
 
     if (newton_bond || i1 < nlocal) {
-      f[i1][0] += f1[0]*a;
-      f[i1][1] += f1[1]*a;
-      f[i1][2] += f1[2]*a;
+      f[i1][0] += f1[0];
+      f[i1][1] += f1[1];
+      f[i1][2] += f1[2];
     }
 
     if (newton_bond || i2 < nlocal) {
-      f[i2][0] += f3[0]*a;
-      f[i2][1] += f3[1]*a;
-      f[i2][2] += f3[2]*a;
+      f[i2][0] += f3[0];
+      f[i2][1] += f3[1];
+      f[i2][2] += f3[2];
     }
 
     if (newton_bond || i3 < nlocal) {
-      f[i3][0] += f2[0]*a;
-      f[i3][1] += f2[1]*a;
-      f[i3][2] += f2[2]*a;
+      f[i3][0] += f2[0];
+      f[i3][1] += f2[1];
+      f[i3][2] += f2[2];
     }
 
     if (newton_bond || i4 < nlocal) {
-      f[i4][0] += f4[0]*a;
-      f[i4][1] += f4[1]*a;
-      f[i4][2] += f4[2]*a;
+      f[i4][0] += f4[0];
+      f[i4][1] += f4[1];
+      f[i4][2] += f4[2];
     }
 
     if (evflag) {
@@ -247,7 +247,7 @@ void ImproperUmbrella::compute(int eflag, int vflag)
       vb3y = x[i4][1] - x[i3][1];
       vb3z = x[i4][2] - x[i3][2];
 
-      ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,
+      ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f2,f4,
                vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
     }
   }
diff --git a/src/USER-MISC/improper_fourier.cpp b/src/USER-MISC/improper_fourier.cpp
index 11478d08eab135f70baa9ab6d69a762b4a302672..927478fa1a9245f640e40562d84033001bd2afba 100644
--- a/src/USER-MISC/improper_fourier.cpp
+++ b/src/USER-MISC/improper_fourier.cpp
@@ -206,17 +206,17 @@ void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int
   dahy = ary-c*hry;
   dahz = arz-c*hrz;
 
-  f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
-  f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
-  f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
+  f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
+  f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
+  f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
 
-  f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
-  f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
-  f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
+  f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
+  f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
+  f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
 
-  f4[0] = dahx*rhr;
-  f4[1] = dahy*rhr;
-  f4[2] = dahz*rhr;
+  f4[0] = dahx*rhr*a;
+  f4[1] = dahy*rhr*a;
+  f4[2] = dahz*rhr*a;
 
   f1[0] = -(f2[0] + f3[0] + f4[0]);
   f1[1] = -(f2[1] + f3[1] + f4[1]);
@@ -225,32 +225,32 @@ void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int
   // apply force to each of 4 atoms
 
   if (newton_bond || i1 < nlocal) {
-    f[i1][0] += f1[0]*a;
-    f[i1][1] += f1[1]*a;
-    f[i1][2] += f1[2]*a;
+    f[i1][0] += f1[0];
+    f[i1][1] += f1[1];
+    f[i1][2] += f1[2];
   }
 
   if (newton_bond || i2 < nlocal) {
-    f[i2][0] += f3[0]*a;
-    f[i2][1] += f3[1]*a;
-    f[i2][2] += f3[2]*a;
+    f[i2][0] += f3[0];
+    f[i2][1] += f3[1];
+    f[i2][2] += f3[2];
   }
 
   if (newton_bond || i3 < nlocal) {
-    f[i3][0] += f2[0]*a;
-    f[i3][1] += f2[1]*a;
-    f[i3][2] += f2[2]*a;
+    f[i3][0] += f2[0];
+    f[i3][1] += f2[1];
+    f[i3][2] += f2[2];
   }
 
   if (newton_bond || i4 < nlocal) {
-    f[i4][0] += f4[0]*a;
-    f[i4][1] += f4[1]*a;
-    f[i4][2] += f4[2]*a;
+    f[i4][0] += f4[0];
+    f[i4][1] += f4[1];
+    f[i4][2] += f4[2];
   }
 
   if (evflag)
-    ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4,
-             vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
+      ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f2,f4,
+               -vb1x,-vb1y,-vb1z,vb2x-vb1x,vb2y-vb1y,vb2z-vb1z,vb3x-vb2x,vb3y-vb2y,vb3z-vb2z);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/USER-OMP/improper_fourier_omp.cpp b/src/USER-OMP/improper_fourier_omp.cpp
index aed04003a5eff70550f4d98830145a7557cc22b5..77dd36b64f62d014bf3335178585bed290794b37 100644
--- a/src/USER-OMP/improper_fourier_omp.cpp
+++ b/src/USER-OMP/improper_fourier_omp.cpp
@@ -239,17 +239,17 @@ void ImproperFourierOMP::add1_thr(const int i1,const int i2,
   dahy = ary-c*hry;
   dahz = arz-c*hrz;
 
-  f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
-  f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
-  f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
+  f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
+  f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
+  f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
 
-  f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
-  f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
-  f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
+  f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
+  f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
+  f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
 
-  f4[0] = dahx*rhr;
-  f4[1] = dahy*rhr;
-  f4[2] = dahz*rhr;
+  f4[0] = dahx*rhr*a;
+  f4[1] = dahy*rhr*a;
+  f4[2] = dahz*rhr*a;
 
   f1[0] = -(f2[0] + f3[0] + f4[0]);
   f1[1] = -(f2[1] + f3[1] + f4[1]);
@@ -258,30 +258,31 @@ void ImproperFourierOMP::add1_thr(const int i1,const int i2,
   // apply force to each of 4 atoms
 
   if (NEWTON_BOND || i1 < nlocal) {
-    f[i1][0] += f1[0]*a;
-    f[i1][1] += f1[1]*a;
-    f[i1][2] += f1[2]*a;
+    f[i1][0] += f1[0];
+    f[i1][1] += f1[1];
+    f[i1][2] += f1[2];
   }
 
   if (NEWTON_BOND || i2 < nlocal) {
-    f[i2][0] += f3[0]*a;
-    f[i2][1] += f3[1]*a;
-    f[i2][2] += f3[2]*a;
+    f[i2][0] += f3[0];
+    f[i2][1] += f3[1];
+    f[i2][2] += f3[2];
   }
 
   if (NEWTON_BOND || i3 < nlocal) {
-    f[i3][0] += f2[0]*a;
-    f[i3][1] += f2[1]*a;
-    f[i3][2] += f2[2]*a;
+    f[i3][0] += f2[0];
+    f[i3][1] += f2[1];
+    f[i3][2] += f2[2];
   }
 
   if (NEWTON_BOND || i4 < nlocal) {
-    f[i4][0] += f4[0]*a;
-    f[i4][1] += f4[1]*a;
-    f[i4][2] += f4[2]*a;
+    f[i4][0] += f4[0];
+    f[i4][1] += f4[1];
+    f[i4][2] += f4[2];
   }
 
   if (EVFLAG)
-    ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f3,f4,
-                 vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,thr);
+    ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f2,f4,
+                 -vb1x,-vb1y,-vb1z,vb2x-vb1x,vb2y-vb1y,vb2z-vb1z,vb3x-vb2x,vb3y-vb2y,vb3z-vb2z,thr);
+
 }
diff --git a/src/USER-OMP/improper_umbrella_omp.cpp b/src/USER-OMP/improper_umbrella_omp.cpp
index dc11f24a4de084b07b666c08b76271d135218115..69d41176c64287c55a3663f4db259df03bd5b319 100644
--- a/src/USER-OMP/improper_umbrella_omp.cpp
+++ b/src/USER-OMP/improper_umbrella_omp.cpp
@@ -212,17 +212,17 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr)
     dahy = ary-c*hry;
     dahz = arz-c*hrz;
 
-    f2[0] = (dhay*vb1z - dhaz*vb1y)*rar;
-    f2[1] = (dhaz*vb1x - dhax*vb1z)*rar;
-    f2[2] = (dhax*vb1y - dhay*vb1x)*rar;
+    f2[0] = (dhay*vb1z - dhaz*vb1y)*rar*a;
+    f2[1] = (dhaz*vb1x - dhax*vb1z)*rar*a;
+    f2[2] = (dhax*vb1y - dhay*vb1x)*rar*a;
 
-    f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar;
-    f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar;
-    f3[2] = (-dhax*vb2y + dhay*vb2x)*rar;
+    f3[0] = (-dhay*vb2z + dhaz*vb2y)*rar*a;
+    f3[1] = (-dhaz*vb2x + dhax*vb2z)*rar*a;
+    f3[2] = (-dhax*vb2y + dhay*vb2x)*rar*a;
 
-    f4[0] = dahx*rhr;
-    f4[1] = dahy*rhr;
-    f4[2] = dahz*rhr;
+    f4[0] = dahx*rhr*a;
+    f4[1] = dahy*rhr*a;
+    f4[2] = dahz*rhr*a;
 
     f1[0] = -(f2[0] + f3[0] + f4[0]);
     f1[1] = -(f2[1] + f3[1] + f4[1]);
@@ -231,27 +231,27 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr)
     // apply force to each of 4 atoms
 
     if (NEWTON_BOND || i1 < nlocal) {
-      f[i1].x += f1[0]*a;
-      f[i1].y += f1[1]*a;
-      f[i1].z += f1[2]*a;
+      f[i1].x += f1[0];
+      f[i1].y += f1[1];
+      f[i1].z += f1[2];
     }
 
     if (NEWTON_BOND || i2 < nlocal) {
-      f[i2].x += f3[0]*a;
-      f[i2].y += f3[1]*a;
-      f[i2].z += f3[2]*a;
+      f[i2].x += f3[0];
+      f[i2].y += f3[1];
+      f[i2].z += f3[2];
     }
 
     if (NEWTON_BOND || i3 < nlocal) {
-      f[i3].x += f2[0]*a;
-      f[i3].y += f2[1]*a;
-      f[i3].z += f2[2]*a;
+      f[i3].x += f2[0];
+      f[i3].y += f2[1];
+      f[i3].z += f2[2];
     }
 
     if (NEWTON_BOND || i4 < nlocal) {
-      f[i4].x += f4[0]*a;
-      f[i4].y += f4[1]*a;
-      f[i4].z += f4[2]*a;
+      f[i4].x += f4[0];
+      f[i4].y += f4[1];
+      f[i4].z += f4[2];
     }
 
     if (EVFLAG) {
@@ -270,7 +270,7 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr)
       vb3y = x[i4].y - x[i3].y;
       vb3z = x[i4].z - x[i3].z;
 
-      ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f3,f4,
+      ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f2,f4,
                    vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,thr);
     }
   }