Coverage for wulffpack/icosahedron.py: 98%
162 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-02 08:51 +0000
1from typing import Dict, Tuple, List
2import numpy as np
3from ase import Atoms
4from .core import BaseParticle
5from .core.form import setup_forms
6from .core.geometry import (get_standardized_structure,
7 get_symmetries,
8 get_angle,
9 get_rotation_matrix,
10 break_symmetry,
11 is_array_in_arrays)
14class Icosahedron(BaseParticle):
15 """
16 An `Icosahedron` object is a generalized Wulff construction
17 of an icosahedral particle.
19 Parameters
20 ----------
21 surface_energies
22 A dictionary with surface energies, where keys are
23 Miller indices and values surface energies (per area)
24 in a unit of choice, such as J/m^2.
25 twin_energy
26 Energy per area for twin boundaries
27 primitive_structure
28 Primitive cell to define the atomic structure used if an atomic
29 structure is requested. By default, an Au FCC structure is used.
30 The crystal has to have cubic symmetry.
31 natoms
32 Together with `lattice_parameter`, this parameter
33 defines the volume of the particle. If an atomic structure
34 is requested, the number of atoms will as closely as possible
35 match this value.
36 symprec
37 Numerical tolerance for symmetry analysis, forwarded to spglib.
38 tol
39 Numerical tolerance parameter.
41 Example
42 -------
43 The following example illustrates some possible uses of an
44 `Icosahedron` object::
46 >>> from wulffpack import Icosahedron
47 >>> from ase.build import bulk
48 >>> from ase.io import write
49 >>> surface_energies = {(1, 1, 1): 1.0, (1, 0, 0): 1.14}
50 >>> particle = Icosahedron(surface_energies,
51 ... twin_energy=0.03,
52 ... primitive_structure=bulk('Au'))
53 >>> particle.view()
54 >>> write('icosahedron.xyz', particle.atoms) # Writes atomic structure
56 """
58 def __init__(self,
59 surface_energies: Dict[tuple, float],
60 twin_energy: float,
61 primitive_structure: Atoms = None,
62 natoms: int = 1000,
63 symprec: float = 1e-5,
64 tol: float = 1e-5):
65 standardized_structure = get_standardized_structure(primitive_structure, symprec=symprec)
66 symmetries = get_symmetries(standardized_structure, symprec=symprec)
67 if len(symmetries) < 48:
68 raise ValueError('An icosahedron can only be created with a '
69 'primitive structure that has cubic symmetry')
71 broken_symmetries = break_symmetry(symmetries, [(1, 1, 1)])
73 if twin_energy > 0.5 * min(surface_energies.values()):
74 raise ValueError('The construction expects a twin energy '
75 'that is smaller than 50 percent of the '
76 'smallest surface energy.')
77 surface_energies = surface_energies.copy()
78 surface_energies['twin'] = twin_energy
80 forms = setup_forms(surface_energies,
81 standardized_structure.cell.T,
82 broken_symmetries,
83 symmetries,
84 twin_boundary=(-1, -1, 1))
86 super().__init__(forms=forms,
87 standardized_structure=standardized_structure,
88 natoms=natoms,
89 ngrains=20,
90 volume_scale=_get_icosahedral_scale_factor()**2,
91 tol=tol)
93 # Duplicate the single tetrahedron to form a complete icosahedron
94 # with 20 tetrahedra
96 # Translate such that tip of the tetrahedron is in the origin
97 min_proj = 1e12
98 for vertex in self._twin_form.facets[0].vertices:
99 proj = np.dot(vertex, (1, 1, 1))
100 if proj < min_proj:
101 min_proj = proj
102 translation = - vertex
103 self.translate_particle(translation)
105 # Back up the vertices in a separate list
106 # (that list is useful if creating an Atoms object)
107 for facet in self._yield_facets():
108 facet.original_vertices = [vertex.copy() for vertex in facet.vertices]
110 # Increase the distance from 111 to fill space
111 middle = np.array([1., 1., 1.])
112 middle /= np.linalg.norm(middle)
113 for facet in self._yield_facets():
114 for i, vertex in enumerate(facet.vertices):
115 if np.allclose(vertex, [0, 0, 0]):
116 continue
117 dist = vertex - np.dot(vertex, middle) * middle
118 facet.vertices[i] = vertex + (_get_icosahedral_scale_factor() - 1) * dist
119 # Tilt the normal
120 previous_normal = facet.normal
121 facet.normal = np.cross(facet.vertices[1] - facet.vertices[0],
122 facet.vertices[2] - facet.vertices[0])
123 if get_angle(facet.normal, previous_normal) > np.pi / 2:
124 facet.normal *= -1
125 facet.normal /= np.linalg.norm(facet.normal)
127 # Make 20 grains
128 symmetries = self._get_symmetry_operations()
129 self._duplicate_particle(symmetries)
131 @property
132 def atoms(self) -> Atoms:
133 """
134 Returns an ASE Atoms object
135 """
136 atoms = self._get_atoms()
138 # Handle fivefold axes and twin faces separately
139 # (they will be duplicated otherwise)
140 max_projections_twin = [-1e9, -1e9, -1e9]
141 twin_normals = [self._twin_form.facets[i].original_normal for i in range(3)]
143 # Find max projections onto twin directions
144 for atom in atoms:
145 for i, normal in enumerate(twin_normals):
146 projection = np.dot(normal, atom.position)
147 if projection > max_projections_twin[i]:
148 max_projections_twin[i] = projection
150 # Now identify all atoms that are
151 # close to the maximum projection
152 fivefold_atoms = Atoms()
153 twin_indices = [set(), set(), set()]
154 for atom in atoms:
155 for i, normal in enumerate(twin_normals):
156 projection = np.dot(normal, atom.position)
157 if abs(projection - max_projections_twin[i]) < 1e-5:
158 twin_indices[i].add(atom.index)
160 # Since they should not be duplicated we will remove them from
161 # the atoms object eventually
162 to_remove = list(twin_indices[0].union(*twin_indices[1:]))
164 # Identify the central atom
165 central_atom = list(twin_indices[0].intersection(*twin_indices[1:]))
166 if central_atom: 166 ↛ 173line 166 didn't jump to line 173, because the condition on line 166 was never false
167 assert len(central_atom) == 1
168 for twin in twin_indices:
169 twin.remove(central_atom[0])
170 central_atom = atoms[central_atom[0]]
172 # Identify the fivefold axes
173 fivefold_axes = []
174 fivefold_atoms = []
175 for i in range(3):
176 fivefold_axes.append(twin_indices[i].intersection(twin_indices[(i + 1) % 3]))
177 for fivefold_axis in fivefold_axes:
178 fivefold_atoms.append(atoms[list(fivefold_axis)])
179 for ind in fivefold_axis:
180 for i in range(3):
181 twin_indices[i].discard(ind)
182 twin_atoms = []
183 for twin in twin_indices:
184 twin_atoms.append(atoms[list(twin)])
185 del atoms[to_remove]
186 tetrahedron_atoms = atoms.copy()
188 # Increase the distance from 111 for each vertex to fill space
189 middle = np.array([1., 1., 1.])
190 middle /= np.linalg.norm(middle)
191 for atoms in [tetrahedron_atoms, *twin_atoms, *fivefold_atoms]:
192 for atom in atoms:
193 if np.allclose(atom.position, [0, 0, 0]): 193 ↛ 194line 193 didn't jump to line 194, because the condition on line 193 was never true
194 continue
195 dist = atom.position - np.dot(atom.position, middle) * middle
196 atom.position = atom.position + (_get_icosahedral_scale_factor() - 1) * dist
198 # Make 20 grains
199 symmetries = self._get_all_symmetry_operations()
200 new_positions = []
202 # Add the "bulk" of each tetrahedron
203 base_tetrahedron = np.array([1., 1., 1])
204 new_positions += _get_unique_coordinates(tetrahedron_atoms, symmetries, base_tetrahedron)
206 # Add the atoms on the 30 twin boundaries
207 base_twin = np.array(sum(atom.position for atom in twin_atoms[0]))
208 new_positions += _get_unique_coordinates(twin_atoms[0], symmetries, base_twin)
210 # Add the atoms along the fivefold axes
211 base_fivefold = np.array(sum(atom.position for atom in fivefold_atoms[0]))
212 new_positions += _get_unique_coordinates(fivefold_atoms[0], symmetries, base_fivefold)
214 atoms = Atoms('{}{}'.format(self.standardized_structure[0].symbol,
215 len(new_positions)),
216 positions=new_positions)
217 atoms.append(central_atom)
218 return atoms
220 def _get_two_fivefold_axes(self) -> Tuple[np.ndarray]:
221 """
222 Identify two fivefold axes as the two vertices which
223 are the furthest away from (1, 1, 1) (which is
224 in the middle of the face)
225 """
226 vertices = np.array(self._twin_form.facets[0].vertices)
227 for i, vertex in enumerate(vertices): 227 ↛ 231line 227 didn't jump to line 231, because the loop on line 227 didn't complete
228 if np.allclose(vertex, [0, 0, 0]): 228 ↛ 227line 228 didn't jump to line 227, because the condition on line 228 was never false
229 vertices = np.delete(vertices, i, 0)
230 break
231 direction = np.array([1, 1, 1])
232 angles = [get_angle(v, direction) for v in vertices]
233 angles, ids = zip(*sorted(zip(angles, list(range(len(angles))))))
234 return vertices[ids[-1]], vertices[ids[-2]]
236 def _get_all_symmetry_operations(self) -> List[np.ndarray]:
237 """
238 Get the 60 icosahedral symmetry operations in the coordinate
239 system defined by the particle as it is currently oriented.
240 """
241 fivefold_1, fivefold_2 = self._get_two_fivefold_axes()
242 R1 = get_rotation_matrix(1 * 2 * np.pi / 5, fivefold_1)
243 R2 = get_rotation_matrix(3 * 2 * np.pi / 5, fivefold_2)
244 symmetries = [np.eye(3)]
246 while len(symmetries) < 60:
247 for S in symmetries:
248 for R in [R1, R2]:
249 S_new = np.dot(R, S)
250 for S in symmetries:
251 if np.allclose(S_new, S):
252 break
253 else:
254 symmetries.append(S_new)
255 if len(symmetries) == 60:
256 break
257 assert len(symmetries) == 60
258 return symmetries
260 def _get_symmetry_operations(self) -> List[np.ndarray]:
261 """
262 Construct the subset of symmetry operations that duplicates
263 a single tetrahedron to all 20 tetrahedra.
264 """
265 fivefold_1, fivefold_2 = self._get_two_fivefold_axes()
266 symmetries = []
267 inversion = - np.eye(3)
268 symmetries.append(inversion)
269 down = get_rotation_matrix(2 * 2 * np.pi / 5, fivefold_1)
270 symmetries.append(down)
271 symmetries.append(np.dot(inversion, down))
272 for i in range(1, 5):
273 R = get_rotation_matrix(i * 2 * np.pi / 5, fivefold_2)
274 symmetries.append(R)
275 symmetries.append(np.dot(inversion, R))
276 symmetries.append(np.dot(R, down))
277 symmetries.append(np.dot(inversion, np.dot(R, down)))
278 assert len(symmetries) == 19
279 return symmetries
281 def get_strain_energy(self, shear_modulus, poissons_ratio):
282 """
283 Return a strain energy as estimated with the formula provided in
284 A. Howie and L. D. Marks in Phil. Mag. A **49**, 95 (1984)
285 [HowMar84]_ (Eq. 23), which assumes an inhomogeneous strain in
286 the particle.
288 Warning
289 -------
290 This value is only approximate. If the icosahedron is
291 heavily truncated, the returned strain energy may be highly
292 inaccurate.
294 Parameters
295 ----------
296 shear_modulus
297 Shear modulus of the material
298 poissons_ratio
299 Poisson's ratio of the material
300 """
301 eps_I = 0.0615
302 strain_energy_density = 2 * shear_modulus * eps_I ** 2 / 9
303 strain_energy_density *= (1 + poissons_ratio) / (1 - poissons_ratio)
304 return strain_energy_density * self.volume
307def _get_unique_coordinates(atoms: Atoms,
308 symmetries: List[np.ndarray],
309 base_element: np.ndarray) -> List[np.ndarray]:
310 """
311 Duplicate atoms with a list of symmetries, but avoid putting
312 atoms on top of each other.
314 atoms
315 Atoms object to duplicate
316 symmetries
317 List of symmetry elements to act on atoms with
318 base_element
319 An vector in Cartesian coordinates. If two symmetry
320 elements carries this vector to the same position,
321 one symmetry element will be skipped.
323 Returns
324 -------
325 A list of Cartesian coordinates with new atomic positions,
326 including the original ones.
327 """
328 unique_coordinates = []
329 symmetrical_elements = []
330 for R in symmetries:
331 element = np.dot(R, base_element)
332 if is_array_in_arrays(element, symmetrical_elements):
333 continue
334 else:
335 symmetrical_elements.append(element)
336 for atom in atoms:
337 unique_coordinates.append(np.dot(R, atom.position))
338 return unique_coordinates
341def _get_icosahedral_scale_factor():
342 k = (5 + np.sqrt(5)) / 8
343 return np.sqrt(2 / (3 * k - 1))