Coverage for wulffpack/core/base_particle.py: 97%
361 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-29 11:40 +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 __str__(self) -> str:
77 width = 38
78 s = []
79 name = self.__class__.__name__
80 s += ['{s:-^{n}}'.format(s=f' {name} ', n=width)]
81 flds = dict(volume='12.4f', area='12.4f',
82 surface_energy='12.2f', average_surface_energy='12.4f',
83 edge_length='12.4f', number_of_corners='12d')
84 for fld, fmt in flds.items():
85 s += ['{fld:22} : {value:{fmt}}'.format(fld=fld, value=getattr(self, fld), fmt=fmt)]
86 s += ['facet_fractions :']
87 for k, v in self.facet_fractions.items():
88 s += [f' {str(k):20} : {v:12.6f}']
89 return '\n'.join(s)
91 def _repr_html_(self) -> str:
92 s = []
93 s += ['<table border="1" class="dataframe"']
94 s += [
95 '<thead><tr><th style="text-align: left;">Property</th><th>Value</th></tr></thead>'
96 ]
97 s += ['<tbody>']
98 flds = dict(volume='12.4f', area='12.4f',
99 surface_energy='12.2f', average_surface_energy='12.4f',
100 edge_length='12.4f', number_of_corners='12d')
101 for fld, fmt in flds.items():
102 s += [
103 f'<tr><td style="text-align: left">{fld:22}</td>'
104 '<td>{value:{fmt}}</td><tr>'.format(value=getattr(self, fld), fmt=fmt)
105 ]
106 s += '<tr><td style="text-align: left">facet_fractions</td><tr>'
107 for k, v in self.facet_fractions.items():
108 s += [
109 '<tr><td style="text-align: right">'
110 f'{str(k):22}</td><td>{v:.6f}</td><tr>'
111 ]
112 s += ['</tbody>']
113 s += ['</table>']
114 return ''.join(s)
116 def _closest_vertices(self, facets):
117 """
118 The closest vertex technique calculates all possible vertices and eliminates
119 vertices outside of the Wulff shape.
120 This works for negative interface energies.
121 """
122 # Calculate intersection point for all possible combinations of three facets
123 facets.sort(key=lambda x: x.energy) # compare with small energies first in inner loop
124 for i, fi in enumerate(facets):
125 ni = fi.normal
126 for j, fj in enumerate(facets[:i]):
127 nj = fj.normal
128 if abs(abs(np.dot(ni, nj)) - 1) < self._SMALL:
129 continue # skip parallel facets
130 for fk in facets[:j]:
131 mat = np.vstack([ni, nj, fk.normal])
132 if np.linalg.det(mat) == 0:
133 continue # skip singular matrices
134 vertex = np.dot(np.linalg.inv(mat), [fi.energy, fj.energy, fk.energy])
136 # check if vertex lies outside of Wulff shape
137 for fl in facets:
138 projection = np.dot(vertex, fl.normal)
139 if projection > fl.energy + self._SMALL:
140 break
141 else: # i.e. loop does not break
142 for f in [fi, fj, fk]:
143 f.add_vertex(vertex.copy())
145 def _dual_wulff_shape(self, facets):
146 """
147 The Wulff construction calculates the vertices of the dual,
148 calculates the convex hull of those vertices using scipy's
149 ConvexHull, and then converts the dual back to the Wulff
150 shape.
151 """
152 # Calculate convex hull of dual vertices
153 duals = []
154 for facet in self._yield_facets():
155 facet_points = facet.normal * facet.energy
156 duals.append(facet_points / np.dot(facet_points, facet_points))
157 hull = ConvexHull(duals)
159 # Calculate vertices of dual again (i.e., back to Wulff shape)
160 for i, j, k in hull.simplices:
161 normal = np.cross(duals[j] - duals[i], duals[k] - duals[i])
162 if np.linalg.norm(normal) < self._SMALL: 162 ↛ 163line 162 didn't jump to line 163 because the condition on line 162 was never true
163 continue # Rare instances of duals on a line
164 if np.isclose(np.dot(normal, facets[i].normal), 0):
165 raise ValueError('The vertex after the calculation of the dual and the'
166 ' original facet are orthogonal. This indicates that the'
167 ' facets that are provided cannot form a closed shape.'
168 ' Please review your input.')
169 normal *= facets[i].energy / np.dot(normal, facets[i].normal)
170 for facet_i in (i, j, k):
171 facets[facet_i].add_vertex(normal.copy())
173 def _scale_size(self, scale_factor: float):
174 """
175 Scale the size of the particle.
177 Parameters
178 ----------
179 scale_factor
180 Factor that all coordinates will be scaled with.
181 """
182 for facet in self._yield_facets():
183 for i, _ in enumerate(facet.vertices):
184 facet.vertices[i] *= scale_factor
185 for i, _ in enumerate(facet.original_vertices):
186 facet.original_vertices[i] *= scale_factor
188 @property
189 def natoms(self) -> List[int]:
190 """
191 The approximate number of atoms in the particle
192 (implicitly defining the volume).
193 """
194 return self._natoms
196 @natoms.setter
197 def natoms(self, natoms: int):
198 """
199 Change the approximate number of atoms and thus the volume.
200 """
201 scale_factor = (natoms / self._natoms) ** (1 / 3.)
202 self._scale_size(scale_factor)
203 self._natoms = natoms
205 @property
206 def forms(self) -> List[Form]:
207 """List of inequivalent forms for the particle"""
208 return self._forms
210 @property
211 def standardized_structure(self) -> Atoms:
212 """
213 The standardized atomic structure that defines the geometry
214 and thus the meaning of the Miller indices. Also forms the building
215 blocks when `particle.atoms` is called.
216 """
217 return self._standardized_structure
219 @property
220 def _twin_form(self) -> Form:
221 """Returns the twin form if there is one, otherwise None."""
222 for form in self.forms:
223 if form.parent_miller_indices == 'twin':
224 return form
225 return None
227 def get_continuous_color_scheme(self,
228 base_colors: dict = None,
229 normalize: bool = False) -> dict:
230 """
231 Returns a dictionary with RGB colors for each form.
232 The colors smoothly interpolate between three
233 base colors, corresponding to (1, 1, 1), (1, 1, 0) and
234 (1, 0, 0). Note that this is sensible primarily for
235 cubic systems.
237 Parameters
238 ----------
239 base_colors
240 User chosen colors for one or several of (1, 1, 1),
241 (1, 1, 0) and (1, 0, 0). To enforce, say, green
242 (1, 1, 1), use ``base_colors={(1, 1, 1): 'g'}``
243 normalize
244 If True, the norm of the RGB vectors will be 1. Note
245 that this may affect the ``base_colors`` too.
246 """
247 from matplotlib.colors import to_rgba
249 # Adapt base_colors dictionary
250 new_base_colors = {}
251 if base_colors:
252 for key in base_colors.keys():
253 adapted_key = tuple(sorted(abs(i) for i in key))
254 new_base_colors[adapted_key] = base_colors[key]
256 # Make sure base_colors is complete
257 default_colors = {(1, 1, 1): '#2980B9',
258 (0, 1, 1): '#E92866',
259 (0, 0, 1): '#FFE82C'}
260 for key, color in default_colors.items():
261 if key not in new_base_colors:
262 new_base_colors[key] = color
263 A = [(1, 1, 1), (0, 1, 1), (0, 0, 1)]
264 base_color_matrix = np.array([to_rgba(new_base_colors[f])[:3] for f in A]).T
266 # Normalize rows in A
267 A = np.array([np.array(miller) / np.linalg.norm(miller) for miller in A]).T
269 # Calculate color for each form
270 colors = {}
271 for form in self.forms:
272 if type(form.parent_miller_indices) is not tuple:
273 continue
274 miller = np.array(sorted(abs(i) for i in form.parent_miller_indices))
275 color = np.linalg.solve(A, miller)
276 color /= np.linalg.norm(color)
277 color = np.dot(base_color_matrix, color)
278 if normalize:
279 color /= np.linalg.norm(color)
280 else:
281 color = [min(1, c) for c in color]
282 colors[form.parent_miller_indices] = color
283 return colors
285 def make_plot(self, ax, alpha: float = 0.85, linewidth: float = 0.3, colors: dict = None):
286 """
287 Plot a particle in an axis object. This function can be used to
288 make customized plots of particles.
290 Parameters
291 ----------
292 ax : matplotlib Axes3DSubplot
293 An axis object with 3d projection
294 alpha
295 Opacity of the faces
296 linewidth
297 Thickness of lines between faces
298 colors
299 Allows custom colors for facets of all or a subset of forms,
300 example `{(1, 1, 1): '#FF0000'}`
302 Example
303 -------
304 In the following example, three different particles are
305 plotted in the same figure::
307 >>> from wulffpack import SingleCrystal, Decahedron, Icosahedron
308 >>> import matplotlib.pyplot as plt
309 >>> from mpl_toolkits.mplot3d import Axes3D
310 >>>
311 >>> surface_energies = {(1, 1, 1): 1.0,
312 ... (1, 0, 0): 1.1,
313 ... (1, 1, 0): 1.15,
314 ... (3, 2, 1): 1.15}
315 >>> twin_energy = 0.05
316 >>>
317 >>> fig = plt.figure(figsize=(3*4.0, 4.0))
318 >>> ax = fig.add_subplot(131, projection='3d')
319 >>> particle = SingleCrystal(surface_energies)
320 >>> particle.make_plot(ax)
321 >>>
322 >>> ax = fig.add_subplot(132, projection='3d')
323 >>> particle = Decahedron(surface_energies,
324 ... twin_energy=0.05)
325 >>> particle.make_plot(ax)
326 >>>
327 >>> ax = fig.add_subplot(133, projection='3d')
328 >>> particle = Icosahedron(surface_energies,
329 ... twin_energy=0.05)
330 >>> particle.make_plot(ax)
331 >>>
332 >>> plt.subplots_adjust(top=1, bottom=0, left=0,
333 ... right=1, wspace=0, hspace=0)
334 >>> plt.savefig('particles.png')
336 """
337 from mpl_toolkits.mplot3d.art3d import (Poly3DCollection,
338 Line3DCollection)
340 # Standardize color dict
341 if colors is None:
342 colors = {}
344 default_colors = ['#D62728',
345 '#E377C2',
346 '#8C564B',
347 '#7F7F7F',
348 '#9467BD',
349 '#BCBD22',
350 '#17BECF',
351 '#AEC7E8',
352 '#FFBB78',
353 '#98DF8A',
354 '#FF9896',
355 '#C5B0D5',
356 '#C49C94',
357 '#F7B6D2',
358 '#C7C7C7',
359 '#DBDB8D',
360 '#9EDAE5']
362 mins = [1e9, 1e9, 1e9]
363 maxs = [-1e9, -1e9, -1e9]
364 poly3d = []
366 used_forms = []
367 color_counter = 0
368 for form in self.forms:
369 if form.parent_miller_indices == 'twin':
370 continue
372 # Determine color
373 if form.parent_miller_indices in colors:
374 color = colors[form.parent_miller_indices]
375 elif form.parent_miller_indices == (1, 1, 1):
376 color = '#4f81f1'
377 colors[form.parent_miller_indices] = color
378 elif tuple(sorted(form.parent_miller_indices)) == (0, 0, 1): 378 ↛ 381line 378 didn't jump to line 381 because the condition on line 378 was always true
379 color = '#f8c73b'
380 colors[form.parent_miller_indices] = color
381 elif tuple(sorted(form.parent_miller_indices)) == (0, 1, 1):
382 color = '#2CA02C'
383 colors[form.parent_miller_indices] = color
384 else:
385 color = default_colors[color_counter % len(default_colors)]
386 colors[form.parent_miller_indices] = color
387 color_counter += 1
389 # Save the used forms to be able to make a legend
390 # (the colors dict may contain forms that are not present
391 # if it was supplied by the user)
392 if form.parent_miller_indices not in used_forms: 392 ↛ 396line 392 didn't jump to line 396 because the condition on line 392 was always true
393 used_forms.append(form.parent_miller_indices)
395 # Plot all facets in the form
396 for facet in form.facets:
397 poly3d = [facet.ordered_vertices]
399 # Facets
400 collection = Poly3DCollection(poly3d, alpha=alpha)
401 collection.set_facecolor(color)
402 ax.add_collection3d(collection)
404 # Lines
405 ax.add_collection3d(Line3DCollection(poly3d, colors='k', linewidths=linewidth))
407 # Find proper ranges for the axes
408 for vertex in facet.ordered_vertices:
409 for i, v in enumerate(vertex):
410 if v > maxs[i]:
411 maxs[i] = v
412 if v < mins[i]:
413 mins[i] = v
415 # Plot fake lines to make a legend
416 for form in used_forms:
417 if isinstance(form, tuple):
418 label = ('{{' + '{}' * len(form) + '}}').format(*form)
419 else:
420 label = form
421 ax.plot([0, 0], [0, 0], color=colors[form], linewidth=8.0, label=label)
423 # Fiddle with the layout
424 extent = max(maxs) - min(mins)
425 bounds = [min(mins) + 0.15 * extent, max(maxs) - 0.15 * extent]
426 ax.set_xlim(*bounds)
427 ax.set_ylim(*bounds)
428 ax.set_zlim(*bounds)
429 ax.set_xticks([])
430 ax.set_yticks([])
431 ax.set_zticks([])
432 ax.axis('off')
433 try:
434 ax.set_box_aspect((1, 1, 1))
435 except AttributeError:
436 # set_box_aspect only available in recent versions of matplotlib
437 pass
439 def view(self, alpha: float = 0.85, linewidth: float = 0.3, colors: dict = None,
440 legend: bool = True, save_as: str = None):
441 """
442 Use matplotlib to view a rendition of the particle.
444 Parameters
445 ----------
446 alpha
447 Opacity of the faces
448 linewidth
449 Thickness of lines between faces
450 colors
451 Allows custom colors for facets of all or a subset of forms,
452 example `{(1, 1, 1): '#FF0000'}`
453 legend
454 Whether or not to show a legend with facet-color definitions
455 save_as
456 Filename to save figure as. If None, show the particle
457 with the GUI instead.
458 """
459 import matplotlib.pyplot as plt
460 from mpl_toolkits.mplot3d import Axes3D # NOQA
461 fig = plt.figure(figsize=(6, 6))
462 ax = fig.add_subplot(111, projection='3d')
463 self.make_plot(ax, alpha=alpha, linewidth=linewidth, colors=colors)
464 plt.subplots_adjust(left=0, bottom=0, right=1, top=1)
465 if legend: 465 ↛ 467line 465 didn't jump to line 467 because the condition on line 465 was always true
466 fig.legend(frameon=False)
467 plt.axis('off')
468 if save_as: 468 ↛ 469line 468 didn't jump to line 469 because the condition on line 468 was never true
469 plt.savefig(save_as, transparent=True)
470 else:
471 plt.show()
473 def write(self, filename: str):
474 """
475 Write particle to file. The file format is derived from the
476 filename. Currently supported fileformats are:
478 * Wavefront .obj
480 Parameters
481 ---------
482 filename
483 Filename of file to write to
484 """
485 supported_file_formats = ['obj']
486 fileformat = filename.split('.')[-1]
487 if fileformat not in supported_file_formats:
488 raise ValueError('File format {} not supported, '.format(fileformat) +
489 'supported formats are: ' + ' '.join(supported_file_formats))
490 with open(filename, 'w') as f:
491 if fileformat == 'obj': 491 ↛ exitline 491 didn't jump to the function exit
492 f.write('# Vertices\n')
493 for form in self.forms:
494 for facet in form.facets:
495 for vertex in facet.ordered_vertices[:-1]:
496 f.write('v {:15.8f} {:15.8f} {:15.8f}\n'.format(*vertex))
498 f.write('# Faces\n')
499 vertex_counter = 1
500 for form in self.forms:
501 # Define a group to be able to color each form separately
502 f.write('g {}\n'.format(form.parent_miller_indices))
503 for facet in form.facets:
504 # Tie the face to the right vertices. Count starts from 1.
505 s = ['{:5d}'.format(v + vertex_counter) for v in range(len(facet.vertices))]
506 vertex_counter += len(facet.vertices)
507 f.write('f ' + ' '.join(s) + '\n')
509 def _yield_facets(self, only_original_grains: bool = False) -> Facet:
510 """Generate all facets that are included in the facets dictionary."""
511 for form in self.forms:
512 for facet in form.facets:
513 if only_original_grains and not facet.original_grain:
514 continue
515 yield facet
517 @property
518 def volume(self) -> float:
519 """Returns the volume of the particle"""
520 volume = 0
521 for form in self.forms:
522 volume += form.volume
523 return volume
525 @property
526 def area(self) -> float:
527 """
528 Returns total area of the surface of the particle (not including
529 twin boundaries).
530 """
531 area = 0
532 for form in self.forms:
533 if form.parent_miller_indices != 'twin':
534 area += form.area
535 return area
537 @property
538 def edge_length(self) -> float:
539 """Returns total edge length of the particle."""
540 edge_length = 0
541 for form in self.forms:
542 if form.parent_miller_indices != 'twin':
543 edge_length += form.edge_length / 2 # Every edge comes twice
544 return edge_length
546 @property
547 def number_of_corners(self) -> float:
548 """Returns the number of corners (vertices) on the particle."""
549 unique_vertices = []
550 for form in self.forms:
551 if form.parent_miller_indices == 'twin':
552 continue
553 for facet in form.facets:
554 for vertex in facet.vertices:
555 if not is_array_in_arrays(vertex, unique_vertices):
556 unique_vertices.append(vertex)
557 return len(unique_vertices)
559 @property
560 def surface_energy(self) -> float:
561 """
562 The total surface energy of the particle (including twin boundaries).
563 """
564 E = 0
565 for form in self.forms:
566 E += form.surface_energy
567 return E
569 @property
570 def facets(self) -> List[Facet]:
571 """The facets of the particle (read-only).
572 """
573 return list(self._yield_facets())
575 @property
576 def facet_fractions(self) -> Dict[tuple, float]:
577 """
578 Returns a dict specifying fraction of each form
579 (not including twin boundaries).
580 """
581 facet_fractions = {}
582 total_area = 0
583 for form in self.forms:
584 if form.parent_miller_indices == 'twin':
585 continue
586 area = form.area
587 total_area += area
588 facet_fractions[form.parent_miller_indices] = facet_fractions.get(
589 form.parent_miller_indices, 0) + area
590 for form in facet_fractions:
591 facet_fractions[form] = facet_fractions[form] / total_area
592 return facet_fractions
594 @property
595 def average_surface_energy(self) -> float:
596 """
597 Average surface energy for the Wulff construction, i.e.,
598 a weighted average over all the facets, where the weights are
599 the area fraction of each facet.
600 """
601 fractions = self.facet_fractions
602 weighted_average = 0
603 used_miller_indices = []
604 for form in self.forms:
605 if form.parent_miller_indices == 'twin' or \
606 form.parent_miller_indices in used_miller_indices:
607 continue
608 weighted_average += form.energy * fractions[form.parent_miller_indices]
609 used_miller_indices.append(form.parent_miller_indices)
610 return weighted_average
612 def _duplicate_particle(self, symmetries):
613 """
614 Duplicate the particle by applying symmetry operations.
616 Useful for making an icosahedron from a single tetrahedron for example.
618 Parameters
619 ----------
620 symmetries: list of NumPy arrays
621 Each array a symmetry operation(operating on Cartesian coordinates)
622 that should create a duplicate.
623 """
624 for form in self.forms:
625 new_facets = []
626 for facet in form.facets:
627 for R in symmetries:
628 new_facet = Facet(normal=np.dot(R, facet.normal),
629 energy=facet.energy,
630 symmetry=None,
631 original_grain=False)
632 for vertex in facet.vertices:
633 new_facet.add_vertex(np.dot(R, vertex))
635 new_facets.append(new_facet)
636 for facet in new_facets:
637 form.facets.append(facet)
639 def translate_particle(self, translation: np.ndarray):
640 """
641 Translate the particle.
643 Parameters
644 ----------
645 translation: list of 3 floats
646 Translation vector
647 """
648 for facet in self._yield_facets():
649 for i, vertex in enumerate(facet.vertices):
650 facet.vertices[i] = vertex + translation
652 def rotate_particle(self, rotation: np.ndarray):
653 """
654 Rotate the particle.
656 Parameters
657 ----------
658 rotation
659 Rotation matrix
660 """
661 if abs(abs(np.linalg.det(rotation)) - 1.0) > self.tol:
662 raise ValueError('Provided matrix is not a rotation matrix '
663 '(its determinant does not equal 1 or -1).')
664 for facet in self._yield_facets():
665 for i, vertex in enumerate(facet.vertices):
666 facet.vertices[i] = np.dot(rotation, vertex)
667 facet.normal = np.dot(rotation, facet.normal)
668 facet.original_normal = np.dot(rotation, facet.original_normal)
669 self.standardized_structure.cell = np.dot(rotation, self.standardized_structure.cell.T).T
670 for atom in self.standardized_structure:
671 atom.position = np.dot(rotation, atom.position)
673 def _get_atoms(self, center_shift: Tuple[float, float, float] = None) -> Atoms:
674 """
675 Returns an ASE Atoms object the atoms of which are all
676 within the facets. In polycrystalline particles, only
677 one grain will be returned.
679 Parameters
680 ----------
681 center_shift
682 Shift of center in Cartesian coordinates.
683 """
684 # Find max and min in x, y and z to bound the volume
685 # where there may be atoms
686 # (using scaled coordinates relative to the origin)
687 mins = [1e9, 1e9, 1e9]
688 maxs = [-1e9, -1e9, -1e9]
689 for facet in self._yield_facets(only_original_grains=True):
690 if facet.original_vertices:
691 vertices = facet.original_vertices
692 else:
693 vertices = facet.vertices
694 for vertex in vertices:
695 scaled_vertex = np.linalg.solve(self.standardized_structure.cell.T, vertex)
696 for i, v in enumerate(scaled_vertex):
697 if v < mins[i]:
698 mins[i] = v
699 if v > maxs[i]:
700 maxs[i] = v
702 # Prepare an unnecessarily large atoms object
703 mins = [int(i) - 1 for i in mins]
704 maxs = [int(np.ceil(i)) + 1 for i in maxs]
705 repeat = [int(np.ceil(i)) - int(j) for i, j in zip(maxs, mins)]
706 atoms = self.standardized_structure.repeat(repeat)
707 atoms.translate(np.dot(self.standardized_structure.cell.T, mins))
709 if center_shift:
710 atoms.translate(center_shift)
712 # Loop over all scaled coordinates and check whether
713 # they are within all facets
714 to_delete = []
715 for atom in atoms:
716 for facet in self._yield_facets(only_original_grains=True):
717 if np.dot(atom.position, facet.original_normal) > \
718 facet.distance_from_origin + 1e-5:
719 to_delete.append(atom.index)
720 break
721 del atoms[to_delete]
722 atoms.set_pbc(False)
723 atoms.set_cell([0., 0., 0.])
724 return atoms