Commit 35e5c535 authored by cmaffeo2's avatar cmaffeo2
Browse files

Changed how crossover dihedrals are written when local_twist is False; fixed...

Changed how crossover dihedrals are written when local_twist is False; fixed routines that walk through helices
parent f4c495a8
......@@ -1476,6 +1476,35 @@ class SegmentModel(ArbdModel):
pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot)
## Get beads above and below
u1,u2 = [b.get_intrahelical_above() for b in (b1,b2)]
d1,d2 = [b.get_intrahelical_below() for b in (b1,b2)]
k_fn = lambda sep: (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
if u1 is not None and u2 is not None:
t0 = 0
k = k_fn( 0.5*(dists[b1][u1]+dists[b2][u2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( u1,b1,b2,u2, pot )
elif d1 is not None and d2 is not None:
t0 = 0
k = k_fn( 0.5*(dists[b1][d1]+dists[b2][d2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( d1,b1,b2,d2, pot )
elif d1 is not None and u2 is not None:
t0 = 180
k = k_fn( 0.5*(dists[b1][d1]+dists[b2][u2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( d1,b1,b2,u2, pot )
elif u1 is not None and d2 is not None:
t0 = 180
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 )
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__:
......@@ -1492,11 +1521,7 @@ class SegmentModel(ArbdModel):
pot = self.get_angle_potential(k,t0)
self.add_angle( b1,b2,o, pot )
## Get beads above and below
u1,u2 = [b.get_intrahelical_above() for b in (b1,b2)]
d1,d2 = [b.get_intrahelical_below() for b in (b1,b2)]
k_fn = lambda sep: (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
t0 = 90
if 'orientation_bead' in b1.__dict__:
o1 = b1.orientation_bead
......@@ -1519,37 +1544,17 @@ class SegmentModel(ArbdModel):
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( o2,b2,b1,d1, pot )
if 'orientation_bead' in b1.__dict__ and 'orientation_bead' in b2.__dict__:
if u1 is not None and u2 is not None:
t0 = 0
k = k_fn( 0.5*(dists[b1][u1]+dists[b2][u2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( u1,b1,b2,u2, pot )
elif d1 is not None and d2 is not None:
t0 = 0
k = k_fn( 0.5*(dists[b1][d1]+dists[b2][d2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( d1,b1,b2,d2, pot )
elif d1 is not None and u2 is not None:
t0 = 180
k = k_fn( 0.5*(dists[b1][d1]+dists[b2][u2]) )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( d1,b1,b2,u2, pot )
elif u1 is not None and d2 is not None:
t0 = 180
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)
lastParticle = None
for cl in crossovers:
c,A,B,d = cl
if A.particle is not lastParticle:
new_cl.append(cl)
lastParticle = A.particle
crossovers = new_cl
for i in range(len(crossovers)-1):
c1,A1,B1,dir1 = crossovers[i]
......@@ -1593,19 +1598,24 @@ class SegmentModel(ArbdModel):
def walk_through_helices(segment, direction=1, processed_segments=None):
"""
First and last segment should be same for circular helices
"""
assert( direction in (1,-1) )
if processed_segments = None:
if processed_segments == None:
processed_segments = set()
def segment_is_new_helix(s):
return isinstance(s0,DoubleStrandedSegment) and s0 not in processed_segments
return isinstance(s,DoubleStrandedSegment) and s not in processed_segments
s = s0 = segment
direction0 = direction
new_s = None
s = segment
## iterate intrahelically connected dsDNA segments
while segment_is_new_helix(s):
## comb through helicalwalk through helix
conn_locs = s.get_contour_sorted_connections_and_locations("intrahelical")[::direction]
processed_segments.add(new_s)
new_s = None
new_dir = None
for i in range(len(conn_locs)):
......@@ -1618,10 +1628,10 @@ class SegmentModel(ArbdModel):
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:
......@@ -1631,6 +1641,7 @@ class SegmentModel(ArbdModel):
def get_consecutive_crossovers(self):
## TODOTODO TEST
crossovers = []
processed_segments = set()
for s1 in self.segments:
......@@ -1643,12 +1654,14 @@ class SegmentModel(ArbdModel):
# 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] )
if s == s0 and len(tmp) > 0:
## end of circular helix, only add first crossover
tmp.append( s.get_contour_sorted_connections_and_locations("crossover")[::d][0] + [d] )
else:
tmp.extend( [cl + [d] for cl in 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
return crossovers
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