import copy
from hoobas.Units import SimulationUnits as SU
import numpy as np
from hoobas.CGBead import Bead as Bead
import warnings
from hoobas import Util
from hoobas.Quaternion import Quat
import inspect
import random
import itertools
import pickle
try:
NETWORKX_FOUND = True
import networkx
except ImportError:
NETWORKX_FOUND = False
print('Networkx package not found, graphs and backbone angle / dihedral building will be disabled')
[docs]class TopologyType(object):
"""
Base bonded force-field topology type which defines the base variables of every bonded type.
The class implements two dictionaries, one for topology constants (energy, distances, etc.) and one for the units of
these constants. a name must also be given to the topology type to distinguish it from other types.
Args:
name: name of the type
topodict: dictionary which defines topology constants
unitdict: dictionary of units of the topology constants with matching keys of topodict
Examples:
This yields a bondname with name of '' and keys 'energy_constant':0.0 (E/LL) and 'distance_constant': 0.0 (L)
>>> bond_type = BondType()
Topology types support initialization by list for harmonic bonds. For instance, this yields a
bondtype named 'BondTypeName', with energy_constant of 1.0 (E/LL) and distance_constant of 5.0 (L)
>>> bond_type = BondType(['BondTypeName', 1.0, 5.0])
This yields a bondtype name 'Anharmonic with 'K2' of 1.0 (E/LL), 'K4' of 1.0 (E/LLLL) and 'D' of 1.0 (L)
>>> bond_type = BondType('Anharmonic', {'K2':1.0, 'K4':1.0, 'D':1.0}, {'K2':'E/LL', 'K4':'E/LLLL', 'D':'L'})
Attributes:
typename: name of the TopologyType ('BondTypeName' or 'Anharmonic' in previous examples)
topology_constants: dictionary with numerical properties of the topology
topology_units: dictionary with unit transformation of the topology
class_name: class of potential associated with TopologyType, reserved for future implementation
"""
def __init__(self, name='', topodict=None, unitdict=None):
self.typename = name
self.class_name = None # class of potential, say harmonic or cosine square
self.topology_constants = {}
self.topology_units = {}
if issubclass(topodict.__class__, dict):
self.topology_constants.update(topodict)
if issubclass(unitdict.__class__, dict):
self.topology_units.update(unitdict)
for (key, val) in self.topology_units.items():
if not issubclass(val.__class__, str):
self.topology_units[key] = str(val) # attempt to stringify the value
def __str__(self):
return self.typename
def __eq__(self, other):
return str(self) == str(other)
def __ne__(self, other):
return str(self) != str(other)
def __getitem__(self, item):
if item in self.topology_constants:
return self.topology_constants[item]
else:
return None
@property
def potential(self):
"""
Potential name associated with bonded parameter types.
:getter: returns the potential name
:setter: sets the potential name
:type: str
"""
return self.class_name
@potential.setter
def potential(self, name):
self.class_name = str.lower(name)
def __nonzero__(self): # override nonzero so we can directly check if bondtype to see if it has a force-field
if len(self.topology_constants) == 0:
return False
z = True
for (k, v) in self.topology_constants.items():
z = z and v == 0.0
return z
[docs] @staticmethod
def distance(obj1, obj2):
"""
Computes a topological distance between force-field parameters, defined as :math:`max_v(v_1 / v_2 - 1.0)`, where
:math:`v` iterates over force-field keys. If named keys of the force-field differ, the distance is set to
:math:`\infty`
:param obj1: First type to compare
:type obj1: TopologyType
:param obj2: Second type to compare
:type obj2: TopologyType
:return: None
"""
for (k, v) in obj1.topology_constants.items():
if k not in obj2.topology_constants:
return float('inf')
for (k, v) in obj2.topology_constants.items():
if k not in obj1.topology_constants:
return float('inf')
mdist = 0.0
for (k, v) in obj1.topology_constants.items():
v2 = obj2.topology_constants[k]
try:
mdist = max(mdist, v/v2 - 1.0, v2/v - 1.0)
if v * v2 < 0:
mdist = float('inf')
except ZeroDivisionError:
if v == v2 == 0.0:
mdist = max(mdist, 0.0)
else:
mdist = float('inf')
return mdist
[docs]class BondType(TopologyType):
"""
A default class for bond types. Has topology properties of energy_constant (E/LL) and distance_constant(L).
"""
def __init__(self, name=None, topodict=None, unitdict=None):
super(BondType, self).__init__(name, topodict, unitdict)
if issubclass(name.__class__, list):
topodict = name
self.name = ''
if name is None:
self.name = ''
if issubclass(topodict.__class__, list): # unparsed by main constructor
self.typename = topodict[0]
self.topology_constants['energy_constant'] = topodict[1]
self.topology_constants['distance_constant'] = topodict[2]
if topodict is None:
self.topology_constants['energy_constant'] = 0.0
self.topology_constants['distance_constant'] = 0.0
if unitdict is None or issubclass(topodict.__class__, list):
self.topology_units['energy_constant'] = 'E/LL'
self.topology_units['distance_constant'] = 'L'
for key in self.topology_constants:
if key not in self.topology_units:
raise ValueError('Dictionary for topology constant does not match dictionary for units')
[docs]class AngleType(TopologyType):
"""
A default class for anglular potentials. Has topology constants of energy_constants (E) and angle_constant (1)
"""
def __init__(self, name='', topodict=None, unitdict=None):
super(AngleType, self).__init__(name, topodict, unitdict)
if issubclass(name.__class__, list):
topodict = name
self.name = ''
if name is None:
self.name = ''
if topodict is None:
self.topology_constants['energy_constant'] = 0.0
self.topology_constants['angle_constant'] = 0.0
if issubclass(topodict.__class__, list):
self.typename = topodict[0]
self.topology_constants['energy_constant'] = topodict[1]
self.topology_constants['angle_constant'] = topodict[2]
if unitdict is None or issubclass(topodict.__class__, list):
self.topology_units['energy_constant'] = 'E'
self.topology_units['angle_constant'] = '1'
for key in self.topology_constants:
if key not in self.topology_units:
raise ValueError('Dictionary for topology constant does not match dictionary for units')
[docs]class DihedralType(TopologyType):
"""
A default class for dihedral types. Has topology constants of energy_constant (E) and angle_constant (1)
"""
def __init__(self, name='', topodict=None, unitdict=None):
super(DihedralType, self).__init__(name, topodict, unitdict)
if issubclass(name.__class__, list):
topodict = name
self.name = ''
if name is None:
self.name = ''
if topodict is None:
self.topology_constants['energy_constant'] = 0.0
self.topology_constants['angle_constant'] = 0.0
if issubclass(topodict.__class__, list):
self.typename = topodict[0]
self.topology_constants['energy_constant'] = topodict[1]
self.topology_constants['angle_constant'] = topodict[2]
if unitdict is None or issubclass(topodict.__class__, list):
self.topology_units['energy_constant'] = 'E'
self.topology_units['angle_constant'] = '1'
for key in self.topology_constants:
if key not in self.topology_units:
raise ValueError('Dictionary for topology constant does not match dictionary for units')
[docs]class ImproperType(DihedralType):
"""
A default class for improper types. Has topology constants of energy_constant (E) and angle_constant (1)
"""
def __init__(self, name='', topodict=None, unitdict=None):
super(ImproperType, self).__init__(name, topodict, unitdict)
if issubclass(name.__class__, list):
topodict = name
self.name = ''
if name is None:
self.name = ''
if topodict is None:
self.topology_constants['energy_constant'] = 0.0
self.topology_constants['angle_constant'] = 0.0
if issubclass(topodict.__class__, list):
self.typename = topodict[0]
self.topology_constants['energy_constant'] = topodict[1]
self.topology_constants['angle_constant'] = topodict[2]
if unitdict is None or issubclass(topodict.__class__, list):
self.topology_units['energy_constant'] = 'E'
self.topology_units['angle_constant'] = '1'
for key in self.topology_constants:
if key not in self.topology_units:
raise ValueError('Dictionary for topology constant does not match dictionary for units')
[docs]class TopologyTypeList(object):
"""
List of topology types (bonded force-field) objects used by CompositeObject
This class behaves similar to dictionary, topology types can be queried by key and can be iterated over.
Additionally, the object supports augmented addition to join multiple lists of force-fields
Example:
>>> bond_type_list = TopologyTypeList(BondType)
>>> a_bond_type = BondType('BondType', [1.0, 1.0])
>>> bond_type_list += a_bond_type
>>> another_type_list = TopologyTypeList(BondType)
>>> bond_type_list += another_type_list
TopologyTypeLists are normally linked to a TopologyList to synchronize remapping and addition of force-fields. Types
can be remapped with the remap method.
Example:
>>> bond_type_list = TopologyTypeList(BondType)
>>> bond_type_list += BondType('NoParams') # a bond type with no parameters
>>> bond_type_list['NoParams'] = BondType('Params', [1.0, 1.0]) # remaps NoParams->Params and sets params
>>> bond_type_list.remap('Params', 'NoParams') # remaps Params->NoParams
>>> bond_type_list['NoParams'] = BondType() # removes coefficients from NoParams
Reduction in number of types can be achived by removing force-fields that are close in parameters, see the make_set
method.
Args:
TopologyTypeClass: class of the held TopologyType
"""
def __init__(self, TopologyTypeClass):
if not issubclass(TopologyTypeClass, TopologyType):
raise SyntaxError('Trying to create a topology type list of non-topology type')
# internals
self.topology_type = TopologyTypeClass
self.topology_list = {}
self.topology_list_link = None
def __iadd__(self, other):
"""
Joins either a list of types of topologies (BondList1 += Bondlist2) or adds a topology type Bondlist += bondtype
:param other: other object to add
:type other: list, list[TopologyType], TopologyType
:return: self
"""
if issubclass(other.__class__, list): # list initialization or list construction
if len(other) == 0: # when an automatic aggregation of types leads to an empty list
return self
# check whether other is of the form ['typename', constant, constant] or [['typename1', c1, c2], ['type'...
if next(other.__iter__()).__class__ != self.topology_type and next(other.__iter__()).__class__ != list:
self += self.topology_type(other) # coerce whatever is in there to a TopologyType
return self
for item in other:
if item.__class__ == self.topology_type: # item in list is already a TopologyType
self += item
else:
self += self.topology_type(item)
return self
if other.__class__ == self.topology_type:
self.topology_list.update({str(other): other})
return self
if self.__class__ == other.__class__:
for topo in other:
self += topo
return self
raise SyntaxError('Trying to add a topology of type: ' + str(other.__class__)
+ ' to a list of type: ' + str(self.topology_type))
def __contains__(self, TopologyType):
return str(TopologyType) in self.topology_list
def __iter__(self):
for (key, topo) in self.topology_list.items():
yield topo
def __getitem__(self, item):
if isinstance(item, str) and item in self.topology_list:
return self.topology_list[item]
elif item.__class__ == self.topology_type:
return self.topology_list[str(item)]
return None
def __setitem__(self, key, value):
if value.__class__ != self.topology_type:
raise SyntaxError('Trying to add a non matching topology type: '
+ str(value.__class__) + ' and :' + str(self.topology_type))
if key in self and str(value) != '':
self.remap(key, str(value))
key = str(value)
if key != str(value): # check whether someone forgot to set the name inside the topology type
value.typename = key
self.topology_list.update({key: value})
[docs] def link(self, link):
"""
Links the topology list to a TopologyTypeList, this is normally done automatically by CompositeObject
:param link: TypeList instance
:type link: TopologyTypeList
:return: None
"""
self.topology_list_link = link
[docs] def unlink(self):
"""
Unlinks the topology type list from the current topology list
:return: None
"""
self.topology_list_link = None
[docs] def make_set(self, tolerance=1e-3, avoid=None): # make an approximate set based on numerical tolerance
"""
Makes a reduced set of topology types by merging types with distance less than the tolerance
:param tolerance: maximum distance between merged types
:type tolerance: float
:param avoid: list of topology types to avoid merging
:type avoid: list
:return:
"""
_copy = copy.deepcopy(self.topology_list)
for (k, v) in _copy.items():
if (avoid is not None and k in avoid) or k not in self.topology_list:
continue
for (k2, v2) in _copy.items():
if (avoid is not None and k2 in avoid) or k2 not in self.topology_list or v2 == v:
continue
d = TopologyType.distance(v, v2)
if d < tolerance:
if self.topology_list_link is not None:
self.topology_list_link.remap(k2, k)
del self.topology_list[k2]
[docs] def remap(self, old, new):
"""
Remaps a topology name to another. If the new name already exists, the old type is destroyed and topologies of
that type remapped to new.
:param old: topology type to remap
:type old: str
:param new: type to map to
:type new: str
:return: None
"""
if new in self and old in self:
warnings.warn('Remap called on two existing TopologyTypes, ' + str(old) +
' will be destroyed; this is normal behavior for reducing topologies types', UserWarning)
self.topology_list.pop(old)
if self.topology_list_link is not None:
self.topology_list_link.remap(old, new)
return
if old not in self:
return
if self.topology_list_link is not None:
self.topology_list_link.remap(old, new)
self[new] = self.topology_list.pop(old)
[docs]class TopologyItem(object):
"""
Defines a topology item, for instance, a bond
Args:
idict: optional dict initializer for the topologies
Note: Subclasses have well defined constructors. Bonds have 2 members and can be list initialized
>>> abond = Bond(['BondType', 0, 1])
>>> another_bond = Bond({'topology_tuple': (1, 2), 'topology_name': 'BondType'})
Attributes:
topology_tuple: indicates the associated particles in the topology. For a bond this would be a 2 membered tuple
topology_name: indicates the name associated with the topology. This is normally linked to a TopologyType name
"""
def __init__(self, idict=None):
self.topology_tuple = None
self.topology_name = ''
if issubclass(idict.__class__, dict):
if 'topology_tuple' in idict:
self.topology_tuple = idict['topology_tuple']
if not issubclass(self.topology_tuple.__class__, tuple):
self.topology_tuple = tuple(self.topology_tuple)
if 'topology_name' in idict:
self.topology_name = idict['topology_name']
if not issubclass(self.topology_name.__class__, str):
self.topology_name = str(self.topology_name)
elif issubclass(idict.__class__, tuple):
self.topology_tuple = idict
def __str__(self):
return self.topology_name
def __getitem__(self, item):
if item == 0:
return self.topology_name
if isinstance(item, int) and len(self.topology_tuple) >= item:
return self.topology_tuple[item - 1]
return None
def __setitem__(self, key, value):
if key == 0:
self.topology_name = str(value)
if isinstance(key, int) and len(self.topology_tuple) >= key:
tlist = [i for i in self.topology_tuple]
tlist[key - 1] = value
self.topology_tuple = tuple(tlist)
def shift_index(self, value):
self.topology_tuple = tuple([i + value for i in self.topology_tuple])
[docs]class Bond(TopologyItem):
def __init__(self, ilist=None):
if issubclass(ilist.__class__, tuple) and len(ilist) != 2:
raise SyntaxError('Trying to create bond that does not have 2 members')
super(Bond, self).__init__(ilist)
if issubclass(ilist.__class__, list):
self.topology_name = ilist[0]
self.topology_tuple = (ilist[1], ilist[2])
[docs]class Angle(TopologyItem):
def __init__(self, ilist=None):
if issubclass(ilist.__class__, tuple) and len(ilist) != 3:
raise SyntaxError('Trying to create angle that does not have 3 members')
super(Angle, self).__init__(ilist)
if issubclass(ilist.__class__, list):
self.topology_name = ilist[0]
self.topology_tuple = (ilist[1], ilist[2], ilist[3])
[docs]class Dihedral(TopologyItem):
def __init__(self, ilist=None):
if issubclass(ilist.__class__, tuple) and len(ilist) != 4:
raise SyntaxError('Trying to create dihedral that does not have 4 members')
super(Dihedral, self).__init__(ilist)
if issubclass(ilist.__class__, list):
self.topology_name = ilist[0]
self.topology_tuple = (ilist[1], ilist[2], ilist[3], ilist[4])
[docs]class Constraint(TopologyItem):
def __init__(self, ilist, distance=0.0): # do not perform checks, constraints could be weird stuff
super(Constraint, self).__init__(ilist)
self.distance = distance
if issubclass(ilist.__class__, list):
self.topology_tuple = (ilist[0], ilist[1])
[docs]class TopologyList(object):
"""
Defines a list of topologies of a given class of topologies.
A TopologyList keeps tracks of all bonded interactions of a given class (bond, angle, dihedral, ...). It essentially
works like a list: topologies can be retrieved by index with [] and iterated over.
Args:
topology_class: this defines the type of topology in the list (to avoid mixing bonds with angles)
Attributes:
topologies: list of topologies of class TopologyClass
topology_cls: class of topologies contained
topology_type_list: link to a list of types of topologies
Note: Adding a new topology (say a bond of type 'A-A') will check if 'A-A' is in topology_type_list and will add it
as a default topology of its type (calling topology_cls constructor) if it is not present. The topology list supports
augmented addition:
>>> BondList = TopologyList(Bond)
>>> ABond = Bond(['ABondName', 0, 1])
>>> BondList += ABond
>>> AnotherBondList = TopologyList(Bond)
>>> BondList += AnotherBondList
"""
def __init__(self, topology_class):
self.topologies = []
self.topology_cls = topology_class
self.topology_type_list = None
if (issubclass(self.topology_cls, Bond) or issubclass(self.topology_cls, Constraint)) and NETWORKX_FOUND:
self.graph = networkx.Graph()
else:
self.graph = None
def __len__(self):
return len(self.topologies)
def __iadd__(self, other):
if issubclass(other.__class__, self.topology_cls): # append a bond to the list
self.topologies.append(other)
if self.topology_type_list is not None:
if other.topology_name not in self.topology_type_list:
ttype = self.topology_type_list.topology_type(name=other.topology_name)
self.topology_type_list += ttype
if self.graph is not None:
self.graph.add_edge(*other.topology_tuple, bond_type=self.topologies[-1].topology_name)
return self
if issubclass(other.__class__, TopologyList) and other.topology_cls == self.topology_cls:
for item in other:
self += item
return self
if issubclass(other.__class__, list): # list initialization
self += self.topology_cls(other)
return self
raise SyntaxError('Trying to add incompatible topology lists')
def shift_index(self, value):
for item in self.topologies:
item.shift_index(value)
if self.graph is not None:
self.graph = networkx.relabel_nodes(self.graph, lambda x: x + value)
def link(self, type_list): # if a link is set, every time a topology is appended, it will append the type too
self.topology_type_list = type_list
def unlink(self):
self.topology_type_list = None
def __iter__(self):
for item in self.topologies:
yield item
def __getitem__(self, item):
return self.topologies[item]
def __contains__(self, item):
names = [str(i) for i in self.topologies]
return item in names
def __isub__(self, other):
if other.__class__ == self.topology_cls.__class__:
index = self.topologies.index(other)
if self.graph is not None:
self.graph.remove_edge(other.topology_tuple[0], other.topology_tuple[1])
self.topologies.pop(index)
return self
elif issubclass(other.__class__, list):
for item in other:
self -= item
return self
def wipe(self):
self.topologies = []
if self.graph is not None:
self.graph = networkx.Graph()
def get_topology(self, topology_tuple, edge=None):
edge = {} if edge is None else edge
topology_name = None if len(edge) == 0 else edge[edge.keys()[0]]
if topology_name is not None:
return self.topology_type_list[topology_name]
for topo in self.topologies:
if set(topology_tuple) == set(topo.topology_tuple) and (topology_name is None or topology_name == topo.topology_name):
return topo
def remap(self, original, final):
for topo in self.topologies:
if str(topo) == original:
topo.topology_name = final
def transform(self, length):
if self.topology_cls != Constraint:
pass
for topo in self.topologies:
topo.distance *= length
[docs]class BeadList(object):
"""
Object to manage and link the list of beads and a numpy array;
Attributes:
_beads: list of bead objects
_positions: N x 3 numpy array of all the bead positions
_references: dictionary of objects which positions are linked with, but are also part of another BeadList
graph: a directed networkx graph that keeps track of which subunit beads belong to
Note: the bead positions (bead[i].position) has the same reference as the position array (self._position[i,:])
Methods:
relink(): resets all the bead position references to the _position array
change_units(new): change all the units of this object to the SimulationUnits new
__iter__(): yields the beads in the beadlist
__getitem__(index): returns the bead at index
__len__: returns the length of the bead array
Properties:
positions(): returns _positions
positions(arr): sets the contents of _positions to arr, but keeps the reference
beads(): returns the list of beads in self._beads
"""
def __init__(self):
self._beads = []
self._positions = np.zeros((0, 3), dtype=np.float)
self._references = []
self.graph = None
self.total_graph = None
def __iadd__(self, other):
ilen = len(self._beads)
if other.__class__ == Bead:
self._beads.append(other)
self._positions = np.append(self._positions, np.array([other.position]), axis=0)
if NETWORKX_FOUND:
root = next(networkx.topological_sort(self.graph))
self.graph.add_edge(root, other)
self.relink()
return self
elif other.__class__ == self.__class__:
self._beads += other.beads
self._positions = np.append(self._positions, other.positions, axis=0)
# link lists together
self._references.append({'object': other,
'start_index': ilen,
'stop_index': ilen + len(other)})
self.relink()
return self
elif issubclass(other.__class__, list):
if NETWORKX_FOUND:
root = next(networkx.topological_sort(self.graph))
tuple_list = []
for item in other:
self._beads.append(item)
self._positions = np.append(self._positions, np.array([item.position]), axis=0)
if NETWORKX_FOUND:
tuple_list.append((root, item))
self.graph.add_edges_from(tuple_list)
self.relink()
return self
def relink(self): # this ensures that the value in bead.position is a slice-reference to self._positions
for ref in self._references:
ref['object']._positions = self._positions[ref['start_index']:ref['stop_index'], :]
ref['object'].relink()
for ibead in range(len(self._beads)):
self._beads[ibead]._position = self._positions[ibead, :]
def __iter__(self):
self.relink()
for bead in self._beads:
yield bead
def __getitem__(self, item):
return self._beads[item]
def __len__(self):
return self.beads.__len__()
def change_units(self, length, energy, mass):
for bead in self._beads:
bead.change_units(length, energy, mass)
@property
def positions(self):
return self._positions
@positions.setter
def positions(self, array):
if array.__class__ != self._positions.__class__ or array.size != self._positions.size or array.ndim != self._positions.ndim:
raise SyntaxError('Position descriptor does not match N x 3 numpy array')
self._positions[:] = array[:]
@property
def beads(self):
return self._beads
def reset(self):
self._beads = []
self._positions = np.array((0, 3), dtype=np.float)
def set_attribute_by_type(self, beadtype, attribute, value):
for bead in self._beads:
if bead.beadtype == beadtype:
setattr(bead, attribute, value)
def remap(self, old, new):
self.set_attribute_by_type(old, 'beadtype', new)
def index(self, bead):
return self.beads.index(bead)
def root_index(self, bead):
root = networkx.topological_sort(self.graph).next()
return root.beads.index(bead)
[docs]class AttachmentSite(dict):
def __init__(self, index, qualifier=None, bound_to=None, orientation=None, restriction=None):
super(AttachmentSite, self).__init__()
self.update({'index': index,
'qualifier': qualifier,
'bound_to': bound_to,
'orientation': orientation,
'restriction': restriction})
def rotate(self, op):
q = Quat(op)
tr = q.transform
if self['orientation'] is not None:
self['orientation'] = tr.dot(self['orientation'])
def __eq__(self, other):
return id(self) == id(other)
class AttachmentSiteList(object):
def __init__(self):
self.free_sites = []
self.bound_sites = []
self.link_to_parent = None
def __iadd__(self, other):
if issubclass(other.__class__, AttachmentSite):
self.free_sites.append(other)
return self
elif issubclass(other.__class__, AttachmentSiteList):
self.free_sites += other.free_sites
self.bound_sites += other.bound_sites
return self
elif issubclass(other.__class__, dict):
self.free_sites.append(AttachmentSite(**other))
return self
elif hasattr(other, '__iter__'):
for obj in other:
self += obj
return self
else:
self += AttachmentSite(other)
return self
def index(self, site):
return self.link_to_parent.p_num.index(site['index'])
def pop(self, index):
return self.free_sites.pop(index)
def pop_site(self, site):
return self.free_sites.pop(self.free_sites.index(site))
def pop_internal_index(self, iidx):
for site in self.free_sites:
if site['index'] == iidx:
return self.free_sites.pop(self.free_sites.index(site))
return IndexError('pop index not in free sites')
def bind_site(self, site, bound_to=None):
site['bound_to'] = bound_to
self.bound_sites.append(self.free_sites.pop(self.free_sites.index(site)))
def shift_indexes(self, offset):
for key in self.free_sites:
key['index'] += offset
for key in self.bound_sites:
key['index'] += offset
def get_qualified(self, qualifier=None):
list_of_qualified = []
for site in self.free_sites:
# parse out qualifications of sites, either None, some object which can be compared or somekind of function
# of beads
if qualifier is None:
list_of_qualified.append(site)
continue
if site['qualifier'] == qualifier:
list_of_qualified.append(site)
continue
return list_of_qualified
def __len__(self):
return len(self.free_sites)
def count_index(self, index):
count = 0
for item in self.bound_sites:
if item['index'] == index:
count += 1
for item in self.free_sites:
if item['index'] == index:
count += 1
return count
def __iter__(self):
it = itertools.chain(self.free_sites, self.bound_sites)
for item in it:
yield item
def rotate(self, op):
for item in self:
item.rotate(op)
def reset_free_sites(self):
self.free_sites = []
@staticmethod
def get_qualified_pair(alist1, aqualifier, alist2, bqualifier):
# get the list of sites from alist1 and alist2 that match their respective qualifiers
qualified_list_a = alist1.get_qualified(aqualifier)
qualified_list_b = alist2.get_qualified(bqualifier)
# if either list is empty, return an empty list
if len(qualified_list_a) == 0 or len(qualified_list_b) == 0:
return []
# build the list of available pairs
pair_list = []
for siteA in qualified_list_a:
for siteB in qualified_list_b:
if siteA['restriction'] is None and siteB['restriction'] is None:
pair_list.append((siteA, siteB))
continue
if siteA['restriction'] is not None and siteB['qualifier'] is None:
continue
if siteB['restriction'] is not None and siteA['qualifier'] is None:
continue
if siteA['restriction'] is not None and siteA['restriction'] == siteB['qualifier'] and siteB['restriction'] is None:
pair_list.append((siteA, siteB))
continue
if siteB['restriction'] is not None and siteB['restriction'] == siteA['qualifier'] and siteA['restriction'] is None:
pair_list.append((siteA, siteB))
continue
if siteA['restriction'] == siteB['qualifier'] and siteB['restriction'] == siteA['qualifier']:
pair_list.append((siteA, siteB))
continue
return pair_list
[docs]def deferred(func):
"""
Decorator to defer function calls until the object is built
:param func: decorated method
:type func: callable
:return: function wrapper
:rtype: callable
"""
def wrapper(*args, **kwargs):
self = args[0]
if self.is_built:
return func(*args, **kwargs)
else:
self.deferred.append({'call': func, 'args': args, 'kwargs': kwargs})
return wrapper
[docs]class CompositeObject(object):
"""
Class that handles topology information of building blocks of hoobas and keeps track of force-field constants
Args:
units: description of the current simulation units, as given by a SimulationUnits object. Defaults to default units
Attributes:
p_num: particle tag index of the CompositeObject used for building
beads: list of beads, managed by the BeadList class
bonds, angles, dihedrals, impropers: list of topologies, Managed by a TopologyList
subunits: Other CompositeObjects that make up the current one
superunit: CompositeObject where this is a subunit
att_list: list of sites where this has been attached
constraints: list of distance constraints in the simulation. Managed by a TopologyList
b_types, a_types, d_types, i_types: list of force-fields related to the bonds. Managed by TopologyTypeList
Properties:
positions: returns a N x 3 numpy array from the BeadList class. The getter returns a reference so that it can
be modified without a set, the setter will not change the underlying reference
"""
_CONSTRAINT_INDEX = 0 # global index needed to translate constraints to bonds
def __init__(self, units=None, **kwargs):
# internal bead representation
self.p_num = []
self.beads = BeadList()
# actual topology
self.bonds = TopologyList(Bond)
self.angles = TopologyList(Angle)
self.dihedrals = TopologyList(Dihedral)
self.impropers = TopologyList(Dihedral)
# subconstituants that have been included, in case a method needs to be called on them
self.subunits = []
self.subunits_graph = None
self._top_level_graph = None
if NETWORKX_FOUND:
self.subunits_graph = networkx.DiGraph()
self.subunits_graph.add_node(self)
self.beads.graph = self.subunits_graph
self._top_level_graph = networkx.DiGraph()
self.beads.total_graph = self._top_level_graph
self.superunit = None
# a list of references to rigid bodies that either belong to this CompositeObject or subunits in the object
self.rigid_body_references = []
self._body = -1
# distance constraints
self.constraints = TopologyList(Constraint)
# force-field information; link it with topology
self.b_types = TopologyTypeList(BondType)
self.bonds.link(self.b_types)
self.b_types.link(self.bonds)
self.a_types = TopologyTypeList(AngleType)
self.angles.link(self.a_types)
self.a_types.link(self.angles)
self.d_types = TopologyTypeList(DihedralType)
self.dihedrals.link(self.d_types)
self.d_types.link(self.dihedrals)
self.i_types = TopologyTypeList(ImproperType)
self.impropers.link(self.i_types)
self.i_types.link(self.impropers)
# attachment topologies
self.att_list = AttachmentSiteList()
self.att_list.link_to_parent = self
# units
if units is None:
self.units = SU()
else:
self.units = copy.deepcopy(units)
# calls that have been deferred to after the object is actually built
self.deferred = []
self.is_built = False
def build(self):
self.is_built = True
# deferred function caller
def build_finalize(self):
for deferred_call in self.deferred:
deferred_call['call'](*deferred_call['args'], **deferred_call['kwargs'])
@property
def top_level_graph(self):
return self._top_level_graph
@top_level_graph.setter
def top_level_graph(self, graph):
self._top_level_graph = graph
self.beads.total_graph = graph
@property
def body(self):
return self._body
@property
def remaining_free_sites(self):
return self.att_list.free_sites
@property
def bound_sites(self):
return self.att_list.bound_sites
@property
def pnum_offset(self):
try:
return self.p_num[0]
except IndexError:
return 0
@pnum_offset.setter
def pnum_offset(self, val):
off = val - self.pnum_offset
for i in range(self.p_num.__len__()):
self.p_num[i] += off
# shift the topology
self.bonds.shift_index(off)
self.angles.shift_index(off)
self.dihedrals.shift_index(off)
self.impropers.shift_index(off)
self.constraints.shift_index(off)
self.att_list.shift_indexes(off)
[docs] def merge(self, other): # merge other in self without modifications to other
"""
merge another composite into this one without changing the other object
:param other: object to merge into this one
:type other: CompositeObject
:return: None
"""
if not other.is_built:
other.build()
other.build_finalize()
if not issubclass(other.__class__, CompositeObject):
raise SyntaxError('Trying to merge a non-CompositeObject into a CompositeObject')
other.change_units(self.units)
index_set = self.p_num[-1] + 1 if len(self.p_num) > 0 else 0
index_shift = index_set - other.pnum_offset
for ref in other.rigid_body_references:
for i in range(len(ref.p_num)):
ref.p_num[i] += index_shift
other.pnum_offset = index_set # TODO: replace pnum with indexes directly from builder.beads by calling top beadlist._bead.index
self.p_num += other.p_num
self.beads += other.beads
# merge topology types first
self.b_types += other.bond_types
self.a_types += other.angle_types
self.d_types += other.dihedral_types
self.i_types += other.improper_types
# merge topology lists
self.bonds += other.bonds
self.angles += other.angles
self.dihedrals += other.dihedrals
self.impropers += other.impropers
self.constraints += other.constraints
# add the merged object into the subunits, grab the rigid body references from it
self.subunits.append(other)
if NETWORKX_FOUND:
self.subunits_graph.add_edges_from(other.subunits_graph.edges)
self.subunits_graph.add_nodes_from(other.subunits_graph.nodes)
self.subunits_graph.add_edge(self, other)
for ref in other.rigid_body_references:
self.rigid_body_references.append(ref)
if other.body != -1:
self.rigid_body_references.append(other)
other.superunit = self
# properties to read types
@property
def bond_types(self):
return self.b_types
@bond_types.setter
def bond_types(self, val):
if self.b_types.topology_list_link == val.topology_list_link and self.bonds.topology_type_list == val:
self.b_types = val
else:
self.bonds.unlink()
self.b_types = val
self.bonds.link(self.b_types)
self.b_types.link(self.bonds)
@property
def angle_types(self):
return self.a_types
@angle_types.setter
def angle_types(self, val):
if self.a_types.topology_list_link == val.topology_list_link and self.angles.topology_type_list == val:
self.a_types = val
else:
self.angles.unlink()
self.a_types = val
self.angles.link(self.a_types)
self.a_types.link(self.angles)
@property
def dihedral_types(self):
return self.d_types
@dihedral_types.setter
def dihedral_types(self, val):
if self.d_types.topology_list_link == val.topology_list_link and self.dihedrals.topology_type_list == val:
self.d_types = val
else:
self.dihedrals.unlink()
self.d_types = val
self.dihedrals.link(self.d_types)
self.d_types.link(self.dihedrals)
@property
def improper_types(self):
return self.i_types
@improper_types.setter
def improper_types(self, val):
if self.i_types.topology_list_link == val.topology_list_link and self.dihedrals.topology_type_list == val:
self.i_types = val
else:
self.impropers.unlink()
self.i_types = val
self.impropers.link(self.i_types)
self.i_types.link(self.impropers)
@deferred
def make_root_topology(self):
self.subunits_graph = networkx.DiGraph()
self.subunits_graph.add_node(self)
self.beads.graph = self.subunits_graph
for bead in self.beads:
self.subunits_graph.add_edge(self, bead)
@deferred
def change_units(self, new):
if not issubclass(new.__class__, SU):
raise SyntaxError('Composite object got non-units, expected : ' + str(SU)
+ ' instead got ' + str(new.__class__))
# calculate transformation
m_len = new.get_length(1.0, self.units.lunit)
m_en = new.get_energy(1.0, self.units.Eunit)
m_mass = new.get_mass(1.0, self.units.munit)
# change topology units
self.beads.change_units(m_len, m_en, m_mass)
self.b_types.transform(m_len, m_en, m_mass)
self.a_types.transform(m_len, m_en, m_mass)
self.d_types.transform(m_len, m_en, m_mass)
self.i_types.transform(m_len, m_en, m_mass)
self.constraints.transform(m_len)
self.units = new
for subitem in self.subunits:
subitem.units = new
@deferred
def reduce_topology(self, tolerance=1e-3, bonds=False, angles=False, dihedrals=False, avoid=None):
"""
Reduces the number of types of topologies by combining types which have similar parameters
:param tolerance: maximum relative distance between topology types
:type tolerance: float
:param bonds: whether to reduce bonds (default False)
:type bonds: bool
:param angles: whether to reduce angles (default False)
:type angles: bool
:param dihedrals: whether to reduce dihedrals and impropers (default False)
:type dihedrals: bool
:param avoid: list of topology names to avoid merging
:type avoid: list(str)
:return: None
"""
if bonds:
self.b_types.make_set(tolerance, avoid)
if angles:
self.a_types.make_set(tolerance, avoid)
if dihedrals:
self.d_types.make_set(tolerance, avoid)
self.i_types.make_set(tolerance, avoid)
@property
def positions(self):
return self.beads.positions
@positions.setter
@deferred
def positions(self, array):
self.beads.positions = array
@property
def center_position(self):
return self.positions[0, :]
@center_position.setter
@deferred
def center_position(self, value):
diff = value - self.center_position
self.positions += diff
@deferred
def remap_bond(self, old, new):
self.b_types.remap(old, new)
@deferred
def remap_angle(self, old, new):
self.a_types.remap(old, new)
@deferred
def remap_dihedral(self, old, new):
self.d_types.remap(old, new)
@deferred
def remap_improper(self, old, new):
self.i_types.remap(old, new)
@deferred
def remap_beadtype(self, old, new):
self.beads.remap(old, new)
@deferred
def set_attribute_by_beadtype(self, beadtype, attribute, value):
self.beads.set_attribute_by_type(beadtype, attribute, value)
@deferred
def set_charge_by_beadtype(self, beadtype, charge):
self.set_attribute_by_beadtype(beadtype, 'charge', charge)
@deferred
def set_diameter_by_beadtype(self, beadtype, diameter):
self.set_attribute_by_beadtype(beadtype, 'diameter', diameter)
@deferred
def add_free_attachment_sites(self, key_search, max_binding_per_bead=1, properties=None):
"""
Deferred handler for adding attachment sites. The function adds attachment sites until the max_binding_per_bead
is achieved on any bead that satisfied the key_search description. Functions passed in max_binding or properties
will query index from the bead index, any other property from bead[index] first and from the object itself if
it is not found in bead[i]
:param key_search: key:value search in bead properties; also allows for 'index' which refers to bead
:type key_search: dict
index. Values may be functions which are evaluated at runtime.
:param max_binding_per_bead: maximum number of binding on a given bead. The function will add
binding sites until that value is matched.
:type max_binding_per_bead: int or callable
:param properties: determines the properties of added attachment sites; see the AttachmentSite class.
callable values passed in the dict will be called in the same fashion as max_binding_per_bead
:type properties: dict
:return: None
"""
for i in range(len(self.beads)):
if hasattr(max_binding_per_bead, '__call__'):
args = inspect.getargspec(max_binding_per_bead).args
arguments = {}
for arg in args:
if arg == 'index':
arguments[arg] = i
elif hasattr(self.beads[i], arg):
arguments[arg] = getattr(self.beads[i], arg)
else:
arguments[arg] = getattr(self, arg)
max_bindings = max_binding_per_bead(**arguments)
else:
max_bindings = max_binding_per_bead
number_of_bindings = self.att_list.count_index(self.p_num[i])
if number_of_bindings >= max_bindings:
continue
ok_to_bind = True
for key, value in key_search.items():
if key != 'index':
if not hasattr(value, '__call__'):
ok_to_bind = ok_to_bind and getattr(self.beads[i], key) == value
else:
args = inspect.getargspec(value).args
if len(args) == 0:
ok_to_bind = ok_to_bind and getattr(self.beads[i], key) == value()
else:
ok_to_bind = ok_to_bind and value(getattr(self.beads[i], key))
else: # check for functions passed in for index
if hasattr(value, '__call__'):
args = inspect.getargspec(value).args
if len(args) == 0:
ok_to_bind = ok_to_bind and value()
else:
ok_to_bind = ok_to_bind and value(i)
else:
ok_to_bind = ok_to_bind and value == i
if ok_to_bind:
number_of_bindings_to_add = max_binding_per_bead - number_of_bindings
for bind in range(number_of_bindings_to_add):
if properties is None:
self.att_list += AttachmentSite(self.p_num[i])
else:
p = copy.deepcopy(properties)
for key, val in p.items(): # check whether a function was passed in
if hasattr(val, '__call__'):
args = inspect.getfullargspec(val).args
if len(args) == 0:
p[key] = val()
else:
d_arg = {}
for arg in args:
if arg == 'index':
d_arg[arg] = i
elif arg == 'bead':
d_arg[arg] = self.beads[i]
elif hasattr(self.beads[i], arg):
d_arg[arg] = getattr(self.beads[i], arg)
else:
d_arg[arg] = getattr(self, arg)
p[key] = val(**d_arg)
self.att_list += AttachmentSite(self.p_num[i], **p)
@deferred
def rotate(self, operation):
"""
Rotates the composite object
:param operation: rotation parsable by the quaternion package
:type operation: Quat, iterable or np.array
:return: None
"""
quat = Quat(operation) # parse the rotation using quaternion
rotation_matrix = quat.transform
for bead in self.beads:
bead.position = rotation_matrix.dot(bead.position)
self.att_list.rotate(operation)
for particle in self.rigid_body_references:
particle.quaternion = quat * particle.quaternion
self.beads[particle.p_num[0]].orientation = particle.quaternion.q_w_ijk
@deferred
def graft_external_objects(self,
external_object,
number=1,
connecting_topology=None,
qualification=None,
qualification_object=None,
merge_attachment_sites=False,
strict_compliance=False):
"""
Deferred handler for grafting objects. Grafts external objects unto this one
:param external_object: which external object to graft
:type external_object: CompositeObject
:param number: how many to graft, functions can take in attributes from the composite object
:type number: int or callable
:param connecting_topology: Which kind of bond to add
:type connecting_topology: str or BondType
:param qualification: which kind of site to graft to; None takes all available sites
:param qualification_object: which kind of site to graft to on the external object; None is unconstrained
:param merge_attachment_sites: Whether to merge the attachment sites of the external into this one
:type merge_attachment_sites: bool
:param strict_compliance: Whether to fully comply with qualifiers and restrictions on external_object
:type strict_compliance: bool
:return: Number of successful grafts
:rtype: int
"""
if hasattr(number, '__call__'):
number = number(**self.resolve_args(number))
for obj_number in range(number):
if hasattr(external_object, '__call__'):
cp = external_object(**self.resolve_args(external_object))
else:
cp = copy.deepcopy(external_object)
sites = self.att_list.get_qualified(qualification) # update the list of qualified sites
if not cp.is_built:
cp.build()
cp.build_finalize()
# get the pair list of bindings
binding_pairs = AttachmentSiteList.get_qualified_pair(self.att_list,
qualification,
cp.att_list,
qualification_object)
if len(binding_pairs) > 0:
binding_site, binding_site_cp = random.choice(binding_pairs)
elif len(sites) > 0 and not strict_compliance:
warnings.warn('Did not find a free binding site in bound object, trying to bind index 0')
binding_site = sites.pop(np.random.randint(0, len(sites)))
binding_site_cp = AttachmentSite(cp.p_num[0], orientation=np.array([0., 0., -1.]))
cp.att_list += binding_site_cp
else:
warnings.warn('Grafted ' + str(obj_number) + ' out of ' + str(number) + ' due to lack of binding sites',
UserWarning)
return obj_number
cp.att_list.bind_site(binding_site_cp, self)
self.att_list.bind_site(binding_site, cp)
self.merge(cp)
self_displacement = cp.positions[cp.att_list.index(binding_site_cp), :]
cp.positions -= self_displacement
# rotate the graft so that the alignment of the two attached beads are pointing towards each other
if binding_site['orientation'] is not None and binding_site_cp['orientation'] is not None:
rot_mat = Util.get_v1_v2_rotmat(binding_site['orientation'], -binding_site_cp['orientation'])
cp.rotate(rot_mat)
cp.positions += self_displacement
# move cp to the grafting site
displacement = self.positions[self.att_list.index(binding_site), :] - cp.positions[
cp.att_list.index(binding_site_cp), :]
cp.positions += displacement
# move the object just slightly away from the binding site
if binding_site['orientation'] is not None:
cp.positions += 0.5 * binding_site['orientation'] / np.linalg.norm(binding_site['orientation'])
if connecting_topology is not None:
if connecting_topology.__class__ == BondType:
self.b_types += connecting_topology
self.bonds += [str(connecting_topology), binding_site['index'], binding_site_cp['index']]
elif connecting_topology.__class__ == str:
self.bonds += [str(connecting_topology), binding_site['index'], binding_site_cp['index']]
else:
raise SyntaxError('Connecting topology is not parsable')
else:
self.bonds += ['Graft:' + str(self)+'-'+str(cp), binding_site['index'], binding_site_cp['index']]
if merge_attachment_sites:
self.att_list += cp.att_list
return number
@deferred
def constraints_to_bonds(self, topodict=None, unitdict=None):
"""
Transforms the constraints in the simulation by bonds. By default will create bonds with an energy_constant of
20000.0.; Constraints are wiped afterwards
:param topodict: topology properties of the created bonds
:param unitdict: units bond properties
:return: None
"""
for constraint in self.constraints:
if topodict is None:
topodict = {'distance_constant': constraint.distance, 'energy_constant': 20000.0}
else:
topodict.update({'distance_constant': constraint.distance})
if unitdict is None:
unitdict = {'distance_constant': 'L', 'energy_constant': 'E/LL'}
self.b_types += BondType('Constraint' + str(CompositeObject._CONSTRAINT_INDEX), topodict, unitdict)
self.bonds += ['Constraint' + str(CompositeObject._CONSTRAINT_INDEX), constraint.topology_tuple[0], constraint.topology_tuple[1]]
CompositeObject._CONSTRAINT_INDEX += 1
self.constraints.wipe()
[docs] def get_subgraph(self, require=None, avoid=None):
"""
Builds a subgraph of the connectivity where nodes attributes match require and does not match avoid. This method
requires the networkx package to work. Objects that are iterable are treated as OR logical requirements.
For instance, require={'beadtype': ['A','B']} will make a subgraph of nodes with beadtypes 'A' or 'B', while
avoid={'beadtype':['A','B']} will avoid both 'A' and 'B'
:param require: properties that are required to be part of the subgraph
:type require: dict
:param avoid: properties that cause rejection of the node from the subgraph
:type avoid: dict
:return: subgraph
:rtype: networkx.Graph
"""
if not NETWORKX_FOUND: # import check
return None
# create a new graph and add edges & nodes from the bonds & constraints
cp_graph = networkx.Graph()
cp_graph.add_edges_from(self.bonds.graph.edges)
cp_graph.add_edges_from(self.constraints.graph.edges)
cp_graph.add_nodes_from(self.bonds.graph.nodes)
cp_graph.add_nodes_from(self.constraints.graph.nodes)
# make a list of nodes to avoid
nodes_to_avoid = []
for node in cp_graph.nodes:
bead_node = self.beads[self.p_num.index(node)]
keep_node = True
if require is not None:
for key, val in require.iteritems():
key_satisfied = False
if hasattr(val, '__iter__'):
for v in val:
key_satisfied |= v == getattr(bead_node, key)
else:
key_satisfied = val == getattr(bead_node, key)
keep_node &= key_satisfied
delete_node = False
if avoid is not None:
for key, val in avoid.iteritems():
key_satisfied = False
if hasattr(val, '__iter__'):
for v in val:
key_satisfied |= v == getattr(bead_node, key)
else:
key_satisfied = val == getattr(bead_node, key)
delete_node |= key_satisfied
if not keep_node or delete_node:
nodes_to_avoid.append(node)
cp_graph.remove_nodes_from(nodes_to_avoid)
return cp_graph
[docs] def resolve_args(self, funct):
"""
Resolve arguments of a function by parsing arg names and gathering them from attributes of the composite object
:param funct: function to resolve
:type funct: callable
:return: argument list
:rtype: dict
"""
args = inspect.getargspec(funct).args
arguments = {}
for arg in args:
if arg != 'index':
try:
arguments[arg] = getattr(self, arg)
except AttributeError:
arguments[arg] = None
return arguments
@deferred
def add_angles(self, require=None, avoid=None, subgraph=None):
"""
Add angles to the system. This works by constructing a subgraph and generating all paths of length 3 in it. The
default mode of operation is to create a subgraph based on get_subgraph(require, avoid). A subgraph generator
can be supplied in the subgraph argument. It can be passed the bonded graph, constrains graph and the directed
graph of composite objects and beads through as 'bonded_graph', 'constraints_graph' and 'composite_graph'
respectively. Other attributes of the CompositeObject are also available.
:param require: List of required bead attributes (see get_subgraph)
:type require: dict
:param avoid: List of bead properties that cause node rejecction (see get_subgraph)
:type avoid: dict
:param subgraph: function that returns a networkx graph of all nodes where angle need to be added
:type subgraph: callable
:return: None
"""
if subgraph is not None and (require is not None or avoid is not None):
warnings.warn('add_dihedrals received a subgraph generator plus require/avoid values. These values will be '
'ignored', UserWarning)
if subgraph is None:
reduced_graph = self.get_subgraph(require, avoid)
else:
args = self.resolve_args(subgraph)
if 'bonded_graph' in args: args['bonded_graph'] = self.bonds.graph
if 'constraints_graph' in args: args['constraints_graph'] = self.constraints.graph
if 'composite_graph' in args: args['composite_graph'] = self.subunits_graph
all_paths = [p for p in networkx.all_pairs_dijkstra_path(reduced_graph, 2)]
for source, paths in all_paths:
for destination, path in paths.items():
if destination > source and len(path) == 3: # bidirectional, i->j & j->i, keep only one
anglename = 'Angle:'+str(self.beads[path[0]]) + '-' + str(self.beads[path[1]]) \
+ '-' + str(self.beads[path[2]])
self.angles += [anglename, path[0], path[1], path[2]]
@deferred
def add_dihedrals(self, require=None, avoid=None, subgraph=None):
"""
Add dihedrals to the system. This works by constructing a subgraph and generating all paths of length 4 in it.
The default mode of operation is to create a subgraph based on get_subgraph(require, avoid). A subgraph
generator can be supplied in the subgraph argument. It can be passed the bonded graph, constrains graph and the
directed graph of composite objects and beads through as 'bonded_graph', 'constraints_graph' and
'composite_graph' respectively. Other attributes of the CompositeObject are also available.
:param require: List of required bead attributes (see get_subgraph)
:type require: dict
:param avoid: List of bead properties that cause node rejecction (see get_subgraph)
:type avoid: dict
:param subgraph: function that returns a networkx graph of all nodes where dihedrals need to be added
:type subgraph: callable
:return: None
"""
if subgraph is not None and (require is not None or avoid is not None):
warnings.warn('add_dihedrals received a subgraph generator plus require/avoid values. These values will be '
'ignored', UserWarning)
if subgraph is None:
reduced_graph = self.get_subgraph(require, avoid)
else:
args = self.resolve_args(subgraph)
if 'bonded_graph' in args: args['bonded_graph'] = self.bonds.graph
if 'constraints_graph' in args: args['constraints_graph'] = self.constraints.graph
if 'composite_graph' in args: args['composite_graph'] = self.subunits_graph
all_paths = [p for p in networkx.all_pairs_dijkstra_path(reduced_graph, 3)]
for source, paths in all_paths:
for destination, path in path.iteritems():
if destination > source and len(path) == 4:
diename = 'Dihedral:'+str(self.beads[path[0]]) + '-' + str(self.beads[path[1]]) + '-' \
+ str(self.beads[path[2]]) + '-' + str(self.beads[path[3]])
self.dihedrals += [diename, path[0], path[1], path[2], path[3]]
@deferred
def dihedral_to_improper(self, type=None):
"""
Replaces dihedrals in the topology by impropers
:param type: which type to replace, None selects all (default None)
:type type: str
:return: None
"""
if type is None:
typelist = copy.deepcopy(self.d_types.topology_list.keys())
for key in typelist:
self.dihedral_to_improper(key)
else:
die_type = self.d_types[type]
self.i_types += ImproperType(name=copy.copy(die_type.typename),
topodict=copy.copy(die_type.topology_constants),
unitdict=copy.copy(die_type.topology_units))
for die in self.dihedrals:
if die.topology_name == type:
self.impropers += die
self.dihedrals.topologies[:] = [item for item in self.dihedrals.topologies if item.topology_name != type]
del self.d_types.topology_list[type]
[docs] def save(self, filepath):
"""
Save the current instance to a file
:param filepath: path to file
:type filepath: str
:return: None
"""
with open(filepath, 'wb') as pickle_file:
pickle.dump(self, pickle_file)
[docs] def override(self, filepath):
"""
Overrides the current object with one taken from a file
:param filepath: file to load
:type filepath: str
:return: None
"""
with open(filepath, 'rb') as pickle_file:
tobj = pickle.load(pickle_file)
self.__dict__.clear()
self.__dict__.update(tobj.__dict__)
[docs] @staticmethod
def load(filepath):
"""
Returns the object contained in a file
:param filepath: path to file
:type filepath: str
:return: object in file
:rtype: CompositeObject
"""
with open(filepath, 'rb') as pickle_file:
t_obj = pickle.load(pickle_file)
return t_obj
[docs] @staticmethod
def import_mbuild(mbuild_object, units=None):
"""
Imports a muild object as a Hoobas composite object
:param mbuild_object: mbuild object to import as Hoobas composite
:type mbuild_object: mbuild.compound
:param units: units of the mbuild object (default None)
:type units: SimulationUnits
:return: Hoobas representation of mbuild_object
:rtype: CompositeObject
"""
composite = CompositeObject(units)
for subcompound in mbuild_object.children:
if not subcompound.port_particle: # avoid ports here
composite.merge(CompositeObject.import_mbuild(subcompound))
# in priciple we could iterate using compound.particles, but structures won't match
if not mbuild_object.children:
composite.beads += Bead(beadtype=str(mbuild_object.name),
position=mbuild_object.pos,
charge=mbuild_object.charge)
if mbuild_object.rigid_id:
composite.beads[0].body = mbuild_object.rigid_id
if mbuild_object.rigid_id:
composite._body = mbuild_object.rigid_id
# copy bonds from parent
if not mbuild_object.parent:
subs = list(mbuild_object.children)
for bond in mbuild_object.bonds():
# sort strings to get consistent string names
bondtype = bond[0].name+'-'+bond[1].name if bond[0].name > bond[1].name else bond[1].name+'-'+bond[0].name
composite.bonds += [bondtype, subs.index(bond[0]), subs.index(bond[1])]
for port in mbuild_object.referenced_ports():
if not port.used:
direction = np.mean(port.xyz_with_ports[0:4, :], axis=0)
direction /= np.linalg.norm(direction)
composite.att_list += AttachmentSite(subs.index(port.anchor), orientation=direction)
composite.is_built = True
return composite