diff --git a/segmentmodel.py b/segmentmodel.py
index 6e9055071b4f1f008982e538283dca92e246fba3..24ba170024846276f29318385904fcbd8a06a323 100644
--- a/segmentmodel.py
+++ b/segmentmodel.py
@@ -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 = []