Made intrahelical dsDNA bonds have geometric term subtracted

mrdna/model/nonbonded.py

@@ -116,7 +116,7 @@ class BasePotential():
         np.savetxt( self.filename(), np.array([r, u]).T, fmt="%f" )
 
 class HarmonicPotential(BasePotential):
-    def __init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix):
+    def __init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix, correct_geometry=False):
         self.k = k
         self.kscale_ = None
         BasePotential.__init__(self, r0, rRange, resolution, maxForce, max_potential, prefix)
@@ -143,10 +143,21 @@ class HarmonicPotential(BasePotential):
 #         self.kscale_ = 1.0
 
 class HarmonicBond(HarmonicPotential):
-    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):
+    def __init__(self, k, r0, rRange=(0,50), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/", correct_geometry=False, temperature=295):
         HarmonicPotential.__init__(self, k, r0, rRange, resolution, maxForce, max_potential, prefix)
-        self.type_ = "bond"
+        self.type_ = "gbond" if correct_geometry else "bond"
         self.kscale_ = 1.0
+        self.correct_geometry = correct_geometry
+        self.temperature = temperature
+
+    def potential(self,dr):
+        u = HarmonicPotential.potential(self,dr)
+        if self.correct_geometry:
+            du = 2*0.58622592*np.log(dr+self.r0) * self.temperature/295
+            du[np.logical_not(np.isfinite(du))] = 0
+            u = u-du
+        return u
+
 
 class HarmonicAngle(HarmonicPotential):
     def __init__(self, k, r0, rRange=(0,181), resolution=0.1, maxForce=None, max_potential=None, prefix="potentials/"):

mrdna/segmentmodel.py

@@ -1619,6 +1619,8 @@ class SegmentModel(ArbdModel):
             assert( kSpring >= 0 )
             if type_ == "bond":
                 self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,1200), max_potential=max_potential)
+            elif type_ == "gbond":
+                self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,1200), max_potential=max_potential, correct_geometry=True, temperature=self.temperature)
             elif type_ == "angle":
                 self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
                 # , resolution = 1, maxForce=0.1)
@@ -1627,9 +1629,10 @@ class SegmentModel(ArbdModel):
             else:
                 raise Exception("Unhandled potential type '%s'" % type_)
         return self._bonded_potential[key]
-    def get_bond_potential(self, kSpring, d):
+    def get_bond_potential(self, kSpring, d, correct_geometry=False):
         assert( d > 0.2 )
-        return self._get_potential("bond", kSpring, d)
+        return self._get_potential("gbond" if correct_geometry else "bond",
+                                   kSpring, d)
     def get_angle_potential(self, kSpring, d):
         return self._get_potential("angle", kSpring, d)
     def get_dihedral_potential(self, kSpring, d, max_potential=None):
@@ -2299,7 +2302,7 @@ class SegmentModel(ArbdModel):
             # dists[[b1,b2]] = dists[[b2,b1]] = sep
             if is_dsdna:
                 k = k_dsdna_bond(d)
-                bond = self.get_bond_potential(k,d)
+                bond = self.get_bond_potential(k,d,correct_geometry=True)
             else:
                 bond = self._get_wlc_sk_bond_potential(d)
             parent.add_bond( b1, b2, bond, exclude=True )