Coverage for hiphive/force_constant_potential.py: 96%
167 statements
« prev ^ index » next coverage.py v7.6.8, created at 2024-11-28 11:20 +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 paramteres
38 metadata : dict
39 metadata dictionary, 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):
89 """ Writes a ForceConstantPotential to file.
91 Parameters
92 ----------
93 filename : str
94 name of file to write 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):
123 """ Reads a ForceConstantPotentialfrom file.
125 Parameters
126 ----------
127 filename : str
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):
224 """ dict : metadata associated with force constant potential """
225 return self._metadata
227 @property
228 def primitive_structure(self):
229 """ ase.Atoms : 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):
265 """ Return the force constants of a compatible structure.
267 Parameters
268 ----------
269 atoms : ase.Atoms
270 input structure
272 Returns
273 -------
274 ForceConstants
275 force constants
276 """
277 return ForceConstantModel(atoms, self).get_force_constants()
279 def __str__(self):
280 orbits = self.orbit_data
281 orbit_counts = Counter([orbit['order'] for orbit in orbits])
282 cluster_counts = Counter()
283 for orbit in orbits:
284 cluster_counts[orbit['order']] += orbit['n_clusters']
286 n = 54
287 s = []
288 s.append(' ForceConstantPotential '.center(n, '='))
289 s.append(f'Spacegroup {self.spacegroup}')
290 s.append(f'Cell:\n{self.primitive_structure.cell}')
291 s.append(f'Basis:\n{self.primitive_structure.basis}')
292 s.append(f'Numbers: {self.primitive_structure.numbers}')
294 s.append(f'Cutoff matrix:\n{self.cs_summary._cutoff_matrix}')
296 for order in sorted(orbit_counts.keys()):
297 n_dofs = self.cs_summary.ndofs_by_order[order]
298 s.append(f'Order {order}, #orbits {orbit_counts[order]}, #cluster '
299 f'{cluster_counts[order]}, #parameters {n_dofs}')
300 s.append(f'Total number of orbits: {len(orbits)}')
301 s.append(f'total number of clusters: {sum(cluster_counts.values())}')
302 s.append(f'total number of parameters: {sum(self.cs_summary.ndofs_by_order.values())}')
303 s.append(''.center(n, '='))
304 return '\n'.join(s)
306 def __repr__(self):
307 return 'ForceConstantPotential(ClusterSpace({!r}, ...), [...])'.format(
308 self.primitive_structure)
310 def _add_default_metadata(self):
311 """Adds default metadata to metadata dict."""
312 import getpass
313 import socket
314 from datetime import datetime
315 from . import __version__ as hiphive_version
317 self._metadata['date_created'] = datetime.now().strftime('%Y-%m-%dT%H:%M:%S')
318 self._metadata['username'] = getpass.getuser()
319 self._metadata['hostname'] = socket.gethostname()
320 self._metadata['hiphive_version'] = hiphive_version