from typing import List, Dict, Union
import numpy as np
from .facet import Facet
from .geometry import (is_array_in_arrays,
where_is_array_in_arrays)
def setup_forms(surface_energies: Dict,
cell: np.ndarray,
symmetries_restricted: List[np.ndarray],
symmetries_full: List[np.ndarray],
twin_boundary: Union[List, tuple] = None,
interface: Union[List, tuple] = None,
tol: float = 1e-6) -> Dict:
"""
Create an adapted dictionary of surface energies based on
the symmetries specified.
This function is relevant for polycrystalline particles, the crystal
structure of which possess higher symmetry than the grains they are made
from. A decahedron, for example, is made of five FCC grains. FCC, being a
cubic crystal lattice, has 48 symmetry elements, but the grain is one of
five and has much fewer symmetry elements. This means that there are
multiple facets that would belong to the {111} form had not the symmetry
been broken. In the decahedral case, there are thus three inequivalent
families of facets that in the cubic case would all have been equivalent
facets belonging to the {111} form. These three families are, respectively,
the twin facets, the usually large facets "on top and underneath" the
particle as well as the re-entrance facets(the "notches"). All of these
will get their own key in the dictionary, describing a representative
member of the form.
Parameters
----------
surface_energies
keys are either tuples with three integers (describing a form in the
cubic case) or the string `'twin'`/`'interface'`, values are
corresponding surface energies
cell
the basis vectors as columns
symmetries_restricted
the matrices for the symmetry elements in the broken symmetry
symmetry_full
the matrices for the symmetry elements in the case of full symmetry
(for example the 48 symmetry elements of the m-3m point group)
twin_boundary
Miller index for a twin boundary if there is one
interface
Miller index for an interface(to a substrate for example)
tol
Numerical tolerance when checking if energies are identical
Returns
-------
dictionary
keys are tuples describing the form(for each inequivalent form in the
broken symmetry case) with three integer, values are corresponding
surface energies
"""
reciprocal_cell = np.linalg.inv(cell).T
forms = []
scaled_normals = [] # To be able to identify over-specified facets
ordered_facets = [] # --"--
if min(surface_energies.values()) < 0:
raise ValueError('Please use only positive '
'surface/twin/interface energies')
for form, energy in surface_energies.items():
inequivalent_normals_scaled = []
if form == 'twin':
if twin_boundary:
forms.append(Form(miller_indices=twin_boundary,
energy=energy / 2,
cell=cell,
symmetries=symmetries_restricted,
parent_miller_indices=form))
elif form == 'interface':
if interface:
forms.append(Form(miller_indices=interface,
energy=energy,
cell=cell,
symmetries=symmetries_restricted,
parent_miller_indices=form))
else:
if len(form) == 4:
miller_indices = convert_bravais_miller_to_miller(form)
else:
miller_indices = form
normal = np.dot(reciprocal_cell, miller_indices) # Cartesian coords
normal_scaled = np.linalg.solve(cell, normal) # scaled coords
# Check that this facet is not symmetrically equivalent to another
scaled_normals.append(normal_scaled)
ordered_facets.append((form, energy))
for R_r in symmetries_full:
transformed_normal_scaled = np.dot(R_r, normal_scaled)
previous_index = where_is_array_in_arrays(transformed_normal_scaled,
scaled_normals[:-1])
if previous_index != -1:
# Then this facet is redundant.
# If it has the same energy, we just skip it,
# otherwise we need to throw an exception.
if abs(energy - ordered_facets[previous_index][1]) > tol:
raise ValueError(
f'The Miller indices {form} '
f'are equivalent to {ordered_facets[previous_index][0]} '
'by symmetry. If your system is such that they should be '
'inequivalent, please specify a primitive structure with '
'the appropriate point group symmetry or provide a '
'custom set of symmetries.')
else:
break
else:
# Then this facet is not symmetrically equivalent to another.
# Now we need to look at all symmetrically equivalent versions
# and check whether broken symmetry makes them inequivalent.
for R_r in symmetries_full:
transformed_normal_scaled = np.dot(R_r, normal_scaled)
for R_f in symmetries_restricted:
if is_array_in_arrays(np.dot(R_f, transformed_normal_scaled),
inequivalent_normals_scaled):
break
else:
transformed_normal = np.dot(cell, transformed_normal_scaled)
miller = np.linalg.solve(reciprocal_cell, transformed_normal)
rounded_miller = np.round(miller)
assert np.allclose(miller, rounded_miller)
rounded_miller = tuple(int(i) for i in rounded_miller)
forms.append(Form(miller_indices=rounded_miller,
energy=energy,
cell=cell,
symmetries=symmetries_restricted,
parent_miller_indices=form))
inequivalent_normals_scaled.append(transformed_normal_scaled)
return forms
def convert_bravais_miller_to_miller(bravais_miller_indices: tuple) -> tuple:
"""
Returns the Miller indices(three integer tuple)
corresponding to Bravais-Miller indices(for integer
tuple).
Paramters
---------
bravais_miller_indices
Four integer tuple
"""
if sum(bravais_miller_indices[:3]) != 0:
raise ValueError('Invalid Bravais-Miller indices (h, k, i, l), '
'h + k + i != 0')
return tuple(bravais_miller_indices[i] for i in [0, 1, 3])