diff --git a/mrdna/coords.py b/mrdna/coords.py index b74906da417e0b7cc787be05fbdb422e24609705..41628e83501c17098094fc12b7d4d4e6594ebcb6 100644 --- a/mrdna/coords.py +++ b/mrdna/coords.py @@ -150,7 +150,7 @@ def quaternion_from_matrix( R ): def rotationAboutAxis(axis,angle, normalizeAxis=True): if normalizeAxis: axis = axis / np.linalg.norm(axis) - angle *= 0.5 * np.pi/180 + angle = angle * 0.5 * np.pi/180 cos = np.cos( angle ) sin = np.sin( angle ) q = [cos] + [sin*x for x in axis] diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py index 068271454fe2eb5c256956f68e5a5f4cf4ae3865..61ddfca9e9ddfa57e7459941503bbdebdbfab7cf 100644 --- a/mrdna/segmentmodel.py +++ b/mrdna/segmentmodel.py @@ -420,7 +420,7 @@ class Segment(ConnectableElement, Group): ids = list(filter(lambda i: position_filter(r[i,:]), ids)) return ids - def translate(self, translation_vector, about=None, position_filter=None, contour_filter=None): + def translate(self, translation_vector, position_filter=None, contour_filter=None): dr = np.array(translation_vector) tck, u = self.position_spline_params r = self.contour_to_position(u) @@ -429,7 +429,7 @@ class Segment(ConnectableElement, Group): ## Translate r[ids,:] = r[ids,:] + dr[np.newaxis,:] - self.set_position_splines(r,u) + self.set_splines(u,r) def rotate(self, rotation_matrix, about=None, position_filter=None, contour_filter=None): tck, u = self.position_spline_params @@ -445,16 +445,16 @@ class Segment(ConnectableElement, Group): ## TODO: do this more efficiently r[ids,:] = np.array([rotation_matrix.dot(r[i,:]-dr) + dr for i in range(r.shape[0])]) - self.set_position_splines(r,u) + self.set_splines(u,r) if self.quaternion_spline_params is not None: ## TODO: performance: don't shift between quaternion and matrix representations so much tck, u = self.quaternion_spline_params - orientations = [self.contour_to_orientation(v)) for v in u] + orientations = [self.contour_to_orientation(v) for v in u] for i in ids: orientations[i,:] = rotation_matrix.dot(orientations[i]) quats = [quaternion_from_matrix(o) for o in orientations] - self.set_orientation_splines(quats,u) + self.set_orientation_splines(u, quats) def _set_splines_from_ends(self): self.quaternion_spline_params = None @@ -1386,6 +1386,7 @@ class SegmentModel(ArbdModel): self._bonded_potential = dict() # cache for bonded potentials self._generate_strands() + self.potentials = [] self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist, escapable_twist) self.useNonbondedScheme( nbDnaScheme ) @@ -1540,6 +1541,26 @@ class SegmentModel(ArbdModel): print("WARNING: called deprecated command '_update_segment_positions; use 'update_splines' instead") return self.update_splines(bead_coordinates) + ## Operations on spline coordinates + def translate(self, translation_vector, position_filter=None): + for s in self.segments: + s.translate(translation_vector, position_filter=position_filter) + + def rotate(self, rotation_matrix, about=None, position_filter=None): + for s in self.segments: + s.rotate(rotation_matrix, about=about, position_filter=position_filter) + + def get_center(self, include_ssdna=False): + if include_ssdna: + segments = self.segments + else: + segments = list(filter(lambda s: isinstance(s,DoubleStrandedSegment), + self.segments)) + centers = [s.get_center() for s in segments] + weights = [s.num_nt*2 if isinstance(s,DoubleStrandedSegment) else s.num_nt for s in segments] + # centers,weights = [np.array(a) for a in (centers,weights)] + return np.average( centers, axis=0, weights=weights) + def update_splines(self, bead_coordinates): """ Set new function for each segments functions contour_to_position and contour_to_orientation """ @@ -1626,14 +1647,8 @@ class SegmentModel(ArbdModel): # print(" {}[{:03f}]: {:0.3f}".format(b.parent.name, b.contour_position, c) ) - ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 ) - tck = ret[0][0] + tck, u = interpolate.splprep( positions, u=contours, s=0, k=1 ) - ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 ) - tck2 = ret[0][0] - - assert( tck[0][0] == tck2[0][0] ) - # if len(beads) < 8: # ret = interpolate.splprep( positions, u=contours, s=0, k=1, full_output=1 ) # tck = ret[0][0] @@ -1652,7 +1667,7 @@ class SegmentModel(ArbdModel): # if ret[2] > 0: # pdb.set_trace() - s.position_spline_params = (tck,ret[1][0]) + s.position_spline_params = (tck,u) """ Get twist """ cb = [e for e in zip(contours,beads) if 'orientation_bead' in e[1].__dict__]