diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp
index 575836fe6d..381629e161 100644
--- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp
+++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp
@@ -79,8 +79,8 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
 {
   int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype;
   double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
-  double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond;
-  double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
+  double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch;
+  double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2,d;
   double fi[3],fj[3],delr1[3],delr2[3];
   double r2inv,r10inv;
   double switch1,switch2;
@@ -175,14 +175,20 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
               powint(c,pm.ap-1)*s;
 
             eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4);
+
+            force_switch=0.0;
+
             if (rsq > pm.cut_innersq) {
               switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
                         (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) /
                         pm.denom_vdw;
               switch2 = 12.0*rsq * (pm.cut_outersq-rsq) *
                         (rsq-pm.cut_innersq) / pm.denom_vdw;
-              force_kernel = force_kernel*switch1 + eng_lj*switch2;
-              eng_lj *= switch1;
+
+              force_kernel *= switch1;
+              force_angle  *= switch1;
+              force_switch  = eng_lj*switch2/rsq;
+              eng_lj       *= switch1;
             }
 
             if (eflag) {
@@ -193,6 +199,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
 
             a = factor_hb*force_angle/s;
             b = factor_hb*force_kernel;
+            d = factor_hb*force_switch;
 
             a11 = a*c / rsq1;
             a12 = -a / (r1*r2);
@@ -205,12 +212,12 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag)
             vz1 = a11*delr1[2] + a12*delr2[2];
             vz2 = a22*delr2[2] + a12*delr1[2];
 
-            fi[0] = vx1 + b*delx;
-            fi[1] = vy1 + b*dely;
-            fi[2] = vz1 + b*delz;
-            fj[0] = vx2 - b*delx;
-            fj[1] = vy2 - b*dely;
-            fj[2] = vz2 - b*delz;
+            fi[0] = vx1 + b*delx + d*delx;
+            fi[1] = vy1 + b*dely + d*dely;
+            fi[2] = vz1 + b*delz + d*delz;
+            fj[0] = vx2 - b*delx - d*delx;
+            fj[1] = vy2 - b*dely - d*dely;
+            fj[2] = vz2 - b*delz - d*delz;
 
             f[i][0] += fi[0];
             f[i][1] += fi[1];
diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp
index 02d6bbba0d..e3c26a5dc2 100644
--- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp
+++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp
@@ -50,8 +50,8 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
 {
   int i,j,k,m,ii,jj,kk,inum,jnum,knum,itype,jtype,ktype;
   double delx,dely,delz,rsq,rsq1,rsq2,r1,r2;
-  double factor_hb,force_angle,force_kernel,evdwl,ehbond;
-  double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
+  double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond;
+  double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2;
   double fi[3],fj[3],delr1[3],delr2[3];
   double r,dr,dexp,eng_morse,switch1,switch2;
   int *ilist,*jlist,*klist,*numneigh,**firstneigh;
@@ -143,6 +143,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
             eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp);
             force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap);
             force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s;
+	    force_switch = 0.0;
 
             if (rsq > pm.cut_innersq) {
               switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) *
@@ -150,8 +151,11 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
                         pm.denom_vdw;
               switch2 = 12.0*rsq * (pm.cut_outersq-rsq) *
                         (rsq-pm.cut_innersq) / pm.denom_vdw;
-              force_kernel = force_kernel*switch1 + eng_morse*switch2;
-              eng_morse *= switch1;
+
+              force_kernel *= switch1;
+	      force_angle  *= switch1;
+	      force_switch  = eng_morse*switch2/rsq;
+              eng_morse    *= switch1;
             }
 
             if (eflag) {
@@ -162,6 +166,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
 
             a = factor_hb*force_angle/s;
             b = factor_hb*force_kernel;
+            d = factor_hb*force_switch;
 
             a11 = a*c / rsq1;
             a12 = -a / (r1*r2);
@@ -174,12 +179,12 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag)
             vz1 = a11*delr1[2] + a12*delr2[2];
             vz2 = a22*delr2[2] + a12*delr1[2];
 
-            fi[0] = vx1 + b*delx;
-            fi[1] = vy1 + b*dely;
-            fi[2] = vz1 + b*delz;
-            fj[0] = vx2 - b*delx;
-            fj[1] = vy2 - b*dely;
-            fj[2] = vz2 - b*delz;
+            fi[0] = vx1 + (b+d)*delx;
+            fi[1] = vy1 + (b+d)*dely;
+            fi[2] = vz1 + (b+d)*delz;
+            fj[0] = vx2 - (b+d)*delx;
+            fj[1] = vy2 - (b+d)*dely;
+            fj[2] = vz2 - (b+d)*delz;
 
             f[i][0] += fi[0];
             f[i][1] += fi[1];