"""
Classes and functions to handle various shapes.
Copyright Schrodinger, LLC. All rights reserved."""
# Contributor: Thomas F. Hughes
import math
from collections import OrderedDict
from past.utils import old_div
import numpy
from scipy import constants
from schrodinger import structure
from schrodinger.structutils import color
from schrodinger.structutils import transform
from .msutils import add_or_update_bond_type
from .nano import xtal
_version = '$Revision 0.0 $'
ORIGIN = [0.0, 0.0, 0.0]
X_AXIS = numpy.array(transform.X_AXIS)
Y_AXIS = numpy.array(transform.Y_AXIS)
Z_AXIS = numpy.array(transform.Z_AXIS)
TEMPLATE_ELEMENT = 'C'
TEMPLATE_BOND_TYPE = structure.BondType.Single
# this parameter handles numerical precision for inf and nan
# cases
INF_NAN_THRESH = 1.0e-12
# see Face.intersectSegmentAndPlane where this parameter
# controls whether a line segment point (starting or ending)
# that lies near a given plane is considered to intersect
# with the plane, if positive then the plane is effectively
# pushed away from the segment and so the decision regarding
# intersection becomes less restrictive meaning that even
# though a formal intersection may be there it is not
# considered to be intersecting, if negative then the plane is
# effectively pulled toward the segment and so the decision
# is more restrictive meaning that even though a formal
# intersection is not there it is still considered to be
# intersecting, if zero the plane is not effectively moved
# and if the point is before the plane it is not intersecting
# while if the point is on or after the plane it is
# intersecting
DISTANCE_THRESH = 0.0001
SQRT2 = math.sqrt(2.0)
INVSQRT2 = old_div(1.0, SQRT2)
GOLDEN = constants.golden
INVGOLDEN = old_div(1.0, GOLDEN)
PI = constants.pi
PLATONIC = 'platonic'
ARCHIMEDEAN = 'archimedean'
PRISM = 'prism'
BASIC = 'basic'
POLYHEDRON_TYPES = [PLATONIC, ARCHIMEDEAN, PRISM]
LENGTH = 'length'
AREA = 'area'
COLINEAR_THRESH = 0.25
[docs]class Vertex(object):
"""
Class to manage vertices.
"""
[docs] def __init__(self, index, coordinates):
"""
Create a vertex instance.
:type index: int
:param index: the index of this vertex
:type coordinates: numpy.array
:param coordinates: the coordinates of this vertex
"""
self.index = index
self.coordinates = coordinates
[docs]class Edge(object):
"""
Class to manage edges.
"""
[docs] def __init__(self, index, vertex_1, vertex_2):
"""
Create an edge instance.
:type index: int
:param index: the index of this edge
:type vertex_1: `Vertex`
:param vertex_1: the first vertex in this edge
:type vertex_2: `Vertex`
:param vertex_2: the second vertex in this edge
"""
self.index = index
self.vertex_1 = vertex_1
self.vertex_2 = vertex_2
self.vector = None
self.length = self.getLength()
[docs] def getLength(self):
"""
Return the edge length.
:rtype: float
:return: the edge length in Ang.
"""
self.vector = self.vertex_2.coordinates - self.vertex_1.coordinates
return transform.get_vector_magnitude(self.vector)
[docs]class Face(object):
"""
Class to manage faces.
"""
[docs] def __init__(self, index, indices, points, num_unique):
"""
Create a face instance.
:type index: int
:param index: the index of this face
:type indices: list
:param indices: the indices of points making up this face
:type points: list of numpy.array
:param points: the points making up this face
:type num_unique: int
:param num_unique: the number of symmetry unique edges per face
"""
self.index = index
self.indices = indices
self.num_unique = num_unique
# in the following, ever_intersected is whether or not
# this face was ever intersected considering all point
# queries, while intersected is whether or not this face
# was intersected on the most recent point query,
# plane_intersected is whether or not the plane that this
# face resides in was intersected on the most recent point
# query
self.ever_intersected = False
self.intersected = False
self.plane_intersected = False
self.vertices = []
for index, point in zip(indices, points):
self.vertices.append(Vertex(index, point))
self.center = get_centroid(
[vertex.coordinates for vertex in self.vertices])
self.normal = None
self.area = None
self.edges = self.getEdges()
self.setNormal()
self.setArea()
self.reference_edges = self.getReferenceEdges(self.edges,
self.num_unique)
[docs] def getEdges(self):
"""
Return a list of Edge.
:rtype: list
:return: contains all Edge for this face
"""
edges = []
num_verts = len(self.vertices)
for i_idx in range(num_verts):
j_idx = (i_idx + 1) % num_verts
edges.append(
Edge(i_idx + 1, self.vertices[i_idx], self.vertices[j_idx]))
return edges
[docs] def setNormal(self):
"""
Set the normal to this face.
"""
# get two vertices in the plane of the face, pick
# the first two so that they are not colinear with
# the center
vertex_a, vertex_b = self.vertices[:2]
center_to_a = vertex_a.coordinates - self.center
center_to_b = vertex_b.coordinates - self.center
# the indexing is counterclockwise from the outside of the
# convex polyhedron from which the face came so cross in that way
# to get the outward normal
self.normal = numpy.cross(center_to_a, center_to_b)
self.normal = transform.get_normalized_vector(self.normal)
[docs] def setArea(self):
"""
Set the area of this face.
"""
# first rotate the polygon from 3D to 2D with the normal along the z-axis,
# translations do not matter here
rotation_matrix = transform.get_alignment_matrix(self.normal, Z_AXIS)
vertices_3d = [vertex.coordinates for vertex in self.vertices]
vertices_2d = []
for vertex_3d in vertices_3d:
vertex_2d = transform.transform_atom_coordinates(numpy.copy(vertex_3d), \
rotation_matrix)
vertices_2d.append(numpy.delete(vertex_2d, 2))
# second get the area
self.area = get_polygon_area(vertices_2d)
[docs] def getReferenceEdges(self, edges, num_unique):
"""
Return a list containing the num_unique number of unique edges.
These will be the reference edges used to orient the reference face
of the polyhedron.
:type edges: list
:param edges: all Edge objects for this face
:type num_unique: int
:param num_unique: the number of symmetry unique edges per face
:rtype: list
:return: unique Edge objects to serve as references
"""
# threshold in Ang. used to determine if two edges are equivalent
# by length
EDGE_THRESHOLD = 0.1
return get_reference_data(edges, LENGTH, num_unique, EDGE_THRESHOLD)
[docs] def intersectSegmentAndPlane(self,
line_start,
line_end,
inf_nan_thresh=INF_NAN_THRESH,
distance_thresh=DISTANCE_THRESH):
"""
Return the intersection point of the specified line segment and the plane
in which this face resides or return None if there is no intersection.
:type line_start: numpy.array
:param line_start: the starting point of the line segment
:type line_end: numpy.array
:param line_end: the ending point of the line segment
:type inf_nan_thresh: float
:param inf_nan_thresh: this parameter handles numerical precision
for inf and nan cases
:type distance_thresh: float
:param distance_thresh: this parameter controls how the intersection of
line segment end-points and a plane are handled for cases where one of the
end-points lies in (or near) the plane (see the comment near the module
level constant DISTANCE_THRESH)
:rtype: numpy.array or None
:return: the single point of intersection along the line segment
or None if there is no intersection
"""
line = line_end - line_start
numerator = numpy.dot(self.center - line_start, self.normal)
denominator = numpy.dot(line, self.normal)
# (1) the distance is infinite, i.e. +-x/0 for non-zero x, if the line
# segment is outside of the plane and parallel to it, this includes a
# line segment of zero length, the sign of the infinity depends on the
# numerator and thus whether the segment is above or below the plane
# (regardless of which way the line segment vector is pointing) (note
# that for use with nanoparticle we will always be talking about the
# positive infinity case because the line segment will be below the
# plane given that one of the end points is the convex polyhedron
# center and the face normals always point outwards), with infinite
# distance the intersection point is at infinity and so we return
# with no intersection
# (2) the distance is NaN, i.e. 0/0, if the line segment is inside the
# plane, i.e. parallel to it and inside and thus intersecting everywhere,
# this includes a line segment of zero length and again regardless
# of which way the line segment vector is pointing, because it
# intersects everywhere by convention we just return the starting point
# (note that for use with nanoparticle we will never see this case
# because one of the line segment end points is the convex polyhedron
# center which necessitates a non-zero numerator)
if abs(denominator) < inf_nan_thresh:
if abs(numerator) < inf_nan_thresh:
return line_start
else:
return None
# the following line can throw a RuntimeWarning, read the comments
# within this method to see how this is being handled (note that for
# use with nanoparticle we will never see a -/+ (distance <= 0) case
# because one of the line segment end points is the convex polyhedron
# center which necessitates that the numerator be positive if the
# denominator is positive, there are however +/+ (distance >= 0),
# -/- (distance >= 0), and +/- (distance <= 0) cases)
distance = old_div(numerator, denominator)
# (3) for use with nanoparticle the distance works like this, take
# a query point (an atom), this is the line_start, that is far above
# the plane (plane containing the polyhedron face in question) (by
# above we mean in the direction of the plane normal), if infinitely
# far away from the plane then the distance is 1 and the intersection
# point is the line_end, i.e. polyhedron center, then as this point
# approaches the plane the distance approaches 0/x, i.e. zero for x
# non-zero, until it is on the plane at which point the distance is
# 0 and the intersection point is the query point position (the line
# start), so for on or above the plane we have a positive distance in
# [0, 1), these are actually the only points that can intersect
#
# nevertheless we will continue the analysis, once our query point is
# below the plane the distance is negative, +/- case, up until the
# line segment becomes perpendicular to the plane normal at which
# point the denominator goes to negative zero (numerator is positive
# non-zero) and the distance is negative infinity (this is caught
# above), once the line segment again becomes out of the plane
# the distance becomes positive infinity, x/+0 for positive non-zero
# x, and then as the point becomes infinitely below the plane we go
# from a distance of infinity to a distance of 1 again, so in summary
# as the point moves from far above the plane, through the plane,
# through the line segment end point, to far below the plane,
# the distances go from 1 to 0 to -infinity to +infinity to 1 and
# intersections can only occur on [0, 1), so we will want to return
# with an intersection there (use a threshold to classify cases near
# the plane)
#
# for use outside of nanoparticle, i.e. the line segment end point
# may be on the plane, the analysis will be done similarly to above
# except swap the starting and ending line segment points so that
# it is the end point that is moving through the plane in the analysis,
# in this case the distances go from 0 to 1 (as you move the ending
# point toward the plane from infinity), then from 1 to +infinity
# (as the point moved to the starting point) then to -infinity (as
# the point moves through the starting point) then from -infinity
# to zero as the end point becomes infinitely below the plane,
#
# note that the case where the line segment end point can be
# on the plane introduces an asymmetry with respect to the
# long range behavior of the distances, if it can't as in
# nanoparticle then the intersecting distances will be [0, 1)
# while if it can the intersecting distances will be (0, 1]
# (4) all other cases do not intersect
# note that for nanoparticle the line_end is the polyhedron center
# which is not on the plane and so the upper distance threshold is
# zero
if abs(numpy.dot(self.normal, line_end - self.center)) < inf_nan_thresh:
lower_distance_thresh = 0
upper_distance_thresh = distance_thresh
else:
lower_distance_thresh = distance_thresh
upper_distance_thresh = 0
if lower_distance_thresh <= distance <= 1 - upper_distance_thresh:
return line_start + distance * line
else:
return None
[docs] def insideFace(self, point):
"""
Return true if the query point lies on or inside the boundaries
of this face, false otherwise.
:type point: numpy.array
:param point: a query point
:rtype: bool
:return: true if the query point lies on or inside the face boundaries,
false otherwise
"""
DISTANCE_THRESHOLD = 0.0000001
ANGLE_THRESHOLD = 0.0001
# handle query points that are close to the vertices
for vertex in self.vertices:
if transform.get_vector_magnitude(vertex.coordinates - point) < \
DISTANCE_THRESHOLD:
return True
# handle all other points including those on edges
first_vec = self.vertices[0].coordinates - point
old_vec = first_vec
revolution = 0.0
for vertex in self.vertices[1:]:
new_vec = vertex.coordinates - point
revolution += transform.get_angle_between_vectors(new_vec, old_vec)
old_vec = new_vec
revolution += transform.get_angle_between_vectors(first_vec, old_vec)
if abs(revolution - 2 * PI) < ANGLE_THRESHOLD:
return True
[docs]class ConvexPolyhedron(object):
"""
Class to manage a convex polyhedron.
"""
[docs] def __init__(self, params, center, ref_face_idx, ref_face_normal_along,
ref_edge_idx, ref_edge_along):
"""
Create an instance.
:type params: list
:param params: list of floating point parameters defining the polyhedron
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
if self.TYPE != PRISM:
self.scale = params[0]
else:
self.scale = 1
self.a_param, self.b_param, self.c_param, \
self.alpha_param, self.beta_param, self.gamma_param = params
# get lattice vectors and original vertices
self.a_vec, self.b_vec, self.c_vec = xtal.get_lattice_vectors(
*params)
origin = numpy.array(xtal.ParserWrapper.ORIGIN)
self.VERTICES = get_parallelepiped_vertices(origin, self.a_vec,
self.b_vec, self.c_vec)
self.center = numpy.array(center)
self.ref_face_idx = ref_face_idx
self.ref_face_normal_along = ref_face_normal_along
self.ref_edge_idx = ref_edge_idx
self.ref_edge_along = ref_edge_along
# get references for alignment
self.vertices = self.getVertices(self.VERTICES, self.scale)
self.faces = self.getFaces(self.vertices, self.INDICES,
self.NUM_UNIQUE_EDGES)
self.reference_faces = self.getReferenceFaces(self.faces, \
self.NUM_UNIQUE_FACES)
reference_face = self.reference_faces[self.ref_face_idx - 1]
self.reference_edges = reference_face.reference_edges
reference_edge = self.reference_edges[self.ref_edge_idx - 1]
# do the alignment and translation to arrive at the final
# scaled, translated, and rotated set of vertices
self.vertices = self.alignPolyhedron(self.vertices, \
reference_face.normal, self.ref_face_normal_along, \
reference_edge.vector, self.ref_edge_along)
self.translateVertices(self.vertices, self.center)
[docs] def updateShape(self, vertices):
"""
Update the shape object using the provided vertices.
:type vertices: list
:param vertices: list of Vertex with new coordinates
"""
self.vertices = vertices
self.center = get_centroid(
[vertex.coordinates for vertex in self.vertices])
# update reference faces and edges (update them rather than redetermine
# them because they involve sorting by double precision data which for
# equivalent data can suffer from numerical precision issues which result
# in different orderings before and after rotation)
# collect reference edge indices for all faces
ref_edge_inds = []
for face in self.faces:
ref_edge_inds.append(
[edge.index - 1 for edge in face.reference_edges])
# collect reference face indices
ref_face_inds = []
for face in self.reference_faces:
ref_face_inds.append(face.index - 1)
# get all faces
self.faces = self.getFaces(self.vertices, self.INDICES,
self.NUM_UNIQUE_EDGES)
# update reference edges for all faces
for face, inds in zip(self.faces, ref_edge_inds):
face.reference_edges = [face.edges[x] for x in inds]
# update reference faces
self.reference_faces = [self.faces[x] for x in ref_face_inds]
# update reference edges
self.reference_edges = self.reference_faces[self.ref_face_idx -
1].reference_edges
# get template
self.template = self.makeTemplate()
[docs] def getVertices(self, vertices, scale=1.0):
"""
Create vertex data for the polyhedron, including an option to
scale all vertices using a multiplicative factor.
:type vertices: list
:param vertices: list of numpy.array triples specifying
the vertices of this polyhedron
:type scale: float
:param scale: multiplicative factor used to scale all vertices
:rtype: list
:return: a list of scaled Vertex
"""
scaled_vertices = []
for index, vertex in enumerate(vertices, 1):
scaled_vertices.append(Vertex(index, scale * vertex))
return scaled_vertices
[docs] def getFaces(self, vertices, all_indices, num_unique):
"""
Create face data for the polyhedron.
:type vertices: list
:param vertices: list of scaled Vertex
:type all_indices: list
:param all_indices: contains sub-lists specifying the vertex
indices for each face
:type num_unique: int
:param num_unique: the number of symmetry unique edges per face
:rtype: list
:return: a list of Face
"""
faces = []
for findex, indices in enumerate(all_indices, 1):
face = [vertices[vindex - 1].coordinates for vindex in indices]
faces.append(Face(findex, indices, face, num_unique))
return faces
[docs] def translateVertices(self, vertices, vector):
"""
Translate the vertices by adding the specified vector.
:type vertices: list
:param vertices: list of scaled Vertex
:type vector: numpy.array triple
:param vector: the vector used to translate the vertices
"""
translated_vertices = []
for vertex in vertices:
translated_vertices.append(Vertex(vertex.index, \
vertex.coordinates + vector))
self.updateShape(translated_vertices)
[docs] def getSurfaceArea(self, faces):
"""
Return the surface area of the polyhedron.
:type faces: list
:param faces: all Face objects for this polyhedron
:rtype: float
:return: the surface area in Ang.^2
"""
return sum((x.area for x in faces))
[docs] def getReferenceFaces(self, faces, num_unique):
"""
Return a list containing the num_unique number of unique faces.
These will be the reference faces used to orient the polyhedron.
:type faces: list
:param faces: all Face objects for this polyhedron
:type num_unique: int
:param num_unique: the number of symmetry unique faces per polyhedron
:rtype: list
:return: unique Face objects to serve as references
"""
# threshold in Ang.^2 used to determine if two faces are equivalent
# by area
AREA_THRESHOLD = 0.1
return get_reference_data(faces, AREA, num_unique, AREA_THRESHOLD)
[docs] def alignPolyhedron(self, vertices, face_vector, face_along, \
edge_vector, edge_along):
"""
Return the polyhedron vertices rotated so as to align the face and edge
vectors.
:type vertices: list
:param vertices: list of Vertex
:type face_vector: numpy.array triple
:param face_vector: the normal of the reference face to be rotated
:type face_along: numpy.array triple
:param face_along: the vector onto which the face normal will be rotated
:type edge_vector: numpy.array triple
:param edge_vector: the vector of the reference edge to be rotated
:type edge_along: numpy.array triple
:param edge_along: the vector onto which the edge vector will be rotated, it
is actually this vector's component that is perpendicular to face_along that
is used in the alignment, this is to safeguard against the case where
edge_along and face_along are not perpendicular, if they are the same
an exception is raised
:rtype: list
:return: list of rotated Vertex
"""
def return_rotated(in_vertices, rot_matrix):
out_vertices = []
for in_vertex in in_vertices:
out_vertex = transform.transform_atom_coordinates(numpy.copy(in_vertex), \
rot_matrix)
out_vertices.append(out_vertex)
return out_vertices
def rotate_collections(verts, face, edge, matrix):
rot_verts = return_rotated(verts, matrix)
rot_face = return_rotated([face], matrix)[0]
rot_edge = return_rotated([edge], matrix)[0]
return rot_verts, rot_face, rot_edge
# for convenience do not use the edge_along vector directly but rather
# use it's component that is perpendicular to the face_along vector
unit_face_along = transform.get_normalized_vector(face_along)
unit_edge_along = transform.get_normalized_vector(edge_along)
perp_edge_along = unit_edge_along - numpy.dot(unit_edge_along, \
unit_face_along) * unit_face_along
if transform.get_vector_magnitude(perp_edge_along) < COLINEAR_THRESH:
no_perp_msg = (
'You have specified that the vector onto which '
'the reference edge vector will be aligned is colinear with '
'the vector onto which the reference face normal will be '
'aligned. These two target vectors must not be colinear.')
raise ValueError(no_perp_msg)
# first rotate the face vector onto the target face vector, we don't
# care about any of the other degrees of freedom just yet
face_matrix = transform.get_alignment_matrix(face_vector, face_along)
# update important parameters
rot_vertices, rot_face_vector, rot_edge_vector = \
rotate_collections([vertex.coordinates for vertex in vertices], \
face_vector, edge_vector, face_matrix)
# second rotate the updated edge vector onto the target edge vector,
# do not use transform.get_alignment_matrix again as it can destroy
# the previous alignment as you aren't guaranteed anything about
# the axis about which the alignment is done (see MATSCI-1971)
angle = transform.get_angle_between_vectors(rot_edge_vector,
perp_edge_along)
unit_rot_face_vector = transform.get_normalized_vector(rot_face_vector)
if angle == 0 or angle == PI:
sign = 1
else:
unit_cross = transform.get_normalized_vector(
numpy.cross(rot_edge_vector, perp_edge_along))
if numpy.dot(unit_cross, unit_rot_face_vector) > 0:
sign = 1
else:
sign = -1
edge_matrix = transform.get_rotation_matrix(sign * unit_rot_face_vector,
angle)
# update important parameters
rot_vertices, rot_face_vector, rot_edge_vector = \
rotate_collections(rot_vertices, rot_face_vector, rot_edge_vector, edge_matrix)
return self.getVertices(rot_vertices)
[docs] def makeTemplate(self):
"""
Create a template, i.e. structure object, for this convex polyhedron.
:rtype: `schrodinger.structure.Structure`
:return: the template structure
"""
template = structure.create_new_structure(num_atoms=len(self.vertices))
template.title = self.NAME
for aatom, vertex in zip(template.atom, self.vertices):
aatom.element = TEMPLATE_ELEMENT
aatom.xyz = vertex.coordinates
for face in self.faces:
first_index = face.vertices[0].index
last_index = first_index
for vertex in face.vertices[1:]:
current_index = vertex.index
add_or_update_bond_type(template, last_index, current_index,
TEMPLATE_BOND_TYPE)
last_index = current_index
add_or_update_bond_type(template, last_index, first_index,
TEMPLATE_BOND_TYPE)
# color the template
color.apply_color_scheme(template, 'element')
return template
[docs] def addPointsToTemplate(self, points):
"""
Add the specified points to this convex polyhedron's template.
:type points: list of numpy.array
:param points: contains the points to be added to the template
for this convex polyhedron
"""
for point in points:
self.template.addAtom(TEMPLATE_ELEMENT, *point)
[docs] def addAlignmentAxesToTemplate(self):
"""
Add unit vectors that mark the primary and secondary alignment axes
to the template.
"""
natoms = self.template.atom_total
face = self.reference_faces[self.ref_face_idx - 1]
unit_normal = transform.get_normalized_vector(face.normal)
self.addPointsToTemplate([face.center, face.center + unit_normal])
self.template.addBond(natoms + 1, natoms + 2, TEMPLATE_BOND_TYPE)
natoms = self.template.atom_total
edge = face.reference_edges[self.ref_edge_idx - 1]
unit_edge = transform.get_normalized_vector(edge.vector)
self.addPointsToTemplate([edge.vertex_2.coordinates + unit_edge])
self.template.addBond(edge.vertex_2.index, natoms + 1,
TEMPLATE_BOND_TYPE)
# color the template
color.apply_color_scheme(self.template, 'element')
[docs] def addNormalsToTemplate(self):
"""
Add the normals of this convex polyhedron to its template.
"""
natoms = self.template.atom_total
for face in self.faces:
self.addPointsToTemplate([face.center, face.center + face.normal])
self.template.addBond(natoms + 1, natoms + 2, TEMPLATE_BOND_TYPE)
natoms += 2
[docs] def getSegmentPlaneIntersections(self, line_start, line_end):
"""
Return a two lists (1) of points where the given line segment intersects
the planes containing this polyhedron's faces and (2) the centers of the faces
whose planes are being intersected.
:type line_start: numpy.array
:param line_start: the start of the line segment
:type line_end: numpy.array
:param line_end: the end of the line segment
:rtype: two lists of numpy.array
:return: the intersection points and the centers of faces whose planes are
being intersected
"""
intersections = []
centers = []
for face in self.faces:
intersection = face.intersectSegmentAndPlane(line_start, line_end)
if intersection is not None:
intersections.append(intersection)
centers.append(face.center)
face.plane_intersected = True
if face.insideFace(intersection):
face.ever_intersected = True
face.intersected = True
else:
face.intersected = False
else:
face.intersected = False
face.plane_intersected = False
return intersections, centers
[docs] def addSegmentPlaneIntersectionsToTemplate(self, line_start, line_end):
"""
Add to this polyhedron's template the intersection points of the specified line
segment and the planes containing the faces. Also draw connections between
these points and the centers of the faces containing the planes that are being
intersected.
:type line_start: numpy.array
:param line_start: the start of the line segment
:type line_end: numpy.array
:param line_end: the end of the line segment
"""
intersections, centers = self.getSegmentPlaneIntersections(
line_start, line_end)
if intersections:
natoms = self.template.atom_total
for intersection, center in zip(intersections, centers):
self.addPointsToTemplate([intersection, center])
self.template.addBond(natoms + 1, natoms + 2,
TEMPLATE_BOND_TYPE)
natoms += 2
[docs] def pointInside(self, point):
"""
Return True if the query point is either on or inside of
this convex polyhedron.
:type point: numpy.array
:param point: the point in question
:rtype: bool
:return: True if the point in question is either on or inside
this polyhedron, False otherwise
"""
# to test the membership of a point in this polyhedron just create a line
# segment from the center of the polyhedron to the query point and test
# to see if any of the planes containing the polyhedrons faces intersect,
# if yes then the point is outside, if no then the point is inside
intersections, centers = self.getSegmentPlaneIntersections(
point, self.center)
return not intersections
[docs] def allFacesIntersected(self):
"""
Return True if all faces, not the planes containing those faces but
the actual faces, of this convex polyhedron have been intersected
at least once given all pointInside queries performed thus far.
:rtype: bool
:return: True if all faces have been intersected, False otherwise
"""
return all([face.ever_intersected for face in self.faces])
# Note the following regarding the computation of volume, surface area,
# etc. properties of the following polyhedral shapes. There is an entire
# field of mathematics exactly for this and our goal is not to provide
# an exhaustive API for this but rather support only those properties
# that we will need in the physical model. For example, the platonic
# solids have general equations for these properties but involve the
# computation of intermediate quantities that are of little physical
# interest. And the properties of other polyhedra, for example the
# cubeoctahedron, must be obtained by deconstructing the shape in terms
# of its platonic ones and then combining the properties of those.
# Because of this we are going to simply use the specific formulas
# for computing these properties.
[docs]class Cube(ConvexPolyhedron):
"""
Class to manage a cube.
"""
NAME = 'cube'
DEFAULT = [5.0]
DESCRIPTION = ['edge length']
TYPE = PLATONIC
# following final VERTICES are for a cube centered
# at the origin and of length 1 and arbitrarily oriented
# as by default it will be rotated so that the normal of
# the largest face lies along the z-axis while the longest
# edge of that face lies along the x-axis
VERTICES = [[1.0, 1.0, 1.0], [1.0, 1.0, -1.0], [1.0, -1.0, 1.0],
[1.0, -1.0, -1.0], [-1.0, 1.0, 1.0], [-1.0, 1.0, -1.0],
[-1.0, -1.0, 1.0], [-1.0, -1.0, -1.0]]
VERTICES = [0.5 * numpy.array(vertex) for vertex in VERTICES]
# here are the indicies of the faces, the order of the sets is irrelevant,
# the order within each set is also irrelevant with the exception that it
# should be counter-clockwise so that the normal point to the outside
# of the cube
INDICES = [[1, 3, 4, 2], [5, 1, 2, 6], [7, 5, 6, 8], [3, 7, 8, 4],
[1, 5, 7, 3], [4, 8, 6, 2]]
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces
NUM_UNIQUE_EDGES = 1
# this is the number of symmetry unique faces per polyhedron
NUM_UNIQUE_FACES = 1
[docs] def __init__(self, scale, center=ORIGIN, ref_face_idx=1, ref_face_normal_along=Z_AXIS, \
ref_edge_idx=1, ref_edge_along=X_AXIS):
"""
Create an instance.
:type scale: float
:param scale: multiplicative scaling factor
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
params = [scale]
super(self.__class__, self).__init__(params, center, ref_face_idx,
ref_face_normal_along,
ref_edge_idx, ref_edge_along)
self.circumradius = self.getCircumRadius(self.scale)
self.volume = self.getVolume(self.scale)
self.surface_area = self.getSurfaceArea(self.faces)
[docs] def getEdgeLength(self, circumradius):
"""
Return the edge length.
:type circumradius: float
:param circumradius: circumradius in Angstrom
:rtype: float
:return: edge length in Angstrom
"""
factor = old_div(math.sqrt(3.0), 2.0)
factor = old_div(1.0, factor)
return factor * circumradius
[docs] def getCircumRadius(self, length):
"""
Return the circumradius.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: circumradius in Angstrom
"""
factor = old_div(math.sqrt(3.0), 2.0)
return factor * length
[docs] def getVolume(self, length):
"""
Return the volume.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: volume in cubic Angstrom
"""
return length**3
[docs]class Tetrahedron(ConvexPolyhedron):
"""
Class to manage a tetrahedron.
"""
NAME = 'tetrahedron'
DEFAULT = [10.0]
DESCRIPTION = ['edge length']
TYPE = PLATONIC
# following final VERTICES are for a tetrahedron centered
# at the origin and of length 1 and arbitrarily oriented
# as by default it will be rotated so that the normal of
# the largest face lies along the z-axis while the longest
# edge of that face lies along the x-axis
VERTICES = [[1.0, 0.0, -1.0 * INVSQRT2], [-1.0, 0.0, -1.0 * INVSQRT2],
[0.0, 1.0, INVSQRT2], [0.0, -1.0, INVSQRT2]]
VERTICES = [0.5 * numpy.array(vertex) for vertex in VERTICES]
# here are the indicies of the faces, the order of the sets is irrelevant,
# the order within each set is also irrelevant with the exception that it
# should be counter-clockwise so that the normal point to the outside
# of the tetrahedron
INDICES = [[1, 4, 2], [3, 4, 1], [2, 4, 3], [1, 2, 3]]
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces
NUM_UNIQUE_EDGES = 1
# this is the number of symmetry unique faces per polyhedron
NUM_UNIQUE_FACES = 1
[docs] def __init__(self, scale, center=ORIGIN, ref_face_idx=1, ref_face_normal_along=Z_AXIS, \
ref_edge_idx=1, ref_edge_along=X_AXIS):
"""
Create an instance.
:type scale: float
:param scale: multiplicative scaling factor
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
params = [scale]
super(self.__class__, self).__init__(params, center, ref_face_idx,
ref_face_normal_along,
ref_edge_idx, ref_edge_along)
self.circumradius = self.getCircumRadius(self.scale)
self.volume = self.getVolume(self.scale)
self.surface_area = self.getSurfaceArea(self.faces)
[docs] def getEdgeLength(self, circumradius):
"""
Return the edge length.
:type circumradius: float
:param circumradius: circumradius in Angstrom
:rtype: float
:return: edge length in Angstrom
"""
factor = math.sqrt(old_div(3.0, 8.0))
factor = old_div(1.0, factor)
return factor * circumradius
[docs] def getCircumRadius(self, length):
"""
Return the circumradius.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: circumradius in Angstrom
"""
factor = math.sqrt(old_div(3.0, 8.0))
return factor * length
[docs] def getVolume(self, length):
"""
Return the volume.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: volume in cubic Angstrom
"""
factor = 6.0 * SQRT2
factor = old_div(1.0, factor)
return factor * length**3
[docs]class Octahedron(ConvexPolyhedron):
"""
Class to manage an octahedron.
"""
NAME = 'octahedron'
DEFAULT = [8.0]
DESCRIPTION = ['edge length']
TYPE = PLATONIC
# following final VERTICES are for an octahedron centered
# at the origin and of length 1 and arbitrarily oriented
# as by default it will be rotated so that the normal of
# the largest face lies along the z-axis while the longest
# edge of that face lies along the x-axis
VERTICES = [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0],
[0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0]]
VERTICES = [INVSQRT2 * numpy.array(vertex) for vertex in VERTICES]
# here are the indicies of the faces, the order of the sets is irrelevant,
# the order within each set is also irrelevant with the exception that it
# should be counter-clockwise so that the normal point to the outside
# of the octahedron
INDICES = [[6, 2, 3], [4, 2, 6], [5, 2, 4], [3, 2, 5], [6, 3, 1], [4, 6, 1],
[5, 4, 1], [3, 5, 1]]
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces
NUM_UNIQUE_EDGES = 1
# this is the number of symmetry unique faces per polyhedron
NUM_UNIQUE_FACES = 1
[docs] def __init__(self, scale, center=ORIGIN, ref_face_idx=1, ref_face_normal_along=Z_AXIS, \
ref_edge_idx=1, ref_edge_along=X_AXIS):
"""
Create an instance.
:type scale: float
:param scale: multiplicative scaling factor
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
params = [scale]
super(self.__class__, self).__init__(params, center, ref_face_idx,
ref_face_normal_along,
ref_edge_idx, ref_edge_along)
self.circumradius = self.getCircumRadius(self.scale)
self.volume = self.getVolume(self.scale)
self.surface_area = self.getSurfaceArea(self.faces)
[docs] def getEdgeLength(self, circumradius):
"""
Return the edge length.
:type circumradius: float
:param circumradius: circumradius in Angstrom
:rtype: float
:return: edge length in Angstrom
"""
factor = old_div(SQRT2, 2.0)
factor = old_div(1.0, factor)
return factor * circumradius
[docs] def getCircumRadius(self, length):
"""
Return the circumradius.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: circumradius in Angstrom
"""
factor = old_div(SQRT2, 2.0)
return factor * length
[docs] def getVolume(self, length):
"""
Return the volume.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: volume in cubic Angstrom
"""
factor = old_div(SQRT2, 3.0)
return factor * length**3
[docs]class Dodecahedron(ConvexPolyhedron):
"""
Class to manage a dodecahedron.
"""
NAME = 'dodecahedron'
DEFAULT = [3.0]
DESCRIPTION = ['edge length']
TYPE = PLATONIC
# following final VERTICES are for a dodecahedron centered
# at the origin and of length 1 and arbitrarily oriented
# as by default it will be rotated so that the normal of
# the largest face lies along the z-axis while the longest
# edge of that face lies along the x-axis
VERTICES = [[1.0, 1.0, 1.0], [1.0, 1.0, -1.0], [1.0, -1.0, 1.0],
[1.0, -1.0, -1.0], [-1.0, 1.0, 1.0], [-1.0, 1.0, -1.0],
[-1.0, -1.0, 1.0], [-1.0, -1.0, -1.0], [0.0, INVGOLDEN, GOLDEN],
[0.0, INVGOLDEN,
-1.0 * GOLDEN], [0.0, -1.0 * INVGOLDEN, GOLDEN],
[0.0, -1.0 * INVGOLDEN, -1.0 * GOLDEN],
[INVGOLDEN, GOLDEN, 0.0], [INVGOLDEN, -1.0 * GOLDEN, 0.0],
[-1.0 * INVGOLDEN, GOLDEN,
0.0], [-1.0 * INVGOLDEN, -1.0 * GOLDEN, 0.0],
[GOLDEN, 0.0, INVGOLDEN], [GOLDEN, 0.0, -1.0 * INVGOLDEN],
[-1.0 * GOLDEN, 0.0, INVGOLDEN],
[-1.0 * GOLDEN, 0.0, -1.0 * INVGOLDEN]]
VERTICES = [0.5 * GOLDEN * numpy.array(vertex) for vertex in VERTICES]
# here are the indicies of the faces, the order of the sets is irrelevant,
# the order within each set is also irrelevant with the exception that it
# should be counter-clockwise so that the normal point to the outside
# of the dodecahedron
INDICES = [[6, 15, 13, 2, 10], [20, 6, 10, 12, 8], [19, 20, 8, 16, 7],
[5, 19, 7, 11, 9], [15, 5, 9, 1, 13], [20, 19, 5, 15, 6],
[10, 2, 18, 4, 12], [8, 12, 4, 14, 16], [7, 16, 14, 3, 11],
[9, 11, 3, 17, 1], [13, 1, 17, 18, 2], [4, 18, 17, 3, 14]]
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces
NUM_UNIQUE_EDGES = 1
# this is the number of symmetry unique faces per polyhedron
NUM_UNIQUE_FACES = 1
[docs] def __init__(self, scale, center=ORIGIN, ref_face_idx=1, ref_face_normal_along=Z_AXIS, \
ref_edge_idx=1, ref_edge_along=X_AXIS):
"""
Create an instance.
:type scale: float
:param scale: multiplicative scaling factor
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
params = [scale]
super(self.__class__, self).__init__(params, center, ref_face_idx,
ref_face_normal_along,
ref_edge_idx, ref_edge_along)
self.circumradius = self.getCircumRadius(self.scale)
self.volume = self.getVolume(self.scale)
self.surface_area = self.getSurfaceArea(self.faces)
[docs] def getEdgeLength(self, circumradius):
"""
Return the edge length.
:type circumradius: float
:param circumradius: circumradius in Angstrom
:rtype: float
:return: edge length in Angstrom
"""
factor = math.sqrt(3.0) * GOLDEN / 2.0
factor = old_div(1.0, factor)
return factor * circumradius
[docs] def getCircumRadius(self, length):
"""
Return the circumradius.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: circumradius in Angstrom
"""
factor = math.sqrt(3.0) * GOLDEN / 2.0
return factor * length
[docs] def getVolume(self, length):
"""
Return the volume.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: volume in cubic Angstrom
"""
factor = old_div((15.0 + 7.0 * math.sqrt(5.0)), 4.0)
return factor * length**3
[docs]class Icosahedron(ConvexPolyhedron):
"""
Class to manage an icosahedron.
"""
NAME = 'icosahedron'
DEFAULT = [5.0]
DESCRIPTION = ['edge length']
TYPE = PLATONIC
# following final VERTICES are for an icosahedron centered
# at the origin and of length 1 and arbitrarily oriented
# as by default it will be rotated so that the normal of
# the largest face lies along the z-axis while the longest
# edge of that face lies along the x-axis
VERTICES = [[0.0, 1.0, GOLDEN], [0.0, 1.0, -1.0 * GOLDEN],
[0.0, -1.0, GOLDEN], [0.0, -1.0, -1.0 * GOLDEN],
[1.0, GOLDEN, 0.0], [1.0, -1.0 * GOLDEN, 0.0],
[-1.0, GOLDEN, 0.0], [-1.0, -1.0 * GOLDEN, 0.0],
[GOLDEN, 0.0, 1.0], [GOLDEN, 0.0, -1.0],
[-1.0 * GOLDEN, 0.0, 1.0], [-1.0 * GOLDEN, 0.0, -1.0]]
VERTICES = [0.5 * numpy.array(vertex) for vertex in VERTICES]
# here are the indicies of the faces, the order of the sets is irrelevant,
# the order within each set is also irrelevant with the exception that it
# should be counter-clockwise so that the normal point to the outside
# of the icosahedron
INDICES = [[2, 4, 12], [2, 12, 7], [2, 7, 5], [2, 5, 10], [2, 10, 4],
[4, 10, 6], [4, 6, 8], [12, 4, 8], [12, 8, 11], [7, 12, 11],
[7, 11, 1], [5, 7, 1], [5, 1, 9], [10, 5, 9], [10, 9, 6],
[8, 6, 3], [11, 8, 3], [1, 11, 3], [9, 1, 3], [6, 9, 3]]
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces
NUM_UNIQUE_EDGES = 1
# this is the number of symmetry unique faces per polyhedron
NUM_UNIQUE_FACES = 1
[docs] def __init__(self, scale, center=ORIGIN, ref_face_idx=1, ref_face_normal_along=Z_AXIS, \
ref_edge_idx=1, ref_edge_along=X_AXIS):
"""
Create an instance.
:type scale: float
:param scale: multiplicative scaling factor
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
params = [scale]
super(self.__class__, self).__init__(params, center, ref_face_idx,
ref_face_normal_along,
ref_edge_idx, ref_edge_along)
self.circumradius = self.getCircumRadius(self.scale)
self.volume = self.getVolume(self.scale)
self.surface_area = self.getSurfaceArea(self.faces)
[docs] def getEdgeLength(self, circumradius):
"""
Return the edge length.
:type circumradius: float
:param circumradius: circumradius in Angstrom
:rtype: float
:return: edge length in Angstrom
"""
factor = old_div(math.sqrt(GOLDEN * math.sqrt(5.0)), 2.0)
factor = old_div(1.0, factor)
return factor * circumradius
[docs] def getCircumRadius(self, length):
"""
Return the circumradius.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: circumradius in Angstrom
"""
factor = old_div(math.sqrt(GOLDEN * math.sqrt(5.0)), 2.0)
return factor * length
[docs] def getVolume(self, length):
"""
Return the volume.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: volume in cubic Angstrom
"""
factor = 5.0 * (3.0 + math.sqrt(5.0)) / 12.0
return factor * length**3
[docs]class Cubeoctahedron(ConvexPolyhedron):
"""
Class to manage a cubeoctahedron.
"""
NAME = 'cubeoctahedron'
DEFAULT = [5.0]
DESCRIPTION = ['edge length']
TYPE = ARCHIMEDEAN
# following final VERTICES are for a cubeoctahedron centered
# at the origin and of length 1 and arbitrarily oriented
# as by default it will be rotated so that the normal of
# the largest face lies along the z-axis while the longest
# edge of that face lies along the x-axis
VERTICES = [[1.0, 1.0, 0.0], [1.0, -1.0, 0.0], [-1.0, 1.0, 0.0],
[-1.0, -1.0, 0.0], [1.0, 0.0, 1.0], [1.0, 0.0, -1.0],
[-1.0, 0.0, 1.0], [-1.0, 0.0, -1.0], [0.0, 1.0, 1.0],
[0.0, 1.0, -1.0], [0.0, -1.0, 1.0], [0.0, -1.0, -1.0]]
VERTICES = [INVSQRT2 * numpy.array(vertex) for vertex in VERTICES]
# here are the indicies of the faces, the order of the sets is irrelevant,
# the order within each set is also irrelevant with the exception that it
# should be counter-clockwise so that the normal point to the outside
# of the cubeoctahedron
INDICES = [[2, 5, 11], [1, 9, 5], [5, 9, 7, 11], [11, 7, 4], [11, 4, 12, 2],
[2, 12, 6], [2, 6, 1, 5], [1, 10, 3, 9], [9, 3, 7], [7, 3, 8, 4],
[4, 8, 12], [12, 8, 10, 6], [6, 10, 1], [10, 8, 3]]
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces
NUM_UNIQUE_EDGES = 1
# this is the number of symmetry unique faces per polyhedron
NUM_UNIQUE_FACES = 2
[docs] def __init__(self, scale, center=ORIGIN, ref_face_idx=1, ref_face_normal_along=Z_AXIS, \
ref_edge_idx=1, ref_edge_along=X_AXIS):
"""
Create an instance.
:type scale: float
:param scale: multiplicative scaling factor
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
params = [scale]
super(self.__class__, self).__init__(params, center, ref_face_idx,
ref_face_normal_along,
ref_edge_idx, ref_edge_along)
self.circumradius = self.getCircumRadius(self.scale)
self.volume = self.getVolume(self.scale)
self.surface_area = self.getSurfaceArea(self.faces)
[docs] def getEdgeLength(self, circumradius):
"""
Return the edge length.
:type circumradius: float
:param circumradius: circumradius in Angstrom
:rtype: float
:return: edge length in Angstrom
"""
return circumradius
[docs] def getCircumRadius(self, length):
"""
Return the circumradius.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: circumradius in Angstrom
"""
return length
[docs] def getVolume(self, length):
"""
Return the volume.
:type length: float
:param length: length in Angstrom
:rtype: float
:return: volume in cubic Angstrom
"""
factor = 5.0 * SQRT2 / 3.0
return factor * length**3
[docs]class Parallelepiped(ConvexPolyhedron):
"""
Class to manage a parallelepiped.
"""
NAME = 'parallelepiped'
DEFAULT = [5.0, 10.0, 12.0, 60.0, 45.0, 80.0]
DESCRIPTION = ['edge a length', 'edge b length', 'edge c length', \
'edge-edge b-c angle', 'edge-edge a-c angle', 'edge-edge a-b angle']
TYPE = PRISM
# here are the indicies of the faces, the order of the sets is irrelevant,
# the order within each set is also irrelevant with the exception that it
# should be counter-clockwise so that the normal point to the outside
# of the parallelepiped
INDICES = [[1, 4, 3, 2], [5, 1, 2, 6], [4, 8, 7, 3], [8, 5, 6, 7],
[8, 4, 1, 5], [6, 2, 3, 7]]
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces
NUM_UNIQUE_EDGES = 2
# this is the number of symmetry unique faces per polyhedron
NUM_UNIQUE_FACES = 3
[docs] def __init__(self, a_param, b_param, c_param, alpha_param, beta_param, gamma_param, \
center=ORIGIN, ref_face_idx=1, ref_face_normal_along=Z_AXIS, \
ref_edge_idx=1, ref_edge_along=X_AXIS):
"""
Create an instance.
:type a_param: float
:param a_param: the length of the parallelepiped along edge a
:type b_param: float
:param b_param: the length of the parallelepiped along edge b
:type c_param: float
:param c_param: the length of the parallelepiped along edge c
:type alpha_param: float
:param alpha_param: the angle between edges b and c
:type beta_param: float
:param beta_param: the angle between edges a and c
:type gamma_param: float
:param gamma_param: the angle between edges a and b
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
params = [
a_param, b_param, c_param, alpha_param, beta_param, gamma_param
]
super(self.__class__, self).__init__(params, center, ref_face_idx,
ref_face_normal_along,
ref_edge_idx, ref_edge_along)
self.volume = self.getVolume(self.a_vec, self.b_vec, self.c_vec)
self.surface_area = self.getSurfaceArea(self.faces)
[docs] def getParallelepipedVertices(self, origin, a_vec, b_vec, c_vec):
"""
Get the vertices of the specified parallelepiped.
:type origin: numpy.array
:param origin: the point of origin of the lattic vectors
: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
:rtype: list of numpy.array
:return: the vertices of the parallelepiped
"""
return get_parallelepiped_vertices(origin, a_vec, b_vec, c_vec)
[docs] def getVolume(self, a_vec, b_vec, c_vec):
"""
Return the volume.
: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
:rtype: float
:return: volume in cubic Angstrom
"""
a_cross_b = numpy.cross(a_vec, b_vec)
volume = abs(numpy.dot(c_vec, a_cross_b))
return volume
[docs]class Slab(Parallelepiped, ConvexPolyhedron):
"""
Class to manage a slab.
"""
NAME = 'slab'
DEFAULT = [5.0, 10.0, 12.0]
DESCRIPTION = ['edge a length', 'edge b length', 'edge c length']
[docs] def __init__(self, a_param, b_param, c_param, \
center=ORIGIN, ref_face_idx=1, ref_face_normal_along=Z_AXIS, \
ref_edge_idx=1, ref_edge_along=X_AXIS):
"""
Create an instance.
:type a_param: float
:param a_param: the length of the parallelepiped along edge a
:type b_param: float
:param b_param: the length of the parallelepiped along edge b
:type c_param: float
:param c_param: the length of the parallelepiped along edge c
:type center: numpy.array
:param center: the center of the polyhedron
:type ref_face_idx: int
:param ref_face_idx: specifies which of this polyhedron's available
reference faces (symmetry unique faces sorted by decreasing area) to use
in the alignment
:type ref_face_normal_along: numpy.array triple
:param ref_face_normal_along: specifies the vector along which the
reference face normal will be aligned
:type ref_edge_idx: int
:param ref_edge_idx: specifies which of the reference faces available
reference edges (symmetry unique edges sorted by decreasing length) to use
in the alignment
:type ref_edge_along: numpy.array triple
:param ref_edge_along: specifies the vector along which the
reference edge will be aligned
"""
params = [a_param, b_param, c_param, 90.0, 90.0, 90.0]
ConvexPolyhedron.__init__(self, params, center, ref_face_idx,
ref_face_normal_along, ref_edge_idx,
ref_edge_along)
self.volume = self.getVolume(self.a_vec, self.b_vec, self.c_vec)
self.surface_area = self.getSurfaceArea(self.faces)
self.template.title = self.NAME
[docs]class Sphere(object):
"""
Class to manage a sphere.
"""
NAME = 'sphere'
DEFAULT = [5.0]
DESCRIPTION = ['radius']
TYPE = BASIC
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces,
# a sphere has no edges
NUM_UNIQUE_EDGES = 0
# this is the number of symmetry unique faces per polyhedron, a
# sphere has no faces
NUM_UNIQUE_FACES = 0
[docs] def __init__(self, radius, center=ORIGIN):
"""
Create an instance.
:type radius: float
:param radius: the radius of the sphere
:type center: numpy.array
:param center: the center of the sphere
"""
self.radius = radius
self.center = numpy.array(center)
self.volume = self.getVolume(self.radius)
self.surface_area = self.getSurfaceArea(self.radius)
[docs] def getVolume(self, radius):
"""
Return the volume.
:type radius: float
:param radius: radius in Angstrom
:rtype: float
:return: volume in cubic Angstrom
"""
factor = 4.0 * PI / 3.0
return factor * radius**3
[docs] def getSurfaceArea(self, radius):
"""
Return the surface area.
:type radius: float
:param radius: radius in Angstrom
:rtype: float
:return: surface area in square Angstrom
"""
factor = 4.0 * PI
return factor * radius**2
[docs] def pointInside(self, point):
"""
Return True if the query point is either on or inside of
this sphere.
:type point: numpy.array
:param point: the point in question
:rtype: bool
:return: True if the point in question is either on or inside
this sphere, False otherwise
"""
return transform.get_vector_magnitude(point -
self.center) <= self.radius
[docs]class Cylinder(object):
"""
Class to manage a cylinder.
"""
NAME = 'cylinder'
DEFAULT = [5.0, 25.0]
DESCRIPTION = ['radius', 'length']
TYPE = BASIC
# this is the number of symmetry unique edges per face, polyhedra
# containing faces with both differing numbers of vertices as well
# as different vertex-to-vertex distances are not supported, so this
# number is set for the polyhedron rather than the individual faces,
# a cylinder has no edges
NUM_UNIQUE_EDGES = 0
# this is the number of symmetry unique faces per polyhedron, a
# cylinder is considered to have a single symmetry unique face
NUM_UNIQUE_FACES = 1
[docs] def __init__(self, radius, length, center=ORIGIN):
"""
Create an instance.
:type radius: float
:param radius: the radius of the cylinder
:type length: float
:param length: the length of the cylinder
:type center: numpy.array
:param center: the center of the cylinder
"""
self.radius = radius
self.length = length
self.center = numpy.array(center)
self.volume = self.getVolume(self.radius, self.length)
self.surface_area = self.getSurfaceArea(self.radius, self.length)
[docs] def getVolume(self, radius, length):
"""
Return the volume.
:type radius: float
:param radius: radius in Angstrom
:type length: float
:param length: length in Angstrom
:rtype: float
:return: volume in cubic Angstrom
"""
return PI * length * radius**2
[docs] def getSurfaceArea(self, radius, length):
"""
Return the surface area.
:type radius: float
:param radius: radius in Angstrom
:type length: float
:param length: length in Angstrom
:rtype: float
:return: surface area in square Angstrom
"""
return 2.0 * PI * radius * (radius + length)
[docs] def pointInside(self, point):
"""
Return True if the query point is either on or inside of
this cylinder.
:type point: numpy.array
:param point: the point in question
:rtype: bool
:return: True if the point in question is either on or inside
this cylinder, False otherwise
"""
axis_unit_vec = X_AXIS
point_vec = point - self.center
length_proj = numpy.dot(point_vec, axis_unit_vec)
radius_proj = transform.get_vector_magnitude(point_vec - \
length_proj * axis_unit_vec)
return radius_proj <= self.radius and abs(length_proj) <= old_div(
self.length, 2.0)
SHAPES_NAMES_TO_OBJECTS_DICT = OrderedDict([
(Cube.NAME, Cube), (Tetrahedron.NAME, Tetrahedron),
(Octahedron.NAME, Octahedron), (Dodecahedron.NAME, Dodecahedron),
(Icosahedron.NAME, Icosahedron), (Cubeoctahedron.NAME, Cubeoctahedron),
(Slab.NAME, Slab), (Parallelepiped.NAME, Parallelepiped),
(Sphere.NAME, Sphere), (Cylinder.NAME, Cylinder)
])
[docs]def get_shape_object_by_name(name):
"""
Return a shape object by name.
:type name: str
:param name: the name of the object wanted
:rtype: object
:return: the shape object
"""
return SHAPES_NAMES_TO_OBJECTS_DICT.get(name)
[docs]def get_polygon_area(vertices):
"""
Return the area of the specified polygon using the shoelace
formula.
:type vertices: list
:param vertices: contains all vertices of the polygon,
each of which is a two dimensional numpy.array, i.e. x and y
:rtype: float
:return: the area of the polygon
"""
nvertices = len(vertices)
area = 0.0
for index in range(nvertices):
index_plus = (index + 1) % nvertices
index_minus = (index - 1) % nvertices
area += vertices[index][0] * (vertices[index_plus][1] - \
vertices[index_minus][1])
return old_div(abs(area), 2.0)
[docs]def get_reference_data(data, attr, num_unique, threshold):
"""
Return a list containing the num_unique number of unique data.
The data will be either a list of Face or a list of Edge characterized
using the attr AREA or LENGTH, respectively. These will be the
reference data used to orient the polyhedron.
:type data: list
:param data: either all Face objects for a given polyhedron or all Edge
objects for a given face
:type attr: str
:param attr: the attribute on which to characterize the data, either
AREA or LENGTH
:type num_unique: int
:param num_unique: the number of symmetry unique data
:type threshold: float
:param threshold: the threshold used to consider if two data are
equivalent by attr (either Ang. or Ang.^2)
:rtype: list
:return: unique data to serve as references
"""
# sort by decreasing attr
data_sorted = sorted(data, key=lambda x: getattr(x, attr), reverse=True)
# in order of decreasing attr collect num_unique unique data
# or run out of data before that in certain cases (see comment below)
primary_data = data_sorted.pop(0)
last_attr = getattr(primary_data, attr)
unique_data = [primary_data]
while len(unique_data) < num_unique and len(data_sorted):
next_data = data_sorted.pop(0)
next_attr = getattr(next_data, attr)
if abs(next_attr - last_attr) > threshold:
unique_data.append(next_data)
last_attr = next_attr
# in certain cases, for example slab and parallelepiped, it is possible
# that with the user input there are redundant faces and edges, i.e. less
# than num_unique unique faces or edges, in this case leave the num_unique
# alone and simply fill out unique_data with references to the last unique
# data, this just means that the user would still be given num_unique choices
# on which to align only that some of those would lead to redundant results
# which is harmless
diff = num_unique - len(unique_data)
if diff:
unique_data.extend([unique_data[-1]] * diff)
return unique_data
[docs]def get_parallelepiped_vertices(origin, a_vec, b_vec, c_vec, center=True):
"""
Get the vertices of the specified parallelepiped.
:type origin: numpy.array
:param origin: the point of origin of the lattic vectors
: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 center: bool
:param center: specifies whether or not to translate the
final vertices so that the centroid is at (0, 0, 0)
:rtype: list of numpy.array
:return: the vertices of the parallelepiped
"""
# by default the final vertices are for a parallelepiped centered
# at (0, 0, 0) and arbitrarily oriented as by default it
# will be rotated so that the normal of the largest face
# lies along the z-axis while the longest edge of that
# face lies along the x-axis
# get vertices
vertices = [
origin, origin + a_vec, origin + a_vec + b_vec, origin + b_vec,
origin + c_vec, origin + c_vec + a_vec, origin + c_vec + a_vec + b_vec,
origin + c_vec + b_vec
]
# we want the parallelepiped center to be at (0, 0, 0)
if center:
centroid = get_centroid(vertices)
vertices = [vertex - centroid for vertex in vertices]
return vertices
[docs]def get_centroid(vertices):
"""
Return the centroid of the provided vertices.
:type vertices: list
:param vertices: numpy.array array of points
:rtype: numpy.array
:return: the centroid
"""
return old_div(sum(vertices), len(vertices))