Coverage for hiphive/core/orbits.py: 97%
147 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-08-01 17:04 +0000
1"""
2Contains the Orbit class which hold onformation about equivalent clusters.
3"""
5import tarfile
6import pickle
7from typing import Union, BinaryIO, TextIO
8import numpy as np
9from .utilities import BiMap
10from .atoms import Atom
11from . import atoms as atoms_module
12from ..input_output.logging_tools import logger
13from ..input_output.read_write_files import add_items_to_tarfile_pickle, read_items_pickle
16logger = logger.getChild('orbits')
19class Orbit:
20 """
21 This class serves as a container for storing data pertaining to an orbit.
23 Attributes
24 ----------
25 orientation_families : list[OrientationFamily]
26 Orientation families of the orbit.
27 eigensymmetries : list[tuple]
28 Each eigensymmetry corresponds to a pair where the first index
29 is the symmetry and the second is the permutation.
30 eigentensors : list[np.ndarray]
31 Decomposition of the force constant into symmetry elements.
32 """
34 # TODO: Make properties of the parameters
35 def __init__(self):
36 self.orientation_families = []
37 self.eigensymmetries = []
38 self.eigentensors = []
39 self.radius = None
40 self.order = None
41 self.maximum_distance = None
43 @property
44 def prototype_index(self) -> int:
45 """ Index of cluster that serves as prototype for this orbit.
46 In the code the first symmetry is always the identity so the first
47 orientation family should always correspond to the prototype.
48 """
49 return self.orientation_families[0].cluster_indices[0]
51 def write(self, f: Union[str, BinaryIO, TextIO]) -> None:
52 """Write an Orbit instance to a file.
54 Parameters
55 ----------
56 f
57 Name of input file (`str`) or stream to write to (file object).
58 """
59 tar_file = tarfile.open(mode='w', fileobj=f)
61 # add all available attributes
62 add_items_to_tarfile_pickle(tar_file, self.__dict__, 'attributes')
64 @staticmethod
65 def read(f: Union[str, BinaryIO, TextIO]):
66 """Load an :class:`Orbit` from file.
68 Parameters
69 ----------
70 f
71 Name of input file (`str`) or stream to load from (file object).
72 """
74 orb = Orbit()
75 tar_file = tarfile.open(mode='r', fileobj=f)
77 # read attributes pickle
78 attributes = read_items_pickle(tar_file, 'attributes')
79 for name, value in attributes.items():
80 orb.__setattr__(name, value)
81 return orb
84class OrientationFamily:
85 """A container for storing information for a "family of orientations".
87 An orbit contains many clusters. Some of the clusters can be translated
88 onto each other and other must first be rotated. A set of clusters in the
89 orbit which can all be translated onto each other are oriented in the same
90 way and belong to the same orientation family. The family is characterized
91 by the symmetry (rotation) which relates it to the prototype structure of
92 the orbit.
93 Since the clusters are generally stored sorted the permutation information
94 must also be stored.
96 Parameters
97 ----------
98 symmetry_index : int
99 The index of the symmetry corresponding to spglibs symmetry.
101 Attributes
102 ----------
103 symmetry_index : int
104 The index of the symmetry corresponding to the symmetry according to spglib.
105 cluster_indices : list[int]
106 The indices of the clusters belonging to this family.
107 permutation_indices : list[int]
108 The indices of the permutation vector.
109 """
111 def __init__(self, symmetry_index=None):
112 self.symmetry_index = symmetry_index
113 self.cluster_indices = []
114 self.permutation_indices = []
116 def write(self, f: Union[str, BinaryIO, TextIO]):
117 """ Write the object to file.
119 Parameters
120 ----------
121 f
122 Name of input file (`str`) or stream to write to (file object).
123 """
124 pickle.dump(self, f)
126 @staticmethod
127 def read(f: Union[str, BinaryIO, TextIO]):
128 """Load an :class:`OrientationFamily` object from a pickle file.
130 Parameters
131 ----------
132 f
133 Name of input file (`str`) or stream to load from (file object).
135 Returns
136 -------
137 OrientationFamily
138 """
139 return pickle.load(f)
142# This is the interface accessible for cluster_space
143def get_orbits(
144 cluster_list: BiMap,
145 atom_list: BiMap,
146 rotation_matrices: list[np.ndarray],
147 translation_vectors: list[np.ndarray],
148 permutations: list,
149 prim: atoms_module.Atoms,
150 symprec: float,
151) -> list[Orbit]:
152 """Generate a list of the orbits for the clusters in a supercell
153 configuration.
155 This method requires as input a list of the clusters in a supercell
156 configuration as well as a set of symmetry operations (rotations and
157 translations). From this information it will generate a list of the orbits,
158 i.e. the set of symmetry inequivalent clusters each associated with its
159 respective set of equivalent clusters.
161 Parameters
162 ----------
163 cluster_list
164 List of clusters.
165 atom_list
166 List of atoms in a supercell.
167 rotation_matrices
168 Rotational symmetries as 3x3 matrices to be imposed (e.g., from spglib).
169 translation_vectors
170 Translational symmetries as 3 vectors to be imposed (e.g., from spglib).
171 permutations
172 Look-up table for permutations.
173 prim
174 Primitive structure.
176 Returns
177 -------
178 Orbits associated with the list of input clusters.
180 """
182 logger.debug('Preparing input...')
183 atoms = prepare_atoms(atom_list)
184 clusters = prepare_clusters(cluster_list)
185 rotations = prepare_rotations(rotation_matrices)
186 translations = prepare_translations(translation_vectors)
187 permutations = prepare_permutations(permutations)
188 cell = prim.cell
189 basis = prim.basis
191 logger.debug('Creating permutation map...')
192 permutation_map, extended_atoms = \
193 get_permutation_map(atoms, rotations, translations, basis, symprec)
195 logger.debug('Creating orbits...')
196 orbits = _get_orbits(permutation_map, extended_atoms, clusters, basis,
197 cell, rotations, permutations)
199 return orbits
202# All prepares are in case we changes some interface stuff
203def prepare_permutations(permutations):
204 perms = BiMap()
205 for p in permutations:
206 perms.append(tuple(p))
207 return perms
210def prepare_atoms(atom_list):
211 atoms = BiMap()
212 for atom in atom_list:
213 atoms.append(atom)
214 return atoms
217def prepare_clusters(cluster_list):
218 clusters = BiMap()
219 for cluster in cluster_list:
220 clusters.append(tuple(cluster))
221 return clusters
224def prepare_rotations(rotation_matrices):
225 return rotation_matrices
228def prepare_translations(translation_vectors):
229 return translation_vectors
232def get_permutation_map(atoms, rotations, translations, basis, symprec):
234 extended_atoms = atoms.copy()
236 permutation_map = np.zeros((len(atoms), len(rotations)), dtype=int)
238 scaled_positions = [atom.spos(basis) for atom in extended_atoms]
240 for sym_index, (R, T) in enumerate(zip(rotations, translations)):
241 for atom_index, spos in enumerate(scaled_positions):
243 new_spos = np.dot(R, spos) + T
244 new_atom = Atom.spos_to_atom(new_spos, basis, symprec)
246 if new_atom not in extended_atoms:
247 extended_atoms.append(new_atom)
249 mapped_atom_index = extended_atoms.index(new_atom)
250 permutation_map[atom_index, sym_index] = mapped_atom_index
252 return permutation_map, extended_atoms
255# The internal function doing the outer loop over orbits
256def _get_orbits(permutation_map, extended_atoms, clusters,
257 basis, cell,
258 rotations, permutations):
259 cluster_is_found = [False] * len(clusters)
260 orbits = []
261 for cluster_index, cluster in enumerate(clusters):
262 if cluster_is_found[cluster_index]:
263 continue
265 orbit = Orbit()
267 cluster_atoms = [extended_atoms[i] for i in cluster]
269 positions = [atom.pos(basis, cell) for atom in cluster_atoms]
270 orbit.radius = get_geometrical_radius(positions)
271 orbit.maximum_distance = get_maximum_distance(positions)
272 orbit.order = len(cluster)
274 populate_orbit(orbit, permutations, clusters,
275 cluster,
276 permutation_map, extended_atoms,
277 cluster_is_found)
278 orbits.append(orbit)
280# bar.tick()
281 return orbits
284# Takes a cluster and generates all equivalent, translated clusters
285def generate_translated_clusters(cluster, extended_atoms):
286 transformed_cluster_atoms = [extended_atoms[i] for i in cluster]
287 tested_offsets = set()
288 for atom in transformed_cluster_atoms:
289 offset = atom.offset
290 if offset in tested_offsets:
291 continue
292 else:
293 tested_offsets.add(offset)
294 translated_cluster = []
295 for atom in transformed_cluster_atoms:
296 new_offset = np.subtract(atom.offset, offset)
297 new_atom = Atom(atom.site, new_offset)
298 translated_cluster.append(extended_atoms.index(new_atom))
299 yield tuple(translated_cluster)
302# Here is the actual categorization
303def populate_orbit(orbit, permutations, clusters,
304 cluster,
305 permutation_map, extended_atoms,
306 cluster_is_found):
307 for sym_index in range(permutation_map.shape[1]):
309 of = OrientationFamily(sym_index)
311 transformed_cluster = permutation_map[cluster, sym_index]
313 for translated_cluster in generate_translated_clusters(
314 transformed_cluster, extended_atoms):
316 argsort = tuple(np.argsort(translated_cluster))
317 perm_index = permutations.index(argsort)
319 translated_cluster = tuple(sorted(translated_cluster))
320 translated_cluster_index = clusters.index(translated_cluster)
322 if cluster == translated_cluster:
323 if (sym_index, perm_index) not in orbit.eigensymmetries: 323 ↛ 326line 323 didn't jump to line 326 because the condition on line 323 was always true
324 orbit.eigensymmetries.append((sym_index, perm_index))
326 if not cluster_is_found[translated_cluster_index]:
327 cluster_is_found[translated_cluster_index] = True
328 of.cluster_indices.append(translated_cluster_index)
329 of.permutation_indices.append(perm_index)
331 if len(of.cluster_indices) > 0:
332 orbit.orientation_families.append(of)
334 return orbit
337def get_geometrical_radius(positions: list[np.ndarray]) -> float:
338 """Compute the geometrical size of a 3-dimensional point cloud. The
339 geometrical size is defined as the average distance to the geometric
340 center.
342 Parameters
343 ----------
344 positions
345 Positions of points in cloud.
347 Returns
348 -------
349 Geometric size of point cloud.
350 """
351 geometric_center = np.mean(positions, axis=0)
352 return np.mean(np.sqrt(np.sum((positions - geometric_center)**2, axis=1)))
355def get_maximum_distance(positions: list[np.ndarray]) -> float:
356 """Compute the maximum distance between any two points in a 3-dimensional
357 point cloud. This is equivalent to the "size" criterion used when imposing
358 a certain (pair) cutoff criterion during construction of a set of clusters.
360 Parameters
361 ----------
362 positions
363 Positions of points in cloud.
365 Returns
366 -------
367 Maximum distance betwee any two points.
368 """
369 if len(positions) == 0: 369 ↛ 370line 369 didn't jump to line 370 because the condition on line 369 was never true
370 return 0.0
371 max_distance = 0.0
372 for pt1 in positions[:-1]:
373 for pt2 in positions[1:]:
374 max_distance = max(max_distance, np.linalg.norm(pt1 - pt2))
375 return max_distance