Coverage for wulffpack/decahedron.py: 100%
114 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-02 08:51 +0000
1from typing import Dict
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)
13class Decahedron(BaseParticle):
14 """
15 A `Decahedron` object is a generalized Wulff construction
16 of a decahedral particle.
18 Parameters
19 ----------
20 surface_energies
21 A dictionary with surface energies, where keys are
22 Miller indices and values surface energies (per area)
23 in a unit of choice, such as J/m^2.
24 twin_energy
25 Energy per area for twin boundaries
26 primitive_structure
27 Primitive cell to define the atomic structure used if an atomic
28 structure is requested. By default, an Au FCC structure is used.
29 The crystal has to have cubic symmetry.
30 natoms
31 Together with `lattice_parameter`, this parameter
32 defines the volume of the particle. If an atomic structure
33 is requested, the number of atoms will as closely as possible
34 match this value.
35 symprec
36 Numerical tolerance for symmetry analysis, forwarded to spglib.
37 tol
38 Numerical tolerance parameter.
40 Example
41 -------
42 The following example illustrates some possible uses of a
43 `Decahedron` object::
45 >>> from wulffpack import Decahedron
46 >>> from ase.build import bulk
47 >>> from ase.io import write
48 >>> surface_energies = {(1, 1, 1): 1.0, (1, 0, 0): 1.14}
49 >>> prim = bulk('Au')
50 >>> particle = Decahedron(surface_energies,
51 ... twin_energy=0.03,
52 ... primitive_structure=bulk('Au'))
53 >>> particle.view()
54 >>> write('decahedron.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)
66 symmetries = get_symmetries(standardized_structure, symprec)
67 if len(symmetries) < 48:
68 raise ValueError('A decahedron can only be created with a '
69 'primitive structure that has cubic symmetry')
71 twin_boundaries = [(-1, 1, 1), (1, -1, 1)]
72 symmetry_axes = [(0, 0, 1), np.cross(*twin_boundaries)]
73 inversion = [False, True]
74 broken_symmetries = break_symmetry(symmetries, symmetry_axes, inversion)
76 if twin_energy > 0.5 * min(surface_energies.values()):
77 raise ValueError('The construction expects a twin energy '
78 'that is smaller than 50 percent of the '
79 'smallest surface energy.')
80 surface_energies = surface_energies.copy()
81 surface_energies['twin'] = twin_energy
82 forms = setup_forms(surface_energies,
83 standardized_structure.cell.T,
84 broken_symmetries,
85 symmetries,
86 twin_boundary=twin_boundaries[0])
88 super().__init__(forms=forms,
89 standardized_structure=standardized_structure,
90 natoms=natoms,
91 ngrains=5,
92 volume_scale=_get_decahedral_scale_factor(),
93 tol=tol)
95 # Duplicate a single tetrahedron to form a complete decahedron
96 # with five tetrahedra
98 # Rotate such that fivefold axis is aligned with z axis
99 fivefold_vector = self.fivefold_axis_vector
100 rotation_axis = np.cross(fivefold_vector, [0, 0, 1])
101 angle = np.arccos(fivefold_vector[2] / np.linalg.norm(fivefold_vector))
102 R = get_rotation_matrix(angle, rotation_axis)
103 self.rotate_particle(R)
105 # Translate such that fivefold axis in x, y = 0, 0.
106 # Dangerous to use max and min in z because the particle
107 # may be truncated such that are several with the same value.
108 # Instead check which vertices are furthest away from the Ino face.
109 vertices = np.array(self._twin_form.facets[0].vertices)
110 # hard code because there may not be an Ino facet...
111 direction = np.array([1, 1, 0])
112 projections = [np.dot(v, direction) for v in vertices]
113 projections, ids = zip(*sorted(zip(projections,
114 list(range(len(projections))))))
115 min_vertex = vertices[ids[0]]
116 max_vertex = vertices[ids[1]]
118 # Check that rotation did its job
119 assert np.linalg.norm(max_vertex[:2] - min_vertex[:2]) < 1e-5
120 translation = np.array([-min_vertex[0], -min_vertex[1],
121 -(min_vertex[2] + max_vertex[2]) / 2])
122 self.translate_particle(translation)
124 # Rotate such that Ino facet aligns with a Cartesian vector
125 current = self._twin_form.facets[0].normal + self._twin_form.facets[1].normal
126 assert abs(current[2]) < 1e-5
127 target = (0, -1, 0)
128 angle = get_angle(current, target)
129 R = get_rotation_matrix(angle, [0, 0, 1])
130 self.rotate_particle(R)
132 # Back up the vertices in a separate list
133 # (that list is useful if creating an Atoms object)
134 for facet in self._yield_facets():
135 facet.original_vertices = [vertex.copy() for vertex in facet.vertices]
137 # Strain the particle to fill space
138 scale = _get_decahedral_scale_factor()
139 for facet in self._yield_facets():
140 for i in range(len(facet.vertices)):
141 facet.vertices[i][0] *= _get_decahedral_scale_factor()
142 # Scale normal too
143 facet.normal[1] *= scale
144 facet.normal[2] *= scale
145 facet.normal /= np.linalg.norm(facet.normal)
147 # Make five grains
148 rotations = []
149 for i in range(1, 5):
150 rotations.append(get_rotation_matrix(i * 2 * np.pi / 5, [0, 0, 1]))
151 self._duplicate_particle(rotations)
153 @property
154 def atoms(self) -> Atoms:
155 """
156 Returns an ASE Atoms object
157 """
158 atoms = self._get_atoms()
160 # We delete one atomic layer on a twin facet because we will
161 # get it back when we make five grains.
162 # But we need to save the atoms on the five-fold axis.
163 # Begin with identifying the maximum projection on one of the
164 # twin facets.
165 max_projection = -1e9
166 twin_direction = self._twin_form.facets[0].original_normal
167 for atom in atoms:
168 projection = np.dot(twin_direction, atom.position)
169 if projection > max_projection:
170 max_projection = projection
171 # Now delete all atoms that are
172 # close to the maximum projection
173 fivefold_atoms = Atoms()
174 twin_indices = []
175 for atom in atoms:
176 projection = np.dot(twin_direction, atom.position)
177 if abs(projection - max_projection) < 1e-5:
178 twin_indices.append(atom.index)
179 if abs(atom.position[0]) < 1e-5 and abs(atom.position[1]) < 1e-5:
180 fivefold_atoms.append(atom)
181 del atoms[twin_indices]
183 # Increase the distance from the x=0 plane by 1.8 % for each atom
184 for atom in atoms:
185 atom.position[0] = atom.position[0] * _get_decahedral_scale_factor()
187 # Make five grains
188 rotations = []
189 for i in range(1, 5):
190 rotations.append(get_rotation_matrix(i * 2 * np.pi / 5, [0, 0, 1]))
191 new_positions = []
192 for R in rotations:
193 for atom in atoms:
194 new_positions.append(np.dot(R, atom.position))
195 atoms += Atoms('{}{}'.format(self.standardized_structure[0].symbol,
196 len(new_positions)),
197 positions=new_positions)
198 atoms += fivefold_atoms
200 return atoms
202 @property
203 def fivefold_axis_vector(self) -> np.ndarray:
204 """
205 Returns a vector pointing in the
206 direction of the five-fold axis.
207 """
208 twins = self._twin_form.facets
209 v = np.cross(twins[0].normal,
210 twins[1].normal)
211 # Normalize
212 v /= np.linalg.norm(v)
213 return v
215 def get_strain_energy(self, shear_modulus: float, poissons_ratio: float) -> float:
216 """
217 Return strain energy as estimated with the formula provided in
218 A. Howie and L. D. Marks in Phil. Mag. A **49**, 95 (1984)
219 [HowMar84]_ (Eq. 10), which assumes an inhomogeneous strain in
220 the particle.
222 Warning
223 -------
224 This value is only approximate. If the decahedron is
225 heavily truncated, the returned strain energy may be highly
226 inaccurate.
228 Parameters
229 ----------
230 shear_modulus
231 Shear modulus of the material
232 poissons_ratio
233 Poisson's ratio of the material
234 """
235 eps_D = 0.0205
236 strain_energy_density = shear_modulus * eps_D ** 2 / 4
237 strain_energy_density /= (1 - poissons_ratio)
238 return strain_energy_density * self.volume
240 @property
241 def aspect_ratio(self) -> float:
242 """
243 Returns the aspect ratio of the decahedron, here defined as
244 the ratio between the longest distance between two vertices
245 projected on the fivefold axis versus the longest distance
246 between two vertices projected on the plane perpendicular
247 to the fivefold axis.
248 """
249 z_coords = []
250 xy_coords = []
251 for facet in self._yield_facets():
252 for vertex in facet.vertices:
253 z_coords.append(vertex[2])
254 xy_coords.append(vertex[:2])
256 # Longest distance between two points along z axis
257 # (because the particle has been rotated so that
258 # fivefold axis aligns with z)
259 height = max(z_coords) - min(z_coords)
261 # Find longest distance between two points in xy plane
262 width = 0
263 for i, xy_coord_i in enumerate(xy_coords):
264 for xy_coord_j in xy_coords[i + 1:]:
265 dist = ((xy_coord_i - xy_coord_j)**2).sum()
266 if dist > width:
267 width = dist
268 return height / width**0.5
271def _get_decahedral_scale_factor() -> float:
272 return np.sqrt(10 - 4 * np.sqrt(5))