Commit f7cd72ff authored by cmaffeo2's avatar cmaffeo2
Browse files

Improved orientation potentials at junctions when local_twist is enabled; now...

Improved orientation potentials at junctions when local_twist is enabled; now dihedrals are more local and work in more situations
parent 6395dcac
......@@ -100,17 +100,21 @@ class Parent():
x.parent = None
def add_bond(self, i,j, bond, exclude=False):
assert( i is not j )
## TODO: how to handle duplicating and cloning bonds
# beads = [b for b in self]
# for b in (i,j): assert(b in beads)
self.bonds.append( (i,j, bond, exclude) )
def add_angle(self, i,j,k, angle):
assert( len(set((i,j,k))) == 3 )
# beads = [b for b in self]
# for b in (i,j,k): assert(b in beads)
self.angles.append( (i,j,k, angle) )
def add_dihedral(self, i,j,k,l, dihedral):
assert( len(set((i,j,k,l))) == 4 )
# beads = [b for b in self]
# for b in (i,j,k,l): assert(b in beads)
self.dihedrals.append( (i,j,k,l, dihedral) )
......
......@@ -306,7 +306,10 @@ def convert_maya_to_segments(helices):
botPos = 0.5 * (bot.get_position() + bot.basepair.get_position())
topPos = 0.5 * (top.get_position() + top.basepair.get_position())
segments[hn] = DoubleStrandedSegment(hn, num_bp = len(fwd_bases),
segname = hn
if segname[:6] == "helix_":
segname = 'D'+segname[6:]
segments[hn] = DoubleStrandedSegment(segname, num_bp = len(fwd_bases),
start_position = botPos,
end_position = topPos)
......
......@@ -213,19 +213,21 @@ class SegmentParticle(PointParticle):
self.other_neighbors = []
self.locations = []
def get_intrahelical_above(self):
def get_intrahelical_above(self, all_types=True):
""" Returns bead directly above self """
# assert( len(self.intrahelical_neighbors) <= 2 )
for b in self.intrahelical_neighbors:
if b.get_contour_position(self.parent, self.contour_position) > self.contour_position:
return b
if all_types or isinstance(b,type(self)):
return b
def get_intrahelical_below(self):
def get_intrahelical_below(self, all_types=True):
""" Returns bead directly below self """
# assert( len(self.intrahelical_neighbors) <= 2 )
for b in self.intrahelical_neighbors:
if b.get_contour_position(self.parent, self.contour_position) < self.contour_position:
return b
if all_types or isinstance(b,type(self)):
return b
def _neighbor_should_be_added(self,b):
if type(self.parent) != type(b.parent):
......@@ -2225,6 +2227,10 @@ class SegmentModel(ArbdModel):
## TODO: get better numbers our ssDNA model
elastic_modulus_times_area = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
d = 5*sep
if b1.type_.name[0] != b2.type_.name[0]:
""" Add a small extra distance to junction """
d += 3
k = conversion*elastic_modulus_times_area/d
# print(sep,d,k)
......@@ -2291,9 +2297,12 @@ class SegmentModel(ArbdModel):
parent.add_dihedral(b1,b2,b2.orientation_bead,b3, dihed)
else:
elif b1.type_.name[0] == "S" and b2.type_.name[0] == "S" and b3.type_.name[0] == "S":
## TODO: get correct number from ssDNA model
k = angle_spring_from_lp(sep,2)
else:
## Considered as a sscrossover below
continue
angle = self.get_angle_potential(k,180)
parent.add_angle( b1, b2, b3, angle )
......@@ -2314,7 +2323,24 @@ class SegmentModel(ArbdModel):
exclusions = set()
for b1 in beads:
exclusions.update( [(b1,b) for b in _recursively_get_beads_within(b1, 20, done=[b1])] )
""" In addition to bond exclusiosn, only add exclusions
within like-type segments (i.e. dsDNA or ssDNA, not
junctions between the two) """
t = type(b1.parent)
if t is DoubleStrandedSegment:
cutoff = 20
elif t is SingleStrandedSegment:
cutoff = 5
else:
raise ValueError("Unexpected polymer segment type")
for b in _recursively_get_beads_within(b1, cutoff, done=[b1]):
if isinstance(b.parent,t):
exclusions.add((b1,b))
else:
break
# exclusions.update( tmp )
if self.DEBUG: print("Adding %d exclusions" % len(exclusions))
for b1,b2 in exclusions:
......@@ -2363,112 +2389,81 @@ class SegmentModel(ArbdModel):
def k_xover_angle(sep):
return 0.5 * k_angle(sep)
def add_local_crossover_orientation_potential(b1,b2,is_parallel=True):
u1,u2 = [b.get_intrahelical_above() for b in (b1,b2)]
d1,d2 = [b.get_intrahelical_below() for b in (b1,b2)]
def add_local_crossover_strand_orientation_potential(b1,b2, b1_on_fwd_strand):
k = k_xover_angle(sep=1)
pot = self.get_angle_potential(k,120)
""" Adds a dihedral angle potential so bead b2 at opposite
end of crossover stays on correct side of helix of b1 """
if 'orientation_bead' in b1.__dict__:
o = b1.orientation_bead
self.add_angle( o,b1,b2, pot )
if 'orientation_bead' in b2.__dict__:
o = b2.orientation_bead
self.add_angle( b1,b2,o, pot )
u1 = b1.get_intrahelical_above(all_types=False)
d1 = b1.get_intrahelical_below(all_types=False)
if 'orientation_bead' in b1.__dict__:
t0 = -90 if A.on_fwd_strand else 90
if not is_parallel: t0 *= -1
sign = 1 if b1_on_fwd_strand else -1
o1 = b1.orientation_bead
if u2 is not None and isinstance(u2.parent,DoubleStrandedSegment):
k = k_xover_angle( dists[b2][u2] )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( o1,b1,b2,u2, pot )
elif d2 is not None and isinstance(d2.parent,DoubleStrandedSegment):
k = k_xover_angle( dists[b2][d2] )
pot = self.get_dihedral_potential(k,-t0)
self.add_dihedral( o1,b1,b2,d2, pot )
if 'orientation_bead' in b2.__dict__:
t0 = -90 if B.on_fwd_strand else 90
if not is_parallel: t0 *= -1
o2 = b2.orientation_bead
if u1 is not None and isinstance(u1.parent,DoubleStrandedSegment):
k = k_xover_angle( dists[b1][u1] )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( o2,b2,b1,u1, pot )
elif d1 is not None and isinstance(d1.parent,DoubleStrandedSegment):
k = k_xover_angle( dists[b1][d1] )
pot = self.get_dihedral_potential(k,-t0)
self.add_dihedral( o2,b2,b1,d1, pot )
# if b1.parent.name == "D72" or b2.parent.name == "D72":
# print()
# print(b1.parent.name, b2.parent.name, b1_on_fwd_strand)
# import pdb
# pdb.set_trace()
a,b,c = b2,b1,d1
if c is None or c is a:
c = u1
sign *= -1
if c is None or c is a: return
try:
d = b1.orientation_bead
except:
return
""" Add connection potentials """
for c,A,B in self.get_connections("terminal_crossover"):
## TODO: use a better description here
b1,b2 = [loc.particle for loc in (c.A,c.B)]
# pot = self.get_bond_potential(4,18.5)
pot = self.get_bond_potential(4,12)
self.add_bond(b1,b2, pot)
k = k_xover_angle(sep=1) # TODO
pot = self.get_dihedral_potential(k, sign*120)
self.add_dihedral( a,b,c,d, pot )
dotProduct = b1.parent.contour_to_tangent(b1.contour_position).dot(
b2.parent.contour_to_tangent(b2.contour_position) )
def add_local_tee_orientation_potential(b1,b2, b1_on_fwd_strand, b2_on_fwd_strand):
## Add potential to provide a particular orinetation
if local_twist:
add_local_crossover_orientation_potential(b1,b2, is_parallel=dotProduct > 0)
""" b1 is the end of a helix, b2 is in the middle This
adds a dihedral angle potential so helix of b1 is oriented
properly relative to strand on b2 """
u1,u2 = [b.get_intrahelical_above(all_types=False) for b in (b1,b2)]
d1,d2 = [b.get_intrahelical_below(all_types=False) for b in (b1,b2)]
""" Add connection potentials """
for c,A,B in self.get_connections("sscrossover"):
b1,b2 = [loc.particle for loc in (c.A,c.B)]
## TODO: improve parameters
pot = self.get_bond_potential(4,6)
self.add_bond(b1,b2, pot)
angle = 150
if not b2_on_fwd_strand: angle -= 180
## Add potentials to provide a sensible orientation
## TODO refine potentials against all-atom simulation data
"""
u1,u2 = [b.get_intrahelical_above() for b in (b1,b2)]
d1,d2 = [b.get_intrahelical_below() for b in (b1,b2)]
k_fn = k_angle
if isinstance(b1.parent,DoubleStrandedSegment):
a,b,c = (d1,b1,b2) if d1.parent is b1.parent else (u1,b1,b2)
else:
a,b,c = (d2,b2,b1) if d2.parent is b2.parent else (u2,b2,b1)
"""
a,b,c = u2,b2,b1
if a is None:
a = d2
angle -= 180
try:
d = b1.orientation_bead
except:
d = None
angle -= 120
# pot = self.get_angle_potential( k_fn(dists[a][b]), 120 ) # a little taller than 90 degrees
# self.add_angle( a,b,c, pot )
# o = b.orientation_bead
# angle=3
# pot = self.get_angle_potential( 1, 0 ) # a little taller than 90 degrees
# self.add_angle( a,b,c, pot )
while angle > 180:
angle -= 360
while angle < -180:
angle += 360
dotProduct = b1.parent.contour_to_tangent(b1.contour_position).dot(
b2.parent.contour_to_tangent(b2.contour_position) )
k = k_xover_angle(sep=1) # TODO
if a is not None and d is not None:
pot = self.get_dihedral_potential(k,angle)
self.add_dihedral( a,b,c,d, pot )
if local_twist:
add_local_crossover_orientation_potential(b1,b2, is_parallel=dotProduct > 0)
## Add 180 degree angle potential
a,b,c = b2,b1,u1
if c is None: c = d1
crossover_bead_pots = set()
for c,A,B in self.get_connections("crossover"):
b1,b2 = [loc.particle for loc in (A,B)]
## Avoid double-counting
if (b1,b2,A.on_fwd_strand,B.on_fwd_strand) in crossover_bead_pots: continue
crossover_bead_pots.add((b1,b2,A.on_fwd_strand,B.on_fwd_strand))
crossover_bead_pots.add((b2,b1,B.on_fwd_strand,A.on_fwd_strand))
pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot)
if c is not None:
pot = self.get_angle_potential(0.5*k,180)
self.add_angle( a,b,c, pot )
def add_parallel_crossover_potential(b1,b2):
## 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)]
u1,u2 = [b.get_intrahelical_above(all_types=False) for b in (b1,b2)]
d1,d2 = [b.get_intrahelical_below(all_types=False) for b in (b1,b2)]
dotProduct = b1.parent.contour_to_tangent(b1.contour_position).dot(
b2.parent.contour_to_tangent(b2.contour_position) )
if dotProduct < 0:
......@@ -2490,16 +2485,121 @@ class SegmentModel(ArbdModel):
t0 = 180
a,b,c,d = (u1,b1,b2,d2)
## TODO?: Check length-dependence of this potential
if a is not None:
k = k_xover_angle( dists[b][a]+dists[c][d] )
pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( a,b,c,d, pot )
...
""" Functions for adding crossover potentials """
def add_ss_crossover_potentials(connection,A,B, add_bond=True):
b1,b2 = [loc.particle for loc in (c.A,c.B)]
if (b1,b2,A.on_fwd_strand,B.on_fwd_strand) in processed_crossovers:
return
processed_crossovers.add((b1,b2,A.on_fwd_strand,B.on_fwd_strand))
processed_crossovers.add((b2,b1,B.on_fwd_strand,A.on_fwd_strand))
if b1 is b2:
""" Catch attempts to add "crossover potentials" at
intrahelical junctions between ds and ssDNA """
if A.container is not b1.parent:
b1 = A.container.get_nearest_bead(A.address)
if B.container is not b2.parent:
b2 = B.container.get_nearest_bead(B.address)
if b1 is b2:
return
## TODO: improve parameters
if add_bond:
pot = self.get_bond_potential(4,12)
self.add_bond(b1,b2, pot)
## Add potentials to provide a sensible orientation
## TODO refine potentials against all-atom simulation data
if local_twist:
add_local_crossover_strand_orientation_potential(b1,b2, A.on_fwd_strand)
add_local_crossover_strand_orientation_potential(b2,b1, B.on_fwd_strand)
def add_crossover_potentials(connection,A,B):
## TODO: use a better description here
b1,b2 = [loc.particle for loc in (c.A,c.B)]
if (b1,b2,A.on_fwd_strand,B.on_fwd_strand) in processed_crossovers:
return
processed_crossovers.add((b1,b2,A.on_fwd_strand,B.on_fwd_strand))
processed_crossovers.add((b2,b1,B.on_fwd_strand,A.on_fwd_strand))
if b1 is b2:
""" Catch attempts to add "crossover potentials" at
intrahelical junctions between ds and ssDNA """
return
""" Add bond potential """
pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot)
""" Add parallel helices potential, possibly """
## Add potential to provide a particular orinetation
nt1,nt2 = [l.get_nt_pos() for l in (c.A,c.B)]
is_end1, is_end2 = [nt in (0,l.container.num_nt-1) for nt,l in zip((nt1,nt2),(c.A,c.B))]
is_T_junction = (is_end1 and not is_end2) or (is_end2 and not is_end1)
if (not is_end1) and (not is_end2):
## TODO?: Only apply this potential if not local_twist
add_parallel_crossover_potential(b1,b2)
# dotProduct = b1.parent.contour_to_tangent(b1.contour_position).dot(
# b2.parent.contour_to_tangent(b2.contour_position) )
if local_twist:
add_local_crossover_orientation_potential(b1,b2, is_parallel=dotProduct > 0)
if is_T_junction:
""" Special case: one helix extends away from another in T-shaped junction """
if is_end1:
b1_forward = c.A.on_fwd_strand if nt1 == 0 else not c.A.on_fwd_strand
add_local_tee_orientation_potential(b1,b2, b1_forward, c.B.on_fwd_strand)
else:
add_local_crossover_strand_orientation_potential(b1,b2, c.A.on_fwd_strand)
if is_end2:
b2_forward = c.B.on_fwd_strand if nt2 == 0 else not c.B.on_fwd_strand
add_local_tee_orientation_potential(b2,b1, b2_forward, c.A.on_fwd_strand)
else:
add_local_crossover_strand_orientation_potential(b2,b1, c.B.on_fwd_strand)
else:
""" Normal case: add orientation potential """
add_local_crossover_strand_orientation_potential(b1,b2, c.A.on_fwd_strand)
add_local_crossover_strand_orientation_potential(b2,b1, c.B.on_fwd_strand)
## TODOTODO check that this works
""" Add connection potentials """
processed_crossovers = set()
# pdb.set_trace()
for c,A,B in self.get_connections("sscrossover"):
p1,p2 = [loc.container for loc in (c.A,c.B)]
assert(any([isinstance(p,SingleStrandedSegment) for p in (p1,p2)]))
add_ss_crossover_potentials(c,A,B)
for c,A,B in self.get_connections("intrahelical"):
ps = [loc.container for loc in (c.A,c.B)]
if any([isinstance(p,SingleStrandedSegment) for p in ps]) and \
any([isinstance(p,DoubleStrandedSegment) for p in ps]):
add_ss_crossover_potentials(c,A,B, add_bond=False)
for c,A,B in sum([self.get_connections(term) for term in ("crossover","terminal_crossover")],[]):
p1,p2 = [loc.container for loc in (c.A,c.B)]
if any([isinstance(p,SingleStrandedSegment) for p in (p1,p2)]):
add_ss_crossover_potentials(c,A,B)
else:
add_crossover_potentials(c,A,B)
## todotodo check that this works
for crossovers in self.get_consecutive_crossovers():
if local_twist: break
## filter crossovers
......@@ -2542,7 +2642,8 @@ class SegmentModel(ArbdModel):
# if n2.idx == 0:
# print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
if sep == 0 and n1 is not n4:
# if sep == 0 and n1 is not n4:
if sep == 0:
# pot = self.get_angle_potential(k,t0)
# self.add_angle( n1,n2,n4, pot )
pass
......
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