diff --git a/src/MOLECULE/bond_fene.cpp b/src/MOLECULE/bond_fene.cpp
index 41d81f5e17..d70a76734e 100644
--- a/src/MOLECULE/bond_fene.cpp
+++ b/src/MOLECULE/bond_fene.cpp
@@ -231,7 +231,8 @@ void BondFENE::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double BondFENE::single(int type, double rsq, int i, int j)
+double BondFENE::single(int type, double rsq, int i, int j, 
+                        double &fforce)
 {
   double r0sq = r0[type] * r0[type];
   double rlogarg = 1.0 - rsq/r0sq;
@@ -250,11 +251,13 @@ double BondFENE::single(int type, double rsq, int i, int j)
   }
 
   double eng = -0.5 * k[type]*r0sq*log(rlogarg);
+  fforce = -k[type]/rlogarg;
   if (rsq < TWO_1_3*sigma[type]*sigma[type]) {
     double sr2,sr6;
     sr2 = sigma[type]*sigma[type]/rsq;
     sr6 = sr2*sr2*sr2;
     eng += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
+    fforce += 48.0*epsilon[type]*sr6*(sr6-0.5)/rsq;
   }
 
   return eng;
diff --git a/src/MOLECULE/bond_fene.h b/src/MOLECULE/bond_fene.h
index 55216b05ad..34c9e37db4 100644
--- a/src/MOLECULE/bond_fene.h
+++ b/src/MOLECULE/bond_fene.h
@@ -35,7 +35,7 @@ class BondFENE : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  double single(int, double, int, int);
+  double single(int, double, int, int, double &);
 
  protected:
   double TWO_1_3;
diff --git a/src/MOLECULE/bond_fene_expand.cpp b/src/MOLECULE/bond_fene_expand.cpp
index a8051468ac..dc9bfa9763 100644
--- a/src/MOLECULE/bond_fene_expand.cpp
+++ b/src/MOLECULE/bond_fene_expand.cpp
@@ -242,7 +242,8 @@ void BondFENEExpand::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double BondFENEExpand::single(int type, double rsq, int i, int j)
+double BondFENEExpand::single(int type, double rsq, int i, int j,
+                        double &fforce)
 {
   double r = sqrt(rsq);
   double rshift = r - shift[type];
@@ -264,11 +265,13 @@ double BondFENEExpand::single(int type, double rsq, int i, int j)
   }
 
   double eng = -0.5 * k[type]*r0sq*log(rlogarg);
+  fforce = -k[type]*rshift/rlogarg/r;
   if (rshiftsq < TWO_1_3*sigma[type]*sigma[type]) {
     double sr2,sr6;
     sr2 = sigma[type]*sigma[type]/rshiftsq;
     sr6 = sr2*sr2*sr2;
     eng += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
+    fforce += 48.0*epsilon[type]*sr6*(sr6-0.5)/rshift/r;
   }
 
   return eng;
diff --git a/src/MOLECULE/bond_fene_expand.h b/src/MOLECULE/bond_fene_expand.h
index 1d8e3a3ede..9177f415c4 100644
--- a/src/MOLECULE/bond_fene_expand.h
+++ b/src/MOLECULE/bond_fene_expand.h
@@ -35,7 +35,7 @@ class BondFENEExpand : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  double single(int, double, int, int);
+  double single(int, double, int, int, double &);
 
  protected:
   double TWO_1_3;
diff --git a/src/MOLECULE/bond_harmonic.cpp b/src/MOLECULE/bond_harmonic.cpp
index 4fb18d13ca..1de76d90b9 100644
--- a/src/MOLECULE/bond_harmonic.cpp
+++ b/src/MOLECULE/bond_harmonic.cpp
@@ -176,10 +176,13 @@ void BondHarmonic::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double BondHarmonic::single(int type, double rsq, int i, int j)
+double BondHarmonic::single(int type, double rsq, int i, int j,
+                        double &fforce)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double rk = k[type] * dr;
+  fforce = 0;
+  if (r > 0.0) fforce = -2.0*rk/r;
   return rk*dr;
 }
diff --git a/src/MOLECULE/bond_harmonic.h b/src/MOLECULE/bond_harmonic.h
index 90bc0e1fb7..173ca5edbc 100644
--- a/src/MOLECULE/bond_harmonic.h
+++ b/src/MOLECULE/bond_harmonic.h
@@ -34,7 +34,7 @@ class BondHarmonic : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  double single(int, double, int, int);
+  double single(int, double, int, int, double &);
 
  protected:
   double *k,*r0;
diff --git a/src/MOLECULE/bond_morse.cpp b/src/MOLECULE/bond_morse.cpp
index 22d2909244..d1e8ffe02d 100644
--- a/src/MOLECULE/bond_morse.cpp
+++ b/src/MOLECULE/bond_morse.cpp
@@ -186,10 +186,13 @@ void BondMorse::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double BondMorse::single(int type, double rsq, int i, int j)
+double BondMorse::single(int type, double rsq, int i, int j,
+                         double &fforce)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double ralpha = exp(-alpha[type]*dr);
+  fforce = 0;
+  if (r > 0.0) fforce = -2.0*d0[type]*alpha[type]*(1-ralpha)*ralpha/r;
   return d0[type]*(1-ralpha)*(1-ralpha);
 }
diff --git a/src/MOLECULE/bond_morse.h b/src/MOLECULE/bond_morse.h
index 87303e4bc7..8e24480eef 100644
--- a/src/MOLECULE/bond_morse.h
+++ b/src/MOLECULE/bond_morse.h
@@ -34,7 +34,7 @@ class BondMorse : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  double single(int, double, int, int);
+  double single(int, double, int, int, double &);
 
  protected:
   double *d0,*alpha,*r0;
diff --git a/src/MOLECULE/bond_nonlinear.cpp b/src/MOLECULE/bond_nonlinear.cpp
index 6577d4380c..b2c689fd5f 100644
--- a/src/MOLECULE/bond_nonlinear.cpp
+++ b/src/MOLECULE/bond_nonlinear.cpp
@@ -181,12 +181,15 @@ void BondNonlinear::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double BondNonlinear::single(int type, double rsq, int i, int j)
+double BondNonlinear::single(int type, double rsq, int i, int j,
+                             double &fforce)
 {
   double r = sqrt(rsq);
   double dr = r - r0[type];
   double drsq = dr*dr;
   double lamdasq = lamda[type]*lamda[type];
   double denom = lamdasq - drsq;
+  double denomsq = denom*denom;
+  fforce = -epsilon[type]/r * 2.0*dr*lamdasq/denomsq;
   return epsilon[type] * drsq / denom;
 }
diff --git a/src/MOLECULE/bond_nonlinear.h b/src/MOLECULE/bond_nonlinear.h
index 02efaa969f..0c8a99fc5b 100644
--- a/src/MOLECULE/bond_nonlinear.h
+++ b/src/MOLECULE/bond_nonlinear.h
@@ -34,7 +34,7 @@ class BondNonlinear : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  double single(int, double, int, int);
+  double single(int, double, int, int, double &);
 
  protected:
   double *epsilon,*r0,*lamda;
diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp
index 9c49a5e510..78e97d986b 100755
--- a/src/MOLECULE/bond_quartic.cpp
+++ b/src/MOLECULE/bond_quartic.cpp
@@ -296,7 +296,8 @@ void BondQuartic::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double BondQuartic::single(int type, double rsq, int i, int j)
+double BondQuartic::single(int type, double rsq, int i, int j,
+                           double &fforce)
 {
   double r,dr,r2,ra,rb,sr2,sr6;
 
@@ -325,10 +326,13 @@ double BondQuartic::single(int type, double rsq, int i, int j)
   rb = dr - b2[type];
 
   eng += k[type]*r2*ra*rb + u0[type];
+  fforce = -k[type]/r * (r2*(ra+rb) + 2.0*dr*ra*rb);
+
   if (rsq < TWO_1_3) {
     sr2 = 1.0/rsq;
     sr6 = sr2*sr2*sr2;
     eng += 4.0*sr6*(sr6-1.0) + 1.0;
+    fforce += 48.0*sr6*(sr6-0.5)/rsq;
   }
 
   return eng;
diff --git a/src/MOLECULE/bond_quartic.h b/src/MOLECULE/bond_quartic.h
index 4fb144c655..1507b4cf98 100644
--- a/src/MOLECULE/bond_quartic.h
+++ b/src/MOLECULE/bond_quartic.h
@@ -35,7 +35,7 @@ class BondQuartic : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  double single(int, double, int, int);
+  double single(int, double, int, int, double &);
 
  protected:
   double TWO_1_3;
diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp
index dc6a637d8b..b5e7a6a8d3 100644
--- a/src/MOLECULE/bond_table.cpp
+++ b/src/MOLECULE/bond_table.cpp
@@ -241,11 +241,14 @@ void BondTable::read_restart(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double BondTable::single(int type, double rsq, int i, int j)
+double BondTable::single(int type, double rsq, int i, int j,
+                         double &fforce)
 {
   double r = sqrt(rsq);
   double u;
-  u_lookup(type,r,u);
+  double mdu;
+  uf_lookup(type,r,u,mdu);
+  fforce = mdu/r;
   return u;
 }
 
diff --git a/src/MOLECULE/bond_table.h b/src/MOLECULE/bond_table.h
index 3a6f75d89b..17db413c87 100644
--- a/src/MOLECULE/bond_table.h
+++ b/src/MOLECULE/bond_table.h
@@ -35,7 +35,7 @@ class BondTable : public Bond {
   double equilibrium_distance(int);
   void write_restart(FILE *);
   void read_restart(FILE *);
-  double single(int, double, int, int);
+  double single(int, double, int, int, double &);
 
  protected:
   int tabstyle,tablength;