import numpy as np
_kT = 0.58622522 # kcal/mol
_k = np.logspace(-8,3,1000)
def _integrate( fn ):
t = np.linspace(0,np.pi,10000)
return np.trapz( fn(t), t[np.newaxis,:], axis = -1 )
_integral = _integrate( lambda t: np.cos(t[np.newaxis,:])*np.sin(t[np.newaxis,:])* np.exp((-0.5*_k[:,np.newaxis]*(t[np.newaxis,:])**2 )/_kT) ) / \
_integrate( lambda t: np.sin(t[np.newaxis,:]) * np.exp((-0.5*_k[:,np.newaxis]*(t[np.newaxis,:])**2 )/_kT) )
assert( (np.diff(_integral) <= 0).sum() == 0 )
[docs]
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)
[docs]
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]))