From 81ab7c6dacdac19b9c8029c02505d2d613d625b8 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Fri, 20 Jul 2018 15:10:50 -0500
Subject: [PATCH] Fixed spring constant for dihedral that spans crossovers,
keeping helices parallel (note this differs from prototyped web service where
90 degree angle potentials were used)
Changed dna bead type clustering to have higher resolution
Streamlined some small parts of the code
---
segmentmodel.py | 56 ++++++++++++++++++++++++++++++-------------------
1 file changed, 35 insertions(+), 21 deletions(-)
diff --git a/segmentmodel.py b/segmentmodel.py
index b00a29f..1b7300e 100644
--- a/segmentmodel.py
+++ b/segmentmodel.py
@@ -1654,10 +1654,9 @@ class SegmentModel(ArbdModel):
data = np.array([b.num_nt for b in beads])[:,np.newaxis]
order = int(2-np.log10(2*max_basepairs_per_bead)//1)
try:
- clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/80, criterion="distance")
+ clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/500, criterion="distance")
cluster_size = [np.mean(data[clusters == i]) for i in np.unique(clusters)]
except:
- ## TODO debug
clusters = data
cluster_size = np.arange(len(data))+1
for b,c in zip(beads,clusters):
@@ -1667,15 +1666,30 @@ class SegmentModel(ArbdModel):
data = np.array([b.num_nt for b in beads])[:,np.newaxis]
order = int(2-np.log10(max_nucleotides_per_bead)//1)
try:
- clusters = hcluster.fclusterdata(data, float(max_basepairs_per_bead)/40, criterion="distance")
+ clusters = hcluster.fclusterdata(data, float(max_nucleotides_per_bead)/500, criterion="distance")
cluster_size = [np.mean(data[clusters == i]) for i in np.unique(clusters)]
except:
- ## TODO debug
clusters = data
cluster_size = np.arange(len(data))+1
for b,c in zip(beads,clusters):
_assign_bead_type(b, cluster_size[c-1], decimals=order)
+
+ # for bead in [b for s in segments for b in s]:
+ # num_nt0 = bead.num_nt
+ # # bead.num_nt = np.around( np.float32(num_nt), decimals=decimals )
+ # key = (bead.type_.name[0].upper(), bead.num_nt)
+ # if key in beadtype_s:
+ # bead.type_ = beadtype_s[key]
+ # else:
+ # t = deepcopy(bead.type_)
+ # t.__dict__["nts"] = bead.num_nt*2 if t.name[0].upper() in ("D","O") else bead.num_nt
+ # # t.name = t.name + "%03d" % (t.nts*10**decimals)
+ # t.name = t.name + "%.16f" % (t.nts)
+ # print( "{} --> {} ({})".format(num_nt0, bead.num_nt, t.name) )
+ # beadtype_s[key] = bead.type_ = t
+
+
""" Update bead indices """
self._countParticleTypes() # probably not needed here
self._updateParticleOrder()
@@ -1853,30 +1867,28 @@ class SegmentModel(ArbdModel):
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
- 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 )
+ a,b,c,d = (u1,b1,b2,u2)
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 )
+ a,b,c,d = (d1,b1,b2,d2 )
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 )
+ a,b,c,d = (d1,b1,b2,u2)
elif u1 is not None and d2 is not None:
t0 = 180
- k = k_fn( 0.5*(dists[b1][u1]+dists[b2][d2]) )
+ a,b,c,d = (u1,b1,b2,d2)
+
+ if a is not None:
+ k = k_fn( dists[b][a]+dists[c][d] )
pot = self.get_dihedral_potential(k,t0)
- self.add_dihedral( u1,b1,b2,d2, pot )
+ self.add_dihedral( a,b,c,d, 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
+ k = k_fn(1)
pot = self.get_angle_potential(k,120)
if 'orientation_bead' in b1.__dict__:
o = b1.orientation_bead
@@ -1933,11 +1945,13 @@ class SegmentModel(ArbdModel):
n1,n2,n3,n4 = (B1.particle, A1.particle, A2.particle, B2.particle)
- ## <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
- ## Assume A is small
- ## int[B_] := Normal[Integrate[ Series[Cos[x] Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]/
- ## Integrate[Series[Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]]
- ## Actually, without assumptions I get fitFun below
+ """
+ <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
+ Assume A is small
+ int[B_] := Normal[Integrate[ Series[Cos[x] Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]/
+ Integrate[Series[Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]]
+ Actually, without assumptions I get fitFun below
+ """
## From http://www.annualreviews.org/doi/pdf/10.1146/annurev.bb.17.060188.001405
## units "3e-19 erg cm/ 295 k K" "nm" =~ 73
--
GitLab