From 671f16053f14a6cb1e75ed64a3f618385f039b68 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Mon, 17 Sep 2018 15:45:44 -0500
Subject: [PATCH] Fixed multiple issues; 3-prime ends of nucleotides now
terminate correctly; strands not cleared unless atoms exist; strands
regenerated with atomic model; cadnano models initialized with more correct
azimuthal angle
---
mrdna/model/arbdmodel.py | 13 ++++++-------
mrdna/model/nbPot.py | 22 +++++++++++++++-------
mrdna/readers/cadnano_segments.py | 5 +++--
mrdna/segmentmodel.py | 21 +++++++++++++++------
4 files changed, 39 insertions(+), 22 deletions(-)
diff --git a/mrdna/model/arbdmodel.py b/mrdna/model/arbdmodel.py
index cd21983..6fb1c03 100644
--- a/mrdna/model/arbdmodel.py
+++ b/mrdna/model/arbdmodel.py
@@ -1105,13 +1105,12 @@ cellBasisVector3 0 0 {cellZ}
if {{$nLast == 0}} {{
temperature 300
- fixedAtoms on
- fixedAtomsForces on
- fixedAtomsFile $prefix.pdb
- fixedAtomsCol B
- minimize 2400
-
- fixedAtoms off
+ # fixedAtoms on
+ # fixedAtomsForces on
+ # fixedAtomsFile $prefix.pdb
+ # fixedAtomsCol B
+ # minimize 2400
+ # fixedAtoms off
minimize 2400
}} else {{
bincoordinates {output_directory}/$prefix-$nLast.restart.coor
diff --git a/mrdna/model/nbPot.py b/mrdna/model/nbPot.py
index 1b9c145..56e46b9 100644
--- a/mrdna/model/nbPot.py
+++ b/mrdna/model/nbPot.py
@@ -3,8 +3,8 @@ import scipy.optimize as opt
from scipy.interpolate import interp1d
from scipy.signal import savgol_filter as savgol
-from .nonbonded import NonbondedScheme
-from .. import get_resource_path
+from mrdna.model.nonbonded import NonbondedScheme
+from mrdna import get_resource_path
dxParam = -1
@@ -33,7 +33,7 @@ y0 = d[:,0] # kcal/mol for 10bp with 10bp, according to JY
y0 = y0 - np.mean(y0[x0>(38.5 + dxParam)]) # set 0 to tail of potential
y0 = savgol(y0, int(5.0/(x0[1]-x0[0])), 3) # smooth potential
-def pot(x,*parms):
+def parametric_potential(x,*parms):
x = np.array(x)
u = interp1d(x0, np.array(parms),
bounds_error = False, fill_value="extrapolate",
@@ -41,7 +41,7 @@ def pot(x,*parms):
u[x>np.max(x0)] = 0
return u
-def potForBp(potfn,x,numBp,*parms):
+def integrated_potential(potfn,x,numBp,*parms):
u = np.zeros(np.shape(x))
for h in range(-numBp,numBp+1):
r = np.sqrt(x**2 + (3.4*h)**2)
@@ -50,7 +50,8 @@ def potForBp(potfn,x,numBp,*parms):
return u
def fitFun(x,*parms):
- return 10 * potForBp(pot,x,20,*parms)
+ ## 1-turn iteracting with 4 turns (more is not needed)
+ return 10 * integrated_potential(parametric_potential,x,20,*parms)
xinterp = np.linspace(np.min(x0),np.max(x0),1000) # interpolate potentialb
yinterp = interp1d(x0,y0)(xinterp)
@@ -62,9 +63,9 @@ def nbPot(x,bps1,bps2):
if larger >= 2:
largerInt = int(0.5*larger)
smaller *= larger/(2*largerInt)
- y = smaller * potForBp(pot,x,largerInt,*popt)
+ y = smaller * integrated_potential(parametric_potential,x,largerInt,*popt)
else:
- y = smaller * larger * pot(x,*popt)
+ y = smaller * larger * parametric_potential(x,*popt)
y = maxForce(x,y, 100)
y = savgol(y, int(5.0/(x0[1]-x0[0])), 3)
return y
@@ -99,3 +100,10 @@ class nbDnaScheme(NonbondedScheme):
u = nbPot(r, 0.5*typeA.nts, 0.5*typeB.nts)
return u
nbDnaScheme = nbDnaScheme()
+
+if __name__ == "__main__":
+ ## run some tests
+ np.savetxt("jj-filtered.dat", np.array((x0,y0)).T)
+ y = fitFun(x,*popt)
+ np.savetxt("1-bp-fit.dat", np.array((x,y)).T),
+
diff --git a/mrdna/readers/cadnano_segments.py b/mrdna/readers/cadnano_segments.py
index 849de60..59fb2fe 100644
--- a/mrdna/readers/cadnano_segments.py
+++ b/mrdna/readers/cadnano_segments.py
@@ -5,7 +5,7 @@ import os,sys
from glob import glob
import re
-from ..coords import readArbdCoords, readAvgArbdCoords
+from ..coords import readArbdCoords, readAvgArbdCoords, rotationAboutAxis
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
from ..model.dna_sequence import m13 as m13seq
@@ -311,7 +311,8 @@ class cadnano_part(SegmentModel):
## TODO get sequence from cadnano api
if zMid in strandOccupancies[0] and zMid in strandOccupancies[1]:
kwargs['num_bp'] = numBps
- seg = DoubleStrandedSegment(**kwargs,**posargs1)
+ start_orientation = rotationAboutAxis(np.array((0,0,1)), (zid1-3.25)*34)
+ seg = DoubleStrandedSegment(**kwargs,**posargs1, start_orientation = start_orientation)
elif zMid in strandOccupancies[0]:
kwargs['num_nt'] = numBps
seg = SingleStrandedSegment(**kwargs,**posargs1)
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index f260055..c8efc6a 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -506,6 +506,9 @@ class Segment(ConnectableElement, Group):
else:
orientation0 = np.eye(3) if axis.dot(zAxis) > 0 else \
rotationAboutAxis( np.array((1,0,0)), 180, normalizeAxis=False )
+ if self.start_orientation is not None:
+ orientation0 = orientation0.dot(self.start_orientation)
+
orientation = rotationAboutAxis( axis, self.twist_per_nt*self.contour_to_nt_pos(s), normalizeAxis=False )
orientation = orientation.dot(orientation0)
else:
@@ -572,7 +575,7 @@ class Segment(ConnectableElement, Group):
if scale is not None and scale != 1:
for a in atoms:
a.position = scale*a.position
- a.beta = 0
+ a.fixed = 1
atoms.position = pos - atoms.atoms_by_name["C1'"].collapsedPosition()
else:
if scale is not None and scale != 1:
@@ -583,7 +586,7 @@ class Segment(ConnectableElement, Group):
for a in atoms:
if a.name[-1] in ("'","P","T"):
a.position = scale*(a.position-r0) + r0
- a.beta = 0
+ a.fixed = 1
atoms.position = pos
atoms.contour_position = contour_position
@@ -1321,14 +1324,15 @@ class Strand(Group):
# pdb.set_trace()
# assert(s.end != s.start)
assert( s.num_nt == 1 or (np.linalg.norm( seg.contour_to_position(contour[-1]) - seg.contour_to_position(contour[0]) ) > 0.1) )
+ nucleotide_count = 0
for c,seq in zip(contour,s.get_sequence()):
+ nucleotide_count += 1
if last is None and not self.is_circular:
seq = "5"+seq
- if strand_segment_count == len(s.strand_segments) and c == 1 and not self.is_circular:
+ if strand_segment_count == len(s.strand_segments) and nucleotide_count == s.num_nt and not self.is_circular:
seq = seq+"3"
nt = seg._generate_atomic_nucleotide( c, s.is_fwd, seq, scale, s )
- s.add(nt)
## Join last basepairs
if last is not None:
@@ -1515,6 +1519,7 @@ class SegmentModel(ArbdModel):
for seg in self.segments:
for d in ('fwd','rev'):
seg.strand_pieces[d] = []
+ self._generate_strands()
def clear_beads(self):
return self._clear_beads()
@@ -1525,7 +1530,8 @@ class SegmentModel(ArbdModel):
s.clear_all()
self.clear_all(keep_children=True)
- self.clear_atomic()
+ if len(self.strands[0].children[0].children) > 0:
+ self.clear_atomic()
## Check that it worked
assert( len([b for b in self]) == 0 )
@@ -2355,6 +2361,9 @@ class SegmentModel(ArbdModel):
try:
for s in self.strands:
self.children.remove(s)
+ for seg in self.segments:
+ for d in ('fwd','rev'):
+ seg.strand_pieces[d] = []
except:
...
self.strands = strands = []
@@ -2674,7 +2683,7 @@ class SegmentModel(ArbdModel):
# print("WARNING: could not find 'P' atom in {}:{} or {}:{}".format( segI, ntI_idx, segJ, ntJ_idx ))
...
- print("PUSH BONDS:", len(push_bonds))
+ # print("PUSH BONDS:", len(push_bonds))
if not self.useTclForces:
with open("%s.exb" % output_name, 'a') as fh:
--
GitLab