diff --git a/dnarbd/segmentmodel.py b/dnarbd/segmentmodel.py
index cfffd9dfeee8bd1dba8fa26d8b8f65baee5d03f0..784d792c24853dc4441efb7c67ac6ebfac56056d 100644
--- a/dnarbd/segmentmodel.py
+++ b/dnarbd/segmentmodel.py
@@ -1882,13 +1882,69 @@ class SegmentModel(ArbdModel):
pot = self.get_dihedral_potential(k,angle)
parent.add_dihedral(o1,b1,b2,o2, pot)
+ def k_angle(sep):
+ return 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
+ 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)]
+
+ k = k_xover_angle(sep=1)
+ pot = self.get_angle_potential(k,120)
+
+ 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 )
+
+ if 'orientation_bead' in b1.__dict__:
+ t0 = -90 if A.on_fwd_strand else 90
+ if not is_parallel: t0 *= -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 )
+
+
""" 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,18.5)
+ pot = self.get_bond_potential(4,12)
self.add_bond(b1,b2, pot)
+ dotProduct = b1.parent.contour_to_tangent(b1.contour_position).dot(
+ b2.parent.contour_to_tangent(b2.contour_position) )
+
+ ## Add potential to provide a particular orinetation
+ if local_twist:
+ add_local_crossover_orientation_potential(b1,b2, is_parallel=dotProduct > 0)
+
+
""" Add connection potentials """
for c,A,B in self.get_connections("sscrossover"):
b1,b2 = [loc.particle for loc in (c.A,c.B)]
@@ -1896,6 +1952,30 @@ class SegmentModel(ArbdModel):
pot = self.get_bond_potential(4,6)
self.add_bond(b1,b2, pot)
+ ## 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)
+ """
+
+ # 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 )
+
+ 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)
crossover_bead_pots = set()
@@ -1910,13 +1990,16 @@ 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)]
+ dotProduct = b1.parent.contour_to_tangent(b1.contour_position).dot(
+ b2.parent.contour_to_tangent(b2.contour_position) )
+ if dotProduct < 0:
+ tmp = d2
+ d2 = u2
+ u2 = tmp
-
- k_fn = lambda sep: (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
a = None
if u1 is not None and u2 is not None:
t0 = 0
@@ -1932,46 +2015,12 @@ class SegmentModel(ArbdModel):
a,b,c,d = (u1,b1,b2,d2)
if a is not None:
- k = k_fn( dists[b][a]+dists[c][d] )
+ 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 )
-
if local_twist:
- k = k_fn(1)
- pot = self.get_angle_potential(k,120)
- 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 )
-
-
- if 'orientation_bead' in b1.__dict__:
- t0 = 90
- if A.on_fwd_strand: t0 = -90
- o1 = b1.orientation_bead
- if u2 is not None:
- k = k_fn( dists[b2][u2] )
- pot = self.get_dihedral_potential(k,t0)
- self.add_dihedral( o1,b1,b2,u2, pot )
- elif d2 is not None:
- k = k_fn( 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: t0 = -90
- o2 = b2.orientation_bead
- if u1 is not None:
- k = k_fn( dists[b1][u1] )
- pot = self.get_dihedral_potential(k,t0)
- self.add_dihedral( o2,b2,b1,u1, pot )
- elif d1 is not None:
- k = k_fn( dists[b1][d1] )
- pot = self.get_dihedral_potential(k,-t0)
- self.add_dihedral( o2,b2,b1,d1, pot )
+ add_local_crossover_orientation_potential(b1,b2, is_parallel=dotProduct > 0)
## TODOTODO check that this works