import copy
from hoobas.Units import SimulationUnits
from hoobas import LinearChain
import numpy as np
import warnings
import math
import fractions
from hoobas import Util
from hoobas.Quaternion import Quat
from hoobas import Composite
from hoobas import SimulationDomain
[docs]class Multilayer(Composite.CompositeObject):
"""
Makes multilayered structures of :py:class:`hoobas.Composite.CompositeObject`
Args:
NLayers: how many layers should be built
alternating_directions (bool): whether objects should be reversed between layers
"""
def __init__(self, NLayers=1, units=None, alternating_directions=True):
super(Multilayer, self).__init__(units=units)
self.midplane = 0.0
self.NLayers = NLayers
self.squish_factor = 0.6 # the typical lipid bilayer is highly squished: tails of opposing sides are touching
# a function of density can be passed number of chains to build in the layer; since we only have a pointer to a
# function, we need to keep track of surface density units
self.surface_conversion_factor = 1.0
# list of every type of thing to build in the layer
self.species = []
# flags for building
self.is_built = False
self.enforced_thickness = False
# thickness of a single layer
self.layer_thickness = 0
# check if layers are alternating
self.altern_orientations = alternating_directions
# XY transverse scale and scaled coordinates
self.xscale = np.array([1., 0., 0.])
self.yscale = np.array([0., 1., 0.])
self.zscale = np.array([0., 0., 1.])
self.scaled_coordinates = np.zeros((0, 3), dtype=np.float)
@property
def squish(self):
return self.squish_factor
@squish.setter
def squish(self, squish):
self.squish_factor = squish
@property
def thickness(self):
return self.layer_thickness * self.NLayers * self.squish_factor
@property
def reference_plane(self):
return self.midplane
@reference_plane.setter
def reference_plane(self, value):
diff = (value - self.midplane) * self.zscale
for bead in self.beads:
bead.position += diff
self.midplane = value
[docs] def add_species(self, species_type, species_per_layer, squish=None):
"""
Add a species of :py:class:`hoobas.Composite.CompositeObjects` to the layer
:param species_type: Object to add
:type species_type: Composite.CompositeObject
:param species_per_layer: Number of objects on each layer
:type species_per_layer: int, callable
:param squish: relative depth of the object
:type squish: float
:return: None
"""
# check if we can figure a thickness to the species
species_type.change_units(self.units)
if squish is None:
squish = 1.0
try:
thickness = species_type.thickness
except AttributeError:
if issubclass(species_type.__class__, LinearChain.LinearChain):
thickness = species_type.nmono * species_type.lmono
else:
thickness = 1.0
warnings.warn('Could not find thickness of object of type L ' + str(species_type.__class__)
+ ', assuming a thickness of 1.0 in units supplied to multilayer class')
# set layer thickness to the largest thickness of all composing ojects
if not self.enforced_thickness and thickness > self.layer_thickness:
self.layer_thickness = thickness
self.species.append({'object': copy.deepcopy(species_type),
'number': species_per_layer,
'thickness': thickness,
'squish': squish})
if hasattr(species_per_layer, '__call__'): # if a density is set
self.species[-1]['density function'] = species_per_layer
def change_units(self, new_units):
# convert every chain of the layer to new units
for specie in self.species:
if hasattr(specie['object'], '__iter__'):
for sub_specie in specie['object']:
sub_specie.change_units(self.units)
else:
specie['object'].change_units(self.units)
# length conversion factor (this is a purely geometrical construct)
l_convert_factor = new_units.get_length(1.0, self.units.lunit)
# units of length
self.layer_thickness *= l_convert_factor
for specie in self.species:
if hasattr(specie['thickness'], '__iter__'):
for sub_th in specie['thickness']:
sub_th *= l_convert_factor
else:
specie['thickness'] *= l_convert_factor
super(Multilayer, self).change_units(new_units)
# units of number / area
self.surface_conversion_factor /= l_convert_factor * l_convert_factor
self.units = new_units
def build(self):
if self.is_built:
warnings.warn('Trying to build an already built layer of type ' + str(self.__class__.__name__)
+ '; check results')
self.is_built = True
# calculate position of the midplanes
if self.NLayers % 2 == 0:
midplanes = np.arange(-0.5 - (self.NLayers/2 - 1), 0.5 + (self.NLayers/2 - 1) + 0.02)
else:
midplanes = np.arange(-self.NLayers/2, self.NLayers/2 + 0.02)
midplanes *= self.layer_thickness * self.squish_factor
for layer_index in range(self.NLayers): # for every layer
# check if layer needs a vertical flip
if layer_index % 2 == 1 and self.altern_orientations:
flipv = True
else:
flipv = False
total_objects = 0
for species_dict in self.species:
total_objects += species_dict['number']
# make a grid that contains all numbers and add offset it by (0.25, 0.25) * distance to avoid overlaps
grid_edge_length = int(math.ceil(float(total_objects)**0.5))
xg, yg = np.meshgrid(np.linspace(-0.5, 0.5, grid_edge_length, endpoint=False),
np.linspace(-0.5, 0.5, grid_edge_length, endpoint=False))
xg = xg.flatten()
yg = yg.flatten()
offsetx = np.random.uniform(-0.5, 0.5, len(xg)) * 0.50 / grid_edge_length
offsety = np.random.uniform(-0.5, 0.5, len(yg)) * 0.50 / grid_edge_length
if layer_index % 2 == 1: # add an offset to the grid half a grid size and fold back into first unit cell
offsetx += 0.5 / grid_edge_length
offsety += 0.5 / grid_edge_length
xg += offsetx
yg += offsety
div = np.ones(len(xg)) * 0.5
xg = xg - np.modf(xg / div)[1]
yg = yg - np.modf(yg / div)[1]
xg = list(xg)
yg = list(yg)
for species_dict in self.species: # for every species in that layer
species_number = species_dict['number']
for NObjects in range(species_number): # build the number of object of that species in that layer
# old pure random mode
# xp = np.random.uniform(-0.5, 0.5)
# yp = np.random.uniform(-0.5, 0.5)
cint = np.random.randint(0, len(xg))
xp = xg.pop(cint)
yp = yg.pop(cint)
c_obj = copy.deepcopy(species_dict['object'])
self.merge(c_obj)
squished = midplanes[layer_index] * species_dict['squish']
# flip the chain if required
if flipv:
c_obj.rotate(Quat([0., 1., 0., 0.]))
if hasattr(c_obj, 'center_position'):
c_obj.center_position = np.array([0.0, 0.0, squished])
elif hasattr(c_obj, 'position'):
c_obj.position = np.array([0.0, 0.0, squished])
else:
raise SyntaxError('Supplied object of class ' + str(c_obj.__class__)
+ ' doesn''t have a method to move positions')
setattr(c_obj, 'scaled_coordinates', np.array([xp, yp, 0.0]))
def _recalculate_numbers_from_area(self, area):
for specie in self.species:
if hasattr(specie['object'], '__iter__'):
for index in range(len(specie['object'])):
if hasattr(specie['density function'][index], '__call__'):
specie['number'][index] = int(specie['density function'][index](area) * self.surface_conversion_factor)
else:
if 'density function' in specie:
specie['number'] = int(specie['density function'](area) * self.surface_conversion_factor)
@Composite.deferred
def __transform_map(self, lattice, position_map=None, normal_map=None):
for item in self.subunits:
local_position = np.array(position_map(item.scaled_coordinates[0], item.scaled_coordinates[1]))
local_normal = np.array(normal_map(item.scaled_coordinates[0], item.scaled_coordinates[1]))
local_normal = np.dot(local_normal, lattice)
metric = 1.0 / np.linalg.norm(local_normal)
rotation = Util.get_rotation_matrix(local_normal * metric)
item.rotate(rotation)
item.positions += local_position
return
@Composite.deferred
def __transform_hkl(self, lattice, hkl=None, shift=None, inplane_vector_one=None, inplane_vector_two=None):
real_normal = np.cross(inplane_vector_one, inplane_vector_two)
unit_normal = real_normal / np.linalg.norm(real_normal)
rotation = Util.get_rotation_matrix(unit_normal)
# determine the number of planes in the unit cell. A plane like (111) looks like it has 2 planes but the second
# is actually a fold of the first one. Combined together, the two planes of (111) form a parallologram
nplanes = 1
for i in hkl:
if i != 0:
nplanes = nplanes * i // math.gcd(i, nplanes)
matrix_transform = np.transpose([inplane_vector_one, inplane_vector_two, unit_normal])
for item in self.subunits:
item.rotate(rotation)
item.positions += np.dot(matrix_transform, item.scaled_coordinates)
dhkl = SimulationDomain.Lattice.d_hkl(lattice, hkl)
# duplicate the plane nplanes -1 times, shifting it by dhkl everytime
single_plane_sub = copy.deepcopy(self.subunits)
for iplane in range(nplanes - 1):
for item in single_plane_sub:
self.subunits.append(copy.deepcopy(item))
item.positions += dhkl * unit_normal
if shift is not None:
for item in self.subunits:
item.positions += shift * dhkl * unit_normal
[docs] def force_layer_thickness(self, thickness):
"""
Forces the individual layers to have a specific thickness
:param thickness: thickness of the layers
:type thickness: float
:return: None
"""
self.enforced_thickness = True
self.layer_thickness = thickness
[docs]class LayeredTiling(Multilayer):
"""
class meant to emulate liquid-liquid phase separation by pre-separating components. The layers will be N x M tiles
of different chains. The chains will be fully phase separated; if a ternary mixture A-B-C needs to be built in tiles
of A-B + C, this can be made by using a list of component and numbers in the species
Args:
nx (int): number of tiles in X direction
ny (int): number of tiles in Y direction
alternating_directions (bool): whether to reverse objects between layers
NLayers (int): number of layers to build
"""
def __init__(self, nx=1, ny=1, **kwargs):
super(LayeredTiling, self).__init__(**kwargs)
# strides in x, y and total number of patches per layer
self.x_period = nx
self.y_period = ny
self.N_patches = nx * ny
self.layer_number_multiplier = 1.0 / self.N_patches
[docs] def add_species(self, species_type, species_per_layer, squish=None):
"""
Adds :py:class:`hoobas.Composite.CompositeObject` species to the layers.
:param species_type: type of species to add
:type species_type: Composite.CompositeObject, list[Composite.CompositeObject]
:param species_per_layer: number of objects per layer
:type species_per_layer: int, list[int]
:param squish: how much individual object types should be squished
:type squish: float, list[float]
:return: None
"""
# check if we can figure a thickness to the species
if hasattr(species_type, '__iter__'):
dict_to_append = {'object': [], 'number': [], 'thickness': [], 'density function': [], 'squish': []}
for specie in species_type:
specie.change_units(self.units)
try:
thickness = specie.thickness
except AttributeError:
if issubclass(specie.__class__, LinearChain.LinearChain):
thickness = specie.nmono * specie.lmono
else:
thickness = 1.0
warnings.warn('Could not find thickness of object of type L ' + str(specie.__class__)
+ ', assuming a thickness of 1.0 in units supplied to multilayer class')
if not self.enforced_thickness and thickness > self.layer_thickness:
self.layer_thickness = thickness
dict_to_append['object'].append(copy.deepcopy(specie))
dict_to_append['thickness'].append(thickness)
if squish is None:
dict_to_append['squish'].append(1.0)
else:
dict_to_append['squish'].append(squish[species_type.index(specie)])
if hasattr(species_per_layer, '__iter__') and len(species_per_layer) != len(species_type):
raise SyntaxError('Number of objects per layer must either be a scalar (function or number) or '
'an array of length equal to the number of types')
# check whether we're passed a scalar or array of numbers and parse whether we have functions passed in
if hasattr(species_per_layer, '__iter__'):
for num in species_per_layer:
dict_to_append['number'].append(num)
if hasattr(num, '__call__'):
dict_to_append['density function'].append(num)
else:
dict_to_append['density function'].append([])
else: # extend the numbers to every type of object
dict_to_append['number'] = [species_per_layer] * len(species_type)
if hasattr(species_per_layer, '__call__'):
dict_to_append['density function'] = [species_per_layer] * len(species_type)
else:
dict_to_append['density function'] = [[]] * len(species_type)
self.species.append(dict_to_append) # append to internal dictionaries
else:
species_type.change_units(self.units)
if squish is None:
squish = 1.0
try:
thickness = species_type.thickness
except AttributeError:
if issubclass(species_type.__class__, LinearChain.LinearChain):
thickness = species_type.nmono * species_type.lmono
else:
thickness = 1.0
warnings.warn('Could not find thickness of object of type L ' + str(species_type.__class__)
+ ', assuming a thickness of 1.0 in units supplied to multilayer class')
# set layer thickness to the largest thickness of all composing ojects
if not self.enforced_thickness and thickness > self.layer_thickness:
self.layer_thickness = thickness
self.species.append({'object': copy.deepcopy(species_type),
'number': species_per_layer,
'thickness': thickness,
'squish': squish})
if hasattr(species_per_layer, '__call__'): # if a density is set
self.species[-1]['density function'] = species_per_layer
def build(self):
if self.is_built:
warnings.warn('Trying to build an already built layer of type ' + str(self.__class__.__name__)
+ '; check results')
self.is_built = True
if self.x_period % len(self.species):
warnings.warn('Number of periods in x non-commensurate with number of species; check results', UserWarning)
if self.y_period % len(self.species):
warnings.warn('Number of periods in y non-commensurate with number of species; check results', UserWarning)
# calculate position of the midplanes
if self.NLayers % 2 == 0:
midplanes = np.arange(-0.5 - (self.NLayers/2 - 1), 0.5 + (self.NLayers/2 - 1) + 0.02)
else:
midplanes = np.arange(-self.NLayers/2, self.NLayers/2 + 0.02)
midplanes *= self.layer_thickness * self.squish_factor
for layer_index in range(self.NLayers): # for every layer
# check if layer needs a vertical flip
if layer_index % 2 == 1 and self.altern_orientations:
flipv = True
else:
flipv = False
# scaled coordinates of the current patch on the surface
xcoords = np.linspace(-0.5, 0.5, num=self.x_period+1, endpoint=True)
ycoords = np.linspace(-0.5, 0.5, num=self.y_period+1, endpoint=True)
for patch_index in range(self.N_patches):
(xidx, yidx) = np.unravel_index(patch_index, (self.x_period, self.y_period))
if len(self.species) > self.x_period:
species_dict = self.species[patch_index % len(self.species)]
else:
species_dict = self.species[(xidx + yidx) % len(self.species)] # one specie per patch; use cyclic index
if hasattr(species_dict['object'], '__iter__'):
specie_range = len(species_dict['object'])
else:
specie_range = 1
# calculate the total number of objects
patch_number_objects = 0
for specie_index in range(specie_range):
try:
patch_number_objects += species_dict['number'][specie_index] / self.N_patches
except TypeError:
patch_number_objects += species_dict['number'] / self.N_patches
# make a grid over the patch, add a random displament vector in range of contact circle and fold back
grid_edge_length = int(math.ceil(float(patch_number_objects) ** 0.5))
xg = np.linspace(xcoords[xidx], xcoords[xidx + 1], grid_edge_length+2)
xg = xg[1:-1]
yg = np.linspace(ycoords[yidx], ycoords[yidx + 1], grid_edge_length + 2)
yg = yg[1:-1]
xg, yg = np.meshgrid(xg, yg)
xg = xg.flatten()
yg = yg.flatten()
offsetx = np.random.uniform(-0.25, 0.25, len(xg)) * 0.5 / grid_edge_length
offsety = np.random.uniform(-0.25, 0.25, len(yg)) * 0.5 / grid_edge_length
if layer_index % 2 == 1: # add an offset to the grid half a grid size and fold into first unit cell
offsetx += 0.5 / grid_edge_length
offsety += 0.5 / grid_edge_length
xg += offsetx
yg += offsety
div = np.ones(len(xg)) * 0.5
xg = xg - np.modf(xg / div)[1]
yg = yg - np.modf(yg / div)[1]
xg = list(xg)
yg = list(yg)
for specie_index in range(specie_range):
try:
local_specie_number = species_dict['number'][specie_index] / self.N_patches # each patch has a lower number of chains
except TypeError:
local_specie_number = species_dict['number'] / self.N_patches
for NObjects in range(int(local_specie_number)):
cint = np.random.randint(0, len(xg))
xp = xg.pop(cint)
yp = yg.pop(cint)
try:
c_obj = copy.deepcopy(species_dict['object'][specie_index])
except TypeError:
c_obj = copy.deepcopy(species_dict['object'])
# flip the chain if required
if flipv:
c_obj.rotate(Quat([0., 1., 0., 0.]))
self.merge(c_obj)
if hasattr(c_obj, 'center_position'):
c_obj.center_position = np.array([0.0, 0.0, midplanes[layer_index]])
elif hasattr(c_obj, 'position'):
c_obj.position = np.array([0.0, 0.0, midplanes[layer_index]])
else:
raise SyntaxError('Supplied object of class ' + str(c_obj.__class__)
+ ' doesn''t have a method to move positions')
setattr(c_obj, 'scaled_coordinates', np.array([xp, yp, 0.0]))