Updated method for spring constant

mrdna/segmentmodel.py

@@ -15,6 +15,8 @@ from scipy import interpolate
 from .model.CanonicalNucleotideAtoms import canonicalNtFwd, canonicalNtRev, seqComplement
 from .model.CanonicalNucleotideAtoms import enmTemplateHC, enmTemplateSQ, enmCorrectionsHC
 
+from .model.spring_from_lp import k_angle as angle_spring_from_lp
+
 # import pdb
 """
 TODO:
@@ -2069,10 +2071,7 @@ 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":
-                ## <cos(q)> = exp(-s/Lp) = integrate( x^4 exp(-A x^2) / 2, {x, 0, pi} ) / integrate( x^2 exp(-A x^2), {x, 0, pi} )
-                ## <cos(q)> ~ 1 - 3/4A                                                                                            
-                ## where A = k_spring / (2 kT)                                                                                    
-                k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
+                k = angle_spring_from_lp(sep,147)
                 if local_twist:
                     ## TODO optimize this paramter
                     k *= 0.5    # halve because orientation beads have similar springs
@@ -2082,7 +2081,7 @@ class SegmentModel(ArbdModel):
 
             else:
                 ## TODO: get correct number from ssDNA model
-                k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2
+                k = angle_spring_from_lp(sep,2)
 
             angle = self.get_angle_potential(k,180)
             parent.add_angle( b1, b2, b3, angle )
@@ -2126,7 +2125,7 @@ class SegmentModel(ArbdModel):
                     parent = self._getParent( b1, b2 )
 
                     """ Add heuristic 90 degree potential to keep orientation bead orthogonal """
-                    k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
+                    k = 0.5*angle_spring_from_lp(sep,147)
                     pot = self.get_angle_potential(k,90)
                     parent.add_angle(o1,b1,b2, pot)
                     parent.add_angle(b1,b2,o2, pot)
@@ -2146,7 +2145,8 @@ class SegmentModel(ArbdModel):
                     parent.add_dihedral(o1,b1,b2,o2, pot)
 
         def k_angle(sep):
-            return 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
+            return  0.5*angle_spring_from_lp(sep,147)
+
         def k_xover_angle(sep):
             return 0.5 * k_angle(sep)