#coding: utf-8
"""
Genshape provides classes to represent colloid shapes in hoobas.
"""
from math import *
from fractions import gcd
import warnings
import numpy as np
import copy
from hoobas import Composite
from hoobas.Util import get_rotation_matrix
from hoobas.Util import get_inverse_rotation_matrix
from hoobas.Quaternion import Quat
from hoobas.Units import SimulationUnits as SimUnits
[docs]class shape(object):
"""
Generic object describing the shape of a colloid in hoobas. The object has a shape.table array which contains
points representing the shape. New shapes can easily be added by adding new objects inheriting the shape class
self.flags is a dictionary for building directives, which can be set by using shape.set_properties
properties is a dict object for additional properties useful for non-trivial cases (read xyz / pdb parsed files):
'normalized' : set to either True or False, defines if the shape is a geometric representation (eg cube) or a finite
size object (eg protein)
'surface' : used by normalized shapes to calculate surface of an object
'density' : used by normalized shapes to calculate mass
'volume' : used by normalized shapes to calculate volume properties (eg mass)
'mass' : mass of the object, mainly used for proteins, overrides the density*volume value
'size' : scaling factor for normalized shapes
Constructor :
:param surf_plane : rotation to give every object to match given z axis crystal axis
:param lattice : used to specify crystal lattice in conjuction with surf_plane
:param properties : dictionary used to set flags
:param units: units to be used by the shape.
General methods :
set_properties : sets the additional properties required for building such as density, mass, etc.
Other internals :
Bonded types :
self.internal_bonds is a list of length N, containing N bonds in the form of [i, j, r0, type]
self.internal_angles is a list of length N, containing N angle potentials in the form [i, j, k, theta0, type]
Rigid body characteristics :
self.Itensor is the inertia tensor of the shape, in the form of 1 x 3 values for principal axis
self.quaternion is the quaternion of the shape, relating the supplied table to the inertia tensor principal axis
"""
def __init__(self, properties=None, units=None):
self.table = np.zeros((0, 3), dtype=np.float32)
self.BuildMethod = lambda: None
# additional DOF for moment corrections
self.additional_points = np.zeros((1, 3), dtype=np.float32) # add the 0, 0, 0 position
self.masses = np.zeros(0, dtype=np.float32)
self.type_suffix = ['']
# internal DOF for soft shells
self.internal_bonds = []
self.internal_angles = []
self.grafts = []
# set shape units
if units is not None and issubclass(units.__class__, SimUnits):
self.units = units
else:
self.units = SimUnits()
if properties is None:
self.flags = {}
else:
self.flags = properties
self.flags['hard_core_safe_dist'] = 0
self.keys = {}
self._pdb = []
# keying atomic symbols to masses (in amu)
self.mkeys = {'C': 12.011, 'O': 15.999, 'H': 1.008, 'N': 14.007, 'S': 32.06}
self.quaternion = Quat(np.eye(3))
self.Itensor = (0.0, 0.0, 0.0)
@property
def num_surf(self):
return self.table.__len__()
@property
def volume(self):
try:
return self.flags['volume']
except KeyError:
return None
@property
def pos(self):
return self.table
@property
def unit_shape(self):
try:
return self.flags['normalized']
except KeyError:
return None
def set_properties(self, properties=None):
if not (properties is None):
self.flags.update(properties)
self.flags['pdb_object'] = True
def remove_duplicates(self, tol=1e-4):
DupRows = []
for i in range(self.table.__len__() - 1):
for j in range(i + 1, self.table.__len__()):
if (self.table[i, 0] - self.table[j, 0]) ** 2 + (self.table[i, 1] - self.table[j, 1]) ** 2 + (
self.table[i, 2] - self.table[j, 2]) ** 2 < tol * tol:
DupRows.append(j)
DupRows = np.unique(DupRows)
self.table = np.delete(self.table, DupRows, axis=0)
def add_grafts(self, graft, number=1, connecting_topology=None, qualifier=None):
if not issubclass(graft.__class__, Composite.CompositeObject):
raise SyntaxError('Cannot graft objects which are not of type CompositeObject. Perform an import first')
self.grafts.append({'object': graft,
'number': number,
'connecting_topology': connecting_topology,
'qualifier': qualifier})
def rotate_object(self, operator):
if operator is None:
operator = np.eye(3)
_qop = Quat(operator)
self.quaternion = _qop * self.quaternion
[docs] def set_geometric_quaternion(self, tol_eps_rel=2.0):
"""
calculates the tensor of the structure in self.table; uses diagonalization methods and has numerical precision
issues
:return:
"""
inertia_tensor = np.zeros((3, 3), dtype=np.longfloat)
for point in self.table:
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]
#####################################
# numerical fixing of some low values
#####################################
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]) < tol_eps_rel * np.finfo(np.float32).eps * _maxel:
inertia_tensor[dim1, dim2] = 0.0
w, v = np.linalg.eigh(inertia_tensor * 0.6 / self.table.__len__())
# 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]
self.quaternion = Quat(v)
self.Itensor = np.array([w[0], w[1], w[2]], dtype=np.float32)
self.quat_to_hoomd()
[docs] def get_N_idx(self, N, **sortargs):
"""
gets the N smallest indexes of norm values in self.table calculates offsets so whatever is defined in
self.additional_points is not taken into account.
This is made for geometric shapes, but can be used for other things.
:param N number of largest norm values to return
:return list of indexes
"""
_norm = np.zeros(self.table.__len__(), dtype=np.float32)
for idx in range(self.table.__len__()):
_norm[idx] = self.table[idx, 0] ** 2.0 + self.table[idx, 1] ** 2.0 + self.table[idx, 2] ** 2.0
args = np.argsort(_norm, **sortargs)
args = args[0: N] + self.additional_points.__len__()
return args.tolist()
def change_units(self, new_units):
# change coordinates, moments of inertia, properties of bead lists, mass Vector,
if new_units is None:
return
if not issubclass(new_units.__class__, SimUnits):
raise SyntaxError('Units supplied to GenShape instance is not a subclass of SimulationUnits')
if new_units.lunit != self.units.lunit:
m_len = new_units.get_length(1.0, self.units.lunit)
if 'normalized' in self.flags and self.flags['normalized'] is False:
self.table *= m_len
if 'size' in self.flags:
self.flags['size'] *= m_len
if 'density' in self.flags:
self.flags['density'] /= m_len * m_len * m_len
if 'I_tensor' in self.flags:
self.flags['I_tensor'] *= m_len * m_len
self.Itensor *= m_len * m_len
else:
m_len = 1.0
if new_units.munit != self.units.munit:
m_mass = new_units.get_mass(1.0, self.units.munit)
if 'normalized' in self.flags and self.flags['normalized'] is False:
self.Itensor *= m_mass
if 'density' in self.flags:
self.flags['density'] *= m_mass
if 'mass' in self.flags:
self.flags['mass'] *= m_mass
else:
m_mass = 1.0
if new_units.Eunit != self.units.Eunit:
E_mul = new_units.get_energy(1.0, self.units.Eunit)
else:
E_mul = 1.0
self.units = new_units
def quat_to_hoomd(self):
# this sets the table to match hoomd internals
for idx in range(len(self.table)):
lv = np.array(copy.deepcopy(self.table[idx]))
lv = self.quaternion.inv().transform.dot(lv)
self.table[idx] = lv
def surface_normal(self, v):
return np.array([v[0], v[1], v[2]]) / np.linalg.norm(np.array([v[0], v[1], v[2]]))
[docs]class Cube(shape):
def __init__(self, Num, properties=None, Radius=None):
shape.__init__(self, properties)
if Radius is None:
Radius = 0.0
NumSide = int(round((Num * (1-Radius)**2/6)**0.5))
#Create 6 faces, with range (-1+Radius) to (1-Radius), NumSide**2 points on each
Gridbase = np.linspace(-1+Radius, 1-Radius, NumSide)
xGrid, yGrid = np.meshgrid(Gridbase,Gridbase)
Tablexy = np.zeros((0,3))
for i in range(NumSide):
for j in range(NumSide):
Tablexy = np.append(Tablexy, np.array([[xGrid[i,j], yGrid[i,j],1]]),axis=0)
# numside is the proportionality constant, pi * R /4 / (2*(1-R)) is the ratio of the side length to quarter cylinder length
PtEdge = NumSide**2 * pi * Radius / (4*2*(1-Radius))
PtAngle = int(PtEdge/NumSide)
#PtL = int((PtEdge*4*(1-Radius)/pi)**0.5)
AngleRange = np.linspace(0,pi/2,PtAngle+2)
TableEdge = np.zeros((0,3))
#XY-YZ edge
for i in range(PtAngle):
for j in range(NumSide):
TableEdge = np.append(TableEdge, np.array([[(1-Radius)+Radius*cos(AngleRange[i+1]),Gridbase[j], (1-Radius)+Radius*sin(AngleRange[i+1])]]),axis =0)
#XYZ vertice
PtVertice = NumSide**2 *(4*pi*Radius**2/8) / (2*(1-Radius))**2
TableVert = np.zeros((0,3))
PtVertice = int(PtVertice*8)
gold_ang = pi*(3-5**0.5)
th = gold_ang*np.arange(PtVertice)
if PtVertice>1:
z = np.linspace(1-1.0/PtVertice, 1.0/PtVertice -1, PtVertice)
for i in range(0, PtVertice):
if 0 < cos(th[i]) and sin(th[i]) > 0 and z[i] > 0:
TableVert = np.append(TableVert, np.array([[(1-Radius)+Radius*(1-z[i]**2)**0.5*cos(th[i]), (1-Radius)+Radius*(1-z[i]**2)**0.5*sin(th[i]),(1-Radius)+Radius*z[i]]]),axis=0)
elif PtVertice == 1 :
TableVert = np.append(TableVert, np.array([[(1-Radius)+Radius*3**0.5, (1-Radius)+Radius*3**0.5,(1-Radius)+Radius*3**0.5]]),axis=0)
#Write 6 faces
for i in range(Tablexy.__len__()):
#XY +/-
self.table = np.append(self.table, [[ Tablexy[i,0], Tablexy[i,1], Tablexy[i,2]]], axis=0)
self.table = np.append(self.table, [[ Tablexy[i,0], Tablexy[i,1],-Tablexy[i,2]]], axis=0)
#YZ +/-
self.table = np.append(self.table, [[ Tablexy[i,2], Tablexy[i,0], Tablexy[i,1]]], axis=0)
self.table = np.append(self.table, [[-Tablexy[i,2], Tablexy[i,0], Tablexy[i,1]]], axis=0)
#ZX
self.table = np.append(self.table, [[ Tablexy[i,1], Tablexy[i,2], Tablexy[i,0]]], axis=0)
self.table = np.append(self.table, [[ Tablexy[i,1],-Tablexy[i,2], Tablexy[i,0]]], axis=0)
# Write 12 Edges
for i in range(TableEdge.__len__()):
# Y kind
self.table = np.append(self.table, [[ TableEdge[i, 0], TableEdge[i, 1], TableEdge[i, 2]]], axis=0)
self.table = np.append(self.table, [[-TableEdge[i, 0], TableEdge[i, 1], TableEdge[i, 2]]], axis=0)
self.table = np.append(self.table, [[-TableEdge[i, 0], TableEdge[i, 1],-TableEdge[i, 2]]], axis=0)
self.table = np.append(self.table, [[ TableEdge[i, 0], TableEdge[i, 1],-TableEdge[i, 2]]], axis=0)
self.table = np.append(self.table, [[ TableEdge[i, 2], TableEdge[i, 0], TableEdge[i, 1]]], axis=0)
self.table = np.append(self.table, [[-TableEdge[i, 2], TableEdge[i, 0], TableEdge[i, 1]]], axis=0)
self.table = np.append(self.table, [[-TableEdge[i, 2],-TableEdge[i, 0], TableEdge[i, 1]]], axis=0)
self.table = np.append(self.table, [[ TableEdge[i, 2],-TableEdge[i, 0], TableEdge[i, 1]]], axis=0)
self.table = np.append(self.table, [[ TableEdge[i, 1], TableEdge[i, 2], TableEdge[i, 0]]], axis=0)
self.table = np.append(self.table, [[ TableEdge[i, 1],-TableEdge[i, 2], TableEdge[i, 0]]], axis=0)
self.table = np.append(self.table, [[ TableEdge[i, 1],-TableEdge[i, 2],-TableEdge[i, 0]]], axis=0)
self.table = np.append(self.table, [[ TableEdge[i, 1], TableEdge[i, 2],-TableEdge[i, 0]]], axis=0)
# Write 8 vertices
for i in range(TableVert.__len__()):
self.table = np.append(self.table, [[ TableVert[i, 0], TableVert[i, 1], TableVert[i, 2]]], axis=0)
self.table = np.append(self.table, [[-TableVert[i, 0], TableVert[i, 1], TableVert[i, 2]]], axis=0)
self.table = np.append(self.table, [[ TableVert[i, 0],-TableVert[i, 1], TableVert[i, 2]]], axis=0)
self.table = np.append(self.table, [[ TableVert[i, 0], TableVert[i, 1],-TableVert[i, 2]]], axis=0)
self.table = np.append(self.table, [[ TableVert[i, 0],-TableVert[i, 1],-TableVert[i, 2]]], axis=0)
self.table = np.append(self.table, [[-TableVert[i, 0], TableVert[i, 1],-TableVert[i, 2]]], axis=0)
self.table = np.append(self.table, [[-TableVert[i, 0],-TableVert[i, 1], TableVert[i, 2]]], axis=0)
self.table = np.append(self.table, [[-TableVert[i, 0],-TableVert[i, 1],-TableVert[i, 2]]], axis=0)
self.remove_duplicates()
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = 3**0.5
self.flags['volume'] = 2**3
self.flags['surface'] = 6*2**2
self.flags['simple_I_tensor'] = True
self.flags['call'] = 'cube'
self.set_geometric_quaternion()
[docs]class Octahedron(shape):
def __init__(self, Num, properties=None):
shape.__init__(self, properties)
#Creating faces; Generate triangular lattice with N points on the lower side. Triangular points are (1,sqrt(3)),(0,0), (2, 0), Limit points are
FacePoint = Num / 8
WhP = 1
Table = np.array([[0, 0]])
TableAdd = np.zeros((1, 2))
tol = 1e-5
while WhP < FacePoint:
TableAdd = Table + [[1, 0]]
TableAdd = np.append(TableAdd,Table+ [[0.5, 3**0.5/2]], axis = 0)
Table = np.append(Table, TableAdd, axis = 0)
DupRows = []
for i in range(Table.__len__()-1):
for j in range(i+1, Table.__len__()):
if (Table[i,0]-Table[j,0])**2 +(Table[i,1]-Table[j,1])**2 < tol:
DupRows.append(j)
DupRows = np.unique(DupRows)
Table = np.delete(Table, DupRows, axis = 0)
WhP = Table.__len__()
FacePoint = Table.__len__()
delta =0
N = max(Table[:,0])+1
q = (2- 2*3**0.5*delta)/ (N-1)
Table = (((Table*q) + [[delta*3**0.5, delta]])-[[1, 3**0.5 * 1/3.0]])*0.5
TableFace = np.zeros((Table.__len__(), 3))
FaAng = (pi - acos(1/3.0))*0.5
for i in range(TableFace.__len__()):
TableFace[i,:] = np.array([[Table[i,0], (Table[i,1])*cos(FaAng),sin(FaAng)*(Table[i,1])]])
TableFace = TableFace + [[0, -3**0.5 / 3.0 *cos(FaAng), 3**0.5 / 6.0 *sin(FaAng)]]
TotalTable = np.zeros((0,3))
#Append 8 faces by symmetry
for i in range(TableFace.__len__()):
TotalTable = np.append(TotalTable, np.array([[TableFace[i,0], TableFace[i,1], TableFace[i,2]]]), axis = 0)
TotalTable = np.append(TotalTable, np.array([[TableFace[i,0], -TableFace[i,1], TableFace[i,2]]]), axis = 0)
TotalTable = np.append(TotalTable, np.array([[TableFace[i,0], -TableFace[i,1], -TableFace[i,2]]]), axis = 0)
TotalTable = np.append(TotalTable, np.array([[TableFace[i,0], TableFace[i,1], -TableFace[i,2]]]), axis = 0)
TotalTable = np.append(TotalTable, np.array([[TableFace[i,1], TableFace[i,0], TableFace[i,2]]]), axis = 0)
TotalTable = np.append(TotalTable, np.array([[-TableFace[i,1], TableFace[i,0], TableFace[i,2]]]), axis = 0)
TotalTable = np.append(TotalTable, np.array([[-TableFace[i,1], TableFace[i,0], -TableFace[i,2]]]), axis = 0)
TotalTable = np.append(TotalTable, np.array([[TableFace[i,1], TableFace[i,0], -TableFace[i,2]]]), axis = 0)
DupRows = []
for i in range(TotalTable.__len__()):
for j in range(i+1, TotalTable.__len__()):
if (TotalTable[i,0]-TotalTable[j,0])**2 +(TotalTable[i,1]-TotalTable[j,1])**2 +(TotalTable[i,2]-TotalTable[j,2])**2 < tol:
DupRows.append(j)
DupRows = np.unique(DupRows)
TotalTable = np.delete(TotalTable, DupRows, axis=0)
self.table = 2*TotalTable
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = 3**0.5
self.flags['volume'] = 2**0.5 * 2**3 / 3
self.flags['surface'] = 8*3**0.5*2**2/4
self.flags['simple_I_tensor'] = True
self.flags['call'] = 'oct'
self.set_geometric_quaternion()
[docs]class ShellShape(shape):
def __init__(self, properties=None, units=None):
super(ShellShape, self).__init__(properties=properties, units=units)
self.shells = [] # list of dictionaries
self.is_built = False
def build(self):
pass
def add_shell(self, key, qualifier=None):
pass
[docs]class PdbProtein(ShellShape):
"""
Provides a class to import PDB protein files into hoobas in a simple fashion
Methods :
parse_pdb_protein(filename = None)
will load a protein from a pdb file. Calculates mass, moment of inertia. If filename is left to none, assumes legacy filename behavior.
add_shell(key, shell_name = None)
parses the keys in the pdb list and adds a shell representation of the parsed keys. Giving it a name allows to bind stuff onto the shell
set_ext_shell_grafts(ext_obj, num = None, linker_bond_type = None, shell_name = None)
binds num ext_objs onto the shell with name shell_name. A bond of linker_bond_type is added between the shell and external object
pdb_build_table()
needs to be called once all shells / grafts are set
"""
def __init__(self, surf_plane=None, lattice=None, properties=None, filename=None, units=None):
# default PDB units are angstroms, this can be overridden by setting units to that of the PDB file
if units is not None and isinstance(units, SimUnits):
_units = units
else:
_units = SimUnits()
_units.set_length('Ang')
_units.set_mass('amu')
shape.__init__(self, surf_plane, lattice, properties, units=_units)
if filename is not None:
self.parse_pdb_protein(filename)
# TODO : add checking against occupancy
[docs] def parse_pdb_protein(self, filename):
"""
Parses a pdbml protein and stores the file in self._pdb
:param filename:
:return:
"""
with open(filename) as f:
self._pdb = f.readlines()
del f
_m = 0.0
# inertia Vector : ixx, iyy, izz, ixy, ixz, iyz
_i_v = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0], dtype=np.longfloat) # catch some stuff close to zero
_com = np.array([0.0, 0.0, 0.0])
_geo_c = np.array([0.0, 0.0, 0.0])
_cnt = 0
for _line in self._pdb:
_s = _line.strip().split()
if _s[0] == 'ATOM':
_m_l = self.mkeys[_s[-1]]
_m += _m_l
posx = 6
posy = 7
posz = 8
if len(_s) == 10:
posx = 5
posy = 6
posz = 7
_p = np.array([float(_s[posx]), float(_s[posy]), float(_s[posz])])
_com += _p * _m_l
_geo_c += _p
_cnt += 1
_com /= _m
_geo_c /= _cnt
for _line in self._pdb:
_s = _line.strip().split()
if _s[0] == 'ATOM':
_m_l = self.mkeys[_s[-1]]
posx = 6
posy = 7
posz = 8
if len(_s) == 10:
posx = 5
posy = 6
posz = 7
_p = np.array([float(_s[posx]), float(_s[posy]), float(_s[posz])])-_com
_i_v += np.array([_p[1]**2 + _p[2]**2, _p[0]**2 + _p[2]**2, _p[0]**2 + _p[1]**2, -_p[0]*_p[1], -_p[0]*_p[2], -_p[1]*_p[2]])*_m_l
_i_v = _i_v.astype(np.float)
########
# numerical precision issues before diagonalization
#######
_maxel = np.max(_i_v.flatten())
for el_idx in range(_i_v.__len__()):
if abs(_i_v[el_idx]) < np.finfo(np.float).eps * 2.0 * _maxel:
_i_v[el_idx] = 0.0
self.flags['normalized'] = False
self.flags['mass'] = _m
self.flags['simple_I_tensor'] = False
self.flags['I_tensor'] = _i_v
self.flags['center_of_mass'] = _com
self.flags['pdb_object'] = True
pdbITensor = np.array([[_i_v[0], _i_v[3], _i_v[4]], [_i_v[3], _i_v[1], _i_v[5]], [_i_v[4], _i_v[5], _i_v[2]]])
w, v = np.linalg.eigh(pdbITensor)
idx = w.argsort(kind='mergesort')[::-1]
w = w[idx]
v = v[:, idx]
# invert one axis
if np.linalg.det(v) < 0:
v[:, 0] = -v[:, 0]
self.quaternion = Quat(v)
self.Itensor = np.array([w[0], w[1], w[2]], dtype=np.float32)
self.__check_build_order()
[docs] def build(self):
"""
creates the list for build, keeping only keyed particles.
:return:
"""
com = self.flags['center_of_mass']
self.table = np.zeros((0, 3))
_t = []
_cnt = 0
if 'shell' in self.keys:
for shl in self.keys['shell']:
for idx in shl[1]:
self.table = np.append(self.table, [np.array([idx[0], idx[1], idx[2]]) - com], axis=0)
_cnt += 1
_t.append(_cnt)
if _t.__len__() >= 1:
self.flags['multiple_surface_types'] = _t
self.is_built = True
self.quat_to_hoomd()
[docs] def add_shell(self, key, qualifier=None):
"""
Adds a shell to the protein object. Able to key multiple commands onto the same shell
self.flags['shell'] is defined as a list of [ shell_name, positions, [dictionaries] ], this backtracks all previous shells
to make sure that no particle is added twice. It is possible to create an empty shell
:param key:
:param qualifier: qualifier to describe the current set of shell beads, this can be any class that is then compared
with a supplied qualifier during the add_grafts if required. The add_graft can also be supplied with a lambda
function that has a key arguement 'qualifier' that will receive this paramter.
:return:
"""
pdb_form = {'HEAD': 0, 'RES': 5, 'RESNAME': 3, 'TYPE': -1, 'CHAIN': 4, 'OCC': 10, 'ATOM': 2}
_l_list = []
for _line in self._pdb:
_s = _line.strip().split()
if _s[0] == 'ATOM':
_l = True
for _k in key.iterkeys():
if hasattr(key[_k], '__iter__'): # argument passed is an iterable, e.g. 'CHAIN':[12, 14, 15]
inkey = False
for elem in key[_k]:
inkey = inkey or _s[pdb_form[_k]] == str(elem)
_l = _l and inkey
else:
_l = _l and _s[pdb_form[_k]] == key[_k]
for _shell_idx in range(self.keys['shell'].__len__()):
if not _l:
break
for shell in self.shells:
for position in shell['position']:
if position[0] == float(_s[6]) and position[1] == float(_s[7]) and position[2] == float(_s[8]):
_l = False
break
if _l:
posx = 6
posy = 7
posz = 8
if len(_s) == 10:
posx=5
posy=6
posz=7
_l_list.append([float(_s[posx]), float(_s[posy]), float(_s[posz])])
# find whether we need to append this to an existing shell or to create a new one
_new_shl = True
for index, shell in enumerate(self.shells): # i in range(self.keys['shell'].__len__()):
if shell['qualifier'].__class__ == qualifier.__class__ and shell['qualifier'] == qualifier:
_new_shl = False
_shl_idx = index
break
if _new_shl:
_shl_idx = -1
self.shells.append({'qualifier': qualifier, 'positions':_l_list})
else:
self.shells[_shl_idx]['positions'] += _l_list
self.__check_build_order()
def __check_build_order(self):
"""
checks that the object wasn't built when directives are passed
"""
if self.is_built:
warnings.warn('Directives were set after the shape was built. Unbuilding the shape', UserWarning)
[docs]class ProteinResidueCG(shape):
_ALL_RESIDUES = ['ARG', 'LYS', 'ASP', 'GLU', 'GLN', 'ASN', 'HIS', 'SER', 'THR', 'TYR', 'CYS',
'MET', 'TRP', 'ALA', 'ILE', 'LEU', 'PHE', 'VAL', 'PRO', 'GLY']
_ALL_RESIDUES_MASS = {'ARG': 1.0}
_ALL_ATOM_MASS = {'N': 14.0}
def __init__(self, surf_plane=None, lattice=None, properties=None, filename=None, filters=None, CG_model=None):
self.BuiltFlag = False
# define PDB in kJ/mol / Angstrom / amu
_units = SimUnits()
_units.set_energy('kJ/mol')
_units.set_length('A')
_units.set_mass('amu')
super(ProteinResidueCG, self).__init__(units= _units, properties= properties)
if filename is None:
raise SyntaxError('Valid PDB file must be supplied to this shape')
self.BuildMethod = self.pdb_build_table
# list of [residue index, residue type, positions = np.array((atom, 3)), [atom indexes], chain#]
self.residue_table = []
self.filtered_residue = []
self.removed_residues = []
# list of connectivity of the protein, [residue #1, residue #2, atomic bond index 0, atomic bond index 1]
self.connectivity_table = [] # full connectity of all residues
self.filtered_connectivity = [] # bonds that are fully in the filtered protein
self.removed_connectivity = [] # bonds that are fully removed
# chain headers listed in PDB
self.chain_head = []
self.chain_size = {}
self.chain_residues = {}
self.chain_ATOM_missing = {}
self.chain_start_index = {}
self.filtered_chain_residues = {}
self.chain_connectivity = {}
self.chain_filtered_connectivity = {}
self.chain_half_connectivity = {}
self.soft_sequences = []
# parse the protein
self.parse_pdb(filename)
self.parse_remark465()
self.parse_residues()
self.parse_connect()
self.rigid_mass = 0.0
self.atom_types = []
self.filters = []
# check type of filters, TODO : check this as filters may require chain lists
if filters is not None:
if hasattr(filters, '__iter__'):
for _filter in filters:
if not isinstance(_filter, ResidueFilter):
_lf = ResidueFilter(_filter)
else:
_lf = _filter
_lf.residues = self.residue_table
_lf.connect = self.connectivity_table
self.filters.append(_lf)
elif not isinstance(filters, ResidueFilter):
self.filters = ResidueFilter(filters)
else:
self.filters = filters
self.__apply_filtering()
self.__apply_filtering()
def parse_pdb(self, filename):
with open(filename) as f:
self._pdb = f.readlines()
@property
def residue_filter(self):
return self.filters
@residue_filter.setter
def residue_filter(self, filters):
if hasattr(filters, '__iter__'):
for _filter in filters:
if not isinstance(_filter, ResidueFilter):
_lf = ResidueFilter(_filter)
else:
_lf = _filter
self.filters.append(_lf)
elif not isinstance(filters, ResidueFilter):
self.filters = ResidueFilter(filters)
else:
self.filters = filters
self.__apply_filtering()
@staticmethod
def parse_ATOM_line(line):
if line[0:6].strip() != 'ATOM':
return ['']
return [line[0:6].strip(), int(line[6:11].strip()), line[12:16].strip(), line[16], line[17:20].strip(), line[21],
int(line[22:26].strip()), line[26], float(line[30:38].strip()), float(line[38:46].strip()), float(line[46:54].strip()),
float(line[54:60].strip()), line[60:66].strip(), line[76:78].strip(), line[78:80].strip()]
def parse_residues(self):
prev_res_number = None
# load one residue at a time
lidx = 0
while(lidx < self._pdb.__len__()):#for lidx in range(self._pdb.__len__()):
line = self._pdb[lidx]
_s = ProteinResidueCG.parse_ATOM_line(line)
if _s[0] != 'ATOM' or _s[6] == prev_res_number:
lidx += 1
continue
c_atom_list = []
c_residue = _s[4]
alt_locations = False
if _s[4] not in ProteinResidueCG._ALL_RESIDUES:
print([str(item) for item in _s])
raise AssertionError('Got unknown residue: ' + str(_s[4]) + ' in the chain')
if _s[3] != ' ':
alt_locations = True
c_residue = _s[4][1:]
cres_number = _s[6]
# append [serial, atom name, chainID, res#, x, y, z, occupancy, element+charge
c_atom_list.append([_s[1], _s[2], _s[5], _s[6], _s[8], _s[9], _s[10], _s[11], _s[12]])
# update the chain list
if _s[5] not in self.chain_head:
self.chain_head.append(_s[5])
self.chain_residues.update({_s[5]: []})
self.filtered_chain_residues.update({_s[5]: []})
try:
nval = int(_s[6])
except ValueError:
try:
nval = int(_s[6][0:-1])
except ValueError:
nval = 1
warnings.warn('Could not properly convert residue number : ' + _s[6] + ' in chain '+ str(_s[5]), RuntimeWarning)
if _s[5] not in self.chain_start_index:
self.chain_start_index.update({_s[5]: nval})
else:
self.chain_start_index[_s[5]] = min(nval, self.chain_start_index[_s[5]])
# load the rest of the residue
while lidx + 1 < self._pdb.__len__() \
and ProteinResidueCG.parse_ATOM_line(self._pdb[lidx+1])[0] != '' \
and ProteinResidueCG.parse_ATOM_line(self._pdb[lidx + 1])[6] == cres_number:
lidx += 1
_sn = ProteinResidueCG.parse_ATOM_line(self._pdb[lidx])
c_atom_list.append([_sn[1], _sn[2], _sn[5], _sn[6], _sn[8], _sn[9], _sn[10], _sn[11], _sn[12]])
# if we have multiple locations, we keep the highest occupancy factor only
if(alt_locations):
_occlist = np.array([float(el[6]) for el in c_atom_list])
_midx = _occlist.flat[abs(_occlist).argmax()]
idx_to_remove = []
for atom_idx in range(c_atom_list.__len__()):
if c_atom_list[atom_idx][5] != _midx:
idx_to_remove.append(atom_idx)
c_atom_list = [item for item in c_atom_list if c_atom_list.index(item) not in idx_to_remove]
self.residue_table.append({'pdb_res_num': cres_number,
'res_type': c_residue,
'atom_pos': np.array([[item[4], item[5], item[6]] for item in c_atom_list]),
'serial_num': [item[0] for item in c_atom_list],
'Chain': c_atom_list[0][2],
'atom_type': [item[1] for item in c_atom_list],
'atom_charge': [item[8] for item in c_atom_list]
})
self.chain_residues[c_atom_list[0][2]].append(self.residue_table[-1])
lidx += 1
# validate that residue numbers are increasing
for chain in self.chain_residues:
for residue_index in range(self.chain_residues[chain].__len__() - 1):
if self.chain_residues[chain][residue_index]['pdb_res_num'] > self.chain_residues[chain][residue_index]['pdb_res_num']:
raise warnings.warn('Potentially non-ordered tags found in residue numbers of chain : ' + str(chain)
+ ' expect bad behavior if residues are missing from ATOM', RuntimeWarning)
def parse_remark465(self):
for line in self._pdb:
_s = line.strip().split()
if _s[0] != 'REMARK':
continue
if _s[1] != '465' or _s.__len__() <= 2:
continue
if _s[2] in ProteinResidueCG._ALL_RESIDUES:
if _s[3] not in self.chain_ATOM_missing:
self.chain_ATOM_missing.update({_s[3]: []})
self.chain_start_index.update({_s[3]: 1})
self.chain_ATOM_missing[_s[3]].append(_s[4])
try:
val = int(_s[4])
self.chain_start_index[_s[3]] = min(val, self.chain_start_index[_s[3]])
except ValueError: # bad indexes in there
val = 0
self.chain_start_index[_s[3]] = min(val, self.chain_start_index[_s[3]])
if val <= 0:
warnings.warn('Found residue numbers smaller than zero in chain ' + str(_s[3]) +
', expect bad behavior if topology is not simple', RuntimeWarning)
def parse_connect(self):
prev_num = 1
prev_c = '0'
idx_a = 0
idx_b = 1
for line in self._pdb:
_s = line.strip().split()
if _s[0] != 'SEQRES':
continue
if int(_s[1]) == prev_num + 1 and prev_c == _s[2]:
self.connectivity_table.append({'Chain': _s[2], 'A': self.connectivity_table[-1]['B'], 'B': _s[4],
'A_i': idx_a, 'B_i': idx_b})
idx_a += 1
idx_b += 1
else:
idx_a = 0
idx_b = 1
prev_num = int(_s[1])
prev_c = _s[2]
if _s[2] not in self.chain_size:
self.chain_size.update({_s[2]: int(_s[3])})
for item in range(4, _s.__len__()-1):
self.connectivity_table.append({'Chain': _s[2], 'A': _s[item], 'B': _s[item+1], 'A_i': idx_a, 'B_i': idx_b})
idx_a += 1
idx_b += 1
# apply the filters in self.filters
def __apply_filtering(self):
# reset filtered lists
for dict_item in self.filtered_chain_residues:
self.filtered_chain_residues[dict_item] = []
# first pop the sequences that are in SEQRES but for which we dont have positions
# TODO : guard against bad author numbering PDB sequencing by using string match instead
for chain in range(self.chain_head.__len__()):
offset = self.chain_start_index[self.chain_head[chain]]
for residue_index in range(offset, self.chain_size[self.chain_head[chain]] + offset):
position_resolved = False
for residue in self.chain_residues[self.chain_head[chain]]:
if residue_index == residue['pdb_res_num']:
position_resolved = True
if position_resolved:
pass
else:
# find the residue type
R_type = 'NULL'
for item in self.connectivity_table:
if item['Chain'] != self.chain_head[chain]:
continue
if item['A_i'] + offset == residue_index:
R_type = item['A']
break
elif item['B_i'] + offset == residue_index:
R_type = item['B']
break
if(R_type == 'NULL'):
raise AssertionError('Could not find a residue in the chain for residue index : '
+ str(residue_index) + ' located in chain ' + str(self.chain_head[chain]))
self.filtered_chain_residues[self.chain_head[chain]].append({'pdb_res_num': residue_index,
'res_type': R_type,
'Chain': self.chain_head[chain]})
# at this point the list of filtered chain residues + chain residue equal to the whole protein
for r_filter in self.filters:
for chain_item in self.chain_residues:
for residue in self.chain_residues[chain_item]:
if r_filter(residue) and r_filter.mode == 'pop':
# pop the item and append it to the filtered list
self.filtered_chain_residues[chain_item].append(self.chain_residues[chain_item].pop(self.chain_residues[chain_item.index(residue)]))
if r_filter(residue) and r_filter.mode == 'remove':
self.chain_residues[chain_item].pop(self.chain_residues[chain_item].index(residue))
for residue in self.filtered_chain_residues[chain_item]:
if r_filter(residue) and r_filter.mode == 'remove':
self.filtered_chain_residues[chain_item].pop(self.filtered_chain_residues[chain_item].index(residue))
# rebuild the filtered connectivity
self.chain_connectivity = {}
self.chain_filtered_connectivity = {}
self.chain_half_connectivity = {}
for bond in self.connectivity_table:
# 4 types of bonds :
# type 1 : rigid-rigid bond
# type 2 : rigid-filtered bond
# type 3 : filter-filter bond
# type 4 : removed particle
type = 0
rigid_t1 = False
rigid_t2 = False
soft_t1 = False
soft_t2 = False
chain = bond['Chain']
if chain not in self.chain_connectivity:
self.chain_connectivity.update({chain: []})
self.chain_filtered_connectivity.update({chain: []})
self.chain_half_connectivity.update({chain: []})
for ritem in self.chain_residues[chain]:
if bond['A_i'] == ritem['pdb_res_num']:
rigid_t1 = True
if bond['B_i'] == ritem['pdb_res_num']:
rigid_t2 = True
if rigid_t1 and rigid_t2:
type = 1
break
if type == 0:
for ritem in self.chain_filtered_connectivity[chain]:
if bond['A_i'] == ritem['pdb_res_num']:
soft_t1 = True
if bond['B_i'] == ritem['pdb_res_num']:
soft_t2 = True
if (rigid_t1 and soft_t2):
type = 2
break
if(rigid_t2 and soft_t1):
type = 2
break
if(soft_t1 and soft_t2):
type = 3
break
if type == 0:
type = 4
continue
if type == 1:
self.chain_connectivity[chain].append({'A': bond['A'], 'B': bond['B'], 'A_i': bond['A_i'], 'B_i': bond['B_i']})
elif type == 2 and rigid_t1:
self.chain_half_connectivity[chain].append({'A': bond['A'], 'B': bond['B'], 'A_i': bond['A_i'], 'B_i': bond['B_i'], 'rigid': 'A'})
elif type == 2 and rigid_t2:
self.chain_half_connectivity[chain].append({'A': bond['A'], 'B': bond['B'], 'A_i': bond['A_i'], 'B_i': bond['B_i'], 'rigid': 'B'})
else:
self.chain_filtered_connectivity[chain].append({'A': bond['A'], 'B': bond['B'], 'A_i': bond['A_i'], 'B_i': bond['B_i']})
# yields the residue sequences that have been filtered
def get_filtered_sequences(self):
for seq in self.soft_sequences:
yield seq
def pdb_build_table(self):
# keep only the chain_residue stuff and build the table one residue at a time. These are non-bonded, but indexes
# must be recalculated first to have correct absolute hoomd indexes
self.table = np.array((0, 3), dtype=np.float32)
# center of mass should be set to 0
com = np.array([0, 0, 0], dtype=np.float32)
total_mass = 0.0
mI = np.array((3, 3), dtype=np.float32)
for chain in self.chain_head:
for residue in self.chain_residues[chain]:
for atom in range(residue['atom_pos'].__len__()):
com += residue['atom_pos'][atom] * ProteinResidueCG._ALL_ATOM_MASS[residue['atom_charge'][atom]]
total_mass += ProteinResidueCG._ALL_ATOM_MASS[residue['atom_charge'][atom]]
com /= total_mass
self.rigid_mass = total_mass
for chain in self.chain_head:
for residue in self.chain_residues[chain]:
for atom in range(residue['atom_pos'].__len__()):
residue['atom_pos'][atom] -= com
mI[0, 0] += (residue['atom_pos'][atom][1] * residue['atom_pos'][atom][1] + residue['atom_pos'][atom][2] * residue['atom_pos'][atom][2]) * ProteinResidueCG._ALL_ATOM_MASS[residue['atom_type'][atom]]
mI[1, 1] += (residue['atom_pos'][atom][0] * residue['atom_pos'][atom][0] + residue['atom_pos'][atom][2] * residue['atom_pos'][atom][2]) * ProteinResidueCG._ALL_ATOM_MASS[residue['atom_type'][atom]]
mI[2, 2] += (residue['atom_pos'][atom][1] * residue['atom_pos'][atom][1] + residue['atom_pos'][atom][0] * residue['atom_pos'][atom][0]) * ProteinResidueCG._ALL_ATOM_MASS[residue['atom_type'][atom]]
for elx in range(3):
for ely in range(elx + 1, 3):
mI[elx, ely] -=residue['atom_pos'][atom][elx] * residue['atom_pos'][atom][ely] * ProteinResidueCG._ALL_ATOM_MASS[residue['atom_type'][atom]]
mI[ely, elx] -= residue['atom_pos'][atom][elx] * residue['atom_pos'][atom][ely] * ProteinResidueCG._ALL_ATOM_MASS[residue['atom_type'][atom]]
########
# numerical precision issues before diagonalization
#######
_maxel = np.max(mI.flatten())
for el_idx in range(3):
for el_idx2 in range(3):
if abs(mI[el_idx, el_idx2]) < np.finfo(np.float).eps * 2.0 * _maxel:
mI[el_idx, el_idx2] = 0.0
self.flags['normalized'] = False
self.flags['mass'] = total_mass
self.flags['simple_I_tensor'] = False
self.flags['I_tensor'] = np.array([mI[0,0], mI[1,1], mI[2,2], mI[0,1], mI[0,2], mI[1,2]], dtype=np.float32)
self.flags['center_of_mass'] = com
self.flags['pdb_object'] = True
w, v = np.linalg.eigh(mI)
idx = w.argsort(kind='mergesort')[::-1]
w = w[idx]
v = v[:, idx]
# invert one axis
if np.linalg.det(v) < 0:
v[:, 0] = -v[:, 0]
self.quaternion = Quat(v)
self.Itensor = np.array([w[0], w[1], w[2]], dtype=np.float32)
# rebuild the index table
idx_change_lookup = {}
for chain in self.chain_head:
if chain not in idx_change_lookup:
idx_change_lookup.update({chain: []})
for residue in self.chain_residues[chain]:
for atom in residue['atom_pos']:
self.table = np.append(self.table, atom, axis=0)
for m in residue['atom_type']:
self.masses = np.append(self.masses, float(ProteinResidueCG._ALL_ATOM_MASS[m]))
self.atom_types.append(m)
idx_change_lookup[chain].append(self.masses.__len__() - 1)
# build the soft sequences, TOTO : add position parsing into atom_pos
for chain in self.chain_head:
for middle_bonds in self.chain_filtered_connectivity:
# check if one of the bond is already in the soft sequences
exists = False
for seq in self.soft_sequences:
if seq['chain'] == chain and (middle_bonds['A_i'] in seq['pdb_num'] or middle_bonds['B_i'] in seq['pdb_num']):
# find which end of the sequence to append to
if middle_bonds['A_i'] in seq['pdb_num']:
if middle_bonds['A_i'] == seq['pdb_num'][-1]:
# push at the end
seq['pdb_num'].append(middle_bonds['B_i'])
seq['res_type'].append(middle_bonds['B'])
elif middle_bonds['A_i'] == seq['pdb_num'][0]:
# push at front
seq['pdb_num'].insert(0, middle_bonds['B_i'])
seq['res_type'].insert(0, middle_bonds['B'])
else:
raise AssertionError('Multiple bonds connecting a residue chain ? ')
elif middle_bonds['B_i'] in seq['pdb_num']:
if middle_bonds['B_i'] == seq['pdb_num'][-1]:
# push at the end
seq['pdb_num'].append(middle_bonds['A_i'])
seq['res_type'].append(middle_bonds['A'])
elif middle_bonds['B_i'] == seq['pdb_num'][0]:
# push at front
seq['pdb_num'].insert(0, middle_bonds['A_i'])
seq['res_type'].insert(0, middle_bonds['A'])
else:
raise AssertionError('Multiple bonds connecting a residue chain ? ')
break
else:
self.soft_sequences.append({'chain': chain,
'pdb_num': [middle_bonds['A_i'], middle_bonds['B_i']],
'res_type': [middle_bonds['A'], middle_bonds['B']],
})
for end_bonds in self.chain_half_connectivity:
if end_bonds['rigid'] == 'A':
rigid_num = idx_change_lookup[end_bonds['A_i']]
init_soft_n = end_bonds['B_i']
soft_t = end_bonds['B']
else:
rigid_num = idx_change_lookup[end_bonds['B_i']]
init_soft_n = end_bonds['A_i']
soft_t = end_bonds['A']
# check if the end bond is connecting a structure like Rigid - soft - Rigid
exists = False
for seq in self.soft_sequences:
if seq['chain'] == chain and init_soft_n in seq['pdb_num']:
exists = True
if init_soft_n == seq['pdb_num'][0]:
seq.update[{'left_connect': rigid_num}]
elif init_soft_n == seq['pdb_num'][-1]:
seq.update[{'right_connect': rigid_num}]
else:
raise AssertionError('soft sequences connected at the middle ? ')
break
if not exists:
self.soft_sequences.append({'chain': chain,
'pdb_num': [init_soft_n],
'res_type': [soft_t],
'left_connect': rigid_num})
def get_chain_index(self, chain):
idx_l = []
for item in self.residue_table:
if item['Chain'] == chain:
idx_l.append(item['pdb_res_num'])
return idx_l
# TODO : change this to take in external models instead of only one
def apply_CG(self):
pass
# specialized dict class
class ResidueFilter(object):
def __init__(self, *args, **kwargs):
self.filter_dict = dict(**kwargs)
self.residue_list = []
self.connectivity = []
@property
def connect(self):
return self.connectivity
@connect.setter
def connect(self, val):
self.connectivity = val
@property
def residues(self):
return self.residue_list
@residues.setter
def residues(self, val):
self.residue_list = val
def __call__(self, val): # checks if the filter applies to the supplied residue / index
if isinstance(val, dict): # residues are supplied as dicts
pass
else: # otherwise assume it's by index as the filter may depend on previous or next residues
pass
[docs]class Sphere(shape):
def __init__(self, Num, properties=None, **kwargs):
shape.__init__(self, properties, **kwargs)
# a sphere of r = 1
self.table = np.zeros((0, 3), dtype=np.float32)
#gold_ang = pi*(3-5**0.5)
#th = gold_ang*np.arange(Num)
gold = (5.0**0.5 + 1)*0.5
th = np.arange(Num) / (gold)
th = 2 * pi * (th - np.floor(th))
if Num > 0:
z = np.linspace(1 - 1.0/Num, 1.0/Num - 1, Num)
else:
z = np.array([])
for i in range(0, Num):
self.table = np.append(self.table, np.array([[(1-z[i]**2)**0.5*cos(th[i]), (1-z[i]**2)**0.5*sin(th[i]),z[i]]]),axis=0)
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = 1
self.flags['surface'] = 4*pi
self.flags['volume'] = (4.0/3.0)*pi
self.flags['simple_I_tensor'] = True
if Num > 0:
self.set_geometric_quaternion()
else:
self.quaternion = Quat(np.array([1.0, 0.0, 0.0, 0.0]))
self.Itensor = np.array([0.0, 0.0, 0.0], dtype=np.float32)
[docs]class RhombicDodecahedron(shape):
def __init__(self, Num, properties=None):
shape.__init__(self, properties)
# Generate one face
ang = acos(1./3.)
Numline = int(round(sqrt(Num/12)))
xline = np.linspace(0.0, 1.0, Numline, dtype=np.float32) - (1.0+cos(ang))/2.0
yline = np.linspace(0.0, 1.0 * sin(ang), Numline, dtype=np.float32) - sin(ang) / 2.0
linedst = xline[1]-xline[0]
TableFace = np.zeros((0, 3))
for i in range(xline.__len__()):
for j in range(yline.__len__()):
TableFace = np.append(TableFace, [[xline[i] + cos(ang)/(Numline-1) * j, yline[j], 0]], axis=0)
arot = -ang
for i in range(TableFace.__len__()):
TableFace[i,:] = np.array([[TableFace[i,0]*cos(arot)-TableFace[i,1]*sin(arot), TableFace[i,0]*sin(arot)+TableFace[i,1]*cos(arot),0]])
TableFace = TableFace - [[0, 0, 6.0**0.5/3.0]]
TableRotP = np.copy(TableFace)
TableRotM = np.copy(TableFace)
for i in range(TableFace.__len__()):
# Rotation 60 deg / x+y axis
ux = 1
uy = 0
L = 0*-sin(pi/2 - ang)
an = pi/3
TableRotP[i,:] = np.array([[(cos(an)+ux**2*(1-cos(an)))*TableRotP[i,0]+ux*uy*(1-cos(an))*TableRotP[i,1]+uy*sin(an)*TableRotP[i,2]+L, uy*ux*(1-cos(an))*TableRotP[i,0] + (cos(an) + uy**2*(1-cos(an)))*TableRotP[i,1] - ux*sin(an)*TableRotP[i,2], -uy*sin(an)*TableRotP[i,0] +ux*sin(an)*TableRotP[i,1]+cos(an)*TableRotP[i,2]]])
ux = cos(ang/2)
uy = sin(ang/2)
vl = ux*TableRotP[i,0] + uy*TableRotP[i,1]
TableRotM[i,:] = np.array([[2*vl*ux - TableRotP[i,0], 2*vl*uy - TableRotP[i,1], TableRotP[i,2]]])
TableTot = np.zeros((0,3))
for i in range(TableFace.__len__()):
TableTot = np.append(TableTot, [[TableFace[i, 0], TableFace[i, 1], TableFace[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[TableFace[i, 0], TableFace[i, 1], -TableFace[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[-TableRotP[i, 0], TableRotP[i,1], TableRotP[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[-TableRotP[i, 0], TableRotP[i,1], -TableRotP[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[TableRotP[i, 0], -TableRotP[i,1], -TableRotP[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[TableRotP[i, 0], -TableRotP[i,1], TableRotP[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[-TableRotM[i, 0], TableRotM[i,1], TableRotM[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[-TableRotM[i, 0], TableRotM[i,1], -TableRotM[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[TableRotM[i, 0], -TableRotM[i,1], -TableRotM[i,2]]],axis = 0)
TableTot = np.append(TableTot, [[TableRotM[i, 0], -TableRotM[i,1], TableRotM[i,2]]],axis = 0)
ux = sin(ang/2)
uy = cos(ang/2)
uz = 0
TableRotS1 = np.copy(TableFace)
for i in range(TableFace.__len__()):
an = pi/2
TableRotS1[i,:] = np.array([[(cos(an)+ux**2*(1-cos(an)))*TableRotS1[i,0]+(ux*uy*(1-cos(an))-uz*sin(an))*TableRotS1[i,1]+(ux*uz*(1-cos(an))+uy*sin(an))*TableRotS1[i,2], (uy*ux*(1-cos(an))+uz*sin(an))*TableRotS1[i,0] + (cos(an) + uy**2*(1-cos(an)))*TableRotS1[i,1] +(uy*uz*(1-cos(an))- ux*sin(an))*TableRotS1[i,2], (ux*uz*(1-cos(an))-uy*sin(an))*TableRotS1[i,0] +(uz*uy*(1-cos(an))+ux*sin(an))*TableRotS1[i,1]+(cos(an)+uz**2*(1-cos(an)))*TableRotS1[i,2]]])
for i in range(TableFace.__len__()):
TableTot = np.append(TableTot, [[TableRotS1[i,0], TableRotS1[i,1], TableRotS1[i,2]]], axis = 0)
TableTot = np.append(TableTot, [[-TableRotS1[i,0], -TableRotS1[i,1], -TableRotS1[i,2]]], axis = 0)
self.table = 2*TableTot
self.remove_duplicates()
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = 3**0.5
self.flags['volume'] = 16.0 / 9.0 *3**0.5 * 2.0**3
self.flags['surface'] = 8*2**0.5 * 2**2
self.flags['simple_I_tensor'] = True
self.flags['call'] = 'rh_dodec'
self.set_geometric_quaternion()
[docs]class Tetrahedron(shape):
def __init__(self, Num, properties=None):
shape.__init__(self, properties)
FacePoint = Num / 4
WhP = 1
Table =np.array([[0, 0, 0]], dtype=np.float32)
TableAdd = np.zeros((1,3), dtype=np.float32)
tol = 1e-5
while WhP < FacePoint:
TableAdd = Table + [[1, 0, 0]]
TableAdd = np.append(TableAdd,Table+ [[0.5, 3**0.5/2, 0]], axis = 0)
Table = np.append(Table, TableAdd, axis = 0)
DupRows = []
for i in range(Table.__len__()-1):
for j in range(i+1, Table.__len__()):
if (Table[i,0]-Table[j,0])**2 +(Table[i,1]-Table[j,1])**2 < tol:
DupRows.append(j)
DupRows = np.unique(DupRows)
Table = np.delete(Table, DupRows, axis=0)
WhP = Table.__len__()
FacePoint = Table.__len__()
Table /= np.max(Table)
_ang = acos(1.0 / 3.0)
_r_mat = get_rotation_matrix([0.0, sin(_ang), cos(_ang)])
TableZ = np.copy(Table)
for i in range(TableZ.__len__()):
TableZ[i,:] = np.dot(TableZ[i,:],_r_mat)
self.table = np.append(self.table, Table - [0.5, sqrt(3) / 6.0, 0.0], axis=0)
TableZ -= [0.5, sqrt(3) / 6.0, 0.0]
self.table = np.append(self.table, TableZ, axis=0)
_ang = 2*pi/3.0
_r_mat = np.array([[cos(_ang), sin(_ang), 0.0], [-sin(_ang), cos(_ang), 0.0], [0.0, 0.0, 1.0]])
for i in range(TableZ.__len__()):
TableZ[i,:] = np.dot(TableZ[i,:], _r_mat)
self.table = np.append(self.table, TableZ, axis=0)
for i in range(TableZ.__len__()):
TableZ[i,:] = np.dot(TableZ[i,:], _r_mat)
self.table = np.append(self.table, TableZ, axis=0)
DupRows = []
for i in range(self.table.__len__()):
for j in range(i+1, self.table.__len__()):
if (self.table[i, 0]-self.table[j, 0])**2 + (self.table[i, 1]-self.table[j, 1])**2 + \
(self.table[i, 2] - self.table[j, 2])**2 < tol:
DupRows.append(j)
DupRows = np.unique(DupRows)
self.table = np.delete(self.table, DupRows, axis=0)
# shift the center to origin
for i in range(self.table.__len__()):
self.table[i, 2] -= 6**0.5 / 12.0
self.table *= 2.0
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = 2.0 * 3**0.5 / 8.0**0.5
self.flags['volume'] = 2.0**3.0 / (sqrt(2.0) * 6.0)
self.flags['surface'] = 2.0**2 * sqrt(3.0)
self.flags['simple_I_tensor'] = True
self.flags['call'] = 'Tetrahedron'
self.set_geometric_quaternion()
[docs]class Cylinder(shape):
def __init__(self, Num, properties=None, RadiusOverLength=None):
shape.__init__(self, properties)
self.BuildMethod = self.BuildCylinder
self.RadiusOverLength = RadiusOverLength
self.Num = Num
# Check if Radius have been parsed
if RadiusOverLength is None:
self.BuiltFlag = False
else:
self.BuiltFlag = True
self.BuildCylinder()
def BuildCylinder(self):
if self.RadiusOverLength is None:
return
if not isinstance(self.RadiusOverLength, float):
if hasattr(self.RadiusOverLength, '__call__'):
self.RadiusOverLength = self.RadiusOverLength()
else:
self.RadiusOverLength = float(self.RadiusOverLength)
if not isinstance(self.RadiusOverLength, float):
raise RuntimeError('Cannot coerce RadiusOverLength into a floating point type')
nbr_beads_circle = int(self.Num * self.RadiusOverLength * 2 * np.pi)
# azimuthal angle given number of beads
delta_angle = 2. * np.pi / nbr_beads_circle
# array of angle to loop around
angles = np.array(range(0, nbr_beads_circle)) * delta_angle
# shifted angles to loop around when ring of beads is shifted
shift_angles = angles + 0.5 * delta_angle
# distance between neighboring beads in ring around circunferece
dS = self.RadiusOverLength * delta_angle
# number of nbr_rings to fit given the neighboring distance between beads in a ring
nbr_rings = int(1.0 / dS)
# Since we have shifted and non-shifted nbr_rings, we need to have an even number of nbr_rings to fit along the box
if nbr_rings % 2 == 1:
nbr_rings = nbr_rings + 1
# ****Note: when fitting the number of rings in the simulation along the z-direction we are
# creating and anisotropic surface density along the symmetry axis of Cylinder
# shift along Z-directions between layers
dZ = float(1.0 / nbr_rings)
self.surf_ratio = dS / dZ
# loop to Create Cylinder
for i in range(nbr_rings):
# loop to add shifted/non-shifted rings
if i % 2 == 0:
# non-shifted ring
_ang = np.array(angles)
else:
# shifted ring
_ang = np.array(shift_angles)
# translate ring along Z-directions
z = dZ * i - 0.5
for j in _ang:
# loop through beads of ring
x, y = self.RadiusOverLength * cos(j), self.RadiusOverLength * sin(j)
# Add rings coordinates to global container
self.table = np.append(self.table, [[x, y, z]], axis=0)
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = self.RadiusOverLength
self.flags['volume'] = np.pi * self.RadiusOverLength**2.0
self.flags['surface'] = 2 * np.pi * self.RadiusOverLength
self.flags['simple_I_tensor'] = True
self.flags['call'] = 'Cylinder'
self.set_geometric_quaternion()
def surface_normal(self, v):
return np.array([v[0], v[1], 0.0])
[docs]class SquareRod(shape):
def __init__(self, Num, surf_plane=None, lattice=None, properties=None, SideOverLength=None):
shape.__init__(self, surf_plane, lattice, properties)
#if SideOverLength is None:
# raise RuntimeError("Didn't set side size")
self.SideOverLength = SideOverLength
self.est_number = Num
self.BuildMethod = self.buildRod
if not isinstance(self.SideOverLength, float):
self.BuiltFlag = False
else:
self.BuiltFlag = True
self.buildRod()
def buildRod(self):
if self.SideOverLength is None:
return
if not isinstance(self.SideOverLength, float):
self.SideOverLength = self.SideOverLength()
if not isinstance(self.SideOverLength, float):
raise SyntaxError("Cannot coerce SideOverLength ratio to a floating point number")
# overall number of beads ~ 4 NS * L
Ns = int(np.sqrt(self.est_number / 4.0 * self.SideOverLength))
# vertice based on symmetric center
x0 = np.linspace(-self.SideOverLength * 0.5, self.SideOverLength * 0.5, Ns)
y0 = np.linspace(-self.SideOverLength * 0.5, self.SideOverLength * 0.5, Ns)
y1 = y0[1:-1]
# horizonals
hL, hH = np.zeros((len(x0), 3)), np.zeros((len(x0), 3))
hL[:, 0], hH[:, 0] = x0, x0
hL[:, 1], hH[:, 1] = np.ones(len(x0)) * y0[0], np.ones(len(x0)) * y0[-1]
# verticals with horizontals vertices
vL, vR = np.zeros((len(y1), 3)), np.zeros((len(y1), 3))
vL[:, 1], vR[:, 1] = y1, y1
vL[:, 0], vR[:, 0] = np.ones(len(y1)) * x0[0], np.ones(len(y1)) * x0[-1]
# square perimeter
sq = np.append(hL, hH, axis=0)
sq = np.append(sq, vL, axis=0)
sq = np.append(sq, vR, axis=0)
# square height and surface
nh = int(1.0 / (x0[1] - x0[0]))
z0 = np.linspace(-0.5, 0.5, nh)
sq_x = np.shape(sq)[0]
for k in range(1, len(z0)):
sq[:, 2] = np.ones(sq_x) * z0[k]
self.table = np.append(self.table, sq, axis=0)
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = 2 ** 0.5 * self.SideOverLength
self.flags['volume'] = self.SideOverLength ** 2.0
self.flags['surface'] = self.SideOverLength * 4.0
self.flags['simple_I_tensor'] = True
self.flags['call'] = 'Square'
self.set_geometric_quaternion()
def surface_normal(self, v):
return np.array([np.sign(v[0]) if abs(v[0]) >= abs(v[1]) else 0.0, np.sign(v[1]) if abs(v[1]) >= abs(v[0]) else 0.0, 0.0])
[docs]class TruncatedTetrahedron(shape):
"""
Creates a truncated tetrahedron
V4 V1
. . . .
V5....V6 . V4 .
. . . . . .
. . . V5 .... V6 .
Side View-> V2............V3 Top View-> V2 . . . . . . . .V3
"""
def __init__(self, ratio=0.5, spacing=0.05, properties=None):
shape.__init__(self, properties)
# define vertices
# bottom three vertices of a (truncated) tetrahedron
v1 = np.array([0, sqrt(3)/3, 0], dtype=np.float32)
v2 = np.array([-0.5, -sqrt(3)/6, 0], dtype=np.float32)
v3 = np.array([0.5, -sqrt(3)/6, 0], dtype=np.float32)
# top vertex of a tetrahedron
v0 = np.array([0, 0, sqrt(6)/3])
# top three vertices of a truncated tetrahedron
# ratio is defined as the height ratio of truncated and regular tetrahedra. For example, ratio = 1 means not truncated.
v4 = ratio*v0 + (1-ratio)*v1
v5 = ratio*v0 + (1-ratio)*v2
v6 = ratio*v0 + (1-ratio)*v3
# define edges
v12 = v2 - v1
v14 = v4 - v1
v13 = v3 - v1
v25 = v5 - v2
v23 = v3 - v2
v45 = v5 - v4
v46 = v6 - v4
num_per_edge = int(1.0/spacing)
num_per_top_edge = int(1.0*(1-ratio)/spacing)
num_per_truncated_edge = int(1.0*ratio/spacing)
# edge units
v12 = v12/num_per_edge
v13 = v13/num_per_edge
v23 = v23/num_per_edge
if num_per_truncated_edge != 0:
v14 = v14/num_per_truncated_edge
v25 = v25/num_per_truncated_edge
if num_per_top_edge != 0:
v45 = v45/num_per_top_edge
v46 = v46/num_per_top_edge
positions = []
# add 3 vertexes V1, V2, V3 and V0
#positions.append(v0)
positions.append(v1)
positions.append(v2)
positions.append(v3)
# create 3 side surfaces
for i in range(num_per_truncated_edge):
for j in range(num_per_edge-i):
pos1 = (0.5+i)*v14 + (0.5+j)*v12 + v1
pos2 = (0.5+i)*v25 + (0.5+j)*v23 + v2
pos3 = (0.5+i)*v14 + (0.5+j)*v13 + v1
positions.append(pos1)
positions.append(pos2)
positions.append(pos3)
# create bottom surface
for i in range(num_per_edge):
for j in range(num_per_edge-i):
pos4 = (0.5+i)*v12 + (0.5+j)*v13 + v1
positions.append(pos4)
# create top surface
for i in range(num_per_top_edge):
for j in range(num_per_top_edge-i):
pos5 = (0.5+i)*v45 + (0.5+j)*v46 + v4
positions.append(pos5)
self.table = np.array(positions, dtype=np.float32)
# shift the center to origin
for i in range(self.table.__len__()):
self.table[i, 2] -= 0.5 * sqrt(6)/3 * ratio
self.table *= 2.0
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = 2.0 * 3**0.5 / 8.0**0.5
self.flags['volume'] = 2.0**3.0 / (sqrt(2.0) * 6.0) * (1 - (1-ratio)**3)
self.flags['surface'] = 2.0**2*sqrt(3.0) * (1 - 0.5 * (1-ratio)**2)
self.flags['simple_I_tensor'] = True
self.flags['call'] = 'TruncatedTetrahedron'
self.set_geometric_quaternion()
[docs]class ElongatedRhombicDodecahedron(RhombicDodecahedron):
def __init__(self, Num, ratio=None, properties=None):
RhombicDodecahedron.__init__(self, Num, properties)
if ratio and ratio > 2.0 / 3.0**0.5:
# a cube centered in origin with side length of 2.
dl = (ratio * 3**0.5 - 2.0) * 0.5
# define vertices
v1 = np.array([2+dl, 0, 0], dtype=np.float32)
v2 = np.array([1+dl, -1, 1], dtype=np.float32)
v3 = np.array([1+dl, 1, 1], dtype=np.float32)
v4 = np.array([1+dl, 1, -1], dtype=np.float32)
v5 = np.array([1+dl, -1, -1], dtype=np.float32)
v6 = np.array([dl, 0, 2], dtype=np.float32)
v7 = np.array([dl, 2, 0], dtype=np.float32)
v8 = np.array([dl, 0, -2], dtype=np.float32)
v9 = np.array([dl, -2, 0], dtype=np.float32)
# by symmetry
v11 = np.array([-(2+dl), 0, 0], dtype=np.float32)
v12 = np.array([-(1+dl), -1, 1], dtype=np.float32)
v13 = np.array([-(1+dl), 1, 1], dtype=np.float32)
v14 = np.array([-(1+dl), 1, -1], dtype=np.float32)
v15 = np.array([-(1+dl), -1, -1], dtype=np.float32)
v16 = np.array([-dl, 0, 2], dtype=np.float32)
v17 = np.array([-dl, 2, 0], dtype=np.float32)
v18 = np.array([-dl, 0, -2], dtype=np.float32)
v19 = np.array([-dl, -2, 0], dtype=np.float32)
# define faces (12 total)
# rhombus (8)
# e.g. v1, v2, v3, v6 (v1 + i * v12 + j * v13)
Numline = int(round(sqrt(Num/12)))
positions = []
for i in range(Numline+1):
for j in range(Numline+1):
# v1, v2, v3, v6
point = v1 + 1.0 * i / Numline * (v2-v1) + 1.0 * j / Numline * (v3-v1)
positions.append(point)
# by symmetry mirror about yz plane
point2 = np.copy(point)
point2[0] = -point2[0]
positions.append(point2)
# v1, v3, v4, v7
point = v1 + 1.0 * i / Numline * (v3-v1) + 1.0 * j / Numline * (v4-v1)
positions.append(point)
# by symmetry mirror about yz plane
point2 = np.copy(point)
point2[0] = -point2[0]
positions.append(point2)
# v1, v4, v5, v8
point = v1 + 1.0 * i / Numline * (v4-v1) + 1.0 * j / Numline * (v5-v1)
positions.append(point)
# by symmetry mirror about yz plane
point2 = np.copy(point)
point2[0] = -point2[0]
positions.append(point2)
# v1, v5, v2, v9
point = v1 + 1.0 * i / Numline * (v5-v1) + 1.0 * j / Numline * (v2-v1)
positions.append(point)
# by symmetry mirror about yz plane
point2 = np.copy(point)
point2[0] = -point2[0]
positions.append(point2)
# six edge faces (4)
# e.g. v2, v6, v9, v12, v16, v19
for i in range(Numline+1):
for j in range(Numline+1-i):
# v2, v6, v9
point = v2 + 1.0 * i / Numline * (v6-v2) + 1.0 * j / Numline * (v9-v2)
positions.append(point)
# by symmetry mirror about yz plane
point2 = np.copy(point)
point2[0] = -point2[0]
positions.append(point2)
# v3, v6, v7
point = v3 + 1.0 * i / Numline * (v6-v3) + 1.0 * j / Numline * (v7-v3)
positions.append(point)
# by symmetry mirror about yz plane
point2 = np.copy(point)
point2[0] = -point2[0]
positions.append(point2)
# v4, v7, v8
point = v4 + 1.0 * i / Numline * (v7-v4) + 1.0 * j / Numline * (v8-v4)
positions.append(point)
# by symmetry mirror about yz plane
point2 = np.copy(point)
point2[0] = -point2[0]
positions.append(point2)
# v5, v8, v9
point = v5 + 1.0 * i / Numline * (v8-v5) + 1.0 * j / Numline * (v9-v5)
positions.append(point)
# by symmetry mirror about yz plane
point2 = np.copy(point)
point2[0] = -point2[0]
positions.append(point2)
# gap is ractangle faces
# update Numline2 on 12/28/2017 see notes
Numline2 = int(round(Numline * dl * 2.0))
#Numline2 = int(round(Numline * dl * 2.0 / 3**0.5)) # old definition
for i in range(Numline+1):
for j in range(Numline2+1):
# v6, v7, v16, v17
point = v6 + 1.0 * i / Numline * (v7-v6) + 1.0 * j / Numline2 * (v16-v6)
positions.append(point)
# v7, v8
point = v7 + 1.0 * i / Numline * (v8-v7) + 1.0 * j / Numline2 * (v17-v7)
positions.append(point)
# v8, v9
point = v8 + 1.0 * i / Numline * (v9-v8) + 1.0 * j / Numline2 * (v18-v8)
positions.append(point)
# v9, v6
point = v6 + 1.0 * i / Numline * (v9-v6) + 1.0 * j / Numline2 * (v16-v6)
positions.append(point)
positions.append(v1)
positions.append(v2)
positions.append(v3)
positions.append(v4)
positions.append(v5)
positions.append(v6)
positions.append(v7)
positions.append(v8)
positions.append(v9)
positions.append(v11)
positions.append(v12)
positions.append(v13)
positions.append(v14)
positions.append(v15)
positions.append(v16)
positions.append(v17)
positions.append(v18)
positions.append(v19)
self.table = np.array(positions, dtype=np.float32)
self.remove_duplicates()
self.flags['normalized'] = True
self.flags['hard_core_safe_dist'] = 3**0.5
self.flags['volume'] = 16.0 / 9.0 *3**0.5 * 2.0**3
self.flags['surface'] = 8*2**0.5 * 2**2
self.flags['simple_I_tensor'] = True
self.flags['call'] = 'el_rh_dodec'
self.set_geometric_quaternion()
else:
warnings.warn("Ratio either not set or too small, reverting to standard RD", RuntimeWarning)