From c84c708dd9b341a085d4841ef33fd76fcf221f1a Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Sun, 2 Dec 2018 10:16:11 -0600
Subject: [PATCH] Updated persistence length to be a little bit shorter

---
 mrdna/segmentmodel.py | 22 ++++++++++++++++++++--
 mrdna/version.py      |  2 +-
 2 files changed, 21 insertions(+), 3 deletions(-)

diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 40ab067..be162ec 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -2072,6 +2072,22 @@ class SegmentModel(ArbdModel):
         #     print("Helix {}: bps {}, beads {}, separation {}".format(s.name, s.num_nt, beadsum, sepsum))
 
         """ Add intrahelical angle potentials """
+        def get_effective_dsDNA_Lp(sep):
+
+            """ The persistence length of our model was found to be a
+            little off (probably due to NB interactions). This
+            attempts to compensate """
+
+            ## For 1 bp, Lp=559, for 25 Lp = 524
+            beads_per_bp = sep/2
+            Lp0 = 147
+            # return 0.93457944*Lp0 ;# factor1
+            return 0.97*Lp0 ;# factor2
+            # factor = bead_per_bp * (0.954-0.8944
+            # return Lp0 * bead_per_bp
+
+        empirical_compensation_factor = max_basepairs_per_bead
+        
         if self.DEBUG: print("Adding intrahelical angle potentials")
         for b1,b2,b3 in self._get_intrahelical_angle_beads():
             ## TODO: could be slightly smarter about sep
@@ -2080,7 +2096,8 @@ class SegmentModel(ArbdModel):
 
             kT = 0.58622522         # kcal/mol
             if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":
-                k = angle_spring_from_lp(sep,147)
+                Lp = get_effective_dsDNA_Lp(sep)
+                k = angle_spring_from_lp(sep,Lp)
                 if local_twist:
                     k_dihed = 0.25*k
                     k *= 0.75    # reduce because orientation beads impose similar springs
@@ -2134,7 +2151,8 @@ class SegmentModel(ArbdModel):
                     parent = self._getParent( b1, b2 )
 
                     """ Add heuristic 90 degree potential to keep orientation bead orthogonal """
-                    k = 0.5*angle_spring_from_lp(sep,147)
+                    Lp = get_effective_dsDNA_Lp(sep)
+                    k = 0.5*angle_spring_from_lp(sep,Lp)
                     pot = self.get_angle_potential(k,90)
                     parent.add_angle(o1,b1,b2, pot)
                     parent.add_angle(b1,b2,o2, pot)
diff --git a/mrdna/version.py b/mrdna/version.py
index 80ef677..b12c139 100644
--- a/mrdna/version.py
+++ b/mrdna/version.py
@@ -133,7 +133,7 @@ def get_version(abbrev=7):
 
 
 if __name__ == "__main__":
-    print( get_git_version() )
+    print( get_version() )
 
 def read_version_file(filename):
     with open(filename) as fh:
-- 
GitLab