diff --git a/src/MOLECULE/bond_fene.cpp b/src/MOLECULE/bond_fene.cpp
index 0f279fccb6..7f55c89296 100644
--- a/src/MOLECULE/bond_fene.cpp
+++ b/src/MOLECULE/bond_fene.cpp
@@ -30,6 +30,13 @@ using MathConst::MY_CUBEROOT2;
 
 /* ---------------------------------------------------------------------- */
 
+BondFENE::BondFENE(LAMMPS *_lmp) : Bond(_lmp)
+{
+  born_matrix_enable = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
 BondFENE::~BondFENE()
 {
   if (allocated && !copymode) {
@@ -262,6 +269,28 @@ double BondFENE::single(int type, double rsq, int /*i*/, int /*j*/, double &ffor
 
 /* ---------------------------------------------------------------------- */
 
+void BondFENE::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2)
+{
+  double r = sqrt(rsq);
+  double r0sq = r0[type] * r0[type];
+  double rlogarg = 1.0 - rsq / r0sq;
+
+  // Contribution from the attractive term
+  du = k[type] * r / rlogarg;
+  du2 = k[type] * (1.0 + rsq / r0sq) / (rlogarg * rlogarg);
+
+  // Contribution from the repulsive Lennard-Jones term
+  if (rsq < MY_CUBEROOT2 * sigma[type] * sigma[type]) {
+    double sr2 = sigma[type] * sigma[type] / rsq;
+    double sr6 = sr2 * sr2 * sr2;
+
+    du += 48.0 * epsilon[type] * sr6 * (0.5 - sr6) / r;
+    du2 += 48.0 * epsilon[type] * sr6 * (13.0 * sr6 - 3.5) / rsq;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void *BondFENE::extract(const char *str, int &dim)
 {
   dim = 1;