Commit c16fbbed authored by cmaffeo2's avatar cmaffeo2
Browse files

Made beadModelTwist omit beadModel-style crossover bonds and angles

parent 2366e578
......@@ -257,7 +257,7 @@ class beadModelTwist():
# self._dihedralParamFiles = set()
self.twistPersistenceLength = twistPersistenceLength
self.apply_extra_crossover_potentials = False
self._buildModel(part, maxBpsPerDNode, maxNtsPerSNode)
......@@ -1255,12 +1255,13 @@ component "data" value 3
self.addAngle( *args )
self.addAngle( n1,n2,n3,Angle(k,180) )
a,d = self._getCrossoverAnglesAndDihedrals()
for nodes,sep in a:
n1,n2,n3 = nodes
k = (1.0/2) * 1.5 * kT * (1.0 / (1-exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
self.addAngle( n1,n2,n3,Angle(k,90) )
if self.apply_extra_crossover_potentials:
a,d = self._getCrossoverAnglesAndDihedrals()
for nodes,sep in a:
n1,n2,n3 = nodes
k = (1.0/2) * 1.5 * kT * (1.0 / (1-exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
self.addAngle( n1,n2,n3,Angle(k,90) )
## Intrahelical 90 deg orientation angles
for nodes,seps in self._getOrientationAngles():
......@@ -1289,40 +1290,42 @@ component "data" value 3
def _buildDihedrals(self, prefix, directory="potentials"):
kT = 0.58622522 # kcal/mol
a,d = self._getCrossoverAnglesAndDihedrals()
for nodes,sep,isFwd1,isFwd2 in d:
n1,n2,n3,n4 = nodes
## <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 = self.twistPersistenceLength/0.34 # set semi-arbitrarily as there is a large spread in literature
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) - exp(-sep/Lp)
k = opt.leastsq( fitFun, x0=exp(-sep/Lp) )
k = k[0][0] * 2*kT*0.00030461742
# intrinsicDegrees=30
# fitFun = lambda x: (1.0/(2*x) - 2*np.sqrt(np.pi)*np.exp(-4*np.pi**2*x) / (np.sqrt(x)*erf(2*np.pi*np.sqrt(x))) ) - \
# ( (intrinsicDegrees*np.pi/180)**2 + 2*(1-exp(-sep/Lp)) )
# k = opt.leastsq( fitFun, x0=1/(1-exp(-sep/Lp)) )
# k = k[0][0] * 2*kT*0.00030461742
t0 = sep*(360.0/10.5)
# pdb.set_trace()
if isFwd1[0]: t0 -= 120
if isFwd2[0]: t0 += 120
t0 = t0 % 360
# if n2.idx == 0:
# print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
self.addDihedral( n1,n2,n3,n4,Dihedral(k,t0) )
if self.apply_extra_crossover_potentials:
a,d = self._getCrossoverAnglesAndDihedrals()
for nodes,sep,isFwd1,isFwd2 in d:
n1,n2,n3,n4 = nodes
## <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 = self.twistPersistenceLength/0.34 # set semi-arbitrarily as there is a large spread in literature
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) - exp(-sep/Lp)
k = opt.leastsq( fitFun, x0=exp(-sep/Lp) )
k = k[0][0] * 2*kT*0.00030461742
# intrinsicDegrees=30
# fitFun = lambda x: (1.0/(2*x) - 2*np.sqrt(np.pi)*np.exp(-4*np.pi**2*x) / (np.sqrt(x)*erf(2*np.pi*np.sqrt(x))) ) - \
# ( (intrinsicDegrees*np.pi/180)**2 + 2*(1-exp(-sep/Lp)) )
# k = opt.leastsq( fitFun, x0=1/(1-exp(-sep/Lp)) )
# k = k[0][0] * 2*kT*0.00030461742
t0 = sep*(360.0/10.5)
# pdb.set_trace()
if isFwd1[0]: t0 -= 120
if isFwd2[0]: t0 += 120
t0 = t0 % 360
# if n2.idx == 0:
# print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
self.addDihedral( n1,n2,n3,n4,Dihedral(k,t0) )
for nodes,seps in self._getOrientationDihedrals():
n1,n2,n3,n4 = nodes
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment