# coding: utf-8
import copy
import warnings
import numpy as np
from hoobas import CGBead
from hoobas.Quaternion import Quat
import inspect
from hoobas import Composite
from hoobas.Units import SimulationUnits as SimUnits
from typing import Union
from typing import List
from typing import Callable
from hoobas.GenShape import shape as Shape
[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: str = 'W', surface_type: str = 'P', shape: Union[None, 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) -> str:
return self.s_type
@property
def center_type(self) -> str:
return self.beads[0].beadtype
@property
def shape_class(self):
return type(self.shape).__name__
@property
def shape_flag(self) -> dict:
return self.shape.flags
@property
def body(self) -> Union[None, int]:
return self._body
@body.setter
def body(self, val: Union[None, int]) -> None:
self._body = val
for rigid_bead in self.body_beads:
rigid_bead.body = val
@property
def body_typelist(self) -> List[str]:
return [bead.beadtype for bead in self.body_beads]
[docs] def relative_positions(self) -> List[tuple]:
"""
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) -> None:
_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) -> None:
"""
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: Union[float, Callable[..., float]], **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: Callable, **kwargs):
super(MultiDisperseColloid, self).__init__(**kwargs)
self.shape_function = shape_property_function
self.shape_args = kwargs
self.size = 0.0
def build(self) -> None:
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.getfullargspec(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: Composite.CompositeObject, units: SimUnits = 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) -> Union[None, int]:
return self._body
@body.setter
def body(self, value: Union[None, int]) -> None:
for bead in self.body_beads:
bead.body = value
self._body = value
@property
def body_typelist(self) -> List[str]:
return [bead.beadtype for bead in self.body_beads]
def build(self) -> None:
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) -> List[tuple]:
"""
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