Source code for atomate2.common.schemas.phonons

"""Schemas for phonon documents."""

import copy
import logging
from pathlib import Path
from typing import Optional, Union

import numpy as np
from emmet.core.math import Matrix3D
from emmet.core.structure import StructureMetadata
from monty.json import MSONable
from phonopy import Phonopy
from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections
from phonopy.structure.symmetry import symmetrize_borns_and_epsilon
from phonopy.units import VaspToTHz
from pydantic import BaseModel, Field
from pymatgen.core import Structure
from pymatgen.io.phonopy import (
    get_ph_bs_symm_line,
    get_ph_dos,
    get_phonopy_structure,
    get_pmg_structure,
)
from pymatgen.io.vasp import Kpoints
from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine
from pymatgen.phonon.dos import PhononDos
from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter
from pymatgen.symmetry.bandstructure import HighSymmKpath
from pymatgen.symmetry.kpath import KPathSeek
from typing_extensions import Self

from atomate2.aims.utils.units import omegaToTHz

logger = logging.getLogger(__name__)


[docs] def get_factor(code: str) -> float: """ Get the frequency conversion factor to THz for each code. Parameters ---------- code: str The code to get the conversion factor for Returns ------- float The correct conversion factor Raises ------ ValueError If code is not defined """ if code in ["forcefields", "vasp"]: return VaspToTHz if code == "aims": return omegaToTHz # Based on CODATA 2002 raise ValueError(f"Frequency conversion factor for code ({code}) not defined.")
[docs] class PhononComputationalSettings(BaseModel): """Collection to store computational settings for the phonon computation.""" # could be optional and implemented at a later stage? npoints_band: int = Field("number of points for band structure computation") kpath_scheme: str = Field("indicates the kpath scheme") kpoint_density_dos: int = Field( "number of points for computation of free energies and densities of states", )
[docs] class ThermalDisplacementData(BaseModel): """Collection to store information on the thermal displacement matrices.""" freq_min_thermal_displacements: float = Field( "cutoff frequency in THz to avoid numerical issues in the " "computation of the thermal displacement parameters" ) thermal_displacement_matrix_cif: Optional[list[list[Matrix3D]]] = Field( None, description="field including thermal displacement matrices in CIF format" ) thermal_displacement_matrix: Optional[list[list[Matrix3D]]] = Field( None, description="field including thermal displacement matrices in Cartesian " "coordinate system", ) temperatures_thermal_displacements: Optional[list[int]] = Field( None, description="temperatures at which the thermal displacement matrices" "have been computed", )
[docs] class PhononUUIDs(BaseModel): """Collection to save all uuids connected to the phonon run.""" optimization_run_uuid: Optional[str] = Field( None, description="optimization run uuid" ) displacements_uuids: Optional[list[str]] = Field( None, description="The uuids of the displacement jobs." ) static_run_uuid: Optional[str] = Field(None, description="static run uuid") born_run_uuid: Optional[str] = Field(None, description="born run uuid")
[docs] class ForceConstants(MSONable): """A force constants class.""" def __init__(self, force_constants: list[list[Matrix3D]]) -> None: self.force_constants = force_constants
[docs] class PhononJobDirs(BaseModel): """Collection to save all job directories relevant for the phonon run.""" displacements_job_dirs: Optional[list[Optional[str]]] = Field( None, description="The directories where the displacement jobs were run." ) static_run_job_dir: Optional[Optional[str]] = Field( None, description="Directory where static run was performed." ) born_run_job_dir: Optional[str] = Field( None, description="Directory where born run was performed." ) optimization_run_job_dir: Optional[str] = Field( None, description="Directory where optimization run was performed." ) taskdoc_run_job_dir: Optional[str] = Field( None, description="Directory where task doc was generated." )
[docs] class PhononBSDOSDoc(StructureMetadata, extra="allow"): # type: ignore[call-arg] """Collection of all data produced by the phonon workflow.""" structure: Optional[Structure] = Field( None, description="Structure of Materials Project." ) phonon_bandstructure: Optional[PhononBandStructureSymmLine] = Field( None, description="Phonon band structure object.", ) phonon_dos: Optional[PhononDos] = Field( None, description="Phonon density of states object.", ) free_energies: Optional[list[float]] = Field( None, description="vibrational part of the free energies in J/mol per " "formula unit for temperatures in temperature_list", ) heat_capacities: Optional[list[float]] = Field( None, description="heat capacities in J/K/mol per " "formula unit for temperatures in temperature_list", ) internal_energies: Optional[list[float]] = Field( None, description="internal energies in J/mol per " "formula unit for temperatures in temperature_list", ) entropies: Optional[list[float]] = Field( None, description="entropies in J/(K*mol) per formula unit" "for temperatures in temperature_list ", ) temperatures: Optional[list[int]] = Field( None, description="temperatures at which the vibrational" " part of the free energies" " and other properties have been computed", ) total_dft_energy: Optional[float] = Field("total DFT energy per formula unit in eV") has_imaginary_modes: Optional[bool] = Field( None, description="if true, structure has imaginary modes" ) # needed, e.g. to compute Grueneisen parameter etc force_constants: Optional[ForceConstants] = Field( None, description="Force constants between every pair of atoms in the structure" ) born: Optional[list[Matrix3D]] = Field( None, description="born charges as computed from phonopy. Only for symmetrically " "different atoms", ) epsilon_static: Optional[Matrix3D] = Field( None, description="The high-frequency dielectric constant" ) supercell_matrix: Matrix3D = Field("matrix describing the supercell") primitive_matrix: Matrix3D = Field( "matrix describing relationship to primitive cell" ) code: str = Field("String describing the code for the computation") phonopy_settings: PhononComputationalSettings = Field( "Field including settings for Phonopy" ) thermal_displacement_data: Optional[ThermalDisplacementData] = Field( "Includes all data of the computation of the thermal displacements" ) jobdirs: Optional[PhononJobDirs] = Field( "Field including all relevant job directories" ) uuids: Optional[PhononUUIDs] = Field("Field including all relevant uuids")
[docs] @classmethod def from_forces_born( cls, structure: Structure, supercell_matrix: np.array, displacement: float, sym_reduce: bool, symprec: float, use_symmetrized_structure: Union[str, None], kpath_scheme: str, code: str, displacement_data: dict[str, list], total_dft_energy: float, epsilon_static: Matrix3D = None, born: Matrix3D = None, **kwargs, ) -> Self: """Generate collection of phonon data. Parameters ---------- structure: Structure object supercell_matrix: numpy array describing the supercell displacement: float size of displacement in angstrom sym_reduce: bool if True, phonopy will use symmetry symprec: float precision to determine kpaths, primitive cells and symmetry in phonopy and pymatgen use_symmetrized_structure: str primitive, conventional or None kpath_scheme: str kpath scheme to generate phonon band structure code: str which code was used for computation displacement_data: output of the displacement data total_dft_energy: float total energy in eV per cell epsilon_static: Matrix3D The high-frequency dielectric constant born: Matrix3D born charges **kwargs: additional arguments """ factor = get_factor(code) # This opens the opportunity to add support for other codes # that are supported by phonopy cell = get_phonopy_structure(structure) if use_symmetrized_structure == "primitive": primitive_matrix: Union[np.ndarray, str] = np.eye(3) else: primitive_matrix = "auto" # TARP: THIS IS BAD! Including for discussions sake if cell.magnetic_moments is not None and primitive_matrix == "auto": if np.any(cell.magnetic_moments != 0.0): raise ValueError( "For materials with magnetic moments, " "use_symmetrized_structure must be 'primitive'" ) cell.magnetic_moments = None phonon = Phonopy( cell, supercell_matrix, primitive_matrix=primitive_matrix, factor=factor, symprec=symprec, is_symmetry=sym_reduce, ) phonon.generate_displacements(distance=displacement) set_of_forces = [np.array(forces) for forces in displacement_data["forces"]] if born is not None and epsilon_static is not None: if len(structure) == len(born): borns, epsilon = symmetrize_borns_and_epsilon( ucell=phonon.unitcell, borns=np.array(born), epsilon=np.array(epsilon_static), symprec=symprec, primitive_matrix=phonon.primitive_matrix, supercell_matrix=phonon.supercell_matrix, is_symmetry=kwargs.get("symmetrize_born", True), ) else: raise ValueError( "Number of born charges does not agree with number of atoms" ) if code == "vasp" and not np.all(np.isclose(borns, 0.0)): phonon.nac_params = { "born": borns, "dielectric": epsilon, "factor": 14.399652, } # Other codes could be added here else: borns = None epsilon = None # Produces all force constants phonon.produce_force_constants(forces=set_of_forces) # with phonopy.load("phonopy.yaml") the phonopy API can be used phonon.save("phonopy.yaml") # get phonon band structure kpath_dict, kpath_concrete = PhononBSDOSDoc.get_kpath( structure=get_pmg_structure(phonon.primitive), kpath_scheme=kpath_scheme, symprec=symprec, ) npoints_band = kwargs.get("npoints_band", 101) qpoints, connections = get_band_qpoints_and_path_connections( kpath_concrete, npoints=kwargs.get("npoints_band", 101) ) # phonon band structures will always be computed filename_band_yaml = "phonon_band_structure.yaml" # TODO: potentially add kwargs to avoid computation of eigenvectors phonon.run_band_structure( qpoints, path_connections=connections, with_eigenvectors=kwargs.get("band_structure_eigenvectors", False), is_band_connection=kwargs.get("band_structure_eigenvectors", False), ) phonon.write_yaml_band_structure(filename=filename_band_yaml) bs_symm_line = get_ph_bs_symm_line( filename_band_yaml, labels_dict=kpath_dict, has_nac=born is not None ) new_plotter = PhononBSPlotter(bs=bs_symm_line) new_plotter.save_plot( filename=kwargs.get("filename_bs", "phonon_band_structure.pdf"), units=kwargs.get("units", "THz"), ) # will determine if imaginary modes are present in the structure imaginary_modes = bs_symm_line.has_imaginary_freq( tol=kwargs.get("tol_imaginary_modes", 1e-5) ) # gets data for visualization on website - yaml is also enough if kwargs.get("band_structure_eigenvectors"): bs_symm_line.write_phononwebsite("phonon_website.json") # get phonon density of states filename_dos_yaml = "phonon_dos.yaml" kpoint_density_dos = kwargs.get("kpoint_density_dos", 7_000) kpoint = Kpoints.automatic_density( structure=get_pmg_structure(phonon.primitive), kppa=kpoint_density_dos, force_gamma=True, ) phonon.run_mesh(kpoint.kpts[0]) phonon.run_total_dos() phonon.write_total_dos(filename=filename_dos_yaml) dos = get_ph_dos(filename_dos_yaml) new_plotter_dos = PhononDosPlotter() new_plotter_dos.add_dos(label="total", dos=dos) new_plotter_dos.save_plot( filename=kwargs.get("filename_dos", "phonon_dos.pdf"), units=kwargs.get("units", "THz"), ) # compute vibrational part of free energies per formula unit temperature_range = np.arange( kwargs.get("tmin", 0), kwargs.get("tmax", 500), kwargs.get("tstep", 10) ) free_energies = [ dos.helmholtz_free_energy( temp=temp, structure=get_pmg_structure(phonon.primitive) ) for temp in temperature_range ] entropies = [ dos.entropy(temp=temp, structure=get_pmg_structure(phonon.primitive)) for temp in temperature_range ] internal_energies = [ dos.internal_energy( temp=temp, structure=get_pmg_structure(phonon.primitive) ) for temp in temperature_range ] heat_capacities = [ dos.cv(temp=temp, structure=get_pmg_structure(phonon.primitive)) for temp in temperature_range ] # will compute thermal displacement matrices # for the primitive cell (phonon.primitive!) # only this is available in phonopy if kwargs.get("create_thermal_displacements"): phonon.run_mesh( kpoint.kpts[0], with_eigenvectors=True, is_mesh_symmetry=False ) freq_min_thermal_displacements = kwargs.get( "freq_min_thermal_displacements", 0.0 ) phonon.run_thermal_displacement_matrices( t_min=kwargs.get("tmin_thermal_displacements", 0), t_max=kwargs.get("tmax_thermal_displacements", 500), t_step=kwargs.get("tstep_thermal_displacements", 100), freq_min=freq_min_thermal_displacements, ) temperature_range_thermal_displacements = np.arange( kwargs.get("tmin_thermal_displacements", 0), kwargs.get("tmax_thermal_displacements", 500), kwargs.get("tstep_thermal_displacements", 100), ) for idx, temp in enumerate(temperature_range_thermal_displacements): phonon.thermal_displacement_matrices.write_cif( phonon.primitive, idx, filename=f"tdispmat_{temp}K.cif" ) _disp_mat = phonon._thermal_displacement_matrices tdisp_mat = _disp_mat.thermal_displacement_matrices.tolist() tdisp_mat_cif = _disp_mat.thermal_displacement_matrices_cif.tolist() else: tdisp_mat = None tdisp_mat_cif = None formula_units = ( structure.composition.num_atoms / structure.composition.reduced_composition.num_atoms ) total_dft_energy_per_formula_unit = ( total_dft_energy / formula_units if total_dft_energy is not None else None ) return cls.from_structure( structure=structure, meta_structure=structure, phonon_bandstructure=bs_symm_line, phonon_dos=dos, free_energies=free_energies, internal_energies=internal_energies, heat_capacities=heat_capacities, entropies=entropies, temperatures=temperature_range.tolist(), total_dft_energy=total_dft_energy_per_formula_unit, has_imaginary_modes=imaginary_modes, force_constants={"force_constants": phonon.force_constants.tolist()} if kwargs["store_force_constants"] else None, born=borns.tolist() if borns is not None else None, epsilon_static=epsilon.tolist() if epsilon is not None else None, supercell_matrix=phonon.supercell_matrix.tolist(), primitive_matrix=phonon.primitive_matrix.tolist(), code=code, thermal_displacement_data={ "temperatures_thermal_displacements": temperature_range_thermal_displacements.tolist(), # noqa: E501 "thermal_displacement_matrix_cif": tdisp_mat_cif, "thermal_displacement_matrix": tdisp_mat, "freq_min_thermal_displacements": freq_min_thermal_displacements, } if kwargs.get("create_thermal_displacements") else None, jobdirs={ "displacements_job_dirs": displacement_data["dirs"], "static_run_job_dir": kwargs["static_run_job_dir"], "born_run_job_dir": kwargs["born_run_job_dir"], "optimization_run_job_dir": kwargs["optimization_run_job_dir"], "taskdoc_run_job_dir": str(Path.cwd()), }, uuids={ "displacements_uuids": displacement_data["uuids"], "born_run_uuid": kwargs["born_run_uuid"], "optimization_run_uuid": kwargs["optimization_run_uuid"], "static_run_uuid": kwargs["static_run_uuid"], }, phonopy_settings={ "npoints_band": npoints_band, "kpath_scheme": kpath_scheme, "kpoint_density_dos": kpoint_density_dos, }, )
[docs] @staticmethod def get_kpath( structure: Structure, kpath_scheme: str, symprec: float, **kpath_kwargs ) -> tuple: """Get high-symmetry points in k-space in phonopy format. Parameters ---------- structure: Structure Object kpath_scheme: str string describing kpath symprec: float precision for symmetry determination **kpath_kwargs: additional parameters that can be passed to this method as a dict """ if kpath_scheme in ("setyawan_curtarolo", "latimer_munro", "hinuma"): high_symm_kpath = HighSymmKpath( structure, path_type=kpath_scheme, symprec=symprec, **kpath_kwargs ) kpath = high_symm_kpath.kpath elif kpath_scheme == "seekpath": high_symm_kpath = KPathSeek(structure, symprec=symprec, **kpath_kwargs) kpath = high_symm_kpath._kpath path = copy.deepcopy(kpath["path"]) for set_idx, label_set in enumerate(kpath["path"]): for lbl_idx, label in enumerate(label_set): path[set_idx][lbl_idx] = kpath["kpoints"][label] return kpath["kpoints"], path