# coding: utf-8
import copy
import random
from itertools import chain
import warnings
import numpy as np
from hoobas import CGBead
from hoobas.Util import get_rotation_matrix
from hoobas.Quaternion import Quat
import inspect
from hoobas import Composite
from hoobas.Units import SimulationUnits as SimUnits
[docs]class Colloid(Composite.CompositeObject):
"""
Object meant for holding CompositeObject comprised of both rigid and flexible parts.
Args:
center_type (str): Additional rigid beadtype prefixed to the colloid
surface_type (str): Beadtype to use when creating bead lists from the shape
shape (GenShape.Shape): Shape to use for constructing the rigid body
"""
def __init__(self, center_type='W', surface_type='P', shape=None, **kwargs):
super(Colloid, self).__init__(**kwargs)
self.c_type = center_type
self.s_type = surface_type
self.flags = {}
self.surf_tags = []
# associated shape, contains list of surface, and directives
self.shape = copy.deepcopy(shape)
self.shape.change_units(self.units)
self.CONST_SHAPE_TABLE = copy.deepcopy(shape.table)
# check if the shape has a supplementary building method that has to be run but wasnt
if hasattr(self.shape, 'BuiltFlag') and self.shape.BuiltFlag is False:
self.shape.BuildMethod()
warnings.warn('Colloid found a non-built shape object and had to build it', UserWarning)
# link diagI with the bead properties
self.diagI = np.array(self.shape.Itensor)
# colloid quaternion is shared with local object quaternion
self.quaternion = self.shape.quaternion
# figure out the mass of the colloid
if 'mass' in kwargs: # mass override
self.mass = kwargs['mass']
elif 'mass' in self.shape.flags:
self.mass = self.shape.flags['mass']
elif 'density' in kwargs and 'volume' in self.shape.flags: # this is for normalized shapes
self.mass = kwargs['density'] * self.shape.flags['volume']
else:
warnings.warn('Unable to determine rigid mass, using value of 10, something is wrong here')
self.mass = 10.0
#####################
# body properties
#####################
self.body_num = 0
self.body_beads = []
self.body_mass = 0.0
self.__initial_rotation()
def build(self):
if hasattr(self.shape, 'BuiltFlag') and self.shape.BuiltFlag is False:
self.shape.BuildMethod()
super(Colloid, self).build()
self._body = 0
@property
def surface_type(self):
return self.s_type
@property
def center_type(self):
return self.beads[0].beadtype
@property
def shape_class(self):
return type(self.shape).__name__
@property
def shape_flag(self):
return self.shape.flags
@property
def body(self):
return self._body
@body.setter
def body(self, val):
self._body = val
for rigid_bead in self.body_beads:
rigid_bead.body = val
@property
def body_typelist(self):
return [bead.beadtype for bead in self.body_beads]
[docs] def relative_positions(self):
"""
Retrieve relative positions with respect to particle center. This is meant to produce tuples for hoomd rigid
body constraints command
:return: list[tuple]
"""
_ = []
for position in self.CONST_SHAPE_TABLE:
_.append((float(position[0]), float(position[1]), float(position[2])))
return _
# determine the initial rotation the colloid has based on its orientation
def __initial_rotation(self):
_mat = self.quaternion.transform
for pidx in range(self.shape.table.__len__()):
lvec = copy.deepcopy(self.shape.table[pidx])
lvec = _mat.dot(lvec)
self.shape.table[pidx] = lvec
@Composite.deferred
def rotate(self, operation):
"""
rotates the colloid
:param operation: transform operation
:return:
"""
translation = copy.copy(self.positions[0, :])
self.positions += -translation
q_op = Quat(operation)
super(Colloid, self).rotate(q_op)
self.positions += translation
self.quaternion = q_op * self.quaternion
self.beads[0].orientation = self.quaternion.q_w_ijk
[docs]class SimpleColloid(Colloid):
"""
This is the class used for simple rigid bodies, i.e., polyhedra, with a single surface atom type. A single
polydispersity can be supplied, which is obtained by passing a function for the size
"""
def __init__(self, size, **kwargs):
super(SimpleColloid, self).__init__(**kwargs)
self.size = size
def build(self):
super(SimpleColloid, self).build()
if hasattr(self.size, '__call__'):
self.size = self.size(**self.resolve_args(self.size))
if not isinstance(self.size, float):
self.size = float(self.size)
warnings.warn('Had to coerce size to floating point value')
# mass in shape class is normalized; specifying mass is an override
if 'mass' not in self.shape.flags:
self.mass *= self.size ** 3
# set the surface mass to be 3/5 of the overall mass to fix the rigid body intertia, this is not used anymore
# in hoomd
if self.shape.num_surf > 0:
self.s_mass = self.mass * 3.0 / 5.0 / self.shape.num_surf
else:
self.s_mass = 0.0
self.body_mass = self.mass
# non-rotating frame for aggregation of rigid properties
self.CONST_SHAPE_TABLE *= self.size
self.p_num = [0]
self.diagI *= self.mass * (self.size ** 2.0)
# the center particle holds the rigid body structure, make sure that mass is not a ref to self.mass for units
self.beads += CGBead.Bead(position=np.array([0.0, 0.0, 0.0], dtype=np.float32), beadtype=self.c_type,
body=0, mass=self.mass, quaternion=self.quaternion,
moment_inertia=self.diagI)
self.body_beads.append(self.beads[-1])
for i in range(len(self.shape.pos)):
self.p_num.append(self.p_num[-1] + 1)
self.beads += CGBead.Bead(position=copy.deepcopy(self.shape.pos[i] * self.size),
beadtype=self.s_type, body=0, mass=self.s_mass)
self.body_beads.append(self.beads[-1])
@Composite.deferred
def additional_directives(self, directive):
if not issubclass(directive.__class__, AdditionalDirective):
raise SyntaxError('Directive supplied not a subclass of the additional directive class')
# apply directive to all points in apply_to_types
_norm = 0.0
for bead in self.beads:
if directive.apply_to_types is None or bead.beadtype in directive.apply_to_types:
setattr(bead, directive.attribute,
directive.function(bead.position[0], bead.position[1], bead.position[2]))
_norm += getattr(bead, directive.attribute)
# if a normalization has been supplied, enforce it
if directive.normalization is not None:
for bead in self.beads:
if directive.apply_to_types is None or bead.beadtype in directive.apply_to_types:
setattr(bead, directive.attribute,
getattr(bead, directive.attribute) * directive.normalization / _norm)
[docs]class MultiDisperseColloid(Colloid):
"""
Class for colloids requiring dispersity in multiple variables. For instance rods with aspect ratios and sizes, which
may have non-zero covariance.
Args:
shape_property_function: function which yields properties of the shape in a dict, must at least return a key
for size
"""
def __init__(self, shape_property_function, **kwargs):
super(MultiDisperseColloid, self).__init__(**kwargs)
self.shape_function = shape_property_function
self.shape_args = kwargs
self.size = 0.0
def build(self):
super(MultiDisperseColloid, self).build()
if not hasattr(self.shape, 'BuiltFlag') or self.shape.BuiltFlag is True:
raise AssertionError("Shape is already built while specifying runtime building parameters")
if not hasattr(self.shape_function, '__call__'):
raise SyntaxError('Cannot call function shape_function')
# do some parsing on the shape function to figure out its arguments
if len(inspect.getargspec(self.shape_function)[0]) == 0:
shape_prop = self.shape_function()
else:
shape_prop = self.shape_function(**self.shape_args)
if 'size' not in shape_prop:
raise SyntaxError('size is not defined')
self.size = shape_prop['size']
for prop in shape_prop:
setattr(self.shape, prop, shape_prop[prop])
self.shape.BuildMethod()
self.CONST_SHAPE_TABLE = copy.deepcopy(self.shape.table)
self.orientation = self.shape.n_plane
self.diagI = self.shape.Itensor
# colloid quaternion is shared with local object quaternion
self.quaternion = self.shape.quaternion
# get mass from shape object
try:
self.mass = self.shape.flags['mass']
except KeyError:
try:
self.mass = self.shape.flags['density'] * self.shape.flags['volume'] * (self.size ** 3)
except KeyError:
warnings.warn('Unable to determine solid body mass, using value of 10.0, something is wrong here',
UserWarning)
self.mass = 10.0
# all colloids have well defined orientations with respect to the base shape function
if self.orientation is None:
self.orientation = [0, 0, 1]
# set the surface mass to be 3/5 of the overall mass to fix the rigid body intertia
if self.shape.num_surf > 0:
self.s_mass = self.mass * 3.0 / 5.0 / self.shape.num_surf
else:
self.s_mass = 0.0
self.body_mass = self.mass
# non-rotating frame for aggregation of rigid properties
self.CONST_SHAPE_TABLE *= self.size
# the center particle holds the rigid body structure
self.diagI *= self.mass * (self.size ** 2.0)
# mass in shape class is normalized
if 'mass' not in self.shape.flags:
self.mass *= self.size ** 3
self.p_num = [0]
self.beads += CGBead.Bead(position=np.array([0.0, 0.0, 0.0], dtype=np.float32), beadtype=self.c_type,
body=0, mass=self.mass, quaternion=self.quaternion,
moment_inertia=self.diagI)
self.body_beads += self.beads[-1]
for i in range(len(self.shape.pos)):
self.p_num += self.p_num[-1] + 1
self.beads += CGBead.Bead(position=copy.deepcopy(self.shape.pos[i] * self.size),
beadtype=self.s_type, body=0, mass=self.s_mass)
self.body_beads += self.beads[i]
@Composite.deferred
def additional_directives(self, directive):
if not issubclass(directive.__class__, AdditionalDirective):
raise SyntaxError('Directive supplied not a subclass of the additional directive class')
# apply directive to all points in apply_to_types
_norm = 0.0
for bead in self.beads:
if directive.apply_to_types is None or bead.beadtype in directive.apply_to_types:
setattr(bead, directive.attribute,
directive.function(bead.position[0], bead.position[1], bead.position[2]))
_norm += getattr(bead, directive.attribute)
# if a normalization has been supplied, enforce it
if directive.normalization is not None:
for bead in self.beads:
if directive.apply_to_types is None or bead.beadtype in directive.apply_to_types:
setattr(bead, directive.attribute,
getattr(bead, directive.attribute) * directive.normalization / _norm)
[docs]class ComplexColloid(Colloid):
"""
Provides an object to build complex colloids defined by shells. This is generally used in conjunction with
:py:class:`hoobas.GenShape.PDBProtein`.
"""
def __init__(self, **kwargs):
super(ComplexColloid, self).__init__(**kwargs)
self.nshells = 0
self.s_mass = 0.0
def build(self):
super(ComplexColloid, self).build()
# mass is set by the I_fixer from the base shape
self.s_mass = 1.0
self.body_mass = self.mass
# A complex colloid should have a moment of inertia defined by the I_fixer
self.beads += CGBead.Bead(position=np.array([0.0, 0.0, 0.0]), beadtype=self.c_type, body=0,
mass=self.body_mass, quaternion=self.quaternion,
moment_inertia=self.diagI)
# in case we have more stuff defined (note additional point[0] is the center)
for i in range(1, self.shape.additional_points.__len__()):
self.beads += CGBead.Bead(position=self.shape.additional_points[i],
beadtype=self.c_type + self.shape.type_suffix[i],
body=0, mass=self.shape.masses[i])
self.p_num.append(self.p_num[-1] + 1)
# check for the number of shells in the shape
self.nshells = len(self.shape.shells)
# the base beads are part of the body
self.body_beads += self.beads
# this colloid is defined by shells. For instance shell with name 'shell1' can correspond to certain grafting
# points. This is used to enforce binding of ligands on specific points for instance.
for shell in self.shape.shells:
for index, item in enumerate(shell['positions']):
self.p_num.append(self.p_num[-1] + 1)
self.beads += CGBead.Bead(position=item, beadtype=self.s_type, body=self.body, mass=self.s_mass)
self.att_list += Composite.AttachmentSite(index=index, qualifier=shell['qualifier'])
self.body_beads.append(self.beads[-1])
[docs]class FromComposite(Composite.CompositeObject):
"""
Creates a rigid body out of a :py:class:`hoobas.Composite.CompositeObject` prototype
"""
def __init__(self, composite, units=None):
super(FromComposite, self).__init__(units=units)
self._composite = composite
self.CONST_SHAPE_TABLE = np.zeros((0, 3))
self.body_beads = []
self.center_type = ''
self.mass = 0.0
self.body = 0
@property
def body(self):
return self._body
@body.setter
def body(self, value):
for bead in self.body_beads:
bead.body = value
self._body = value
@property
def body_typelist(self):
return [bead.beadtype for bead in self.body_beads]
def build(self):
self._composite.build()
self._composite.build_finalize()
self.merge(self._composite)
for bead in self.beads:
if bead.body != -1:
self.body_beads.append(bead)
self.center_type = self.body_beads[0].beadtype
# compute moment of inertia and rotation
inertia_tensor = np.zeros((3, 3), dtype=np.float)
for bead in self.body_beads:
point = bead.position - self.body_beads[0].position
inertia_tensor[0, 0] += point[1] * point[1] + point[2] * point[2]
inertia_tensor[1, 1] += point[0] * point[0] + point[2] * point[2]
inertia_tensor[2, 2] += point[0] * point[0] + point[1] * point[1]
inertia_tensor[0, 1] -= point[0] * point[1]
inertia_tensor[0, 2] -= point[0] * point[2]
inertia_tensor[1, 2] -= point[1] * point[2]
inertia_tensor[1, 0] = inertia_tensor[0, 1]
inertia_tensor[2, 0] = inertia_tensor[0, 2]
inertia_tensor[2, 1] = inertia_tensor[1, 2]
inertia_tensor = inertia_tensor.astype(dtype=np.float32)
_maxel = np.max(inertia_tensor)
for dim1 in range(3):
for dim2 in range(3):
if abs(inertia_tensor[dim1, dim2]) < 2.0 * np.finfo(np.float32).eps * _maxel:
inertia_tensor[dim1, dim2] = 0.0
w, v = np.linalg.eigh(inertia_tensor * 72.0)
# sort the eigenvalues
idx = w.argsort(kind='mergesort')[::-1]
w = w[idx]
v = v[:, idx]
# invert one axis if the coordinate system has become left-handed
if np.linalg.det(v) < 0:
v[:, 0] = -v[:, 0]
for bead in self.body_beads:
self.mass += bead.mass
self.body_beads[0].mass = self.mass
self.rotate(Quat(v).inv())
self.beads[0].moment_inertia = np.array([w[0], w[1], w[2]])
self.CONST_SHAPE_TABLE = np.zeros((len(self.body_beads)-1, 3))
# set the const table
for i, bead in enumerate(self.body_beads):
if i == 0:
continue
self.CONST_SHAPE_TABLE[i - 1, :] = bead.position - self.beads[0].position
super(FromComposite, self).build()
[docs] def relative_positions(self):
"""
Retrieve relative positions with respect to particle center. This is meant to produce tuples for hoomd rigid
body constraints command
:return: list[tuple]
"""
_ = []
for position in self.CONST_SHAPE_TABLE:
_.append((float(position[0]), float(position[1]), float(position[2])))
return _
def change_units(self, new):
m_len = new.get_length(1.0, self.units.lunit)
self.CONST_SHAPE_TABLE *= m_len
super(FromComposite, self).change_units(new)
class AdditionalDirective(object):
def __init__(self):
self.function = lambda self, x, y, z: 1
self.attribute = 'charge'
self.apply_to_types = None
self.normalization = None