diff --git a/lib/reax/reax_connect.F b/lib/reax/reax_connect.F
index 886a68b8f3..958605c2e0 100644
--- a/lib/reax/reax_connect.F
+++ b/lib/reax/reax_connect.F
@@ -186,7 +186,7 @@ c$$$*     end if
 ********************************************************************** 
 ********************************************************************** 
 
-      subroutine srtbon1(lprune,lhb,hbcut_in)
+      subroutine srtbon1(lprune,lhb,hbcut_in,lhbnew_in)
 
 ********************************************************************** 
 #include "cbka.blk"
@@ -230,8 +230,9 @@ c$$$      call timer(65)
 c$$$      close (65)
 c$$$      end if
 
-c Transfer hbcut from C++ calling function
+c Transfer hbcut and lhbnew from C++ calling function
       hbcut = hbcut_in
+      lhbnew = lhbnew_in
 
       do i1=1,na
       abo(i1)=0.0d0
@@ -1528,6 +1529,7 @@ c
       ihb(nhb,8)=k3
 *     write (64,*)nhb,ihb(nhb,1),j3,j1,j2,nbohb,k1,k2,k3,bo(nbohb),
 *    $dishb
+
       end if
 
    10 continue
diff --git a/lib/reax/reax_poten.F b/lib/reax/reax_poten.F
index a407a0b8ca..cff08de890 100644
--- a/lib/reax/reax_poten.F
+++ b/lib/reax/reax_poten.F
@@ -2473,7 +2473,11 @@ c$$$      end if
       sin2=sinhu*sinhu
       exphu1=exp(-vhb1(ityhb)*boa)
       exphu2=exp(-vhb2(ityhb)*(rhu1+rhu2-2.0))
-      ehbh=(1.0-exphu1)*dehb(ityhb)*exphu2*sin2*sin2*sin2*sin2
+      if (lhbnew .eq. 0) then
+         ehbh=(1.0-exphu1)*dehb(ityhb)*exphu2*sin2*sin2*sin2*sin2
+      else
+         ehbh=(1.0-exphu1)*dehb(ityhb)*exphu2*sin2*sin2
+      endif
       ehb=ehb+ehbh
       estrain(j(2))=estrain(j(2))+ehbh !2nd atom energy
 
@@ -2482,12 +2486,20 @@ c$$$      end if
 *     Calculate first derivatives                                    *
 *                                                                    *
 ********************************************************************** 
-      dehbdbo=vhb1(ityhb)*exphu1*dehb(ityhb)*exphu2*sin2*sin2*
-     $sin2*sin2
-      dehbdv=(1.0-exphu1)*dehb(ityhb)*exphu2*
-     $4.0*sin2*sin2*sin2*sinhu*cos(hhb(i1)/2.0)
-      dehbdrda=(1.0-exphu1)*dehb(ityhb)*sin2*sin2*sin2*sin2*
-     $vhb2(ityhb)*(rhb(ityhb)/(rda*rda)-1.0/rhb(ityhb))*exphu2
+      if (lhbnew .eq. 0) then
+         dehbdbo=vhb1(ityhb)*exphu1*dehb(ityhb)*exphu2*sin2*sin2*
+     $        sin2*sin2
+         dehbdv=(1.0-exphu1)*dehb(ityhb)*exphu2*
+     $        4.0*sin2*sin2*sin2*sinhu*cos(hhb(i1)/2.0)
+         dehbdrda=(1.0-exphu1)*dehb(ityhb)*sin2*sin2*sin2*sin2*
+     $        vhb2(ityhb)*(rhb(ityhb)/(rda*rda)-1.0/rhb(ityhb))*exphu2
+      else
+         dehbdbo=vhb1(ityhb)*exphu1*dehb(ityhb)*exphu2*sin2*sin2
+         dehbdv=(1.0-exphu1)*dehb(ityhb)*exphu2*
+     $        2.0*sin2*sinhu*cos(hhb(i1)/2.0)
+         dehbdrda=(1.0-exphu1)*dehb(ityhb)*sin2*sin2*
+     $        vhb2(ityhb)*(rhb(ityhb)/(rda*rda)-1.0/rhb(ityhb))*exphu2
+      endif
 
       if (Lvirial.eq.1) then
          do k1=1,3