Source code for arbdmodel.polymer

from abc import abstractmethod, ABCMeta

from .logger import logger

from pathlib import Path
from copy import deepcopy

import numpy as np
from scipy import interpolate

from .model import ArbdModel
from .coords import rotationAboutAxis, quaternion_from_matrix, quaternion_to_matrix
from . import Group, PointParticle
from .interactions import HarmonicBond

"""
TODO:
 - test for large systems
 - rework Location class 
 - remove recursive calls
 - document
 - develop unit test suite
"""

# class ParticleNotConnectedError(Exception):
#     pass

[docs] class Location(): """ Site for connection within an object """ def __init__(self, parent, contour_position, type_, on_fwd_strand = True): ## TODO: remove cyclic references(?) assert( isinstance(parent, ConnectableElement) ) self.parent = parent self.contour_position = contour_position # represents position along contour length in polymer self.type_ = type_ # self.particle = None self.connection = None
[docs] def get_connected_location(self): if self.connection is None: return None else: return self.connection.other(self)
[docs] def set_connection(self, connection): self.connection = connection # TODO weakref?
[docs] def get_monomer_index(self): try: idx = self.parent.contour_to_monomer_index(self.contour_position, nearest_monomer=True) except: if self.contour_position == 0: idx = 0 elif self.contour_position == 1: idx = self.parent.num_monomers-1 else: raise return idx
def __repr__(self): if self.on_fwd_strand: on_fwd = "on_fwd_strand" else: on_fwd = "on_rev_strand" return "<Location {}.{}[{:.2f},{:d}]>".format( self.parent.name, self.type_, self.contour_position, self.on_fwd_strand)
[docs] class Connection(): """ Class for connecting two elements """ def __init__(self, A, B, type_ = None): assert( isinstance(A,Location) ) assert( isinstance(B,Location) ) self.A = A self.B = B self.type_ = type_
[docs] def other(self, location): if location is self.A: return self.B elif location is self.B: return self.A else: raise ValueError
[docs] def delete(self): self.A.parent.connections.remove(self) if self.B.parent is not self.A.parent: self.B.parent.connections.remove(self) self.A.connection = None self.B.connection = None
def __repr__(self): return "<Connection {}--{}--{}]>".format( self.A, self.type_, self.B )
[docs] class ConnectableElement(): """ ConnectableElement is an abstract base class that represents an object that can be connected to other objects. This class provides functionality for managing connections between objects, where connections are typically made between 'Location' objects that are associated with the ConnectableElement. Attributes: locations (list): A list of Location objects associated with this element. Alias for connection_locations. connection_locations (list): A list of Location objects associated with this element. connections (list): A list of Connection objects that this element is a part of. Methods: get_locations(type_=None, exclude=()): Returns locations of a specified type, excluding certain types. get_location_at(contour_position, on_fwd_strand=True, new_type="crossover"): Returns a location at a specific contour position and strand orientation, creating one if necessary. get_connections_and_locations(connection_type=None, exclude=()): Returns a list of [connection, location_in_self, location_in_other] for each connection. Note: - The class is designed to be subclassed. - It appears to be part of a system for modeling polymer structures, possibly DNA origami or similar structures. - The class enforces that objects cannot be connected to themselves directly. """ def __init__(self, connection_locations=None, connections=None): if connection_locations is None: connection_locations = [] if connections is None: connections = [] ## TODO decide on names self.locations = self.connection_locations = connection_locations self.connections = connections
[docs] def get_locations(self, type_=None, exclude=()): locs = [l for l in self.connection_locations if (type_ is None or l.type_ == type_) and l.type_ not in exclude] counter = dict() for l in locs: if l in counter: counter[l] += 1 else: counter[l] = 1 assert( np.all( [counter[l] == 1 for l in locs] ) ) return locs
[docs] def get_location_at(self, contour_position, on_fwd_strand=True, new_type="crossover"): loc = None if (self.num_monomers == 1): ## Assumes that intrahelical connections have been made before crossovers for l in self.locations: if l.on_fwd_strand == on_fwd_strand and l.connection is None: assert(loc is None) loc = l # assert( loc is not None ) else: for l in self.locations: if l.contour_position == contour_position and l.on_fwd_strand == on_fwd_strand: assert(loc is None) loc = l if loc is None: loc = Location( self, contour_position=contour_position, type_=new_type, on_fwd_strand=on_fwd_strand ) return loc
[docs] def get_connections_and_locations(self, connection_type=None, exclude=()): """ Returns a list with each entry of the form: connection, location_in_self, location_in_other """ type_ = connection_type ret = [] for c in self.connections: if (type_ is None or c.type_ == type_) and c.type_ not in exclude: if c.A.parent is self: ret.append( [c, c.A, c.B] ) elif c.B.parent is self: ret.append( [c, c.B, c.A] ) else: import pdb pdb.set_trace() raise Exception("Object contains connection that fails to refer to object") return ret
def _connect(self, other, connection): ## TODO fix circular references A,B = [connection.A, connection.B] A.connection = B.connection = connection self.connections.append(connection) if other is not self: other.connections.append(connection) else: raise NotImplementedError("Objects cannot yet be connected to themselves; if you are attempting to make a circular object, try breaking the object into multiple components") l = A.parent.locations if A not in l: l.append(A) l = B.parent.locations if B not in l: l.append(B)
[docs] class PolymerSection(ConnectableElement): """ Base class that describes a linear section of a polymer. A PolymerSection represents a continuous section of a polymer, defined by the number of monomers and their spatial properties. It handles spatial transformations, contour-to-position mapping, and manages connections to other polymer elements. Attributes: segname (str): Name of the polymer section. num_monomers (int): Number of monomers in this polymer section. monomer_length (float): Length of each monomer. start_position (ndarray): 3D coordinates of the section's starting point. end_position (ndarray): 3D coordinates of the section's ending point. start_orientation (ndarray, optional): Initial orientation matrix. position_spline_params (tuple): Spline parameters for position interpolation. quaternion_spline_params (tuple, optional): Spline parameters for orientation interpolation. sequence (list, optional): Sequence of monomers if applicable. """ def __init__(self, name, num_monomers, monomer_length, start_position = None, end_position = None, parent = None, **kwargs): ConnectableElement.__init__(self, connection_locations=[], connections=[]) if 'segname' not in kwargs: self.segname = name for key,val in kwargs.items(): self.__dict__[key] = val self.num_monomers = int(num_monomers) self.monomer_length = monomer_length if start_position is None: start_position = np.array((0,0,0)) if end_position is None: end_position = np.array((0,0,self.monomer_length*num_monomers)) + start_position self.start_position = start_position self.end_position = end_position self.start_orientation = None ## Set up interpolation for positions self._set_splines_from_ends() # self.sequence = None def __repr__(self): return "<{} {}[{:d}]>".format( type(self), self.segname, self.num_monomers )
[docs] def set_splines(self, contours, coords): """ Sets up splines for the polymer's position. This method creates cubic splines that interpolate the polymer's position coordinates along its contour. These splines can be used for smooth interpolation between discrete polymer positions. Parameters ---------- contours : array_like Contour parameter values corresponding to each coordinate point. Typically represents positions along the polymer chain. coords : ndarray Coordinates of the polymer, shape (n_points, n_dimensions). Each row represents the position of a polymer segment in space. Note ----- The method uses scipy.interpolate.splprep with linear splines (k=1) and no smoothing (s=0). The resulting spline parameters are stored in self.position_spline_params as a tuple containing the knot points, coefficients, and degree of the spline (tck) and the parameter values (u). """ tck, u = interpolate.splprep( coords.T, u=contours, s=0, k=1) self.position_spline_params = (tck,u)
[docs] def set_orientation_splines(self, contours, quaternions): """ Set spline interpolation parameters for orientation. Creates a linear interpolation spline (k=1) for quaternion orientations along contour points. The spline parameters are stored in the self.quaternion_spline_params attribute. Parameters ---------- contours : array_like Parametric positions along the polymer contour where orientations are defined. quaternions : array_like Array of quaternions defining orientations at each contour point. Should be of shape (n, 4) where n is the number of contour points. Note ----- Uses scipy.interpolate.splprep with smoothing factor s=0 (exact interpolation) and k=1 (linear interpolation). """ tck, u = interpolate.splprep( quaternions.T, u=contours, s=0, k=1) self.quaternion_spline_params = (tck,u)
[docs] def get_center(self): tck, u = self.position_spline_params return np.mean(self.contour_to_position(u), axis=0)
def _get_location_positions(self): return [self.contour_to_monomer_index(l.contour_position) for l in self.locations]
[docs] def insert_monomers(self, at_monomer: int, num_monomer: int, seq=tuple()): assert(np.isclose(np.around(num_monomer),num_monomer)) if at_monomer < 0: raise ValueError("Attempted to insert DNA into {} at a negative location".format(self)) if at_monomer > self.num_monomer-1: raise ValueError("Attempted to insert DNA into {} at beyond the end of the Polymer".format(self)) if num_monomers < 0: raise ValueError("Attempted to insert DNA a negative amount of DNA into {}".format(self)) num_monomers = np.around(num_monomer) monomer_positions = self._get_location_positions() new_monomer_positions = [p if p <= at_monomer else p+num_monomers for p in monomer_positions] ## TODO: handle sequence self.num_monomers = self.num_monomer+num_monomer for l,p in zip(self.locations, new_monomer_positions): l.contour_position = self.monomer_index_to_contour(p)
[docs] def remove_monomers(self, first_monomer: int, last_monomer: int): """ Removes nucleotides between first_monomer and last_monomer, inclusive """ assert(np.isclose(np.around(first_monomer),first_monomer)) assert(np.isclose(np.around(last_monomer),last_monomer)) tmp = min((first_monomer,last_monomer)) last_monomer = max((first_monomer,last_monomer)) fist_monomer = tmp if first_monomer < 0 or first_monomer > self.num_monomer-2: raise ValueError("Attempted to remove DNA from {} starting at an invalid location {}".format(self, first_monomer)) if last_monomer < 1 or last_monomer > self.num_monomer-1: raise ValueError("Attempted to remove DNA from {} ending at an invalid location {}".format(self, last_monomer)) if first_monomer == last_monomer: return first_monomer = np.around(first_monomer) last_monomer = np.around(last_monomer) monomer_positions = self._get_location_positions() bad_locations = list(filter(lambda p: p >= first_monomer and p <= last_monomer, monomer_positions)) if len(bad_locations) > 0: raise Exception("Attempted to remove DNA containing locations {} from {} between {} and {}".format(bad_locations,self,first_monomer,last_monomer)) removed_monomer = last_monomer-first_monomer+1 new_monomer_positions = [p if p <= last_monomer else p-removed_monomer for p in monomer_positions] num_monomers = self.num_monomer-removed_monomer if self.sequence is not None and len(self.sequence) == self.num_monomer: self.sequence = [s for s,i in zip(self.sequence,range(self.num_monomer)) if i < first_monomer or i > last_monomer] assert( len(self.sequence) == num_monomers ) self.num_monomers = num_monomers for l,p in zip(self.locations, new_monomer_positions): l.contour_position = self.monomer_position_to_contour(p)
def __filter_contours(contours, positions, position_filter, contour_filter): u = contours r = positions ## Filter ids = list(range(len(u))) if contour_filter is not None: ids = list(filter(lambda i: contour_filter(u[i]), ids)) if position_filter is not None: ids = list(filter(lambda i: position_filter(r[i,:]), ids)) return ids
[docs] 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) ids = PolymerSection.__filter_contours(u, r, position_filter, contour_filter) if len(ids) == 0: return ## Translate r[ids,:] = r[ids,:] + dr[np.newaxis,:] self.set_splines(u,r)
[docs] def rotate(self, rotation_matrix, about=None, position_filter=None, contour_filter=None): tck, u = self.position_spline_params r = self.contour_to_position(u) ids = PolymerSection.__filter_contours(u, r, position_filter, contour_filter) if len(ids) == 0: return if about is None: ## TODO: do this more efficiently r[ids,:] = np.array([rotation_matrix.dot(r[i,:]) for i in ids]) else: dr = np.array(about) ## TODO: do this more efficiently r[ids,:] = np.array([rotation_matrix.dot(r[i,:]-dr) + dr for i in ids]) 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] for i in ids: orientations[i,:] = rotation_matrix.dot(orientations[i]) quats = [quaternion_from_matrix(o) for o in orientations] self.set_orientation_splines(u, quats)
def _set_splines_from_ends(self, resolution=4): self.quaternion_spline_params = None r0 = np.array(self.start_position)[np.newaxis,:] r1 = np.array(self.end_position)[np.newaxis,:] u = np.linspace(0,1, max(3,self.num_monomers//int(resolution))) s = u[:,np.newaxis] coords = (1-s)*r0 + s*r1 self.set_splines(u, coords)
[docs] def contour_to_monomer_index(self, contour_pos, nearest_monomer=False): idx = contour_pos*(self.num_monomer) - 0.5 if nearest_monomer: assert( np.isclose(np.around(idx),idx) ) idx = np.around(idx) return idx
[docs] def monomer_index_to_contour(self,monomer_pos): return (monomer_pos+0.5)/(self.num_monomers)
[docs] def contour_to_position(self,s): p = interpolate.splev( s, self.position_spline_params[0] ) if len(p) > 1: p = np.array(p).T return p
[docs] def contour_to_tangent(self,s): t = interpolate.splev( s, self.position_spline_params[0], der=1 ) t = (t / np.linalg.norm(t,axis=0)) return t.T
[docs] def contour_to_orientation(self,s): assert( isinstance(s,float) or isinstance(s,int) or len(s) == 1 ) # TODO make vectorized version if self.quaternion_spline_params is None: axis = self.contour_to_tangent(s) axis = axis / np.linalg.norm(axis) rotAxis = np.cross(axis,np.array((0,0,1))) rotAxisL = np.linalg.norm(rotAxis) zAxis = np.array((0,0,1)) if rotAxisL > 0.001: theta = np.arcsin(rotAxisL) * 180/np.pi if axis.dot(zAxis) < 0: theta = 180-theta orientation0 = rotationAboutAxis( rotAxis/rotAxisL, theta, normalizeAxis=False ).T 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_monomer_index(s), normalizeAxis=False ) # orientation = orientation.dot(orientation0) orientation = orientation0 else: q = interpolate.splev( s, self.quaternion_spline_params[0] ) if len(q) > 1: q = np.array(q).T # TODO: is this needed? orientation = quaternion_to_matrix(q) return orientation
[docs] def get_contour_sorted_connections_and_locations(self,type_): sort_fn = lambda c: c[1].contour_position cl = self.get_connections_and_locations(type_) return sorted(cl, key=sort_fn)
[docs] def add_location(self, nt, type_, on_fwd_strand=True): ## Create location if needed, add to polymer c = self.monomer_index_to_contour(nt) assert(c >= 0 and c <= 1) # TODO? loc = self.Location( contour_position=c, type_=type_, on_fwd_strand=is_fwd ) loc = Location( self, contour_position=c, type_=type_, on_fwd_strand=on_fwd_strand ) self.locations.append(loc)
[docs] def iterate_connections_and_locations(self, reverse=False): ## connections to other polymers cl = self.get_contour_sorted_connections_and_locations() if reverse: cl = cl[::-1] for c in cl: yield c
[docs] class PolymerGroup(): def __init__(self, polymers=[], **kwargs): """ Simulation system parameters """ ## TODO decide whether to move these elsewhere # self.dimensions = dimensions # self.temperature = temperature # self.timestep = timestep # self.cutoff = cutoff # self.decomposition_period = decomposition_period # self.pairlist_distance = pairlist_distance # self.nonbonded_resolution = nonbonded_resolution self.polymers = polymers self.grid_potentials = [] # self.useNonbondedScheme( nbDnaScheme )
[docs] def get_connections(self,type_=None,exclude=()): """ Find all connections in model, without double-counting """ added=set() ret=[] for s in self.polymers: items = [e for e in s.get_connections_and_locations(type_,exclude=exclude) if e[0] not in added] added.update([e[0] for e in items]) ret.extend( list(sorted(items,key=lambda x: x[1].contour_position)) ) return ret
[docs] def extend(self, other, copy=True, include_strands=False): assert( isinstance(other, PolymerSectionModel) ) if copy: for s in other.polymers: self.polymers.append(deepcopy(s)) if include_strands: for s in other.strands: self.strands.append(deepcopy(s)) else: for s in other.polymers: self.polymers.append(s) if include_strands: for s in other.strands: self.strands.append(s) self._clear_beads()
[docs] def update(self, polymer , copy=False): assert( isinstance(polymer, PolymerSection) ) if copy: polymer = deepcopy(polymer) self.polymers.append(polymer)
""" Operations on spline coordinates """
[docs] def translate(self, translation_vector, position_filter=None): for s in self.polymers: s.translate(translation_vector, position_filter=position_filter)
[docs] def rotate(self, rotation_matrix, about=None, position_filter=None): for s in self.polymers: s.rotate(rotation_matrix, about=about, position_filter=position_filter)
[docs] def get_center(self, include_ssdna=False): if include_ssdna: polymers = self.polymers else: polymers = list(filter(lambda s: isinstance(s,DoubleStrandedPolymerSection), self.polymers)) centers = [s.get_center() for s in polymers] weights = [s.num_monomers*2 if isinstance(s,DoubleStrandedPolymerSection) else s.num_monomers for s in polymers] # centers,weights = [np.array(a) for a in (centers,weights)] return np.average( centers, axis=0, weights=weights)
[docs] def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ): positions = [] for s in self.polymers: positions.append(s.contour_to_position(0)) positions.append(s.contour_to_position(0.5)) positions.append(s.contour_to_position(1)) positions = np.array(positions) dx,dy,dz = [(np.max(positions[:,i])-np.min(positions[:,i])+30)*padding_factor for i in range(3)] if isotropic: dx = dy = dz = max((dx,dy,dz)) return [dx,dy,dz]
[docs] def add_grid_potential(self, grid_file, scale=1, per_monomer=True): grid_file = Path(grid_file) if not grid_file.is_file(): raise ValueError("Grid file {} does not exist".format(grid_file)) if not grid_file.is_absolute(): grid_file = Path.cwd() / grid_file self.grid_potentials.append((grid_file,scale,per_,pmp,er))
[docs] def vmd_tube_tcl(self, file_name="drawTubes.tcl", radius=5): with open(file_name, 'w') as tclFile: tclFile.write("## beginning TCL script \n") def draw_tube(polymer,radius_value=10, color="cyan", resolution=5): tclFile.write("## Tube being drawn... \n") contours = np.linspace(0,1, max(2,1+polymer.num_monomers//resolution) ) rs = [polymer.contour_to_position(c) for c in contours] radius_value = str(radius_value) tclFile.write("graphics top color {} \n".format(str(color))) for i in range(len(rs)-2): r0 = rs[i] r1 = rs[i+1] filled = "yes" if i in (0,len(rs)-2) else "no" tclFile.write("graphics top cylinder {{ {} {} {} }} {{ {} {} {} }} radius {} resolution 30 filled {} \n".format(r0[0], r0[1], r0[2], r1[0], r1[1], r1[2], str(radius_value), filled)) tclFile.write("graphics top sphere {{ {} {} {} }} radius {} resolution 30\n".format(r1[0], r1[1], r1[2], str(radius_value))) r0 = rs[-2] r0 = rs[-1] tclFile.write("graphics top cylinder {{ {} {} {} }} {{ {} {} {} }} radius {} resolution 30 filled yes \n".format(r0[0], r0[1], r0[2], r1[0], r1[1], r1[2], str(radius_value))) ## material tclFile.write("graphics top materials on \n") tclFile.write("graphics top material AOEdgy \n") ## iterate through the model polymers for s in self.polymers: draw_tube(s,radius,"cyan")
[docs] def vmd_cylinder_tcl(self, file_name="drawCylinders.tcl", radius=5): with open(file_name, 'w') as tclFile: tclFile.write("## beginning TCL script \n") def draw_cylinder(polymer,radius_value=10,color="cyan"): tclFile.write("## cylinder being drawn... \n") r0 = polymer.contour_to_position(0) r1 = polymer.contour_to_position(1) radius_value = str(radius_value) color = str(color) tclFile.write("graphics top color {} \n".format(color)) tclFile.write("graphics top cylinder {{ {} {} {} }} {{ {} {} {} }} radius {} resolution 30 filled yes \n".format(r0[0], r0[1], r0[2], r1[0], r1[1], r1[2], radius_value)) ## material tclFile.write("graphics top materials on \n") tclFile.write("graphics top material AOEdgy \n") ## iterate through the model polymers for s in self.polymers: draw_cylinder(s,radius,"cyan")
[docs] class PolymerBeads(Group, metaclass=ABCMeta): """ Abstract class for bead-based models that are built from Polymer objects """ def __init__(self, polymer, sequence=None, monomers_per_bead_group = 1, rest_length = None, **kwargs): self.polymer = polymer self.sequence = sequence # self.spring_constant = spring_constant self.monomers_per_bead_group = monomers_per_bead_group if rest_length is None: rest_length = polymer.monomer_length * self.monomers_per_bead_group self.rest_length = rest_length for prop in ('segname','chain'): if prop not in kwargs: try: self.__dict__[prop] = polymer.__dict__[prop] except: pass Group.__init__(self, **kwargs) @property def monomers_per_bead_group(self): return self.__monomers_per_bead_group @property def num_bead_groups(self): return self.__num_bead_groups @monomers_per_bead_group.setter def monomers_per_bead_group(self,value): if value <= 0: raise ValueError("monomers_per_bead_group must be positive") self.num_bead_groups = int(np.ceil(self.polymer.num_monomers/value)) @num_bead_groups.setter def num_bead_groups(self,value): if value <= 0: raise ValueError("num_bead_groups must be positive") self.__num_bead_groups = value self.__monomers_per_bead_group = self.polymer.num_monomers/value # get exact ratio def _clear_beads(self): ## self.children = [] self.clear_all() @abstractmethod def _generate_ith_bead_group(self, index, position, orientation): """ Normally a bead, but sometimes a group of beads """ bead = PointParticle(type_, position, resid = index+1) pass @abstractmethod def _join_adjacent_bead_groups(self, ids): if len(ids) == 2: b1,b2 = [self.children[i] for i in ids] """ units "10 kJ/N_A" kcal_mol """ bond = HarmonicBond(k = self.spring_constant, r0 = self.rest_length, rRange = (0,500), resolution = 0.01, maxForce = 50) self.add_bond( i=b1, j=b2, bond = bond, exclude=True ) else: pass def _generate_beads(self): for i in range(self.num_bead_groups): c = self.polymer.monomer_index_to_contour(i * self.monomers_per_bead_group + (self.monomers_per_bead_group-1)/2) r = self.polymer.contour_to_position(c) o = self.polymer.contour_to_orientation(c) obj = self._generate_ith_bead_group(i, r, o) self.add(obj) for di in range(1,4): for i in range(len(self.children)-di): ids = tuple(i+_di for _di in range(di+1)) self._join_adjacent_bead_groups(ids)
[docs] class PolymerModel(ArbdModel, metaclass=ABCMeta): def __init__(self, polymers, sequences = None, monomers_per_bead_group = 1, **kwargs): """ Assign sequences """ if sequences is None: sequences = [None for i in range(len(polymers))] self.polymer_group = PolymerGroup(polymers) self.sequences = sequences self.monomers_per_bead_group = monomers_per_bead_group ArbdModel.__init__(self, [], **kwargs) """ Generate beads """ self.generate_beads()
[docs] def update_splines(self, coords): i = 0 for p,c in zip(self.polymer_group.polymers,self.children): n = len(c.children) p.set_splines(np.linspace(0,1,n), coords[i:i+n]) i += n
@abstractmethod def _generate_polymer_beads(self, polymer, sequence, polymer_index = None): pass
[docs] def generate_beads(self): self.clear_all() self.peptides = [self._generate_polymer_beads(p,s,i) for i,(p,s) in enumerate(zip(self.polymer_group.polymers,self.sequences))] for s in self.peptides: self.add(s) s._generate_beads()