Coverage for wulffpack/core/form.py: 98%
82 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-05-06 13:03 +0000
1from typing import List, Dict, Union
2import numpy as np
3from .facet import Facet
4from .geometry import (is_array_in_arrays,
5 where_is_array_in_arrays)
8class Form():
9 """
10 A `Form` object contains all facets that are equivalent by the
11 symmetry of the system under consideration. For a cubic system,
12 this means that the form with Miller indices {100} will contain
13 the facets with normals [100], [-100], [010] etc.
15 Parameters
16 ----------
17 miller_indices
18 Miller indices for a representative member of the form.
19 energy
20 The surface energy per area for a facet in this form.
21 cell
22 The cell, with the basis vectors as columns.
23 symmetries
24 Symmetry elements of the system. If symmetry is broken
25 by, e.g., the presence of an interface, this list
26 should only contain the symmetries surviving.
27 parent_miller_indices
28 If symmetry is broken, it may still be of interest to
29 know what its form would be had not the symmetry been
30 broken. This attribute contains the Miller indices for
31 a representative of such a form.
32 """
34 def __init__(self, miller_indices: tuple,
35 energy: float,
36 cell: np.ndarray,
37 symmetries: List[np.ndarray],
38 parent_miller_indices: tuple):
39 self.miller_indices = miller_indices
40 self.energy = energy
41 self.parent_miller_indices = parent_miller_indices
43 # Create all facets belonging to this form
44 self.facets = []
45 used_normals = []
46 reciprocal_cell = np.linalg.inv(cell).T
47 normal = np.dot(reciprocal_cell, self.miller_indices)
48 normal_scaled = np.linalg.solve(cell, normal)
49 for R in symmetries:
50 new_normal_scaled = np.dot(R, normal_scaled)
51 new_normal = np.dot(cell, new_normal_scaled)
52 array_index = where_is_array_in_arrays(new_normal, used_normals)
53 if array_index == -1:
54 # Then we found a new facet
55 self.facets.append(Facet(normal=tuple(new_normal),
56 energy=energy,
57 symmetry=R))
58 used_normals.append(new_normal)
59 else:
60 # If the facet was already there, save instead the
61 # symmetry operation that took us there
62 facet = self.facets[array_index]
63 assert np.allclose(facet.normal, new_normal / np.linalg.norm(new_normal))
64 facet.symmetries.append(R)
66 @property
67 def area(self) -> float:
68 """Returns the total area of the facets in this form."""
69 return len(self.facets) * self.facets[0].area
71 @property
72 def surface_energy(self) -> float:
73 """Returns the total surface energy of the facets in this form."""
74 return len(self.facets) * self.facets[0].surface_energy
76 @property
77 def volume(self) -> float:
78 """
79 Returns the total volume formed by pyramids formed by each
80 facet and the origin.
81 """
82 return len(self.facets) * self.facets[0].get_volume()
84 @property
85 def edge_length(self) -> float:
86 """Returns the length of all facets in the form."""
87 return len(self.facets) * self.facets[0].perimeter_length
90def setup_forms(surface_energies: Dict,
91 cell: np.ndarray,
92 symmetries_restricted: List[np.ndarray],
93 symmetries_full: List[np.ndarray],
94 twin_boundary: Union[List, tuple] = None,
95 interface: Union[List, tuple] = None,
96 tol: float = 1e-6) -> Dict:
97 """
98 Create an adapted dictionary of surface energies based on
99 the symmetries specified.
101 This function is relevant for polycrystalline particles, the crystal
102 structure of which possess higher symmetry than the grains they are made
103 from. A decahedron, for example, is made of five FCC grains. FCC, being a
104 cubic crystal lattice, has 48 symmetry elements, but the grain is one of
105 five and has much fewer symmetry elements. This means that there are
106 multiple facets that would belong to the {111} form had not the symmetry
107 been broken. In the decahedral case, there are thus three inequivalent
108 families of facets that in the cubic case would all have been equivalent
109 facets belonging to the {111} form. These three families are, respectively,
110 the twin facets, the usually large facets "on top and underneath" the
111 particle as well as the re-entrance facets(the "notches"). All of these
112 will get their own key in the dictionary, describing a representative
113 member of the form.
115 Parameters
116 ----------
117 surface_energies
118 keys are either tuples with three integers (describing a form in the
119 cubic case) or the string `'twin'`/`'interface'`, values are
120 corresponding surface energies
121 cell
122 the basis vectors as columns
123 symmetries_restricted
124 the matrices for the symmetry elements in the broken symmetry
125 symmetry_full
126 the matrices for the symmetry elements in the case of full symmetry
127 (for example the 48 symmetry elements of the m-3m point group)
128 twin_boundary
129 Miller index for a twin boundary if there is one
130 interface
131 Miller index for an interface(to a substrate for example)
132 tol
133 Numerical tolerance when checking if energies are identical
135 Returns
136 -------
137 dictionary
138 keys are tuples describing the form(for each inequivalent form in the
139 broken symmetry case) with three integer, values are corresponding
140 surface energies
141 """
142 reciprocal_cell = np.linalg.inv(cell).T
143 forms = []
144 scaled_normals = [] # To be able to identify over-specified facets
145 ordered_facets = [] # --"--
146 if min(e for form, e in surface_energies.items() if form != 'interface') < 0:
147 raise ValueError('Please use only positive surface/twin energies')
149 for form, energy in surface_energies.items():
150 inequivalent_normals_scaled = []
151 if form == 'twin':
152 if twin_boundary:
153 forms.append(Form(miller_indices=twin_boundary,
154 energy=energy / 2,
155 cell=cell,
156 symmetries=symmetries_restricted,
157 parent_miller_indices=form))
158 elif form == 'interface':
159 if interface:
160 forms.append(Form(miller_indices=interface,
161 energy=energy,
162 cell=cell,
163 symmetries=symmetries_restricted,
164 parent_miller_indices=form))
165 else:
166 if len(form) == 4:
167 miller_indices = convert_bravais_miller_to_miller(form)
168 else:
169 miller_indices = form
170 normal = np.dot(reciprocal_cell, miller_indices) # Cartesian coords
171 normal_scaled = np.linalg.solve(cell, normal) # scaled coords
173 # Check that this facet is not symmetrically equivalent to another
174 scaled_normals.append(normal_scaled)
175 ordered_facets.append((form, energy))
176 for R_r in symmetries_full:
177 transformed_normal_scaled = np.dot(R_r, normal_scaled)
179 previous_index = where_is_array_in_arrays(transformed_normal_scaled,
180 scaled_normals[:-1])
181 if previous_index != -1:
182 # Then this facet is redundant.
183 # If it has the same energy, we just skip it,
184 # otherwise we need to throw an exception.
185 if abs(energy - ordered_facets[previous_index][1]) > tol:
186 raise ValueError(
187 f'The Miller indices {form} '
188 f'are equivalent to {ordered_facets[previous_index][0]} '
189 'by symmetry. If your system is such that they should be '
190 'inequivalent, please specify a primitive structure with '
191 'the appropriate point group symmetry or provide a '
192 'custom set of symmetries.')
193 else:
194 break
195 else:
196 # Then this facet is not symmetrically equivalent to another.
197 # Now we need to look at all symmetrically equivalent versions
198 # and check whether broken symmetry makes them inequivalent.
199 for R_r in symmetries_full:
200 transformed_normal_scaled = np.dot(R_r, normal_scaled)
201 for R_f in symmetries_restricted:
202 if is_array_in_arrays(np.dot(R_f, transformed_normal_scaled),
203 inequivalent_normals_scaled):
204 break
205 else:
206 transformed_normal = np.dot(cell, transformed_normal_scaled)
207 miller = np.linalg.solve(reciprocal_cell, transformed_normal)
208 rounded_miller = np.round(miller)
209 assert np.allclose(miller, rounded_miller)
210 rounded_miller = tuple(int(i) for i in rounded_miller)
211 forms.append(Form(miller_indices=rounded_miller,
212 energy=energy,
213 cell=cell,
214 symmetries=symmetries_restricted,
215 parent_miller_indices=form))
216 inequivalent_normals_scaled.append(transformed_normal_scaled)
217 return forms
220def convert_bravais_miller_to_miller(bravais_miller_indices: tuple) -> tuple:
221 """
222 Returns the Miller indices(three integer tuple)
223 corresponding to Bravais-Miller indices(for integer
224 tuple).
226 Paramters
227 ---------
228 bravais_miller_indices
229 Four integer tuple
230 """
231 if sum(bravais_miller_indices[:3]) != 0: 231 ↛ 232line 231 didn't jump to line 232 because the condition on line 231 was never true
232 raise ValueError('Invalid Bravais-Miller indices (h, k, i, l), '
233 'h + k + i != 0')
234 return tuple(bravais_miller_indices[i] for i in [0, 1, 3])