diff --git a/mrdna/model/arbdmodel.py b/mrdna/model/arbdmodel.py
index e2a1f90e477eff0b6b5d613bea18ee2974c8029e..88d22b581f515e61f8c97f226e166834e6366708 100644
--- a/mrdna/model/arbdmodel.py
+++ b/mrdna/model/arbdmodel.py
@@ -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) )
diff --git a/mrdna/readers/polygon_mesh.py b/mrdna/readers/polygon_mesh.py
index 4e2027c63ac5cedb09bb5a70ffcaf638f92151e9..700ad09fd34295430323d766254195bd4e6b4519 100644
--- a/mrdna/readers/polygon_mesh.py
+++ b/mrdna/readers/polygon_mesh.py
@@ -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)
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index afbc29588a4693b87893bee064eb84ebf2351c6d..62502502fed8c7bf3fb23798ce5e3bba28ed2a14 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -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