Coverage for icet/core/cluster_space.py: 97%
312 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1"""
2This module provides the :class:`ClusterSpace` class.
3"""
5import os
6import copy
7import itertools
8import pickle
9import tarfile
10import tempfile
11from collections.abc import Iterable
12from math import log10, floor
13from typing import Dict, List, Union, Tuple
15import numpy as np
16import spglib
18from _icet import ClusterSpace as _ClusterSpace
19from ase import Atoms
20from ase.io import read as ase_read
21from ase.io import write as ase_write
22from icet.core.orbit_list import OrbitList
23from icet.core.structure import Structure
24from icet.core.sublattices import Sublattices
25from icet.tools.geometry import (ase_atoms_to_spglib_cell,
26 get_occupied_primitive_structure)
27from pandas import DataFrame
30class ClusterSpace(_ClusterSpace):
31 """This class provides functionality for generating and maintaining
32 cluster spaces.
34 Note
35 ----
36 In :program:`icet` all :class:`Atoms <ase.Atoms>` objects must have
37 periodic boundary conditions. When constructing cluster expansions
38 for surfaces and nanoparticles it is therefore recommended to
39 surround the structure with vacuum and use periodic boundary
40 conditions. This can be achieved by using :func:`Atoms.center <ase.Atoms.center>`.
42 Parameters
43 ----------
44 structure
45 Atomic structure.
46 cutoffs
47 Cutoff radii per order that define the cluster space.
49 Cutoffs are specified in units of Ångstrom and refer to the
50 longest distance between two atoms in the cluster. The first
51 element refers to pairs, the second to triplets, the third
52 to quadruplets, and so on. :attr:`cutoffs=[7.0, 4.5]` thus implies
53 that all pairs distanced 7 Å or less will be included,
54 as well as all triplets among which the longest distance is no
55 longer than 4.5 Å.
56 chemical_symbols
57 List of chemical symbols, each of which must map to an element
58 of the periodic table.
60 If a list of chemical symbols is provided, all sites on the
61 lattice will have the same allowed occupations as the input
62 list.
64 If a list of list of chemical symbols is provided then the
65 outer list must be the same length as the :attr:`structure` object and
66 :attr:`chemical_symbols[i]` will correspond to the allowed species
67 on lattice site ``i``.
68 symprec
69 Tolerance imposed when analyzing the symmetry using spglib.
70 position_tolerance
71 Tolerance applied when comparing positions in Cartesian coordinates.
73 Examples
74 --------
75 The following snippets illustrate several common situations::
77 >>> from ase.build import bulk
78 >>> from ase.io import read
79 >>> from icet import ClusterSpace
81 >>> # AgPd alloy with pairs up to 7.0 A and triplets up to 4.5 A
82 >>> prim = bulk('Ag')
83 >>> cs = ClusterSpace(structure=prim, cutoffs=[7.0, 4.5],
84 ... chemical_symbols=[['Ag', 'Pd']])
85 >>> print(cs)
87 >>> # (Mg,Zn)O alloy on rocksalt lattice with pairs up to 8.0 A
88 >>> prim = bulk('MgO', crystalstructure='rocksalt', a=6.0)
89 >>> cs = ClusterSpace(structure=prim, cutoffs=[8.0],
90 ... chemical_symbols=[['Mg', 'Zn'], ['O']])
91 >>> print(cs)
93 >>> # (Ga,Al)(As,Sb) alloy with pairs, triplets, and quadruplets
94 >>> prim = bulk('GaAs', crystalstructure='zincblende', a=6.5)
95 >>> cs = ClusterSpace(structure=prim, cutoffs=[7.0, 6.0, 5.0],
96 ... chemical_symbols=[['Ga', 'Al'], ['As', 'Sb']])
97 >>> print(cs)
99 >>> # PdCuAu alloy with pairs and triplets
100 >>> prim = bulk('Pd')
101 >>> cs = ClusterSpace(structure=prim, cutoffs=[7.0, 5.0],
102 ... chemical_symbols=[['Au', 'Cu', 'Pd']])
103 >>> print(cs)
105 """
107 def __init__(self,
108 structure: Atoms,
109 cutoffs: List[float],
110 chemical_symbols: Union[List[str], List[List[str]]],
111 symprec: float = 1e-5,
112 position_tolerance: float = None) -> None:
114 if not isinstance(structure, Atoms): 114 ↛ 115line 114 didn't jump to line 115, because the condition on line 114 was never true
115 raise TypeError('Input configuration must be an ASE Atoms object'
116 f', not type {type(structure)}.')
117 if not all(structure.pbc):
118 raise ValueError('Input structure must be periodic.')
119 if symprec <= 0:
120 raise ValueError('symprec must be a positive number.')
122 self._config = {'symprec': symprec}
123 self._cutoffs = cutoffs.copy()
124 self._input_structure = structure.copy()
125 self._input_chemical_symbols = copy.deepcopy(chemical_symbols)
126 chemical_symbols = self._get_chemical_symbols()
128 self._pruning_history: List[tuple] = []
130 # set up primitive
131 occupied_primitive, primitive_chemical_symbols = get_occupied_primitive_structure(
132 self._input_structure, chemical_symbols, symprec=self.symprec)
133 self._primitive_chemical_symbols = primitive_chemical_symbols
134 assert len(occupied_primitive) == len(primitive_chemical_symbols)
136 # derived tolerances
137 if position_tolerance is None:
138 self._config['position_tolerance'] = symprec
139 else:
140 if position_tolerance <= 0:
141 raise ValueError('position_tolerance must be a positive number')
142 self._config['position_tolerance'] = position_tolerance
143 effective_box_size = abs(np.linalg.det(occupied_primitive.cell)) ** (1 / 3)
144 tol = self.position_tolerance / effective_box_size
145 tol = min(tol, self._config['position_tolerance'] / 5)
146 self._config['fractional_position_tolerance'] = round(tol, -int(floor(log10(abs(tol)))))
148 # set up orbit list
149 self._orbit_list = OrbitList(
150 structure=occupied_primitive,
151 cutoffs=self._cutoffs,
152 chemical_symbols=self._primitive_chemical_symbols,
153 symprec=self.symprec,
154 position_tolerance=self.position_tolerance,
155 fractional_position_tolerance=self.fractional_position_tolerance)
156 self._orbit_list.remove_orbits_with_inactive_sites()
158 # call (base) C++ constructor
159 _ClusterSpace.__init__(self,
160 orbit_list=self._orbit_list,
161 position_tolerance=self.position_tolerance,
162 fractional_position_tolerance=self.fractional_position_tolerance)
164 def _get_chemical_symbols(self):
165 """ Returns chemical symbols using input structure and
166 chemical symbols. Carries out multiple sanity checks. """
168 # setup chemical symbols as List[List[str]]
169 if all(isinstance(i, str) for i in self._input_chemical_symbols):
170 chemical_symbols = [self._input_chemical_symbols] * len(self._input_structure)
171 # also accept tuples and other iterables but not, e.g., List[List, str]
172 # (need to check for str explicitly here because str is an Iterable)
173 elif not all(isinstance(i, Iterable) and not isinstance(i, str)
174 for i in self._input_chemical_symbols):
175 raise TypeError('chemical_symbols must be List[str] or List[List[str]], not {}'.format(
176 type(self._input_chemical_symbols)))
177 elif len(self._input_chemical_symbols) != len(self._input_structure):
178 msg = 'chemical_symbols must have same length as structure. '
179 msg += 'len(chemical_symbols) = {}, len(structure)= {}'.format(
180 len(self._input_chemical_symbols), len(self._input_structure))
181 raise ValueError(msg)
182 else:
183 chemical_symbols = copy.deepcopy(self._input_chemical_symbols)
185 for i, symbols in enumerate(chemical_symbols):
186 if len(symbols) != len(set(symbols)):
187 raise ValueError(
188 'Found duplicates of allowed chemical symbols on site {}.'
189 ' allowed species on site {}= {}'.format(i, i, symbols))
191 if len([tuple(sorted(s)) for s in chemical_symbols if len(s) > 1]) == 0:
192 raise ValueError('No active sites found')
194 return chemical_symbols
196 def _get_chemical_symbol_representation(self):
197 """Returns a str version of the chemical symbols that is
198 easier on the eyes.
199 """
200 sublattices = self.get_sublattices(self.primitive_structure)
201 nice_str = []
202 for sublattice in sublattices.active_sublattices:
203 sublattice_symbol = sublattice.symbol
204 nice_str.append('{} (sublattice {})'.format(
205 list(sublattice.chemical_symbols), sublattice_symbol))
206 return ', '.join(nice_str)
208 def _get_string_representation(self,
209 print_threshold: int = None,
210 print_minimum: int = 10) -> str:
211 """
212 String representation of the cluster space that provides an overview of
213 the orbits (order, radius, multiplicity etc) that constitute the space.
215 Parameters
216 ----------
217 print_threshold
218 if the number of orbits exceeds this number print dots
219 print_minimum
220 number of lines printed from the top and the bottom of the orbit
221 list if `print_threshold` is exceeded
223 Returns
224 -------
225 multi-line string
226 string representation of the cluster space.
227 """
229 def repr_orbit(orbit, header=False):
230 formats = {'order': '{:2}',
231 'radius': '{:8.4f}',
232 'multiplicity': '{:4}',
233 'index': '{:4}',
234 'orbit_index': '{:4}',
235 'multicomponent_vector': '{:}',
236 'sublattices': '{:}'}
237 s = []
238 for name, value in orbit.items():
239 if name == 'sublattices':
240 str_repr = formats[name].format('-'.join(value))
241 else:
242 str_repr = formats[name].format(value)
243 n = max(len(name), len(str_repr))
244 if header:
245 s += ['{s:^{n}}'.format(s=name, n=n)]
246 else:
247 s += ['{s:^{n}}'.format(s=str_repr, n=n)]
248 return ' | '.join(s)
250 # basic information
251 # (use largest orbit to obtain maximum line length)
252 prototype_orbit = self.as_list[-1]
253 width = len(repr_orbit(prototype_orbit))
254 s = []
255 s += ['{s:=^{n}}'.format(s=' Cluster Space ', n=width)]
256 s += [' {:38} : {}'.format('space group', self.space_group)]
257 s += [' {:38} : {}'
258 .format('chemical species', self._get_chemical_symbol_representation())]
259 s += [' {:38} : {}'.format('cutoffs',
260 ' '.join(['{:.4f}'.format(c) for c in self.cutoffs]))]
261 s += [' {:38} : {}'.format('total number of parameters', len(self))]
262 t = ['{}= {}'.format(k, c)
263 for k, c in self.number_of_orbits_by_order.items()]
264 s += [' {:38} : {}'.format('number of parameters by order', ' '.join(t))]
265 for key, value in sorted(self._config.items()):
266 s += [' {:38} : {}'.format(key, value)]
268 # table header
269 s += [''.center(width, '-')]
270 s += [repr_orbit(prototype_orbit, header=True)]
271 s += [''.center(width, '-')]
273 # table body
274 index = 0
275 orbit_list_info = self.as_list
276 while index < len(orbit_list_info):
277 if (print_threshold is not None and
278 len(self) > print_threshold and
279 index >= print_minimum and
280 index <= len(self) - print_minimum):
281 index = len(self) - print_minimum
282 s += [' ...']
283 s += [repr_orbit(orbit_list_info[index])]
284 index += 1
285 s += [''.center(width, '=')]
287 return '\n'.join(s)
289 def __str__(self) -> str:
290 """ String representation. """
291 return self._get_string_representation(print_threshold=50)
293 def _repr_html_(self) -> str:
294 """ HTML representation. Used, e.g., in jupyter notebooks. """
295 s = ['<h4>Cluster Space</h4>']
296 s += ['<table border="1" class="dataframe">']
297 s += ['<thead><tr><th style="text-align: left;">Field</th><th>Value</th></tr></thead>']
298 s += ['<tbody>']
299 s += [f'<tr><td style="text-align: left;">Space group</td><td>{self.space_group}</td></tr>']
300 for sl in self.get_sublattices(self.primitive_structure).active_sublattices:
301 s += [f'<tr><td style="text-align: left;">Sublattice {sl.symbol}</td>'
302 f'<td>{sl.chemical_symbols}</td></tr>']
303 s += [f'<tr><td style="text-align: left;">Cutoffs</td><td>{self.cutoffs}</td></tr>']
304 s += ['<tr><td style="text-align: left;">Total number of parameters</td>'
305 f'<td>{len(self)}</td></tr>']
306 for k, n in self.number_of_orbits_by_order.items():
307 s += [f'<tr><td style="text-align: left;">Number of parameters of order {k}</td>'
308 f'<td>{n}</td></tr>']
309 for key, value in sorted(self._config.items()):
310 s += [f'<tr><td style="text-align: left;">{key}</td><td>{value}</td></tr>']
311 s += ['</tbody>']
312 s += ['</table>']
313 return ''.join(s)
315 def __repr__(self) -> str:
316 """ Representation. """
317 s = type(self).__name__ + '('
318 s += f'structure={self.primitive_structure.__repr__()}'
319 s += f', cutoffs={self._cutoffs.__repr__()}'
320 s += f', chemical_symbols={self._input_chemical_symbols.__repr__()}'
321 s += f', position_tolerance={self._config["position_tolerance"]}'
322 s += ')'
323 return s
325 @property
326 def symprec(self) -> float:
327 """ Tolerance imposed when analyzing the symmetry using spglib. """
328 return self._config['symprec']
330 @property
331 def position_tolerance(self) -> float:
332 """ Tolerance applied when comparing positions in Cartesian coordinates. """
333 return self._config['position_tolerance']
335 @property
336 def fractional_position_tolerance(self) -> float:
337 """ Tolerance applied when comparing positions in fractional coordinates. """
338 return self._config['fractional_position_tolerance']
340 @property
341 def space_group(self) -> str:
342 """ Space group of the primitive structure in international notion (via spglib). """
343 structure_as_tuple = ase_atoms_to_spglib_cell(self.primitive_structure)
344 return spglib.get_spacegroup(structure_as_tuple, symprec=self._config['symprec'])
346 @property
347 def as_list(self) -> List[dict]:
348 """Representation of cluster space as list with information regarding
349 order, radius, multiplicity etc.
350 """
351 data = []
352 zerolet = dict(
353 index=0,
354 order=0,
355 radius=0,
356 multiplicity=1,
357 orbit_index=-1,
358 multicomponent_vector='.',
359 sublattices='.',
360 )
361 data.append(zerolet)
363 sublattices = self.get_sublattices(self.primitive_structure)
364 index = 0
365 for orbit_index in range(len(self.orbit_list)):
366 orbit = self.orbit_list.get_orbit(orbit_index)
367 representative_cluster = orbit.representative_cluster
368 orbit_sublattices = [
369 sublattices[sublattices.get_sublattice_index_from_site_index(ls.index)].symbol
370 for ls in representative_cluster.lattice_sites]
371 for cv_element in orbit.cluster_vector_elements:
372 index += 1
373 data.append(dict(
374 index=index,
375 order=representative_cluster.order,
376 radius=representative_cluster.radius,
377 multiplicity=cv_element['multiplicity'],
378 orbit_index=orbit_index,
379 multicomponent_vector=cv_element['multicomponent_vector'],
380 sublattices=orbit_sublattices
381 ))
382 return data
384 def to_dataframe(self) -> DataFrame:
385 """ Returns a representation of the cluster space as a DataFrame. """
386 df = DataFrame.from_dict(self.as_list)
387 del df['index']
388 return df
390 @property
391 def number_of_orbits_by_order(self) -> dict:
392 """ Number of orbits by order in the form of a dictionary
393 where keys and values represent order and number of orbits,
394 respectively.
395 """
396 count_orbits: Dict[int, int] = {}
397 for orbit in self.as_list:
398 k = orbit['order']
399 count_orbits[k] = count_orbits.get(k, 0) + 1
400 return dict(sorted(count_orbits.items()))
402 def get_cluster_vector(self, structure: Atoms) -> np.ndarray:
403 """
404 Returns the cluster vector for a structure.
406 Parameters
407 ----------
408 structure
409 Atomic configuration.
410 """
411 if not isinstance(structure, Atoms): 411 ↛ 412line 411 didn't jump to line 412, because the condition on line 411 was never true
412 raise TypeError('Input structure must be an ASE Atoms object')
414 try:
415 cv = _ClusterSpace.get_cluster_vector(
416 self,
417 structure=Structure.from_atoms(structure),
418 fractional_position_tolerance=self.fractional_position_tolerance)
419 except Exception as e:
420 self.assert_structure_compatibility(structure)
421 raise Exception(str(e))
422 return cv
424 def get_coordinates_of_representative_cluster(self, orbit_index: int) -> List[Tuple[float]]:
425 """
426 Returns the positions of the sites in the representative cluster of the selected orbit.
428 Parameters
429 ----------
430 orbit_index
431 Index of the orbit for which to return the positions of the sites.
432 """
433 # Raise exception if chosen orbit index not in current list of orbit indices
434 if orbit_index not in range(len(self._orbit_list)):
435 raise ValueError('The input orbit index is not in the list of possible values.')
436 return self._orbit_list.get_orbit(orbit_index).representative_cluster.positions
438 def _remove_orbits(self, indices: List[int]) -> None:
439 """
440 Removes orbits.
442 Parameters
443 ----------
444 indices
445 Indices to all orbits to be removed.
446 """
447 size_before = len(self._orbit_list)
449 # Since we remove orbits, orbit indices will change,
450 # so we run over the orbits in reverse order.
451 for ind in reversed(sorted(indices)):
452 self._orbit_list.remove_orbit(ind)
454 size_after = len(self._orbit_list)
455 assert size_before - len(indices) == size_after
457 def _prune_orbit_list(self, indices: List[int]) -> None:
458 """
459 Prunes the internal orbit list and maintains the history.
461 Parameters
462 ----------
463 indices
464 Indices to all orbits to be removed.
465 """
466 self._remove_orbits(indices)
467 self._pruning_history.append(('prune', indices))
469 @property
470 def primitive_structure(self) -> Atoms:
471 """ Primitive structure on which cluster space is based. """
472 structure = self._get_primitive_structure().to_atoms()
473 # Decorate with the "real" symbols (instead of H, He, Li etc)
474 for atom, symbols in zip(structure, self._primitive_chemical_symbols):
475 atom.symbol = min(symbols)
476 return structure
478 @property
479 def chemical_symbols(self) -> List[List[str]]:
480 """ Species identified by their chemical symbols. """
481 return self._primitive_chemical_symbols.copy()
483 @property
484 def cutoffs(self) -> List[float]:
485 """
486 Cutoffs for different n-body clusters. The cutoff radius (in
487 Ångstroms) defines the largest interatomic distance in a
488 cluster.
489 """
490 return self._cutoffs
492 @property
493 def orbit_list(self):
494 """ Orbit list that defines the cluster in the cluster space. """
495 return self._orbit_list
497 def get_possible_orbit_occupations(self, orbit_index: int) -> List[List[str]]:
498 """ Returns possible occupations of the orbit.
500 Parameters
501 ----------
502 orbit_index
503 Index of orbit of interest.
504 """
505 orbit = self.orbit_list.orbits[orbit_index]
506 indices = [ls.index for ls in orbit.representative_cluster.lattice_sites]
507 allowed_species = [self.chemical_symbols[index] for index in indices]
508 return list(itertools.product(*allowed_species))
510 def get_sublattices(self, structure: Atoms) -> Sublattices:
511 """ Returns the sublattices of the input structure.
513 Parameters
514 ----------
515 structure
516 Atomic structure the sublattices are based on.
517 """
518 sl = Sublattices(self.chemical_symbols,
519 self.primitive_structure,
520 structure,
521 fractional_position_tolerance=self.fractional_position_tolerance)
522 return sl
524 def assert_structure_compatibility(self, structure: Atoms, vol_tol: float = 1e-5) -> None:
525 """ Raises error if structure is not compatible with this cluster space.
527 Parameters
528 ----------
529 structure
530 Structure to check for compatibility with cluster space.
531 vol_tol
532 Tolerance imposed when comparing volumes.
533 """
534 # check volume
535 vol1 = self.primitive_structure.get_volume() / len(self.primitive_structure)
536 vol2 = structure.get_volume() / len(structure)
537 if abs(vol1 - vol2) > vol_tol:
538 raise ValueError(f'Volume per atom of structure ({vol1}) does not match the volume of'
539 f' the primitive structure ({vol2}; vol_tol= {vol_tol}).')
541 # check occupations
542 sublattices = self.get_sublattices(structure)
543 sublattices.assert_occupation_is_allowed(structure.get_chemical_symbols())
545 # check pbc
546 if not all(structure.pbc):
547 raise ValueError('Input structure must be periodic.')
549 def merge_orbits(self,
550 equivalent_orbits: Dict[int, List[int]],
551 ignore_permutations: bool = False) -> None:
552 """ Combines several orbits into one. This allows one to make custom
553 cluster spaces by manually declaring the clusters in two or more
554 orbits to be equivalent. This is a powerful approach for simplifying
555 the cluster spaces of low-dimensional structures such as
556 surfaces or nanoparticles.
558 The procedure works in principle for any number of components. Note,
559 however, that in the case of more than two components the outcome of
560 the merging procedure inherits the treatment of the multi-component
561 vectors of the orbit chosen as the representative one.
563 Parameters
564 ----------
565 equivalent_orbits
566 The keys of this dictionary denote the indices of the orbit into
567 which to merge. The values are the indices of the orbits that are
568 supposed to be merged into the orbit denoted by the key.
569 ignore_permutations
570 If ``True`` orbits will be merged even if their multi-component
571 vectors and/or site permutations differ. While the object will
572 still be functional, the cluster space may not be properly spanned
573 by the resulting cluster vectors.
575 Note
576 ----
577 The orbit index should not be confused with the index shown when
578 printing the cluster space.
580 Examples
581 --------
582 The following snippet illustrates the use of this method to create a
583 cluster space for a (111) FCC surface, in which only the singlets for
584 the first and second layer are distinct as well as the in-plane pair
585 interaction in the topmost layer. All other singlets and pairs are
586 respectively merged into one orbit. After merging there aree only 3
587 singlets and 2 pairs left with correspondingly higher multiplicities.
589 >>> from icet import ClusterSpace
590 >>> from ase.build import fcc111
591 >>>
592 >>> # Create primitive surface unit cell
593 >>> structure = fcc111('Au', size=(1, 1, 8), a=4.1, vacuum=10, periodic=True)
594 >>>
595 >>> # Set up initial cluster space
596 >>> cs = ClusterSpace(structure=structure, cutoffs=[3.8], chemical_symbols=['Au', 'Ag'])
597 >>>
598 >>> # At this point, one can inspect the orbits in the cluster space by printing the
599 >>> # ClusterSpace object and accessing the individial orbits.
600 >>> # There will be 4 singlets and 8 pairs.
601 >>>
602 >>> # Merge singlets for the third and fourth layers as well as all pairs except for
603 >>> # the one corresponding to the in-plane interaction in the topmost surface
604 >>> # layer.
605 >>> cs.merge_orbits({2: [3], 4: [6, 7, 8, 9, 10, 11]})
606 """
608 self._pruning_history.append(('merge', equivalent_orbits))
609 orbits_to_delete = []
610 for k1, orbit_indices in equivalent_orbits.items():
611 orbit1 = self.orbit_list.get_orbit(k1)
613 for k2 in orbit_indices:
615 # sanity checks
616 if k1 == k2:
617 raise ValueError(f'Cannot merge orbit {k1} with itself.')
618 if k2 in orbits_to_delete:
619 raise ValueError(f'Orbit {k2} cannot be merged into orbit {k1}'
620 ' since it was already merged with another orbit.')
621 orbit2 = self.orbit_list.get_orbit(k2)
622 if orbit1.order != orbit2.order:
623 raise ValueError(f'The order of orbit {k1} ({orbit1.order}) does not'
624 f' match the order of orbit {k2} ({orbit2.order}).')
626 if not ignore_permutations:
627 # compare site permutations
628 permutations1 = [el['site_permutations']
629 for el in orbit1.cluster_vector_elements]
630 permutations2 = [el['site_permutations']
631 for el in orbit2.cluster_vector_elements]
632 for vec_group1, vec_group2 in zip(permutations1, permutations2):
633 if len(vec_group1) != len(vec_group2) or \
634 not np.allclose(np.array(vec_group1), np.array(vec_group2)):
635 raise ValueError(f'Orbit {k1} and orbit {k2} have different '
636 'site permutations.')
638 # compare multi-component vectors (maybe this is redundant because
639 # site permutations always differ if multi-component vectors differ?)
640 mc_vectors1 = [el['multicomponent_vector']
641 for el in orbit1.cluster_vector_elements]
642 mc_vectors2 = [el['multicomponent_vector']
643 for el in orbit2.cluster_vector_elements]
644 if not all(np.allclose(vec1, vec2) 644 ↛ 646line 644 didn't jump to line 646, because the condition on line 644 was never true
645 for vec1, vec2 in zip(mc_vectors1, mc_vectors2)):
646 raise ValueError(f'Orbit {k1} and orbit {k2} have different '
647 'multi-component vectors.')
649 # merge
650 self._merge_orbit(k1, k2)
651 orbits_to_delete.append(k2)
653 # update merge/prune history
654 self._remove_orbits(orbits_to_delete)
656 def is_supercell_self_interacting(self, structure: Atoms) -> bool:
657 """
658 Checks whether a structure has self-interactions via periodic
659 boundary conditions.
660 Returns ``True`` if the structure contains self-interactions via periodic
661 boundary conditions, otherwise ``False``.
663 Parameters
664 ----------
665 structure
666 Structure to be tested.
667 """
668 ol = self.orbit_list.get_supercell_orbit_list(
669 structure=structure,
670 fractional_position_tolerance=self.fractional_position_tolerance)
671 orbit_indices = set()
672 for orbit in ol.orbits:
673 for cluster in orbit.clusters:
674 indices = tuple(sorted([site.index for site in cluster.lattice_sites]))
675 if indices in orbit_indices:
676 return True
677 else:
678 orbit_indices.add(indices)
679 return False
681 def write(self, filename: str) -> None:
682 """
683 Saves cluster space to a file.
685 Parameters
686 ---------
687 filename
688 Name of file to which to write.
689 """
691 with tarfile.open(name=filename, mode='w') as tar_file:
693 # write items
694 items = dict(cutoffs=self._cutoffs,
695 chemical_symbols=self._input_chemical_symbols,
696 pruning_history=self._pruning_history,
697 symprec=self.symprec,
698 position_tolerance=self.position_tolerance)
699 temp_file = tempfile.TemporaryFile()
700 pickle.dump(items, temp_file)
701 temp_file.seek(0)
702 tar_info = tar_file.gettarinfo(arcname='items', fileobj=temp_file)
703 tar_file.addfile(tar_info, temp_file)
704 temp_file.close()
706 # write structure
707 temp_file = tempfile.NamedTemporaryFile(delete=False)
708 temp_file.close()
709 ase_write(temp_file.name, self._input_structure, format='json')
710 with open(temp_file.name, 'rb') as tt:
711 tar_info = tar_file.gettarinfo(arcname='atoms', fileobj=tt)
712 tar_file.addfile(tar_info, tt)
713 os.remove(temp_file.name)
715 @staticmethod
716 def read(filename: str):
717 """
718 Reads cluster space from file and returns :attr:`ClusterSpace` object.
720 Parameters
721 ---------
722 filename
723 Name of file from which to read cluster space.
724 """
725 if isinstance(filename, str):
726 tar_file = tarfile.open(mode='r', name=filename)
727 else:
728 tar_file = tarfile.open(mode='r', fileobj=filename)
730 # read items
731 items = pickle.load(tar_file.extractfile('items'))
733 # read structure
734 temp_file = tempfile.NamedTemporaryFile(delete=False)
735 temp_file.write(tar_file.extractfile('atoms').read())
736 temp_file.close()
737 structure = ase_read(temp_file.name, format='json')
738 os.remove(temp_file.name)
740 tar_file.close()
742 # ensure backward compatibility
743 if 'symprec' not in items: # pragma: no cover
744 items['symprec'] = 1e-5
745 if 'position_tolerance' not in items: # pragma: no cover
746 items['position_tolerance'] = items['symprec']
748 cs = ClusterSpace(structure=structure,
749 cutoffs=items['cutoffs'],
750 chemical_symbols=items['chemical_symbols'],
751 symprec=items['symprec'],
752 position_tolerance=items['position_tolerance'])
753 if len(items['pruning_history']) > 0:
754 if isinstance(items['pruning_history'][0], tuple): 754 ↛ 763line 754 didn't jump to line 763, because the condition on line 754 was never false
755 for key, value in items['pruning_history']:
756 if key == 'prune':
757 cs._prune_orbit_list(value)
758 elif key == 'merge': 758 ↛ 755line 758 didn't jump to line 755, because the condition on line 758 was never false
759 # It is safe to ignore permutations here because otherwise
760 # the orbits could not have been merged in the first place.
761 cs.merge_orbits(value, ignore_permutations=True)
762 else: # for backwards compatibility
763 for value in items['pruning_history']:
764 cs._prune_orbit_list(value)
766 return cs
768 def copy(self):
769 """ Returns copy of :class:`ClusterSpace` instance. """
770 cs_copy = ClusterSpace(structure=self._input_structure,
771 cutoffs=self.cutoffs,
772 chemical_symbols=self._input_chemical_symbols,
773 symprec=self.symprec,
774 position_tolerance=self.position_tolerance)
776 for key, value in self._pruning_history:
777 if key == 'prune':
778 cs_copy._prune_orbit_list(value)
779 elif key == 'merge': 779 ↛ 776line 779 didn't jump to line 776, because the condition on line 779 was never false
780 # It is safe to ignore permutations here because otherwise
781 # the orbits could not have been merged in the first place.
782 cs_copy.merge_orbits(value, ignore_permutations=True)
783 return cs_copy