"""
This module introduces the ForceConstantPotential object which acts as the
finalized force constant model.
"""
import pickle
import tarfile
import numpy as np
from collections import Counter
from typing import Any, Dict, List
from .core.atoms import Atoms
from .force_constant_model import ForceConstantModel
from .core.orbits import Orbit
from .core.orbits import OrientationFamily
from .core.tensors import rotation_to_cart_coord, rotate_tensor
from .input_output.read_write_files import (add_items_to_tarfile_pickle,
add_items_to_tarfile_custom,
add_list_to_tarfile_custom,
read_items_pickle,
read_list_custom)
from .input_output.logging_tools import logger
logger = logger.getChild('ForceConstantPotential')
[docs]
class ForceConstantPotential:
""" A finalized force constant model. Can produce force constants for any
structure compatible with the structure for which the model was set up.
Parameters
----------
cs : ClusterSpace
The cluster space the model is based upon.
parameters : numpy.ndarray
The fitted parameters.
metadata : dict
Meta data; will be pickled when object is written to file.
"""
def __init__(self, cs, parameters, metadata=None):
self._prim = cs.primitive_structure.copy()
self.cluster_list = cs.cluster_list.copy()
self.atom_list = cs.atom_list.copy()
self.orbits = []
self.spacegroup = cs.spacegroup
self._config = cs._config
self.cs_summary = cs.summary
# add metadata
if metadata is None:
metadata = dict()
self._metadata = metadata
self._add_default_metadata()
# Extract the eigentensors from the cluster space and use the paramters
# to construct the finalized force constants
parameters = cs._map_parameters(parameters)
p = 0
for orb in cs.orbits:
new_orbit = Orbit()
fc = np.zeros(orb.eigentensors[0].shape)
for et, a in zip(orb.eigentensors, parameters[p:]):
fc += et * a
new_orbit.force_constant = fc
new_orbit.order = orb.order
new_orbit.radius = orb.radius
new_orbit.maximum_distance = orb.maximum_distance
for of in orb.orientation_families:
new_of = OrientationFamily()
new_of.cluster_indices = of.cluster_indices.copy()
sym_ind = of.symmetry_index
R = rotation_to_cart_coord(cs.rotation_matrices[sym_ind],
self.primitive_structure.cell)
fc = rotate_tensor(new_orbit.force_constant, R.T)
perm = cs.permutations[of.permutation_indices[0]]
new_of.force_constant = fc.transpose(perm)
new_orbit.orientation_families.append(new_of)
self.orbits.append(new_orbit)
p += len(orb.eigentensors)
@property
def symprec(self):
return self._config['symprec']
[docs]
def write(self, filename: str):
""" Writes a ForceConstantPotential to file.
Parameters
----------
filename
Name of file to write :class:`ForceConstantPotential` to.
"""
# Create a tar archive
if isinstance(filename, str):
tar_file = tarfile.open(name=filename, mode='w')
else:
raise ValueError('filename must be str')
# objects with custom write
add_list_to_tarfile_custom(tar_file, self.orbits, 'orbits')
# prim with its builtin write/read functions. Note prim is a hiphive.core.Atoms object
items_custom = {'_prim': self._prim}
add_items_to_tarfile_custom(tar_file, items_custom)
# Attributes in pickle format
pickle_attributes = ['_config', '_metadata', 'spacegroup', 'atom_list',
'cluster_list', 'cs_summary']
items_pickle = dict()
for attribute in pickle_attributes:
items_pickle[attribute] = self.__getattribute__(attribute)
add_items_to_tarfile_pickle(tar_file, items_pickle, 'attributes')
# Done!
tar_file.close()
[docs]
@staticmethod
def read(filename: str):
""" Reads a :class:`ForceConstantPotential` object from file.
Parameters
----------
filename
Name of input file to load ForceConstantPotential from.
Returns
-------
ForceConstantPotential
The original object as stored in the file.
"""
# Try usage of old read format, Remove in hiphive 1.0
try:
old_fcp = ForceConstantPotential._read_old(filename)
logger.warning('This fcp was written with a version <1.0. Please rewrite it.')
return old_fcp
except Exception:
pass
# Instantiate empty cs obj.
fcp = ForceConstantPotential.__new__(ForceConstantPotential)
# Load from file on disk
if type(filename) is str:
tar_file = tarfile.open(mode='r', name=filename)
else:
raise ValueError('filename must be str')
# Attributes with custom read
fileobj = tar_file.extractfile('_prim')
fcp._prim = Atoms.read(fileobj)
fcp.orbits = read_list_custom(tar_file, 'orbits', Orbit.read)
# Attributes
attributes = read_items_pickle(tar_file, 'attributes')
for name, value in attributes.items():
fcp.__setattr__(name, value)
# Done!
tar_file.close()
return fcp
# TODO: Remove?
def _write_old(self, f):
"""Writes a force constant potential to file using old format.
Parameters
----------
f : str or file object
name of input file (str) or stream to write to (file object)
"""
if isinstance(f, str):
with open(f, 'wb') as fobj:
pickle.dump(self, fobj)
else:
try:
pickle.dump(self, f)
except Exception:
raise Exception('Failed writing to file.')
@staticmethod
def _read_old(f):
"""Reads a force constant potential from file in old format.
Parameters
----------
f : str or file object
name of input file (str) or stream to load from (file object)
Returns
-------
ForceConstantPotential
the original object as stored in the file
"""
if isinstance(f, str):
with open(f, 'rb') as fobj:
# This allows for reading FCPs with ASE-3.17 and 3.18
fcp = pickle.load(fobj)
_prim = fcp._prim
if hasattr(_prim, '_cell'): # 3.17
cell = _prim._cell
else: # 3.18
cell = _prim.cell[:]
# assume PBC True (as it has to be True in hiphive)
pbc = [True, True, True]
new_prim = Atoms(
symbols=_prim.symbols, positions=_prim.positions, cell=cell, pbc=pbc)
fcp._prim = new_prim
return fcp
else:
return pickle.load(f)
@property
def metadata(self) -> dict[str, Any]:
""" Metadata associated with force constant potential. """
return self._metadata
@property
def primitive_structure(self) -> Atoms:
""" Atomic structure. """
return self._prim.copy()
@property
def orbit_data(self) -> List[Dict[str, Any]]:
"""List of dictionaries containing detailed information for each
orbit, e.g., cluster radius and force constant.
"""
data = []
for orbit_index, orbit in enumerate(self.orbits):
d = {}
d['index'] = orbit_index
d['order'] = orbit.order
d['radius'] = orbit.radius
d['maximum_distance'] = orbit.maximum_distance
d['n_clusters'] = len(orbit.orientation_families)
types = []
for atom_ind in self.cluster_list[orbit.prototype_index]:
types.append(self.primitive_structure.numbers[
self.atom_list[atom_ind].site])
d['prototype_cluster'] = self.cluster_list[orbit.prototype_index]
d['prototype_atom_types'] = types
d['geometrical_order'] = len(set(d['prototype_cluster']))
d['force_constant'] = orbit.force_constant
d['force_constant_norm'] = np.linalg.norm(orbit.force_constant)
data.append(d)
return data
[docs]
def print_tables(self):
""" Prints information concerning the underlying cluster space to stdout,including,
e.g., the number of cluster, orbits, and parameters by order and number of bodies. """
self.cs_summary.print_tables()
[docs]
def get_force_constants(self, atoms: Atoms) -> ForceConstantModel:
""" Return the force constants of a compatible structure.
Parameters
----------
atoms
Input structure.
Returns
-------
Force constants in format for the input structure.
"""
return ForceConstantModel(atoms, self).get_force_constants()
def __str__(self):
orbits = self.orbit_data
orbit_counts = Counter([orbit['order'] for orbit in orbits])
cluster_counts = Counter()
for orbit in orbits:
cluster_counts[orbit['order']] += orbit['n_clusters']
n = 54
s = []
s.append(' ForceConstantPotential '.center(n, '='))
s.append(f'Spacegroup {self.spacegroup}')
s.append(f'Cell:\n{self.primitive_structure.cell}')
s.append(f'Basis:\n{self.primitive_structure.basis}')
s.append(f'Numbers: {self.primitive_structure.numbers}')
s.append(f'Cutoff matrix:\n{self.cs_summary._cutoff_matrix}')
for order in sorted(orbit_counts.keys()):
n_dofs = self.cs_summary.ndofs_by_order[order]
s.append(f'Order {order}, #orbits {orbit_counts[order]}, #cluster '
f'{cluster_counts[order]}, #parameters {n_dofs}')
s.append(f'Total number of orbits: {len(orbits)}')
s.append(f'total number of clusters: {sum(cluster_counts.values())}')
s.append(f'total number of parameters: {sum(self.cs_summary.ndofs_by_order.values())}')
s.append(''.center(n, '='))
return '\n'.join(s)
def __repr__(self):
return 'ForceConstantPotential(ClusterSpace({!r}, ...), [...])'.format(
self.primitive_structure)
def _add_default_metadata(self):
"""Adds default metadata to metadata dict."""
import getpass
import socket
from datetime import datetime
from . import __version__ as hiphive_version
self._metadata['date_created'] = datetime.now().strftime('%Y-%m-%dT%H:%M:%S')
self._metadata['username'] = getpass.getuser()
self._metadata['hostname'] = socket.gethostname()
self._metadata['hiphive_version'] = hiphive_version
