diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index cafdc811bb..6720e8980f 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -67,6 +67,7 @@ jobs: mkdir -p ~/.abinit/pseudos cp -r tests/test_data/abinit/pseudos/ONCVPSP-PBE-SR-PDv0.4 ~/.abinit/pseudos uv pip install .[strict,strict-forcefields,tests,abinit] + micromamba install -n a2 -c conda-forge phono3py --yes # Added uv pip install torch-runstats uv pip install --no-deps nequip==0.5.6 diff --git a/pyproject.toml b/pyproject.toml index d2384e1689..6892c1c2eb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,6 +43,16 @@ amset = ["amset>=0.4.15", "pydash"] cclib = ["cclib"] mp = ["mp-api>=0.37.5"] phonons = ["phonopy>=1.10.8", "seekpath>=2.0.0"] +hiphive = [ + "phonopy==2.21.2", + "hiphive @ git+https://gitlab.com/jsyony37/hiphive.git@personal#egg=hiphive", + "f90nml==1.4.4", + "spglib>=2.2.0", + "ase>=3.22.1", + "phono3py==2.8.0", + "numpy==1.26.4", + "scipy==1.14.1", + ] lobster = ["ijson>=3.2.2", "lobsterpy>=0.4.0"] defects = [ "dscribe>=1.2.0", @@ -117,6 +127,10 @@ strict = [ "seekpath==2.1.0", "tblite==0.3.0; python_version < '3.12'", "typing-extensions==4.12.2", + "hiphive @ git+https://gitlab.com/jsyony37/hiphive.git@personal#egg=hiphive", + "f90nml==1.4.4", + "spglib>=2.2.0", + "scipy==1.14.1", ] strict-forcefields = [ "calorine==3.0", @@ -223,7 +237,7 @@ isort.known-first-party = ["atomate2"] isort.split-on-trailing-comma = false [tool.ruff.format] -docstring-code-format = true +# docstring-code-format = true [tool.ruff.lint.per-file-ignores] "__init__.py" = ["F401"] diff --git a/src/atomate2/common/flows/hiphive.py b/src/atomate2/common/flows/hiphive.py new file mode 100644 index 0000000000..f93316740d --- /dev/null +++ b/src/atomate2/common/flows/hiphive.py @@ -0,0 +1,398 @@ +"""Common Flow for calculating harmonic & anharmonic props of phonon.""" + +# Basic Python packages +from __future__ import annotations + +import logging +from abc import ABC, abstractmethod +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, ClassVar + +from jobflow import Flow, Maker + +from atomate2.common.jobs.hiphive import ( + generate_frequencies_eigenvectors, + generate_phonon_displacements, +) + +# Jobflow packages +from atomate2.common.jobs.phonons import ( + get_supercell_size, + get_total_energy_per_cell, + run_phonon_displacements, +) +from atomate2.common.jobs.utils import structure_to_conventional, structure_to_primitive +from atomate2.forcefields.jobs import ( + CHGNetStaticMaker, + ForceFieldRelaxMaker, + ForceFieldStaticMaker, +) + +# Atomate2 packages +from atomate2.vasp.jobs.phonons import PhononDisplacementMaker +from atomate2.vasp.sets.core import StaticSetGenerator + +if TYPE_CHECKING: + from pathlib import Path + + from emmet.core.math import Matrix3D + from pymatgen.core.structure import Structure # type: ignore # noqa: PGH003 + + from atomate2.vasp.flows.core import DoubleRelaxMaker + from atomate2.vasp.jobs.base import BaseVaspMaker + +SUPPORTED_CODES = ["vasp", "forcefields"] + +logger = logging.getLogger(__name__) + +__all__ = ["BaseHiphiveMaker"] + +__author__ = "Alex Ganose, Junsoo Park, Zhuoying Zhu, Hrushikesh Sahasrabuddhe" +__email__ = "aganose@lbl.gov, jsyony37@lbl.gov, zyzhu@lbl.gov, hpsahasrabuddhe@lbl.gov" + + +@dataclass +class BaseHiphiveMaker(Maker, ABC): + """ + Workflow to calc. interatomic force constants and vibrational props using hiPhive. + + A summary of the workflow is as follows: + 1. Structure relaxtion + 2. Calculate a supercell transformation matrix that brings the + structure as close as cubic as possible, with all lattice lengths + greater than 5 nearest neighbor distances. Then, perturb the atomic sites + for each supercell using a Fixed displacement rattle procedure. The atoms are + perturbed roughly according to a normal deviation around the average value. + A number of standard deviation perturbation distances are included. Multiple + supercells may be generated for each perturbation distance. Then, run static + VASP calculations on each perturbed supercell to calculate atomic forces. + Then, aggregate the forces and the perturbed structures. + 3. Conduct the fit atomic force constants using the regression schemes in hiPhive. + 4. Perform phonon renormalization at finite temperature - useful when unstable + modes exist + 5. Output the interatomic force constants, and phonon band structure and density of + states to the database + 6. Solve the lattice thermal conductivity using ShengBTE and output to the database. + + Args + ---------- + name : str + Name of the flows produced by this maker. + bulk_relax_maker (BaseVaspMaker | None): + The VASP input generator for bulk relaxation, + default is DoubleRelaxMaker using TightRelaxMaker. + phonon_displacement_maker (BaseVaspMaker | None): + The VASP input generator for phonon displacement calculations, + default is PhononDisplacementMaker. + ff_displacement_maker (BaseVaspMaker | None): + The force field displacement maker, default is CHGNetStaticMaker. + min_length (float): + Minimum length of supercell lattice vectors in Angstroms, default is 13.0. + max_length (float): + Minimum length of supercell lattice vectors in Angstroms, default is 25.0. + prefer_90_degrees (bool): + Whether to prefer 90 degree angles in supercell matrix, + default is True. + supercell_matrix_kwargs (dict): + Keyword arguments for supercell matrix calculation, default is {}. + IMAGINARY_TOL (float): + Imaginary frequency tolerance in THz, default is 0.025. + MESH_DENSITY (float): + Mesh density for phonon calculations, default is 100.0. + T_QHA (list): + Temperatures for phonopy thermodynamic calculations, + default is [0, 100, 200, ..., 2000]. + T_RENORM (list): + Temperatures for renormalization calculations, default is [1500]. + T_KLAT (int): + Temperature for lattice thermal conductivity calculation, default is 300. + FIT_METHOD (str): + Method for fitting force constants, default is "rfe". + sym_reduce (bool): + Whether to reduce the number of deformations using symmetry, default is True. + symprec (float): + Symmetry precision to use in the reduction of symmetry to find the + primitive/conventional cell and to handle all symmetry-related tasks in phonopy, + default is 1e-4. + displacement (float): + Displacement distance for phonons, default is 0.01. + get_supercell_size_kwargs (dict): + Keyword arguments for get_supercell_size, default is {}. + use_symmetrized_structure (str | None): + Whether to use a symmetrized structure, default is None. + static_energy_maker (BaseVaspMaker | None): + The VASP input generator for static calculations, default is None. + born_maker (BaseVaspMaker | None): + The VASP input generator for BORN charge calculations, default is None. + create_thermal_displacements (bool): + Whether to create thermal displacements, default is True. + generate_frequencies_eigenvectors_kwargs (dict): + Keyword arguments for generate_frequencies_eigenvectors, default is {}. + kpath_scheme (str): + Scheme to generate kpoints, default is "seekpath". + code (str): + The calculator (VASP, MLFF etc) to use for the calculations, default is None. + store_force_constants (bool): + Whether to store the force constants, default is True. + socket (bool): + Whether to use the socket interface for VASP, default is False. + """ + + name: str = "Lattice-Dynamics" + bulk_relax_maker: DoubleRelaxMaker | ForceFieldRelaxMaker | None = None + phonon_displacement_maker: BaseVaspMaker | ForceFieldStaticMaker | None = field( + default_factory=lambda: PhononDisplacementMaker( + input_set_generator=StaticSetGenerator(auto_lreal=True) + ) + ) + ff_displacement_maker: ForceFieldStaticMaker | None = field( + default_factory=CHGNetStaticMaker + ) + min_length: float | None = 13.0 + max_length: float | None = 25.0 + prefer_90_degrees: bool = True + supercell_matrix_kwargs: dict = field(default_factory=dict) + IMAGINARY_TOL = 0.025 # in THz + MESH_DENSITY = 100.0 # should always be a float + T_QHA: ClassVar[list[int]] = [ + i * 100 for i in range(21) + ] # Temp. for phonopy calc. of thermo. properties (free energy etc.) + FIT_METHOD = "least-squares" # least-squares #omp #rfe #elasticnet + sym_reduce: bool = True + symprec: float = 1e-4 + displacement: float = 0.01 + get_supercell_size_kwargs: dict = field(default_factory=dict) + use_symmetrized_structure: str | None = None + static_energy_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None + born_maker: ForceFieldStaticMaker | BaseVaspMaker | None = None + create_thermal_displacements: bool = True + generate_frequencies_eigenvectors_kwargs: dict = field(default_factory=dict) + kpath_scheme: str = "seekpath" + code: str = None + store_force_constants: bool = True + socket: bool = False + + def make( + self, + mpid: str, + structure: Structure, + bulk_modulus: float, + supercell_matrix: Matrix3D | None = None, + disp_cut: float | None = None, + cutoffs: list[list[float]] | None = None, + prev_dir: str | Path | None = None, + born: list[Matrix3D] | None = None, + epsilon_static: Matrix3D | None = None, + total_dft_energy_per_formula_unit: float | None = None, + ) -> Flow: + """ + Make flow to calculate the harmonic & anharmonic properties of phonon. + + Parameters + ---------- + mpid (str): + The Materials Project ID (MPID) of the material. + structure (Structure): + The A pymatgen structure of the material. + bulk_modulus (float): + Bulk modulus of the material in GPa. + supercell_matrix (list[list[int]], optional): + Supercell transformation matrix, default is None. + disp_cut (float, optional): + Cutoff distance for displacements in Angstroms, default is None. + cutoffs (List[List[float]], optional): + List of cutoff distances for different force constants fitting, + default is None. + prev_dir (str | Path | None, optional): + Previous RELAX calculation directory to use for copying outputs., + default is None. + born (List[Matrix3D], optional): + Born charges for the material, default is None. + epsilon_static (Matrix3D, optional): + Static dielectric tensor for the material, default is None. + total_dft_energy_per_formula_unit (float, optional): + Total DFT energy per formula unit, default is None. + """ + use_symmetrized_structure = self.use_symmetrized_structure + kpath_scheme = self.kpath_scheme + valid_structs = (None, "primitive", "conventional") + if use_symmetrized_structure not in valid_structs: + raise ValueError( + f"Invalid {use_symmetrized_structure=}, use one of {valid_structs}" + ) + + if use_symmetrized_structure != "primitive" and kpath_scheme != "seekpath": + raise ValueError( + f"You can't use {kpath_scheme=} with the primitive standard " + "structure, please use seekpath" + ) + + valid_schemes = ("seekpath", "hinuma", "setyawan_curtarolo", "latimer_munro") + if kpath_scheme not in valid_schemes: + raise ValueError( + f"{kpath_scheme=} is not implemented, use one of {valid_schemes}" + ) + + if self.code is None or self.code not in SUPPORTED_CODES: + raise ValueError( + "The code variable must be passed and it must be a supported code." + f" Supported codes are: {SUPPORTED_CODES}" + ) + + jobs = [] + + # TODO: should this be after or before structural optimization as the + # optimization could change the symmetry we could add a tutorial and point out + # that the structure should be nearly optimized before the phonon workflow + if self.use_symmetrized_structure == "primitive": + # These structures are compatible with many + # of the kpath algorithms that are used for Materials Project + prim_job = structure_to_primitive(structure, self.symprec) + jobs.append(prim_job) + structure = prim_job.output + elif self.use_symmetrized_structure == "conventional": + # it could be beneficial to use conventional standard structures to arrive + # faster at supercells with right angles + conv_job = structure_to_conventional(structure, self.symprec) + jobs.append(conv_job) + structure = conv_job.output + + optimization_run_job_dir = None + optimization_run_uuid = None + + if self.bulk_relax_maker is not None: + # optionally relax the structure + bulk_kwargs = {} + if self.prev_calc_dir_argname is not None: + bulk_kwargs[self.prev_calc_dir_argname] = prev_dir + bulk = self.bulk_relax_maker.make(structure, **bulk_kwargs) + jobs.append(bulk) + structure = bulk.output.structure + prev_dir = bulk.output.dir_name + optimization_run_job_dir = bulk.output.dir_name + optimization_run_uuid = bulk.output.uuid + + # if supercell_matrix is None, supercell size will be determined after relax + # maker to ensure that cell lengths are really larger than threshold + if supercell_matrix is None: + supercell_job = get_supercell_size( + structure, + self.min_length, + self.max_length, + self.prefer_90_degrees, + **self.get_supercell_size_kwargs, + ) + jobs.append(supercell_job) + supercell_matrix = supercell_job.output + + # Computation of static energy + total_dft_energy = None + static_run_job_dir = None + static_run_uuid = None + if (self.static_energy_maker is not None) and ( + total_dft_energy_per_formula_unit is None + ): + static_job_kwargs = {} + if self.prev_calc_dir_argname is not None: + static_job_kwargs[self.prev_calc_dir_argname] = prev_dir + static_job = self.static_energy_maker.make( + structure=structure, **static_job_kwargs + ) + jobs.append(static_job) + total_dft_energy = static_job.output.output.energy + static_run_job_dir = static_job.output.dir_name + static_run_uuid = static_job.output.uuid + prev_dir = static_job.output.dir_name + elif total_dft_energy_per_formula_unit is not None: + # to make sure that one can reuse results from Doc + compute_total_energy_job = get_total_energy_per_cell( + total_dft_energy_per_formula_unit, structure + ) + jobs.append(compute_total_energy_job) + total_dft_energy = compute_total_energy_job.output + + # get a phonon object from phonopy + displacements = generate_phonon_displacements( + structure=structure, + supercell_matrix=supercell_matrix, + fixed_displs=[0.01, 0.03, 0.08, 0.1], + ) + jobs.append(displacements) + + # perform the phonon displacement calculations + displacement_calcs = run_phonon_displacements( + displacements=displacements.output, + structure=structure, + supercell_matrix=supercell_matrix, + phonon_maker=self.phonon_displacement_maker, + socket=self.socket, + prev_dir_argname=self.prev_calc_dir_argname, + prev_dir=prev_dir, + ) + jobs.append(displacement_calcs) + + # Computation of BORN charges + born_run_job_dir = None + born_run_uuid = None + if self.born_maker is not None and (born is None or epsilon_static is None): + born_kwargs = {} + if self.prev_calc_dir_argname is not None: + born_kwargs[self.prev_calc_dir_argname] = prev_dir + born_job = self.born_maker.make(structure, **born_kwargs) + jobs.append(born_job) + + # I am not happy how we currently access "born" charges + # This is very vasp specific code aims and forcefields + # do not support this at the moment, if this changes we have + # to update this section + epsilon_static = born_job.output.calcs_reversed[0].output.epsilon_static + born = born_job.output.calcs_reversed[0].output.outcar["born"] + born_run_job_dir = born_job.output.dir_name + born_run_uuid = born_job.output.uuid + + logger.info("Generating phonon frequencies and eigenvectors") + phonon_collect = generate_frequencies_eigenvectors( + supercell_matrix=supercell_matrix, + displacement=self.displacement, + sym_reduce=self.sym_reduce, + symprec=self.symprec, + use_symmetrized_structure=self.use_symmetrized_structure, + kpath_scheme=self.kpath_scheme, + code=self.code, + structure=structure, + displacement_data=displacement_calcs.output, + epsilon_static=epsilon_static, + born=born, + total_dft_energy=total_dft_energy, + static_run_job_dir=static_run_job_dir, + static_run_uuid=static_run_uuid, + born_run_job_dir=born_run_job_dir, + born_run_uuid=born_run_uuid, + optimization_run_job_dir=optimization_run_job_dir, + optimization_run_uuid=optimization_run_uuid, + create_thermal_displacements=self.create_thermal_displacements, + store_force_constants=self.store_force_constants, + bulk_modulus=bulk_modulus, + disp_cut=disp_cut, + **self.generate_frequencies_eigenvectors_kwargs, + ) + + jobs.append(phonon_collect) + + return Flow( + jobs=jobs, + output=phonon_collect.output, + name=f"{mpid}_" f"{disp_cut}_" f"{cutoffs}_" f"{self.name}", + ) + + @property + @abstractmethod + def prev_calc_dir_argname(self) -> str | None: + """Name of argument informing static maker of previous calculation directory. + + As this differs between different DFT codes (e.g., VASP, CP2K), it + has been left as a property to be implemented by the inheriting class. + + Note: this is only applicable if a relax_maker is specified; i.e., two + calculations are performed for each ordering (relax -> static) + """ diff --git a/src/atomate2/common/jobs/hiphive.py b/src/atomate2/common/jobs/hiphive.py new file mode 100644 index 0000000000..7076a105e2 --- /dev/null +++ b/src/atomate2/common/jobs/hiphive.py @@ -0,0 +1,243 @@ +"""Jobs for calculating harmonic & anharmonic props of phonon using hiPhive.""" + +# Basic Python packages +from __future__ import annotations + +import warnings +from typing import TYPE_CHECKING + +import numpy as np +import scipy as sp + +# Jobflow packages +from jobflow import job +from monty.serialization import dumpfn + +# Phonopy & Phono3py +from pymatgen.core.structure import Structure +from pymatgen.io.ase import AseAtomsAdaptor +from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine +from pymatgen.phonon.dos import PhononDos +from pymatgen.transformations.standard_transformations import SupercellTransformation + +from atomate2.common.schemas.hiphive import ForceConstants, PhononBSDOSDoc +from atomate2.utils.log import initialize_logger + +warnings.filterwarnings("ignore", module="pymatgen") +warnings.filterwarnings("ignore", module="ase") + +if TYPE_CHECKING: + from emmet.core.math import Matrix3D + + +logger = initialize_logger(level=3) + +IMAGINARY_TOL = 0.1 # in THz # changed from 0.025 +FIT_METHOD = "rfe" + +ev2j = sp.constants.elementary_charge +hbar = sp.constants.hbar # J-s +kb = sp.constants.Boltzmann # J/K + +__all__ = [ + "generate_phonon_displacements", + "generate_frequencies_eigenvectors", +] + +__author__ = "Alex Ganose, Junsoo Park, Zhuoying Zhu, Hrushikesh Sahasrabuddhe" +__email__ = "aganose@lbl.gov, jsyony37@lbl.gov, zyzhu@lbl.gov, hpsahasrabuddhe@lbl.gov" + + +@job(data=[Structure]) +def generate_phonon_displacements( + structure: Structure, + supercell_matrix: np.array, + fixed_displs: list[float], + nconfigs: int = 1, +) -> list[Structure]: + """ + Generate displaced structures with phonopy. + + Parameters + ---------- + structure: Structure object + Fully optimized input structure for phonon run + supercell_matrix: np.array + array to describe supercell matrix + displacement: float + displacement in Angstrom + sym_reduce: bool + if True, symmetry will be used to generate displacements + symprec: float + precision to determine symmetry + use_symmetrized_structure: str or None + primitive, conventional or None + kpath_scheme: str + scheme to generate kpath + code: + code to perform the computations + """ + logger.info(f"supercell_matrix = {supercell_matrix}") + supercell_structure = SupercellTransformation( + scaling_matrix=supercell_matrix + ).apply_transformation(structure) + logger.info(f"supercell_structure = {supercell_structure}") + structure_data = { + "structure": structure, + "supercell_structure": supercell_structure, + "supercell_matrix": supercell_matrix, + } + + dumpfn(structure_data, "structure_data.json") + + rattled_structures = [ + supercell_structure for _ in range(nconfigs * len(fixed_displs)) + ] + + atoms_list = [] + for idx, rattled_structure in enumerate(rattled_structures): + logger.info(f"iter number = {idx}") + atoms = AseAtomsAdaptor.get_atoms(rattled_structure) + atoms_tmp = atoms.copy() + + # number of unique displs + n_displs = len(fixed_displs) + + # map the index to distance from fixed_displs list + distance_mapping = {i: fixed_displs[i] for i in range(n_displs)} + + # Calculate the distance based on the index pattern + distance = distance_mapping[idx % n_displs] + logger.info(f"distance_{idx % n_displs} = {distance}") + + # set the random seed for reproducibility + # 6 is the number of fixed_displs + rng = np.random.default_rng(seed=(idx // n_displs)) # idx // 6 % 6 + + total_inds = list(enumerate(atoms)) + logger.info(f"total_inds = {total_inds}") + logger.info(f"len(total_inds) = {len(total_inds)}") + + def generate_normal_displacement( + distance: float, n: int, rng: np.random.Generator + ) -> np.ndarray: + directions = rng.normal(size=(n, 3)) + normalizer = np.linalg.norm(directions, axis=1, keepdims=True) + distance_normal_distribution = rng.normal( + distance, distance / 5, size=(n, 1) + ) + displacements = distance_normal_distribution * directions / normalizer + logger.info(f"displacements = {displacements}") + return displacements + + # Generate displacements + disp_normal = generate_normal_displacement( + distance, len(total_inds), rng + ) # Uncomment this later + mean_displacements = np.linalg.norm( + disp_normal, axis=1 + ).mean() # Uncomment this later + + logger.info(f"mean_displacements = {mean_displacements}") + + atoms_tmp = atoms.copy() + + # add the disp_normal to the all the atoms in the structure + for i in range(len(total_inds)): + atoms_tmp.positions[i] += disp_normal[i] + + atoms_list.append(atoms_tmp) + + # Convert back to pymatgen structure + structures_pymatgen = [] + for atoms_ase in range(len(atoms_list)): + logger.info(f"atoms: {atoms_ase}") + logger.info(f"structures_ase_all[atoms]: {atoms_list[atoms_ase]}") + structure_i = AseAtomsAdaptor.get_structure(atoms_list[atoms_ase]) + structures_pymatgen.append(structure_i) + + for i in range(len(structures_pymatgen)): + structures_pymatgen[i].to(f"POSCAR_{i}", "poscar") + + dumpfn(structures_pymatgen, "perturbed_structures.json") + + return structures_pymatgen + + +@job( + output_schema=PhononBSDOSDoc, + data=[PhononDos, PhononBandStructureSymmLine, ForceConstants], +) +def generate_frequencies_eigenvectors( + structure: Structure, + supercell_matrix: np.array, + displacement: float, + sym_reduce: bool, + symprec: float, + use_symmetrized_structure: str | None, + kpath_scheme: str, + code: str, + displacement_data: dict[str, list], + total_dft_energy: float, + epsilon_static: Matrix3D = None, + born: Matrix3D = None, + bulk_modulus: float = None, + disp_cut: float = 0.05, + **kwargs, +) -> PhononBSDOSDoc: + """ + Analyze the phonon runs and summarize the results. + + Parameters + ---------- + structure: Structure object + Fully optimized structure used for phonon runs + supercell_matrix: np.array + array to describe supercell + displacement: float + displacement in Angstrom used for supercell computation + sym_reduce: bool + if True, symmetry will be used in phonopy + symprec: float + precision to determine symmetry + use_symmetrized_structure: str + primitive, conventional, None are allowed + kpath_scheme: str + kpath scheme for phonon band structure computation + code: str + code to run computations + displacement_data: dict + outputs from displacements + total_dft_energy: float + total DFT energy in eV per cell + epsilon_static: Matrix3D + The high-frequency dielectric constant + born: Matrix3D + Born charges + bulk_modulus: float + Bulk modulus in GPa + disp_cut: float + Displacement cut-off in Angstrom + kwargs: dict + Additional parameters that are passed to PhononBSDOSDoc.from_forces_born + """ + logger.info("Starting generate_frequencies_eigenvectors()") + return PhononBSDOSDoc.from_forces_born( + structure=structure.remove_site_property(property_name="magmom") + if "magmom" in structure.site_properties + else structure, + supercell_matrix=supercell_matrix, + displacement=displacement, + sym_reduce=sym_reduce, + symprec=symprec, + use_symmetrized_structure=use_symmetrized_structure, + kpath_scheme=kpath_scheme, + code=code, + displacement_data=displacement_data, + total_dft_energy=total_dft_energy, + epsilon_static=epsilon_static, + born=born, + bulk_modulus=bulk_modulus, + disp_cut=disp_cut, + **kwargs, + ) diff --git a/src/atomate2/common/jobs/phonons.py b/src/atomate2/common/jobs/phonons.py index 39457c8c9f..a12fd7359f 100644 --- a/src/atomate2/common/jobs/phonons.py +++ b/src/atomate2/common/jobs/phonons.py @@ -195,6 +195,15 @@ def generate_phonon_displacements( supercells = phonon.supercells_with_displacements + from monty.serialization import dumpfn + + # Convert back to pymatgen structure + structures_pymatgen = [get_pmg_structure(cell) for cell in supercells] + for i in range(len(structures_pymatgen)): + structures_pymatgen[i].to(f"POSCAR_{i}", "poscar") + + dumpfn(structures_pymatgen, "perturbed_structures.json") + return [get_pmg_structure(cell) for cell in supercells] @@ -307,6 +316,7 @@ def run_phonon_displacements( "forces": [], "uuids": [], "dirs": [], + "structure": [], } phonon_job_kwargs = {} if prev_dir is not None and prev_dir_argname is not None: @@ -322,11 +332,13 @@ def run_phonon_displacements( phonon_job.update_maker_kwargs( {"_set": {"write_additional_data->phonon_info:json": info}}, dict_mod=True ) + phonon_jobs.append(phonon_job) outputs["displacement_number"] = list(range(len(displacements))) outputs["uuids"] = [phonon_job.output.uuid] * len(displacements) outputs["dirs"] = [phonon_job.output.dir_name] * len(displacements) outputs["forces"] = phonon_job.output.output.all_forces + outputs["structure"].append(phonon_job.output.output.structure) else: for idx, displacement in enumerate(displacements): if prev_dir is not None: @@ -352,6 +364,7 @@ def run_phonon_displacements( outputs["uuids"].append(phonon_job.output.uuid) outputs["dirs"].append(phonon_job.output.dir_name) outputs["forces"].append(phonon_job.output.output.forces) + outputs["structure"].append(phonon_job.output.output.structure) displacement_flow = Flow(phonon_jobs, outputs) return Response(replace=displacement_flow) diff --git a/src/atomate2/common/schemas/hiphive.py b/src/atomate2/common/schemas/hiphive.py new file mode 100644 index 0000000000..14e22fff12 --- /dev/null +++ b/src/atomate2/common/schemas/hiphive.py @@ -0,0 +1,2130 @@ +"""Schemas for phonon documents.""" + +from __future__ import annotations + +import contextlib +import copy +import json +import logging +import os +import warnings +from itertools import product +from pathlib import Path +from typing import TYPE_CHECKING, Any, Optional, Union + +import matplotlib.pyplot as plt +import numpy as np +import psutil +import scipy as sp +from emmet.core.math import Matrix3D +from emmet.core.structure import StructureMetadata + +# Hiphive packages +from hiphive import ( + ClusterSpace, + ForceConstantPotential, + StructureContainer, + enforce_rotational_sum_rules, +) +from hiphive import ForceConstants as HiphiveForceConstants +from hiphive.cutoffs import estimate_maximum_cutoff, is_cutoff_allowed +from hiphive.fitting import Optimizer +from hiphive.utilities import get_displacements +from monty.json import MSONable +from monty.serialization import dumpfn +from phono3py.phonon3.gruneisen import Gruneisen + +# Phonopy & Phono3py +from phonopy import Phonopy +from phonopy.file_IO import parse_FORCE_CONSTANTS +from phonopy.interface.hiphive_interface import phonopy_atoms_to_ase +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 + +# Pymatgen packages +from pymatgen.core import Structure +from pymatgen.io.ase import AseAtomsAdaptor +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.analyzer import SpacegroupAnalyzer +from pymatgen.symmetry.bandstructure import HighSymmKpath +from pymatgen.symmetry.kpath import KPathSeek +from pymatgen.transformations.standard_transformations import SupercellTransformation +from typing_extensions import Self + +from atomate2.aims.utils.units import omegaToTHz + +warnings.filterwarnings("ignore", module="pymatgen") +warnings.filterwarnings("ignore", module="ase") + +if TYPE_CHECKING: + from ase.atoms import Atoms + from emmet.core.math import Matrix3D + from pymatgen.core.structure import Structure + + +logger = logging.getLogger(__name__) + +ev2j = sp.constants.elementary_charge +hbar = sp.constants.hbar # J-s +kb = sp.constants.Boltzmann # J/K + + +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.") + + +def get_cutoffs(supercell_structure: Structure) -> list[list[float]]: + """ + Determine the trial cutoffs based on a supercell structure. + + This function calculates and returns the best cutoffs for 2nd, 3rd, and 4th order + interactions for a given supercell structure. It performs linear interpolation to + find cutoffs for specific target degrees of freedom (DOFs), generates combinations + of cutoffs, and filters them to retain unique DOFs. + + Args: + supercell_structure: A structure. + + Returns + ------- + list[list[float]]: A list of lists where each sublist contains the best cutoffs + for 2nd, 3rd, and 4th order interactions. + """ + # Initialize lists and variables + cutoffs_2nd_list = [] + cutoffs_3rd_list = [] + cutoffs_4th_list = [] + n_doffs_2nd_list: list[int] = [] + n_doffs_3rd_list: list[int] = [] + n_doffs_4th_list: list[int] = [] + # create a list of cutoffs to check starting from 2 to 11, with a step size of 0.1 + supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + max_cutoff = estimate_maximum_cutoff(supercell_atoms) + cutoffs_2nd_to_check = np.arange(3, max_cutoff, 0.3) + cutoffs_3rd_to_check = np.arange(3, max_cutoff, 0.1) + cutoffs_4th_to_check = np.arange(3, max_cutoff, 0.1) + + # Assume that supercell_atoms is defined before + supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + + def get_best_cutoff( + cutoffs_to_check: np.ndarray, target_n_doff: int, order: int + ) -> tuple[float, list[int], np.ndarray]: + """ + Find the best cutoff value for a given order of interaction. + + This function iterates over a range of cutoff values, evaluates the DOFs, and + finds the best cutoff that meets the target number of DOFs. It performs + iterative refinement to narrow down the best cutoff. + + Args: + cutoffs_to_check (numpy.ndarray): An array of cutoff values to check. + target_n_doff (int): The target number of degrees of freedom (DOFs). + order (int): The order of interaction (2 for 2nd order, 3 for 3rd order, + 4 for 4th order). + + Returns + ------- + tuple: A tuple containing: + - best_cutoff (float): The best cutoff value. + - n_doffs_list (list[int]): A list of DOFs corresponding to the cutoffs + checked. + - cutoffs_to_check (numpy.ndarray): The array of cutoff values checked. + + """ + n_doffs_list: list[int] = [] + best_cutoff = np.inf + + if order == 2: + increment_size = 0.2 + elif order == 3: + increment_size = 0.25 + elif order == 4: + increment_size = 0.05 + + for cutoff in cutoffs_to_check: + # check if order == 2 and if no element in n_doffs_list is greater than 1000 + if order == 2 and all(n_doff < 1000 for n_doff in n_doffs_list): + with contextlib.suppress(ValueError): + cs = ClusterSpace( + supercell_atoms, [cutoff], symprec=1e-3, acoustic_sum_rules=True + ) + elif order == 3 and all(n_doff < 2200 for n_doff in n_doffs_list): + with contextlib.suppress(ValueError): + cs = ClusterSpace( + supercell_atoms, + [2, cutoff], + symprec=1e-3, + acoustic_sum_rules=True, + ) + elif order == 4 and all(n_doff < 2200 for n_doff in n_doffs_list): + with contextlib.suppress(ValueError): + cs = ClusterSpace( + supercell_atoms, + [2, 3, cutoff], + symprec=1e-3, + acoustic_sum_rules=True, + ) + + try: + n_doff = cs.get_n_dofs_by_order(order) + if ( + (order == 2 and n_doff < 1000) + or (order == 3 and n_doff < 2200) + or (order == 4 and n_doff < 2200) + ): + logger.info(f"adding n_doff = {n_doff} to the n_doff list") + n_doffs_list.append(n_doff) + elif ( + (order == 2 and n_doff > 1000) + or (order == 3 and n_doff > 2200) + or (order == 4 and n_doff > 2200) + ): + # remove all the cutoffs from cutoffs_to_check from the current + # cutoff once we are inside this block + cutoff_index = np.where(cutoffs_to_check == cutoff)[0][0] + cutoffs_to_check = cutoffs_to_check[:cutoff_index] + # do the same for n_doffs_list + n_doffs_list = n_doffs_list[: cutoff_index + 1] + break + except UnboundLocalError: + logger.info(f"UnboundLocalError for cutoff = {cutoff}") + # find the index of the cutoff in the cutoffs_to_check list + cutoff_index = np.where(cutoffs_to_check == cutoff)[0][0] + logger.info(cutoff_index) + # remove only the cutoff corresponding to the index from the list + cutoffs_to_check = np.delete(cutoffs_to_check, cutoff_index) + continue + + # Find the closest cutoff to the target + closest_index = np.argmin(np.abs(np.array(n_doffs_list) - target_n_doff)) + best_cutoff = cutoffs_to_check[closest_index] + + while increment_size >= 0.01: + if order == 2: + with contextlib.suppress(ValueError): + cs = ClusterSpace( + supercell_atoms, + [best_cutoff], + symprec=1e-3, + acoustic_sum_rules=True, + ) + elif order == 3: + with contextlib.suppress(ValueError): + cs = ClusterSpace( + supercell_atoms, + [3, best_cutoff], + symprec=1e-3, + acoustic_sum_rules=True, + ) + elif order == 4: + with contextlib.suppress(ValueError): + cs = ClusterSpace( + supercell_atoms, + [3, 3, best_cutoff], + symprec=1e-3, + acoustic_sum_rules=True, + ) + + n_doff = cs.get_n_dofs_by_order(order) + + if n_doff > target_n_doff: + best_cutoff -= increment_size / 2 + else: + best_cutoff += increment_size / 2 + + increment_size /= 4 + + return best_cutoff, n_doffs_list, cutoffs_to_check + + # Get best cutoffs for 2nd, 3rd, and 4th order + logger.info("Getting best cutoffs for 2nd order") + best_2nd_order_cutoff, n_doffs_2nd_list, cutoffs_2nd_to_check = get_best_cutoff( + cutoffs_2nd_to_check, 100, 2 + ) + logger.info("Getting best cutoffs for 3rd order") + best_3rd_order_cutoff, n_doffs_3rd_list, cutoffs_3rd_to_check = get_best_cutoff( + cutoffs_3rd_to_check, 1000, 3 + ) + logger.info("Getting best cutoffs for 4th order") + best_4th_order_cutoff, n_doffs_4th_list, cutoffs_4th_to_check = get_best_cutoff( + cutoffs_4th_to_check, 1000, 4 + ) + + cutoffs_2nd_list.append(best_2nd_order_cutoff) + cutoffs_3rd_list.append(best_3rd_order_cutoff) + cutoffs_4th_list.append(best_4th_order_cutoff) + + best_cutoff_list = [ + [best_2nd_order_cutoff, best_3rd_order_cutoff, best_4th_order_cutoff] + ] + logger.info(f"best_cutoff_list = {best_cutoff_list}") + + # Linear interpolation to find cutoffs for targets + def interpolate_and_find_cutoffs( + cutoffs_to_check: np.ndarray, n_doffs_list: list[int], targets: list[int] + ) -> np.ndarray: + """ + Perform linear interpolation to find cutoff values for specific target DOFs. + + This function interpolates the given cutoff values and their corresponding DOFs + to find the cutoff values that achieve the specified target DOFs. + + Args: + cutoffs_to_check (numpy.ndarray): An array of cutoff values to check. + n_doffs_list (list[int]): A list of DOFs corresponding to the cutoffs checked. + targets (list[int]): A list of target DOFs for which to find the cutoff values. + + Returns + ------- + numpy.ndarray: An array of interpolated cutoff values corresponding to the + target DOFs. + """ + logger.info(f"cutoffs_to_check = {cutoffs_to_check}") + logger.info(f"n_doffs_list = {n_doffs_list}") + logger.info(f"targets = {targets}") + logger.info(f"len(cutoffs_to_check) = {len(cutoffs_to_check)}") + logger.info(f"len(n_doffs_list) = {len(n_doffs_list)}") + cutoffs_for_targets = np.interp(targets, n_doffs_list, cutoffs_to_check) + logger.info(f"cutoffs_for_targets = {cutoffs_for_targets}") + return cutoffs_for_targets + + cutoffs_2nd_list.extend( + interpolate_and_find_cutoffs(cutoffs_2nd_to_check, n_doffs_2nd_list, [70, 220]) + ) + cutoffs_3rd_list.extend( + interpolate_and_find_cutoffs( + cutoffs_3rd_to_check, n_doffs_3rd_list, [1500, 2000] + ) + ) + cutoffs_4th_list.extend( + interpolate_and_find_cutoffs( + cutoffs_4th_to_check, n_doffs_4th_list, [1500, 2000] + ) + ) + + # Sort the lists + cutoffs_2nd_list = sorted(cutoffs_2nd_list) + cutoffs_3rd_list = sorted(cutoffs_3rd_list) + cutoffs_4th_list = sorted(cutoffs_4th_list) + + # Generate combinations of cutoffs + cutoffs = np.array( + list(map(list, product(cutoffs_2nd_list, cutoffs_3rd_list, cutoffs_4th_list))) + ) + logger.info(f"cutoffs = {cutoffs}") + logger.info(f"len(cutoffs) = {len(cutoffs)}") + logger.info(f"cutoffs[0] = {cutoffs[0]}") + dofs_list = [] + for cutoff_value in cutoffs: + logger.info(f"cutoff inside the for loop = {cutoff_value}") + cutoff = cutoff_value.tolist() + logger.info(f"cutoff.tolist() = {cutoff}") + with contextlib.suppress(ValueError): + cs = ClusterSpace( + supercell_atoms, cutoff, symprec=1e-3, acoustic_sum_rules=True + ) + n_doff_2 = cs.get_n_dofs_by_order(2) + n_doff_3 = cs.get_n_dofs_by_order(3) + n_doff_4 = cs.get_n_dofs_by_order(4) + dofs_list.append([n_doff_2, n_doff_3, n_doff_4]) + logger.info(f"dofs_list = {dofs_list}") + + # Save the plots for cutoffs vs n_doffs + def save_plot( + cutoffs_to_check: np.ndarray, n_doffs_list: list[int], order: int + ) -> None: + """Save the plot for cutoffs vs n_doffs.""" + plt.figure() + plt.scatter(cutoffs_to_check, n_doffs_list, color="blue", label="Data Points") + plt.plot( + cutoffs_to_check, n_doffs_list, color="red", label="Linear Interpolation" + ) + plt.xlabel(f"Cutoffs {order} Order") + plt.ylabel(f"n_doffs_{order}") + plt.title(f"Linear Interpolation for n_doffs_{order} vs Cutoffs {order} Order") + plt.legend() + plt.savefig(f"cutoffs_{order}_vs_n_doffs_{order}.png") + + save_plot(cutoffs_2nd_to_check, n_doffs_2nd_list, 2) + save_plot(cutoffs_3rd_to_check, n_doffs_3rd_list, 3) + save_plot(cutoffs_4th_to_check, n_doffs_4th_list, 4) + + logger.info("We have completed the cutoffs calculation.") + logger.info(f"cutoffs = {cutoffs}") + + max_cutoff = estimate_maximum_cutoff(AseAtomsAdaptor.get_atoms(supercell_structure)) + cutoffs[cutoffs > max_cutoff] = max_cutoff + logger.info(f"CUTOFFS \n {cutoffs}") + logger.info(f"MAX_CUTOFF \n {max_cutoff}") + good_cutoffs = np.all(cutoffs < max_cutoff - 0.1, axis=1) + logger.info(f"GOOD CUTOFFS \n{good_cutoffs}") + + cutoffs = cutoffs[good_cutoffs].tolist() + logger.info(f"cutoffs_used = {cutoffs}") + + def filter_cutoffs( + cutoffs: list[list[int]], dofs_list: list[list[int]] + ) -> list[list[int]]: + """Filter cutoffs based on unique DOFs.""" + # Map cutoffs to dofs_list + cutoffs_to_dofs = { + tuple(c): tuple(d) for c, d in zip(cutoffs, dofs_list, strict=False) + } + logger.info(f"Cutoffs to DOFs mapping: {cutoffs_to_dofs}") + + # Track seen dofs and keep only the first occurrence of each unique dofs + seen_dofs = set() + new_cutoffs = [] + + for c, d in zip(cutoffs, dofs_list, strict=False): + d_tuple = tuple(d) + if d_tuple not in seen_dofs: + seen_dofs.add(d_tuple) + new_cutoffs.append(c) + + logger.info(f"Unique DOFs: {seen_dofs}") + logger.info(f"New cutoffs: {new_cutoffs}") + + return new_cutoffs + + cutoffs = filter_cutoffs(cutoffs, dofs_list) + logger.info(f"Filtered cutoffs: {cutoffs}") + + # Round each order of the cutoff to the first decimal place + rounded_cutoffs = [[round(order, 1) for order in cutoff] for cutoff in cutoffs] + logger.info(f"Rounded cutoffs: {rounded_cutoffs}") + + return rounded_cutoffs + + +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", + ) + + +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", + ) + + +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") + + +class ForceConstants(MSONable): + """A force constants class.""" + + def __init__(self, force_constants: list[list[Matrix3D]]) -> None: + self.force_constants = force_constants + + +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." + ) + + +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") + + @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, + bulk_modulus: float, + epsilon_static: Matrix3D = None, + born: Matrix3D = None, + fit_method: str | None = "rfe", + disp_cut: float | None = None, + temperature_qha: float | list | dict | None = None, + imaginary_tol: float = 0.025, # in THz + cutoffs: list[list[float]] | None = 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 + """ + logger.info("Starting from_forces_born.") + factor = get_factor(code) + # This opens the opportunity to add support for other codes + # that are supported by phonopy + + if temperature_qha is None: + temperature_qha = [i * 100 for i in range(21)] + + 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) + # add the def run_hiphive() here and get the 2nd order FCs out. + # then use JZ's way to add FC to the phonon object. + + logger.info("Starting run_hiphive.") + # 3. Hiphive Fitting of FCPs upto 4th order + PhononBSDOSDoc.run_hiphive( + parent_structure=structure, + fit_method=fit_method, + disp_cut=disp_cut, + bulk_modulus=bulk_modulus, + temperature_qha=temperature_qha, + imaginary_tol=imaginary_tol, + # prev_dir_json_saver=static_calcs.output["current_dir"], + cutoffs=cutoffs, + displacement_data=displacement_data, + supercell_matrix=supercell_matrix, + t_qha=[ + i * 100 for i in range(21) + ], # Temp. for phonopy calc. of thermo. properties (free energy etc.) + ) + logger.info("Completed run_hiphive.") + + # Read the force constants from the output file of pheasy code + force_constants = parse_FORCE_CONSTANTS( + filename="FORCE_CONSTANTS_2ND" + ) # FORCE_CONSTANTS_2ND FORCE_CONSTANTS + phonon.force_constants = force_constants + # symmetrize the force constants to make them physically correct based + # on the space group symmetry of the crystal structure. + phonon.symmetrize_force_constants() + + # 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", 1000), 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 # noqa: SLF001 + 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 + ) + + logger.info("Finished from_forces_born.") + + 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, + }, + ) + + @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 # noqa: SLF001 + else: + raise ValueError(f"Unexpected {kpath_scheme=}") + + 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 + + @staticmethod + def run_hiphive( + parent_structure: Structure, + displacement_data: dict[str, list], + cutoffs: list[list] | None = None, + fit_method: str | None = "rfe", + disp_cut: float | None = None, + bulk_modulus: float | None = None, + temperature_qha: float | list | dict | None = None, + imaginary_tol: float | None = None, + # prev_dir_json_saver: str | None = None, + supercell_structure: Structure | None = None, + supercell_matrix: np.array | None = None, + t_qha: float | list | dict | None = None, + ) -> dict: + """ + Fit force constants using hiPhive. + + Requires "perturbed_structures.json", "perturbed_forces.json", and + "structure_data.json" files to be present in the current working directory. + + Args: + cutoffs (Optional[list[list]]): A list of cutoffs to trial. If None, + a set of trial cutoffs will be generated based on the structure + (default). + separate_fit: If True, harmonic and anharmonic force constants are fit + separately and sequentially, harmonic first then anharmonic. If + False, then they are all fit in one go. Default is False. + disp_cut: if separate_fit=True, determines the mean displacement of + perturbed structure to be included in harmonic (<) or + anharmonic (>) fitting + imaginary_tol (float): Tolerance used to decide if a phonon mode + is imaginary, in THz. + fit_method (str): Method used for fitting force constants. This can + be any of the values allowed by the hiphive ``Optimizer`` class. + """ + logger.info(f"cutoffs = {cutoffs}") + logger.info(f"disp_cut is {disp_cut}") + logger.info(f"fit_method is {fit_method}") + + supercell_structure = SupercellTransformation( + scaling_matrix=supercell_matrix + ).apply_transformation(parent_structure) + + if cutoffs is None: + # cutoffs = get_cutoffs(supercell_structure) + cutoffs = [[4, 3, 2]] # Si + logger.info(f"cutoffs is {cutoffs}") + else: + pass + + t_qha = temperature_qha if temperature_qha is not None else t_qha + if isinstance(t_qha, list): + t_qha.sort() + else: + # Handle the case where T_qha is not a list + # You can choose to raise an exception or handle it differently + # For example, if T_qha is a single temperature value, you can + # convert it to a list + t_qha = [t_qha] + + logger.info(f"t_qha is {t_qha}") + + structures = [] + logger.info(f"supercell_structure is {supercell_structure}") + + set_of_forces = [np.array(forces) for forces in displacement_data["forces"]] + supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + for structure, forces in zip( + displacement_data["structure"], set_of_forces, strict=False + ): + logger.info(f"structure is {structure}") + atoms = AseAtomsAdaptor.get_atoms(structure) + displacements = get_displacements(atoms, supercell_atoms) + atoms.new_array("displacements", displacements) + atoms.new_array("forces", forces) + atoms.positions = supercell_atoms.get_positions() + structures.append(atoms) + + # Calculate mean displacements + mean_displacements = np.linalg.norm(displacements, axis=1).mean() + logger.info( + f"Mean displacements while reading individual displacements: " + f"{mean_displacements}" + ) + # Calculate standard deviation of displacements + std_displacements = np.linalg.norm(displacements, axis=1).std() + logger.info( + "Standard deviation of displacements while reading individual" + "displacements: " + f"{std_displacements}" + ) + + all_cutoffs = cutoffs + logger.info(f"all_cutoffs is {all_cutoffs}") + + fcs, param, cs, fitting_data, fcp = PhononBSDOSDoc.fit_force_constants( + parent_structure=parent_structure, + supercell_matrix=supercell_matrix, + supercell_structure=supercell_structure, + structures=structures, + all_cutoffs=all_cutoffs, + disp_cut=disp_cut, + imaginary_tol=imaginary_tol, + fit_method=fit_method, + ) + + logger.info("Saving Harmonic props") + thermal_data, phonopy = PhononBSDOSDoc.harmonic_properties( + parent_structure, supercell_matrix, fcs, t_qha, imaginary_tol + ) + + anharmonic_data = PhononBSDOSDoc.anharmonic_properties( + phonopy, + fcs, + t_qha, + thermal_data["heat_capacity"], + thermal_data["n_imaginary"], + bulk_modulus, + ) + + if fcs is None: + raise RuntimeError("Could not find a force constant solution") + + if isinstance(fcs, ForceConstants): + logger.info("Writing force_constants") + fcs.write("force_constants.fcs") + else: + logger.info("fcs is not an instance of ForceConstants") + + if isinstance(fcp, ForceConstantPotential): + logger.info("Writing force_constants_potential") + fcp.write("force_constants_potential.fcp") + + logger.info("Saving parameters") + np.savetxt("parameters.txt", param) + + if isinstance(cs, ClusterSpace): + logger.info("Writing cluster_space") + cs.write("cluster_space.cs") + logger.info("cluster_space writing is complete") + else: + logger.info("cs is not an instance of ClusterSpace") + + logger.info("Saving phonopy_params") + phonopy.save("phonopy_params.yaml") + fitting_data["best_n_imaginary"] = thermal_data.pop("n_imaginary") + thermal_data.update(anharmonic_data) + logger.info("Saving fitting_data") + dumpfn(fitting_data, "fitting_data.json") + logger.info("Saving thermal_data") + dumpfn(thermal_data, "thermal_data.json") + + logger.info("Writing cluster space and force_constants") + logger.info(f"{type(fcs)}") + + logger.info( + f"phonopy_atoms_to_ase(phonopy.supercell) = " + f"{phonopy_atoms_to_ase(phonopy.supercell)}" + ) + logger.info(f"supercell_atoms = {supercell_atoms}") + from pymatgen.io.ase import MSONAtoms + + structure_data_phonopy_atoms = { + "supercell_structure_phonopy_atoms": MSONAtoms( + phonopy_atoms_to_ase(phonopy.supercell) + ), + } + structure_data_pymatgen_atoms = { + "supercell_structure_pymatgen_atoms": supercell_atoms, + } + + dumpfn(structure_data_phonopy_atoms, "structure_data_phonopy_atoms.json") + dumpfn(structure_data_pymatgen_atoms, "structure_data_pymatgen_atoms.json") + + primitive_atoms_phonopy = phonopy_atoms_to_ase(phonopy.primitive) + primitive_atoms_pymatgen = AseAtomsAdaptor.get_atoms(parent_structure) + + prim_structure_data_phonopy_atoms = { + "primitive_structure_phonopy_atoms": MSONAtoms(primitive_atoms_phonopy), + } + prim_structure_data_pymatgen_atoms = { + "primitive_structure_pymatgen_atoms": primitive_atoms_pymatgen, + } + # dumpfn primitive and supercell atoms + dumpfn( + prim_structure_data_phonopy_atoms, "prim_structure_data_phonopy_atoms.json" + ) + dumpfn( + prim_structure_data_pymatgen_atoms, + "prim_structure_data_pymatgen_atoms.json", + ) + + primitive_structure_phonopy = AseAtomsAdaptor.get_structure( + primitive_atoms_phonopy + ) + primitive_structure_pymatgen = AseAtomsAdaptor.get_structure( + primitive_atoms_pymatgen + ) + + # dumpfn primitive and supercell structures + dumpfn(primitive_structure_phonopy, "primitive_structure_phonopy.json") + dumpfn(primitive_structure_pymatgen, "primitive_structure_pymatgen.json") + + primitive_structure_phonopy_new = SpacegroupAnalyzer( + primitive_structure_phonopy + ).find_primitive() + primitive_structure_pymatgen_new = SpacegroupAnalyzer( + primitive_structure_pymatgen + ).find_primitive() + + # dumpfn primitive and supercell structures + dumpfn(primitive_structure_phonopy_new, "primitive_structure_phonopy_new.json") + dumpfn( + primitive_structure_pymatgen_new, "primitive_structure_pymatgen_new.json" + ) + + if fitting_data["best_n_imaginary"] == 0: + # if True: + logger.info("No imaginary modes! Writing ShengBTE files") + atoms = AseAtomsAdaptor.get_atoms(parent_structure) + # # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD_1", atoms, order=3) + # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD_from_fcs_pymatgen_struct", + # atoms) + fcs.write_to_phonopy("FORCE_CONSTANTS_2ND", format="text") + + # primitive_atoms_phonopy = phonopy_atoms_to_ase(phonopy.primitive) + # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD_from_fcs_phonopy_struct", + # primitive_atoms_phonopy) + + HiphiveForceConstants.write_to_phonopy( + fcs, "fc2.hdf5", "hdf5" + ) # ForceConstants + # ForceConstants.write_to_phono3py(fcs, "fc3.hdf5", "hdf5") + HiphiveForceConstants.write_to_phono3py( + fcs, "fc3.hdf5", order=3 + ) # ForceConstants + + ### detour from hdf5 + supercell_atoms_phonopy = phonopy_atoms_to_ase(phonopy.supercell) + supercell_atoms_pymatgen = AseAtomsAdaptor.get_atoms(supercell_structure) + + # dumpfn(supercell_atoms_phonopy, "supercell_atoms_phonopy.json") + dumpfn(supercell_atoms_pymatgen, "supercell_atoms_pymatgen.json") + + # check if the supercell_atoms are the same + if supercell_atoms_phonopy == supercell_atoms_pymatgen: + logger.info("supercell_atoms are the same") + else: + logger.info("supercell_atoms are different") + + supercell_atoms = supercell_atoms_phonopy + # supercell_atoms = supercell_atoms_pymatgen + # fcs = ForceConstants.read_phono3py(supercell_atoms, "fc3.hdf5", order=3) + fcs = HiphiveForceConstants.read_phono3py( + supercell_atoms, "fc3.hdf5", order=3 + ) # ForceConstants + # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, order=3, fc_tol=1e-4) + # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD", atoms, fc_tol=1e-4) + + # # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD_from_hdf5_phonopy_struct", + # atoms) # this was the original way of writing shengBTE files + # fcs.write_to_shengBTE("FORCE_CONSTANTS_3RD_from_hdf5_phonopy_struct", + # phonopy_atoms_to_ase(phonopy.primitive)) + + supercell_atoms = supercell_atoms_pymatgen + fcs = HiphiveForceConstants.read_phono3py( + supercell_atoms, "fc3.hdf5", order=3 + ) # ForceConstants + fcs.write_to_shengBTE( + "FORCE_CONSTANTS_3RD_from_hdf5_pymatgen_struct", atoms, order=3 + ) + + else: + logger.info(f"best_n_imaginary = {fitting_data['best_n_imaginary']}") + logger.info("ShengBTE files not written due to imaginary modes.") + logger.info("You may want to perform phonon renormalization.") + + current_dir = os.getcwd() + + outputs: dict[str, Any] = { + "thermal_data": thermal_data, + "anharmonic_data": anharmonic_data, + "fitting_data": fitting_data, + "param": param, + "current_dir": current_dir, + } + + return outputs + + @staticmethod + def fit_force_constants( + parent_structure: Structure, + supercell_matrix: np.ndarray, + supercell_structure: Structure, + structures: list[Atoms], + all_cutoffs: list[list[float]], + disp_cut: float | None = 0.055, + imaginary_tol: float = 0.025, # in THz + fit_method: str | None = "rfe", + n_jobs: int | None = -1, + fit_kwargs: dict | None = None, + ) -> tuple: + """ + Fit FC using hiphive. + + Fit force constants using hiphive and select the optimum cutoff values. + The optimum cutoffs will be determined according to: + 1. Number imaginary modes < ``max_n_imaginary``. + 2. Most negative imaginary frequency < ``max_imaginary_freq``. + 3. Least number of imaginary modes. + 4. Lowest free energy at 300 K. + If criteria 1 and 2 are not satisfied, None will be returned as the + force constants. + Args: + parent_structure: Initial input structure. + supercell_matrix: Supercell transformation matrix. + structures: A list of ase atoms objects with "forces" and + "displacements" arrays added, as required by hiPhive. + all_cutoffs: A nested list of cutoff values to trial. Each set of + cutoffs contains the radii for different orders starting with second + order. + disp_cut: determines the mean displacement of perturbed + structure to be included in harmonic (<) or anharmonic (>) fitting + imaginary_tol: Tolerance used to decide if a phonon mode is imaginary, + in THz. + max_n_imaginary: Maximum number of imaginary modes allowed in the + the final fitted force constant solution. If this criteria is not + reached by any cutoff combination this FireTask will fizzle. + max_imaginary_freq: Maximum allowed imaginary frequency in the + final fitted force constant solution. If this criteria is not + reached by any cutoff combination this FireTask will fizzle. + fit_method: Method used for fitting force constants. This can be + any of the values allowed by the hiphive ``Optimizer`` class. + n_jobs: Number of processors to use for fitting coefficients. + -1 means use all processors. + fit_kwargs: Additional arguments passed to the hiphive force constant + optimizer. + + Returns + ------- + A tuple of the best fitted force constants as a hiphive + ``SortedForceConstants`` object, array of parameters, cluster space, + and a dictionary of information on the fitting results. + """ + logger.info("Starting force constant fitting.") + logger.info(f"disp_cut is {disp_cut}") + logger.info(f"fit_method is {fit_method}") + + fitting_data: dict[str, Any] = { + "cutoffs": [], + "rmse_test": [], + "imaginary": [], + "cs_dofs": [], + "n_imaginary": [], + "parent_structure": parent_structure, + "supercell_structure": supercell_structure, + "supercell_matrix": supercell_matrix, + "fit_method": fit_method, + "disp_cut": disp_cut, + "imaginary_tol": imaginary_tol, + "best_cutoff": None, + "best_rmse": np.inf, + } + + best_fit = { + "n_imaginary": np.inf, + "rmse_test": np.inf, + "imaginary": None, + "cs_dofs": [None, None, None], + "cluster_space": None, + "force_constants": None, + "parameters": None, + "cutoffs": None, + "force_constants_potential": None, + } + # all_cutoffs = all_cutoffs[0] #later change it back to all_cutoffs + n_cutoffs = len(all_cutoffs) + logger.info(f"len_cutoffs={n_cutoffs}") + + fit_kwargs = fit_kwargs if fit_kwargs else {} + if fit_method == "rfe" and n_jobs == -1: + fit_kwargs["n_jobs"] = 1 + elif fit_method == "lasso": + fit_kwargs["lasso"] = dict(max_iter=1000) + elif fit_method == "elasticnet": + fit_kwargs = {"max_iter": 100000} + + logger.info(f"CPU COUNT: {os.cpu_count()}") + + logger.info("We are starting Joblib_s parallelized jobs") + + cutoff_results = [] + for i, cutoffs in enumerate(all_cutoffs): + result = PhononBSDOSDoc._run_cutoffs( + i=i, + cutoffs=cutoffs, + n_cutoffs=n_cutoffs, + parent_structure=parent_structure, + supercell_structure=supercell_structure, + structures=structures, + supercell_matrix=supercell_matrix, + fit_method=fit_method, + disp_cut=disp_cut, + imaginary_tol=imaginary_tol, + fit_kwargs=fit_kwargs, + ) + cutoff_results.append(result) + + for result in cutoff_results: + if result is None: + logger.info("result is None") + continue + logger.info(f"result = {result}") + if result != {}: + if "cutoffs" in result: + fitting_data["cutoffs"].append(result["cutoffs"]) + else: + logger.info("Cutoffs not found in result") + continue + if "rmse_test" in result: + fitting_data["rmse_test"].append(result["rmse_test"]) + if "n_imaginary" in result: + fitting_data["n_imaginary"].append(result["n_imaginary"]) + if "cs_dofs" in result: + fitting_data["cs_dofs"].append(result["cs_dofs"]) + if "imaginary" in result: + fitting_data["imaginary"].append(result["imaginary"]) + if "rmse_test" in result and ( + result["rmse_test"] < best_fit["rmse_test"] + ): + best_fit.update(result) + fitting_data["best_cutoff"] = result["cutoffs"] + fitting_data["best_rmse"] = result["rmse_test"] + else: + logger.info("result is an empty dictionary") + + logger.info(f"best_fit = {best_fit}") + logger.info(f"fitting_data = {fitting_data}") + logger.info("Finished fitting force constants.") + + return ( + best_fit["force_constants"], + best_fit["parameters"], + best_fit["cluster_space"], + fitting_data, + best_fit["force_constants_potential"], + ) + + @staticmethod + def _run_cutoffs( + i: int, + cutoffs: list[float], + n_cutoffs: int, + parent_structure: Structure, + supercell_structure: Structure, + structures: list[Atoms], + supercell_matrix: np.ndarray, + fit_method: str | None, + disp_cut: float | None, + fit_kwargs: dict[str, Any], + imaginary_tol: float = 0.025, # in THz + ) -> dict[str, Any]: + logger.info(f"Testing cutoffs {i+1} out of {n_cutoffs}: {cutoffs}") + logger.info(f"fit_method is {fit_method}") + supercell_atoms = AseAtomsAdaptor.get_atoms(supercell_structure) + supercell_atoms = structures[0] + logger.info(f"supercell_atoms = {supercell_atoms}") + + if not is_cutoff_allowed(supercell_atoms, max(cutoffs)): + logger.info( + "Skipping cutoff due as it is not commensurate with supercell size" + ) + return {} + + try: + cs = ClusterSpace( + supercell_atoms, cutoffs, symprec=1e-3, acoustic_sum_rules=True + ) + except ValueError: + logger.info("ValueError encountered, moving to the next cutoff") + return {} + + logger.debug(cs.__repr__()) + logger.info(cs) + cs_2_dofs = cs.get_n_dofs_by_order(2) + cs_3_dofs = cs.get_n_dofs_by_order(3) + cs_4_dofs = cs.get_n_dofs_by_order(4) + cs_dofs = [cs_2_dofs, cs_3_dofs, cs_4_dofs] + logger.info(cs_dofs) + n2nd = cs.get_n_dofs_by_order(2) + nall = cs.n_dofs + + logger.info("Fitting harmonic force constants separately") + separate_fit = True # Change this back to true + logger.info(f"disp_cut = {disp_cut}") + + sc = PhononBSDOSDoc.get_structure_container( + cs, structures, separate_fit, disp_cut, ncut=n2nd, param2=None + ) + opt = Optimizer(sc.get_fit_data(), fit_method, [0, n2nd], **fit_kwargs) + try: + opt.train() + except Exception: + logger.exception( + f"Error occurred in opt.train() for cutoff: {i}, {cutoffs}" + ) + + return {} + param_harmonic = opt.parameters # harmonic force constant parameters + param_tmp = np.concatenate( + (param_harmonic, np.zeros(cs.n_dofs - len(param_harmonic))) + ) + + parent_phonopy = get_phonopy_structure(parent_structure) + phonopy = Phonopy(parent_phonopy, supercell_matrix=supercell_matrix) + phonopy.primitive.get_number_of_atoms() + + # Ensure supercell_matrix is a numpy array + supercell_matrix = np.array(supercell_matrix) + mesh = supercell_matrix.diagonal() * 2 + mesh = [10, 10, 10] # TODO change this later + + fcp = ForceConstantPotential(cs, param_tmp) + # supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) + logger.info(f"supercell atoms = {supercell_atoms}") + fcs = fcp.get_force_constants(supercell_atoms) + logger.info("Did you get the large Condition number error?") + + phonopy.set_force_constants(fcs.get_fc_array(2)) + phonopy.set_mesh(mesh, is_eigenvectors=True, is_mesh_symmetry=True) + phonopy.run_mesh(mesh, with_eigenvectors=True, is_mesh_symmetry=True) + omega = phonopy.mesh.frequencies # THz + omega = np.sort(omega.flatten()) + logger.info(f"omega_one_shot_fit = {omega}") + imaginary = np.any(omega < -1e-3) + logger.info(f"imaginary_one_shot_fit = {imaginary}") + n_imaginary = int(np.sum(omega < -np.abs(imaginary_tol))) + + if imaginary: + logger.info( + "Imaginary modes found! Fitting anharmonic force constants separately" + ) + sc = PhononBSDOSDoc.get_structure_container( + cs, structures, separate_fit, disp_cut, ncut=n2nd, param2=param_harmonic + ) + opt = Optimizer(sc.get_fit_data(), fit_method, [n2nd, nall], **fit_kwargs) + + try: + opt.train() + except Exception: + logger.exception( + f"Error occurred in opt.train() for cutoff: {i}, {cutoffs}" + ) + return {} + param_anharmonic = opt.parameters # anharmonic force constant parameters + + parameters = np.concatenate((param_harmonic, param_anharmonic)) # combine + + if nall != len(parameters): + raise ValueError("nall is not equal to the length of parameters.") + logger.info(f"Training complete for cutoff: {i}, {cutoffs}") + + parent_phonopy = get_phonopy_structure(parent_structure) + phonopy = Phonopy(parent_phonopy, supercell_matrix=supercell_matrix) + phonopy.primitive.get_number_of_atoms() + mesh = supercell_matrix.diagonal() * 2 + mesh = [10, 10, 10] # TODO change this later + + fcp = ForceConstantPotential(cs, parameters) + # supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) + logger.info(f"supercell atoms = {supercell_atoms}") + fcs = fcp.get_force_constants(supercell_atoms) + logger.info("Did you get the large Condition number error?") + + phonopy.set_force_constants(fcs.get_fc_array(2)) + phonopy.set_mesh( + mesh, is_eigenvectors=False, is_mesh_symmetry=False + ) # run_mesh(is_gamma_center=True) + phonopy.run_mesh( + mesh=100.0, with_eigenvectors=False, is_mesh_symmetry=False + ) + omega = phonopy.mesh.frequencies # THz + omega = np.sort(omega.flatten()) + logger.info(f"omega_seperate_fit = {omega}") + imaginary = np.any(omega < -1e-3) + logger.info(f"imaginary_seperate_fit = {imaginary}") + + # Phonopy's way of calculating phonon frequencies + structure_phonopy = get_phonopy_structure(parent_structure) + phonon = Phonopy(structure_phonopy, supercell_matrix=supercell_matrix) + phonon.set_force_constants(fcs.get_fc_array(2)) + phonon.run_mesh(mesh, is_mesh_symmetry=False, is_gamma_center=True) + mesh = phonon.get_mesh_dict() + omega = mesh["frequencies"] + omega = np.sort(omega.flatten()) + logger.info(f"omega_phonopy_seperate_fit = {omega}") + # imaginary = np.any(omega < -1e-3) + # logger.info(f"imaginary_phonopy_seperate_fit = {imaginary}") + imaginary = np.any(omega < -imaginary_tol) + logger.info(f"imaginary_phonopy_seperate_fit = {imaginary}") + n_imaginary = int(np.sum(omega < -np.abs(imaginary_tol))) + + else: + logger.info("No imaginary modes! Fitting all force constants in one shot") + separate_fit = False + sc = PhononBSDOSDoc.get_structure_container( + cs, structures, separate_fit, disp_cut=None, ncut=None, param2=None + ) + opt = Optimizer(sc.get_fit_data(), fit_method, [0, nall], **fit_kwargs) + + try: + opt.train() + except Exception: + logger.exception( + f"Error occurred in opt.train() for cutoff: {i}, {cutoffs}" + ) + return {} + parameters = opt.parameters + logger.info(f"Training complete for cutoff: {i}, {cutoffs}") + + logger.info(f"parameters before enforcing sum rules {parameters}") + logger.info(f"Memory use: {psutil.virtual_memory().percent} %") + parameters = enforce_rotational_sum_rules( + cs, parameters, ["Huang", "Born-Huang"] + ) + fcp = ForceConstantPotential(cs, parameters) + # supercell_atoms = phonopy_atoms_to_ase(phonopy.supercell) + fcs = fcp.get_force_constants(supercell_atoms) + logger.info(f"FCS generated for cutoff {i}, {cutoffs}") + + return { + "cutoffs": cutoffs, + "rmse_test": opt.rmse_test, + "cluster_space": sc.cluster_space, + "parameters": parameters, + "force_constants": fcs, + "force_constants_potential": fcp, + "imaginary": imaginary, + "cs_dofs": cs_dofs, + "n_imaginary": n_imaginary, + } + + @staticmethod + def get_structure_container( + cs: ClusterSpace, + structures: list[Atoms], + separate_fit: bool, + disp_cut: float, + ncut: int, + param2: np.ndarray, + ) -> StructureContainer: + """Get a hiPhive StructureContainer from cutoffs and a list of atoms objects. + + Args: + cs: ClusterSpace + structures: A list of ase atoms objects with the "forces" and + "displacements" arrays included. + separate_fit: Boolean to determine whether harmonic and anharmonic fitting + are to be done separately (True) or in one shot (False) + disp_cut: if separate_fit true, determines the mean displacement of + perturbed structure to be included in harmonic (<) or anharmonic (>) + fitting + ncut: the parameter index where fitting separation occurs + param2: previously fit parameter array (harmonic only for now, hence 2). + + Returns + ------- + A hiPhive StructureContainer. + """ + sc = StructureContainer(cs) + logger.info(f"sc = {sc}") + logger.info(f"initial shape of fit matrix = {sc.data_shape}") + saved_structures = [] + for _, structure in enumerate(structures): + displacements = structure.get_array("displacements") + forces = structure.get_array("forces") + # Calculate mean displacements + mean_displacements = np.linalg.norm(displacements, axis=1).mean() + logger.info(f"Mean displacements: {mean_displacements}") + # Calculate mean forces + mean_forces = np.linalg.norm(forces, axis=1).mean() + logger.info(f"Mean Forces: {mean_forces}") + # Calculate standard deviation of displacements + std_displacements = np.linalg.norm(displacements, axis=1).std() + logger.info(f"Standard deviation of displacements: " f"{std_displacements}") + # Calculate standard deviation of forces + std_forces = np.linalg.norm(forces, axis=1).std() + logger.info(f"Standard deviation of forces: " f"{std_forces}") + if not separate_fit: # fit all + sc.add_structure(structure) + # for harmonic fitting + elif separate_fit and param2 is None and mean_displacements < disp_cut: + logger.info("We are in harmonic fitting if statement") + # logger.info(f"mean_disp = {mean_displacements}") + logger.info(f"mean_forces = {mean_forces}") + sc.add_structure(structure) + # for anharmonic fitting + elif separate_fit and param2 is not None and mean_displacements >= disp_cut: + logger.info("We are in anharmonic fitting if statement") + # logger.info(f"mean_disp = {mean_displacements}") + logger.info(f"mean_forces = {mean_forces}") + sc.add_structure(structure) + saved_structures.append(structure) + + logger.info( + "final shape of fit matrix (total # of atoms in all added" + f"supercells, n_dofs) = (rows, columns) = {sc.data_shape}" + ) + logger.info("We have completed adding structures") + logger.info(f"sc.get_fit_data() = {sc.get_fit_data()}") + + if separate_fit and param2 is not None: # do after anharmonic fitting + a_mat = sc.get_fit_data()[0] # displacement matrix + f_vec = sc.get_fit_data()[1] # force vector + f_vec -= np.dot(a_mat[:, :ncut], param2) # subtract harmonic forces + sc.delete_all_structures() + for i, structure in enumerate(saved_structures): + natoms = structure.get_global_number_of_atoms() + ndisp = natoms * 3 + structure.set_array( + "forces", f_vec[i * ndisp : (i + 1) * ndisp].reshape(natoms, 3) + ) + sc.add_structure(structure) + + logger.debug(sc.__repr__()) + + return sc + + @staticmethod + def harmonic_properties( + structure: Structure, + supercell_matrix: np.ndarray, + fcs: ForceConstants, + temperature: list, + imaginary_tol: float = 0.025, # in THz + mesh: list = None, + ) -> tuple[dict, Phonopy]: + """ + Thermodynamic (harmonic) phonon props calculated using the force constants. + + Args: + structure: The parent structure. + supercell_matrix: The supercell transformation matrix. + force_constants: The force constants in numpy format. + imaginary_tol: Tolerance used to decide if a phonon mode is imaginary, + in THz. + + Returns + ------- + A tuple of the number of imaginary modes at Gamma, the minimum phonon + frequency at Gamma, and the free energy, entropy, and heat capacity + """ + logger.info("Evaluating harmonic properties...") + fcs2 = fcs.get_fc_array(2) + logger.info("fcs2 & fcs3 read...") + logger.info(f"fcs2 = {fcs2}") + parent_phonopy = get_phonopy_structure(structure) + phonopy = Phonopy(parent_phonopy, supercell_matrix=supercell_matrix) + # natom = phonopy.primitive.get_number_of_atoms() + # Ensure supercell_matrix is a numpy array + supercell_matrix = np.array(supercell_matrix) + mesh = supercell_matrix.diagonal() * 2 + mesh = [10, 10, 10] # TODO change this later + logger.info(f"Mesh: {mesh}") + + phonopy.set_force_constants(fcs2) + phonopy.set_mesh(mesh, is_eigenvectors=True, is_mesh_symmetry=False) + phonopy.run_thermal_properties(temperatures=temperature) + logger.info("Thermal properties successfully run!") + + _, free_energy, entropy, heat_capacity = phonopy.get_thermal_properties() + + ## Use the following lines to convert the units to eV/atom + # free_energy *= 1000/sp.constants.Avogadro/eV2J/natom # kJ/mol to eV/atom + # entropy *= 1/sp.constants.Avogadro/eV2J/natom # J/K/mol to eV/K/atom + # heat_capacity *= 1/sp.constants.Avogadro/ev2j/natom # J/K/mol to eV/K/atom + + freq = phonopy.mesh.frequencies # in THz + logger.info(f"Frequencies: {freq}") + logger.info(f"freq_flatten = {np.sort(freq.flatten())}") + + n_imaginary = int(np.sum(freq < -np.abs(imaginary_tol))) + + if n_imaginary == 0: + logger.info("No imaginary modes!") + else: + logger.warning("Imaginary modes found!") + + ### Added from Phonon workflow + logger.info("starting the thermal prop calculation as per the phonon workflow") + from pymatgen.io.phonopy import ( + get_ph_bs_symm_line, + get_ph_dos, + get_pmg_structure, + ) + + from atomate2.common.schemas.phonons import PhononBSDOSDoc + + kpath_scheme = "seekpath" + symprec = 1e-4 + from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections + from pymatgen.io.vasp import Kpoints + from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter + + # with phonopy.load("phonopy.yaml") the phonopy API can be used + phonopy.save("phonopy.yaml") + + # get phonon band structure + kpath_dict, kpath_concrete = PhononBSDOSDoc.get_kpath( + structure=get_pmg_structure(phonopy.primitive), + kpath_scheme=kpath_scheme, + symprec=symprec, + ) + + # npoints_band = 101 + qpoints, connections = get_band_qpoints_and_path_connections( + kpath_concrete, + npoints=101, # changed from 101 + ) + + # phonon band structures will always be computed + filename_band_yaml = "phonon_band_structure.yaml" + + # TODO: potentially add kwargs to avoid computation of eigenvectors + phonopy.run_band_structure( + qpoints, + path_connections=connections, + with_eigenvectors=False, + is_band_connection=False, + ) + phonopy.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=False + ) + new_plotter = PhononBSPlotter(bs=bs_symm_line) + new_plotter.save_plot( + filename="phonon_band_structure.pdf", + units="THz", + ) + + # get phonon density of states + filename_dos_yaml = "phonon_dos.yaml" + + kpoint_density_dos = 7000 + kpoint = Kpoints.automatic_density( + structure=get_pmg_structure(phonopy.primitive), + kppa=kpoint_density_dos, + force_gamma=True, + ) + phonopy.run_mesh(kpoint.kpts[0]) + phonopy.run_total_dos() + phonopy.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="phonon_dos.pdf", + units="THz", + ) + + # compute vibrational part of free energies per formula unit + temperature_range = np.arange(0, 1000, 10) + + free_energies = [ + dos.helmholtz_free_energy( + temp=temp, structure=get_pmg_structure(phonopy.primitive) + ) + for temp in temperature_range + ] + + entropies = [ + dos.entropy(temp=temp, structure=get_pmg_structure(phonopy.primitive)) + for temp in temperature_range + ] + + internal_energies = [ + dos.internal_energy( + temp=temp, structure=get_pmg_structure(phonopy.primitive) + ) + for temp in temperature_range + ] + + heat_capacities = [ + dos.cv(temp=temp, structure=get_pmg_structure(phonopy.primitive)) + for temp in temperature_range + ] + + # save a json file of temperature_range, free_energies, entropies, + # internal_energies, heat_capacities + # Organize data into a dictionary + data = { + "temperature_range": temperature_range.tolist(), + "free_energies": free_energies, + "entropies": entropies, + "internal_energies": internal_energies, + "heat_capacities": heat_capacities, + } + + # Define the output file path + output_file_path = "thermodynamic_properties.json" + + # Write the dictionary to a JSON file + with open(output_file_path, "w") as f: + json.dump(data, f, indent=4) # `indent=4` for pretty printing + + logger.info("ending the thermal prop calculation as per the phonon workflow") + + ### Added from Phonon workflow + + return { + "temperature": temperature, + "free_energy": free_energy, + "entropy": entropy, + "heat_capacity": heat_capacity, + "n_imaginary": n_imaginary, + }, phonopy + + @staticmethod + def anharmonic_properties( + phonopy: Phonopy, + fcs: ForceConstants, + temperature: list, + heat_capacity: np.ndarray, + n_imaginary: float, + bulk_modulus: float | None = None, + ) -> dict: + """ + Anharmonic properties calculated using the force constants. + + Args: + phonopy: The phonopy object. + fcs: The force constants in numpy format. + temperature: The temperature range. + heat_capacity: The heat capacity in J/K/mol. + n_imaginary: The number of imaginary modes. + bulk_modulus: The bulk modulus. + + Returns + ------- + A dictionary of the Gruneisen parameter, thermal expansion, and + expansion fraction. + """ + if n_imaginary == 0: + logger.info("Evaluating anharmonic properties...") + fcs2 = fcs.get_fc_array(2) + fcs3 = fcs.get_fc_array(3) + grun, cte, dlfrac = PhononBSDOSDoc.gruneisen( + phonopy, + fcs2, + fcs3, + temperature, + heat_capacity, + bulk_modulus=bulk_modulus, + ) + else: # do not calculate these if imaginary modes exist + logger.warning( + "Gruneisen and thermal expansion cannot be calculated with" + "imaginary modes. All set to 0." + ) + grun = np.zeros((len(temperature), 3)) + cte = np.zeros((len(temperature), 3)) + dlfrac = np.zeros((len(temperature), 3)) + + return { + "gruneisen": grun, + "thermal_expansion": cte, + "expansion_fraction": dlfrac, + } + + @staticmethod + def get_total_grun( + omega: np.ndarray, grun: np.ndarray, kweight: np.ndarray, t: float + ) -> np.ndarray: + """ + Calculate the total Gruneisen parameter. + + Args: + omega: The phonon frequencies. + grun: The Gruneisen parameters. + kweight: The k-point weights. + t: The temperature. + + Returns + ------- + The total Gruneisen parameter. + """ + # total = 0 + weight = 0 + nptk = omega.shape[0] + nbands = omega.shape[1] + omega = abs(omega) * 1e12 * 2 * np.pi + if t == 0: + total = np.zeros((3, 3)) + grun_total_diag = np.zeros(3) + else: + total = np.zeros((3, 3)) + for i in range(nptk): + for j in range(nbands): + x = hbar * omega[i, j] / (2.0 * kb * t) + dbe = (x / np.sinh(x)) ** 2 + weight += dbe * kweight[i] + total += dbe * kweight[i] * grun[i, j] + total = total / weight + grun_total_diag = np.array([total[0, 2], total[1, 1], total[2, 0]]) + + def percent_diff(a: float, b: float) -> float: + return abs((a - b) / b) + + # This process preserves cell symmetry upon thermal expansion, i.e., + # it prevents symmetry-identical directions from inadvertently expanding + # by different ratios when/if the Gruneisen routine returns slightly + # different ratios for those directions + avg012 = np.mean( + (grun_total_diag[0], grun_total_diag[1], grun_total_diag[2]) + ) + avg01 = np.mean((grun_total_diag[0], grun_total_diag[1])) + avg02 = np.mean((grun_total_diag[0], grun_total_diag[2])) + avg12 = np.mean((grun_total_diag[1], grun_total_diag[2])) + if percent_diff(grun_total_diag[0], avg012) < 0.1: + if percent_diff(grun_total_diag[1], avg012) < 0.1: + if percent_diff(grun_total_diag[2], avg012) < 0.1: # all similar + grun_total_diag[0] = avg012 + grun_total_diag[1] = avg012 + grun_total_diag[2] = avg012 + elif ( + percent_diff(grun_total_diag[2], avg02) < 0.1 + ): # 0 and 2 similar + grun_total_diag[0] = avg02 + grun_total_diag[2] = avg02 + elif ( + percent_diff(grun_total_diag[2], avg12) < 0.1 + ): # 1 and 2 similar + grun_total_diag[1] = avg12 + grun_total_diag[2] = avg12 + else: + pass + elif percent_diff(grun_total_diag[1], avg01) < 0.1: # 0 and 1 similar + grun_total_diag[0] = avg01 + grun_total_diag[1] = avg01 + elif percent_diff(grun_total_diag[1], avg12) < 0.1: # 1 and 2 similar + grun_total_diag[1] = avg12 + grun_total_diag[2] = avg12 + else: + pass + elif percent_diff(grun_total_diag[0], avg01) < 0.1: # 0 and 1 similar + grun_total_diag[0] = avg01 + grun_total_diag[1] = avg01 + elif percent_diff(grun_total_diag[0], avg02) < 0.1: # 0 and 2 similar + grun_total_diag[0] = avg02 + grun_total_diag[2] = avg02 + else: # nothing similar + pass + + return grun_total_diag + + @staticmethod + def gruneisen( + phonopy: Phonopy, + fcs2: np.ndarray, + fcs3: np.ndarray, + temperature: list, + heat_capacity: np.ndarray, # in J/K/mol + bulk_modulus: float | None = None, # in GPa + ) -> tuple[list, list, list]: + """ + Calculate the Gruneisen parameter. + + Args: + phonopy: The phonopy object. + fcs2: The second order force constants. + fcs3: The third order force constants. + temperature: The temperature range. + heat_capacity: The heat capacity in J/K/mol. + bulk_modulus: The bulk modulus in GPa. + + Returns + ------- + A tuple of the Gruneisen parameter, linear thermal expansion coefficient, + and linear expansion fraction. + """ + gruneisen = Gruneisen(fcs2, fcs3, phonopy.supercell, phonopy.primitive) + gruneisen.set_sampling_mesh(phonopy.mesh_numbers, is_gamma_center=True) + gruneisen.run() + grun = gruneisen.get_gruneisen_parameters() # (nptk,nmode,3,3) + omega = gruneisen._frequencies # noqa: SLF001 + # qp = gruneisen._qpoints + kweight = gruneisen._weights # noqa: SLF001 + + grun_tot = [ + PhononBSDOSDoc.get_total_grun(omega, grun, kweight, temp) + for temp in temperature + ] + grun_tot = np.nan_to_num(np.array(grun_tot)) + + # linear thermal expansion coefficient and fraction + if bulk_modulus is None: + cte = None + dlfrac = None + else: + # heat_capacity *= eV2J*phonopy.primitive.get_number_of_atoms() + # # eV/K/atom to J/K + heat_capacity *= 1 / sp.constants.Avogadro # J/K/mol to J/K + # to convert from J/K/atom multiply by + # phonopy.primitive.get_number_of_atoms() + # Convert heat_capacity to an array if it's a scalar + # heat_capacity = np.array([heat_capacity]) + logger.info(f"heat capacity = {heat_capacity}") + vol = phonopy.primitive.get_volume() + + logger.info(f"grun_tot: {grun_tot}") + # logger.info(f"grun_tot shape: {grun_tot.shape}") + logger.info(f"heat_capacity shape: {heat_capacity.shape}") + logger.info(f"heat_capacity: {heat_capacity}") + logger.info(f"vol: {vol}") + logger.info(f"bulk_modulus: {bulk_modulus}") + # cte = grun_tot*heat_capacity.repeat(3)/(vol/10**30)/(bulk_modulus*10**9)/3 + cte = ( + grun_tot + * heat_capacity.repeat(3).reshape(len(heat_capacity), 3) + / (vol / 10**30) + / (bulk_modulus * 10**9) + / 3 + ) + cte = np.nan_to_num(cte) + if len(temperature) == 1: + dlfrac = cte * temperature + else: + dlfrac = PhononBSDOSDoc.thermal_expansion(temperature, cte) + logger.info(f"Gruneisen: \n {grun_tot}") + logger.info(f"Coefficient of Thermal Expansion: \n {cte}") + logger.info(f"Linear Expansion Fraction: \n {dlfrac}") + + return grun_tot, cte, dlfrac + + @staticmethod + def thermal_expansion( + temperature: list, + cte: np.array, + ) -> np.ndarray: + """ + Calculate the linear thermal expansion fraction. + + Args: + temperature: The temperature range. + cte: The linear thermal expansion coefficient. + + Returns + ------- + The linear thermal expansion fraction. + """ + if len(temperature) != len(cte): + raise ValueError("Length of temperature and cte lists must be equal.") + if 0 not in temperature: + temperature = [0, *temperature] + cte = np.array([np.array([0, 0, 0]), *list(cte)]) + temperature = np.array(temperature) + ind = np.argsort(temperature) + temperature = temperature[ind] + cte = np.array(cte)[ind] + # linear expansion fraction + dlfrac = copy.copy(cte) + for t in range(len(temperature)): + dlfrac[t, :] = np.trapezoid(cte[: t + 1, :], temperature[: t + 1], axis=0) + dlfrac = np.nan_to_num(dlfrac) + logger.info(f"Linear Expansion Fraction: \n {dlfrac}") + return dlfrac diff --git a/src/atomate2/forcefields/flows/hiphive.py b/src/atomate2/forcefields/flows/hiphive.py new file mode 100644 index 0000000000..67e9a213c5 --- /dev/null +++ b/src/atomate2/forcefields/flows/hiphive.py @@ -0,0 +1,96 @@ +"""Define the ForceField HiphiveMaker. + +Uses hiPhive, phono3py, phonopy & alma/shengbte for calculating harmonic & anharmonic +props of phonon. +""" + +# Basic Python packages +from __future__ import annotations + +import logging +from dataclasses import dataclass, field + +from atomate2.common.flows.hiphive import BaseHiphiveMaker +from atomate2.forcefields.jobs import ( + ForceFieldRelaxMaker, + ForceFieldStaticMaker, + MACERelaxMaker, + MACEStaticMaker, +) + +logger = logging.getLogger(__name__) + +__all__ = ["HiphiveMaker"] + +__author__ = "Alex Ganose, Junsoo Park, Zhuoying Zhu, Hrushikesh Sahasrabuddhe" +__email__ = "aganose@lbl.gov, jsyony37@lbl.gov, zyzhu@lbl.gov, hpsahasrabuddhe@lbl.gov" + + +@dataclass +class HiphiveMaker(BaseHiphiveMaker): + """ + Workflow to calc. interatomic force constants and vibrational props using hiPhive. + + A summary of the workflow is as follows: + 1. Structure relaxtion + 2. Calculate a supercell transformation matrix that brings the + structure as close as cubic as possible, with all lattice lengths + greater than 5 nearest neighbor distances. + 3. Perturb the atomic sites for each supercell using a Monte Carlo + rattle procedure. The atoms are perturbed roughly according to a + normal deviation. A number of standard deviation perturbation distances + are included. Multiple supercells may be generated for each perturbation + distance. + 4. Run static VASP calculations on each perturbed supercell to calculate + atomic forces. + 5. Aggregate the forces and conduct the fit atomic force constants using + the minimization schemes in hiPhive. + 6. Output the interatomic force constants, and phonon band structure and + density of states to the database. + 7. Optional: Perform phonon renormalization at finite temperature - useful + when unstable modes exist + 8. Optional: Solve the lattice thermal conductivity using ShengBTE and + output to the database + + Args + ---------- + name : str + Name of the flows produced by this maker. + static_energy_maker (BaseVaspMaker): + The VASP input generator for static calculations, default is StaticMaker. + bulk_relax_maker (BaseVaspMaker | None): + The VASP input generator for bulk relaxation, + default is DoubleRelaxMaker using TightRelaxMaker. + phonon_displacement_maker (BaseVaspMaker | None): + The VASP input generator for phonon displacement calculations, + default is PhononDisplacementMaker. + """ + + name: str = "Lattice-Dynamics-FORCE_FIELD" + static_energy_maker: ForceFieldStaticMaker | None = field( + default_factory=MACEStaticMaker + ) + bulk_relax_maker: ForceFieldRelaxMaker = field( + default_factory=lambda: MACERelaxMaker(relax_kwargs={"fmax": 0.00001}) + ) + phonon_displacement_maker: ForceFieldStaticMaker = field( + default_factory=MACEStaticMaker + ) + # ff_displacement_maker: ForceFieldStaticMaker = field( + # default_factory=CHGNetStaticMaker + # ) + # phonon_displacement_maker: ForceFieldStaticMaker = field( + # default_factory=lambda: CHGNetRelaxMaker(relax_kwargs={"fmax": 2}) + # ) + + @property + def prev_calc_dir_argname(self) -> None: + """Name of argument informing static maker of previous calculation directory. + + As this differs between different DFT codes (e.g., VASP, CP2K), it + has been left as a property to be implemented by the inheriting class. + + Note: this is only applicable if a relax_maker is specified; i.e., two + calculations are performed for each ordering (relax -> static) + """ + return diff --git a/src/atomate2/settings.py b/src/atomate2/settings.py index 912013fc27..bc696fb9c1 100644 --- a/src/atomate2/settings.py +++ b/src/atomate2/settings.py @@ -243,6 +243,9 @@ class Atomate2Settings(BaseSettings): "parsing QChem directories useful for storing duplicate of FW.json", ) + # Phono3py settings + PHONO3PY_CMD: str = Field("phono3py", description="Command to run phono3py.") + @model_validator(mode="before") @classmethod def load_default_settings(cls, values: dict[str, Any]) -> dict[str, Any]: diff --git a/src/atomate2/vasp/files.py b/src/atomate2/vasp/files.py index d24cefd5c7..e5ac79cffd 100644 --- a/src/atomate2/vasp/files.py +++ b/src/atomate2/vasp/files.py @@ -19,6 +19,8 @@ from atomate2.vasp.sets.base import VaspInputGenerator +__all__ = ["copy_vasp_outputs", "get_largest_relax_extension", "copy_hiphive_outputs"] + logger = logging.getLogger(__name__) @@ -37,7 +39,7 @@ def copy_vasp_outputs( For folders containing multiple calculations (e.g., suffixed with relax1, relax2, etc), this function will only copy the files with the highest numbered suffix and - the suffix will be removed. Additional vasp files will be also be copied with the + the suffix will be removed. Additional vasp files will be also be copied with the same suffix applied. Lastly, this function will gunzip any gzipped files. Parameters @@ -104,10 +106,8 @@ def copy_vasp_outputs( if relax_ext: all_files = optional_files + required_files files_to_rename = { - file.name.replace(".gz", ""): file.name.replace(relax_ext, "").replace( - ".gz", "" - ) - for file in all_files + k.name.replace(".gz", ""): k.name.replace(relax_ext, "").replace(".gz", "") + for k in all_files } rename_files(files_to_rename, allow_missing=True, file_client=file_client) @@ -202,3 +202,66 @@ def write_vasp_input_set( logger.info("Writing VASP input set.") vis.write_input(directory, **kwargs) + + +@auto_fileclient +def copy_hiphive_outputs( + src_dir: Path | str, + src_host: str | None = None, + file_client: FileClient | None = None, +) -> None: + """ + Copy Non-VASP output files to the current directory. + + Parameters + ---------- + src_dir : str or Path + The source directory. + src_host : str or None + The source hostname used to specify a remote filesystem. Can be given as + either "username@remote_host" or just "remote_host" in which case the username + will be inferred from the current user. If ``None``, the local filesystem will + be used as the source. + file_client : .FileClient + A file client to use for performing file operations. + """ + src_dir = strip_hostname(src_dir) # TODO: Handle hostnames properly. + + logger.info(f"Copying Non-VASP inputs from {src_dir}") + + directory_listing = file_client.listdir(src_dir, host=src_host) + optional_files = [] + + for file in [ + "cluster_space.cs", + "force_constants.fcs", + "force_constants_potential.fcp", + "parameters.txt", + "fitting_data.json", + "phonopy_params.yaml", + "thermal_data.json", + "renorm_thermal_data.json", + "structure.json", + "relaxed_structure.json", + "structure_data.json", + "perturbed_structures.json", + "perturbed_forces.json", + "perturbed_forces_new.json", + "fc2.hdf5", + "fc3.hdf5", + "FORCE_CONSTANTS_2ND", + "FORCE_CONSTANTS_3RD", + ]: + found_file = get_zfile(directory_listing, file, allow_missing=True) + if found_file is not None: + optional_files.append(found_file) + + copy_files( + src_dir, + src_host=src_host, + # include_files=required_files + optional_files, + include_files=optional_files, + file_client=file_client, + ) + + logger.info("Finished copying Non-VASP inputs") diff --git a/src/atomate2/vasp/flows/hiphive.py b/src/atomate2/vasp/flows/hiphive.py new file mode 100644 index 0000000000..745c6fe0d6 --- /dev/null +++ b/src/atomate2/vasp/flows/hiphive.py @@ -0,0 +1,146 @@ +"""Define the VASP HiphiveMaker. + +Uses hiPhive, phono3py, phonopy & alma/shengbte for calculating harmonic & anharmonic +props of phonon. +""" + +# Basic Python packages +from __future__ import annotations + +import logging +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +from atomate2.common.flows.hiphive import BaseHiphiveMaker +from atomate2.vasp.flows.core import DoubleRelaxMaker + +# Atomate2 packages +from atomate2.vasp.jobs.core import TightRelaxMaker +from atomate2.vasp.jobs.phonons import PhononDisplacementMaker +from atomate2.vasp.sets.core import StaticSetGenerator, TightRelaxSetGenerator + +if TYPE_CHECKING: + from atomate2.vasp.jobs.base import BaseVaspMaker + +logger = logging.getLogger(__name__) + +__all__ = ["HiphiveMaker"] + +__author__ = "Alex Ganose, Junsoo Park, Zhuoying Zhu, Hrushikesh Sahasrabuddhe" +__email__ = "aganose@lbl.gov, jsyony37@lbl.gov, zyzhu@lbl.gov, hpsahasrabuddhe@lbl.gov" + + +@dataclass +class HiphiveMaker(BaseHiphiveMaker): + """ + Workflow to calc. interatomic force constants and vibrational props using hiPhive. + + A summary of the workflow is as follows: + 1. Structure relaxtion + 2. Calculate a supercell transformation matrix that brings the + structure as close as cubic as possible, with all lattice lengths + greater than 5 nearest neighbor distances. + 3. Perturb the atomic sites for each supercell using a Monte Carlo + rattle procedure. The atoms are perturbed roughly according to a + normal deviation. A number of standard deviation perturbation distances + are included. Multiple supercells may be generated for each perturbation + distance + 4. Run static VASP calculations on each perturbed supercell to calculate + atomic forces. + 5. Aggregate the forces and conduct the fit atomic force constants using + the minimization schemes in hiPhive. + 6. Output the interatomic force constants, and phonon band structure and + density of states to the database. + 7. Optional: Perform phonon renormalization at finite temperature - useful + when unstable modes exist + 8. Optional: Solve the lattice thermal conductivity using ShengBTE and + output to the database. + + Args + ---------- + name : str + Name of the flows produced by this maker. + static_energy_maker (BaseVaspMaker): + The VASP input generator for static calculations, default is StaticMaker. + bulk_relax_maker (BaseVaspMaker | None): + The VASP input generator for bulk relaxation, + default is DoubleRelaxMaker using TightRelaxMaker. + phonon_displacement_maker (BaseVaspMaker | None): + The VASP input generator for phonon displacement calculations, + default is PhononDisplacementMaker. + """ + + name: str = "Lattice-Dynamics-VASP" + bulk_relax_maker: DoubleRelaxMaker = field( + default_factory=lambda: DoubleRelaxMaker.from_relax_maker( + TightRelaxMaker( + input_set_generator=TightRelaxSetGenerator( + user_incar_settings={ + "PREC": "Accurate", + "GGA": "PS", + "IBRION": 2, + "NSW": 20, + "NELMIN": 5, + "ISIF": 3, + # "ENCUT": 648.744200, + "EDIFF": 1.000000e-06, + # "EDIFFG": -1.000000e-08, + "ISMEAR": 0, + "SIGMA": 1.000000e-02, + "IALGO": 38, + "LREAL": ".FALSE.", + "ADDGRID": ".TRUE.", + "LWAVE": ".FALSE.", + "LCHARG": ".FALSE.", + "NPAR": 4, + } + ) + ) + ) + ) + phonon_displacement_maker: BaseVaspMaker | None = field( + default_factory=lambda: PhononDisplacementMaker( + input_set_generator=StaticSetGenerator( + user_kpoints_settings={"reciprocal_density": 100}, + user_incar_settings={ + "ADDGRID": True, + "ALGO": "Normal", + "EDIFF": 1e-06, + # "EDIFFG": -1.000000e-08, + "ENCUT": 600, + "GGA": "PS", + "IBRION": -1, + "ISIF": 3, + "ISMEAR": 0, + "ISPIN": 2, + "LAECHG": False, + "LASPH": True, + "LCHARG": False, + "LORBIT": 11, + "LREAL": "Auto", + "LVHAR": False, + "LVTOT": False, + "LWAVE": False, + # "MAGMOM": 250*0.6, + "NCORE": 6, + "NELM": 100, + "NSW": 0, + "PREC": "Accurate", + "SIGMA": 0.1, # changed from 0.1 + }, + # auto_ispin=True, + ) + ) + ) + + @property + def prev_calc_dir_argname(self) -> str: + """Name of argument informing static maker of previous calculation directory. + + As this differs between different DFT codes (e.g., VASP, CP2K), it + has been left as a property to be implemented by the inheriting class. + + Note: this is only applicable if a relax_maker is specified; i.e., two + calculations are performed for each ordering (relax -> static) + """ + return "prev_dir" diff --git a/src/atomate2/vasp/schemas/hiphive.py b/src/atomate2/vasp/schemas/hiphive.py new file mode 100644 index 0000000000..d64ba7d6ba --- /dev/null +++ b/src/atomate2/vasp/schemas/hiphive.py @@ -0,0 +1,7 @@ +"""Schemas for hiphive documents.""" + +from pydantic import BaseModel + + +class LTCDoc(BaseModel): + """Collection to store Lattice thermal conductivity.""" diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/INCAR.gz new file mode 100644 index 0000000000..022da53496 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..140caace7f Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/KPOINTS.gz new file mode 100644 index 0000000000..5553fe153d Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..f18c035bd9 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POSCAR.gz new file mode 100644 index 0000000000..2bb621f58d Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..3bc34861cb Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/CONTCAR.gz new file mode 100644 index 0000000000..ffa710f5e2 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/INCAR.gz new file mode 100644 index 0000000000..022da53496 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..140caace7f Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/KPOINTS.gz new file mode 100644 index 0000000000..5553fe153d Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..f18c035bd9 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/OUTCAR.gz new file mode 100644 index 0000000000..63cea0af08 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POSCAR.gz new file mode 100644 index 0000000000..2bb621f58d Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..3bc34861cb Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..5c10a0a90c Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_1_4/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/INCAR.gz new file mode 100644 index 0000000000..510730ceb3 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..a4052bb72d Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/KPOINTS.gz new file mode 100644 index 0000000000..61dc0530bf Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..990d991ad9 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POSCAR.gz new file mode 100644 index 0000000000..810397d92a Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..91af663b67 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/CONTCAR.gz new file mode 100644 index 0000000000..3e3dd01843 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/INCAR.gz new file mode 100644 index 0000000000..510730ceb3 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..a4052bb72d Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/KPOINTS.gz new file mode 100644 index 0000000000..61dc0530bf Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..990d991ad9 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/OUTCAR.gz new file mode 100644 index 0000000000..f0540bd9de Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POSCAR.gz new file mode 100644 index 0000000000..810397d92a Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..91af663b67 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..46a2e71e38 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_2_4/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/INCAR.gz new file mode 100644 index 0000000000..5526729680 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..3e18728e70 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/KPOINTS.gz new file mode 100644 index 0000000000..85bec82587 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..ed22c9311b Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POSCAR.gz new file mode 100644 index 0000000000..901623a3de Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..0f8f2ac662 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/CONTCAR.gz new file mode 100644 index 0000000000..9e3577fe9e Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/INCAR.gz new file mode 100644 index 0000000000..5526729680 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..3e18728e70 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/KPOINTS.gz new file mode 100644 index 0000000000..85bec82587 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..ed22c9311b Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/OUTCAR.gz new file mode 100644 index 0000000000..516277e02b Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POSCAR.gz new file mode 100644 index 0000000000..901623a3de Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..0f8f2ac662 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..c948be6bc0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_3_4/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/INCAR.gz new file mode 100644 index 0000000000..e161e7a9a4 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..5c7cd84c10 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/KPOINTS.gz new file mode 100644 index 0000000000..51eda54170 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..b2b64dc3f9 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POSCAR.gz new file mode 100644 index 0000000000..4ad5a9bde0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..ad72e020b5 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/CONTCAR.gz new file mode 100644 index 0000000000..50b17c5a9c Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/INCAR.gz new file mode 100644 index 0000000000..e161e7a9a4 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..5c7cd84c10 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/KPOINTS.gz new file mode 100644 index 0000000000..51eda54170 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..b2b64dc3f9 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/OUTCAR.gz new file mode 100644 index 0000000000..908676336a Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POSCAR.gz new file mode 100644 index 0000000000..4ad5a9bde0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..ad72e020b5 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..af00e73621 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/phonon_static_4_4/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/INCAR.gz new file mode 100644 index 0000000000..4c038e26d2 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..c6c97072a3 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/KPOINTS.gz new file mode 100644 index 0000000000..2a5d966fc3 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..17d4303ba0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POSCAR.gz new file mode 100644 index 0000000000..9eb6f9c88b Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..30bc977ee5 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/CONTCAR.gz new file mode 100644 index 0000000000..8b2a0781db Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/INCAR.gz new file mode 100644 index 0000000000..4c038e26d2 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..c6c97072a3 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/KPOINTS.gz new file mode 100644 index 0000000000..2a5d966fc3 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..17d4303ba0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/OUTCAR.gz new file mode 100644 index 0000000000..d4eb0e7936 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POSCAR.gz new file mode 100644 index 0000000000..9eb6f9c88b Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..30bc977ee5 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..cc164476b8 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_1/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/INCAR.gz new file mode 100644 index 0000000000..b3ad87ba62 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..c906e38ae0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/KPOINTS.gz new file mode 100644 index 0000000000..bc7149f63c Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..762bcf5dd6 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POSCAR.gz new file mode 100644 index 0000000000..f3344eeb85 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..f4db069962 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/CONTCAR.gz new file mode 100644 index 0000000000..6e18be0e72 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/INCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/INCAR.gz new file mode 100644 index 0000000000..b3ad87ba62 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..c906e38ae0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/KPOINTS.gz new file mode 100644 index 0000000000..bc7149f63c Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..762bcf5dd6 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/OUTCAR.gz new file mode 100644 index 0000000000..c574fef800 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POSCAR.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POSCAR.gz new file mode 100644 index 0000000000..f3344eeb85 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..f4db069962 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..990daa22ac Binary files /dev/null and b/tests/test_data/vasp/Si_hiphive/tight_relax_2/outputs/vasprun.xml.gz differ diff --git a/tests/vasp/flows/test_hiphive.py b/tests/vasp/flows/test_hiphive.py new file mode 100644 index 0000000000..583306b002 --- /dev/null +++ b/tests/vasp/flows/test_hiphive.py @@ -0,0 +1,667 @@ +from jobflow import run_locally +from numpy.testing import assert_allclose +from pymatgen.core.structure import Structure +from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine +from pymatgen.phonon.dos import PhononDos + +from atomate2.common.flows.hiphive import BaseHiphiveMaker +from atomate2.common.flows.phonons import BasePhononMaker +from atomate2.common.powerups import add_metadata_to_flow +from atomate2.common.schemas.hiphive import ( + ForceConstants, + PhononBSDOSDoc, + PhononComputationalSettings, + PhononJobDirs, + PhononUUIDs, + ThermalDisplacementData, +) +from atomate2.vasp.flows.core import DoubleRelaxMaker +from atomate2.vasp.flows.hiphive import HiphiveMaker +from atomate2.vasp.jobs.base import BaseVaspMaker +from atomate2.vasp.jobs.core import TightRelaxMaker +from atomate2.vasp.jobs.phonons import PhononDisplacementMaker + + +def test_hiphive_wf_vasp(mock_vasp, clean_dir, test_dir): + # mapping from job name to directory containing test files + ref_paths = { + "tight relax 1": "Si_hiphive/tight_relax_1", + "tight relax 2": "Si_hiphive/tight_relax_2", + "phonon static 1/4": "Si_hiphive/phonon_static_1_4", + "phonon static 2/4": "Si_hiphive/phonon_static_2_4", + "phonon static 3/4": "Si_hiphive/phonon_static_3_4", + "phonon static 4/4": "Si_hiphive/phonon_static_4_4", + } + + # settings passed to fake_run_vasp; adjust these to check for certain INCAR settings + fake_run_vasp_kwargs = { + "tight relax 1": {"incar_settings": ["NSW", "ISMEAR"]}, + "tight relax 2": {"incar_settings": ["NSW", "ISMEAR"]}, + "phonon static 1/4": {"incar_settings": ["NSW", "ISMEAR"]}, + "phonon static 2/4": {"incar_settings": ["NSW", "ISMEAR"]}, + "phonon static 3/4": {"incar_settings": ["NSW", "ISMEAR"]}, + "phonon static 4/4": {"incar_settings": ["NSW", "ISMEAR"]}, + } + + # automatically use fake VASP and write POTCAR.spec dulsring the test + mock_vasp(ref_paths, fake_run_vasp_kwargs) + + si_struct = Structure.from_file( + test_dir / "vasp/Si_hiphive/tight_relax_1/inputs/POSCAR.gz" + ) + + min_length = 6 + disp_cut = 0.05 + min_atoms = 30 + max_atoms = 50 + supercell_matrix_kwargs = { + "min_atoms": min_atoms, + "max_atoms": max_atoms, + "force_diagonal": True, + } + + job = HiphiveMaker( + supercell_matrix_kwargs=supercell_matrix_kwargs, + min_length=min_length, + code="vasp", + ).make( + mpid="mp-149", + structure=si_struct, + bulk_modulus=89, + disp_cut=disp_cut, + ) + + job = add_metadata_to_flow( + flow=job, + additional_fields={"project": "hiphive_unit_test"}, + class_filter=( + DoubleRelaxMaker, + TightRelaxMaker, + BaseVaspMaker, + PhononDisplacementMaker, + BasePhononMaker, + HiphiveMaker, + BaseHiphiveMaker, + ), + ) + + # run the flow or job and ensure that it finished running successfully + responses = run_locally(job, create_folders=True, ensure_success=True) + + # validate the outputs + assert isinstance(responses[job.jobs[-1].uuid][1].output, PhononBSDOSDoc) + assert_allclose( + responses[job.jobs[-1].uuid][1].output.free_energies, + [ + 5996.475041164254, + 5996.474199726997, + 5996.453005834313, + 5996.25654877896, + 5995.353731735742, + 5992.817833713253, + 5987.568222978423, + 5978.570835629409, + 5964.919147881639, + 5945.847051748049, + 5920.718784006014, + 5889.015101127187, + 5850.320531539404, + 5804.31181004804, + 5750.746924263767, + 5689.454540956843, + 5620.323896494623, + 5543.29535227478, + 5458.351796556117, + 5365.511000296318, + 5264.818956639393, + 5156.344172743387, + 5040.1728431841275, + 4916.404813348597, + 4785.150234051645, + 4646.526810312552, + 4500.657554082521, + 4347.668960096908, + 4187.6895342360895, + 4020.8486137953378, + 3847.2754283317154, + 3667.098358019808, + 3480.444353636197, + 3287.438488437673, + 3088.203617387223, + 2882.8601235242654, + 2671.5257348863015, + 2454.3153983769434, + 2231.341199439378, + 2002.712318421595, + 1768.5350161846118, + 1528.912642871017, + 1283.9456648708326, + 1033.7317059395778, + 778.3655991752091, + 517.9394471767135, + 252.5426882119234, + -17.737833364351307, + -292.81779386174793, + -572.6153153175517, + -857.0508897863285, + -1146.047303584171, + -1439.5295620632312, + -1737.4248153672058, + -2039.6622855164374, + -2346.17319508792, + -2656.890697687518, + -2971.7498103565536, + -3290.687348009972, + -3613.641859967181, + -3940.553568607455, + -4271.364310158468, + -4606.017477607981, + -4944.457965714119, + -5286.632118078246, + -5632.487676235766, + -5981.973730713599, + -6335.040673998304, + -6691.640155355467, + -7051.725037438825, + -7415.249354626328, + -7782.168273019967, + -8152.438052046258, + -8526.016007594952, + -8902.860476634527, + -9282.930783244366, + -9666.18720600498, + -10052.590946689326, + -10442.104100200186, + -10834.689625700265, + -11230.311318883712, + -11628.933785339732, + -12030.522414960777, + -12435.043357349836, + -12842.46349818321, + -13252.750436487013, + -13665.872462787453, + -14081.798538096675, + -14500.498273697767, + -14921.941911693975, + -15346.10030628897, + -15772.944905766402, + -16202.447735138441, + -16634.58137943451, + -17069.318967602594, + -17506.634156996963, + -17946.501118427204, + -18388.8945217448, + -18833.789521944404, + -19281.161745758254, + ], + rtol=1e-5, + atol=1e-5, + ) + + assert isinstance( + responses[job.jobs[-1].uuid][1].output.phonon_bandstructure, + PhononBandStructureSymmLine, + ) + assert isinstance(responses[job.jobs[-1].uuid][1].output.phonon_dos, PhononDos) + assert isinstance( + responses[job.jobs[-1].uuid][1].output.thermal_displacement_data, + ThermalDisplacementData, + ) + assert isinstance(responses[job.jobs[-1].uuid][1].output.structure, Structure) + assert_allclose( + responses[job.jobs[-1].uuid][1].output.temperatures, + [ + 0, + 10, + 20, + 30, + 40, + 50, + 60, + 70, + 80, + 90, + 100, + 110, + 120, + 130, + 140, + 150, + 160, + 170, + 180, + 190, + 200, + 210, + 220, + 230, + 240, + 250, + 260, + 270, + 280, + 290, + 300, + 310, + 320, + 330, + 340, + 350, + 360, + 370, + 380, + 390, + 400, + 410, + 420, + 430, + 440, + 450, + 460, + 470, + 480, + 490, + 500, + 510, + 520, + 530, + 540, + 550, + 560, + 570, + 580, + 590, + 600, + 610, + 620, + 630, + 640, + 650, + 660, + 670, + 680, + 690, + 700, + 710, + 720, + 730, + 740, + 750, + 760, + 770, + 780, + 790, + 800, + 810, + 820, + 830, + 840, + 850, + 860, + 870, + 880, + 890, + 900, + 910, + 920, + 930, + 940, + 950, + 960, + 970, + 980, + 990, + ], + ) + assert responses[job.jobs[-1].uuid][1].output.has_imaginary_modes is False + assert isinstance( + responses[job.jobs[-1].uuid][1].output.force_constants, ForceConstants + ) + assert isinstance(responses[job.jobs[-1].uuid][1].output.jobdirs, PhononJobDirs) + assert isinstance(responses[job.jobs[-1].uuid][1].output.uuids, PhononUUIDs) + assert_allclose( + responses[job.jobs[-1].uuid][1].output.supercell_matrix, + [[-2.0, -2.0, 2.0], [-2.0, 2.0, -2.0], [2.0, -2.0, -2.0]], + ) + assert_allclose( + responses[job.jobs[-1].uuid][1].output.primitive_matrix, + [[0.0, 0.0, -1.0], [0.0, -1.0, 0.0], [-1.0, 0.0, 0.0]], + rtol=1e-5, + atol=1e-10, + ) + assert responses[job.jobs[-1].uuid][1].output.code == "vasp" + assert isinstance( + responses[job.jobs[-1].uuid][1].output.phonopy_settings, + PhononComputationalSettings, + ) + assert responses[job.jobs[-1].uuid][1].output.phonopy_settings.npoints_band == 101 + assert ( + responses[job.jobs[-1].uuid][1].output.phonopy_settings.kpath_scheme + == "seekpath" + ) + assert ( + responses[job.jobs[-1].uuid][1].output.phonopy_settings.kpoint_density_dos + == 7000 + ) + assert_allclose( + responses[job.jobs[-1].uuid][1].output.entropies, + [ + 0.0, + 0.0003514422735110269, + 0.005925725827513983, + 0.0421930242943324, + 0.15432043988535035, + 0.3712097344657012, + 0.696064829429383, + 1.118532005016821, + 1.6245532858262144, + 2.2004381823394468, + 2.8338482900251196, + 3.5137936977700504, + 4.230489766226699, + 4.97527041996249, + 5.74054538827298, + 6.519761547651541, + 7.3073469656429015, + 8.098634327473317, + 8.889769541131486, + 9.677613662407056, + 10.459645261027463, + 11.23386823809233, + 11.998728031644424, + 12.753037524746905, + 13.495912858140203, + 14.226718667047802, + 14.945021895033998, + 15.650553183379555, + 16.343174812151844, + 17.02285422247485, + 17.6896422410225, + 18.34365523402357, + 18.98506052524345, + 19.614064512954503, + 20.230903011233345, + 20.835833419815362, + 21.429128394328497, + 22.01107074586041, + 22.581949346638353, + 23.142055858356123, + 23.6916821325635, + 24.231118159630707, + 24.760650465075617, + 25.280560870320738, + 25.791125549932108, + 26.292614329668734, + 26.78529017972682, + 27.269408865799768, + 27.74521872732361, + 28.21296055780888, + 28.67286756669662, + 29.125165405897548, + 29.57007224723024, + 30.00779889948354, + 30.438548955892493, + 30.862518964513292, + 31.27989861537905, + 31.690870939466592, + 32.095612515450156, + 32.49429368099625, + 32.8870787459943, + 33.274126205644905, + 33.655588951760954, + 34.03161448099267, + 34.40234509898039, + 34.76791811967774, + 35.128466059284285, + 35.484116824386234, + 35.83499389403375, + 36.18121649558882, + 36.52289977426247, + 36.86015495632818, + 37.193089506052495, + 37.52180727642644, + 37.84640865381455, + 38.166990696663746, + 38.483647268433174, + 38.79646916492007, + 39.105544236166025, + 39.410957503134135, + 39.712791269350575, + 40.01112522770538, + 40.30603656260583, + 40.59760004767405, + 40.88588813917651, + 41.170971065369265, + 41.45291691193741, + 41.731791703702086, + 42.007659482762385, + 42.2805823832338, + 42.5506207027384, + 42.81783297079617, + 43.0822760142604, + 43.344005019934286, + 43.60307359449979, + 43.85953382188391, + 44.11343631818204, + 44.36483028425238, + 44.61376355609009, + 44.86028265308503, + ], + atol=1e-6, + ) + assert_allclose( + responses[job.jobs[-1].uuid][1].output.heat_capacities, + [ + 0.0, + 0.0011868003840771405, + 0.028034154256417927, + 0.2009551099977339, + 0.6485781063419351, + 1.3595431369205297, + 2.255468264450785, + 3.266343406515264, + 4.345030847460179, + 5.459276810546103, + 6.584049337246838, + 7.698421094730239, + 8.785103204310895, + 9.83066401753151, + 10.825552552619895, + 11.763768034484103, + 12.64229266802984, + 13.460442497125165, + 14.219247130533086, + 14.920916094983887, + 15.568410803255496, + 16.165118665117106, + 16.7146154410529, + 17.220498811453748, + 17.686276808699812, + 18.11529706060478, + 18.510705522893172, + 18.87542595398425, + 19.212153576086067, + 19.523358122597617, + 19.81129282657783, + 20.0780069227314, + 20.32535998475757, + 20.555036963001697, + 20.768563175541495, + 20.96731877972719, + 21.15255244185952, + 21.325394053424603, + 21.486866430302257, + 21.63789598916472, + 21.779322431952757, + 21.911907491358473, + 22.03634280229226, + 22.15325696964408, + 22.263221903570138, + 22.36675849165218, + 22.464341673696424, + 22.556404980421142, + 22.643344592332866, + 22.7255229700477, + 22.803272102387663, + 22.876896413903097, + 22.94667536911254, + 23.012865806745612, + 23.07570403363029, + 23.13540770457687, + 23.19217751165858, + 23.246198703651245, + 23.29764245404487, + 23.346667093953943, + 23.393419224402656, + 23.438034720823556, + 23.480639641159293, + 23.521351047676305, + 23.560277751467613, + 23.597520987622158, + 23.633175028154362, + 23.667327739007018, + 23.700061086750253, + 23.731451599988663, + 23.761570789948593, + 23.790485534238893, + 23.81825842735472, + 23.844948101117776, + 23.870609517913, + 23.89529423928512, + 23.919050672195343, + 23.941924295003787, + 23.96395786503493, + 23.98519160939703, + 24.00566340056069, + 24.025408918053603, + 24.044461797496005, + 24.062853768083, + 24.080614779513724, + 24.09777311927245, + 24.11435552108127, + 24.13038726526754, + 24.145892271720452, + 24.160893186049147, + 24.175411459499095, + 24.18946742313331, + 24.20308035673952, + 24.21626855288366, + 24.229049376493077, + 24.241439320319376, + 24.253454056600646, + 24.26510848521541, + 24.276416778595795, + 24.287392423644828, + ], + rtol=1e-5, + atol=1e-5, + ) + + assert_allclose( + responses[job.jobs[-1].uuid][1].output.internal_energies, + [ + 5996.475041164254, + 5996.477713192104, + 5996.571519393219, + 5997.522338549997, + 6001.526548372723, + 6011.378319476532, + 6029.332111781312, + 6056.868075013315, + 6094.883409774392, + 6143.886487177431, + 6204.103612017741, + 6275.532406879701, + 6357.97930247125, + 6451.096963612934, + 6554.423277575255, + 6667.418772039798, + 6789.499409913215, + 6920.063186840122, + 7058.510312832553, + 7204.257595003152, + 7356.748007670027, + 7515.456501542569, + 7679.893008919433, + 7849.603442786817, + 8024.169318723846, + 8203.206475764464, + 8386.363245452072, + 8573.318318240243, + 8763.778480239047, + 8957.476336882552, + 9154.16809917656, + 9353.631479073356, + 9555.663720188079, + 9760.079776153992, + 9966.710639614892, + 10175.401818834645, + 10386.01195518593, + 10598.411572652747, + 10812.48194943522, + 11028.114101419309, + 11245.20786741417, + 11463.671086488877, + 11683.41885833677, + 11904.372878276388, + 12126.460839208761, + 12349.615893555434, + 12573.77616887825, + 12798.884331517582, + 13024.887193173528, + 13251.735355892522, + 13479.382891409346, + 13707.787051234469, + 13936.908004270792, + 14166.708599096673, + 14397.154148366308, + 14628.212233058292, + 14859.85252455166, + 15092.04662272923, + 15324.76790850378, + 15557.991409336022, + 15791.693676467215, + 16025.852672725627, + 16260.447669887044, + 16495.459154676966, + 16730.868742597315, + 16966.65909884522, + 17202.81386566677, + 17439.317595555447, + 17676.15568976463, + 17913.31434165674, + 18150.780484458755, + 18388.541743036425, + 18626.586389336906, + 18864.90330118366, + 19103.481924137446, + 19342.312236164504, + 19581.384714877116, + 19820.690307133795, + 20060.22040080573, + 20299.966798533835, + 20539.921693316566, + 20780.077645783087, + 21020.427563019086, + 21260.96467882429, + 21501.682535291275, + 21742.574965604603, + 21983.636077967978, + 22224.86024057487, + 22466.242067545318, + 22707.77640575772, + 22949.4583225106, + 23191.283093954506, + 23433.24619423888, + 23675.34328532347, + 23917.570207407614, + 24159.922969934454, + 24402.397743130532, + 24644.990850044236, + 24887.69875904936, + 25130.51807678263, + ], + rtol=1e-5, + atol=1e-5, + ) + assert responses[job.jobs[-1].uuid][1].output.chemsys == "Si"