Coverage for wulffpack/core/form.py: 98%
82 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-02 08:51 +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(surface_energies.values()) < 0:
147 raise ValueError('Please use only positive '
148 'surface/twin/interface energies')
150 for form, energy in surface_energies.items():
151 inequivalent_normals_scaled = []
152 if form == 'twin':
153 if twin_boundary:
154 forms.append(Form(miller_indices=twin_boundary,
155 energy=energy / 2,
156 cell=cell,
157 symmetries=symmetries_restricted,
158 parent_miller_indices=form))
159 elif form == 'interface':
160 if interface:
161 forms.append(Form(miller_indices=interface,
162 energy=energy,
163 cell=cell,
164 symmetries=symmetries_restricted,
165 parent_miller_indices=form))
166 else:
167 if len(form) == 4:
168 miller_indices = convert_bravais_miller_to_miller(form)
169 else:
170 miller_indices = form
171 normal = np.dot(reciprocal_cell, miller_indices) # Cartesian coords
172 normal_scaled = np.linalg.solve(cell, normal) # scaled coords
174 # Check that this facet is not symmetrically equivalent to another
175 scaled_normals.append(normal_scaled)
176 ordered_facets.append((form, energy))
177 for R_r in symmetries_full:
178 transformed_normal_scaled = np.dot(R_r, normal_scaled)
180 previous_index = where_is_array_in_arrays(transformed_normal_scaled,
181 scaled_normals[:-1])
182 if previous_index != -1:
183 # Then this facet is redundant.
184 # If it has the same energy, we just skip it,
185 # otherwise we need to throw an exception.
186 if abs(energy - ordered_facets[previous_index][1]) > tol:
187 raise ValueError(
188 f'The Miller indices {form} '
189 f'are equivalent to {ordered_facets[previous_index][0]} '
190 'by symmetry. If your system is such that they should be '
191 'inequivalent, please specify a primitive structure with '
192 'the appropriate point group symmetry or provide a '
193 'custom set of symmetries.')
194 else:
195 break
196 else:
197 # Then this facet is not symmetrically equivalent to another.
198 # Now we need to look at all symmetrically equivalent versions
199 # and check whether broken symmetry makes them inequivalent.
200 for R_r in symmetries_full:
201 transformed_normal_scaled = np.dot(R_r, normal_scaled)
202 for R_f in symmetries_restricted:
203 if is_array_in_arrays(np.dot(R_f, transformed_normal_scaled),
204 inequivalent_normals_scaled):
205 break
206 else:
207 transformed_normal = np.dot(cell, transformed_normal_scaled)
208 miller = np.linalg.solve(reciprocal_cell, transformed_normal)
209 rounded_miller = np.round(miller)
210 assert np.allclose(miller, rounded_miller)
211 rounded_miller = tuple(int(i) for i in rounded_miller)
212 forms.append(Form(miller_indices=rounded_miller,
213 energy=energy,
214 cell=cell,
215 symmetries=symmetries_restricted,
216 parent_miller_indices=form))
217 inequivalent_normals_scaled.append(transformed_normal_scaled)
218 return forms
221def convert_bravais_miller_to_miller(bravais_miller_indices: tuple) -> tuple:
222 """
223 Returns the Miller indices(three integer tuple)
224 corresponding to Bravais-Miller indices(for integer
225 tuple).
227 Paramters
228 ---------
229 bravais_miller_indices
230 Four integer tuple
231 """
232 if sum(bravais_miller_indices[:3]) != 0: 232 ↛ 233line 232 didn't jump to line 233, because the condition on line 232 was never true
233 raise ValueError('Invalid Bravais-Miller indices (h, k, i, l), '
234 'h + k + i != 0')
235 return tuple(bravais_miller_indices[i] for i in [0, 1, 3])