import posixpath
from typing import Any, Optional
import numpy as np
from ase.atoms import Atoms
from phonopy import Phonopy
from phonopy.file_IO import write_FORCE_CONSTANTS
from atomistics.shared.output import OutputPhonons, OutputThermodynamic
from atomistics.workflows.interface import Workflow
from atomistics.workflows.phonons.helper import (
analyse_results_for_harmonic_approximation,
get_band_structure,
get_hesse_matrix,
get_tasks_for_harmonic_approximation,
get_thermal_properties_for_harmonic_approximation,
plot_band_structure,
plot_dos,
)
[docs]
class PhonopyWorkflow(Workflow):
"""
Phonopy wrapper for the calculation of free energy in the framework of quasi harmonic approximation.
Get output via `get_thermal_properties()`.
Note:
- This class does not consider the thermal expansion. For this, use `QuasiHarmonicJob` (find more in its
docstring)
- Depending on the value given in `job.input['interaction_range']`, this class automatically changes the number of
atoms. The input parameters of the reference job might have to be set appropriately (e.g. use `k_mesh_density`
for DFT instead of setting k-points directly).
- The structure used in the reference job should be a relaxed structure.
- Theory behind it: https://en.wikipedia.org/wiki/Quasi-harmonic_approximation
"""
[docs]
def __init__(
self,
structure: Atoms,
interaction_range: float = 10.0,
displacement: float = 0.01,
dos_mesh: int = 20,
primitive_matrix: Optional[np.ndarray] = None,
number_of_snapshots: Optional[int] = None,
):
"""
Initialize the PhonopyWorkflow.
Args:
structure (Atoms): The structure used in the reference job.
interaction_range (float, optional): The interaction range. Defaults to 10.0.
displacement (float, optional): The displacement. Defaults to 0.01.
dos_mesh (int, optional): The DOS mesh. Defaults to 20.
primitive_matrix (np.ndarray, optional): The primitive matrix. Defaults to None.
number_of_snapshots (int, optional): The number of snapshots. Defaults to None.
"""
self._interaction_range = interaction_range
self._displacement = displacement
self._dos_mesh = dos_mesh
self._number_of_snapshots = number_of_snapshots
self.structure = structure
self._primitive_matrix = primitive_matrix
self.phonopy: Phonopy | None = None
self._phonopy_dict: dict[str, Any] = {}
def _get_phonopy(self) -> Phonopy:
if self.phonopy is None:
raise ValueError(
"Please call generate_structures() before accessing phonon results."
)
return self.phonopy
[docs]
def generate_structures(self) -> dict:
"""
Generate structures.
Returns:
dict: The generated structures.
"""
task_dict, self.phonopy = get_tasks_for_harmonic_approximation(
structure=self.structure,
primitive_matrix=self._primitive_matrix,
displacement=self._displacement,
number_of_snapshots=self._number_of_snapshots,
interaction_range=self._interaction_range,
)
return task_dict
[docs]
def analyse_structures(
self, output_dict: dict, output_keys: tuple[str] = OutputPhonons.keys()
) -> dict:
"""
Analyse structures.
Args:
output_dict (dict): The output dictionary.
output_keys (tuple[str], optional): The output keys. Defaults to OutputPhonons.keys().
Returns:
dict: The analysed structures.
"""
self._phonopy_dict = analyse_results_for_harmonic_approximation(
phonopy=self._get_phonopy(),
output_dict=output_dict,
dos_mesh=self._dos_mesh,
number_of_snapshots=self._number_of_snapshots,
output_keys=output_keys,
)
return self._phonopy_dict
[docs]
def get_thermal_properties(
self,
t_min: float = 1.0,
t_max: float = 1500.0,
t_step: float = 50.0,
temperatures: Optional[np.ndarray] = None,
cutoff_frequency: Optional[float] = None,
pretend_real: bool = False,
band_indices: Optional[np.ndarray] = None,
is_projection: bool = False,
output_keys: tuple[str] = OutputThermodynamic.keys(),
) -> dict:
"""
Get thermal properties.
Args:
t_min (float, optional): The minimum sample temperature. Defaults to 1.0.
t_max (float, optional): The maximum sample temperature. Defaults to 1500.0.
t_step (float, optional): The temperature sample interval. Defaults to 50.0.
temperatures (np.ndarray, optional): Custom array of temperature samples. Defaults to None.
cutoff_frequency (float, optional): The cutoff frequency. Defaults to None.
pretend_real (bool, optional): Whether to pretend real. Defaults to False.
band_indices (np.ndarray, optional): The band indices. Defaults to None.
is_projection (bool, optional): Whether it is a projection. Defaults to False.
output_keys (tuple[str], optional): The output keys. Defaults to OutputThermodynamic.keys().
Returns:
dict: The thermal properties.
"""
return get_thermal_properties_for_harmonic_approximation(
phonopy=self._get_phonopy(),
t_min=t_min,
t_max=t_max,
t_step=t_step,
temperatures=temperatures,
cutoff_frequency=cutoff_frequency,
pretend_real=pretend_real,
band_indices=band_indices,
is_projection=is_projection,
output_keys=output_keys,
)
[docs]
def get_dynamical_matrix(self, npoints: int = 101) -> np.ndarray:
"""
Get the dynamical matrix.
Args:
npoints (int, optional): The number of points. Defaults to 101.
Returns:
np.ndarray: The dynamical matrix.
"""
phonopy = self._get_phonopy()
phonopy.auto_band_structure(
npoints=npoints,
with_eigenvectors=False,
with_group_velocities=False,
plot=False,
write_yaml=False,
filename="band.yaml",
)
dynamical_matrix = phonopy.dynamical_matrix
if dynamical_matrix is None or dynamical_matrix.dynamical_matrix is None:
raise ValueError("Phonopy did not produce a dynamical matrix.")
return np.real_if_close(dynamical_matrix.dynamical_matrix)
[docs]
def dynamical_matrix_at_q(self, q: np.ndarray) -> np.ndarray:
"""
Get the dynamical matrix at a given q.
Args:
q (np.ndarray): The q value.
Returns:
np.ndarray: The dynamical matrix.
"""
return np.real_if_close(self._get_phonopy().get_dynamical_matrix_at_q(q))
[docs]
def write_phonopy_force_constants(
self, file_name: str = "FORCE_CONSTANTS", cwd: Optional[str] = None
):
"""
Write the Phonopy force constants.
Args:
file_name (str, optional): The file name. Defaults to "FORCE_CONSTANTS".
cwd (str, optional): The current working directory. Defaults to None.
"""
if cwd is not None:
file_name = posixpath.join(cwd, file_name)
phonopy = self._get_phonopy()
if phonopy.force_constants is None:
raise ValueError("Phonopy force constants are not available.")
write_FORCE_CONSTANTS(
force_constants=phonopy.force_constants, filename=file_name
)
[docs]
def get_hesse_matrix(self) -> np.ndarray:
"""
Get the Hesse matrix.
Returns:
np.ndarray: The Hesse matrix.
"""
return get_hesse_matrix(phonopy=self._get_phonopy())
[docs]
def get_band_structure(
self,
npoints: int = 101,
with_eigenvectors: bool = False,
with_group_velocities: bool = False,
):
"""
Get the band structure.
Args:
npoints (int, optional): The number of points. Defaults to 101.
with_eigenvectors (bool, optional): Whether to include eigenvectors. Defaults to False.
with_group_velocities (bool, optional): Whether to include group velocities. Defaults to False.
Returns:
[type]: [description]
"""
return get_band_structure(
phonopy=self._get_phonopy(),
npoints=npoints,
with_eigenvectors=with_eigenvectors,
with_group_velocities=with_group_velocities,
)
[docs]
def plot_band_structure(
self, axis=None, *args, label: Optional[str] = None, **kwargs
):
"""
Plot the band structure.
Args:
axis ([type], optional): The axis. Defaults to None.
label (str, optional): The label. Defaults to None.
Returns:
[type]: [description]
"""
return plot_band_structure(
self._get_phonopy(),
axis,
*args,
label=label,
**kwargs,
)
[docs]
def plot_dos(self, *args, axis=None, **kwargs):
"""
Plot the DOS.
Args:
axis ([type], optional): The axis. Defaults to None.
Returns:
[type]: [description]
"""
return plot_dos(
phonopy_dict=self._phonopy_dict,
*args,
axis=axis,
**kwargs,
)