"""
This module provides functionality for storing and handling of force constants.
"""
import numpy as np
import tarfile
from abc import ABC, abstractmethod
from itertools import product
from string import ascii_lowercase, ascii_uppercase
from typing import List, Tuple
from scipy.linalg import eig
from ase import Atoms, units
from .input_output import shengBTE as shengBTEIO
from .input_output import phonopy as phonopyIO
from .input_output.read_write_files import add_items_to_tarfile_pickle, read_items_pickle,\
add_ase_atoms_to_tarfile, read_ase_atoms_from_tarfile
from .input_output import gpumd as gpumdIO
[docs]class ForceConstants(ABC):
""" Base class for force constants """
def __init__(self, supercell: Atoms):
if not all(supercell.pbc):
raise ValueError('configuration must have periodic boundary conditions')
self._supercell = supercell.copy()
@abstractmethod
def __getitem__(self):
pass
@property
def n_atoms(self) -> int:
""" number of atoms """
return len(self.supercell)
@property
def supercell(self) -> Atoms:
""" supercell associated with force constants """
return self._supercell.copy()
@property
def clusters(self) -> list:
""" sorted list of clusters """
return sorted(self._fc_dict.keys(), key=lambda c: (len(c), c))
def __len__(self) -> int:
""" int : number of clusters (or force constants) """
return len(self._fc_dict)
def __add__(self, fcs_other) -> None:
""" ForceConstants : addition of two force constants """
if type(self) != type(fcs_other):
raise ValueError('ForceConstants objects are of different types')
# if they have overlapping orders, raise
if any(order in self.orders for order in fcs_other.orders):
raise ValueError('ForceConstants objects share at least one order')
# check that they share the same supercell
if self.supercell != fcs_other.supercell:
raise ValueError('supercells do not match')
# add force constants
fc_merged = {**self._fc_dict, **fcs_other._fc_dict}
return self.__class__(fc_merged, self.supercell)
[docs] def get_fc_dict(self, order: int = None) -> dict:
""" Returns force constant dictionary for one specific order.
The returned dict may be sparse or may be dense depending on the
underlying force constants.
Parameters
----------
order
force constants returned for this order
Returns
-------
dictionary with keys corresponding to clusters and values to
respective force constant
"""
if order is None:
return self._fc_dict
if order not in self.orders:
raise ValueError('Order {} not in ForceConstants'.format(order))
fc_order = {}
for c, fc in self._fc_dict.items():
if len(c) == order:
fc_order[c] = fc
return fc_order
[docs] def get_fc_array(self, order: int, format: str = 'phonopy') -> np.ndarray:
""" Returns force constants in array format for specified order.
Parameters
----------
order
force constants for this order will be returned
format
specify which format (shape) the NumPy array should have,
possible values are `phonopy` and `ase`
Returns
-------
NumPy array with shape `(N,)*order + (3,)*order` where `N` is
the number of atoms
"""
if order not in self.orders:
raise ValueError('Order not in orders')
if format not in ['ase', 'phonopy']:
raise ValueError('Format must be either ase or phonopy')
# generate fc array for phonopy format
fc_array = np.zeros((self.n_atoms, ) * order + (3, ) * order)
for cluster in product(range(self.n_atoms), repeat=order):
fc_array[cluster] = self[cluster]
if format == 'ase':
if order != 2:
raise ValueError('ASE format works only for order 2')
return fc_array.transpose([0, 2, 1, 3]).reshape(
self.n_atoms * 3, self.n_atoms * 3)
else:
return fc_array
[docs] def compute_gamma_frequencies(self) -> np.ndarray:
""" Returns the Gamma frequencies in THz using the second-order force
constants. """
fc2_array = self.get_fc_array(order=2)
masses = self.supercell.get_masses()
return _compute_gamma_frequencies(fc2_array, masses)
[docs] def assert_acoustic_sum_rules(self, order: int = None, tol: float = 1e-6):
""" Asserts that force constants obey acoustic sum rules.
Parameters
----------
order
specifies which order to check, if None all are checked
tol
numeric tolerance for checking sum rules
Raises
------
AssertionError
if acoustic sum rules are violated
"""
# set up orders
if order is None:
orders = self.orders
else:
if order not in self.orders:
raise ValueError('Order not available in FCS')
orders = [order]
atomic_indices = range(self.n_atoms)
for order in orders:
assert_msg = 'Acoustic sum rule for order {} violated for atom'
assert_msg += ' {}' * (order - 1) + ' x'
for ijk in product(atomic_indices, repeat=order-1):
fc_sum_ijk = np.zeros((3, )*order)
for m in atomic_indices:
cluster = ijk + (m, )
fc_sum_ijk += self[cluster]
assert np.linalg.norm(fc_sum_ijk) < tol, assert_msg.format(order, *ijk)
[docs] def print_force_constant(self, cluster: Tuple[int]) -> None:
"""
Prints force constants for a cluster in a nice format.
Parameters
----------
cluster
sites belonging to the cluster
"""
print(self._repr_fc(cluster))
def __eq__(self, other):
# check supercells are the same
if not len(self.supercell) == len(other.supercell):
return False
if not np.allclose(self.supercell.positions, other.supercell.positions):
return False
if not np.allclose(self.supercell.cell, other.supercell.cell):
return False
if not all(self.supercell.numbers == other.supercell.numbers):
return False
# check orders and clusters are the same
if self.orders != other.orders:
return False
if self.clusters != other.clusters:
return False
# check force constants
for c in self.clusters:
if not np.allclose(self[c], other[c]):
return False
return True
def __str__(self) -> str:
s = []
s.append(' ForceConstants '.center(54, '='))
s.append('Orders: {}'.format(self.orders))
s.append('Atoms: {}'.format(self.supercell))
s.append('')
if len(self) > 10:
for c in self.clusters[:3]:
s.append(self._repr_fc(c)+'\n')
s.append('...\n')
for c in self.clusters[-3:]:
s.append(self._repr_fc(c)+'\n')
else:
for c in self.clusters:
s.append(self._repr_fc(c)+'\n')
return '\n'.join(s)
def __repr__(self) -> str:
fc_dict_str = '{{{}: {}, ..., {}: {}}}'.format(
self.clusters[0], self[self.clusters[0]].round(5),
self.clusters[-1], self[self.clusters[-1]].round(5))
return ('ForceConstants(fc_dict={}, atoms={!r})'
.format(fc_dict_str, self.supercell))
def _repr_fc(self, cluster: Tuple[int],
precision: float = 5, suppress: bool = True) -> str:
"""
Representation for single cluster and its force constant.
Parameters
----------
cluster
tuple of ints indicating the sites belonging to the cluster
"""
s = []
s.append('Cluster: {}'.format(cluster))
for atom_index in cluster:
s.append(str(self.supercell[atom_index]))
s.append('Force constant:')
s.append(np.array_str(self[cluster], precision=precision,
suppress_small=suppress))
return '\n'.join(s)
def _sanity_check_dict(self, fc_dict: dict) -> None:
""" Checks that all indices in clusters are between 0 and number of
atoms. """
# Use flat numpy array for fast check
cluster_indices = np.concatenate([cluster for cluster in fc_dict.keys()])
if not np.all((0 <= cluster_indices) & (cluster_indices < self.n_atoms)):
# If the check failed, do slower list processing to find the erring cluster
for cluster in fc_dict.keys():
if not all(0 <= i < self.n_atoms for i in cluster):
raise ValueError('Cluster {} not in supercell'.format(cluster))
[docs] @classmethod
def from_arrays(cls, supercell: Atoms,
fc2_array: np.ndarray = None, fc3_array: np.ndarray = None):
""" Constructs FCs from numpy arrays.
One or both of fc2_array and fc3_array must not be None
Parameters
----------
supercell
supercell structure
fc2_array
second-order force constant in phonopy format, i.e. must have shape (N, N, 3, 3)
fc3_array
third-order force constant in phonopy format, i.e. must have shape (N, N, N, 3, 3, 3)
"""
if fc2_array is None and fc3_array is None:
raise ValueError('Please provide force constant arrays.')
n_atoms = len(supercell)
if fc2_array is None:
fc2_dict = dict()
else:
if fc2_array.shape != (n_atoms, n_atoms, 3, 3):
raise ValueError('fc2 array has wrong shape')
fc2_dict = array_to_dense_dict(fc2_array)
if fc3_array is None:
fc3_dict = dict()
else:
if fc3_array.shape != (n_atoms, n_atoms, n_atoms, 3, 3, 3):
raise ValueError('fc2 array has wrong shape')
fc3_dict = array_to_dense_dict(fc3_array)
fc_dict = {**fc2_dict, **fc3_dict}
return cls.from_dense_dict(fc_dict, supercell)
[docs] @classmethod
def from_sparse_dict(cls, fc_dict: dict, supercell: Atoms):
""" Assumes label symmetries, meaning only one cluster for each
permuation should be included
Parameters
----------
fc_dict
keys corresponding to clusters and values to the force constants
supercell
atomic configuration
"""
return SortedForceConstants(fc_dict, supercell=supercell)
[docs] @classmethod
def from_dense_dict(cls, fc_dict: dict, supercell: Atoms):
""" All permutations of clusters that are not zero must be listed,
if label symmetries are fullfilled will return a SortedForceConstants
Parameters
----------
fc_dict
keys corresponding to clusters and values to the force constants
supercell
atomic configuration
"""
if check_label_symmetries(fc_dict):
fc_sparse_dict = dense_dict_to_sparse_dict(fc_dict)
return SortedForceConstants(fc_sparse_dict, supercell=supercell)
else:
return RawForceConstants(fc_dict, supercell=supercell)
[docs] @classmethod
def read_phonopy(cls, supercell: Atoms, fname: str, format: str = None):
""" Reads force constants from a phonopy calculation.
Parameters
----------
supercell
supercell structure (`SPOSCAR`)
fname
name of second-order force constant file
format
format for second-order force constants;
possible values: "text", "hdf5"
"""
fc2_array = phonopyIO.read_phonopy_fc2(fname, format=format)
return cls.from_arrays(supercell, fc2_array)
[docs] @classmethod
def read_phono3py(cls, supercell: Atoms, fname: str):
""" Reads force constants from a phono3py calculation.
Parameters
----------
supercell
supercell structure (`SPOSCAR`)
fname
name of third-order force constant file
"""
fc3_array = phonopyIO.read_phonopy_fc3(fname)
return cls.from_arrays(supercell, fc3_array=fc3_array)
[docs] @classmethod
def read_shengBTE(cls, supercell: Atoms, fname: str, prim: Atoms, symprec=1e-5):
""" Reads third order force constants from a shengBTE calculation.
shengBTE force constants will be mapped onto a supercell.
Parameters
----------
supercell
supercell structure
fname
name of third-order force constant file
prim
primitive configuration (must be equivalent to structure used in
the shengBTE calculation)
symprec
structural symmetry tolerance
"""
fcs = shengBTEIO.read_shengBTE_fc3(fname, prim, supercell, symprec)
return fcs
[docs] @classmethod
def read(cls, fname: str):
""" Reads ForceConstants from file.
Parameters
----------
fname
name of file from which to read
"""
tar_file = tarfile.open(mode='r', name=fname)
items = read_items_pickle(tar_file, 'fc_dict')
fc_dict = items['fc_dict']
fcs_type = items['fcs_type']
supercell = read_ase_atoms_from_tarfile(tar_file, 'supercell')
tar_file.close()
if fcs_type == 'SortedForceConstants':
return SortedForceConstants(fc_dict, supercell)
elif fcs_type == 'RawForceConstants':
return RawForceConstants(fc_dict, supercell)
else:
raise ValueError('FCS type not recongnized')
[docs] def write(self, fname: str) -> None:
""" Writes entire ForceConstants object to file.
Parameters
----------
fname
name of file to which to write
"""
tar_file = tarfile.open(name=fname, mode='w')
items_pickle = dict(fc_dict=self._fc_dict, fcs_type=self.__class__.__name__)
add_items_to_tarfile_pickle(tar_file, items_pickle, 'fc_dict')
add_ase_atoms_to_tarfile(tar_file, self.supercell, 'supercell')
tar_file.close()
[docs] def write_to_phonopy(self, fname: str, format: str = None) -> None:
"""
Writes force constants in phonopy format.
Parameters
----------
fname
name of file to which to write second-order force constant
format
format for second-order force constants;
possible values: "text", "hdf5"
"""
phonopyIO.write_phonopy_fc2(fname, self, format=format)
[docs] def write_to_phono3py(self, fname: str) -> None:
"""
Writes force constants in phono3py format.
Parameters
----------
fname
name of file to which to write third-order force constant
"""
phonopyIO.write_phonopy_fc3(fname, self)
[docs] def write_to_shengBTE(self, fname: str, prim: Atoms, **kwargs) -> None:
"""
Writes third order force constants in shengBTE format.
Parameters
----------
fname
name of file to which to write third-order force constant
prim
primitive configuration (must be equivalent to structure used in
the shengBTE calculation)
"""
shengBTEIO.write_shengBTE_fc3(fname, self, prim, **kwargs)
[docs]class SortedForceConstants(ForceConstants):
""" Force constants with label symmetries.
Parameters
----------
fc_dict : dict
keys corresponding to clusters and values to the force constants,
should only contain sorted clusters
supercell : ase.Atoms
"""
def __init__(self, fc_dict: dict, supercell: Atoms) -> None:
super().__init__(supercell)
self._sanity_check_dict(fc_dict)
self._fc_dict = fc_dict
self._orders = sorted(set(len(c) for c in self._fc_dict.keys()))
def __getitem__(self, cluster: Tuple[int]):
sorted_cluster = tuple(sorted(cluster))
# cluster not in fcs
if sorted_cluster not in self._fc_dict.keys():
return np.zeros((3,)*len(cluster))
# return fc for the unsorted cluster
inv_perm = np.argsort(np.argsort(cluster))
return self._fc_dict[sorted_cluster].transpose(inv_perm)
@property
def orders(self) -> List[int]:
""" orders for which force constants exist """
return self._orders.copy()
def _sanity_check_dict(self, fc_dict: dict) -> None:
super()._sanity_check_dict(fc_dict)
# also check clusters are sorted
for cluster in fc_dict.keys():
if cluster != tuple(sorted(cluster)):
raise ValueError('Found unsorted cluster {}'.format(cluster))
[docs] def write_to_GPUMD(self, fname_fc, fname_clusters, order, tol=1e-10):
"""
Writes force constants of the specified order in GPUMD format.
Parameters
----------
fname_fc : str
name of file which contains the lookup force constants
fname_clusters : str
name of file which contains the clusters and the fc lookup index
order : int
force constants for this order will be written to file
tol : float
if the norm of a force constant is less than tol then it is not written.
if two force-constants are within tol; they are considered equal.
"""
gpumdIO.write_fcs_gpumd(fname_fc, fname_clusters, self, order, tol)
[docs]class RawForceConstants(ForceConstants):
""" Force constants without label symmetries.
Parameters
----------
fc_dict : dict
keys corresponding to clusters and values to the force constants,
should contain all clusters with nonzero force constants
supercell : ase.Atoms
"""
def __init__(self, fc_dict: dict, supercell: Atoms) -> None:
super().__init__(supercell)
self._sanity_check_dict(fc_dict)
self._fc_dict = fc_dict
self._orders = sorted(set(len(c) for c in self._fc_dict.keys()))
def __getitem__(self, cluster: Tuple[int]):
if cluster not in self._fc_dict.keys():
return np.zeros((3,)*len(cluster))
return self._fc_dict[cluster]
@property
def orders(self) -> List[int]:
""" orders for which force constants exist """
return self._orders.copy()
# ====================== #
# HELPER FUNCTIONS #
# ====================== #
[docs]def array_to_dense_dict(fc_array: np.ndarray, fc_tol: float = 1e-10) -> dict:
""" Constructs a dense dict from an fc array in phonopy format.
Force constants with norm smaller than fc_tol will be considered zero and
therefore not included in the fc_dict.
Parameters
----------
fc_array
force constant array in phonopy format
fc_tol
tolerance for considering force constants zero or not
"""
# sanity check
n_atoms = fc_array.shape[0]
order = int(len(fc_array.shape) / 2)
if fc_array.shape != (n_atoms, ) * order + (3, ) * order:
raise ValueError('fc array has bad shape')
# construct dense dict
fc_dict = dict()
for cluster in product(range(n_atoms), repeat=order):
fc = fc_array[cluster]
if np.linalg.norm(fc) > fc_tol:
fc_dict[cluster] = fc
return fc_dict
[docs]def check_label_symmetries(fc_dict: dict) -> bool:
""" Checks label symmetries for dense fc dict.
TODO
----
tol, which one to use etc
Parameters
----------
fc_dict
keys corresponding to clusters and values to the force constants
"""
for cluster, fc in fc_dict.items():
inv_perm = np.argsort(np.argsort(cluster))
sorted_cluster = tuple(sorted(cluster))
if sorted_cluster not in fc_dict:
return False
if not np.allclose(fc, fc_dict[sorted_cluster].transpose(inv_perm)):
return False
return True
[docs]def dense_dict_to_sparse_dict(fc_dict: dict) -> dict:
""" Converts dense dict to sparse dict.
This does not check if label symmetry is True, but rather will just keep
the sorted clusters and their force constants.
Parameters
----------
fc_dict
keys corresponding to clusters and values to the force constants
"""
fc_dict_sparse = dict()
for cluster, fc in fc_dict.items():
if cluster == tuple(sorted(cluster)):
fc_dict_sparse[cluster] = fc
return fc_dict_sparse
[docs]def symbolize_force_constant(fc: np.ndarray, tol: float = 1e-10) -> np.ndarray:
"""Carries out a symbolic symmetrization of a force constant tensor.
Parameters
----------
fc
force constant tensor
tol
tolerance used to decide whether two elements are identical
Returns
-------
symbolic representation of force constant matrix
"""
fc_int = np.round(fc / tol).astype(int)
fc_chars = np.empty(fc_int.shape, dtype=object).flatten()
all_chars = ascii_lowercase + ascii_uppercase
lookup_chars = {}
for i, val in enumerate(fc_int.flatten()):
if val == 0:
fc_chars[i] = 0
elif val in lookup_chars.keys():
fc_chars[i] = lookup_chars[val]
elif -val in lookup_chars.keys():
fc_chars[i] = '-{:}'.format(lookup_chars[-val])
else:
lookup_chars[val] = all_chars[len(lookup_chars.keys())]
fc_chars[i] = lookup_chars[val]
return fc_chars.reshape(fc_int.shape)
def _compute_gamma_frequencies(fc2: np.ndarray, masses: List[float]) -> np.ndarray:
""" Computes Gamma frequencies using second-order force constants.
Assumes fc2 is in units of eV/A2.
Parameters
----------
fc2
second-order force constants in phonopy format
masses
mass of each atom
Returns
-------
Gamma frequencies in THz
"""
n_atoms = fc2.shape[0]
if len(masses) != n_atoms:
raise ValueError('Length of masses not compatible with fc2')
mass_matrix = np.sqrt(np.outer(masses, masses))
# divide with mass matrix
fc2_tmp = np.zeros((n_atoms, n_atoms, 3, 3))
for pair in product(range(n_atoms), repeat=2):
fc2_tmp[pair] = fc2[pair] / mass_matrix[pair]
# reshape into matrix and solve eigenvalues
fc2_tmp = fc2_tmp.transpose([0, 2, 1, 3]).reshape(n_atoms * 3, n_atoms * 3)
eigen_vals, _ = eig(fc2_tmp)
eigen_vals *= 1e20 / units.J * units.kg # [1/s**2]
eigen_vals.sort()
# if negative eigenval, set frequency to negative
gamma_frequencies = []
for val in eigen_vals.real:
if val >= 0:
gamma_frequencies.append(np.sqrt(val))
else:
gamma_frequencies.append(-np.sqrt(np.abs(val)))
# Sort and convert to THz
gamma_frequencies = np.array(gamma_frequencies) / 1e12 / (2 * np.pi)
gamma_frequencies.sort()
return gamma_frequencies