diff --git a/arbdmodel.py b/arbdmodel.py
index fb2ed3167114e2e662e2b27b3c3c47b06cf5d552..aab8e96e06f557f5318182509c40a97f80dd22fc 100644
--- a/arbdmodel.py
+++ b/arbdmodel.py
@@ -81,6 +81,16 @@ class Parent():
x.parent = self
self.children.append(x)
+ def clear_all(self, keep_children=False):
+ if keep_children == False:
+ for x in self.children:
+ x.parent = None
+ self.children = []
+ self.bonds = []
+ self.angles = []
+ self.dihedrals = []
+ self.exclusions = []
+
def remove(self,x):
if x in self.children:
self.children.remove(x)
@@ -462,6 +472,13 @@ class ArbdModel(PdbModel):
self.cacheUpToDate = False
+ def clear_all(self, keep_children=False):
+ Parent.clear_all(self, keep_children=keep_children)
+ self.particles = []
+ self.numParticles = 0
+ self.type_counts = None
+ self._nbParamFiles = []
+ self._written_bond_files = dict()
def _getNbScheme(self, typeA, typeB):
scheme = None
diff --git a/coords.py b/coords.py
index 0a4d4fb38c075904801be00b415e92dab06e4fc9..368b32d70d78ad2bc0cc80d018d4064adfbebed7 100644
--- a/coords.py
+++ b/coords.py
@@ -13,7 +13,7 @@ def minimizeRmsd(coordsB, coordsA, weights=None, maxIter=100):
while tol > 1e-6:
q,cB,comA = _minimizeRmsd(cNext,coordsA, weights)
- R = R.dot(quaternionToMatrix3(q))
+ R = R.dot(quaternion_to_matrix(q))
assert( np.all(np.isreal( R )) )
comB += cB
@@ -32,7 +32,7 @@ def minimizeRmsd(coordsB, coordsA, weights=None, maxIter=100):
def minimizeRmsd(coordsB, coordsA, weights=None):
q,comA,comB = _minimizeRmsd(coordsB, coordsA, weights)
assert( np.all(np.isreal( q )) )
- return quaternionToMatrix3(q),comA,comB
+ return quaternion_to_matrix(q),comA,comB
## http://onlinelibrary.wiley.com/doi/10.1002/jcc.21439/full
@@ -92,7 +92,7 @@ def _minimizeRmsd(coordsB, coordsA, weights=None):
# print("CHECK:", K.dot(q)-vals[i]*q )
return q, comB, comA
-def quaternionToMatrix3(q):
+def quaternion_to_matrix(q):
assert(len(q) == 4)
## It looks like the wikipedia article I used employed a less common convention for q (see below
@@ -109,10 +109,48 @@ def quaternionToMatrix3(q):
return np.array(R)
+def quaternion_from_matrix( R ):
+ e1 = R[0]
+ e2 = R[1]
+ e3 = R[2]
+
+ # d1 = 0.5 * np.sqrt( 1+R[0,0]+R[1,1]+R[2,2] )
+ # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
+ # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
+ # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] )
+
+ d1 = 1+R[0,0]+R[1,1]+R[2,2]
+ d2 = 1+R[0,0]-R[1,1]-R[2,2]
+ d2 = 1+R[0,0]-R[1,1]-R[2,2]
+ d2 = 1+R[0,0]-R[1,1]-R[2,2]
+
+ # minD = min((d1,d2,d3,d3))
+ maxD = max((d1,d2))
+
+ ## TODO: add more approaches TypeError: 0 < idim < 11 must hold
+ d = 0.5 / np.sqrt(maxD)
+ if d1 == maxD:
+ return np.array(( 1.0/(4*d),
+ d * (R[2,1]-R[1,2]),
+ d * (R[0,2]-R[2,0]),
+ d * (R[1,0]-R[0,1]) ))
+ elif d2 == maxD:
+ return np.array(( d * (R[2,1]-R[1,2]),
+ 1.0/(4*d),
+ d * (R[0,1]+R[1,0]),
+ d * (R[0,2]+R[2,0]) ))
+
def rotationAboutAxis(axis,angle, normalizeAxis=True):
if normalizeAxis: axis = axis / np.linalg.norm(axis)
angle *= 0.5 * np.pi/180
cos = np.cos( angle )
sin = np.sin( angle )
q = [cos] + [sin*x for x in axis]
- return quaternionToMatrix3(q)
+ return quaternion_to_matrix(q)
+
+def readArbdCoords(fname):
+ coords = []
+ with open(fname) as fh:
+ for line in fh:
+ coords.append([float(x) for x in line.split()[1:]])
+ return np.array(coords)
diff --git a/segmentmodel.py b/segmentmodel.py
index a67bfdd2f5e59f90a6f4b495baa34bd8b3f88794..b3c14caee58329f4a71f2ba48759b9108d7e5b0e 100644
--- a/segmentmodel.py
+++ b/segmentmodel.py
@@ -1,13 +1,15 @@
import numpy as np
from arbdmodel import PointParticle, ParticleType, Group, ArbdModel
-from coords import rotationAboutAxis
+from coords import rotationAboutAxis, quaternion_from_matrix, quaternion_to_matrix
from nonbonded import *
from copy import copy, deepcopy
from nbPot import nbDnaScheme
from scipy.special import erf
import scipy.optimize as opt
+from scipy import interpolate
+import types
"""
TODO:
@@ -94,6 +96,41 @@ class Segment(ConnectableElement, Group):
self.start_position = start_position
self.end_position = end_position
+ ## Set up interpolation for positions
+ a = np.array([self.start_position,self.end_position]).T
+ tck, u = interpolate.splprep( a, u=[0,1], s=0, k=1)
+ self.position_spline_params = tck
+
+
+ def contour_to_position(self,s):
+ p = interpolate.splev( s, self.position_spline_params )
+ if len(p) > 1: p = np.array(p).T
+ return p
+
+ def contour_to_tangent(self,s):
+ t = interpolate.splev( s, self.position_spline_params, der=1 )
+ if len(t) > 1: t = np.array(t).T
+ return t
+
+
+ def contour_to_orientation(self,s):
+ if self.start_orientation is not None:
+ # axis = self.start_orientation.dot( np.array((0,0,1)) )
+ if self.quaternion_spline_params is None:
+ axis = self.contour_to_tangent(s)
+ # orientation = rotationAboutAxis( axis, s*self.twist_per_nt*self.num_nts, normalizeAxis=False )
+ orientation = rotationAboutAxis( axis, s*self.twist_per_nt*self.num_nts, normalizeAxis=False )
+ ## TODO: ensure this is correct
+ # orientation = self.start_orientation.dot(orientation) # .dot( self.start_orientation )
+ orientation = orientation.dot( self.start_orientation )
+ else:
+ q = interpolate.splev(s, self.quaternion_spline_params)
+ orientation = quaternion_to_matrix(q)
+ else:
+ orientation = None
+ return orientation
+
+
def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead ):
raise NotImplementedError
@@ -143,13 +180,11 @@ class Segment(ConnectableElement, Group):
""" Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """
## TODO: decide whether to remove bead_model argument
+ ## (currently unused)
# self._bead_model_generation += 1
# self._bead_model_max_nts_per_bead = max_nts_per_bead
- direction = self.end_position - self.start_position
-
-
num_beads = self._get_num_beads( max_basepairs_per_bead, max_nucleotides_per_bead )
nts_per_bead = float(self.num_nts)/num_beads
twist_per_bead = nts_per_bead * self.twist_per_nt
@@ -159,7 +194,8 @@ class Segment(ConnectableElement, Group):
if num_beads <= 2:
## not yet implemented for dsDNA
assert( isinstance(self, SingleStrandedSegment) )
- b = self._generate_one_bead(self.start_position,
+ b = self._generate_one_bead(0.5,
+ self.start_position,
self.num_nts,
self.start_orientation)
self.children.append(b)
@@ -167,7 +203,6 @@ class Segment(ConnectableElement, Group):
self._assign_particles_to_locations()
return
- assert( np.linalg.norm(direction) > 0 )
for i in range(num_beads+1):
nts = nts_per_bead
@@ -175,7 +210,7 @@ class Segment(ConnectableElement, Group):
nts *= 0.5
s = i*float(nts_per_bead)/(self.num_nts) # contour
- pos = direction * s + self.start_position
+ pos = self.contour_to_position(s)
if self.start_orientation is not None:
axis = self.start_orientation.dot( np.array((0,0,1)) )
orientation = rotationAboutAxis( axis, s*self.twist_per_nt*self.num_nts, normalizeAxis=False )
@@ -185,7 +220,7 @@ class Segment(ConnectableElement, Group):
else:
orientation = None
- b = self._generate_one_bead(pos,nts,orientation)
+ b = self._generate_one_bead(s,pos,nts,orientation)
self.children.append(b)
self.beads.append(b) # don't add orientation bead
@@ -242,6 +277,15 @@ class DoubleStrandedSegment(Segment):
self.end5 = Location( self, address=-1, type_= "end5" )
self.end3 = Location( self, address=-1, type_ = "end3" )
+ ## Set up interpolation for azimuthal angles
+ a = np.array([self.start_position,self.end_position]).T
+ tck, u = interpolate.splprep( a, u=[0,1], s=0, k=1)
+ self.position_spline_params = tck
+
+ ## TODO: initialize sensible spline for orientation
+ self.quaternion_spline_params = None
+
+
## Convenience methods
def connect_start5(self, end3, type_="intrahelical", force_connection=False):
if isinstance(end3, SingleStrandedSegment):
@@ -275,7 +319,7 @@ class DoubleStrandedSegment(Segment):
def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
return (self.num_nts // max_basepairs_per_bead) + 1
- def _generate_one_bead(self, pos, nts, orientation=None):
+ def _generate_one_bead(self, contour_position, pos, nts, orientation=None):
if self.local_twist:
assert(orientation is not None)
# opos = pos + np.array((2,0,0)).dot(orientation)
@@ -283,11 +327,14 @@ class DoubleStrandedSegment(Segment):
o = PointParticle( Segment.orientation_particle, opos, nts,
num_nts=nts, parent=self )
bead = PointParticle( Segment.dsDNA_particle, pos, nts,
- num_nts=nts, parent=self, orientation_bead=o )
+ num_nts=nts, parent=self,
+ orientation_bead=o,
+ contour_position=contour_position )
else:
bead = PointParticle( Segment.dsDNA_particle, pos, nts,
- num_nts=nts, parent=self )
+ num_nts=nts, parent=self,
+ contour_position=contour_position )
return bead
@@ -359,9 +406,10 @@ class SingleStrandedSegment(Segment):
# pdb.set_trace()
return (self.num_nts // max_nucleotides_per_bead) + 1
- def _generate_one_bead(self, pos, nts, orientation=None):
+ def _generate_one_bead(self, contour_position, pos, nts, orientation=None):
return PointParticle( Segment.ssDNA_particle, pos, nts,
- num_nts=nts, parent=self )
+ num_nts=nts, parent=self,
+ contour_position=contour_position )
def _assign_particles_to_locations(self):
self.start.particle = self.children[0]
@@ -392,7 +440,7 @@ class SegmentModel(ArbdModel):
self._bonded_potential = dict() # cache bonded potentials
- self._generate_bead_model(segments, max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
+ self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
self.useNonbondedScheme( nbDnaScheme )
@@ -526,14 +574,72 @@ class SegmentModel(ArbdModel):
fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - np.exp(-sep/Lp)
k = opt.leastsq( fitFun, x0=np.exp(-sep/Lp) )
return k[0][0] * 2*kT*0.00030461742
+
+ # def _update_segment_positions(self, bead_coordinates):
+ # """ Set new function for each segments functions
+ # contour_to_position and contour_to_orientation """
+
+ # dsDnaHelixNeighborDist=50
+ # dsDnaAllNeighborDist=30
+ # ssDnaHelixNeighborDist=25
+ # ssDnaAllNeighborDist=25
+ # beads = [b in s.beads for s in self.segments]
+ # positions = np.array([b.position for b in beads])
+ # neighborhood = dict()
+
+ # ## Assign neighborhood to each bead
+ # for b in beads:
+ # dists = b.position[np.newaxis,:] - positions
+ # dists = np.linalg.norm(dists, axis=-1)
+ # neighborhood[b] = np.where( dists < 50 )
+
+ """ Mapping between different resolution models """
+ def _clear_beads(self):
+ for s in self.segments:
+ s.clear_all()
+ s.beads = []
+ self.clear_all(keep_children=True)
+
+ def _update_segment_positions(self, bead_coordinates):
+ """ Set new function for each segments functions
+ contour_to_position and contour_to_orientation """
- def _generate_bead_model(self, segments,
+ for s in self.segments:
+ beads = [b for b in s.beads]
+ ids = [b.idx for b in beads]
+
+ """ Get positions """
+ positions = bead_coordinates[ids,:].T
+ contours = [b.contour_position for b in beads]
+ tck, u = interpolate.splprep( positions, u=contours, s=0, )
+
+ s.position_spline_params = tck
+
+ """ Get twist """
+ if 'orientation_bead' in beads[0].__dict__:
+ tangents = s.contour_to_tangent(contours)
+ quats = []
+ for b,t in zip(beads,tangents):
+ o = b.orientation_bead
+ angleVec = o.position - b.position
+ angleVec = angleVec - angleVec.dot(t)*t
+ angleVec = angleVec/np.linalg.norm(angleVec)
+ y = np.cross(t,angleVec)
+ quats.append( quaternion_from_matrix( np.array([t,y,angleVec])) )
+ quats = np.array(quats)
+ tck, u = interpolate.splprep( quats.T, u=contours, s=0, )
+ s.quaternion_spline_params = tck
+
+
+ ## TODO: set twist
+
+ def _generate_bead_model(self,
max_basepairs_per_bead = 7,
max_nucleotides_per_bead = 4,
- local_twist=False
- ):
+ local_twist=False):
+ segments = self.segments
""" Generate beads """
if self.DEBUG: print("Generating beads")
@@ -715,6 +821,9 @@ class SegmentModel(ArbdModel):
pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot)
+ self._updateParticleOrder()
+
+
# def get_bead(self, location):
# if type(location.container) is not list:
# s = self.segments.index(location.container)