# coding: utf-8
from math import *
import random
import copy
import types
import warnings
from itertools import chain
from hoobas import SimulationDomain
import numpy as np
from hoobas import CGBead
from hoobas import Composite
from hoobas import Layers
from hoobas.Util import get_rotation_matrix
from hoobas.Util import gen_random_mat
from hoobas import Util
from hoobas.Units import SimulationUnits as SimUnits
import inspect
import functools
import os
[docs]class Builder(Composite.CompositeObject):
"""
General building methods of Hoobas. This class only assembles the system and related commands. It is composed with
the following classes for specialized post-processing:
Electrostatics: charge normalization (permittivity)
Bias: adding biasing potentials in the form of bonds, fixed Einstein lattices, etc.
FileWriter: writing simulation information to files
Args:
domain: Simulation domain used
"""
_LICENSE_DISPLAY = True
@staticmethod
def display_license():
this_path = os.path.dirname(os.path.realpath(__file__))
with open(os.path.join(this_path, '..', 'LICENSE'), 'r') as f:
lic = f.read()
print(lic)
f.close()
def __init__(self, domain, **kwargs):
if 'units' in kwargs and issubclass(kwargs['units'].__class__, SimUnits):
units = kwargs['units']
else:
units = SimUnits()
super(Builder, self).__init__(units)
#################################
# used for internal representation
#################################
self.__types = []
self.__obj_index = []
self.__obj_end_index = []
# DEFAULT ARGUMENTS
self.z_multiplier = 1.0
self.filename = 'HOOBAS_FILE'
# list of properties defined in the Hoomd xml formats. First list refers to particle properties,
# second to bonded interaction types
self.properties = ['velocity', 'acceleration', 'diameter', 'charge', 'body', 'orientation', 'angmom',
'moment_inertia', 'image']
self.topology_properties = ['bonds', 'angles', 'dihedrals', 'impropers']
# Display the license if it hasn't been displayed already
if Builder._LICENSE_DISPLAY:
Builder.display_license()
Builder._LICENSE_DISPLAY = False
# overriding defaults
for key in kwargs:
setattr(self, key, kwargs.get(key))
# set units
domain.change_units(self.units)
self.impose_box = []
self.domain = copy.deepcopy(domain)
self.ref_crystal = []
###########################################
# Data for salt dielectric constant values
###########################################
self.c_salt = 0.0
# types that comes from rigid bodies that will be created on the fly by hoomd
self.ext_rigid_types = []
# list of beads of atoms that are set to be built by hoomd
self.ext_atom_types = []
self.ext_layers = []
self.is_built = True # builders are top level objects and are always built or builing other things
# dealing with additional modifications
self.Electrostatics = Electrostatics(self.beads, self.units)
self.FileWriter = FileWriter(self)
self.Bias = Bias(self)
# do the building
self.__build_domain()
# this is a name overload
@property
def sys_box(self):
"""
Queries the simulation domain of the builder
:return: box
:rtype: iterable
"""
return self.current_box()
@property
def bead_types(self):
"""
Beadtypes in the simulation
:return: beadtypes
:rtype: list
"""
p_typelist = copy.deepcopy(self.ext_rigid_types)
for bead in self.beads:
p_typelist += [bead.type]
return list(set(p_typelist))
@property
def center_types(self):
"""
Beadtypes of centers of rigid bodies in the simulation
:return: beadtypes
:rtype: list
"""
_ct = []
for i in range(self.rigid_body_references.__len__()):
if not isinstance(self.rigid_body_references[i].center_type, types.StringTypes):
_ct.append(self.rigid_body_references[i].center_type)
else:
_ct.append([self.rigid_body_references[i].center_type])
return list(set(list(chain.from_iterable(_ct))))
@property
def surface_types(self):
"""
Tags of the surfaces of rigid bodies in the simulation
:return: tags
:rtype: list
"""
_st = []
for i in range(self.rigid_body_references.__len__()):
if not isinstance(self.rigid_body_references[i].surface_type, types.StringTypes):
_st.append(self.rigid_body_references[i].surface_type)
else:
_st.append([self.rigid_body_references[i].surface_type])
return list(set(list(chain.from_iterable(_st))))
@property
def center_tags(self):
"""
Tags of the center of rigid bodies in the simulation
:return: tags
:rtype: list
"""
_tags = []
for i in range(self.rigid_body_references.__len__()):
_tags.append(self.rigid_body_references[i].pnum_offset)
return _tags
@property
def num_beads(self):
"""
Number of beads in the simulation
:return: number of beads
:rtype: int
"""
return len(self.beads)
@property
def charge_list(self):
"""
Returns the list of charged beadtypes in the simulation
:return: charged beadtypes
:rtype: list
"""
_tlist = []
for bead in self.beads:
if abs(bead.charge) > 1e-5:
_tlist.append([bead.beadtype])
return list(set(list(chain.from_iterable(_tlist))))
@property
def obj_start_index(self):
"""
Queries the list of starting indexes of each object added to the simulation. This can be used to create groups
in HOOMD-blue for instance, in conjunction with ending indexes
:return: Starting indexes
:rtype: list
"""
return self.__obj_index
@property
def obj_end_index(self):
"""
Queries the list of last indexes of each object added to the simulation. This can be used to create groups in
HOOMD-blue, when used in conjunction with the starting indexes
:return: Ending indexes
:rtype: list
"""
return self.__obj_end_index
[docs] def shift_pos(self, displacement):
"""
shift positions of all beads by a constant displacement
:param displacement: Displacement to apply to all beads
:return:
"""
for bead in self.beads:
bead.position += displacement
[docs] def set_rotation_function(self, mode=None, key_search=None):
"""
Sets the rotation of rigid bodies in the domain. Note that all rigid bodies will be rotated according to their
orientations, which includes rigid objects that are grafted on other objects. For these cases, the key_search
option must be employed.
:param mode: (string, None), either 'none' or 'random'
:param key_search: (string, dict), used for random list of properties the rigid body must have in order to be rotated.
:return:
"""
if mode is None or mode == 'none':
for ref in self.rigid_body_references:
ref.rotate(get_rotation_matrix(ref.orientation))
if mode == 'random' and key_search is None:
for item in self.subunits:
item.rotate(gen_random_mat())
if mode == 'random' and key_search == 'rigid':
for ref in self.rigid_body_references:
ref.rotate(gen_random_mat())
if mode == 'random' and hasattr(key_search, '__iter__'):
if 'rigid' in key_search and key_search['rigid']:
for ref in self.rigid_body_references:
ok_to_rotate = True
for prop in key_search:
if hasattr(ref, prop):
ok_to_rotate &= getattr(ref, prop) == key_search[prop]
else:
ok_to_rotate = False
warnings.warn('Build : attribute ' + str(prop) + ' not found', SyntaxWarning)
if ok_to_rotate:
ref.rotate(gen_random_mat())
else:
for ref in self.subunits:
ok_to_rotate = True
for prop in key_search:
if hasattr(ref, prop):
ok_to_rotate &= getattr(ref, prop) == key_search[prop]
else:
ok_to_rotate = False
warnings.warn('Build : attribute ' + str(prop) + ' not found', SyntaxWarning)
if ok_to_rotate:
ref.rotate(gen_random_mat())
def __build_domain(self):
"""
performs the building based on the supplied simulation domain
:return:
"""
self.__obj_index = []
self.__obj_end_index = []
for item in self.domain.table:
center_position = item['position']
composite = copy.deepcopy(item['object'])
self.merge(composite)
composite.change_units(self.units)
composite.center_position = center_position
if composite.body != -1:
composite.body = composite.pnum_offset
self.ref_crystal.append(composite)
self.__obj_index.append({'index': composite.pnum_offset, 'object': composite})
self.__obj_end_index.append({'index': self.p_num[-1], 'object': composite})
self.__types.append(str(composite))
def current_box(self):
if self.impose_box.__len__() == 0:
vx, vy, vz = self.domain.direct_lattice_vectors
if 'vertical_slice' in self.domain.flags and self.domain.flags['vertical_slice']:
L = Util.c_hoomd_box([vx, vy, vz], self.domain.bounds, z_multi=5.0)
else:
L = Util.c_hoomd_box([vx, vy, vz], self.domain.bounds)
return L
else:
return self.impose_box
[docs] def copy_ext_obj(self, ext_obj):
"""
Directly imports a composite object into the simulation domain without performing any changes
:param ext_obj: object to import
:return:
"""
composite_copy = copy.deepcopy(ext_obj)
self.merge(composite_copy)
[docs] def add_N_ext_obj(self, ext_obj, N):
"""
Adds N composite objects to the system positioned at random locations, with random rotation
:param ext_obj: external object to be copied.
:param N: number of objects
:return:
"""
# create the current box, set all centers of build-like objects inside, set everything by periodicBC
L = self.current_box()
for n in range(N):
composite_copy = copy.deepcopy(ext_obj)
if hasattr(composite_copy, 'do_on_copy') and callable(getattr(composite_copy, 'do_on_copy')):
composite_copy.do_on_copy()
self.merge(composite_copy)
composite_copy.rotate(Util.gen_random_mat())
random_position = np.array([random.uniform(-L[0]/2.0, L[0]/2.0),
random.uniform(-L[1]/2.0, L[1]/2.0),
random.uniform(-L[2]/2.0, L[2]/2.0)])
composite_copy.center_position = random_position
[docs] def add_rho_molar_ions(self, rho, qtype='ion', ion_mass=1.0, q=1.0, ion_diam=1.0):
"""
Adds a volumetric density of rho ions to the simulation domain
:param rho: molar density (mol / L)
:type rho: float
:param qtype: beadtype for added ions
:type qtype: basestring
:param ion_mass: mass of added ions
:type ion_mass: float
:param ion_diam: size of ions
:type ion_diam: float
:return: None
"""
# rho *= (6.02*10**23) / (10**3*(2*10**-9)**3)**-1
L = self.current_box()
if L.__len__() == 6:
_ = np.array([[L[0], L[1] * L[3], L[2] * L[4]], [0.0, L[1], L[2] * L[5]], [0.0, 0.0, L[2]]],
dtype=np.float32)
elif L.__len__() == 3:
_ = np.array([[L[0], 0.0, 0.0], [0.0, L[1], 0.0], [0.0, 0.0, L[2]]], dtype=np.float32)
else:
_ = np.array([[L[0], 0.0, 0.0], [0.0, L[0], 0.0], [0.0, 0.0, L[0]]], dtype=np.float32)
_v = 0.0
for col in self.rigid_body_references:
if 'volume' in col.shape_flag:
_v += col.shape_flag['volume'] * col.size**3.0
V = np.linalg.det(_) - _v # volume of the box - excluded volume of solid particles
# 4.81 is magic number of units for mol / L to units of l = 2 nm
N = int(rho * 4.81 * (self.units.get_length(1.0, 'R')**-3.0) * V)
self.add_ions(N=N, ion_mass=ion_mass, ion_diam=ion_diam, q=q, qtype=qtype)
[docs] def add_ions(self, N, qtype='ion', ion_mass=1.0, q=1.0, ion_diam=1.0):
"""
Adds N ions in the simulation, randomly positioned
:param N: Number of ions
:type N: int
:param qtype: ion beadtype
:type qtype: str
:param ion_mass: mass of ions
:type ion_mass: float
:param q: charge of ions
:type q: float
:param ion_diam: diameter of ions
:type ion_diam: float
:return: None
"""
L = self.current_box()
beadlist = []
for i in range(N):
_rej_check = True
_gen_pos = np.array([random.uniform(-L[0] / 2.01, L[0] / 2.01), random.uniform(-L[1] / 2.01, L[1] / 2.01),
random.uniform(-L[2] / 2.01, L[2] / 2.01)])
while _rej_check:
_rej_check = False
for j in range(self.rigid_body_references.__len__()):
try:
if np.linalg.norm(self.rigid_body_references[j].center_position - _gen_pos) < self.rigid_body_references[j]._sh.flags['hard_core_safe_dist'] * self.rigid_body_references[j]._sh.flags['size']:
_gen_pos = np.array(
[random.uniform(-L[0] / 2.01, L[0] / 2.01), random.uniform(-L[1] / 2.01, L[1] / 2.01),
random.uniform(-L[2] / 2.01, L[2] / 2.01)])
_rej_check = True
except KeyError:
pass
beadlist.append(CGBead.Bead(position=_gen_pos, beadtype=qtype, mass=ion_mass, charge=q, diameter=ion_diam))
try:
self.p_num.append(self.p_num[-1]+1)
except IndexError:
self.p_num.append(0)
self.beads += beadlist
[docs] def fix_remaining_charge(self, ptype='ion', ntype='ion', pion_mass=1.0, nion_mass=1.0,
qp=1.0, qn=-1.0, pion_diam=1.0, nion_diam=1.0, isrerun=False):
"""
Enforces electroneutrality of the system
:param ptype: type name of cation
:param ntype: type name of anion
:param pion_mass: mass of cation
:param nion_mass: mass of anion
:param qp: charge of cation
:param qn: charge of anion
:param pion_diam: diameter of cation
:param nion_diam: diameter of anion
:param isrerun: internal parameter; should not be used from outside
:return: None
"""
_internal_charge = 0.0
for i in range(self.beads.__len__()):
_internal_charge += self.beads[i].charge
if _internal_charge == 0:
return
elif _internal_charge > 0:
_rerun_check = (_internal_charge % qn == 0)
self.add_ions(N=int(ceil(-_internal_charge / float(qn))), qtype=ntype, ion_mass=nion_mass, q=qn,
ion_diam=nion_diam)
else:
_rerun_check = (_internal_charge % qp == 0)
self.add_ions(N=int(ceil(-_internal_charge / float(qp))), qtype=ptype, ion_mass=pion_mass, q=qp,
ion_diam=pion_diam)
if not isrerun and _rerun_check:
self.fix_remaining_charge(ptype, ntype, pion_mass, nion_mass, qp, qn, pion_diam, nion_diam, isrerun=True)
elif _rerun_check:
warnings.warn('Cannot fix the intrinsic charge in the system', UserWarning)
[docs] def add_layer(self, multilayer, midplane=0.0, displace_elements=True):
"""
Adds a 2D object in the simulation box. It is automatically added in the (001) plane, crossing the z=0 plane.
For most simulation boxes, this is a flat plane laying in the XY plane. By default, the box size is increased
and the simulation contents that would had overlapped are pushed away. This behavior is controlled by
displace_elements
:param multilayer: 2D layer to add
:type multilayer: Layers.Multilayer
:param midplane: location of the plane
:type midplane: float
:param displace_elements: whether to displace existing contents away from the layer
:type displace_elements: bool
:return: None
"""
# if displace is set to True, make an opening in whatever is currently in by compressing the whole box;
# the layer should be periodic only in 2D so push along vz
#length of box along vz is 2 * bounds[2] * vz so we can just calculate thickness and increase bounds
# by a correct amount
if not issubclass(multilayer.__class__, Layers.Multilayer):
warnings.warn('Supplied layer object is not of layers class, results may vary', UserWarning)
ml = copy.deepcopy(multilayer)
ml.change_units(self.units)
ml._recalculate_numbers_from_area(np.linalg.norm(np.cross(self.domain.vx, self.domain.vy)))
if not ml.is_built:
ml.build()
ml.transform(np.array([np.array(self.domain.vx) * 2.0 * self.domain.bounds[0],
np.array(self.domain.vy) * 2.0 * self.domain.bounds[1],
np.array(self.domain.vz) * 2.0 * self.domain.bounds[2]]),
hkl=[0, 0, 1])
ml.reference_plane = midplane
if displace_elements:
# vz needs to be increased by thickness / ( 2 * |vz|)
print('Current vz value : ' + str(self.domain.vz) + ', with multiplier : ' + str(self.domain.bounds[2]))
vz_half_increase = ml.thickness / (4.0 * np.linalg.norm(self.domain.vz))
self.domain.bounds[2] += 2.0 * vz_half_increase
print('Calculated vz increase of : ' + str(2.0 * vz_half_increase) + ' to accomodate layer')
# get box vector for lattice decomposition
a1 = np.array(self.domain.vx)
a2 = np.array(self.domain.vy)
a3 = np.array(self.domain.vz)
lattice_matrix = np.array([[a1[0], a2[0], a3[0]],
[a1[1], a2[1], a3[1]],
[a1[2], a2[2], a3[2]]], dtype=np.float32)
# displace everything else
for bead in self.beads:
fractional_position = list(np.linalg.solve(lattice_matrix, np.array(bead.position)))
if fractional_position[2] >= midplane:
fractional_position[2] += vz_half_increase * 2.0
else:
fractional_position[2] -= vz_half_increase * 2.0
bead.position = np.dot(lattice_matrix, fractional_position)
self.merge(ml)
self.ext_layers.append(ml)
[docs] def enforce_PBC(self, z_box_multi=None):
"""
enforces periodic boundary conditions on particles. The last vector length can be multiplied by a constant to
create empty space
:param z_box_multi: z multiplier for PBC
:type z_box_multi: float
:return:
"""
a1 = np.array(self.domain.vx)
a2 = np.array(self.domain.vy)
a3 = np.array(self.domain.vz)
if z_box_multi is not None:
a3 *= z_box_multi
lattice_matrix = np.array([[a1[0], a2[0], a3[0]],
[a1[1], a2[1], a3[1]],
[a1[2], a2[2], a3[2]]], dtype=np.float32)
for bead in self.beads:
fractional_positions = list(np.linalg.solve(lattice_matrix, np.array(bead.position)))
for i in range(fractional_positions.__len__()):
while fractional_positions[i] > self.domain.bounds[i]:
fractional_positions[i] -= 2*self.domain.bounds[i]
bead.image[i] += 1
while fractional_positions[i] < -self.domain.bounds[i]:
fractional_positions[i] += 2*self.domain.bounds[i]
bead.image[i] -= 1
bead.position = np.dot(lattice_matrix, fractional_positions)
[docs] def enforce_XYPBC(self):
"""
Enforces a special 2D periodic conditions in the XY plane while leaving Z free. (x,y) components of the last
box vectors are assumed to be zero, while the z components of the first two vectors are assumed to be zero. This
is called when making crystal slabs.
:return: None
"""
a1 = np.array(self.domain.vx)
a2 = np.array(self.domain.vy)
xy_matrix = np.array([[a1[0], a2[0]], [a1[1], a2[1]]], dtype=np.float32)
if abs(a1[2]) > 1e-5:
warnings.warn('XY bound Vector (a1) has non-zero z component; check results')
if abs(a2[2]) > 1e-5:
warnings.warn('XY bound Vector (a2) has non-zero z component; check results')
for bead in self.beads:
fractional_position = list(np.linalg.solve(xy_matrix, np.array(bead.position[0:2])))
for i in range(fractional_position.__len__()):
while fractional_position[i] > self.domain.bounds[i]:
fractional_position[i] -= 2 * self.domain.bounds[i]
bead.image[i] += 1
while fractional_position[i] < -self.domain.bounds[i]:
fractional_position[i] += 2 * self.domain.bounds[i]
bead.image[i] -= 1
real_positions = np.dot(xy_matrix, fractional_position)
bead.position[0] = real_positions[0]
bead.position[1] = real_positions[1]
[docs]class Electrostatics(object):
"""
Class dealing with charge normalization in the system and permittivity values.
Args:
beads: list of beads for which charge is managed
units: units to which beads correspond
"""
def __init__(self, beads, units):
self.beads = beads
self.normalized = False
self.units = units
self._normalization_constant = 1.0
# values of permittivity taken from :
# J-W Shen, C. Li, N. F. A. Vegt, C. Peter,
# Transferability of coarse grained potentials : implicit solver models for hydrated ions,
# J. Chem. Th. and Comp., 2011
self.XsaltCPeter = np.array([0.0, 0.4865, 0.6846, 1.0, 2.0, 3.0, 4.0, 5.0], dtype=np.float32)
self.YsaltCPeter = np.array([71.7457, 62.5617, 60.1328, 56.5655, 47.8368, 40.2467, 35.8444, 32.2770],
dtype=np.float32)
@property
def permittivity(self):
"""
Defines the relative permittivity of the system. For atomistic simulations, this value should be 1. For CG
systems, this depends on the force-field. For Martini, this is 10.0.
:getter: Returns the permittivity value
:setter: Sets the permittivity value
:type: float
"""
return (self._normalization_constant * self.units.elem_charge())**2.0
@permittivity.setter
def permittivity(self, value):
if self.normalized:
self.denormalize()
self._normalization_constant = value ** 0.5 * 1.0 / self.units.elem_charge()
@property
def normalization(self):
"""
The normalization constant of charges in the system. This is equivalent to :math:`\epsilon^{1/2} / \mathcal{q}`,
where :math:`\epsilon` is the permittivity and :math:`\mathcal{q}` the charge unit
:getter: returns the normalization constant
:setter: sets the normalization constant
:type: float
"""
return self._normalization_constant
@normalization.setter
def normalization(self, value):
if self.normalized:
self.denormalize()
self._normalization_constant = value
[docs] def normalize(self):
"""
Normalizes the charges in the simulation. Normally, this should be done automatically by builders
:return: None
"""
if not self.normalized:
for bead in self.beads:
bead.charge /= self._normalization_constant
self.normalized = True
[docs] def denormalize(self):
"""
Removes the normalization of charges in the system
:return: None
"""
if self.normalized:
for bead in self.beads:
bead.charge *= self._normalization_constant
self.normalized = False
[docs] def set_eps_to_salt_CPeter(self, salt_concentration):
"""
Sets the permittivity of the simulation to known values for water-NaCl mixtures at a given ion concentration.
These values are taken from: J.-W. Shen, C. Li, N. F. A. Vegt, C. Peter, "Transferability of coarse-grained
potentials: implicit solver models for hydrated ions", J. Chem. Th. and Comp., 2011
:param salt_concentration: Salt concentration in water (M)
:type salt_concentration: float
:return: None
"""
self.c_salt = salt_concentration
# curve source :
# J-W Shen, C. Li, N. F. A. Vegt, C. Peter,
# Transferability of coarse grained potentials : implicit solver models for hydrated ions,
# J. Chem. Th. and Comp., 2011
self.permittivity = np.interp(salt_concentration, self.XsaltCPeter, self.YsaltCPeter)
print('Using values from an article, please cite : J.W. Chan, C. Li, N. F. A. Vegt, C. Peter, Transferability '
'of coarse grained potentials : implicit solver models for hydrated ions, J. Chem. Th. and Comp., 2011')
def set_eps_to_salt_DePablo(self, val, temp):
# epsilon source :
# D. M. Hinckley, J. J. De Pablo
# Coarse-Grained Ions for Nucleic Acid Modeling
# J. Chem. Th. Comp. 2015
self.c_salt = val
self.permittivity = 249.4 - 0.788 * temp + 7.20e-4 * temp**2.0
[docs] def get_type_charge_product(self, typeA, typeB):
"""
returns the charge product of two types of beads in the lists. Will not return correct values if multiple beads
of a given type can have different charges
:param typeA: first beadtype
:type typeA: str
:param typeB: second beadtype
:type typeB: str
:return: charge product
:rtype: float
"""
_qa = 0.0
_qb = 0.0
# define some bools so we dont iterate for no reason
_qaf = False
_qbf = False
for bead in self.beads:
if bead.beadtype == typeA:
_qa = bead.charge * 1.0 if self.normalized else self._normalization_constant
_qaf = True
if bead.beadtype == typeB:
_qb = bead.charge * 1.0 if self.normalized else self._normalization_constant
_qbf = True
if _qaf and _qbf:
break
if not _qaf:
warnings.warn(
'Build : get_type_charge_product : Did not find type ' + typeA +
' in the lists, returned 0 for charge product')
if not _qbf:
warnings.warn(
'Build : get_type_charge_product : Did not find type ' + typeB +
' in the lists, returned 0 for charge product')
return _qa * _qb
[docs]class Bias(object):
"""
Object to manage adding biases to the system in the form of bonds, e.g. for umbrella sampling or thermodynamic
integration
"""
def __init__(self, parent):
self.parent = parent
[docs] def add_einstein_crystal_lattice(self, additional_offsets=None, beadopts=None, bond_name=None, colloid_test=None):
"""
Creates additional beads on each center of the colloids. Default type is 'EC' with a default mass of 1.0. Bonds
are appended between this particle and the colloid center for thermodynamic integration. Bond name is 'EC-bond'
by default. Additional points to append can be supplied by additional_offsets, which are offsets from the center
of the colloid by index.
Conditions to evaluate on the properties of the colloid can be specified by a colloid dict.
:param additional_offsets: additional points where to add beads
:param beadopts: options to be passed to the bead constructor
:param bond_name: name of the bond to append
:param colloid_test: dict of properties that has to be specified
:return: None
"""
# manage default arguments
if beadopts is None:
beadopts = {'beadtype': 'EC', 'mass': 1.0}
if bond_name is None:
bond_name = 'EC-bond'
if colloid_test is None:
colloid_test = {}
self.parent.__types.append(beadopts['beadtype'])
for i in range(self.rigid_body_references.__len__()):
validates = True
# test whether we satisfy requirements of the colloid specifications
try:
for key in colloid_test:
if not getattr(self.rigid_body_references[i], key) == colloid_test[key]:
validates = False
# some property is missing from the colloid
except AttributeError:
validates = False
if not validates:
continue
# append the beads / bonds into the lists
self.parent.beads += CGBead.Bead(position=self.parent.rigid_body_references[i].center_position, **beadopts)
self.parent.p_num.append(self.parent.p_num[-1] + 1)
self.parent.bonds += [bond_name, self.parent.rigid_body_references[i].pnum_offset, self.parent.p_num[-1]]
if additional_offsets is not None:
for offset in additional_offsets:
self.parent.beads += CGBead.Bead(position=self.parent.beads[self.parent.rigid_body_references[i].pnum_offset + offset], **beadopts)
self.parent.p_num.append(self.parent.p_num[-1] + 1)
self.parent.bonds += [bond_name, self.parent.rigid_body_references[i].pnum_offset + offset, self.parent.p_num[-1]]
[docs] def add_bonds_colloids(self, bond_name=None, colloid_test=None):
"""
Add bonds between every center of colloid objects in the domain
:param bond_name: Name of the bond type created
:param colloid_test: Dictionary of properties the colloid must satisfy
:type colloid_test: dict
:return:
"""
if colloid_test is None:
colloid_test = {}
bond_lists = []
if bond_name is None:
bond_name = 'Colloid-bond'
for col_idx in range(self.parent.rigid_body_references.__len__()):
validates = True
try:
for key in colloid_test:
if not getattr(self.parent.rigid_body_references[col_idx], key) == colloid_test[key]:
validates = False
except AttributeError:
validates = False
if validates:
bond_lists.append(col_idx)
for idx_i in range(bond_lists.__len__()):
for idx_j in range(idx_i + 1, bond_lists.__len__()):
self.parent.bonds += [bond_name + '-' + str(idx_i) + '-' + str(idx_j),
self.parent.rigid_body_references[idx_i].pnum_offset,
self.parent.rigid_body_references[idx_j].pnum_offset]
[docs]class FileWriter(object):
"""
Handles file writing from builder objects
This is automatically made by the builder class and methods can be accessed by builder.FileWriter.method
Args:
parent: Parent builder object
"""
def __init__(self, parent):
self.parent = parent
[docs] def export_pdb(self, file_path): # well this thing is ugly but works fine
"""
Writes a very short pdb file
:param file_path: file path to write to
:type file_path: str
:return: None
"""
pdb_file = open(file_path, 'w')
for idx in range(len(self.parent.beads)):
stringarray = [' ' for x in range(80)]
stringarray[0:4] = list('ATOM')
idxstring = str(idx)
stringarray[11 - len(idxstring):11] = list(idxstring)
beadstring = copy.copy(self.parent.beads[idx].beadtype[0:4])
stringarray[16 - len(beadstring):16] = list(beadstring)
if self.parent.beads[idx].residue is not None:
resstring = copy.copy(self.parent.beads[idx].residue[0:3])
else:
resstring = 'UDF'
stringarray[20 - len(resstring):20] = list(resstring)
stringarray[21] = 'A'
stringarray[34:35] = ['-', '3']
stringarray[30:37] = list('%8.3f' % self.parent.beads[idx].position[0])
stringarray[38:45] = list('%8.3f' % self.parent.beads[idx].position[1])
stringarray[46:53] = list('%8.3f' % self.parent.beads[idx].position[2])
stringarray[54:59] = list('%6.3f' % 1.0)
stringarray[60:65] = list('%6.3f' % 1.0)
stringarray[77] = 'H'
if self.parent.beads[idx].charge > 0:
stringarray[78] = '+'
elif self.parent.beads[idx].charge < 0:
stringarray[78] = '-'
else:
stringarray[78] = ' '
pstring = ''.join(stringarray)+'\n'
pdb_file.write(pstring)
[docs]class HOOMDBuilder(Builder):
"""
Building methods specialized for the HOOMD package
Args:
domain (SimulationDomain.Domain): simulation domain used by the builder
shapes (list of GenShape objects): list of shapes to build colloids with
kwargs:
units (SimulationUnits): units used by the builder
"""
def __init__(self, domain, **kwargs):
super(HOOMDBuilder, self).__init__(domain, **kwargs)
self.correct_pnum_body_lists()
[docs] def set_snapshot(self, snapshot, PBC=True):
"""
Sets a given HOOMD snapshots to match the contents of the HOOMDBuilder object
:param snapshot: empty snapshot
:type snapshot: hoomd.data.snapshot
:param PBC: Whether to enforce PBC conditions
:type PBC: bool
:return: None
"""
if not hasattr(snapshot, 'particles') or not hasattr(getattr(snapshot, 'particles'), 'N'):
raise AssertionError('Not a valid hoomd snapshot')
if snapshot.particles.N != self.beads.__len__():
warnings.warn('snapshot size doesnt match with internal list size, resizing', UserWarning)
snapshot.resize(self.beads.__len__())
# enforce PBC & normalize charge
if PBC and 'vertical_slice' in self.domain.flags and self.domain.flags['vertical_slice']:
self.enforce_XYPBC()
elif PBC:
self.enforce_PBC(z_box_multi=1.0)
if len(self.impose_box) == 0:
self.match_hoomd_box_conventions()
self.Electrostatics.normalize()
# list all unique topology type in the system
p_typelist = self.bead_types
bond_typelist = []
angle_typelist = []
dihedral_typelist = []
improper_typelist = []
for bond in self.bonds:
bond_typelist += [bond[0]]
bond_typelist = list(set(bond_typelist))
for angle in self.angles:
angle_typelist += [angle[0]]
angle_typelist = list(set(angle_typelist))
for dihedral in self.dihedrals:
dihedral_typelist += [dihedral[0]]
dihedral_typelist = list(set(dihedral_typelist))
for improper in self.impropers:
improper_typelist += [improper[0]]
improper_typelist = list(set(improper_typelist))
for idx in range(0, len(self.beads)):
for prop in self.beads[idx]:
getattr(snapshot.particles, prop[0])[idx] = prop[1]
snapshot.particles.typeid[idx] = p_typelist.index(self.beads[idx].type)
# set topology in the snapshot
snapshot.bonds.resize(len(self.bonds))
snapshot.bonds.types = bond_typelist
for idx in range(self.bonds.__len__()):
snapshot.bonds.group[idx][0] = self.bonds[idx][1]
snapshot.bonds.group[idx][1] = self.bonds[idx][2]
snapshot.bonds.typeid[idx] = bond_typelist.index(self.bonds[idx][0])
snapshot.angles.resize(len(self.angles))
snapshot.angles.types = angle_typelist
for idx in range(self.angles.__len__()):
snapshot.angles.group[idx][0] = self.angles[idx][1]
snapshot.angles.group[idx][1] = self.angles[idx][2]
snapshot.angles.group[idx][2] = self.angles[idx][3]
snapshot.angles.typeid[idx] = angle_typelist.index(self.angles[idx][0])
snapshot.dihedrals.resize(len(self.dihedrals))
snapshot.dihedrals.types = dihedral_typelist
for idx in range(self.dihedrals.__len__()):
snapshot.dihedrals.group[idx][0] = self.dihedrals[idx][1]
snapshot.dihedrals.group[idx][1] = self.dihedrals[idx][2]
snapshot.dihedrals.group[idx][2] = self.dihedrals[idx][3]
snapshot.dihedrals.group[idx][3] = self.dihedrals[idx][4]
snapshot.dihedrals.typeid[idx] = dihedral_typelist.index(self.dihedrals[idx][0])
snapshot.impropers.resize(len(self.impropers))
snapshot.impropers.types = improper_typelist
for idx in range(self.impropers.__len__()):
snapshot.impropers.group[idx][0] = self.impropers[idx][1]
snapshot.impropers.group[idx][1] = self.impropers[idx][2]
snapshot.impropers.group[idx][2] = self.impropers[idx][3]
snapshot.impropers.group[idx][3] = self.impropers[idx][4]
snapshot.impropers.typeid[idx] = improper_typelist.index(self.impropers[idx][0])
snapshot.constraints.resize(len(self.constraints))
for idx, constraint in enumerate(self.constraints):
snapshot.constraints.group[idx][0] = constraint.topology_tuple[0]
snapshot.constraints.group[idx][1] = constraint.topology_tuple[1]
snapshot.constraints.value[idx] = constraint.distance
[docs] def charge_to_types_per_rigid_body(self, types):
"""
Sets a charge to every bead listed in types. This charge only depends on the rigid body number. This was used
in O'Brien et al., "Exploring the zone of anisotropy and broken symmetries in DNA-mediated nanoparticle
self-assembly", PNAS 2016
:param types: beadtypes to apply charge on
:return: None
"""
for index, ref in enumerate(self.rigid_body_references):
for bead in ref.beads:
if bead.beadtype in types:
bead.charge = index + 1
[docs] def aggregate_rigid_tuples(self):
"""
Provides a method to obtain arguments for constrain.rigid(). This is required for rigid body dynamics
:return: parameters to pass to hoomd
:rtype: list[tuple]
"""
# first create the full list then only keep uniques of these lists (only a single of the three have to be
# different for a particle to be added
# build full lists
_r_c_typelist = []
_r_typelist = []
_r_tuple_list = []
for particle in self.rigid_body_references:
if hasattr(particle, 'relative_positions'):
_r_c_typelist.append(particle.center_type)
_r_typelist.append(particle.body_typelist[1:])
_r_tuple_list.append(particle.relative_positions())
# check non-particle data, inserts of rigid bodies should keep lists ordered
for idx in range(self.beads.__len__()):
# body to add
if (self.beads[idx].body != -1
and self.beads[idx].beadtype not in _r_c_typelist
and self.beads[idx].body == idx):
idx_diff = 1
loc_rc = str(self.beads[idx])
loc_rel_type = []
loc_rel_list = []
try:
# assume that they're not cut by edges
while self.beads[idx + idx_diff].body == self.beads[idx].body:
loc_rel_type += self.beads[idx+idx_diff].beadtype
loc_rel_list += (self.beads[idx+idx_diff].position[0] - self.beads[idx].position[0],
self.beads[idx+idx_diff].position[1] - self.beads[idx].position[1],
self.beads[idx+idx_diff].position[2] - self.beads[idx].position[2])
idx_diff += 1
except IndexError:
pass
_r_c_typelist += [loc_rc]
_r_typelist += [loc_rel_type]
_r_tuple_list += [loc_rel_list]
# at least add the first element of the list
try:
_r_c_out = [_r_c_typelist[0]]
_r_type_out = [_r_typelist[0]]
_r_tuple_out = [_r_tuple_list[0]]
except IndexError:
# some particle is defined as rigid but is empty
warnings.warn('Some colloids have empty lists of positions / types, unable to return rigid structures'
'or all colloids are soft')
return [], [], []
# find all uniques particles
for idx in range(1, _r_c_typelist.__len__()):
# all compare refers to the full comparison with all previously added rigid particles it is set to true if
# the three lists are different to all previously added three lists
_all_compare_bool = True
# compare refers to the current comparison in the already added lists
_compare_bool = False
for idx_compare in range(0, _r_c_out.__len__()):
if _r_c_out[idx_compare] != _r_c_typelist[idx]:
_compare_bool = True
elif _r_type_out[idx_compare] != _r_typelist[idx]:
_compare_bool = True
elif _r_tuple_out[idx_compare] != _r_tuple_list[idx]:
_compare_bool = True
_all_compare_bool &= _compare_bool
if not _all_compare_bool:
break
if _all_compare_bool:
_r_c_out.append(_r_c_typelist[idx])
_r_type_out.append(_r_typelist[idx])
_r_tuple_out.append(_r_tuple_list[idx])
# aggregate into tuples
_out = []
for idx in range(_r_c_out.__len__()):
_out.append((_r_c_out[idx], _r_type_out[idx], _r_tuple_out[idx]))
return _out
[docs] def correct_pnum_body_lists(self):
"""
Fixes the body attribute of beads to their particle tag number to fit HOOMD blue standards
:return: None
"""
for ref in self.rigid_body_references:
ref.body = ref.pnum_offset
[docs] def match_hoomd_box_conventions(self):
"""
Shifts the base vectors used in the build to hoomd standards
:return: None
"""
if not type(self.domain).__name__ == 'Lattice':
return
vx, vy, vz = self.domain.direct_lattice_vectors
_v_a = vx
_t = np.dot(_v_a, _v_a)**0.5
for i in range(_v_a.__len__()):
_v_a[i] /= _t
_v_b = [1, 0, 0]
_v = np.cross(_v_a, _v_b)
_s = (np.dot(_v, _v))**0.5
_c = np.dot(_v_a, _v_b)
_id = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
_vx = [[0, -_v[2], _v[1]], [_v[2], 0, -_v[0]], [-_v[1], _v[0], 0]]
if _s == 0:
_hoomd_mat = _id
else:
_hoomd_mat = _id + _vx + np.linalg.matrix_power(_vx, 2) * (1 - _c) / (_s ** 2)
self.rotate(_hoomd_mat)
[docs] def import_snapshot(self, snapshot):
"""
Imports a hoomd snapshot as a HOOMDBuilder object.
:param snapshot: snapshot to import
:type snapshot: hoomd.data.snapshot
:return: None
"""
warnings.warn('Import function of HOOMDBuilder does not import force-fields', UserWarning)
snap_box = snapshot.box
dom = [[snap_box.Lx, 0.0, 0.0],
[snap_box.Ly * snap_box.xy, snap_box.Ly, 0.0],
[snap_box.Lz * snap_box.xz, snap_box.Ly * snap_box.yz, snap_box.Lz]]
domain = SimulationDomain.EmptyBox(dom, units=self.units)
self.__init__(domain, units=self.units)
typeid = snapshot.particles.types
for p in snapshot.particles:
self.beads += CGBead.Bead(position=p.position,
beadtype=typeid[p.typeid],
charge=p.charge,
mass=p.mass,
moment_inertia=p.moment_inertia,
diameter=p.diameter,
body=p.body,
image=p.image,
quaternion=p.orientation)
try:
self.p_num.append(self.p_num[-1]+1)
except IndexError:
self.p_num = [0]
for bond in snapshot.bonds:
self.bonds += Composite.Bond([bond.type, bond.a, bond.b])
for angle in snapshot.angles:
self.angles += Composite.Angle([angle.type, angle.a, angle.b, angle.c])
for die in snapshot.dihedrals:
self.dihedrals += Composite.Dihedral([die.type, die.a, die.b, die.c, die.d])
for imp in snapshot.impropers:
self.impropers += Composite.Dihedral([imp.type, imp.a, imp.b, imp.c, imp.d])
for cons in snapshot.constraints:
self.constraints += Composite.Constraint([cons.a, cons.b], cons.d)
[docs]class openMMBuilder(Builder):
"""
Building methods specialized for the openMM simulation package
Args:
domain (SimulationDomain.Domain): domain to build
kwargs: units to build this
"""
def __init__(self, domain, **kwargs):
super(openMMBuilder, self).__init__(domain, **kwargs)
[docs] def match_box_conventions(self):
"""
shifts the base vectors used in the build to openMM standards
:return: None
"""
if not type(self.domain).__name__ == 'Lattice':
return
vx, vy, vz = self.domain.direct_lattice_vectors
_v_a = vx
_t = np.dot(_v_a, _v_a)**0.5
for i in range(_v_a.__len__()):
_v_a[i] /= _t
_v_b = [1, 0, 0]
_v = np.cross(_v_a, _v_b)
_s = (np.dot(_v, _v))**0.5
_c = np.dot(_v_a, _v_b)
_id = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
_vx = [[0, -_v[2], _v[1]], [_v[2], 0, -_v[0]], [-_v[1], _v[0], 0]]
if _s == 0:
_mat = _id
else:
_mat = _id + _vx + np.linalg.matrix_power(_vx, 2) * (1 - _c) / (_s ** 2)
self.rotate(_mat)
[docs] def write(self, topology):
"""
Write the openMMBuilder contents to a openMM topology class
:param topology: openMM topology object
:return: None
"""
self.change_units(SimUnits(length='nm', energy='kJ/mol', mass='amu'))
self.match_box_conventions()
v = self.current_box()
topology.setPeriodicBoxVectors([v[0], 0.0, 0.0], [v[1] * v[3], v[1], 0.0], [v[2] * v[4], v[2] * v[5], v[5]])
for composite in self.subunits:
ch = topology.addChain()
prev_res = None
prev_res_name = None
default_residue = None
for bead in composite.beads:
if bead.residue is None and prev_res is None:
default_residue = topology.addResidue(chain=ch, name='UNK')
prev_res = default_residue
elif bead.residue is None and prev_res != default_residue:
default_residue = topology.addResidue(chain=ch, name='UNK')
prev_res = default_residue
elif bead.residue == prev_res_name:
pass
else:
prev_res = topology.addResidue(chain=ch, name=bead.residue)
topology.addAtom(name=str(bead), element=str(bead), residue=prev_res)
n_leftovers = len(self.beads) - self.subunits[-1].p_num[-1] - 1
if n_leftovers > 0:
ch = topology.addChain()
res = topology.addResidue(chain=ch, name='UNK')
for bead_i in range(n_leftovers):
i = self.subunits[-1].p_num[-1] + bead_i
topology.addAtoms(name=str(self.beads[i]), element=str(self.beads[i]), residue=res)
for bond in self.bonds:
topology.addBond(topology.atoms[bond.topology_tuple[0]], topology.atoms[bond.topology_tuple[1]])
[docs]class ESPResSOppBuilder(Builder):
"""
Building methods specialized for ESPResSO++ simulation package. Imports yield mangled names as the simulation package
has no internal string representation for bonded types.
Args:
domain (SimulationDomain.Domain): domain to build
units (SimulationUnits): units to use
"""
def __init__(self, domain, **kwargs):
super(ESPResSOppBuilder, self).__init__(domain, **kwargs)
if self.domain.direct_lattice_vectors[0][1] != 0.0 or self.domain.direct_lattice_vectors[0][2] != 0.0 or \
self.domain.direct_lattice_vectors[1][0] != 0.0 or self.domain.direct_lattice_vectors[1][2] != 0.0 or \
self.domain.direct_lattice_vectors[2][0] != 0.0 or self.domain.direct_lattice_vectors[2][1] != 0.0:
warnings.warn('ESPResSO++ does not handle non-orthorhombic domains', UserWarning)
[docs] def write(self, system, espressopp):
"""
Writes the contents of ESPResSOppBuilder to an espresso++ system. The package itself is an argument to this
method
:param system: espresso++ system to set
:type system: espressopp.system
:param espressopp: espresso++ package
:return: None
"""
storage = system.storage
storage.removeAllParticles()
for bead in self.beads:
storage.addParticles([[self.beads.index(bead), espressopp.Real3D(*bead.position), bead.mass, bead.charge]], 'id', 'pos', 'mass', 'q')
[docs] def pairlist_by_type(self, bondtype):
"""
In espresso++, a pairlist represents a bond interaction type. This returns all bonds corresponding to the
supplied bondtype argument
:param bondtype: builder bondtype representation
:type bondtype: str
:return: bonds
:rtype: list[Composite.Bond]
"""
bonds = []
for bond in self.bonds:
if bond.topology_name == bondtype:
bonds.append(bond.topology_tuple)
return bonds
[docs] def triplelist_by_type(self, angletype):
"""
In espresso++, a triple list represents an angle interaction type. This returns all angles corresponding to the
supplied angletype argument
:param angletype: builder angletype representation
:type angletype: str
:return: angles
:rtype: list[Composite.Angle]
"""
angles = []
for angle in self.angles:
if angle.topology_name == angletype:
angles.append(angle.topology_tuple)
return angles
[docs] def quadruplelist_by_type(self, dihedral_type):
"""
In espresso++, a triple list represents an angle interaction type. This returns all angles corresponding to the
supplied angletype argument
:param dihedral_type: builder dihedral representation
:type dihedral_type: str
:return: dihedrals
:rtype: list[Composite.Dihedral]
"""
dies = []
for die in self.dihedrals:
if die.topology_name == dihedral_type:
dies.append(die.topology_tuple)
return dies
[docs] def read(self, system):
"""
Reads an espresso++ system into the builder object. Bonded force-field names obtained are mangled
:param system: espresso++ system
:type system: espressopp.system
:return: None
"""
dom = SimulationDomain.EmptyBox([system.bc.boxL[0], system.bc.boxL[1], system.bc.boxL[2]])
self.__init__(dom)
try:
i = 0
while True:
p = system.storage.getParticle(i)
self.beads += CGBead.Bead(beadtype=str(p.type),
position=p.position,
charge=p.q,
mass=p.mass)
i += 1
except:
pass
for i in range(system.getNumberOfInteractions()):
inter = system.getInteraction(i)
if hasattr(inter, 'getFixedPairList'):
ibond = list(chain(*inter.getFixedPairList().getBonds()))
ibond = [tpl for s in ibond for tpl in s]
for btuple in ibond:
self.bonds += ['ESPResSOPPBond-'+str(i)] + list(btuple)
if hasattr(inter, 'getFixedTripleList'):
iangle = list(chain(*inter.getFixedTripleList.getAngles()))
for atuple in iangle:
self.angles += ['ESPResSOPPAngle-'+str(i)] + list(atuple)
if hasattr(inter, 'getFixedQuadrupleList'):
idie = list(chain(*inter.getFixedQuadrupleList().getDihedrals()))
for ituple in idie:
self.dihedrals += ['ESPResSOPPDihedral-'+str(i)] + list(ituple)
[docs]class CrossBuilder(Builder):
"""
This class provides functionality for cross MD package functionalities. This inherits from all the classes defined
in MetaBuilder.BuilderTypes and also copies all class methods in another namespace, such that if a class MD_Build
has a method called read, it located in CrossBuilder.read and CrossBuilder.MD_Build.read. This is meant to avoid
overriden methods by external builders added to BuilderTypes.
Args:
domain: Simulation domain passed to Builder
units(default = None): units of simulation. Defaults to reduced units
"""
__metaclass__ = MetaBuilder
def __init__(self, domain, **kwargs):
for cls in MetaBuilder.BuilderTypes:
cls.__init__(self, SimulationDomain.EmptyCube(1.0), **kwargs)
setattr(self, cls.__name__, type(cls.__name__, (object,), {})())
# do some inspection on the methods of cls to copy methods to another namespace, this is meant to avoid
# having compatibility issues if the MetaBuilder class were to be expanded upon by many different people
# overriding each other methods.
attrs = inspect.classify_class_attrs(cls)
for attr in attrs:
if attr.defining_class == cls and attr.kind == 'method':
setattr(getattr(self, cls.__name__), attr.name, functools.partial(attr[3], self))
Builder.__init__(self, domain, **kwargs)