Skip to content

GitLab

Commit a36bc78c authored by cmaffeo2's avatar cmaffeo2
Browse files

Initial implementation of splined positions/orientations in segmentModel

parent 6457307f
Hide whitespace changes
Inline Side-by-side
arbdmodel.py
View file @ a36bc78c
......@@ -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
......
coords.py
View file @ a36bc78c
......@@ -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)
segmentmodel.py
View file @ a36bc78c
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)
......
Markdown is supported
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