"""
Classes and functions for crystal planes.
Copyright Schrodinger, LLC. All rights reserved."""
# Contributor: Thomas F. Hughes
from collections import OrderedDict
from past.utils import old_div
from functools import reduce
from math import gcd
import numpy
from scipy.spatial import ConvexHull
from schrodinger.application.matsci import msutils
from schrodinger.application.matsci import shapes
from schrodinger.application.matsci.nano import xtal
from schrodinger.structutils import transform
_version = '$Revision 0.0 $'
[docs]def ext_gcd(a, b):
"""
Solve ax + by = gcd(a, b) using the Extended Euclidean Algorithm.
Return (1) the greatest common divisor (gcd) of integers a and b
and (2) the integer Bezout coefficients x and y.
:type a: int
:param a: the a coefficient
:type b: int
:param b: the b coefficient
:rtype: int, int, int
:return: gcd(a, b), Bezout coefficient x, and Bezout coefficient y
"""
iterate = lambda r_prev, q, r: (r_prev - q * r, r)
r_prev, r = a, b
x_prev, x = 1, 0
y_prev, y = 0, 1
while r != 0:
q = old_div(r_prev, r)
(r, r_prev) = iterate(r_prev, q, r)
(x, x_prev) = iterate(x_prev, q, x)
(y, y_prev) = iterate(y_prev, q, y)
return r_prev, x_prev, y_prev
[docs]def reduce_hkl(hkl):
"""
Reduce hkl to the smallest set of indices
param list hkl: Miller indices
:retype: list, int
:return: Reduced Miller indices, divisor
"""
divisor = abs(reduce(gcd, hkl))
hkl = [int(index / divisor) for index in hkl]
return hkl, divisor
[docs]class CrystalPlane(object):
"""
Manage a crystal plane object.
"""
SQUARE = OrderedDict(
[(1, numpy.array([-1.0, 1.0, 0.0])), (2, numpy.array([-1.0, -1.0,
0.0])),
(3, numpy.array([1.0, -1.0, 0.0])), (4, numpy.array([1.0, 1.0, 0.0]))])
DISTANCE_THRESH = -0.001
SAME_VECTOR_THRESH = 0.0001
SLAB_THRESHOLD = 0.0000001
[docs] def __init__(self,
h_index,
k_index,
l_index,
a_vec,
b_vec,
c_vec,
origin=None,
logger=None):
"""
Create an instance.
:type h_index: int
:param h_index: the h Miller index
:type k_index: int
:param k_index: the k Miller index
:type l_index: int
:param l_index: the l Miller index
:type a_vec: numpy.array
:param a_vec: the a lattice vector
:type b_vec: numpy.array
:param b_vec: the b lattice vector
:type c_vec: numpy.array
:param c_vec: the c lattice vector
:type origin: numpy.array
:param origin: the origin of the lattice vectors in Angstrom
:type logger: logging.getLogger
:param logger: output logger
"""
self.h_index = h_index
self.k_index = k_index
self.l_index = l_index
self.a_vec = a_vec
self.b_vec = b_vec
self.c_vec = c_vec
if origin is None:
self.origin = numpy.array(xtal.ParserWrapper.ORIGIN)
else:
self.origin = origin
self.logger = logger
self.ra_vec = None
self.rb_vec = None
self.rc_vec = None
self.inter_planar_separation = None
self.normal_vec = None
self.rotate_to_z = None
self.inv_rotate_to_z = None
self.checkMillerIndices()
self.getReciprocals()
self.getNormal()
self.u_vec = None
self.v_vec = None
self.w_vec = None
[docs] def checkMillerIndices(self):
"""
Check the user provided Miller indices.
"""
bad_indices_msg = (
'You have provided the Miller indices '
'(000). At least a single Miller index must be non-zero.')
if self.h_index == self.k_index == self.l_index == 0:
if self.logger:
self.logger.error(bad_indices_msg)
raise ValueError(bad_indices_msg)
[docs] def getReciprocals(self):
"""
Return the reciprocal lattice vectors.
:rtype: three numpy.array
:return: the three reciprocal lattice vectors.
"""
self.ra_vec, self.rb_vec, self.rc_vec = \
xtal.get_reciprocal_lattice_vectors(self.a_vec, self.b_vec, self.c_vec)
return self.ra_vec, self.rb_vec, self.rc_vec
[docs] def getNormal(self):
"""
Return the normal vector.
:rtype: numpy.array
:return: the normal vector for this plane
"""
rnormal_vec = self.h_index * self.ra_vec + self.k_index * self.rb_vec + \
self.l_index * self.rc_vec
self.inter_planar_separation = old_div(
1.0, transform.get_vector_magnitude(rnormal_vec))
self.normal_vec = \
self.inter_planar_separation * transform.get_normalized_vector(rnormal_vec)
return self.normal_vec
[docs] def getLinDepPlaneVectors(self):
"""
Return three typically used plane vectors that are linearly dependent.
:rtype: numpy.array, numpy.array, numpy.array
:return: typically used plane vectors that are linearly dependent
"""
a_vec_prime = self.k_index * self.a_vec - self.h_index * self.b_vec
b_vec_prime = self.l_index * self.a_vec - self.h_index * self.c_vec
c_vec_prime = self.l_index * self.b_vec - self.k_index * self.c_vec
return a_vec_prime, b_vec_prime, c_vec_prime
[docs] def getSimpleSlabVectors(self):
"""
This sets the simple, i.e. two of the Miller indices
are zero, transformation matrix into self.basis and
sets the simple slab vectors.
"""
miller_indices = [self.h_index, self.k_index, self.l_index]
num_zero = miller_indices.count(0)
if num_zero != 2:
raise ValueError('Number of null Miller indices is not two.')
if self.h_index:
self.basis = numpy.array([(0, 0, 1), (0, 1, 0), (-1, 0, 0)])
elif self.k_index:
self.basis = numpy.array([(1, 0, 0), (0, 0, 1), (0, -1, 0)])
elif self.l_index:
self.basis = numpy.array([(1, 0, 0), (0, 1, 0), (0, 0, 1)])
self.u_vec, self.v_vec, self.w_vec = \
self.transformVectors(self.a_vec, self.b_vec, self.c_vec)
[docs] def getSlabVectors(self):
"""
This sets the transformation matrix into self.basis. Basis
vectors are chosen such that a and b axes are in the plane
of the Miller plane, and c-axis is out of this plane (NOT
necessarily normal to it). Also sets the slab vectors.
"""
# if this is a simple case, i.e. two of the Miller
# indices are zero, then just handle it now
try:
self.getSimpleSlabVectors()
return
except ValueError:
pass
# the algorithm used here follows from that given in
# https://wiki.fysik.dtu.dk/ase/_downloads/general_surface.pdf
# to summarize, three plane vectors, a', b', and c', are
# formed from linear combinations of the three lattice vectors,
# a, b, and c, using the Miller indices:
# a' = k*a - h*b
# b' = l*a - h*c
# c' = l*b - k*c
# one of these vectors is picked to be one of the two plane
# vectors while a linear combination of the other two vectors
# is taken to form the last plane vector. this linear combination
# is taken so as to both minimize the planar area spanned by the two
# planar vectors and to make them as orthogonal as possible. these
# are determined subject to their cross product having a length of one
# inter-planar separation
# get some typcially used linearly dependent vectors that
# span the plane
a_vec_prime, b_vec_prime, c_vec_prime = self.getLinDepPlaneVectors()
# set the first plane vector with linear combinations of
# a_vec_prime and b_vec_prime, the linear combination to
# minimize the area requires solving pk + ql = 1
gcd_kl, p_coef, q_coef = ext_gcd(self.k_index, self.l_index)
u_vec = p_coef * a_vec_prime + q_coef * b_vec_prime
# if we can make the planar vectors more orthogonal then go ahead,
# in the original document cited above they mention the need for
# a threshold due to numerical precision however I have not yet
# seen a case where this is necessary
numerator = numpy.dot(u_vec, c_vec_prime)
denominator = \
numpy.dot(self.l_index * a_vec_prime - self.k_index * b_vec_prime, \
c_vec_prime)
if abs(denominator) > self.SLAB_THRESHOLD:
coef = -1 * int(msutils.roundup(numerator / denominator))
p_coef += coef * self.l_index
q_coef -= coef * self.k_index
# Set transformation matrix into self.basis
vec1 = p_coef * numpy.array((self.k_index, -self.h_index, 0)) + \
q_coef * numpy.array((self.l_index, 0, -self.h_index))
vec2 = old_div(numpy.array((0, self.l_index, -self.k_index)), \
abs(xtal.gcd(self.l_index, self.k_index)))
gcd_coef, a_coef, b_coef = ext_gcd(
p_coef * self.k_index + q_coef * self.l_index, self.h_index)
vec3 = (b_coef, a_coef * p_coef, a_coef * q_coef)
self.basis = numpy.array([vec1, vec2, vec3]).T
self.u_vec, self.v_vec, self.w_vec = \
self.transformVectors(self.a_vec, self.b_vec, self.c_vec)
[docs] def getSpanningVectors(self, st):
"""
Return the spanning vectors of this bounding box.
:type st: schrodinger.structure.Structure
:param st: the structure
:rtype: list of numpy.array
:return: contains vectors spanning the parallelepiped and its sides
"""
try:
vecs = xtal.get_vectors_from_chorus(st)
except ValueError:
params = xtal.get_lattice_param_properties(st)
vecs = xtal.get_lattice_vectors(*params)
a_vec, b_vec, c_vec = vecs
return [
a_vec, b_vec, c_vec, a_vec + b_vec, a_vec + c_vec, b_vec + c_vec,
a_vec + b_vec + c_vec
]
[docs] def getBestSpanningVector(self, st):
"""
Return the spanning vector with the largest projection onto the
plane normal vector.
:type st: schrodinger.structure.Structure
:param st: the structure
:rtype: numpy.array
:return: the best spanning vector
"""
spanning_vecs = self.getSpanningVectors(st)
unit_normal_vec = transform.get_normalized_vector(self.normal_vec)
pairs = [(x, abs(numpy.dot(x, unit_normal_vec))) for x in spanning_vecs]
return max(pairs, key=lambda x: x[1])[0]
[docs] def getNumPlanes(self, st):
"""
Return the number of planes that will fit inside the bounding box.
:type st: schrodinger.structure.Structure
:param st: the structure
:rtype: int
:return: the number of planes that will fit inside the bounding box
"""
unit_normal_vec = transform.get_normalized_vector(self.normal_vec)
spanning_vec = self.getBestSpanningVector(st)
spanning_normal_vec = abs(numpy.dot(spanning_vec,
unit_normal_vec)) * unit_normal_vec
nplanes = int(
round(
old_div(
transform.get_vector_magnitude(spanning_normal_vec),
self.inter_planar_separation)))
return nplanes
[docs] def getInterPlanarSeparation(self):
"""
Return the inter-planar separation in Angstrom.
:rtype: float
:return: the inter-planar separation in Angstrom
"""
return self.inter_planar_separation
[docs] def getRotationToZ(self):
"""
Return the rotation matrix needed to rotate this plane to the XY-plane
as well as its inverse.
:rtype: two numpy.array
:return: the rotation matrix that rotates this plane to the XY-plane and
its inverse.
"""
self.rotate_to_z = transform.get_alignment_matrix(
self.normal_vec, numpy.array(transform.Z_AXIS))
self.inv_rotate_to_z = numpy.array(numpy.matrix(self.rotate_to_z).I)
return self.rotate_to_z, self.inv_rotate_to_z
[docs] def getSquareVertices(self):
"""
Return the vertices of a square that lies in this plane. The square has
and edge-length of 2 Angstrom. It is rotated from the XY-plane, centered on
origin, into this plane.
:rtype: list of numpy.array
:return: the vertices of the squre that lies in this plane
"""
self.getRotationToZ()
vertices = []
for vertex in self.SQUARE.values():
vertices.append(
transform.transform_atom_coordinates(
numpy.copy(vertex), self.inv_rotate_to_z))
return vertices
[docs] def getParallelepipedLineSegments(self, st):
"""
Return the line segments that make this bounding box.
:type st: schrodinger.structure.Structure
:param st: the structure
:rtype: list of tuples of heads and tails of 12 line segments.
:return: the line segments that make this bounding box
"""
a_vec, b_vec, c_vec, ab_vec, ac_vec, bc_vec, abc_vec = \
self.getSpanningVectors(st)
# (start, end) line segment pairs
# yapf: disable
line_segments = [
(self.origin, self.origin + a_vec),
(self.origin, self.origin + b_vec),
(self.origin, self.origin + c_vec),
(self.origin + a_vec, self.origin + ab_vec),
(self.origin + a_vec, self.origin + ac_vec),
(self.origin + b_vec, self.origin + ab_vec),
(self.origin + b_vec, self.origin + bc_vec),
(self.origin + ab_vec, self.origin + abc_vec),
(self.origin + c_vec, self.origin + ac_vec),
(self.origin + c_vec, self.origin + bc_vec),
(self.origin + ac_vec, self.origin + abc_vec),
(self.origin + bc_vec, self.origin + abc_vec)
]
# yapf: enable
return line_segments
[docs] def getPlaneBoxIntersections(self, st, vertices):
"""
Return the points where the plane containing the specified
vertices intersects the parallelepiped.
:type st: schrodinger.structure.Structure
:param st: the structure
:type vertices: list of numpy.array
:param vertices: the vertices of the square that lies in this plane
:rtype: list of numpy.array
:return: the points of intersection
"""
num_unique = 1
aface = shapes.Face(1, list(self.SQUARE), vertices, num_unique)
intersections = []
for line_segment in self.getParallelepipedLineSegments(st):
start, end = line_segment
intersection = aface.intersectSegmentAndPlane(start, end, \
distance_thresh=self.DISTANCE_THRESH)
if intersection is not None:
if not intersections:
intersections.append(intersection)
else:
redundant = False
for point in intersections:
if transform.get_vector_magnitude(intersection - point) < \
self.SAME_VECTOR_THRESH:
redundant = True
break
if not redundant:
intersections.append(intersection)
return intersections
[docs] def getOrderedIntersections(self, intersections):
"""
Return the provided list of planar points in counter-clockwise
order.
:type intersections: list of numpy.array
:param intersections: some intersection points in a plane
:rtype: list of numpy.array
:return: those planar intersections in counter-clockwise order
"""
intersections += [old_div(sum(intersections), len(intersections))]
intersections_xy = []
for intersection in intersections:
intersection_xy = transform.transform_atom_coordinates(
numpy.copy(intersection) - self.origin, self.rotate_to_z)
z_coord = intersection_xy[2]
intersections_xy.append(numpy.delete(intersection_xy, 2))
indices = ConvexHull(numpy.matrix(intersections_xy)).vertices
intersections_xy_with_z = [
numpy.append(intersections_xy[index], z_coord) for index in indices
]
intersections = []
for intersection in intersections_xy_with_z:
point = transform.transform_atom_coordinates(
numpy.copy(intersection), self.inv_rotate_to_z) + self.origin
intersections.append(point)
return intersections
[docs] def getVertices(self, normal_vec, idx, draw_location=0, thickness=1):
"""
Return a list of numpy.array containing vertices of a plane
with the given index.
:type normal_vec: numpy.array
:param normal_vec: the normal vector
:type idx: int
:param idx: the plane index, the sign of the index controls
whether behind or ahead of the normal vector origin
:type draw_location: float
:param draw_location: specifies the starting location at which planes
will be drawn in terms of a fraction of the normal vector
which has a length of one inter-planar spacing
:type thickness: float
:param thickness: specifies the thickness or distance between
consecutive planes in terms of a fraction of the normal vector
which has a length of one inter-planar spacing
:rtype: list of numpy.array
:return: the plane vertices
"""
# get square vertices of the hkl plane that passes through (0, 0, 0)
square_vertices = self.getSquareVertices()
# define a vector from (0, 0, 0) to the given draw_location
shift_vec = self.origin + draw_location * normal_vec
# define an offset vector
offset_vec = idx * thickness * normal_vec
return [v + shift_vec + offset_vec for v in square_vertices]
[docs] def getVerticesOfAllPlanes(self,
st,
draw_location=0,
thickness=1,
also_draw_planes_behind=True):
"""
Return a list of lists of points where the set of planes
intersect the parallelepiped.
:type st: schrodinger.structure.Structure
:param st: the structure
:type draw_location: float
:param draw_location: specifies the starting location at which planes
will be drawn in terms of a fraction of the normal vector
which has a length of one inter-planar spacing
:type thickness: float
:param thickness: specifies the thickness or distance between
consecutive planes in terms of a fraction of the normal vector
which has a length of one inter-planar spacing
:type also_draw_planes_behind: bool
:param also_draw_planes_behind: whether to also draw planes behind
the specified draw location, this is in addition to always
drawing planes ahead of the draw location
:rtype: list of list of numpy.array
:return: where the planes intersect the parallelepiped
"""
# number of planes that can fit inside the largest bounding box
nplanes = self.getNumPlanes(st)
# get normal vector of proper sign
spanning_vec = self.getBestSpanningVector(st)
normal_vec = numpy.sign(
numpy.dot(spanning_vec,
transform.get_normalized_vector(
self.normal_vec))) * self.normal_vec
# if the zero plane (idx_plane == 0) doesn't intersect the parallelepiped
# then create one outside the box at the parallelepiped origin at the start
# of the normal vector, the size of the plane will be proportional to the
# length of the normal vector
if also_draw_planes_behind:
idxs_plane = range(-nplanes, nplanes + 1)
else:
idxs_plane = range(nplanes + 1)
vertices_of_all_planes = []
for idx_plane in idxs_plane:
vertices = self.getVertices(
normal_vec,
idx_plane,
draw_location=draw_location,
thickness=thickness)
intersections = self.getPlaneBoxIntersections(st, vertices)
if len(intersections) >= 3:
if len(intersections) >= 4:
intersections = self.getOrderedIntersections(intersections)
vertices_of_all_planes.append(intersections)
elif not idx_plane:
normal_len = transform.get_vector_magnitude(normal_vec)
vertices = [
self.origin + normal_len *
transform.get_normalized_vector(v - self.origin)
for v in vertices
]
vertices_of_all_planes.append(vertices)
return vertices_of_all_planes