Coverage for wulffpack/core/base_particle.py: 97%
330 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-05-06 13:03 +0000
1from typing import Dict, List, Tuple
2import numpy as np
3from ase import Atoms
4from .facet import Facet
5from .form import Form
6from scipy.spatial import ConvexHull
7from .geometry import is_array_in_arrays
10class BaseParticle():
11 """
12 Base class for Wulff constructions, useful for systems of
13 (almost) any symmetry.
15 Parameters
16 ----------
17 forms
18 List of forms belonging to the system at hand
19 standardized_structure
20 The structural unit implicitly carrying information about
21 basis vectors etc.
22 natoms
23 Used as a target number of atoms if an atomistic structure
24 is requested. Together with `standardized_structure`,
25 `natoms` also defines the volume of the particle.
26 ngrains
27 If only a part of a complete particle is to be constructed,
28 for example a fifth of a decahedron, this should be reflected
29 by this parameter.
30 volume_scale
31 Decahedral and icosahedral particle need to be strained to
32 fill space. With this parameter, the volume is scaled such
33 that the volume of the particle is correct *after* this
34 straining.
35 tolerance
36 Numerical tolerance
37 """
39 def __init__(self, forms: List[Form],
40 standardized_structure: Atoms,
41 natoms: int,
42 ngrains: int = 1,
43 volume_scale: int = 1.0,
44 tol: float = 1e-5):
45 self._natoms = natoms
46 self._standardized_structure = standardized_structure
47 self._forms = forms
48 self.tol = tol
49 self._SMALL = 1e-12
51 facets = list(self._yield_facets())
52 if min(f.energy for f in facets) < 0: # interface energy can be negative
53 self._closest_vertices(facets)
54 else:
55 self._dual_wulff_shape(facets)
57 # Remove vertices that are on the line between two other vertices
58 # (relevant in the Winterbottom case and some corner cases)
59 for facet in self._yield_facets():
60 facet.remove_redundant_vertices()
62 # Delete facets with less than three vertices
63 to_pop = []
64 for i, form in enumerate(self.forms):
65 if len(form.facets[0].vertices) < 3:
66 to_pop.append(i)
67 for i in reversed(to_pop):
68 self.forms.pop(i)
70 # Scale everyting such that the volume matches the specified one
71 target_volume = self.natoms * standardized_structure.get_volume() / \
72 len(standardized_structure)
73 scale_factor = (target_volume / (ngrains * self.volume * volume_scale)) ** (1 / 3)
74 self._scale_size(scale_factor)
76 def _closest_vertices(self, facets):
77 """
78 The closest vertex technique calculates all possible vertices and eliminates
79 vertices outside of the Wulff shape.
80 This works for negative interface energies.
81 """
82 # Calculate intersection point for all possible combinations of three facets
83 facets.sort(key=lambda x: x.energy) # compare with small energies first in inner loop
84 for i, fi in enumerate(facets):
85 ni = fi.normal
86 for j, fj in enumerate(facets[:i]):
87 nj = fj.normal
88 if abs(abs(np.dot(ni, nj)) - 1) < self._SMALL:
89 continue # skip parallel facets
90 for fk in facets[:j]:
91 mat = np.vstack([ni, nj, fk.normal])
92 if np.linalg.det(mat) == 0:
93 continue # skip singular matrices
94 vertex = np.dot(np.linalg.inv(mat), [fi.energy, fj.energy, fk.energy])
96 # check if vertex lies outside of Wulff shape
97 for fl in facets:
98 projection = np.dot(vertex, fl.normal)
99 if projection > fl.energy + self._SMALL:
100 break
101 else: # i.e. loop does not break
102 for f in [fi, fj, fk]:
103 f.add_vertex(vertex.copy())
105 def _dual_wulff_shape(self, facets):
106 """
107 The Wulff construction calculates the vertices of the dual,
108 calculates the convex hull of those vertices using scipy's
109 ConvexHull, and then converts the dual back to the Wulff
110 shape.
111 """
112 # Calculate convex hull of dual vertices
113 duals = []
114 for facet in self._yield_facets():
115 facet_points = facet.normal * facet.energy
116 duals.append(facet_points / np.dot(facet_points, facet_points))
117 hull = ConvexHull(duals)
119 # Calculate vertices of dual again (i.e., back to Wulff shape)
120 for i, j, k in hull.simplices:
121 normal = np.cross(duals[j] - duals[i], duals[k] - duals[i])
122 if np.linalg.norm(normal) < self._SMALL: 122 ↛ 123line 122 didn't jump to line 123 because the condition on line 122 was never true
123 continue # Rare instances of duals on a line
124 normal *= facets[i].energy / np.dot(normal, facets[i].normal)
125 for facet_i in (i, j, k):
126 facets[facet_i].add_vertex(normal.copy())
128 def _scale_size(self, scale_factor: float):
129 """
130 Scale the size of the particle.
132 Parameters
133 ----------
134 scale_factor
135 Factor that all coordinates will be scaled with.
136 """
137 for facet in self._yield_facets():
138 for i, _ in enumerate(facet.vertices):
139 facet.vertices[i] *= scale_factor
140 for i, _ in enumerate(facet.original_vertices):
141 facet.original_vertices[i] *= scale_factor
143 @property
144 def natoms(self) -> List[int]:
145 """
146 The approximate number of atoms in the particle
147 (implicitly defining the volume).
148 """
149 return self._natoms
151 @natoms.setter
152 def natoms(self, natoms: int):
153 """
154 Change the approximate number of atoms and thus the volume.
155 """
156 scale_factor = (natoms / self._natoms) ** (1 / 3.)
157 self._scale_size(scale_factor)
158 self._natoms = natoms
160 @property
161 def forms(self) -> List[Form]:
162 """List of inequivalent forms for the particle"""
163 return self._forms
165 @property
166 def standardized_structure(self) -> Atoms:
167 """
168 The standardized atomic structure that defines the geometry
169 and thus the meaning of the Miller indices. Also forms the building
170 blocks when `particle.atoms` is called.
171 """
172 return self._standardized_structure
174 @property
175 def _twin_form(self) -> Form:
176 """Returns the twin form if there is one, otherwise None."""
177 for form in self.forms:
178 if form.parent_miller_indices == 'twin':
179 return form
180 return None
182 def get_continuous_color_scheme(self,
183 base_colors: dict = None,
184 normalize: bool = False) -> dict:
185 """
186 Returns a dictionary with RGB colors for each form.
187 The colors smoothly interpolate between three
188 base colors, corresponding to (1, 1, 1), (1, 1, 0) and
189 (1, 0, 0). Note that this is sensible primarily for
190 cubic systems.
192 Parameters
193 ----------
194 base_colors
195 User chosen colors for one or several of (1, 1, 1),
196 (1, 1, 0) and (1, 0, 0). To enforce, say, green
197 (1, 1, 1), use ``base_colors={(1, 1, 1): 'g'}``
198 normalize
199 If True, the norm of the RGB vectors will be 1. Note
200 that this may affect the ``base_colors`` too.
201 """
202 from matplotlib.colors import to_rgba
204 # Adapt base_colors dictionary
205 new_base_colors = {}
206 if base_colors:
207 for key in base_colors.keys():
208 adapted_key = tuple(sorted(abs(i) for i in key))
209 new_base_colors[adapted_key] = base_colors[key]
211 # Make sure base_colors is complete
212 default_colors = {(1, 1, 1): '#2980B9',
213 (0, 1, 1): '#E92866',
214 (0, 0, 1): '#FFE82C'}
215 for key, color in default_colors.items():
216 if key not in new_base_colors:
217 new_base_colors[key] = color
218 A = [(1, 1, 1), (0, 1, 1), (0, 0, 1)]
219 base_color_matrix = np.array([to_rgba(new_base_colors[f])[:3] for f in A]).T
221 # Normalize rows in A
222 A = np.array([np.array(miller) / np.linalg.norm(miller) for miller in A]).T
224 # Calculate color for each form
225 colors = {}
226 for form in self.forms:
227 if type(form.parent_miller_indices) is not tuple:
228 continue
229 miller = np.array(sorted(abs(i) for i in form.parent_miller_indices))
230 color = np.linalg.solve(A, miller)
231 color /= np.linalg.norm(color)
232 color = np.dot(base_color_matrix, color)
233 if normalize:
234 color /= np.linalg.norm(color)
235 else:
236 color = [min(1, c) for c in color]
237 colors[form.parent_miller_indices] = color
238 return colors
240 def make_plot(self, ax, alpha: float = 0.85, linewidth: float = 0.3, colors: dict = None):
241 """
242 Plot a particle in an axis object. This function can be used to
243 make customized plots of particles.
245 Parameters
246 ----------
247 ax : matplotlib Axes3DSubplot
248 An axis object with 3d projection
249 alpha
250 Opacity of the faces
251 linewidth
252 Thickness of lines between faces
253 colors
254 Allows custom colors for facets of all or a subset of forms,
255 example `{(1, 1, 1): '#FF0000'}`
257 Example
258 -------
259 In the following example, three different particles are
260 plotted in the same figure::
262 >>> from wulffpack import SingleCrystal, Decahedron, Icosahedron
263 >>> import matplotlib.pyplot as plt
264 >>> from mpl_toolkits.mplot3d import Axes3D
265 >>>
266 >>> surface_energies = {(1, 1, 1): 1.0,
267 ... (1, 0, 0): 1.1,
268 ... (1, 1, 0): 1.15,
269 ... (3, 2, 1): 1.15}
270 >>> twin_energy = 0.05
271 >>>
272 >>> fig = plt.figure(figsize=(3*4.0, 4.0))
273 >>> ax = fig.add_subplot(131, projection='3d')
274 >>> particle = SingleCrystal(surface_energies)
275 >>> particle.make_plot(ax)
276 >>>
277 >>> ax = fig.add_subplot(132, projection='3d')
278 >>> particle = Decahedron(surface_energies,
279 ... twin_energy=0.05)
280 >>> particle.make_plot(ax)
281 >>>
282 >>> ax = fig.add_subplot(133, projection='3d')
283 >>> particle = Icosahedron(surface_energies,
284 ... twin_energy=0.05)
285 >>> particle.make_plot(ax)
286 >>>
287 >>> plt.subplots_adjust(top=1, bottom=0, left=0,
288 ... right=1, wspace=0, hspace=0)
289 >>> plt.savefig('particles.png')
291 """
292 from mpl_toolkits.mplot3d.art3d import (Poly3DCollection,
293 Line3DCollection)
295 # Standardize color dict
296 if colors is None:
297 colors = {}
299 default_colors = ['#D62728',
300 '#E377C2',
301 '#8C564B',
302 '#7F7F7F',
303 '#9467BD',
304 '#BCBD22',
305 '#17BECF',
306 '#AEC7E8',
307 '#FFBB78',
308 '#98DF8A',
309 '#FF9896',
310 '#C5B0D5',
311 '#C49C94',
312 '#F7B6D2',
313 '#C7C7C7',
314 '#DBDB8D',
315 '#9EDAE5']
317 mins = [1e9, 1e9, 1e9]
318 maxs = [-1e9, -1e9, -1e9]
319 poly3d = []
321 used_forms = []
322 color_counter = 0
323 for form in self.forms:
324 if form.parent_miller_indices == 'twin':
325 continue
327 # Determine color
328 if form.parent_miller_indices in colors:
329 color = colors[form.parent_miller_indices]
330 elif form.parent_miller_indices == (1, 1, 1):
331 color = '#4f81f1'
332 colors[form.parent_miller_indices] = color
333 elif tuple(sorted(form.parent_miller_indices)) == (0, 0, 1): 333 ↛ 336line 333 didn't jump to line 336 because the condition on line 333 was always true
334 color = '#f8c73b'
335 colors[form.parent_miller_indices] = color
336 elif tuple(sorted(form.parent_miller_indices)) == (0, 1, 1):
337 color = '#2CA02C'
338 colors[form.parent_miller_indices] = color
339 else:
340 color = default_colors[color_counter % len(default_colors)]
341 colors[form.parent_miller_indices] = color
342 color_counter += 1
344 # Save the used forms to be able to make a legend
345 # (the colors dict may contain forms that are not present
346 # if it was supplied by the user)
347 if form.parent_miller_indices not in used_forms: 347 ↛ 351line 347 didn't jump to line 351 because the condition on line 347 was always true
348 used_forms.append(form.parent_miller_indices)
350 # Plot all facets in the form
351 for facet in form.facets:
352 poly3d = [facet.ordered_vertices]
354 # Facets
355 collection = Poly3DCollection(poly3d, alpha=alpha)
356 collection.set_facecolor(color)
357 ax.add_collection3d(collection)
359 # Lines
360 ax.add_collection3d(Line3DCollection(poly3d, colors='k', linewidths=linewidth))
362 # Find proper ranges for the axes
363 for vertex in facet.ordered_vertices:
364 for i, v in enumerate(vertex):
365 if v > maxs[i]:
366 maxs[i] = v
367 if v < mins[i]:
368 mins[i] = v
370 # Plot fake lines to make a legend
371 for form in used_forms:
372 if isinstance(form, tuple):
373 label = ('{{' + '{}' * len(form) + '}}').format(*form)
374 else:
375 label = form
376 ax.plot([0, 0], [0, 0], color=colors[form], linewidth=8.0, label=label)
378 # Fiddle with the layout
379 extent = max(maxs) - min(mins)
380 bounds = [min(mins) + 0.15 * extent, max(maxs) - 0.15 * extent]
381 ax.set_xlim(*bounds)
382 ax.set_ylim(*bounds)
383 ax.set_zlim(*bounds)
384 ax.set_xticks([])
385 ax.set_yticks([])
386 ax.set_zticks([])
387 ax.axis('off')
388 try:
389 ax.set_box_aspect((1, 1, 1))
390 except AttributeError:
391 # set_box_aspect only available in recent versions of matplotlib
392 pass
394 def view(self, alpha: float = 0.85, linewidth: float = 0.3, colors: dict = None,
395 legend: bool = True, save_as: str = None):
396 """
397 Use matplotlib to view a rendition of the particle.
399 Parameters
400 ----------
401 alpha
402 Opacity of the faces
403 linewidth
404 Thickness of lines between faces
405 colors
406 Allows custom colors for facets of all or a subset of forms,
407 example `{(1, 1, 1): '#FF0000'}`
408 legend
409 Whether or not to show a legend with facet-color definitions
410 save_as
411 Filename to save figure as. If None, show the particle
412 with the GUI instead.
413 """
414 import matplotlib.pyplot as plt
415 from mpl_toolkits.mplot3d import Axes3D # NOQA
416 fig = plt.figure(figsize=(6, 6))
417 ax = fig.add_subplot(111, projection='3d')
418 self.make_plot(ax, alpha=alpha, linewidth=linewidth, colors=colors)
419 plt.subplots_adjust(left=0, bottom=0, right=1, top=1)
420 if legend: 420 ↛ 422line 420 didn't jump to line 422 because the condition on line 420 was always true
421 fig.legend(frameon=False)
422 plt.axis('off')
423 if save_as: 423 ↛ 424line 423 didn't jump to line 424 because the condition on line 423 was never true
424 plt.savefig(save_as, transparent=True)
425 else:
426 plt.show()
428 def write(self, filename: str):
429 """
430 Write particle to file. The file format is derived from the
431 filename. Currently supported fileformats are:
433 * Wavefront .obj
435 Parameters
436 ---------
437 filename
438 Filename of file to write to
439 """
440 supported_file_formats = ['obj']
441 fileformat = filename.split('.')[-1]
442 if fileformat not in supported_file_formats:
443 raise ValueError('File format {} not supported, '.format(fileformat) +
444 'supported formats are: ' + ' '.join(supported_file_formats))
445 with open(filename, 'w') as f:
446 if fileformat == 'obj': 446 ↛ exitline 446 didn't jump to the function exit
447 f.write('# Vertices\n')
448 for form in self.forms:
449 for facet in form.facets:
450 for vertex in facet.ordered_vertices[:-1]:
451 f.write('v {:15.8f} {:15.8f} {:15.8f}\n'.format(*vertex))
453 f.write('# Faces\n')
454 vertex_counter = 1
455 for form in self.forms:
456 # Define a group to be able to color each form separately
457 f.write('g {}\n'.format(form.parent_miller_indices))
458 for facet in form.facets:
459 # Tie the face to the right vertices. Count starts from 1.
460 s = ['{:5d}'.format(v + vertex_counter) for v in range(len(facet.vertices))]
461 vertex_counter += len(facet.vertices)
462 f.write('f ' + ' '.join(s) + '\n')
464 def _yield_facets(self, only_original_grains: bool = False) -> Facet:
465 """Generate all facets that are included in the facets dictionary."""
466 for form in self.forms:
467 for facet in form.facets:
468 if only_original_grains and not facet.original_grain:
469 continue
470 yield facet
472 @property
473 def volume(self) -> float:
474 """Returns the volume of the particle"""
475 volume = 0
476 for form in self.forms:
477 volume += form.volume
478 return volume
480 @property
481 def area(self) -> float:
482 """
483 Returns total area of the surface of the particle (not including
484 twin boundaries).
485 """
486 area = 0
487 for form in self.forms:
488 if form.parent_miller_indices != 'twin':
489 area += form.area
490 return area
492 @property
493 def edge_length(self) -> float:
494 """Returns total edge length of the particle."""
495 edge_length = 0
496 for form in self.forms:
497 if form.parent_miller_indices != 'twin':
498 edge_length += form.edge_length / 2 # Every edge comes twice
499 return edge_length
501 @property
502 def number_of_corners(self) -> float:
503 """Returns the number of corners (vertices) on the particle."""
504 unique_vertices = []
505 for form in self.forms:
506 if form.parent_miller_indices == 'twin':
507 continue
508 for facet in form.facets:
509 for vertex in facet.vertices:
510 if not is_array_in_arrays(vertex, unique_vertices):
511 unique_vertices.append(vertex)
512 return len(unique_vertices)
514 @property
515 def surface_energy(self) -> float:
516 """
517 The total surface energy of the particle (including twin boundaries).
518 """
519 E = 0
520 for form in self.forms:
521 E += form.surface_energy
522 return E
524 @property
525 def facet_fractions(self) -> Dict[tuple, float]:
526 """
527 Returns a dict specifying fraction of each form
528 (not including twin boundaries).
529 """
530 facet_fractions = {}
531 total_area = 0
532 for form in self.forms:
533 if form.parent_miller_indices == 'twin':
534 continue
535 area = form.area
536 total_area += area
537 facet_fractions[form.parent_miller_indices] = facet_fractions.get(
538 form.parent_miller_indices, 0) + area
539 for form in facet_fractions:
540 facet_fractions[form] = facet_fractions[form] / total_area
541 return facet_fractions
543 @property
544 def average_surface_energy(self) -> float:
545 """
546 Average surface energy for the Wulff construction, i.e.,
547 a weighted average over all the facets, where the weights are
548 the area fraction of each facet.
549 """
550 fractions = self.facet_fractions
551 weighted_average = 0
552 used_miller_indices = []
553 for form in self.forms:
554 if form.parent_miller_indices == 'twin' or \
555 form.parent_miller_indices in used_miller_indices:
556 continue
557 weighted_average += form.energy * fractions[form.parent_miller_indices]
558 used_miller_indices.append(form.parent_miller_indices)
559 return weighted_average
561 def _duplicate_particle(self, symmetries):
562 """
563 Duplicate the particle by applying symmetry operations.
565 Useful for making an icosahedron from a single tetrahedron for example.
567 Parameters
568 ----------
569 symmetries: list of NumPy arrays
570 Each array a symmetry operation(operating on Cartesian coordinates)
571 that should create a duplicate.
572 """
573 for form in self.forms:
574 new_facets = []
575 for facet in form.facets:
576 for R in symmetries:
577 new_facet = Facet(normal=np.dot(R, facet.normal),
578 energy=facet.energy,
579 symmetry=None,
580 original_grain=False)
581 for vertex in facet.vertices:
582 new_facet.add_vertex(np.dot(R, vertex))
584 new_facets.append(new_facet)
585 for facet in new_facets:
586 form.facets.append(facet)
588 def translate_particle(self, translation: np.ndarray):
589 """
590 Translate the particle.
592 Parameters
593 ----------
594 translation: list of 3 floats
595 Translation vector
596 """
597 for facet in self._yield_facets():
598 for i, vertex in enumerate(facet.vertices):
599 facet.vertices[i] = vertex + translation
601 def rotate_particle(self, rotation: np.ndarray):
602 """
603 Rotate the particle.
605 Parameters
606 ----------
607 rotation
608 Rotation matrix
609 """
610 if abs(abs(np.linalg.det(rotation)) - 1.0) > self.tol:
611 raise ValueError('Provided matrix is not a rotation matrix '
612 '(its determinant does not equal 1 or -1).')
613 for facet in self._yield_facets():
614 for i, vertex in enumerate(facet.vertices):
615 facet.vertices[i] = np.dot(rotation, vertex)
616 facet.normal = np.dot(rotation, facet.normal)
617 facet.original_normal = np.dot(rotation, facet.original_normal)
618 self.standardized_structure.cell = np.dot(rotation, self.standardized_structure.cell.T).T
619 for atom in self.standardized_structure:
620 atom.position = np.dot(rotation, atom.position)
622 def _get_atoms(self, center_shift: Tuple[float, float, float] = None) -> Atoms:
623 """
624 Returns an ASE Atoms object the atoms of which are all
625 within the facets. In polycrystalline particles, only
626 one grain will be returned.
628 Parameters
629 ----------
630 center_shift
631 Shift of center in Cartesian coordinates.
632 """
633 # Find max and min in x, y and z to bound the volume
634 # where there may be atoms
635 # (using scaled coordinates relative to the origin)
636 mins = [1e9, 1e9, 1e9]
637 maxs = [-1e9, -1e9, -1e9]
638 for facet in self._yield_facets(only_original_grains=True):
639 if facet.original_vertices:
640 vertices = facet.original_vertices
641 else:
642 vertices = facet.vertices
643 for vertex in vertices:
644 scaled_vertex = np.linalg.solve(self.standardized_structure.cell.T, vertex)
645 for i, v in enumerate(scaled_vertex):
646 if v < mins[i]:
647 mins[i] = v
648 if v > maxs[i]:
649 maxs[i] = v
651 # Prepare an unnecessarily large atoms object
652 mins = [int(i) - 1 for i in mins]
653 maxs = [int(np.ceil(i)) + 1 for i in maxs]
654 repeat = [int(np.ceil(i)) - int(j) for i, j in zip(maxs, mins)]
655 atoms = self.standardized_structure.repeat(repeat)
656 atoms.translate(np.dot(self.standardized_structure.cell.T, mins))
658 if center_shift:
659 atoms.translate(center_shift)
661 # Loop over all scaled coordinates and check whether
662 # they are within all facets
663 to_delete = []
664 for atom in atoms:
665 for facet in self._yield_facets(only_original_grains=True):
666 if np.dot(atom.position, facet.original_normal) > \
667 facet.distance_from_origin + 1e-5:
668 to_delete.append(atom.index)
669 break
670 del atoms[to_delete]
671 atoms.set_pbc(False)
672 atoms.set_cell([0., 0., 0.])
673 return atoms