Skip to content
Snippets Groups Projects
Commit 7212bdce authored by cmaffeo2's avatar cmaffeo2
Browse files

Updated bead model potentials to work with antiparallel DNA; also made... 

Updated bead model potentials to work with antiparallel DNA; also made terminal_crossovers bond shorter at 1.2 nm (vs 1.85)
parent cec3b33e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 89 additions and 40 deletions
dnarbd/segmentmodel.py
+ 89
− 40
View file @ 7212bdce
...@@ -1882,13 +1882,69 @@ class SegmentModel(ArbdModel): ...@@ -1882,13 +1882,69 @@ class SegmentModel(ArbdModel):
pot = self.get_dihedral_potential(k,angle) pot = self.get_dihedral_potential(k,angle)
parent.add_dihedral(o1,b1,b2,o2, pot) 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 """ """ Add connection potentials """
for c,A,B in self.get_connections("terminal_crossover"): for c,A,B in self.get_connections("terminal_crossover"):
## TODO: use a better description here ## TODO: use a better description here
b1,b2 = [loc.particle for loc in (c.A,c.B)] 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) 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 """ """ Add connection potentials """
for c,A,B in self.get_connections("sscrossover"): for c,A,B in self.get_connections("sscrossover"):
b1,b2 = [loc.particle for loc in (c.A,c.B)] b1,b2 = [loc.particle for loc in (c.A,c.B)]
...@@ -1896,6 +1952,30 @@ class SegmentModel(ArbdModel): ...@@ -1896,6 +1952,30 @@ class SegmentModel(ArbdModel):
pot = self.get_bond_potential(4,6) pot = self.get_bond_potential(4,6)
self.add_bond(b1,b2, pot) 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() crossover_bead_pots = set()
...@@ -1910,13 +1990,16 @@ class SegmentModel(ArbdModel): ...@@ -1910,13 +1990,16 @@ class SegmentModel(ArbdModel):
pot = self.get_bond_potential(4,18.5) pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot) self.add_bond(b1,b2, pot)
## Get beads above and below ## Get beads above and below
u1,u2 = [b.get_intrahelical_above() for b in (b1,b2)] u1,u2 = [b.get_intrahelical_above() for b in (b1,b2)]
d1,d2 = [b.get_intrahelical_below() 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 a = None
if u1 is not None and u2 is not None: if u1 is not None and u2 is not None:
t0 = 0 t0 = 0
...@@ -1932,46 +2015,12 @@ class SegmentModel(ArbdModel): ...@@ -1932,46 +2015,12 @@ class SegmentModel(ArbdModel):
a,b,c,d = (u1,b1,b2,d2) a,b,c,d = (u1,b1,b2,d2)
if a is not None: 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) pot = self.get_dihedral_potential(k,t0)
self.add_dihedral( a,b,c,d, pot ) self.add_dihedral( a,b,c,d, pot )
if local_twist: if local_twist:
k = k_fn(1) add_local_crossover_orientation_potential(b1,b2, is_parallel=dotProduct > 0)
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 )
## TODOTODO check that this works ## TODOTODO check that this works
......
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