Commit f4c495a8 authored by cmaffeo2's avatar cmaffeo2
Browse files

Updated crossover dihedrals when local_twist is False

parent 1eb40b0b
......@@ -33,6 +33,7 @@ TODO:
- rework Location class
- remove recursive calls
- document
- add unit test of helices connected to themselves
"""
class Location():
......@@ -1475,11 +1476,7 @@ class SegmentModel(ArbdModel):
pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot)
if not local_twist:
## TODOTODO
...
else:
if local_twist:
k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
if 'orientation_bead' in b1.__dict__:
# t0 = 90 + 60
......@@ -1543,15 +1540,115 @@ class SegmentModel(ArbdModel):
k = k_fn( 0.5*(dists[b1][u1]+dists[b2][d2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( u1,b1,b2,d2, pot )
## TODOTODO check that this works
for crossovers in self.get_consecutive_crossovers():
## filter crossovers
new_cl = []
particles = set()
for c,A,B in crossovers:
if A.particle not in particles:
new_cl.append(c)
particles.add(A.particle)
crossovers = new_cl
for i in range(len(crossovers)-1):
c1,A1,B1,dir1 = crossovers[i]
c2,A2,B2,dir2 = crossovers[j]
s1,s2 = [l.container for l in (A1,A2)]
sep = A1.get_nt_position(s1) - A2.get_nt_position(s2)
sep = np.abs(sep)
n1,n2,n3,n4 = (B1.particle, A1.particle, A2.particle, B2.particle)
## <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 = s1.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
t0 = sep*(360.0/s1.twist_per_nt)
# pdb.set_trace()
if c1.on_fwd_strand: t0 -= 120
if dir1 != dir2:
c2_on_fwd = not c2.on_fwd_strand
else:
c2_on_fwd = c2.on_fwd_strand
if c2_on_fwd: t0 += 120
# t0 = (t0 % 360
# if n2.idx == 0:
# print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
pot = self.get_dihedral_potential(k,t0)
dihedrals.add( n1,n2,n3,n4, pot )
def get_consecutive_crossovers():
def walk_through_helices(segment, direction=1, processed_segments=None):
assert( direction in (1,-1) )
if processed_segments = None:
processed_segments = set()
for s in self.segments:
if s in processed_segments: continue
for s in self.segments:
for cl in s.get_contour_sorted_connections_and_locations("crossover"):
...
def segment_is_new_helix(s):
return isinstance(s0,DoubleStrandedSegment) and s0 not in processed_segments
s = s0 = segment
direction0 = direction
while segment_is_new_helix(s):
## comb through helicalwalk through helix
conn_locs = s.get_contour_sorted_connections_and_locations("intrahelical")[::direction]
new_s = None
new_dir = None
for i in range(len(conn_locs)):
c,A,B = conn_locs[i]
## TODO: handle change of direction
address = 1*(direction==-1)
if A.address == address and segment_is_new_helix(B.container):
new_s = B.container
assert(B.address in (0,1))
new_dir = 2*(B.address == 0) - 1
break
processed_segments.add(s)
yield s,direction
s = new_s # will break if None
direction = new_dir
# if new_s is None:
# break
# else:
# s = new_s
# yield s
## return s
def get_consecutive_crossovers(self):
crossovers = []
processed_segments = set()
for s1 in self.segments:
if not isinstance(s1,DoubleStrandedSegment):
continue
if s1 in processed_segments: continue
s0,d0 = list(SegmentModel.walk_through_helices(s1,direction=-1))[-1]
# s,direction = get_start_of_helices()
tmp = []
for s,d in SegmentModel.walk_through_helices(s0,-d0):
tmp.extend( s.get_contour_sorted_connections_and_locations("crossover")[::d] )
processed_segments.add(s)
## TODOTODO
## if s and s0 are connected, do something
crossovers.append(tmp)
return crosovers
def _generate_strands(self):
self.strands = strands = []
......
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