Coverage for hiphive/force_constant_potential.py: 96%
167 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-08-01 17:04 +0000
1"""
2This module introduces the ForceConstantPotential object which acts as the
3finalized force constant model.
4"""
6import pickle
7import tarfile
8import numpy as np
10from collections import Counter
11from typing import Any, Dict, List
12from .core.atoms import Atoms
13from .force_constant_model import ForceConstantModel
14from .core.orbits import Orbit
15from .core.orbits import OrientationFamily
16from .core.tensors import rotation_to_cart_coord, rotate_tensor
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 .input_output.logging_tools import logger
25logger = logger.getChild('ForceConstantPotential')
28class ForceConstantPotential:
29 """ A finalized force constant model. Can produce force constants for any
30 structure compatible with the structure for which the model was set up.
32 Parameters
33 ----------
34 cs : ClusterSpace
35 The cluster space the model is based upon.
36 parameters : numpy.ndarray
37 The fitted parameters.
38 metadata : dict
39 Meta data; will be pickled when object is written to file.
40 """
42 def __init__(self, cs, parameters, metadata=None):
44 self._prim = cs.primitive_structure.copy()
45 self.cluster_list = cs.cluster_list.copy()
46 self.atom_list = cs.atom_list.copy()
47 self.orbits = []
48 self.spacegroup = cs.spacegroup
49 self._config = cs._config
50 self.cs_summary = cs.summary
52 # add metadata
53 if metadata is None:
54 metadata = dict()
55 self._metadata = metadata
56 self._add_default_metadata()
58 # Extract the eigentensors from the cluster space and use the paramters
59 # to construct the finalized force constants
60 parameters = cs._map_parameters(parameters)
61 p = 0
62 for orb in cs.orbits:
63 new_orbit = Orbit()
64 fc = np.zeros(orb.eigentensors[0].shape)
65 for et, a in zip(orb.eigentensors, parameters[p:]):
66 fc += et * a
67 new_orbit.force_constant = fc
68 new_orbit.order = orb.order
69 new_orbit.radius = orb.radius
70 new_orbit.maximum_distance = orb.maximum_distance
71 for of in orb.orientation_families:
72 new_of = OrientationFamily()
73 new_of.cluster_indices = of.cluster_indices.copy()
74 sym_ind = of.symmetry_index
75 R = rotation_to_cart_coord(cs.rotation_matrices[sym_ind],
76 self.primitive_structure.cell)
77 fc = rotate_tensor(new_orbit.force_constant, R.T)
78 perm = cs.permutations[of.permutation_indices[0]]
79 new_of.force_constant = fc.transpose(perm)
80 new_orbit.orientation_families.append(new_of)
81 self.orbits.append(new_orbit)
82 p += len(orb.eigentensors)
84 @property
85 def symprec(self):
86 return self._config['symprec']
88 def write(self, filename: str):
89 """ Writes a ForceConstantPotential to file.
91 Parameters
92 ----------
93 filename
94 Name of file to write :class:`ForceConstantPotential` to.
95 """
97 # Create a tar archive
98 if isinstance(filename, str):
99 tar_file = tarfile.open(name=filename, mode='w')
100 else:
101 raise ValueError('filename must be str')
103 # objects with custom write
104 add_list_to_tarfile_custom(tar_file, self.orbits, 'orbits')
106 # prim with its builtin write/read functions. Note prim is a hiphive.core.Atoms object
107 items_custom = {'_prim': self._prim}
108 add_items_to_tarfile_custom(tar_file, items_custom)
110 # Attributes in pickle format
111 pickle_attributes = ['_config', '_metadata', 'spacegroup', 'atom_list',
112 'cluster_list', 'cs_summary']
113 items_pickle = dict()
114 for attribute in pickle_attributes:
115 items_pickle[attribute] = self.__getattribute__(attribute)
116 add_items_to_tarfile_pickle(tar_file, items_pickle, 'attributes')
118 # Done!
119 tar_file.close()
121 @staticmethod
122 def read(filename: str):
123 """ Reads a :class:`ForceConstantPotential` object from file.
125 Parameters
126 ----------
127 filename
128 Name of input file to load ForceConstantPotential from.
130 Returns
131 -------
132 ForceConstantPotential
133 The original object as stored in the file.
134 """
136 # Try usage of old read format, Remove in hiphive 1.0
137 try:
138 old_fcp = ForceConstantPotential._read_old(filename)
139 logger.warning('This fcp was written with a version <1.0. Please rewrite it.')
140 return old_fcp
141 except Exception:
142 pass
144 # Instantiate empty cs obj.
145 fcp = ForceConstantPotential.__new__(ForceConstantPotential)
147 # Load from file on disk
148 if type(filename) is str:
149 tar_file = tarfile.open(mode='r', name=filename)
150 else:
151 raise ValueError('filename must be str')
153 # Attributes with custom read
154 fileobj = tar_file.extractfile('_prim')
155 fcp._prim = Atoms.read(fileobj)
157 fcp.orbits = read_list_custom(tar_file, 'orbits', Orbit.read)
159 # Attributes
160 attributes = read_items_pickle(tar_file, 'attributes')
161 for name, value in attributes.items():
162 fcp.__setattr__(name, value)
164 # Done!
165 tar_file.close()
166 return fcp
168 # TODO: Remove?
169 def _write_old(self, f):
170 """Writes a force constant potential to file using old format.
172 Parameters
173 ----------
174 f : str or file object
175 name of input file (str) or stream to write to (file object)
176 """
177 if isinstance(f, str): 177 ↛ 181line 177 didn't jump to line 181 because the condition on line 177 was always true
178 with open(f, 'wb') as fobj:
179 pickle.dump(self, fobj)
180 else:
181 try:
182 pickle.dump(self, f)
183 except Exception:
184 raise Exception('Failed writing to file.')
186 @staticmethod
187 def _read_old(f):
188 """Reads a force constant potential from file in old format.
190 Parameters
191 ----------
192 f : str or file object
193 name of input file (str) or stream to load from (file object)
195 Returns
196 -------
197 ForceConstantPotential
198 the original object as stored in the file
199 """
200 if isinstance(f, str):
201 with open(f, 'rb') as fobj:
203 # This allows for reading FCPs with ASE-3.17 and 3.18
204 fcp = pickle.load(fobj)
205 _prim = fcp._prim
207 if hasattr(_prim, '_cell'): # 3.17 207 ↛ 208line 207 didn't jump to line 208 because the condition on line 207 was never true
208 cell = _prim._cell
209 else: # 3.18
210 cell = _prim.cell[:]
212 # assume PBC True (as it has to be True in hiphive)
213 pbc = [True, True, True]
215 new_prim = Atoms(
216 symbols=_prim.symbols, positions=_prim.positions, cell=cell, pbc=pbc)
217 fcp._prim = new_prim
218 return fcp
219 else:
220 return pickle.load(f)
222 @property
223 def metadata(self) -> dict[str, Any]:
224 """ Metadata associated with force constant potential. """
225 return self._metadata
227 @property
228 def primitive_structure(self) -> Atoms:
229 """ Atomic structure. """
230 return self._prim.copy()
232 @property
233 def orbit_data(self) -> List[Dict[str, Any]]:
234 """List of dictionaries containing detailed information for each
235 orbit, e.g., cluster radius and force constant.
236 """
237 data = []
238 for orbit_index, orbit in enumerate(self.orbits):
239 d = {}
240 d['index'] = orbit_index
241 d['order'] = orbit.order
242 d['radius'] = orbit.radius
243 d['maximum_distance'] = orbit.maximum_distance
244 d['n_clusters'] = len(orbit.orientation_families)
246 types = []
247 for atom_ind in self.cluster_list[orbit.prototype_index]:
248 types.append(self.primitive_structure.numbers[
249 self.atom_list[atom_ind].site])
250 d['prototype_cluster'] = self.cluster_list[orbit.prototype_index]
251 d['prototype_atom_types'] = types
253 d['geometrical_order'] = len(set(d['prototype_cluster']))
254 d['force_constant'] = orbit.force_constant
255 d['force_constant_norm'] = np.linalg.norm(orbit.force_constant)
256 data.append(d)
257 return data
259 def print_tables(self):
260 """ Prints information concerning the underlying cluster space to stdout,including,
261 e.g., the number of cluster, orbits, and parameters by order and number of bodies. """
262 self.cs_summary.print_tables()
264 def get_force_constants(self, atoms: Atoms) -> ForceConstantModel:
265 """ Return the force constants of a compatible structure.
267 Parameters
268 ----------
269 atoms
270 Input structure.
272 Returns
273 -------
274 Force constants in format for the input structure.
275 """
276 return ForceConstantModel(atoms, self).get_force_constants()
278 def __str__(self):
279 orbits = self.orbit_data
280 orbit_counts = Counter([orbit['order'] for orbit in orbits])
281 cluster_counts = Counter()
282 for orbit in orbits:
283 cluster_counts[orbit['order']] += orbit['n_clusters']
285 n = 54
286 s = []
287 s.append(' ForceConstantPotential '.center(n, '='))
288 s.append(f'Spacegroup {self.spacegroup}')
289 s.append(f'Cell:\n{self.primitive_structure.cell}')
290 s.append(f'Basis:\n{self.primitive_structure.basis}')
291 s.append(f'Numbers: {self.primitive_structure.numbers}')
293 s.append(f'Cutoff matrix:\n{self.cs_summary._cutoff_matrix}')
295 for order in sorted(orbit_counts.keys()):
296 n_dofs = self.cs_summary.ndofs_by_order[order]
297 s.append(f'Order {order}, #orbits {orbit_counts[order]}, #cluster '
298 f'{cluster_counts[order]}, #parameters {n_dofs}')
299 s.append(f'Total number of orbits: {len(orbits)}')
300 s.append(f'total number of clusters: {sum(cluster_counts.values())}')
301 s.append(f'total number of parameters: {sum(self.cs_summary.ndofs_by_order.values())}')
302 s.append(''.center(n, '='))
303 return '\n'.join(s)
305 def __repr__(self):
306 return 'ForceConstantPotential(ClusterSpace({!r}, ...), [...])'.format(
307 self.primitive_structure)
309 def _add_default_metadata(self):
310 """Adds default metadata to metadata dict."""
311 import getpass
312 import socket
313 from datetime import datetime
314 from . import __version__ as hiphive_version
316 self._metadata['date_created'] = datetime.now().strftime('%Y-%m-%dT%H:%M:%S')
317 self._metadata['username'] = getpass.getuser()
318 self._metadata['hostname'] = socket.gethostname()
319 self._metadata['hiphive_version'] = hiphive_version