Coverage for hiphive/cluster_space.py: 98%
Shortcuts on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
Shortcuts on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Contains the ClusterSpace object central to hiPhive
3"""
5import numpy as np
6import tarfile
8from ase.data import chemical_symbols
9from collections import OrderedDict
10from copy import deepcopy
12from .core.cluster_space_builder import build_cluster_space
13from .core.atoms import Atoms
14from .core.orbits import Orbit
15from .cluster_space_data import ClusterSpaceData
16from .input_output.logging_tools import logger
17from .input_output.read_write_files import (add_items_to_tarfile_pickle,
18 add_items_to_tarfile_custom,
19 add_list_to_tarfile_custom,
20 read_items_pickle,
21 read_list_custom)
22from .cutoffs import Cutoffs, CutoffMaximumBody, BaseClusterFilter
24from .config import Config
25logger = logger.getChild('ClusterSpace')
28class ClusterSpace:
29 """Primitive object for handling clusters and force constants of a structure.
31 Parameters
32 ----------
33 prototype_structure : ase.Atoms
34 prototype structure; spglib will be used to find a suitable cell based
35 on this structure unless the cell is already a primitive cell.
36 cutoffs : list or Cutoffs
37 cutoff radii for different orders starting with second order
38 cluster_filter : ClusterFilter
39 accepts a subclass of hiphive.filters.BaseClusterFilter to further
40 control which orbits to include.
41 config : Config object
42 a configuration object that holds information on how the cluster space
43 should be built, e.g., values for tolerances and specifications
44 regarding the handling of acoustic sum rules; if ``config`` is
45 not given then the keyword arguments that follow below can be
46 used for configuration.
47 acoustic_sum_rules : bool
48 If True the aucostic sum rules will be enforced by constraining the
49 parameters.
50 symprec : float
51 numerical precision that will be used for analyzing the symmetry (this
52 parameter will be forwarded to
53 `spglib <https://atztogo.github.io/spglib/>`_)
54 length_scale : float
55 This will be used as a normalization constant for the eigentensors
57 Examples
58 --------
60 To instantiate a :class:`ClusterSpace` object one has to specify a
61 prototype structure and cutoff radii for each cluster order that
62 should be included. For example the following snippet will set up
63 a :class:`ClusterSpace` object for a body-centered-cubic (BCC)
64 structure including second order terms up to a distance of 5 A and
65 third order terms up to a distance of 4 A.
67 >>> from ase.build import bulk
68 >>> from hiphive import ClusterSpace
69 >>> prim = bulk('W')
70 >>> cs = ClusterSpace(prim, [5.0, 4.0])
72 """
73 # TODO: This class probably needs some more documentation
75 def __init__(self, prototype_structure, cutoffs, config=None,
76 cluster_filter=None, **kwargs):
78 if not all(prototype_structure.pbc):
79 raise ValueError('prototype_structure must have pbc.')
81 if isinstance(cutoffs, Cutoffs):
82 self._cutoffs = cutoffs
83 elif isinstance(cutoffs, list):
84 self._cutoffs = CutoffMaximumBody(cutoffs, len(cutoffs) + 1)
85 else:
86 raise TypeError('cutoffs must be a list or a Cutoffs object')
88 if config is None:
89 config = Config(**kwargs)
90 else:
91 if not isinstance(config, Config): 91 ↛ 92line 91 didn't jump to line 92, because the condition on line 91 was never true
92 raise TypeError('config kw must be of type {}'.format(Config))
93 if len(kwargs) > 0:
94 raise ValueError('use either Config or kwargs, not both')
95 self._config = config
97 if cluster_filter is None:
98 self._cluster_filter = BaseClusterFilter()
99 else:
100 self._cluster_filter = cluster_filter
102 self._atom_list = None
103 self._cluster_list = None
104 self._symmetry_dataset = None
105 self._permutations = None
106 self._prim = None
107 self._orbits = None
109 self._constraint_vectors = None
110 # TODO: How to handle the constraint matrices? Should they even be
111 # stored?
112 self._constraint_matrices = None
113 # Is this the best way or should the prim be instantiated separately?
115 build_cluster_space(self, prototype_structure)
116 self.summary = ClusterSpaceData(self)
118 # TODO: Should everything here be properties? deepcopy/ref etc.?
119 @property
120 def n_dofs(self):
121 """int : number of free parameters in the model
123 If the sum rules are not enforced the number of DOFs is the same as
124 the total number of eigentensors in all orbits.
125 """
126 return self._get_n_dofs()
128 @property
129 def cutoffs(self):
130 """ Cutoffs : cutoffs used for constructing the cluster space """
131 return deepcopy(self._cutoffs)
133 @property
134 def symprec(self):
135 """ float : symprec value used when constructing the cluster space """
136 return self._config['symprec']
138 @property
139 def acoustic_sum_rules(self):
140 """ bool : True if acoustic sum rules are enforced """
141 return self._config['acoustic_sum_rules']
143 @property
144 def length_scale(self):
145 """ float : normalization constant of the force constants """
146 return self._config['length_scale']
148 @property
149 def primitive_structure(self):
150 """ ase.Atoms : structure of the lattice """
151 return self._prim.copy()
153 @property
154 def spacegroup(self):
155 """ str : space group of the lattice structure obtained from spglib """
156 return self._symmetry_dataset['international'] + ' ({})'.format(
157 self._symmetry_dataset['number'])
159 @property
160 def wyckoff_sites(self):
161 """ list : wyckoff sites in the primitive cell """
162 return self._symmetry_dataset['equivalent_atoms']
164 @property
165 def rotation_matrices(self):
166 """ list(numpy.ndarray) : symmetry elements (`3x3` matrices)
167 representing rotations """
168 return self._symmetry_dataset['rotations'].copy()
170 @property
171 def translation_vectors(self):
172 """ list(numpy.ndarray) : symmetry elements representing
173 translations """
174 # TODO: bug incoming!
175 return (self._symmetry_dataset['translations'] % 1).copy()
177 @property
178 def permutations(self):
179 """ list(numpy.ndarray) : lookup for permutation references """
180 return deepcopy(self._permutations)
182 @property
183 def atom_list(self):
184 """ BiMap : atoms inside the cutoff relative to the of the center cell
185 """
186 return self._atom_list
188 @property
189 def cluster_list(self):
190 """ BiMap : clusters possible within the cutoff """
191 return self._cluster_list
193 @property
194 def orbits(self): # TODO: add __getitem__ method
195 """ list(Orbit) : orbits associated with the lattice structure. """
196 return self._orbits
198 @property
199 def orbit_data(self):
200 """ list(dict) : detailed information for each orbit, e.g., cluster
201 radius and atom types.
202 """
203 data = []
204 p = 0
205 for orbit_index, orbit in enumerate(self.orbits):
206 d = {}
207 d['index'] = orbit_index
208 d['order'] = orbit.order
209 d['radius'] = orbit.radius
210 d['maximum_distance'] = orbit.maximum_distance
211 d['n_clusters'] = len(orbit.orientation_families)
212 d['eigentensors'] = orbit.eigentensors
213 d['n_parameters'] = len(d['eigentensors'])
215 types, wyckoff_sites = [], []
216 for atom_index in self.cluster_list[orbit.prototype_index]:
217 atom = self.atom_list[atom_index]
218 types.append(self.primitive_structure.numbers[atom.site])
219 wyckoff_sites.append(self.wyckoff_sites[atom.site])
220 d['prototype_cluster'] = self.cluster_list[orbit.prototype_index]
221 d['prototype_atom_types'] = tuple(types)
222 d['prototype_wyckoff_sites'] = tuple(wyckoff_sites)
224 d['geometrical_order'] = len(set(d['prototype_cluster']))
225 d['parameter_indices'] = np.arange(p, p + len(orbit.eigentensors))
227 p += len(orbit.eigentensors)
228 data.append(d)
230 return data
232 def get_parameter_indices(self, order):
233 """
234 Returns a list of the parameter indices associated with the requested
235 order.
237 Parameters
238 ----------
239 order : int
240 order for which to return the parameter indices
242 Returns
243 -------
244 list(int)
245 list of parameter indices associated with the requested order
247 Raises
248 ------
249 ValueError
250 if the order is not included in the cluster space
251 """
252 order = int(order)
253 if order not in self.cutoffs.orders: 253 ↛ 254line 253 didn't jump to line 254, because the condition on line 253 was never true
254 raise ValueError('Order must be in {}'.format(self.cutoffs.orders))
255 min_param = 0
256 max_param = 0
257 for orbit in self.orbits:
258 if orbit.order < order:
259 min_param += len(orbit.eigentensors)
260 max_param = min_param
261 elif orbit.order == order:
262 max_param += len(orbit.eigentensors)
263 else:
264 break
265 rows, cols = self._cvs.nonzero()
266 parameters = set()
267 for r, c in zip(rows, cols):
268 if min_param <= r < max_param:
269 parameters.add(c)
270 for r, c in zip(rows, cols):
271 if c in parameters:
272 assert min_param <= r < max_param, 'The internals are broken!'
274 return sorted(parameters)
276 def get_n_dofs_by_order(self, order):
277 """ Returns number of degrees of freedom for the given order.
279 Parameters
280 ----------
281 order : int
282 order for which to return the number of dofs
284 Returns
285 -------
286 int
287 number of degrees of freedom
288 """
289 return len(self.get_parameter_indices(order=order))
291 def _get_n_dofs(self):
292 """ Returns the number of degrees of freedom. """
293 return self._cvs.shape[1]
295 def _map_parameters(self, parameters):
296 """ Maps irreducible parameters to the real parameters associated with
297 the eigentensors.
298 """
299 if len(parameters) != self.n_dofs: 299 ↛ 300line 299 didn't jump to line 300, because the condition on line 299 was never true
300 raise ValueError('Invalid number of parameters, please provide {} '
301 'parameters'.format(self.n_dofs))
302 return self._cvs.dot(parameters)
304 def print_tables(self):
305 """ Prints information concerning the underlying cluster space to stdout,including,
306 e.g., the number of cluster, orbits, and parameters by order and number of bodies. """
307 self.summary.print_tables()
309 def print_orbits(self):
310 """ Prints a list of all orbits. """
311 orbits = self.orbit_data
313 def str_orbit(index, orbit):
314 elements = ' '.join(chemical_symbols[n] for n in
315 orbit['prototype_atom_types'])
316 fields = OrderedDict([
317 ('index', '{:^3}'.format(index)),
318 ('order', '{:^3}'.format(orbit['order'])),
319 ('elements', '{:^18}'.format(elements)),
320 ('radius', '{:^8.4f}'.format(orbit['radius'])),
321 ('prototype', '{:^18}'.format(str(orbit['prototype_cluster']))),
322 ('clusters', '{:^4}'.format(orbit['n_clusters'])),
323 ('parameters', '{:^3}'.format(len(orbit['eigentensors']))),
324 ])
326 s = []
327 for name, value in fields.items():
328 n = max(len(name), len(value))
329 if index < 0:
330 s += ['{s:^{n}}'.format(s=name, n=n)]
331 else:
332 s += ['{s:^{n}}'.format(s=value, n=n)]
333 return ' | '.join(s)
335 # table header
336 width = max(len(str_orbit(-1, orbits[-1])), len(str_orbit(0, orbits[-1])))
337 print(' List of Orbits '.center(width, '='))
338 print(str_orbit(-1, orbits[0]))
339 print(''.center(width, '-'))
341 # table body
342 for i, orbit in enumerate(orbits):
343 print(str_orbit(i, orbit))
344 print(''.center(width, '='))
346 def __str__(self):
348 def str_order(order, header=False):
349 formats = {'order': '{:2}',
350 'n_orbits': '{:5}',
351 'n_clusters': '{:5}'}
352 s = []
353 for name, value in order.items():
354 str_repr = formats[name].format(value)
355 n = max(len(name), len(str_repr))
356 if header:
357 s += ['{s:^{n}}'.format(s=name, n=n)]
358 else:
359 s += ['{s:^{n}}'.format(s=str_repr, n=n)]
360 return ' | '.join(s)
362 # collect data
363 orbits = self.orbit_data
364 orders = self.cutoffs.orders
366 order_data = {o: dict(order=o, n_orbits=0, n_clusters=0) for o in orders}
367 for orbit in orbits:
368 o = orbit['order']
369 order_data[o]['n_orbits'] += 1
370 order_data[o]['n_clusters'] += orbit['n_clusters']
372 # prototype with max order to find column width
373 max_order = max(orders)
374 prototype = order_data[max_order]
375 n = max(len(str_order(prototype)), 54)
377 # basic information
378 s = []
379 s.append(' Cluster Space '.center(n, '='))
380 data = [('Spacegroup', self.spacegroup),
381 ('symprec', self.symprec),
382 ('Sum rules', self.acoustic_sum_rules),
383 ('Length scale', self.length_scale),
384 ('Cutoffs', self.cutoffs),
385 ('Cell', self.primitive_structure.cell),
386 ('Basis', self.primitive_structure.basis),
387 ('Numbers', self.primitive_structure.numbers),
388 ('Total number of orbits', len(orbits)),
389 ('Total number of clusters',
390 sum([order_data[order]['n_clusters'] for order in orders])),
391 ('Total number of parameters', self._get_n_dofs()
392 )]
393 for field, value in data:
394 if str(value).count('\n') > 1:
395 s.append('{:26} :\n{}'.format(field, value))
396 else:
397 s.append('{:26} : {}'.format(field, value))
399 # table header
400 s.append(''.center(n, '-'))
401 s.append(str_order(prototype, header=True))
402 s.append(''.center(n, '-'))
403 for order in orders:
404 s.append(str_order(order_data[order]).rstrip())
405 s.append(''.center(n, '='))
406 return '\n'.join(s)
408 def __repr__(self):
409 s = 'ClusterSpace({!r}, {!r}, {!r})'
410 return s.format(self.primitive_structure, self.cutoffs, self._config)
412 def copy(self):
413 return deepcopy(self)
415 def write(self, fileobj):
416 """ Writes cluster space to file.
418 The instance is saved into a custom format based on tar-files. The
419 resulting file will be a valid tar file and can be browsed by by a tar
420 reader. The included objects are themself either pickles, npz or other
421 tars.
423 Parameters
424 ----------
425 fileobj : str or file-like object
426 If the input is a string a tar archive will be created in the
427 current directory. Otherwise the input must be a valid file
428 like object.
429 """
430 # Create a tar archive
431 if isinstance(fileobj, str):
432 tar_file = tarfile.open(name=fileobj, mode='w')
433 else:
434 tar_file = tarfile.open(fileobj=fileobj, mode='w')
436 # Attributes in pickle format
437 pickle_attributes = ['_config',
438 '_symmetry_dataset', '_permutations',
439 '_atom_list', '_cluster_list']
440 items_pickle = dict()
441 for attribute in pickle_attributes:
442 items_pickle[attribute] = self.__getattribute__(attribute)
443 add_items_to_tarfile_pickle(tar_file, items_pickle, 'attributes')
445 # Constraint matrices and vectors in hdf5 format
446 items = dict(cvs=self._cvs)
447 add_items_to_tarfile_pickle(tar_file, items, 'constraint_vectors')
449 # Cutoffs and prim with their builtin write/read functions
450 items_custom = {'_cutoffs': self._cutoffs, '_prim': self._prim}
451 add_items_to_tarfile_custom(tar_file, items_custom)
453 # Orbits
454 add_list_to_tarfile_custom(tar_file, self._orbits, 'orbits')
455 add_list_to_tarfile_custom(tar_file, self._dropped_orbits,
456 'dropped_orbits')
458 # Done!
459 tar_file.close()
461 def read(f):
462 """ Reads a cluster space from file.
464 Parameters
465 ----------
466 f : str or file object
467 name of input file (str) or stream to load from (file object)
468 """
470 # Instantiate empty cs obj.
471 cs = ClusterSpace.__new__(ClusterSpace)
473 # Load from file on disk or file-like
474 if type(f) is str:
475 tar_file = tarfile.open(mode='r', name=f)
476 else:
477 tar_file = tarfile.open(mode='r', fileobj=f)
479 # Attributes
480 attributes = read_items_pickle(tar_file, 'attributes')
481 for name, value in attributes.items():
482 cs.__setattr__(name, value)
484 # Load the constraint matrices into their dict
485 items = read_items_pickle(tar_file, 'constraint_vectors')
486 cs._cvs = items['cvs']
488 # Cutoffs and prim via custom save funcs
489 fileobj = tar_file.extractfile('_cutoffs')
490 cs._cutoffs = Cutoffs.read(fileobj)
492 fileobj = tar_file.extractfile('_prim')
493 cs._prim = Atoms.read(fileobj)
495 # Orbits are stored in a separate archive
496 cs._orbits = read_list_custom(tar_file, 'orbits', Orbit.read)
497 cs._dropped_orbits = read_list_custom(
498 tar_file, 'dropped_orbits', Orbit.read)
500 tar_file.close()
502 # create summary object based on CS
503 cs.summary = ClusterSpaceData(cs)
505 # Done!
506 return cs