Speed up model construction using lookup table for twist spring constants

mrdna/model/spring_from_lp.py

@@ -14,3 +14,44 @@ assert( (np.diff(_integral) <= 0).sum() == 0 )
 def k_angle(sep,Lp):
     val = np.exp(-sep/Lp)
     return np.interp(val,_integral,_k) * 0.00030461742 # convert to degree^2
+
+""" Twist! """
+from scipy.special import erf
+_k_twist = np.logspace(-8,3,10000) # in units of kT
+with np.errstate(divide='ignore',invalid='ignore'):
+    _integral_twist = np.real(erf( (4*np.pi*_k_twist + 1j)/(2*np.sqrt(_k_twist)) )) * np.exp(-1/(4*_k_twist)) / erf(2*np.sqrt(_k_twist)*np.pi)
+def k_twist(sep,Lp,temperature=295):
+    kT = temperature * 0.0019872065 # kcal/mol
+    val = np.exp(-sep/Lp)
+    return np.interp(val,_integral_twist,_k_twist) * 2 * kT * 0.00030461742 # convert to degree^2
+
+def _TEST_get_spring_constant():
+
+    def __get_twist_spring_constant(sep, temperature=295):
+        import scipy.optimize as opt
+
+        """ sep in nt """
+        kT = temperature * 0.0019872065 # kcal/mol
+        twist_persistence_length = 90  # set semi-arbitrarily as there is a large spread in literature
+        ## <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
+        ##   Assume A is small
+        ## int[B_] :=  Normal[Integrate[ Series[Cos[x] Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]/
+        ##             Integrate[Series[Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]]
+
+        ## Actually, without assumptions I get fitFun below
+        ## From http://www.annualreviews.org/doi/pdf/10.1146/annurev.bb.17.060188.001405
+        ##   units "3e-19 erg cm/ 295 k K" "nm" =~ 73
+        Lp = twist_persistence_length/0.34
+
+        fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - np.exp(-sep/Lp)
+        k = opt.leastsq( fitFun, x0=np.exp(-sep/Lp) )
+        return k[0][0] * 2*kT*0.00030461742
+
+    seps = np.linspace(0.5,25,1000)
+    vals1 = np.array([__get_twist_spring_constant(s) for s in seps])
+    vals2 = k_twist(seps, 90/0.34 )
+
+    abs_diff = np.abs(vals1-vals2)
+    max_diff = np.max(abs_diff)
+    idx = abs_diff==max_diff
+    print("Max diff ({}) at {}, where sep = {}, vals1 = {}, vals2 = {}".format(max_diff, np.where(idx)[0], seps[idx], vals1[idx], vals2[idx]))

mrdna/segmentmodel.py

@@ -16,6 +16,7 @@ from .model.CanonicalNucleotideAtoms import canonicalNtFwd, canonicalNtRev, seqC
 from .model.CanonicalNucleotideAtoms import enmTemplateHC, enmTemplateSQ, enmCorrectionsHC
 
 from .model.spring_from_lp import k_angle as angle_spring_from_lp
+from .model.spring_from_lp import k_twist as twist_spring_from_lp
 
 import csv
 
@@ -1697,21 +1698,11 @@ class SegmentModel(ArbdModel):
 
     def _get_twist_spring_constant(self, sep):
         """ sep in nt """
-        kT = self.temperature * 0.0019872065 # kcal/mol
         twist_persistence_length = 90  # set semi-arbitrarily as there is a large spread in literature
-        ## <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
-        ##   Assume A is small
-        ## int[B_] :=  Normal[Integrate[ Series[Cos[x] Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]/
-        ##             Integrate[Series[Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]]
-
-        ## Actually, without assumptions I get fitFun below
-        ## From http://www.annualreviews.org/doi/pdf/10.1146/annurev.bb.17.060188.001405
         ##   units "3e-19 erg cm/ 295 k K" "nm" =~ 73
         Lp = twist_persistence_length/0.34 
 
-        fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - np.exp(-sep/Lp)
-        k = opt.leastsq( fitFun, x0=np.exp(-sep/Lp) )
-        return k[0][0] * 2*kT*0.00030461742
+        return twist_spring_from_lp(sep, Lp, self.temperature)
 
     def extend(self, other, copy=True, include_strands=False):
         assert( isinstance(other, SegmentModel) )