#coding: utf-8
import random
#from fractions import gcd
from math import *
import warnings
import functools
import numpy as np
import copy
from hoobas.Util import get_inverse_rotation_matrix
from hoobas.Util import uniquetol
from hoobas.Units import SimulationUnits as SimUnits
[docs]class Domain(object):
"""
Class to manage the simulation domain managed by the Builder class. Handles box coordinates. This should not be
called directly but inherited
Args:
units (Units.SimulationUnits): units used for dimensions
Attributes:
table: stores a list of tuples (position, CompositeObject)
"""
def __init__(self, units=None):
self.table = []
self.BoxSize = 0
self.flags = {}
# setup box / lattice vectors
self.vx = np.zeros((3,), dtype=np.float)
self.vy = np.zeros((3,), dtype=np.float)
self.vz = np.zeros((3,), dtype=np.float)
self.lattice = np.zeros((3, 3), dtype=np.float)
# lattice multipliers
self.bounds = np.array([0.5, 0.5, 0.5])
if units is not None and issubclass(units.__class__, SimUnits):
self.units = units
else:
self.units = SimUnits()
@property
def direct_lattice_vectors(self):
"""
returns crystal lattice vectors as a tuple
:return: crystal vectors
:rtype: tuple(np.ndarray)
"""
return self.vx, self.vy, self.vz
[docs] def change_units(self, new_units):
"""
changes the unit system
:param new_units: units to change this object to
:type new_units: Units.SimulationUnits
:return: None
"""
_mul = new_units.get_length(1.0, self.units.lunit)
for el_tab in self.table:
el_tab['position'] *= _mul
for item in self.vx:
item *= _mul
for item in self.vy:
item *= _mul
for item in self.vz:
item *= _mul
self.BoxSize *= _mul * _mul * _mul
self.lattice *= _mul
self.units = new_units
[docs] @staticmethod
def parse_input_lattice(lattice):
"""
try to parse the supplied lattice, which can be a constant for cubic lattices, three numbers for orthorhombic or
nine for other symmetries
:param lattice: lattice to parse
:type lattice: np.ndarray, float, list[float] or list[list[float]]
:return: parsed lattice
:rtype: np.ndarray
"""
_e = IndexError('Domain : Lattice : parse_input_lattice : cannot parse the supplied lattice')
if isinstance(lattice, np.ndarray) and lattice.shape == (3, 3):
return lattice
if isinstance(lattice, int):
lattice = float(lattice)
if isinstance(lattice, float):
return np.array([[lattice, 0, 0], [0, lattice, 0], [0, 0, lattice]], dtype=np.float32)
elif hasattr(lattice, '__iter__'):
if isinstance(lattice[0], float):
return np.array([[lattice[0], 0, 0], [0, lattice[1], 0], [0, 0, lattice[2]]], dtype=np.float32)
elif hasattr(lattice[0], '__iter__'):
return np.array(
[[lattice[0][0], lattice[0][1], lattice[0][2]], [lattice[1][0], lattice[1][1], lattice[1][2]],
[lattice[2][0], lattice[2][1], lattice[2][2]]], dtype=np.float32)
else:
raise _e
else:
raise _e
[docs] def add_random_positions(self, composite, number=0, size=None):
"""
Adds composite objects at random positions inside the volume
:param composite: object to add inside the volume
:type composite: Composite.CompositeObject
:param number: number of objects to add
:type number: int
:param size: excluded volume of the appended object. None avoids calculation of excluded volumes
:type size: float or None
:return: None
"""
added_table = []
# counting variables
Padd = 0 # added particles
current_trial = 0 # number of trials on the current particle
Nfailed = 0 # number of trials of the whole method
while Padd < number:
# generate a random position
gpos = np.zeros(3, dtype=np.float32)
for dim in range(3):
gpos += np.random.uniform(-self.bounds[dim], self.bounds[dim]) * self.lattice[dim, :]
# check whether that position is actually valid
reject = False
for table_tuple in self.table:
dist = np.linalg.norm(table_tuple['position'] - gpos)
this_size = size if size is not None else 0.0
other_size = table_tuple['size'] if table_tuple['size'] is not None else 0.0
if dist < (other_size + this_size) / 2.0:
reject = True
break
if not reject:
Padd += 1
current_trial = 0
added_table.append({'position': gpos, 'object': composite, 'size': size, 'dimension': 0})
else:
current_trial += 1
if current_trial > 1000:
if(Nfailed > 1000):
raise RuntimeError('Domain : Lattice : Unable to insert all requested particles ')
Nfailed += 1
Padd = 0
current_trial = 0
added_table = []
self.table.extend(added_table)
[docs]class Lattice(Domain):
"""
provides a class for building lattices of :py:class:`hoobas.Composite.CompositeObject`. Class provides methods to
add particles inside the unit cell by use of an offset, a rotation method and a cut method
Args:
lattice (float, 1 x 3 iterable, 3 x 3 iterable): crystal unit cell
"""
def __init__(self, lattice, units=None):
if units is not None and issubclass(units.__class__, SimUnits):
_units = units
else:
_units = SimUnits()
Domain.__init__(self, units=_units)
self.flags['vertical_slice'] = False
self.lattice = self.parse_input_lattice(lattice)
self.unit_volume = np.dot(self.lattice[0], np.cross(self.lattice[1], self.lattice[2]))
self.reciprocal = 2*pi*np.transpose(np.linalg.inv(self.lattice))
self.bounds = [0.5, 0.5, 0.5]
self.rotation_matrix = np.eye(3)
[docs] @staticmethod
def d_hkl(lattice, hkl):
"""
Calculates the distance between planes
:param lattice: crystal lattice
:type lattice: np.ndarray
:param hkl: miller indices defining the plane
:type hkl: list[int]
:return: distance
:rtype: float
"""
a = np.linalg.norm(lattice[0, :])
b = np.linalg.norm(lattice[1, :])
c = np.linalg.norm(lattice[2, :])
unit_volume = abs(np.dot(np.cross(lattice[0, :], lattice[1, :]), lattice[2, :]))
alpha = acos(np.dot(lattice[1], lattice[2]) / np.linalg.norm(lattice[1]) / np.linalg.norm(lattice[2]))
beta = acos(np.dot(lattice[0], lattice[2]) / np.linalg.norm(lattice[0]) / np.linalg.norm(lattice[2]))
gamma = acos(np.dot(lattice[0], lattice[1]) / np.linalg.norm(lattice[1]) / np.linalg.norm(lattice[0]))
return (b * b * c * c * sin(alpha) ** 2 * hkl[0] ** 2 + a * a * c * c * sin(beta) ** 2 * hkl[1] ** 2 +
a * a * b * b * sin(gamma) ** 2 * hkl[2] ** 2 +
2.0 * a * b * c * c * (cos(alpha) * cos(beta) - cos(gamma)) * hkl[0] * hkl[1] +
2.0 * a * a * b * c * (cos(beta) * cos(gamma) - cos(alpha)) * hkl[1] * hkl[2] +
2.0 * a * b * b * c * (cos(gamma) * cos(alpha) - cos(beta)) * hkl[0] * hkl[2]
) ** -0.5 * unit_volume
[docs] @staticmethod
def normal_to(lattice, hkl):
"""
Calculates the normal vector to a crystal plane
:param lattice: crystal lattice
:type lattice: np.ndarray
:param hkl: miller indices defining the plane
:type hkl: list[int]
:return: normal vector
:rtype: np.ndarray
"""
recip = 2.0 * pi * np.transpose(np.linalg.inv(lattice))
return np.array([hkl[0] * recip[0, :] + hkl[1] * recip[1, :] + hkl[2] * recip[2, :]])
[docs] def cut_to_dimensions(self, bounds):
"""
rotates the crystal system so that the surface plane face supplied faces the Z direction. The lattice is cut
to match the integer dimensions supplied in bounds.
:param bounds: bounds of the cell
:type bounds: list[int]
:return: none
"""
for item in self.table:
item['position'] = self.rotation_matrix.dot(item['position'])
self._expand_unit_cell([int(np.ceil(i)) for i in bounds]) # this is a # of replicas, not where we cut
b1 = self.lattice[0, :]
b2 = self.lattice[1, :]
b3 = self.lattice[2, :]
b1 = self.rotation_matrix.dot(b1)
b2 = self.rotation_matrix.dot(b2)
b3 = self.rotation_matrix.dot(b3)
# decomposition matrix
decomposition_matrix = np.transpose(np.array([b1, b2, b3]))
objects_to_keep = []
for item in self.table:
item_position = item['position']
# decompose into base vectors b1, b2, b3
a1, a2, a3 = np.linalg.solve(decomposition_matrix, item_position)
# check lattice params, need some tolerance check in here
_tol = 1e-3
if (-bounds[0] + _tol < a1 <= bounds[0] + _tol) \
and (-bounds[1] + _tol < a2 <= bounds[1] + _tol) \
and (-bounds[2] + _tol < a3 <= bounds[2] + _tol):
objects_to_keep.append(item)
self.table = objects_to_keep
self.vx = b1
self.vy = b2
self.vz = b3
self.lattice = decomposition_matrix
self.bounds = bounds
[docs] def generate_crystal_slab(self, hkl, inplane_basis, bounds, make_at=None):
"""
rotates the crystal system so that the surface plane face supplied faces the Z direction. The new crystal axes
are generally trigonal in nature and the crystallinity in Z is not guaranteed unless the rotated axes point
along Z. Decomposition in XY plane is made along plane_vectors which are defined as integer sums of b1, b2, b3.
:param hkl: surface plane to expose, note that this is reduced so that hlk=[0,0,2] is equivalent to hkl=[0,0,1]
:type hkl: list[int]
:param inplane_basis: lattice vectors to use as basis in XY plane. For (111), a good choice is (101) and (110)
:type inplane_basis: list[list[int]]
:param bounds: lattice size in terms of inplane_basis + dhkl
:type bounds: list[int]
:param make_at: Specify where to cut the plane, either lattice point, arbitrary point or origin
:type make_at: int or list[float] or None
:return: None
"""
hkl[:] = (x // functools.reduce(gcd, hkl) for x in hkl)
surface_normal = Lattice.normal_to(self.lattice, hkl)
surface_normal /= np.linalg.norm(surface_normal)
dhkl = Lattice.d_hkl(self.lattice, hkl)
self.rotation_matrix = get_inverse_rotation_matrix(surface_normal.squeeze())
self.flags['vertical_slice'] = True # flag the builder for cut lattices
if make_at is None:
base_position = np.array([0.0, 0.0, 0.0])
elif hasattr(make_at, '__iter__'):
base_position = make_at
else:
base_position = self.table[make_at]['position']
if hkl[0] != 0: # get a point on the plane through miller indices
pt = self.lattice[0, :] * 1.0 / hkl[0]
elif hkl[1] != 0:
pt = self.lattice[1, :] * 1.0 / hkl[1]
else:
pt = self.lattice[2, :] * 1.0 / hkl[2]
d = np.dot(surface_normal, pt) # figure out the plane constant
shift = (np.dot(surface_normal, base_position) + d) / dhkl
b1 = self.rotation_matrix.dot(self.lattice[0, :])
b2 = self.rotation_matrix.dot(self.lattice[1, :])
b3 = self.rotation_matrix.dot(self.lattice[2, :])
b1, b2 = inplane_basis[0][0] * b1 + inplane_basis[0][1] * b2 + inplane_basis[0][2] * b3, \
inplane_basis[1][0] * b1 + inplane_basis[1][1] * b2 + inplane_basis[1][2] * b3
b3 = np.array([0.0, 0.0, dhkl])
b1[2] = b2[2] = 0.0
bv = [b1, b2, b3]
# expand lattice to make sure the bounds fit in
base_bounds = [sum([int(np.ceil(abs(np.dot(self.lattice[i, :], bv[j])/np.linalg.norm(self.lattice[i, :])/np.linalg.norm(bv[j])) * bounds[j])) for j in range(3)], 1) for i in range(3)]
base_bounds[2] += int(np.ceil(base_position[2]))
self._expand_unit_cell(base_bounds)
for item in self.table:
item['position'] = self.rotation_matrix.dot(item['position'])
item['object'].rotate(self.rotation_matrix)
# decomposition matrix
decomposition_matrix = np.transpose(np.array(bv))
objects_to_keep = []
for item in self.table:
current_position = item['position']
# decompose into base vectors b1, b2, b3
a1, a2, a3 = np.linalg.solve(decomposition_matrix, current_position)
# check lattice params, need some tolerance check in here
_tol = 1e-3
if (-bounds[0] + _tol < a1 <= bounds[0] + _tol) \
and (-bounds[1] + _tol < a2 <= bounds[1] + _tol) \
and (-bounds[2] - _tol < a3 + shift <= bounds[2] + _tol):
objects_to_keep.append(item)
self.vx = b1
self.vy = b2
self.vz = b3
self.table = objects_to_keep
self.lattice = decomposition_matrix
self.bounds = bounds
[docs] def add_particles_on_lattice(self, fractional_coordinates, composite, wrap=True, size=None):
"""
Adds particles at given lattice points. Objects will be built in every unit cell
:param fractional_coordinates: iterable describing the fractional crystal coordinates of the object to add
:type fractional_coordinates: list(float)
:param composite: object to add in every unit cell given by the coordinates
:type composite: Composite.CompositeObject
:param wrap: whether fractional coordinates are wrapped into the first unit cell
:type wrap: bool
:param size: excluded volume size if other objects are added
:type size: float or None
:return: None
"""
# check for wrapping back into cell
if wrap:
for dim in range(3):
while fractional_coordinates[dim] > 1.0:
fractional_coordinates[dim] -= 1.0
while fractional_coordinates[dim] < 0.0:
fractional_coordinates[dim] += 1.0
base = fractional_coordinates[0] * self.lattice[0, :] + \
fractional_coordinates[1] * self.lattice[1, :] + \
fractional_coordinates[2] * self.lattice[2, :]
item = {'object': composite,
'size': size,
'position': base,
'dimension': 0}
self.table.append(item)
def _expand_unit_cell(self, dims=None):
if dims is None:
dims = [0, 0, 0]
unit_base = copy.deepcopy(self.table)
for dim1 in range(-dims[0], dims[0]+1):
for dim2 in range(-dims[1], dims[1]+1):
for dim3 in range(-dims[2], dims[2]+1):
if dim1 == dim2 == dim3 == 0:
continue
for item in unit_base:
position = item['position'] + dim1 * self.lattice[0, :] + \
dim2 * self.lattice[1, :] + \
dim3 * self.lattice[2, :]
nitem = {'object': copy.deepcopy(item['object']),
'size': item['size'],
'position': position,
'dimension': item['dimension']}
self.table.append(nitem)
[docs] def add_layer(self, hkl, composite, displacement=None):
"""
Adds a 2D composite object laying on a crystal plane
:param hkl: miller indices of the plane
:type hkl: list[int]
:param composite: 2D object to add
:type composite: Layers.Multilayer
:param displacement: fractional distance between the origin and the plane
:type displacement: float
:return: None
"""
if displacement is None:
displacement = 0.0
# check if displacement is given as a vector or a distance along the normal vector
nv = Lattice.normal_to(self.lattice, hkl)
position = np.squeeze(displacement * Lattice.d_hkl(self.lattice, hkl) * nv / np.linalg.norm(nv))
composite.transform(lattice=self.lattice, hkl=hkl, shift=displacement)
self.table.append({'object': composite, 'size': 0.0, 'position': position, 'dimension': 2})
[docs] def add_line(self, hkl, composite, displacement=None):
"""
Adds a 1D composite object laying on a crystal direction
:param hkl: crystal direction
:type hkl: list[int]
:param composite: 1D object to add
:type composite: Layers.LineLayer
:param displacement: displacement to apply on the line object
:type displacement: list[float]
:return: None
"""
composite.transform(lattice=self.lattice, hkl=hkl)
if displacement is None:
displacement = np.array([0.0, 0.0, 0.0])
else:
displacement = displacement[0] * self.lattice[0, :] + \
displacement[1] * self.lattice[1, :] + \
displacement[2] * self.lattice[2, :]
self.table.append({'object': composite,
'size': 0.0,
'position': displacement,
'dimension': 1})
[docs] def change_units(self, new_units):
"""
changes the unit system
:param new_units: units to change this object to
:type new_units: Units.SimulationUnits
:return: None
"""
self.reciprocal /= new_units.get_length(1.0, self.units.lunit)
super(Lattice, self).change_units(new_units)
[docs]class EmptyBox(Domain):
"""
Creates an empty domain. The system_size used defines symmetry. A single values defines a cube, a 1x3 iterable
defines an orthorhombic box while a 3x3 iterable defines a triclinic symmetry.
Args:
system_size (float, 1x3 iterable, 3x3 iterable): defines the box dimensions
"""
def __init__(self, system_size, units=None):
Domain.__init__(self, units)
self.bounds = [0.5, 0.5, 0.5]
self.rand_flag = True
self.table = []
self.table_size = 0
self.lattice = self.parse_input_lattice(system_size)
self.vx = list(self.lattice[0, :])
self.vy = list(self.lattice[1, :])
self.vz = list(self.lattice[2, :])
self.BoxSize = abs(np.dot(self.lattice[0, :], np.cross(self.lattice[1, :], self.lattice[2, :])))
[docs]class EmptyCube(EmptyBox):
"""
Creates an empty cube
Args:
size: cube side length
"""
def __init__(self, size, units=None):
super(EmptyCube, self).__init__(size, units=units)
[docs]class EmptyHexagonalPrism(EmptyBox):
"""
Creates an empty hexagonal prism, aligned along z
Args:
a: length of the hexagon side
c: height of the prism
"""
def __init__(self, a, c, units=None):
_l = [[a, 0., 0.], [a*cos(radians(120.0)), a*sin(radians(120.)), 0.], [0., 0., c]]
super(EmptyHexagonalPrism, self).__init__(_l, units=units)