from typing import Dict, List, Tuple
import numpy as np
from ase import Atoms
from .facet import Facet
from .form import Form
from scipy.spatial import ConvexHull
from .geometry import is_array_in_arrays
[docs]
class BaseParticle():
"""
Base class for Wulff constructions, useful for systems of
(almost) any symmetry.
The Wulff construction calculates the vertices of the dual,
calculates the convex hull of those vertices using scipy's
ConvexHull, and then converts the dual back to the Wulff
shape.
Parameters
----------
forms
List of forms belonging to the system at hand
standardized_structure
The structural unit implicitly carrying information about
basis vectors etc.
natoms
Used as a target number of atoms if an atomistic structure
is requested. Together with `standardized_structure`,
`natoms` also defines the volume of the particle.
ngrains
If only a part of a complete particle is to be constructed,
for example a fifth of a decahedron, this should be reflected
by this parameter.
volume_scale
Decahedral and icosahedral particle need to be strained to
fill space. With this parameter, the volume is scaled such
that the volume of the particle is correct *after* this
straining.
tolerance
Numerical tolerance
"""
def __init__(self, forms: List[Form],
standardized_structure: Atoms,
natoms: int,
ngrains: int = 1,
volume_scale: int = 1.0,
tol: float = 1e-5):
self._natoms = natoms
self._standardized_structure = standardized_structure
self._forms = forms
self.tol = tol
# Calculate convex hull of dual vertices
duals = []
facets = [facet for facet in self._yield_facets()]
for facet in self._yield_facets():
facet_points = facet.normal * facet.energy
duals.append(facet_points / np.dot(facet_points, facet_points))
hull = ConvexHull(duals)
# Calculate vertices of dual again (i.e., back to Wulff shape)
for i, j, k in hull.simplices:
normal = np.cross(duals[j] - duals[i], duals[k] - duals[i])
if np.linalg.norm(normal) < 1e-12:
continue # Rare instances of duals on a line
normal *= facets[i].energy / np.dot(normal, facets[i].normal)
for facet_i in (i, j, k):
facets[facet_i].add_vertex(normal.copy())
# Remove vertices that are on the line between two other vertices
# (relevant in the Winterbottom case and some corner cases)
for facet in self._yield_facets():
facet.remove_redundant_vertices()
# Delete facets with less than three vertices
to_pop = []
for i, form in enumerate(self.forms):
if len(form.facets[0].vertices) < 3:
to_pop.append(i)
for i in reversed(to_pop):
self.forms.pop(i)
# Scale everyting such that the volume matches the specified one
target_volume = self.natoms * standardized_structure.get_volume() / \
len(standardized_structure)
scale_factor = (target_volume / (ngrains * self.volume * volume_scale)) ** (1 / 3)
self._scale_size(scale_factor)
def _scale_size(self, scale_factor: float):
"""
Scale the size of the particle.
Parameters
----------
scale_factor
Factor that all coordinates will be scaled with.
"""
for facet in self._yield_facets():
for i, _ in enumerate(facet.vertices):
facet.vertices[i] *= scale_factor
for i, _ in enumerate(facet.original_vertices):
facet.original_vertices[i] *= scale_factor
@property
def natoms(self) -> List[int]:
"""
The approximate number of atoms in the particle
(implicitly defining the volume).
"""
return self._natoms
@natoms.setter
def natoms(self, natoms: int):
"""
Change the approximate number of atoms and thus the volume.
"""
scale_factor = (natoms / self._natoms) ** (1 / 3.)
self._scale_size(scale_factor)
self._natoms = natoms
@property
def forms(self) -> List[Form]:
"""List of inequivalent forms for the particle"""
return self._forms
@property
def standardized_structure(self) -> Atoms:
"""
The standardized atomic structure that defines the geometry
and thus the meaning of the Miller indices. Also forms the building
blocks when `particle.atoms` is called.
"""
return self._standardized_structure
@property
def _twin_form(self) -> Form:
"""Returns the twin form if there is one, otherwise None."""
for form in self.forms:
if form.parent_miller_indices == 'twin':
return form
return None
[docs]
def get_continuous_color_scheme(self,
base_colors: dict = None,
normalize: bool = False) -> dict:
"""
Returns a dictionary with RGB colors for each form.
The colors smoothly interpolate between three
base colors, corresponding to (1, 1, 1), (1, 1, 0) and
(1, 0, 0). Note that this is sensible primarily for
cubic systems.
Parameters
----------
base_colors
User chosen colors for one or several of (1, 1, 1),
(1, 1, 0) and (1, 0, 0). To enforce, say, green
(1, 1, 1), use ``base_colors={(1, 1, 1): 'g'}``
normalize
If True, the norm of the RGB vectors will be 1. Note
that this may affect the ``base_colors`` too.
"""
from matplotlib.colors import to_rgba
# Adapt base_colors dictionary
new_base_colors = {}
if base_colors:
for key in base_colors.keys():
adapted_key = tuple(sorted(abs(i) for i in key))
new_base_colors[adapted_key] = base_colors[key]
# Make sure base_colors is complete
default_colors = {(1, 1, 1): '#2980B9',
(0, 1, 1): '#E92866',
(0, 0, 1): '#FFE82C'}
for key, color in default_colors.items():
if key not in new_base_colors:
new_base_colors[key] = color
A = [(1, 1, 1), (0, 1, 1), (0, 0, 1)]
base_color_matrix = np.array([to_rgba(new_base_colors[f])[:3] for f in A]).T
# Normalize rows in A
A = np.array([np.array(miller) / np.linalg.norm(miller) for miller in A]).T
# Calculate color for each form
colors = {}
for form in self.forms:
if type(form.parent_miller_indices) is not tuple:
continue
miller = np.array(sorted(abs(i) for i in form.parent_miller_indices))
color = np.linalg.solve(A, miller)
color /= np.linalg.norm(color)
color = np.dot(base_color_matrix, color)
if normalize:
color /= np.linalg.norm(color)
else:
color = [min(1, c) for c in color]
colors[form.parent_miller_indices] = color
return colors
[docs]
def make_plot(self, ax, alpha: float = 0.85, linewidth: float = 0.3, colors: dict = None):
"""
Plot a particle in an axis object. This function can be used to
make customized plots of particles.
Parameters
----------
ax : matplotlib Axes3DSubplot
An axis object with 3d projection
alpha
Opacity of the faces
linewidth
Thickness of lines between faces
colors
Allows custom colors for facets of all or a subset of forms,
example `{(1, 1, 1): '#FF0000'}`
Example
-------
In the following example, three different particles are
plotted in the same figure::
>>> from wulffpack import SingleCrystal, Decahedron, Icosahedron
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D
>>>
>>> surface_energies = {(1, 1, 1): 1.0,
... (1, 0, 0): 1.1,
... (1, 1, 0): 1.15,
... (3, 2, 1): 1.15}
>>> twin_energy = 0.05
>>>
>>> fig = plt.figure(figsize=(3*4.0, 4.0))
>>> ax = fig.add_subplot(131, projection='3d')
>>> particle = SingleCrystal(surface_energies)
>>> particle.make_plot(ax)
>>>
>>> ax = fig.add_subplot(132, projection='3d')
>>> particle = Decahedron(surface_energies,
... twin_energy=0.05)
>>> particle.make_plot(ax)
>>>
>>> ax = fig.add_subplot(133, projection='3d')
>>> particle = Icosahedron(surface_energies,
... twin_energy=0.05)
>>> particle.make_plot(ax)
>>>
>>> plt.subplots_adjust(top=1, bottom=0, left=0,
... right=1, wspace=0, hspace=0)
>>> plt.savefig('particles.png')
"""
from mpl_toolkits.mplot3d.art3d import (Poly3DCollection,
Line3DCollection)
# Standardize color dict
if colors is None:
colors = {}
default_colors = ['#D62728',
'#E377C2',
'#8C564B',
'#7F7F7F',
'#9467BD',
'#BCBD22',
'#17BECF',
'#AEC7E8',
'#FFBB78',
'#98DF8A',
'#FF9896',
'#C5B0D5',
'#C49C94',
'#F7B6D2',
'#C7C7C7',
'#DBDB8D',
'#9EDAE5']
mins = [1e9, 1e9, 1e9]
maxs = [-1e9, -1e9, -1e9]
poly3d = []
used_forms = []
color_counter = 0
for form in self.forms:
if form.parent_miller_indices == 'twin':
continue
# Determine color
if form.parent_miller_indices in colors:
color = colors[form.parent_miller_indices]
elif form.parent_miller_indices == (1, 1, 1):
color = '#4f81f1'
colors[form.parent_miller_indices] = color
elif tuple(sorted(form.parent_miller_indices)) == (0, 0, 1):
color = '#f8c73b'
colors[form.parent_miller_indices] = color
elif tuple(sorted(form.parent_miller_indices)) == (0, 1, 1):
color = '#2CA02C'
colors[form.parent_miller_indices] = color
else:
color = default_colors[color_counter % len(default_colors)]
colors[form.parent_miller_indices] = color
color_counter += 1
# Save the used forms to be able to make a legend
# (the colors dict may contain forms that are not present
# if it was supplied by the user)
if form.parent_miller_indices not in used_forms:
used_forms.append(form.parent_miller_indices)
# Plot all facets in the form
for facet in form.facets:
poly3d = [facet.ordered_vertices]
# Facets
collection = Poly3DCollection(poly3d, alpha=alpha)
collection.set_facecolor(color)
ax.add_collection3d(collection)
# Lines
ax.add_collection3d(Line3DCollection(poly3d, colors='k', linewidths=linewidth))
# Find proper ranges for the axes
for vertex in facet.ordered_vertices:
for i, v in enumerate(vertex):
if v > maxs[i]:
maxs[i] = v
if v < mins[i]:
mins[i] = v
# Plot fake lines to make a legend
for form in used_forms:
if isinstance(form, tuple):
label = ('{{' + '{}' * len(form) + '}}').format(*form)
else:
label = form
ax.plot([0, 0], [0, 0], color=colors[form], linewidth=8.0, label=label)
# Fiddle with the layout
extent = max(maxs) - min(mins)
bounds = [min(mins) + 0.15 * extent, max(maxs) - 0.15 * extent]
ax.set_xlim(*bounds)
ax.set_ylim(*bounds)
ax.set_zlim(*bounds)
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
ax.axis('off')
try:
ax.set_box_aspect((1, 1, 1))
except AttributeError:
# set_box_aspect only available in recent versions of matplotlib
pass
[docs]
def view(self, alpha: float = 0.85, linewidth: float = 0.3, colors: dict = None,
legend: bool = True, save_as: str = None):
"""
Use matplotlib to view a rendition of the particle.
Parameters
----------
alpha
Opacity of the faces
linewidth
Thickness of lines between faces
colors
Allows custom colors for facets of all or a subset of forms,
example `{(1, 1, 1): '#FF0000'}`
legend
Whether or not to show a legend with facet-color definitions
save_as
Filename to save figure as. If None, show the particle
with the GUI instead.
"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # NOQA
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111, projection='3d')
self.make_plot(ax, alpha=alpha, linewidth=linewidth, colors=colors)
plt.subplots_adjust(left=0, bottom=0, right=1, top=1)
if legend:
fig.legend(frameon=False)
plt.axis('off')
if save_as:
plt.savefig(save_as, transparent=True)
else:
plt.show()
[docs]
def write(self, filename: str):
"""
Write particle to file. The file format is derived from the
filename. Currently supported fileformats are:
* Wavefront .obj
Parameters
---------
filename
Filename of file to write to
"""
supported_file_formats = ['obj']
fileformat = filename.split('.')[-1]
if fileformat not in supported_file_formats:
raise ValueError('File format {} not supported, '.format(fileformat) +
'supported formats are: ' + ' '.join(supported_file_formats))
with open(filename, 'w') as f:
if fileformat == 'obj':
f.write('# Vertices\n')
for form in self.forms:
for facet in form.facets:
for vertex in facet.ordered_vertices[:-1]:
f.write('v {:15.8f} {:15.8f} {:15.8f}\n'.format(*vertex))
f.write('# Faces\n')
vertex_counter = 1
for form in self.forms:
# Define a group to be able to color each form separately
f.write('g {}\n'.format(form.parent_miller_indices))
for facet in form.facets:
# Tie the face to the right vertices. Count starts from 1.
s = ['{:5d}'.format(v + vertex_counter) for v in range(len(facet.vertices))]
vertex_counter += len(facet.vertices)
f.write('f ' + ' '.join(s) + '\n')
def _yield_facets(self, only_original_grains: bool = False) -> Facet:
"""Generate all facets that are included in the facets dictionary."""
for form in self.forms:
for facet in form.facets:
if only_original_grains and not facet.original_grain:
continue
yield facet
@property
def volume(self) -> float:
"""Returns the volume of the particle"""
volume = 0
for form in self.forms:
volume += form.volume
return volume
@property
def area(self) -> float:
"""
Returns total area of the surface of the particle (not including
twin boundaries).
"""
area = 0
for form in self.forms:
if form.parent_miller_indices != 'twin':
area += form.area
return area
@property
def edge_length(self) -> float:
"""Returns total edge length of the particle."""
edge_length = 0
for form in self.forms:
if form.parent_miller_indices != 'twin':
edge_length += form.edge_length / 2 # Every edge comes twice
return edge_length
@property
def number_of_corners(self) -> float:
"""Returns the number of corners (vertices) on the particle."""
unique_vertices = []
for form in self.forms:
if form.parent_miller_indices == 'twin':
continue
for facet in form.facets:
for vertex in facet.vertices:
if not is_array_in_arrays(vertex, unique_vertices):
unique_vertices.append(vertex)
return len(unique_vertices)
@property
def surface_energy(self) -> float:
"""
The total surface energy of the particle (including twin boundaries).
"""
E = 0
for form in self.forms:
E += form.surface_energy
return E
@property
def facet_fractions(self) -> Dict[tuple, float]:
"""
Returns a dict specifying fraction of each form
(not including twin boundaries).
"""
facet_fractions = {}
total_area = 0
for form in self.forms:
if form.parent_miller_indices == 'twin':
continue
area = form.area
total_area += area
facet_fractions[form.parent_miller_indices] = facet_fractions.get(
form.parent_miller_indices, 0) + area
for form in facet_fractions:
facet_fractions[form] = facet_fractions[form] / total_area
return facet_fractions
@property
def average_surface_energy(self) -> float:
"""
Average surface energy for the Wulff construction, i.e.,
a weighted average over all the facets, where the weights are
the area fraction of each facet.
"""
fractions = self.facet_fractions
weighted_average = 0
used_miller_indices = []
for form in self.forms:
if form.parent_miller_indices == 'twin' or \
form.parent_miller_indices in used_miller_indices:
continue
weighted_average += form.energy * fractions[form.parent_miller_indices]
used_miller_indices.append(form.parent_miller_indices)
return weighted_average
def _duplicate_particle(self, symmetries):
"""
Duplicate the particle by applying symmetry operations.
Useful for making an icosahedron from a single tetrahedron for example.
Parameters
----------
symmetries: list of NumPy arrays
Each array a symmetry operation(operating on Cartesian coordinates)
that should create a duplicate.
"""
for form in self.forms:
new_facets = []
for facet in form.facets:
for R in symmetries:
new_facet = Facet(normal=np.dot(R, facet.normal),
energy=facet.energy,
symmetry=None,
original_grain=False)
for vertex in facet.vertices:
new_facet.add_vertex(np.dot(R, vertex))
new_facets.append(new_facet)
for facet in new_facets:
form.facets.append(facet)
[docs]
def translate_particle(self, translation: np.ndarray):
"""
Translate the particle.
Parameters
----------
translation: list of 3 floats
Translation vector
"""
for facet in self._yield_facets():
for i, vertex in enumerate(facet.vertices):
facet.vertices[i] = vertex + translation
[docs]
def rotate_particle(self, rotation: np.ndarray):
"""
Rotate the particle.
Parameters
----------
rotation
Rotation matrix
"""
if abs(abs(np.linalg.det(rotation)) - 1.0) > self.tol:
raise ValueError('Provided matrix is not a rotation matrix '
'(its determinant does not equal 1 or -1).')
for facet in self._yield_facets():
for i, vertex in enumerate(facet.vertices):
facet.vertices[i] = np.dot(rotation, vertex)
facet.normal = np.dot(rotation, facet.normal)
facet.original_normal = np.dot(rotation, facet.original_normal)
self.standardized_structure.cell = np.dot(rotation, self.standardized_structure.cell.T).T
for atom in self.standardized_structure:
atom.position = np.dot(rotation, atom.position)
def _get_atoms(self, center_shift: Tuple[float, float, float] = None) -> Atoms:
"""
Returns an ASE Atoms object the atoms of which are all
within the facets. In polycrystalline particles, only
one grain will be returned.
Parameters
----------
center_shift
Shift of center in Cartesian coordinates.
"""
# Find max and min in x, y and z to bound the volume
# where there may be atoms
# (using scaled coordinates relative to the origin)
mins = [1e9, 1e9, 1e9]
maxs = [-1e9, -1e9, -1e9]
for facet in self._yield_facets(only_original_grains=True):
if facet.original_vertices:
vertices = facet.original_vertices
else:
vertices = facet.vertices
for vertex in vertices:
scaled_vertex = np.linalg.solve(self.standardized_structure.cell.T, vertex)
for i, v in enumerate(scaled_vertex):
if v < mins[i]:
mins[i] = v
if v > maxs[i]:
maxs[i] = v
# Prepare an unnecessarily large atoms object
mins = [int(i) - 1 for i in mins]
maxs = [int(np.ceil(i)) + 1 for i in maxs]
repeat = [int(np.ceil(i)) - int(j) for i, j in zip(maxs, mins)]
atoms = self.standardized_structure.repeat(repeat)
atoms.translate(np.dot(self.standardized_structure.cell.T, mins))
if center_shift:
atoms.translate(center_shift)
# Loop over all scaled coordinates and check whether
# they are within all facets
to_delete = []
for atom in atoms:
for facet in self._yield_facets(only_original_grains=True):
if np.dot(atom.position, facet.original_normal) > \
facet.distance_from_origin + 1e-5:
to_delete.append(atom.index)
break
del atoms[to_delete]
atoms.set_pbc(False)
atoms.set_cell([0., 0., 0.])
return atoms