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