Coverage for wulffpack/core/form.py: 99%
101 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-29 11:40 +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 def __str__(self) -> str:
67 width = 38
68 s = []
69 s += ['{s:-^{n}}'.format(s=' Form ', n=width)]
70 flds = dict(miller_indices='', parent_miller_indices='',
71 area='12.4f', surface_energy='12.2f', edge_length='12.4f')
72 for fld, fmt in flds.items():
73 s += ['{fld:22} : {value:{fmt}}'.format(fld=fld, value=getattr(self, fld), fmt=fmt)]
74 return '\n'.join(s)
76 def _repr_html_(self) -> str:
77 s = []
78 s += ['<table border="1" class="dataframe"']
79 s += [
80 '<thead><tr><th style="text-align: left;">Property</th><th>Value</th></tr></thead>'
81 ]
82 s += ['<tbody>']
83 flds = dict(miller_indices='', parent_miller_indices='',
84 area='12.4f', surface_energy='12.2f', edge_length='12.4f')
85 for fld, fmt in flds.items():
86 s += [
87 f'<tr><td style="text-align: left">{fld:22}</td>'
88 '<td>{value:{fmt}}</td><tr>'.format(value=getattr(self, fld), fmt=fmt)
89 ]
90 s += ['</tbody>']
91 s += ['</table>']
92 return ''.join(s)
94 @property
95 def area(self) -> float:
96 """Returns the total area of the facets in this form."""
97 return len(self.facets) * self.facets[0].area
99 @property
100 def surface_energy(self) -> float:
101 """Returns the total surface energy of the facets in this form."""
102 return len(self.facets) * self.facets[0].surface_energy
104 @property
105 def volume(self) -> float:
106 """
107 Returns the total volume formed by pyramids formed by each
108 facet and the origin.
109 """
110 return len(self.facets) * self.facets[0].get_volume()
112 @property
113 def edge_length(self) -> float:
114 """Returns the length of all facets in the form."""
115 return len(self.facets) * self.facets[0].perimeter_length
118def setup_forms(surface_energies: Dict,
119 cell: np.ndarray,
120 symmetries_restricted: List[np.ndarray],
121 symmetries_full: List[np.ndarray],
122 twin_boundary: Union[List, tuple] = None,
123 interface: Union[List, tuple] = None,
124 tol: float = 1e-6) -> Dict:
125 """
126 Create an adapted dictionary of surface energies based on
127 the symmetries specified.
129 This function is relevant for polycrystalline particles, the crystal
130 structure of which possess higher symmetry than the grains they are made
131 from. A decahedron, for example, is made of five FCC grains. FCC, being a
132 cubic crystal lattice, has 48 symmetry elements, but the grain is one of
133 five and has much fewer symmetry elements. This means that there are
134 multiple facets that would belong to the {111} form had not the symmetry
135 been broken. In the decahedral case, there are thus three inequivalent
136 families of facets that in the cubic case would all have been equivalent
137 facets belonging to the {111} form. These three families are, respectively,
138 the twin facets, the usually large facets "on top and underneath" the
139 particle as well as the re-entrance facets(the "notches"). All of these
140 will get their own key in the dictionary, describing a representative
141 member of the form.
143 Parameters
144 ----------
145 surface_energies
146 keys are either tuples with three integers (describing a form in the
147 cubic case) or the string `'twin'`/`'interface'`, values are
148 corresponding surface energies
149 cell
150 the basis vectors as columns
151 symmetries_restricted
152 the matrices for the symmetry elements in the broken symmetry
153 symmetry_full
154 the matrices for the symmetry elements in the case of full symmetry
155 (for example the 48 symmetry elements of the m-3m point group)
156 twin_boundary
157 Miller index for a twin boundary if there is one
158 interface
159 Miller index for an interface(to a substrate for example)
160 tol
161 Numerical tolerance when checking if energies are identical
163 Returns
164 -------
165 dictionary
166 keys are tuples describing the form(for each inequivalent form in the
167 broken symmetry case) with three integer, values are corresponding
168 surface energies
169 """
170 reciprocal_cell = np.linalg.inv(cell).T
171 forms = []
172 scaled_normals = [] # To be able to identify over-specified facets
173 ordered_facets = [] # --"--
174 if min(e for form, e in surface_energies.items() if form != 'interface') < 0:
175 raise ValueError('Please use only positive surface/twin energies')
177 for form, energy in surface_energies.items():
178 inequivalent_normals_scaled = []
179 if form == 'twin':
180 if twin_boundary:
181 forms.append(Form(miller_indices=twin_boundary,
182 energy=energy / 2,
183 cell=cell,
184 symmetries=symmetries_restricted,
185 parent_miller_indices=form))
186 elif form == 'interface':
187 if interface:
188 forms.append(Form(miller_indices=interface,
189 energy=energy,
190 cell=cell,
191 symmetries=symmetries_restricted,
192 parent_miller_indices=form))
193 else:
194 if len(form) == 4:
195 miller_indices = convert_bravais_miller_to_miller(form)
196 else:
197 miller_indices = form
198 normal = np.dot(reciprocal_cell, miller_indices) # Cartesian coords
199 normal_scaled = np.linalg.solve(cell, normal) # scaled coords
201 # Check that this facet is not symmetrically equivalent to another
202 scaled_normals.append(normal_scaled)
203 ordered_facets.append((form, energy))
204 for R_r in symmetries_full:
205 transformed_normal_scaled = np.dot(R_r, normal_scaled)
207 previous_index = where_is_array_in_arrays(transformed_normal_scaled,
208 scaled_normals[:-1])
209 if previous_index != -1:
210 # Then this facet is redundant.
211 # If it has the same energy, we just skip it,
212 # otherwise we need to throw an exception.
213 if abs(energy - ordered_facets[previous_index][1]) > tol:
214 raise ValueError(
215 f'The Miller indices {form} '
216 f'are equivalent to {ordered_facets[previous_index][0]} '
217 'by symmetry. If your system is such that they should be '
218 'inequivalent, please specify a primitive structure with '
219 'the appropriate point group symmetry or provide a '
220 'custom set of symmetries.')
221 else:
222 break
223 else:
224 # Then this facet is not symmetrically equivalent to another.
225 # Now we need to look at all symmetrically equivalent versions
226 # and check whether broken symmetry makes them inequivalent.
227 for R_r in symmetries_full:
228 transformed_normal_scaled = np.dot(R_r, normal_scaled)
229 for R_f in symmetries_restricted:
230 if is_array_in_arrays(np.dot(R_f, transformed_normal_scaled),
231 inequivalent_normals_scaled):
232 break
233 else:
234 transformed_normal = np.dot(cell, transformed_normal_scaled)
235 miller = np.linalg.solve(reciprocal_cell, transformed_normal)
236 rounded_miller = np.round(miller)
237 assert np.allclose(miller, rounded_miller)
238 rounded_miller = tuple(int(i) for i in rounded_miller)
239 forms.append(Form(miller_indices=rounded_miller,
240 energy=energy,
241 cell=cell,
242 symmetries=symmetries_restricted,
243 parent_miller_indices=form))
244 inequivalent_normals_scaled.append(transformed_normal_scaled)
245 return forms
248def convert_bravais_miller_to_miller(bravais_miller_indices: tuple) -> tuple:
249 """
250 Returns the Miller indices(three integer tuple)
251 corresponding to Bravais-Miller indices(for integer
252 tuple).
254 Paramters
255 ---------
256 bravais_miller_indices
257 Four integer tuple
258 """
259 if sum(bravais_miller_indices[:3]) != 0: 259 ↛ 260line 259 didn't jump to line 260 because the condition on line 259 was never true
260 raise ValueError('Invalid Bravais-Miller indices (h, k, i, l), '
261 'h + k + i != 0')
262 return tuple(bravais_miller_indices[i] for i in [0, 1, 3])