diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 41c7bdf320..692f846393 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -60,16 +60,21 @@ jobs: run: | micromamba install -n a2 -c conda-forge enumlib packmol bader --yes + - name: Install ALA-Mode + run: | + micromamba activate a2 + micromamba install -n a2 -c conda-forge alm + - name: Install dependencies run: | micromamba activate a2 python -m pip install --upgrade pip 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,approxneb,aims] + uv pip install .[strict,strict-forcefields,tests,abinit,approxneb,aims,pheasy] uv pip install torch-runstats torch_dftd uv pip install --no-deps nequip==0.5.6 - + - name: Install pymatgen from master if triggered by pymatgen repo dispatch if: github.event_name == 'repository_dispatch' && github.event.action == 'pymatgen-ci-trigger' run: | diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index 737c252eb4..e723907004 100644 --- a/docs/user/codes/vasp.md +++ b/docs/user/codes/vasp.md @@ -342,9 +342,39 @@ phonon_flow = PhononMaker(min_length=15.0, store_force_constants=False).make( ) ``` +#### Pheasy + +Alternatively, users can accelerate the calculation of interatomic force constants using the machine-learning-based [Pheasy code](https://doi.org/10.48550/arXiv.2508.01020). +`Pheasy` can be installed with `pip install pheasy`. +By design, these workflows have the same basic structure as the harmonic forcefield workflows and use [Phonopy](https://doi.org/10.7566/JPSJ.92.012001) in part to compute the phonon spectrum. +To use `Pheasy` in the previous example, we would replace the import string to `from atomate2.vasp.flows.pheasy import PhononMaker`. + +By default, this workflow does not compute anharmonic force constants, but can be extended to using the `cal_anhar_fcs` kwarg and the `ALAMODE` code. + +To install ALAMODE, see their [installation guidelines](https://alamode.readthedocs.io/en/latest/install.html#). +Linux and MacOS x86-64 users can try to install using conda forge: +``` +conda install -c conda-forge alm +``` + +Windows and MacOS ARM users cannot use the pre-built wheels on conda forge at this time, and should instead try the following installation using `conda`: +``` +conda install -c conda-forge "numpy<=2.2" scipy h5py compilers “libblas=*=*mkl” spglib boost eigen cmake ipython mkl-include openmpi --yes +git clone https://github.com/ttadano/ALM.git +cd ALM/python +python setup.py build +pip install -e . +``` +NB: MacOS users will need to ensure that `gcc` and `g++` are used rather than `clang` - both can be installed with `homebrew`. +Note also that `boost` and `eigen` can be installed via `homebrew`. +For example, using `gcc-15` from `homebrew`, one might set: +``` +export CC=gcc-15 ; CXX=g++-15 ; CXX_FLAGS=-DOPENMP +``` + ### Grüneisen parameter workflow -Calculates mode-dependent Grüneisen parameters with the help of Phonopy. +Calculates mode-dependent Grüneisen parameters with the help of [Phonopy](https://doi.org/10.7566/JPSJ.92.012001). Initially, a tight structural relaxation is performed to obtain a structure without forces on the atoms. The optimized structure (ground state) is further expanded and diff --git a/pyproject.toml b/pyproject.toml index 7c4655901d..95128c7c8b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,6 +45,7 @@ amset = ["amset>=0.4.15", "pydash"] cclib = ["cclib>=1.8.1"] mp = ["mp-api>=0.37.5"] phonons = ["phonopy>=1.10.8,<2.43.3", "seekpath>=2.0.0"] # issue writing CIFs because of bug in phonopy >=2.43.3 +pheasy = ["hiphive==1.3.1", "numpy<=2.2", "pheasy==0.0.2"] lobster = ["ijson>=3.2.2", "lobsterpy>=0.4.0"] defects = [ "dscribe>=1.2.0", diff --git a/src/atomate2/common/flows/anharmonicity.py b/src/atomate2/common/flows/anharmonicity.py index 97899de6cb..d416dda801 100644 --- a/src/atomate2/common/flows/anharmonicity.py +++ b/src/atomate2/common/flows/anharmonicity.py @@ -10,7 +10,7 @@ from typing import TYPE_CHECKING from warnings import warn -from jobflow import Flow, Maker +from jobflow import Flow, Maker, Response, job from atomate2.common.jobs.anharmonicity import ( displace_structure, @@ -160,6 +160,7 @@ def make( jobs = [phonon_flow, anharmon_flow, results] return Flow(jobs, results.output) + @job def make_from_phonon_doc( self, phonon_doc: PhononBSDOSDoc, @@ -266,7 +267,8 @@ def make_from_phonon_doc( jobs.append(sigma_calcs) sigma_a_vals = sigma_calcs.output - return Flow(jobs, sigma_a_vals) + flow = Flow(jobs, sigma_a_vals) + return Response(replace=flow, output=flow.output) @property @abstractmethod diff --git a/src/atomate2/common/flows/pheasy.py b/src/atomate2/common/flows/pheasy.py new file mode 100644 index 0000000000..0d1bd474f6 --- /dev/null +++ b/src/atomate2/common/flows/pheasy.py @@ -0,0 +1,362 @@ +"""Flow for calculating (an)harmonic FCs and phonon energy renorma. with pheasy.""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Literal + +from atomate2.common.flows.phonons import BasePhononMaker as PurePhonopyMaker +from atomate2.common.jobs.pheasy import ( + generate_frequencies_eigenvectors, + generate_phonon_displacements, + get_supercell_size, +) +from atomate2.common.jobs.phonons import run_phonon_displacements + +if TYPE_CHECKING: + from pathlib import Path + + from emmet.core.math import Matrix3D + from jobflow import Flow, Job + from pymatgen.core.structure import Structure + + from atomate2.aims.jobs.base import BaseAimsMaker + from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker + from atomate2.vasp.jobs.base import BaseVaspMaker + +SUPPORTED_CODES = frozenset(("vasp", "aims", "forcefields")) + + +@dataclass +class BasePhononMaker(PurePhonopyMaker, ABC): + """Maker to calculate harmonic phonons with LASSO-based ML code Pheasy. + + Calculate the zero-K harmonic phonons of a material and higher-order FCs. + Initially, a tight structural relaxation is performed to obtain a structure + without forces on the atoms. Subsequently, supercells with all atoms displaced + by a small amplitude (generally using 0.01 A) are generated and accurate forces + are computed for these structures for the second order force constants. With the + help of pheasy (LASSO technique), these forces are then converted into a dynamical + matrix. In this Workflow, we separate the harmonic phonon calculations and + anharmonic force constants calculations. To correct for polarization effects, a + correction of the dynamical matrix based on BORN charges can be performed. Finally, + phonon densities of states, phonon band structures and thermodynamic properties + are computed. For the anharmonic force constants, the supercells with all atoms + displaced by a larger amplitude (generally using 0.08 A) are generated and accurate + forces are computed for these structures. With the help of pheasy (LASSO technique), + the third- and fourth-order force constants are extracted at once. + + .. Note:: + It is heavily recommended to symmetrize the structure before passing it to + this flow. Otherwise, a different space group might be detected and too many + displacement calculations will be required for pheasy phonon calculation. It + is recommended to check the convergence parameters here and adjust them if + necessary. The default might not be strict enough for your specific case. + Additionally, for high-throughoput calculations, it is recommended to calculate + the residual forces on the atoms in the supercell after the relaxation. Then the + forces on displaced supercells can deduct the residual forces to reduce the + error in the dynamical matrix. + + Parameters + ---------- + name : str + Name of the flow produced by this maker. + sym_reduce : bool + Whether to reduce the number of deformations using symmetry. + symprec : float + Symmetry precision to use in the + reduction of symmetry to find the primitive/conventional cell + (use_primitive_standard_structure, use_conventional_standard_structure) + and to handle all symmetry-related tasks in pheasy, we recommend to + use the value of 1e-3. + displacement: float + displacement distance for phonons, for most cases 0.01 A is a good choice, + but it can be increased to 0.02 A for heavier elements. + num_displaced_supercells: int + number of displacements to be generated using a random-displacement approach + for harmonic phonon calculations. The default value is 0 and the number of + displacements is automatically determined by the number of atoms in the + supercell and its space group. + cal_anhar_fcs: bool + if set to True, anharmonic force constants(FCs) up to fourth-order FCs will + be calculated. The default value is False, and only harmonic phonons will + be calculated. + displacement_anhar: float + displacement distance for anharmonic force constants(FCs) up to fourth-order + FCs, for most cases 0.08 A is a good choice, but it can be increased to 0.1 A. + num_disp_anhar: int + number of displacements to be generated using a random-displacement approach + for anharmonic phonon calculations. The default value is 0 and the number of + displacements is automatically determined by the number of atoms in the + supercell, cutoff distance for anharmonic FCs its space group. generally, + 50 large-distance displacements are enough for most cases. + fcs_cutoff_radius: list + cutoff distance for anharmonic force constants(FCs) up to fourth-order FCs. + The default value is [-1, 12, 10], which means that the cutoff distance for + second-order FCs is the Wigner-Seitz cell boundary and the cutoff distance + for third-order FCs is 12 Borh, and the cutoff distance for fourth-order FCs + is 10 Bohr. Generally, the default value is good enough. + min_length: float + minimum length of lattice constants will be used to create the supercell, + the default value is 14.0 A. In most cases, the default value is good + enough, but it can be increased for larger supercells. + prefer_90_degrees: bool + if set to True, supercell algorithm will first try to find a supercell + with 3 90 degree angles. + get_supercell_size_kwargs: dict + kwargs that will be passed to get_supercell_size to determine supercell size + use_symmetrized_structure: str + allowed strings: "primitive", "conventional", None + + - "primitive" will enforce to start the phonon computation + from the primitive standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + This makes it possible to use certain k-path definitions + with this workflow. Otherwise, we must rely on seekpath + - "conventional" will enforce to start the phonon computation + from the conventional standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + We will however use seekpath and primitive structures + as determined by from phonopy to compute the phonon band structure + bulk_relax_maker: .ForceFieldRelaxMaker, .BaseAimsMaker, .BaseVaspMaker, or None + A maker to perform a tight relaxation on the bulk. + Set to ``None`` to skip the + bulk relaxation + static_energy_maker: .ForceFieldRelaxMaker, .BaseAimsMaker, .BaseVaspMaker, or None + A maker to perform the computation of the DFT energy on the bulk. + Set to ``None`` to skip the + static energy computation + born_maker: .ForceFieldStaticMaker, .BaseAsimsMaker, .BaseVaspMaker, or None + Maker to compute the BORN charges. + phonon_displacement_maker: .ForceFieldStaticMaker, .BaseAimsMaker, .BaseVaspMaker + Maker used to compute the forces for a supercell. + generate_frequencies_eigenvectors_kwargs : dict + Keyword arguments passed to :obj:`generate_frequencies_eigenvectors`. + create_thermal_displacements: bool + Bool that determines if thermal_displacement_matrices are computed + kpath_scheme: str + scheme to generate kpoints. Please be aware that + you can only use seekpath with any kind of cell + Otherwise, please use the standard primitive structure + Available schemes are: + "seekpath", "hinuma", "setyawan_curtarolo", "latimer_munro". + "seekpath" and "hinuma" are the same definition but + seekpath can be used with any kind of unit cell as + it relies on phonopy to handle the relationship + to the primitive cell and not pymatgen + code: str + determines the dft or force field code. + store_force_constants: bool + if True, force constants will be stored + socket: bool + If True, use the socket for the calculation + + """ + + name: str = "phonon" + sym_reduce: bool = True + symprec: float = 1e-3 + displacement: float = 0.01 + num_displaced_supercells: int = 0 + cal_anhar_fcs: bool = False + displacement_anhar: float = 0.08 + num_disp_anhar: int = 0 + fcs_cutoff_radius: list = field( + default_factory=lambda: [-1, 12, 10] + ) # units in Bohr + renorm_phonon: bool = False + renorm_temp: list = field(default_factory=lambda: [100, 700, 100]) + cal_ther_cond: bool = False + ther_cond_mesh: list = field(default_factory=lambda: [20, 20, 20]) + ther_cond_temp: list = field(default_factory=lambda: [100, 700, 100]) + min_length: float | None = 8.0 + max_atoms: float | None = 200 + force_90_degrees: bool = True + force_diagonal: bool = True + get_supercell_size_kwargs: dict = field(default_factory=dict) + use_symmetrized_structure: Literal["primitive", "conventional"] | None = None + bulk_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | BaseAimsMaker | None = None + static_energy_maker: ForceFieldRelaxMaker | BaseVaspMaker | BaseAimsMaker | None = ( + None + ) + born_maker: ForceFieldStaticMaker | BaseVaspMaker | None = None + phonon_displacement_maker: ForceFieldStaticMaker | BaseVaspMaker | BaseAimsMaker = ( + None + ) + create_thermal_displacements: bool = False + 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 get_displacements( + self, structure: Structure, supercell_matrix: Matrix3D + ) -> Job | Flow: + """ + Get displaced supercells. + + Parameters + ---------- + structure: Structure + supercell_matrix: Matrix3D + + Returns + ------- + Job|Flow + """ + return generate_phonon_displacements( + structure=structure, + supercell_matrix=supercell_matrix, + displacement=self.displacement, + num_displaced_supercells=self.num_displaced_supercells, + cal_anhar_fcs=self.cal_anhar_fcs, + displacement_anhar=self.displacement_anhar, + num_disp_anhar=self.num_disp_anhar, + fcs_cutoff_radius=self.fcs_cutoff_radius, + sym_reduce=self.sym_reduce, + symprec=self.symprec, + use_symmetrized_structure=self.use_symmetrized_structure, + kpath_scheme=self.kpath_scheme, + code=self.code, + ) + + def run_displacements( + self, + displacements: Job | Flow, + prev_dir: str | Path | None, + structure: Structure, + supercell_matrix: Matrix3D, + ) -> Job | Flow: + """ + Perform displacement calculations. + + Parameters + ---------- + displacements: Job | Flow + prev_dir: str | Path | None + structure: Structure + supercell_matrix: Matrix3D + + Returns + ------- + Job | Flow + """ + # perform the phonon displacement calculations + return 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, + store_displaced_structures=True, + ) + + def get_results( + self, + born: Matrix3D, + born_run_job_dir: str, + born_run_uuid: str, + displacement_calcs: Job | Flow, + epsilon_static: Matrix3D, + optimization_run_job_dir: str, + optimization_run_uuid: str, + static_run_job_dir: str, + static_run_uuid: str, + structure: Structure, + supercell_matrix: Matrix3D | None, + total_dft_energy: float, + ) -> Job | Flow: + """ + Calculate the harmonic phonons etc. + + Parameters + ---------- + born: Matrix3D + born_run_job_dir: str + born_run_uuid: str + displacement_calcs: Job | Flow + epsilon_static: Matrix3D + optimization_run_job_dir:str + optimization_run_uuid:str + static_run_job_dir:str + static_run_uuid:str + structure: Structure + supercell_matrix: Matrix3D + total_dft_energy: float + + Returns + ------- + Job | Flow + """ + return generate_frequencies_eigenvectors( + supercell_matrix=supercell_matrix, + displacement=self.displacement, + cal_anhar_fcs=self.cal_anhar_fcs, + fcs_cutoff_radius=self.fcs_cutoff_radius, + renorm_phonon=self.renorm_phonon, + renorm_temp=self.renorm_temp, + cal_ther_cond=self.cal_ther_cond, + ther_cond_mesh=self.ther_cond_mesh, + ther_cond_temp=self.ther_cond_temp, + 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, + **self.generate_frequencies_eigenvectors_kwargs, + ) + + def get_supercell_matrix(self, structure: Structure) -> Job | Flow: + """ + Get supercell matrix. + + Parameters + ---------- + structure: Structure + + Returns + ------- + Job | Flow + """ + return get_supercell_size( + structure, + self.min_length, + self.max_atoms, + self.force_90_degrees, + self.force_diagonal, + ) + + @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/flows/phonons.py b/src/atomate2/common/flows/phonons.py index 9e6a21c90b..1469ecd48b 100644 --- a/src/atomate2/common/flows/phonons.py +++ b/src/atomate2/common/flows/phonons.py @@ -23,6 +23,7 @@ from typing import Literal from emmet.core.math import Matrix3D + from jobflow import Job from pymatgen.core.structure import Structure from atomate2.aims.jobs.base import BaseAimsMaker @@ -268,14 +269,7 @@ def make( # 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=structure, - min_length=self.min_length, - max_length=self.max_length, - prefer_90_degrees=self.prefer_90_degrees, - allow_orthorhombic=self.allow_orthorhombic, - **self.get_supercell_size_kwargs, - ) + supercell_job = self.get_supercell_matrix(structure) jobs.append(supercell_job) supercell_matrix = supercell_job.output @@ -306,27 +300,12 @@ def make( total_dft_energy = compute_total_energy_job.output # get a phonon object from phonopy - displacements = generate_phonon_displacements( - structure=structure, - 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, - ) + displacements = self.get_displacements(structure, supercell_matrix) 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, + displacement_calcs = self.run_displacements( + displacements, prev_dir, structure, supercell_matrix ) jobs.append(displacement_calcs) @@ -349,7 +328,85 @@ def make( born_run_job_dir = born_job.output.dir_name born_run_uuid = born_job.output.uuid - phonon_collect = generate_frequencies_eigenvectors( + phonon_collect = self.get_results( + born, + born_run_job_dir, + born_run_uuid, + displacement_calcs, + epsilon_static, + optimization_run_job_dir, + optimization_run_uuid, + static_run_job_dir, + static_run_uuid, + structure, + supercell_matrix, + total_dft_energy, + ) + + jobs.append(phonon_collect) + + # create a flow including all jobs for a phonon computation + return Flow(jobs, phonon_collect.output) + + def get_supercell_matrix(self, structure: Structure) -> Job | Flow: + """ + Get supercell matrix. + + Parameters + ---------- + structure: Structure + + Returns + ------- + Job|Flow + """ + return get_supercell_size( + structure=structure, + min_length=self.min_length, + max_length=self.max_length, + prefer_90_degrees=self.prefer_90_degrees, + allow_orthorhombic=self.allow_orthorhombic, + **self.get_supercell_size_kwargs, + ) + + def get_results( + self, + born: Matrix3D, + born_run_job_dir: str, + born_run_uuid: str, + displacement_calcs: Job | Flow, + epsilon_static: Matrix3D, + optimization_run_job_dir: str, + optimization_run_uuid: str, + static_run_job_dir: str, + static_run_uuid: str, + structure: Structure, + supercell_matrix: Matrix3D | None, + total_dft_energy: float, + ) -> Job | Flow: + """ + Calculate the harmonic phonons etc. + + Parameters + ---------- + born: Matrix3D + born_run_job_dir: str + born_run_uuid: str + displacement_calcs: Job | Flow + epsilon_static: Matrix3D + optimization_run_job_dir:str + optimization_run_uuid:str + static_run_job_dir:str + static_run_uuid:str + structure: Structure + supercell_matrix: Matrix3D + total_dft_energy: float + + Returns + ------- + Job | Flow + """ + return generate_frequencies_eigenvectors( supercell_matrix=supercell_matrix, displacement=self.displacement, sym_reduce=self.sym_reduce, @@ -373,10 +430,62 @@ def make( **self.generate_frequencies_eigenvectors_kwargs, ) - jobs.append(phonon_collect) + def run_displacements( + self, + displacements: Job | Flow, + prev_dir: str | Path | None, + structure: Structure, + supercell_matrix: Matrix3D, + ) -> Job | Flow: + """ + Perform displacement calculations. - # create a flow including all jobs for a phonon computation - return Flow(jobs, phonon_collect.output) + Parameters + ---------- + displacements: Job | Flow + prev_dir: str | Path | None + structure: Structure + supercell_matrix: Matrix3D + + Returns + ------- + Job | Flow + """ + return 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, + ) + + def get_displacements( + self, structure: Structure, supercell_matrix: Matrix3D + ) -> Job | Flow: + """ + Get displaced supercells. + + Parameters + ---------- + structure: Structure + supercell_matrix: Matrix3D + + Returns + ------- + Job|Flow + """ + return generate_phonon_displacements( + structure=structure, + 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, + ) @property @abstractmethod diff --git a/src/atomate2/common/jobs/anharmonicity.py b/src/atomate2/common/jobs/anharmonicity.py index 5ebace0e9e..9416f9db2d 100644 --- a/src/atomate2/common/jobs/anharmonicity.py +++ b/src/atomate2/common/jobs/anharmonicity.py @@ -22,6 +22,8 @@ from pathlib import Path from typing import Any + from emmet.core.math import Matrix3D + from atomate2.aims.jobs.base import BaseAimsMaker from atomate2.common.schemas.phonons import ForceConstants, PhononBSDOSDoc from atomate2.forcefields.jobs import ForceFieldStaticMaker @@ -338,7 +340,7 @@ def displace_structure( def build_dynmat( - force_constants: ForceConstants, + force_constants: ForceConstants | list[list[Matrix3D]], structure: Structure, ) -> np.ndarray: """Build the dynamical matrix. @@ -356,7 +358,7 @@ def build_dynmat( The dynamical matrix """ force_constants_2d = ( - np.array(force_constants.force_constants) + np.array(getattr(force_constants, "force_constants", force_constants)) .swapaxes(1, 2) .reshape(2 * (len(structure) * 3,)) ) @@ -449,7 +451,7 @@ def get_sigma_a_per_mode( @job def get_forces( - force_constants: ForceConstants, + force_constants: ForceConstants | list[list[Matrix3D]], phonon_supercell: Structure, displaced_structures: dict[str, list], ) -> list[np.ndarray]: @@ -470,7 +472,7 @@ def get_forces( List of forces in the form [DFT forces, harmonic forces] """ force_constants_2d = np.swapaxes( - force_constants.force_constants, + getattr(force_constants, "force_constants", force_constants), 1, 2, ).reshape(2 * (len(phonon_supercell) * 3,)) diff --git a/src/atomate2/common/jobs/gruneisen.py b/src/atomate2/common/jobs/gruneisen.py index d78fceb636..8a67a9ee31 100644 --- a/src/atomate2/common/jobs/gruneisen.py +++ b/src/atomate2/common/jobs/gruneisen.py @@ -5,6 +5,7 @@ import logging from typing import TYPE_CHECKING +from emmet.core.phonon import PhononBSDOSDoc from jobflow import Flow, Response, job from pymatgen.phonon.gruneisen import ( GruneisenParameter, @@ -14,7 +15,6 @@ from atomate2 import SETTINGS from atomate2.common.schemas.gruneisen import GruneisenParameterDocument -from atomate2.common.schemas.phonons import PhononBSDOSDoc if TYPE_CHECKING: from pathlib import Path @@ -85,8 +85,7 @@ def run_phonon_jobs( set_symmetry = list(set(symmetry)) if len(set_symmetry) == 1: jobs = [] - phonon_yaml_dirs = dict.fromkeys(("ground", "plus", "minus"), None) - phonon_imaginary_modes = dict.fromkeys(("ground", "plus", "minus"), None) + _for_post_process = {} for st, struct in opt_struct.items(): # phonon run for all 3 optimized structures (ground state, expanded, shrunk) phonon_kwargs = {} @@ -107,15 +106,12 @@ def run_phonon_jobs( ) jobs.append(phonon_job) # store each phonon run task doc - phonon_yaml_dirs[st] = phonon_job.output.jobdirs.taskdoc_run_job_dir - phonon_imaginary_modes[st] = phonon_job.output.has_imaginary_modes + _for_post_process[st] = phonon_job.output + processed = get_calc_meta(_for_post_process) return Response( - replace=Flow(jobs), - output={ - "phonon_yaml": phonon_yaml_dirs, - "imaginary_modes": phonon_imaginary_modes, - }, + replace=Flow([*jobs, processed]), + output=processed.output, ) logger.warning( msg="Different space groups were detected for the optimized structures." @@ -124,6 +120,27 @@ def run_phonon_jobs( return Response(output={"error": "different space groups"}, stop_jobflow=True) +@job +def get_calc_meta( + phonon_jobs_output: dict[str, PhononBSDOSDoc], +) -> dict[str, dict[str, Path | bool]]: + """Return the metadata associated with a set of phonon calculations.""" + return { + "phonon_yaml": { + label: next( + iter( + [cm.dir_name for cm in pbd.calc_meta if cm.name == "taskdoc_run"] + or [None] + ) + ) + for label, pbd in phonon_jobs_output.items() + }, + "imaginary_modes": { + label: pbd.has_imaginary_modes for label, pbd in phonon_jobs_output.items() + }, + } + + @job( output_schema=GruneisenParameterDocument, data=[GruneisenParameter, GruneisenPhononBandStructureSymmLine], diff --git a/src/atomate2/common/jobs/pheasy.py b/src/atomate2/common/jobs/pheasy.py new file mode 100644 index 0000000000..d8df5565b8 --- /dev/null +++ b/src/atomate2/common/jobs/pheasy.py @@ -0,0 +1,945 @@ +"""Jobs for running phonon calculations with phonopy and pheasy.""" + +from __future__ import annotations + +import logging +import shlex +import subprocess +from pathlib import Path +from typing import TYPE_CHECKING + +import numpy as np +from ase.io import read as ase_read +from emmet.core import __version__ as _emmet_core_version +from emmet.core.phonon import PhononBSDOSDoc +from hiphive import ClusterSpace, ForceConstantPotential, enforce_rotational_sum_rules +from hiphive import ForceConstants as HiPhiveForceConstants +from hiphive.cutoffs import estimate_maximum_cutoff +from hiphive.utilities import extract_parameters +from jobflow import job +from packaging.version import parse as parse_version +from phonopy.file_IO import parse_FORCE_CONSTANTS, write_force_constants_to_hdf5 +from phonopy.interface.vasp import write_vasp +from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections +from phonopy.structure.symmetry import symmetrize_borns_and_epsilon +from pymatgen.core import Structure +from pymatgen.io.phonopy import ( + get_ph_bs_symm_line, + get_ph_dos, + get_phonopy_structure, + get_pmg_structure, +) +from pymatgen.io.vasp import Kpoints +from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine +from pymatgen.phonon.dos import PhononDos +from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter +from pymatgen.transformations.advanced_transformations import ( + CubicSupercellTransformation, +) + +from atomate2.common.jobs.phonons import _generate_phonon_object, _get_kpath + +if TYPE_CHECKING: + from emmet.core.math import Matrix3D + +logger = logging.getLogger(__name__) + +try: + from alm import ALM +except ImportError: + ALM = None + +_DEFAULT_FILE_PATHS = { + "force_displacements": "dataset_forces.npy", + "displacements": "dataset_disps.npy", + "displacements_folded": "dataset_disps_array_rr.npy", + "phonopy": "phonopy.yaml", + "band_structure": "phonon_band_structure.yaml", + "band_structure_plot": "phonon_band_structure.pdf", + "dos": "phonon_dos.yaml", + "dos_plot": "phonon_dos.pdf", + "force_constants": "FORCE_CONSTANTS", + "harmonic_displacements": "disp_matrix.npy", + "anharmonic_displacements": "disp_matrix_anhar.npy", + "harmonic_force_matrix": "force_matrix.npy", + "anharmonic_force_matrix": "force_matrix_anhar.npy", + "website": "phonon_website.json", +} + + +@job +def get_supercell_size( + structure: Structure, + min_length: float, + max_atoms: int, + force_90_degrees: bool, + force_diagonal: bool, +) -> list[list[float]]: + """ + Determine supercell size with given min_length and max_length. + + Parameters + ---------- + structure: Structure Object + Input structure that will be used to determine supercell + min_length: float + minimum length of cell in Angstrom + max_length: float + maximum length of cell in Angstrom + prefer_90_degrees: bool + if True, the algorithm will try to find a cell with 90 degree angles first + allow_orthorhombic: bool + if True, orthorhombic supercells are allowed + **kwargs: + Additional parameters that can be set. + """ + transformation = CubicSupercellTransformation( + min_length=min_length, + max_atoms=max_atoms, + force_90_degrees=force_90_degrees, + force_diagonal=force_diagonal, + angle_tolerance=1e-2, + allow_orthorhombic=False, + ) + transformation.apply_transformation(structure=structure) + return transformation.transformation_matrix.transpose().tolist() + + +@job(data=[Structure]) +def generate_phonon_displacements( + structure: Structure, + supercell_matrix: np.array, + displacement: float, + num_displaced_supercells: int, + cal_anhar_fcs: bool, + displacement_anhar: float, + num_disp_anhar: int, + fcs_cutoff_radius: list[int], + sym_reduce: bool, + symprec: float, + use_symmetrized_structure: str | None, + kpath_scheme: str, + code: str, + random_seed: int | None = 103, + verbose: bool = False, +) -> list[Structure]: + """Generate small-distance perturbed structures with phonopy based on two ways. + + (we will directly use the pheasy to generate the supercell in the near future) + 1. finite-displacment method (one displaced atom) when the displacement number + is less than 3. 2. random-displacement method (all-displaced atoms) when the + displacement number is more than 3. + + 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 (default: 0.01) + num_displaced_supercells: int + number of displaced supercells defined by users + cal_anhar_fcs: bool + TODO : docstr + displacement_anhar: float + TODO : docstr + 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: str + code to perform the computations + random_seed : int | None = 103 + Random seed to use in generating randomly-displaced structures. + verbose : bool = False + Whether to log warnings. + + """ + # TODO: remove ALMODE dependence for 2nd order force constants + if not ALM: + raise ImportError( + "Error importing ALM. Please ensure the 'alm' library is installed." + ) + phonon = _generate_phonon_object( + structure, + supercell_matrix, + displacement, + sym_reduce, + symprec, + use_symmetrized_structure, + kpath_scheme, + code, + verbose=verbose, + ) + + # 1. the ALM module is used to determine the number of free parameters + # (irreducible force constants) corresponding to the second order + # force constants (FCs) given a supercell. + # 2. Based on the number of free parameters, we can determine how many + # displaced supercells we need to use to extract the second order force + # constants. Generally, the number of free parameters should be less than + # 3 * natom(supercell) * num_displaced_supercells. However, the full rank + # of the matrix can not always guarantee accurate results, you + # may need to displace more random configurations. Use at least one or + # two more configurations based on the suggested number of displacements. + supercell_ph = phonon.supercell + lattice = supercell_ph.cell + positions = supercell_ph.scaled_positions + numbers = supercell_ph.numbers + natom = len(numbers) + + # get the number of free parameters of 2ND FCs from ALM, labeled as n_fp + with ALM(lattice, positions, numbers) as alm: + alm.define(1) + alm.suggest() + n_fp = alm._get_number_of_irred_fc_elements(1) # noqa: SLF001 + + # get the number of displaced supercells based on the number of free parameters + num_disp_sc = int(np.ceil(n_fp / (3.0 * natom))) + + if verbose: + logger.info( + f"There are {n_fp} free parameters for the second-order " + "force constants (FCs)." + f"There are {3 * natom * num_disp_sc} equations used to " + "obtain the second-order FCs." + "CAUTION: you may need to increase the number of " + "displacements in some cases." + "If the number of atoms in the supercell are less than 100 and " + "all lattice constants are less than 10 Å, the user is advised " + "to use 1-2 more randomly-displaced configurations." + ) + + # get the number of displaced supercells from phonopy to compared with the number + # of 3, if the number of displaced supercells is less than 3, we will use the finite + # displacement method to generate the supercells. Otherwise, we will use the random + # displacement method to generate the supercells. + if len(phonon.displacements) > 3: + phonon.generate_displacements( + distance=displacement, + number_of_snapshots=( + num_displaced_supercells + if num_displaced_supercells != 0 + else int(np.ceil(num_disp_sc * 1.8)) + 1 + ), + random_seed=random_seed, + ) + + supercells = phonon.supercells_with_displacements + displacements = [get_pmg_structure(cell) for cell in supercells] + + # Here, the ALAMODE code is used to determine the number of + # third and fourth-order FCs are needed for the supercell + if cal_anhar_fcs: + # Due to the cutoff radius of the force constants use the unit of Borh in ALM, + # we need to convert the cutoff radius from Angstrom to Bohr. + with ALM(lattice * 1.89, positions, numbers) as alm: + # Define the force constants up to fourth order with a list of + # cutoff radius + alm.define(3, fcs_cutoff_radius) + # Perform symmetry analysis and suggest irreducible force constants. + alm.suggest() + # Get the number of irreducible elements for both 3RD- and 4TH-order + # force constants + n_rd_anh = alm._get_number_of_irred_fc_elements( # noqa: SLF001 + 2 + ) + alm._get_number_of_irred_fc_elements(3) # noqa: SLF001 + # we can determine how many displaced supercells we need to use to extract + # the 3rd and 4th order force constants, and we can add a scaling factor + # to reduce the number of displaced supercells due to we use the lasso + # technique. + num_d_anh = int(np.ceil(n_rd_anh / (3.0 * natom))) + num_dis_cells_anhar = num_disp_anhar if num_disp_anhar != 0 else num_d_anh + + num_dis_cells_anhar = 20 + # generate the supercells for anharmonic force constants + phonon.generate_displacements( + distance=displacement_anhar, + number_of_snapshots=num_dis_cells_anhar, + random_seed=random_seed, + ) + supercells = phonon.supercells_with_displacements + displacements += [get_pmg_structure(cell) for cell in supercells] + + # add the equilibrium structure to the list for calculating + # the residual forces. + displacements.append(get_pmg_structure(phonon.supercell)) + return displacements + + +@job( + output_schema=PhononBSDOSDoc, + data=[PhononDos, PhononBandStructureSymmLine, "force_constants"], +) +def generate_frequencies_eigenvectors( + structure: Structure, + supercell_matrix: np.array, + displacement: float, + cal_anhar_fcs: bool, + fcs_cutoff_radius: list[int], + renorm_phonon: bool, + cal_ther_cond: bool, + ther_cond_mesh: list[int], + ther_cond_temp: list[int], + 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, + **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 + verbose : bool = False + Whether to log error messages. + kwargs: dict + Additional parameters that are passed to PhononBSDOSDoc.from_forces_born + """ + phonon = _generate_phonon_object( + structure, + supercell_matrix, + displacement, + sym_reduce, + symprec, + use_symmetrized_structure, + kpath_scheme, + code, + verbose=False, + ) + + # Write the POSCAR and SPOSCAR files for the input of pheasy code + supercell = phonon._supercell # noqa: SLF001 + write_vasp("POSCAR", get_phonopy_structure(structure)) + write_vasp("SPOSCAR", supercell) + + # get the force-displacement dataset from previous calculations + dataset_forces = np.array(displacement_data["forces"]) + np.save(_DEFAULT_FILE_PATHS["force_displacements"], dataset_forces) + + # To deduct the residual forces on an equilibrium structure to eliminate the + # fitting error + dataset_forces_array_rr = dataset_forces - dataset_forces[-1, :, :] + + # force matrix on the displaced structures + dataset_forces_array_disp = dataset_forces_array_rr[:-1, :, :] + + # To handle the large dispalced distance in the dataset + dataset_disps = np.array( + [disps.frac_coords for disps in displacement_data["displaced_structures"]] + ) + np.save(_DEFAULT_FILE_PATHS["displacements"], dataset_disps) + + dataset_disps_array_rr = np.round( + (dataset_disps - supercell.get_scaled_positions()), decimals=16 + ) + np.save(_DEFAULT_FILE_PATHS["displacements_folded"], dataset_disps_array_rr) + + dataset_disps_array_rr = np.where( + dataset_disps_array_rr > 0.5, + dataset_disps_array_rr - 1.0, + dataset_disps_array_rr, + ) + dataset_disps_array_rr = np.where( + dataset_disps_array_rr < -0.5, + dataset_disps_array_rr + 1.0, + dataset_disps_array_rr, + ) + + # Transpose the displacement array on the + # last two axes (atoms and coordinates) + dataset_disps_array_rr_transposed = np.transpose(dataset_disps_array_rr, (0, 2, 1)) + + # Perform matrix multiplication with the transposed supercell.cell + # 'ij' for supercell.cell.T and + # 'nkj' for the transposed dataset_disps_array_rr + dataset_disps_array_rr_cartesian = np.einsum( + "ij,njk->nik", supercell.cell.T, dataset_disps_array_rr_transposed + ) + # Transpose back to the original format + dataset_disps_array_rr_cartesian = np.transpose( + dataset_disps_array_rr_cartesian, (0, 2, 1) + ) + + dataset_disps_array_use = dataset_disps_array_rr_cartesian[:-1, :, :] + + # separate the dataset into harmonic and anharmonic parts + num_har = dataset_disps_array_use.shape[0] + if cal_anhar_fcs: + if not ALM: + raise ImportError( + "Error importing ALM. Please ensure the 'alm' library is installed." + ) + + supercell_ph = phonon.supercell + lattice = supercell_ph.cell + positions = supercell_ph.scaled_positions + numbers = supercell_ph.numbers + natom = len(numbers) + + # get the number of free parameters of 2ND FCs from ALM, labeled as n_fp + with ALM(lattice, positions, numbers) as alm: + alm.define(1) + alm.suggest() + n_fp = alm._get_number_of_irred_fc_elements(1) # noqa: SLF001 + + # get the number of displaced supercells based on the + # number of free parameters + num = int(np.ceil(n_fp / (3.0 * natom))) + + # get the number of displaced supercells from phonopy to compared + # with the number of 3, if the number of displaced supercells is + # less than 3, we will use the finite displacement method to generate + # the supercells. Otherwise, we will use the random displacement + # method to generate the supercells. + num_disp_f = len(phonon.displacements) + num_har = int(np.ceil(num * 1.8)) if num_disp_f > 3 else num_disp_f + + np.save( + _DEFAULT_FILE_PATHS["harmonic_displacements"], + dataset_disps_array_use[:num_har, :, :], + ) + np.save( + _DEFAULT_FILE_PATHS["harmonic_force_matrix"], + dataset_forces_array_disp[:num_har, :, :], + ) + + # get the born charges and dielectric constant + 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 + + prim = ase_read("POSCAR") + supercell = ase_read("SPOSCAR") + + # Create the clusters and orbitals for second order force constants + # For the variables: --w, --nbody, they are used to specify the order of the + # force constants. in the near future, we will add the option to specify the + # order of the force constants. And these two variables can be defined by the + # users. + pheasy_cmd_1 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} " + f"-s -w 2 --symprec {float(symprec)} --nbody 2" + ) + + # Create the null space to further reduce the free parameters for + # specific force constants and make them physically correct. + pheasy_cmd_2 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -c --symprec " + f"{float(symprec)} -w 2" + ) + + # Generate the Compressive Sensing matrix,i.e., displacement matrix + # for the input of machine leaning method.i.e., LASSO, + pheasy_cmd_3 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -w 2 -d " + f"--symprec {float(symprec)} " + f"--ndata {int(num_har)} --disp_file" + ) + + # Here we set a criteria to determine which method to use to generate the + # force constants. If the number of displacements is larger than 3, we + # will use the LASSO method to generate the force constants. Otherwise, + # we will use the least-squred method to generate the force constants. + if len(phonon.displacements) > 3: + # Calculate the force constants using the LASSO method due to the + # random-displacement method Obviously, the rotaional invariance + # constraint, i.e., tag: --rasr BHH, is enforced during the + # fitting process. + pheasy_cmd_4 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -f --full_ifc " + f"-w 2 --symprec {float(symprec)} " + f"-l LASSO --std --rasr BHH --ndata {int(num_har)}" + ) + + else: + # Calculate the force constants using the least-squred method + pheasy_cmd_4 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -f --full_ifc " + f"-w 2 --symprec {float(symprec)} " + f"--rasr BHH --ndata {int(num_har)}" + ) + + logger.info("Start running pheasy in cluster") + + subprocess.call(shlex.split(pheasy_cmd_1)) + subprocess.call(shlex.split(pheasy_cmd_2)) + subprocess.call(shlex.split(pheasy_cmd_3)) + subprocess.call(shlex.split(pheasy_cmd_4)) + + # When this code is run on Github tests, it is failing because it is + # not able to find the FORCE_CONSTANTS file. This is because the file is + # somehow getting generated in some temp directory. Can you fix the bug? + fc_file = Path(_DEFAULT_FILE_PATHS["force_constants"]) + + if cal_anhar_fcs and fc_file.exists(): + np.save( + _DEFAULT_FILE_PATHS["anharmonic_displacements"], + dataset_disps_array_use[num_har:, :, :], + ) + np.save( + _DEFAULT_FILE_PATHS["anharmonic_force_matrix"], + dataset_forces_array_disp[num_har:, :, :], + ) + num_anhar = dataset_forces_array_disp.shape[0] - num_har + + # We next begin to generate the anharmonic force constants up to fourth + # order using the LASSO method + pheasy_cmd_5 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -s -w 4 --symprec " + f"{float(symprec)} " + f"--nbody 2 3 3 --c3 {float(fcs_cutoff_radius[1] / 1.89)} " + f"--c4 {float(fcs_cutoff_radius[2] / 1.89)}" + ) + + pheasy_cmd_6 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -c --symprec " + f"{float(symprec)} -w 4" + ) + pheasy_cmd_7 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -w 4 -d --symprec " + f"{float(symprec)} " + f"--ndata {int(num_anhar)} --disp_file" + ) + pheasy_cmd_8 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -f -w 4 --fix_fc2 " + f"--symprec {float(symprec)} " + f"--ndata {int(num_anhar)} " + ) + + subprocess.call(shlex.split(pheasy_cmd_5)) + subprocess.call(shlex.split(pheasy_cmd_6)) + subprocess.call(shlex.split(pheasy_cmd_7)) + subprocess.call(shlex.split(pheasy_cmd_8)) + + # begin to renormzlize the phonon energies + if renorm_phonon: + pheasy_cmd_9 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -f -w 4 --fix_fc2 " + f"--hdf5 --symprec {float(symprec)} " + f"--ndata {int(num_anhar)}" + ) + + subprocess.call(shlex.split(pheasy_cmd_9)) + + # write the born charges and dielectric constant to the pheasy format + + # begin to convert the force constants to the phonopy and phono3py format + # for the further lattice thermal conductivity calculations + if cal_ther_cond and fc_file.exists(): + # convert the 2ND order force constants to the phonopy format + fc_phonopy_text = parse_FORCE_CONSTANTS(filename=fc_file) + write_force_constants_to_hdf5(fc_phonopy_text, filename="fc2.hdf5") + + # convert the 3RD order force constants to the phonopy format + + prim_hiphive = ase_read("POSCAR") + supercell_hiphive = ase_read("SPOSCAR") + fcs = HiPhiveForceConstants.read_shengBTE( + supercell_hiphive, "FORCE_CONSTANTS_3RD", prim_hiphive + ) + fcs.write_to_phono3py("fc3.hdf5") + + phono3py_cmd = ( + f"phono3py --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} {int(supercell_matrix[2][2])} " + f"--fc2 --fc3 --br --isotope --wigner " + f"--mesh {ther_cond_mesh[0]} {ther_cond_mesh[1]} {ther_cond_mesh[2]} " + f"--tmin {ther_cond_temp[0]} --tmax {ther_cond_temp[1]} " + f"--tstep {ther_cond_temp[2]}" + ) + + subprocess.call(shlex.split(phono3py_cmd)) + + if fc_file.exists(): + # Read the force constants from the output file of pheasy code + force_constants = parse_FORCE_CONSTANTS(filename=fc_file) + 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(_DEFAULT_FILE_PATHS["phonopy"]) + + # get phonon band structure + kpath_dict, kpath_concrete = _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.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_hdf5_band_structure(filename=_DEFAULT_FILE_PATHS["band_structure"]) + phonon.write_yaml_band_structure(filename=_DEFAULT_FILE_PATHS["band_structure"]) + bs_symm_line = get_ph_bs_symm_line( + _DEFAULT_FILE_PATHS["band_structure"], + labels_dict=kpath_dict, + has_nac=born is not None, + ) + + bs_plot_file = kwargs.get("filename_bs", _DEFAULT_FILE_PATHS["band_structure_plot"]) + dos_plot_file = kwargs.get("filename_dos", _DEFAULT_FILE_PATHS["dos_plot"]) + + new_plotter = PhononBSPlotter(bs=bs_symm_line) + new_plotter.save_plot( + filename=bs_plot_file, + 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) + ) + + # If imaginary modes are present, we first use the hiphive code to enforce + # some symmetry constraints to eliminate the imaginary modes (generally work + # for small imaginary modes near Gamma point). If the imaginary modes are + # still present, we will use the pheasy code to generate the force constants + # using a shorter cutoff (10 A) to eliminate the imaginary modes, also we + # just want to remove the imaginary modes near Gamma point. In the future, + # we will only use the pheasy code to do the job. + + if imaginary_modes: + # Define a cluster space using the largest cutoff you can + max_cutoff = estimate_maximum_cutoff(supercell) - 0.01 + cutoffs = [max_cutoff] # only second order needed + cs = ClusterSpace(prim, cutoffs) + + # import the phonopy force constants using the correct supercell also + # provided by phonopy + fcs = HiPhiveForceConstants.read_phonopy(supercell, "FORCE_CONSTANTS") + + # Find the parameters that best fits the force constants given you + # cluster space + parameters = extract_parameters(fcs, cs) + + # Enforce the rotational sum rules + parameters_rot = enforce_rotational_sum_rules( + cs, parameters, ["Huang", "Born-Huang"], alpha=1e-6 + ) + + # use the new parameters to make a fcp and then create the force + # constants and write to a phonopy file + fcp = ForceConstantPotential(cs, parameters_rot) + fcs = fcp.get_force_constants(supercell) + new_fc_file = f"{_DEFAULT_FILE_PATHS['force_constants']}_short_cutoff" + fcs.write_to_phonopy(new_fc_file, format="text") + + force_constants = parse_FORCE_CONSTANTS(filename=new_fc_file) + phonon.force_constants = force_constants + phonon.symmetrize_force_constants() + + phonon.run_band_structure( + qpoints, path_connections=connections, with_eigenvectors=True + ) + phonon.write_yaml_band_structure(filename=_DEFAULT_FILE_PATHS["band_structure"]) + bs_symm_line = get_ph_bs_symm_line( + _DEFAULT_FILE_PATHS["band_structure"], + labels_dict=kpath_dict, + has_nac=born is not None, + ) + + new_plotter = PhononBSPlotter(bs=bs_symm_line) + + new_plotter.save_plot( + filename=bs_plot_file, + units=kwargs.get("units", "THz"), + ) + + imaginary_modes = bs_symm_line.has_imaginary_freq( + tol=kwargs.get("tol_imaginary_modes", 1e-5) + ) + + # Using a shorter cutoff (10 A) to generate the force constants to + # eliminate the imaginary modes near Gamma point in phesay code + if imaginary_modes: + pheasy_cmd_11 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -s -w 2 --c2 " + f"10.0 --symprec {float(symprec)} " + f"--nbody 2" + ) + + pheasy_cmd_12 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -c --symprec " + f"{float(symprec)} --c2 10.0 -w 2" + ) + + pheasy_cmd_13 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -w 2 -d --symprec " + f"{float(symprec)} --c2 10.0 " + f"--ndata {int(num_har)} --disp_file" + ) + + phonon.generate_displacements(distance=displacement) + + if len(phonon.displacements) > 3: + pheasy_cmd_14 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -f --c2 10.0 " + f"--full_ifc -w 2 --symprec {float(symprec)} " + f"-l LASSO --std --rasr BHH --ndata {int(num_har)}" + ) + + else: + pheasy_cmd_14 = ( + f"pheasy --dim {int(supercell_matrix[0][0])} " + f"{int(supercell_matrix[1][1])} " + f"{int(supercell_matrix[2][2])} -f --full_ifc " + f"--c2 10.0 -w 2 --symprec {float(symprec)} " + f"--rasr BHH --ndata {int(num_har)}" + ) + + subprocess.call(shlex.split(pheasy_cmd_11)) + subprocess.call(shlex.split(pheasy_cmd_12)) + subprocess.call(shlex.split(pheasy_cmd_13)) + subprocess.call(shlex.split(pheasy_cmd_14)) + + force_constants = parse_FORCE_CONSTANTS(filename=new_fc_file) + phonon.force_constants = force_constants + phonon.symmetrize_force_constants() + + phonon.save(_DEFAULT_FILE_PATHS["phonopy"]) + + # get phonon band structure + kpath_dict, kpath_concrete = _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 + phonon.run_band_structure( + qpoints, path_connections=connections, with_eigenvectors=True + ) + phonon.write_yaml_band_structure(filename=_DEFAULT_FILE_PATHS["band_structure"]) + bs_symm_line = get_ph_bs_symm_line( + _DEFAULT_FILE_PATHS["band_structure"], + labels_dict=kpath_dict, + has_nac=born is not None, + ) + new_plotter = PhononBSPlotter(bs=bs_symm_line) + + new_plotter.save_plot( + filename=bs_plot_file, + units=kwargs.get("units", "THz"), + ) + + 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(_DEFAULT_FILE_PATHS["website"]) + + # get phonon density of states + 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=_DEFAULT_FILE_PATHS["dos"]) + dos = get_ph_dos(_DEFAULT_FILE_PATHS["dos"]) + new_plotter_dos = PhononDosPlotter() + new_plotter_dos.add_dos(label="total", dos=dos) + new_plotter_dos.save_plot( + filename=dos_plot_file, + units=kwargs.get("units", "THz"), + ) + + # 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 + ) + + cls_constructor = ( + "migrate_fields" + if parse_version(_emmet_core_version) >= parse_version("0.85.1") + else "from_structure" + ) + return getattr(PhononBSDOSDoc, cls_constructor)( + structure=structure, + meta_structure=structure, + phonon_bandstructure=bs_symm_line, + phonon_dos=dos, + total_dft_energy=total_dft_energy_per_formula_unit, + has_imaginary_modes=imaginary_modes, + force_constants=( + {"force_constants": phonon.force_constants.tolist()} + if kwargs.get("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"], + }, + post_process_settings={ + "npoints_band": npoints_band, + "kpath_scheme": kpath_scheme, + "kpoint_density_dos": kpoint_density_dos, + }, + ) diff --git a/src/atomate2/common/jobs/phonons.py b/src/atomate2/common/jobs/phonons.py index 41b7eafb10..2ec5a43f35 100644 --- a/src/atomate2/common/jobs/phonons.py +++ b/src/atomate2/common/jobs/phonons.py @@ -3,24 +3,39 @@ from __future__ import annotations import contextlib +import copy import logging import warnings +from pathlib import Path from typing import TYPE_CHECKING import numpy as np +from emmet.core import __version__ as _emmet_core_version +from emmet.core.phonon import PhononBSDOSDoc from jobflow import Flow, Response, job +from packaging.version import parse as parse_version from phonopy import Phonopy +from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections +from phonopy.structure.symmetry import symmetrize_borns_and_epsilon +from phonopy.units import VaspToTHz from pymatgen.core import Structure -from pymatgen.io.phonopy import get_phonopy_structure, get_pmg_structure +from pymatgen.io.phonopy import ( + get_ph_bs_symm_line, + get_ph_dos, + get_phonopy_structure, + get_pmg_structure, +) +from pymatgen.io.vasp import Kpoints from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine from pymatgen.phonon.dos import PhononDos +from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter +from pymatgen.symmetry.bandstructure import HighSymmKpath +from pymatgen.symmetry.kpath import KPathSeek -from atomate2.common.schemas.phonons import ForceConstants, PhononBSDOSDoc, get_factor +from atomate2.aims.utils.units import omegaToTHz from atomate2.common.utils import get_supercell_matrix if TYPE_CHECKING: - from pathlib import Path - from emmet.core.math import Matrix3D from atomate2.aims.jobs.base import BaseAimsMaker @@ -31,6 +46,67 @@ logger = logging.getLogger(__name__) +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 ["ase", "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_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 + """ + valid_schemes = {"setyawan_curtarolo", "latimer_munro", "hinuma", "seekpath"} + if kpath_scheme in (valid_schemes - {"seekpath"}): + 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=}, must be one of {valid_schemes}") + + 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 + + @job def get_total_energy_per_cell( total_dft_energy_per_formula_unit: float, structure: Structure @@ -90,8 +166,7 @@ def get_supercell_size( ) -@job(data=[Structure]) -def generate_phonon_displacements( +def _generate_phonon_object( structure: Structure, supercell_matrix: np.array, displacement: float, @@ -100,9 +175,9 @@ def generate_phonon_displacements( use_symmetrized_structure: str | None, kpath_scheme: str, code: str, -) -> list[Structure]: - """ - Generate displaced structures with phonopy. + verbose: bool = False, +) -> Phonopy: + """Bundle commonly-used Phonopy object construction. Parameters ---------- @@ -122,31 +197,38 @@ def generate_phonon_displacements( scheme to generate kpath code: code to perform the computations + use_min_dof : bool = False + Whether to use the minimal number of degrees of freedom + in calculating randomly-displaced structures. + Requires ALAMODE. + verbose : bool = False + Whether to log error messages. Returns ------- - List[Structure] - Displaced structures + Phonopy object. """ - warnings.warn( - "Initial magnetic moments will not be considered for the determination " - "of the symmetry of the structure and thus will be removed now.", - stacklevel=2, + if "magmom" in structure.site_properties and verbose: + warnings.warn( + "Initial magnetic moments will not be considered for the determination " + "of the symmetry of the structure and thus will be removed now.", + stacklevel=2, + ) + + cell = get_phonopy_structure( + structure.copy().remove_site_property(property_name="magmom") + if "magmom" in structure.site_properties + else structure.copy() ) - if "magmom" in structure.site_properties: - # remove_site_property is in-place so make a structure copy first - no_mag_struct = structure.copy().remove_site_property(property_name="magmom") - else: - no_mag_struct = structure - cell = get_phonopy_structure(no_mag_struct) factor = get_factor(code) # a bit of code repetition here as I currently # do not see how to pass the phonopy object? - if use_symmetrized_structure == "primitive" and kpath_scheme != "seekpath": - primitive_matrix: np.ndarray | str = np.eye(3) - else: - primitive_matrix = "auto" + primitive_matrix: np.ndarray | str = ( + np.eye(3) + if use_symmetrized_structure == "primitive" and kpath_scheme != "seekpath" + else "auto" + ) # TARP: THIS IS BAD! Including for discussions sake if cell.magnetic_moments is not None and primitive_matrix == "auto": @@ -166,15 +248,72 @@ def generate_phonon_displacements( is_symmetry=sym_reduce, ) phonon.generate_displacements(distance=displacement) + return phonon - supercells = phonon.supercells_with_displacements +@job(data=[Structure]) +def generate_phonon_displacements( + structure: Structure, + supercell_matrix: np.array, + displacement: float, + sym_reduce: bool, + symprec: float, + use_symmetrized_structure: str | None, + kpath_scheme: str, + code: str, + verbose: bool = False, +) -> 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 + use_min_dof : bool = False + Whether to use the minimal number of degrees of freedom + in calculating randomly-displaced structures. + Requires ALAMODE. + verbose : bool = False + Whether to log error messages. + + Returns + ------- + List[Structure] + Displaced structures + """ + phonon = _generate_phonon_object( + structure, + supercell_matrix, + displacement, + sym_reduce, + symprec, + use_symmetrized_structure, + kpath_scheme, + code, + verbose=verbose, + ) + supercells = phonon.supercells_with_displacements return [get_pmg_structure(cell) for cell in supercells] @job( output_schema=PhononBSDOSDoc, - data=[PhononDos, PhononBandStructureSymmLine, ForceConstants], + data=[PhononDos, PhononBandStructureSymmLine, "force_constants"], ) def generate_frequencies_eigenvectors( structure: Structure, @@ -187,8 +326,8 @@ def generate_frequencies_eigenvectors( code: str, displacement_data: dict[str, list], total_dft_energy: float, - epsilon_static: Matrix3D = None, - born: Matrix3D = None, + epsilon_static: Matrix3D | None = None, + born: Matrix3D | None = None, **kwargs, ) -> PhononBSDOSDoc: """ @@ -223,22 +362,238 @@ def generate_frequencies_eigenvectors( kwargs: dict Additional parameters that are passed to PhononBSDOSDoc.from_forces_born """ - 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, + phonon = _generate_phonon_object( + structure, + supercell_matrix, + displacement, + sym_reduce, + symprec, + use_symmetrized_structure, + kpath_scheme, + code, + verbose=False, + ) + 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, # TODO: where is this magic number coming from? + } + # Other codes could be added here + else: + borns = None + epsilon = None + + # Produces all force constants + phonon.produce_force_constants(forces=set_of_forces) + + filename_phonopy_yaml = kwargs.get("filename_phonopy_yaml", "phonopy.yaml") + create_force_constants_file = kwargs.get("create_force_constants_file", False) + force_constants_filename = kwargs.get("force_constants_filename", "FORCE_CONSTANTS") + phonon.save( + filename_phonopy_yaml, + settings={ + "force_constants": kwargs.get( + "store_force_constants", not create_force_constants_file + ) + }, + ) + if create_force_constants_file: + from phonopy.file_IO import write_FORCE_CONSTANTS + + write_FORCE_CONSTANTS( # save force_constants to text file + phonon.force_constants, filename=force_constants_filename + ) + + # get phonon band structure + kpath_dict, kpath_concrete = _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=npoints_band + ) + + # phonon band structures will always be computed + filename_band_yaml = kwargs.get("filename_band_yaml", "phonon_band_structure.yaml") + # 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 = kwargs.get("filename_dos_yaml", "phonon_dos.yaml") + # 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, + ) + + # projected dos + if kwargs.get("calculate_pdos", False): + phonon.run_mesh(kpoint.kpts[0], with_eigenvectors=True, is_mesh_symmetry=False) + phonon_dos_sigma = kwargs.get("phonon_dos_sigma") + dos_use_tetrahedron_method = kwargs.get("dos_use_tetrahedron_method", True) + phonon.run_projected_dos( + sigma=phonon_dos_sigma, + use_tetrahedron_method=dos_use_tetrahedron_method, + ) + phonon.write_projected_dos() + + phonon.run_mesh(kpoint.kpts[0]) + phonon_dos_sigma = kwargs.get("phonon_dos_sigma") + dos_use_tetrahedron_method = kwargs.get("dos_use_tetrahedron_method", True) + phonon.run_total_dos( + sigma=phonon_dos_sigma, use_tetrahedron_method=dos_use_tetrahedron_method + ) + 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"), + ) + + # 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 + ) + + volume_per_formula_unit = structure.volume / formula_units + + cls_constructor = ( + "migrate_fields" + if parse_version(_emmet_core_version) >= parse_version("0.85.1") + else "from_structure" + ) + return getattr(PhononBSDOSDoc, cls_constructor)( + structure=structure, + meta_structure=structure, + phonon_bandstructure=bs_symm_line.as_dict(), + phonon_dos=dos.as_dict(), + total_dft_energy=total_dft_energy_per_formula_unit, + volume_per_formula_unit=volume_per_formula_unit, + formula_units=formula_units, + 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, - displacement_data=displacement_data, - total_dft_energy=total_dft_energy, - epsilon_static=epsilon_static, - born=born, - **kwargs, + 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"], + }, + post_process_settings={ + "npoints_band": npoints_band, + "kpath_scheme": kpath_scheme, + "kpoint_density_dos": kpoint_density_dos, + }, + **kwargs.get("additional_fields", {}), ) @@ -247,10 +602,11 @@ def run_phonon_displacements( displacements: list[Structure], structure: Structure, supercell_matrix: Matrix3D, - phonon_maker: BaseVaspMaker | ForceFieldStaticMaker | BaseAimsMaker = None, - prev_dir: str | Path = None, - prev_dir_argname: str = None, + phonon_maker: BaseVaspMaker | ForceFieldStaticMaker | BaseAimsMaker | None = None, + prev_dir: str | Path | None = None, + prev_dir_argname: str | None = None, socket: bool = False, + store_displaced_structures: bool = False, ) -> Flow: """ Run phonon displacements. @@ -274,14 +630,15 @@ def run_phonon_displacements( argument name for the prev_dir variable socket: bool If True use the socket-io interface to increase performance + store_displaced_structures : bool = False + Whether to also save the displaced structures. """ phonon_jobs = [] - outputs: dict[str, list] = { - "displacement_number": [], - "forces": [], - "uuids": [], - "dirs": [], - } + save_props = {"displacement_number", "forces", "uuids", "dirs"} + if store_displaced_structures: + save_props.add("displaced_structures") + outputs: dict[str, list] = {k: [] for k in save_props} + phonon_job_kwargs = {} if prev_dir is not None and prev_dir_argname is not None: phonon_job_kwargs[prev_dir_argname] = prev_dir @@ -301,12 +658,13 @@ def run_phonon_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 + + # TODO: ensure order is correct. + if store_displaced_structures: + outputs["displaced_structures"] = displacements else: for idx, displacement in enumerate(displacements): - if prev_dir is not None: - phonon_job = phonon_maker.make(displacement, prev_dir=prev_dir) - else: - phonon_job = phonon_maker.make(displacement) + phonon_job = phonon_maker.make(displacement, prev_dir=prev_dir) phonon_job.append_name(f" {idx + 1}/{len(displacements)}") # we will add some meta data @@ -326,6 +684,8 @@ 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) + if store_displaced_structures: + outputs["displaced_structures"].append(displacement) displacement_flow = Flow(phonon_jobs, outputs) return Response(replace=displacement_flow) diff --git a/src/atomate2/common/jobs/qha.py b/src/atomate2/common/jobs/qha.py index 51f657cd7c..76af89564b 100644 --- a/src/atomate2/common/jobs/qha.py +++ b/src/atomate2/common/jobs/qha.py @@ -5,10 +5,11 @@ import logging from typing import TYPE_CHECKING +import numpy as np +from emmet.core.phonon import PhononBSDOSDoc from jobflow import Flow, Response, job -from atomate2.common.schemas.phonons import PhononBSDOSDoc -from atomate2.common.schemas.qha import PhononQHADoc +from atomate2.common.schemas.qha import PhononQHADoc, PhononSummaryData from atomate2.common.utils import get_supercell_matrix if TYPE_CHECKING: @@ -86,7 +87,10 @@ def get_phonon_jobs( phonon_job.append_name(f" eos deformation {istructure + 1}") phonon_jobs.append(phonon_job) outputs.append(phonon_job.output) - replace_flow = Flow(phonon_jobs, outputs) + concat_output_job = calc_thermo_data(outputs) + replace_flow = Flow( + [*phonon_jobs, concat_output_job], output=concat_output_job.output + ) return Response(replace=replace_flow) @@ -95,7 +99,7 @@ def get_phonon_jobs( data=["free_energies", "heat_capacities", "entropies", "helmholtz_volume"], ) def analyze_free_energy( - phonon_outputs: list[PhononBSDOSDoc], + phonon_outputs: list[PhononSummaryData], structure: Structure, t_max: float = None, pressure: float = None, @@ -107,8 +111,8 @@ def analyze_free_energy( Parameters ---------- - phonon_outputs: list[PhononBSDOSDoc] - list of PhononBSDOSDoc objects + phonon_outputs: list[PhononSummaryData] + list of PhononSummaryData objects structure: Structure object Corresponding structure object. t_max: float @@ -148,7 +152,7 @@ def analyze_free_energy( if (not output.has_imaginary_modes) or ignore_imaginary_modes: electronic_energies[itemp].append(output.total_dft_energy) # convert from J/mol in kJ/mol - free_energies[itemp].append(output.free_energies[itemp] / 1000.0) + free_energies[itemp].append(output.free_energies[itemp] * 1e-3) heat_capacities[itemp].append(output.heat_capacities[itemp]) entropies[itemp].append(output.entropies[itemp]) sorted_volume.append(output.volume_per_formula_unit) @@ -174,3 +178,60 @@ def analyze_free_energy( supercell_matrix=supercell_matrix, **kwargs, ) + + +@job(data=[PhononSummaryData]) +def calc_thermo_data( + phonon_docs: list[PhononBSDOSDoc], + t_min: float = 0.0, + t_max: float = 500.0, + t_step: int = 10, +) -> list[PhononSummaryData]: + """Save temperature-depdenent thermodynamic state variables. + + Parameters + ---------- + phonon_docs : list of PhononBSDOSDoc + List of phonon output documents. + t_min: float = 0.0 + Minimum temperature in K to compute data, defaults to 0 K. + t_max : float = 500. + Maximum temperature in K to compute data, defaults to 500 K. + t_step: int = 10 + Increments for temperature data in K, defaults to 10 K. + + Returns + ------- + list of PhononSummaryData containing high-level thermodynamic data + """ + temperature = np.arange(t_min, t_max, t_step) + remap = { + "temperature": "temperatures", + "entropy": "entropies", + "heat_capacity": "heat_capacities", + "internal_energy": "internal_energies", + "free_energy": "free_energies", + "structure": "meta_structure", + } + return [ + PhononSummaryData.from_structure( + **{ + remap.get(k, k): getattr(ph_doc, k, None) + for k in ( + "structure", + "total_dft_energy", + "has_imaginary_modes", + "volume_per_formula_unit", + "formula_units", + "supercell_matrix", + ) + }, + **{ + remap[k]: vals + for k, vals in ph_doc.compute_thermo_quantities( + temperature, normalization=None + ).items() + }, + ) + for ph_doc in phonon_docs + ] diff --git a/src/atomate2/common/schemas/gruneisen.py b/src/atomate2/common/schemas/gruneisen.py index f4ad6f97fd..36560af68d 100644 --- a/src/atomate2/common/schemas/gruneisen.py +++ b/src/atomate2/common/schemas/gruneisen.py @@ -23,7 +23,7 @@ from pymatgen.phonon.plotter import GruneisenPhononBSPlotter, GruneisenPlotter from typing_extensions import Self -from atomate2.common.schemas.phonons import PhononBSDOSDoc +from atomate2.common.jobs.phonons import _get_kpath logger = logging.getLogger(__name__) @@ -201,7 +201,7 @@ def from_phonon_yamls( img_format=compute_gruneisen_param_kwargs.get("img_format", "pdf"), ) # get phonon band structure - kpath_dict, kpath_concrete = PhononBSDOSDoc.get_kpath( + kpath_dict, kpath_concrete = _get_kpath( structure=structure, kpath_scheme=kpath_scheme, symprec=symprec ) qpoints, _connections = get_band_qpoints_and_path_connections( diff --git a/src/atomate2/common/schemas/phonons.py b/src/atomate2/common/schemas/phonons.py index c33e9b4b9f..7ce5222bb4 100644 --- a/src/atomate2/common/schemas/phonons.py +++ b/src/atomate2/common/schemas/phonons.py @@ -1,100 +1,48 @@ -"""Schemas for phonon documents.""" +"""Deprecated stub; use emmet.core.phonon instead.""" -import copy -import logging -from pathlib import Path -from typing import Union +from __future__ import annotations -import numpy as np -from emmet.core.math import Matrix3D -from emmet.core.structure import StructureMetadata +import warnings +from typing import TYPE_CHECKING + +from emmet.core.phonon import ( + PhononBSDOSDoc, + PhononComputationalSettings, + ThermalDisplacementData, +) from monty.json import MSONable -from phonopy import Phonopy -from phonopy.phonon.band_structure import get_band_qpoints_and_path_connections -from phonopy.structure.symmetry import symmetrize_borns_and_epsilon -from phonopy.units import VaspToTHz from pydantic import BaseModel, Field -from pymatgen.core import Structure -from pymatgen.io.phonopy import ( - get_ph_bs_symm_line, - get_ph_dos, - get_phonopy_structure, - get_pmg_structure, -) -from pymatgen.io.vasp import Kpoints -from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine -from pymatgen.phonon.dos import PhononDos -from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter -from pymatgen.symmetry.bandstructure import HighSymmKpath -from pymatgen.symmetry.kpath import KPathSeek -from typing_extensions import Self - -from atomate2.aims.utils.units import omegaToTHz - -logger = logging.getLogger(__name__) - - -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 ["ase", "forcefields", "vasp"]: - return VaspToTHz - if code == "aims": - return omegaToTHz # Based on CODATA 2002 - raise ValueError(f"Frequency conversion factor for code ({code}) not defined.") +if TYPE_CHECKING: + from emmet.core.math import Matrix3D +warnings.warn( + "atomate2.common.schemas.phonons is deprecated " + "and will be removed on 1 January, 2026. " + "Please migrate your code to use emmet.core.phonon", + stacklevel=2, + category=DeprecationWarning, +) -class PhononComputationalSettings(BaseModel): - """Collection to store computational settings for the phonon computation.""" +__all__ = [ + "ForceConstants", + "PhononBSDOSDoc", + "PhononComputationalSettings", + "PhononJobDirs", + "PhononUUIDs", + "ThermalDisplacementData", +] - # 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 ForceConstants(MSONable): + """DEPRECATED: A force constants class.""" -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: list[list[Matrix3D]] | None = Field( - None, description="field including thermal displacement matrices in CIF format" - ) - thermal_displacement_matrix: list[list[Matrix3D]] | None = Field( - None, - description="field including thermal displacement matrices in Cartesian " - "coordinate system", - ) - temperatures_thermal_displacements: list[int] | None = Field( - None, - description="temperatures at which the thermal displacement matrices" - "have been computed", - ) + def __init__(self, force_constants: list[list[Matrix3D]]) -> None: + self.force_constants = force_constants class PhononUUIDs(BaseModel): - """Collection to save all uuids connected to the phonon run.""" + """DEPRECATED: Collection to save all uuids connected to the phonon run.""" optimization_run_uuid: str | None = Field(None, description="optimization run uuid") displacements_uuids: list[str] | None = Field( @@ -104,15 +52,8 @@ class PhononUUIDs(BaseModel): born_run_uuid: str | None = 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.""" + """DEPRECATED: Save all job directories relevant for the phonon run.""" displacements_job_dirs: list[str | None] | None = Field( None, description="The directories where the displacement jobs were run." @@ -129,487 +70,3 @@ class PhononJobDirs(BaseModel): taskdoc_run_job_dir: str | None = 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: Structure | None = Field( - None, description="Structure of Materials Project." - ) - - phonon_bandstructure: PhononBandStructureSymmLine | None = Field( - None, - description="Phonon band structure object.", - ) - - phonon_dos: PhononDos | None = Field( - None, - description="Phonon density of states object.", - ) - - free_energies: list[float] | None = Field( - None, - description="vibrational part of the free energies in J/mol per " - "formula unit for temperatures in temperature_list", - ) - - heat_capacities: list[float] | None = Field( - None, - description="heat capacities in J/K/mol per " - "formula unit for temperatures in temperature_list", - ) - - internal_energies: list[float] | None = Field( - None, - description="internal energies in J/mol per " - "formula unit for temperatures in temperature_list", - ) - entropies: list[float] | None = Field( - None, - description="entropies in J/(K*mol) per formula unit" - "for temperatures in temperature_list ", - ) - - temperatures: list[int] | None = Field( - None, - description="temperatures at which the vibrational" - " part of the free energies" - " and other properties have been computed", - ) - - total_dft_energy: float | None = Field( - None, description="total DFT energy per formula unit in eV" - ) - - volume_per_formula_unit: float | None = Field( - None, description="volume per formula unit in Angstrom**3" - ) - - formula_units: int | None = Field(None, description="Formula units per cell") - - has_imaginary_modes: bool | None = Field( - None, description="if true, structure has imaginary modes" - ) - - # needed, e.g. to compute Grueneisen parameter etc - force_constants: ForceConstants | None = Field( - None, description="Force constants between every pair of atoms in the structure" - ) - - born: list[Matrix3D] | None = Field( - None, - description="Born charges as computed from phonopy. Only for symmetrically " - "different atoms", - ) - - epsilon_static: Matrix3D | None = Field( - None, description="The high-frequency dielectric constant" - ) - - supercell_matrix: Matrix3D | None = Field( - None, description="matrix describing the supercell" - ) - primitive_matrix: Matrix3D | None = Field( - None, description="matrix describing relationship to primitive cell" - ) - - code: str | None = Field( - None, description="String describing the code for the computation" - ) - - phonopy_settings: PhononComputationalSettings | None = Field( - None, description="Field including settings for Phonopy" - ) - - thermal_displacement_data: ThermalDisplacementData | None = Field( - None, - description="Includes all data of the computation of the thermal displacements", - ) - - jobdirs: PhononJobDirs | None = Field( - None, description="Field including all relevant job directories" - ) - - uuids: PhononUUIDs | None = Field( - None, description="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, - epsilon_static: Matrix3D = None, - born: Matrix3D = None, - **kwargs, - ) -> Self: - """Generate collection of phonon data. - - Parameters - ---------- - structure: Structure object - supercell_matrix: numpy array describing the supercell - displacement: float - size of displacement in angstrom - sym_reduce: bool - if True, phonopy will use symmetry - symprec: float - precision to determine kpaths, - primitive cells and symmetry in phonopy and pymatgen - use_symmetrized_structure: str - primitive, conventional or None - kpath_scheme: str - kpath scheme to generate phonon band structure - code: str - which code was used for computation - displacement_data: - output of the displacement data - total_dft_energy: float - total energy in eV per cell - epsilon_static: Matrix3D - The high-frequency dielectric constant - born: Matrix3D - Born charges - **kwargs: - additional arguments - """ - additional_fields = kwargs.get("additional_fields", {}) - factor = get_factor(code) - # This opens the opportunity to add support for other codes - # that are supported by phonopy - - cell = get_phonopy_structure(structure) - - if use_symmetrized_structure == "primitive": - primitive_matrix: Union[np.ndarray, str] = np.eye(3) - else: - primitive_matrix = "auto" - - # TARP: THIS IS BAD! Including for discussions sake - if cell.magnetic_moments is not None and primitive_matrix == "auto": - if np.any(cell.magnetic_moments != 0.0): - raise ValueError( - "For materials with magnetic moments, " - "use_symmetrized_structure must be 'primitive'" - ) - cell.magnetic_moments = None - - phonon = Phonopy( - cell, - supercell_matrix, - primitive_matrix=primitive_matrix, - factor=factor, - symprec=symprec, - is_symmetry=sym_reduce, - ) - phonon.generate_displacements(distance=displacement) - set_of_forces = [np.array(forces) for forces in displacement_data["forces"]] - - if born is not None and epsilon_static is not None: - if len(structure) == len(born): - borns, epsilon = symmetrize_borns_and_epsilon( - ucell=phonon.unitcell, - borns=np.array(born), - epsilon=np.array(epsilon_static), - symprec=symprec, - primitive_matrix=phonon.primitive_matrix, - supercell_matrix=phonon.supercell_matrix, - is_symmetry=kwargs.get("symmetrize_born", True), - ) - else: - raise ValueError( - "Number of Born charges does not agree with number of atoms" - ) - if code == "vasp" and not np.all(np.isclose(borns, 0.0)): - phonon.nac_params = { - "born": borns, - "dielectric": epsilon, - "factor": 14.399652, - } - # Other codes could be added here - else: - borns = None - epsilon = None - - # Produces all force constants - phonon.produce_force_constants(forces=set_of_forces) - - filename_phonopy_yaml = kwargs.get("filename_phonopy_yaml", "phonopy.yaml") - create_force_constants_file = kwargs.get("create_force_constants_file", False) - force_constants_filename = kwargs.get( - "force_constants_filename", "FORCE_CONSTANTS" - ) - # if kwargs.get("filename_phonopy_yaml") is None: - # kwargs["filename_phonopy_yaml"] = "phonopy.yaml" - - # with phonopy.load("phonopy.yaml") the phonopy API can be used - phonon.save( - filename_phonopy_yaml, - settings={ - "force_constants": kwargs.get( - "store_force_constants", not create_force_constants_file - ) - }, - ) - if create_force_constants_file: - from phonopy.file_IO import write_FORCE_CONSTANTS - - write_FORCE_CONSTANTS( # save force_constants to text file - phonon.force_constants, filename=force_constants_filename - ) - - # 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=npoints_band - ) - - # phonon band structures will always be computed - filename_band_yaml = kwargs.get( - "filename_band_yaml", "phonon_band_structure.yaml" - ) - # 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 = kwargs.get("filename_dos_yaml", "phonon_dos.yaml") - # 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, - ) - - # projected dos - if kwargs.get("calculate_pdos", False): - phonon.run_mesh( - kpoint.kpts[0], with_eigenvectors=True, is_mesh_symmetry=False - ) - phonon_dos_sigma = kwargs.get("phonon_dos_sigma") - dos_use_tetrahedron_method = kwargs.get("dos_use_tetrahedron_method", True) - phonon.run_projected_dos( - sigma=phonon_dos_sigma, - use_tetrahedron_method=dos_use_tetrahedron_method, - ) - phonon.write_projected_dos() - - phonon.run_mesh(kpoint.kpts[0]) - phonon_dos_sigma = kwargs.get("phonon_dos_sigma") - dos_use_tetrahedron_method = kwargs.get("dos_use_tetrahedron_method", True) - phonon.run_total_dos( - sigma=phonon_dos_sigma, use_tetrahedron_method=dos_use_tetrahedron_method - ) - phonon.write_total_dos(filename=filename_dos_yaml) - dos = get_ph_dos(filename_dos_yaml) - new_plotter_dos = PhononDosPlotter() - new_plotter_dos.add_dos(label="total", dos=dos) - new_plotter_dos.save_plot( - filename=kwargs.get("filename_dos", "phonon_dos.pdf"), - units=kwargs.get("units", "THz"), - ) - - # compute vibrational part of free energies per formula unit - temperature_range = np.arange( - kwargs.get("tmin", 0), kwargs.get("tmax", 500), kwargs.get("tstep", 10) - ) - - free_energies = [ - dos.helmholtz_free_energy( - temp=temp, structure=get_pmg_structure(phonon.primitive) - ) - for temp in temperature_range - ] - - entropies = [ - dos.entropy(temp=temp, structure=get_pmg_structure(phonon.primitive)) - for temp in temperature_range - ] - - internal_energies = [ - dos.internal_energy( - temp=temp, structure=get_pmg_structure(phonon.primitive) - ) - for temp in temperature_range - ] - - heat_capacities = [ - dos.cv(temp=temp, structure=get_pmg_structure(phonon.primitive)) - for temp in temperature_range - ] - - # will compute thermal displacement matrices - # for the primitive cell (phonon.primitive!) - # only this is available in phonopy - if kwargs.get("create_thermal_displacements"): - phonon.run_mesh( - kpoint.kpts[0], with_eigenvectors=True, is_mesh_symmetry=False - ) - freq_min_thermal_displacements = kwargs.get( - "freq_min_thermal_displacements", 0.0 - ) - phonon.run_thermal_displacement_matrices( - t_min=kwargs.get("tmin_thermal_displacements", 0), - t_max=kwargs.get("tmax_thermal_displacements", 500), - t_step=kwargs.get("tstep_thermal_displacements", 100), - freq_min=freq_min_thermal_displacements, - ) - - temperature_range_thermal_displacements = np.arange( - kwargs.get("tmin_thermal_displacements", 0), - kwargs.get("tmax_thermal_displacements", 500), - kwargs.get("tstep_thermal_displacements", 100), - ) - for idx, temp in enumerate(temperature_range_thermal_displacements): - phonon.thermal_displacement_matrices.write_cif( - phonon.primitive, idx, filename=f"tdispmat_{temp}K.cif" - ) - _disp_mat = phonon._thermal_displacement_matrices # 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 - ) - - volume_per_formula_unit = structure.volume / formula_units - - doc = 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, - volume_per_formula_unit=volume_per_formula_unit, - formula_units=formula_units, - 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, - }, - ) - - return doc.model_copy(update=additional_fields) - - @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 - """ - valid_schemes = {"setyawan_curtarolo", "latimer_munro", "hinuma", "seekpath"} - if kpath_scheme in (valid_schemes - {"seekpath"}): - 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=}, must be one of {valid_schemes}" - ) - - 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 diff --git a/src/atomate2/common/schemas/qha.py b/src/atomate2/common/schemas/qha.py index 6d6fc3e261..7df3901df0 100644 --- a/src/atomate2/common/schemas/qha.py +++ b/src/atomate2/common/schemas/qha.py @@ -14,6 +14,48 @@ logger = logging.getLogger(__name__) +class PhononSummaryData(StructureMetadata): + """Save thermodynamic state variables at a series of temperatures.""" + + structure: Structure | None = Field( + None, + description=( + "Structure associated with the phonon calculation " + "used to generate this data" + ), + ) + + supercell_matrix: Matrix3D | None = Field( + None, description="matrix describing the supercell." + ) + + total_dft_energy: float | None = Field( + None, description="The total DFT energy associated with the structure." + ) + + volume_per_formula_unit: float | None = Field( + None, description="volume per formula unit in Angstrom**3." + ) + + formula_units: int | None = Field(None, description="Formula units per cell.") + + has_imaginary_modes: bool | None = Field( + None, description="Whether the phonon spectrum has imaginary modes." + ) + + temperatures: list[float] | None = Field(None, description="Temperature in K.") + free_energies: list[float] | None = Field( + None, description="Helmholtz free energies in J/mol." + ) + heat_capacities: list[float] | None = Field( + None, description="Heat capacities in J/(K . mol)." + ) + entropies: list[float] | None = Field(None, description="Entropies in J/K.") + internal_energies: list[float] | None = Field( + None, description="Internal energies in J/mol" + ) + + class PhononQHADoc(StructureMetadata, extra="allow"): # type: ignore[call-arg] """Collection of all data produced by the qha workflow.""" diff --git a/src/atomate2/forcefields/flows/eos.py b/src/atomate2/forcefields/flows/eos.py index 4567f357ab..fab91b3d2f 100644 --- a/src/atomate2/forcefields/flows/eos.py +++ b/src/atomate2/forcefields/flows/eos.py @@ -95,6 +95,7 @@ def from_force_field_name( force_field_name=force_field_name, relax_cell=False ), static_maker=None, + **kwargs, ) return cls( name=f"{force_field_name.split('MLFF.')[-1]} EOS Maker", diff --git a/src/atomate2/forcefields/flows/pheasy.py b/src/atomate2/forcefields/flows/pheasy.py new file mode 100644 index 0000000000..d9db491809 --- /dev/null +++ b/src/atomate2/forcefields/flows/pheasy.py @@ -0,0 +1,145 @@ +"""Flows for calculating phonons.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Literal + +from atomate2 import SETTINGS +from atomate2.common.flows.pheasy import BasePhononMaker +from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker + + +@dataclass +class PhononMaker(BasePhononMaker): + """ + Maker to calculate harmonic phonons with a force field. + + Calculate the harmonic phonons of a material. Initially, a tight structural + relaxation is performed to obtain a structure without forces on the atoms. + Subsequently, supercells with one displaced atom are generated and accurate + forces are computed for these structures. With the help of phonopy, these + forces are then converted into a dynamical matrix. To correct for polarization + effects, a correction of the dynamical matrix based on BORN charges can + be performed. The BORN charges can be supplied manually. + Finally, phonon densities of states, phonon band structures + and thermodynamic properties are computed. + + .. Note:: + It is heavily recommended to symmetrize the structure before passing it to + this flow. Otherwise, a different space group might be detected and too + many displacement calculations will be generated. + It is recommended to check the convergence parameters here and + adjust them if necessary. The default might not be strict enough + for your specific case. + + Parameters + ---------- + name : str + Name of the flows produced by this maker. + sym_reduce : bool + Whether to reduce the number of deformations using symmetry. + symprec : float + Symmetry precision to use in the + reduction of symmetry to find the primitive/conventional cell + (use_primitive_standard_structure, use_conventional_standard_structure) + and to handle all symmetry-related tasks in phonopy + displacement: float + displacement distance for phonons + min_length: float + min length of the supercell that will be built + prefer_90_degrees: bool + if set to True, supercell algorithm will first try to find a supercell + with 3 90 degree angles + get_supercell_size_kwargs: dict + kwargs that will be passed to get_supercell_size to determine supercell size + use_symmetrized_structure: str + allowed strings: "primitive", "conventional", None + + - "primitive" will enforce to start the phonon computation + from the primitive standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + This makes it possible to use certain k-path definitions + with this workflow. Otherwise, we must rely on seekpath + - "conventional" will enforce to start the phonon computation + from the conventional standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + We will however use seekpath and primitive structures + as determined by from phonopy to compute the phonon band structure + bulk_relax_maker : .ForceFieldRelaxMaker or None + A maker to perform a tight relaxation on the bulk. + Set to ``None`` to skip the + bulk relaxation + static_energy_maker : .ForceFieldRelaxMaker or None + A maker to perform the computation of the DFT energy on the bulk. + Set to ``None`` to skip the + static energy computation + born_maker: .ForceFieldStaticMaker or None + Maker to compute the BORN charges. + phonon_displacement_maker : .ForceFieldStaticMaker or None + Maker used to compute the forces for a supercell. + generate_frequencies_eigenvectors_kwargs : dict + Keyword arguments passed to :obj:`generate_frequencies_eigenvectors`. + create_thermal_displacements: bool + Bool that determines if thermal_displacement_matrices are computed + kpath_scheme: str + scheme to generate kpoints. Please be aware that + you can only use seekpath with any kind of cell + Otherwise, please use the standard primitive structure + Available schemes are: + "seekpath", "hinuma", "setyawan_curtarolo", "latimer_munro". + "seekpath" and "hinuma" are the same definition but + seekpath can be used with any kind of unit cell as + it relies on phonopy to handle the relationship + to the primitive cell and not pymatgen + code: str + determines the dft or force field code. + store_force_constants: bool + if True, force constants will be stored + socket: bool + If True, use the socket for the calculation + """ + + name: str = "phonon" + sym_reduce: bool = True + symprec: float = SETTINGS.PHONON_SYMPREC + displacement: float = 0.01 + min_length: float | None = 20.0 + prefer_90_degrees: bool = True + get_supercell_size_kwargs: dict = field(default_factory=dict) + use_symmetrized_structure: Literal["primitive", "conventional"] | None = None + bulk_relax_maker: ForceFieldRelaxMaker | None = field( + default_factory=lambda: ForceFieldRelaxMaker( + force_field_name="MACE", relax_kwargs={"fmax": 0.00001} + ) + ) + static_energy_maker: ForceFieldStaticMaker | None = field( + default_factory=lambda: ForceFieldStaticMaker(force_field_name="MACE") + ) + phonon_displacement_maker: ForceFieldStaticMaker = field( + default_factory=lambda: ForceFieldStaticMaker(force_field_name="MACE") + ) + create_thermal_displacements: bool = False + generate_frequencies_eigenvectors_kwargs: dict = field(default_factory=dict) + kpath_scheme: str = "seekpath" + store_force_constants: bool = True + code: str = "forcefields" + born_maker: ForceFieldStaticMaker | None = None + + @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/openff/__init__.py b/src/atomate2/openff/__init__.py index 762a41ab61..a000ae41e2 100644 --- a/src/atomate2/openff/__init__.py +++ b/src/atomate2/openff/__init__.py @@ -1,5 +1,13 @@ """Module for classical md workflows.""" +from importlib.util import find_spec + +if not find_spec("openff"): + raise ImportError( + "openff must be installed via conda forge to use with atomate2:\n" + "conda install -c conda-force openff" + ) + from openff.interchange import Interchange from openff.toolkit.topology import Topology from openff.toolkit.topology.molecule import Molecule diff --git a/src/atomate2/vasp/flows/pheasy.py b/src/atomate2/vasp/flows/pheasy.py new file mode 100644 index 0000000000..fca5039a00 --- /dev/null +++ b/src/atomate2/vasp/flows/pheasy.py @@ -0,0 +1,195 @@ +"""Define the VASP PhononMaker.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Literal + +from atomate2 import SETTINGS +from atomate2.common.flows.pheasy import BasePhononMaker +from atomate2.vasp.flows.core import DoubleRelaxMaker +from atomate2.vasp.jobs.core import DielectricMaker, StaticMaker, TightRelaxMaker +from atomate2.vasp.jobs.phonons import PhononDisplacementMaker +from atomate2.vasp.sets.core import StaticSetGenerator + +if TYPE_CHECKING: + from atomate2.vasp.jobs.base import BaseVaspMaker + + +@dataclass +class PhononMaker(BasePhononMaker): + """Maker to calculate harmonic phonons with LASSO-based ML code Pheasy. + + Calculate the zero-K harmonic phonons of a material and higher-order FCs. + Initially, a tight structural relaxation is performed to obtain a structure + without forces on the atoms. Subsequently, supercells with all atoms displaced + by a small amplitude (generally using 0.01 A) are generated and accurate forces + are computed for these structures for the second order force constants. With the + help of pheasy (LASSO technique), these forces are then converted into a dynamical + matrix. In this Workflow, we separate the harmonic phonon calculations and + anharmonic force constants calculations. To correct for polarization effects, a + correction of the dynamical matrix based on BORN charges can be performed. Finally, + phonon densities of states, phonon band structures and thermodynamic properties + are computed. For the anharmonic force constants, the supercells with all atoms + displaced by a larger amplitude (generally using 0.08 A) are generated and accurate + forces are computed for these structures. With the help of pheasy (LASSO technique), + the third- and fourth-order force constants are extracted at once. + + .. Note:: + It is heavily recommended to symmetrize the structure before passing it to + this flow. Otherwise, a different space group might be detected and too + many displacement calculations will be generated. + It is recommended to check the convergence parameters here and + adjust them if necessary. The default might not be strict enough + for your specific case. + + Parameters + ---------- + name : str + Name of the flow produced by this maker. + sym_reduce : bool + Whether to reduce the number of deformations using symmetry. + symprec : float + Symmetry precision to use in the + reduction of symmetry to find the primitive/conventional cell + (use_primitive_standard_structure, use_conventional_standard_structure) + and to handle all symmetry-related tasks in pheasy, we recommend to + use the value of 1e-3. + displacement: float + displacement distance for phonons, for most cases 0.01 A is a good choice, + but it can be increased to 0.02 A for heavier elements. + num_displaced_supercells: int + number of displacements to be generated using a random-displacement approach + for harmonic phonon calculations. The default value is 0 and the number of + displacements is automatically determined by the number of atoms in the + supercell and its space group. + cal_anhar_fcs: bool + if set to True, anharmonic force constants(FCs) up to fourth-order FCs will + be calculated. The default value is False, and only harmonic phonons will + be calculated. + displacement_anhar: float + displacement distance for anharmonic force constants(FCs) up to fourth-order + FCs, for most cases 0.08 A is a good choice, but it can be increased to 0.1 A. + num_disp_anhar: int + number of displacements to be generated using a random-displacement approach + for anharmonic phonon calculations. The default value is 0 and the number of + displacements is automatically determined by the number of atoms in the + supercell, cutoff distance for anharmonic FCs its space group. generally, + 50 large-distance displacements are enough for most cases. + fcs_cutoff_radius: list + cutoff distance for anharmonic force constants(FCs) up to fourth-order FCs. + The default value is [-1, 12, 10], which means that the cutoff distance for + second-order FCs is the Wigner-Seitz cell boundary and the cutoff distance + for third-order FCs is 12 Borh, and the cutoff distance for fourth-order FCs + is 10 Bohr. Generally, the default value is good enough. + min_length: float + minimum length of lattice constants will be used to create the supercell, + the default value is 14.0 A. In most cases, the default value is good + enough, but it can be increased for larger supercells. + prefer_90_degrees: bool + if set to True, supercell algorithm will first try to find a supercell + with 3 90 degree angles. + get_supercell_size_kwargs: dict + kwargs that will be passed to get_supercell_size to determine supercell size + use_symmetrized_structure: str + allowed strings: "primitive", "conventional", None + + - "primitive" will enforce to start the phonon computation + from the primitive standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + This makes it possible to use certain k-path definitions + with this workflow. Otherwise, we must rely on seekpath + - "conventional" will enforce to start the phonon computation + from the conventional standard structure + according to Setyawan, W., & Curtarolo, S. (2010). + High-throughput electronic band structure calculations: + Challenges and tools. Computational Materials Science, + 49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010. + We will however use seekpath and primitive structures + as determined by from phonopy to compute the phonon band structure + bulk_relax_maker: .BaseVaspMaker, or None + A maker to perform a tight relaxation on the bulk. + Set to ``None`` to skip the + bulk relaxation + static_energy_maker: .BaseVaspMaker, or None + A maker to perform the computation of the DFT energy on the bulk. + Set to ``None`` to skip the + static energy computation + born_maker: .BaseVaspMaker, or None + Maker to compute the BORN charges. + phonon_displacement_maker: .BaseVaspMaker + Maker used to compute the forces for a supercell. + generate_frequencies_eigenvectors_kwargs : dict + Keyword arguments passed to :obj:`generate_frequencies_eigenvectors`. + create_thermal_displacements: bool + Bool that determines if thermal_displacement_matrices are computed + kpath_scheme: str + scheme to generate kpoints. Please be aware that + you can only use seekpath with any kind of cell + Otherwise, please use the standard primitive structure + Available schemes are: + "seekpath", "hinuma", "setyawan_curtarolo", "latimer_munro". + "seekpath" and "hinuma" are the same definition but + seekpath can be used with any kind of unit cell as + it relies on phonopy to handle the relationship + to the primitive cell and not pymatgen + code: str = "vasp" + determines the DFT code. currently only vasp is implemented. + This keyword might enable the implementation of other codes + in the future + store_force_constants: bool + if True, force constants will be stored + socket: bool + If True, use the socket for the calculation + """ + + name: str = "phonon" + sym_reduce: bool = True + symprec: float = SETTINGS.PHONON_SYMPREC + cal_anhar_fcs: bool = False + displacement: float = 0.01 + displacement_anhar: float = 0.08 + num_displaced_supercells: int = 0 + num_disp_anhar: int = 0 + fcs_cutoff_radius: list = field(default_factory=lambda: [-1, 12, 10]) + min_length: float | None = 8.0 + max_atoms: float | None = 200 + force_90_degrees: bool = True + force_diagonal: bool = True + allow_orthorhombic: bool = False + prefer_90_degrees: bool = True + get_supercell_size_kwargs: dict = field(default_factory=dict) + use_symmetrized_structure: Literal["primitive", "conventional"] | None = None + create_thermal_displacements: bool = False + generate_frequencies_eigenvectors_kwargs: dict = field(default_factory=dict) + kpath_scheme: str = "seekpath" + store_force_constants: bool = True + socket: bool = False + code: str = "vasp" + bulk_relax_maker: BaseVaspMaker | None = field( + default_factory=lambda: DoubleRelaxMaker.from_relax_maker(TightRelaxMaker()) + ) + static_energy_maker: BaseVaspMaker | None = field( + default_factory=lambda: StaticMaker( + input_set_generator=StaticSetGenerator(auto_ispin=True) + ) + ) + born_maker: BaseVaspMaker | None = field(default_factory=DielectricMaker) + phonon_displacement_maker: BaseVaspMaker = field( + default_factory=PhononDisplacementMaker + ) + + @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/tests/aims/test_flows/test_phonon_workflow.py b/tests/aims/test_flows/test_phonon_workflow.py index f6477a9702..f90ce3e83c 100644 --- a/tests/aims/test_flows/test_phonon_workflow.py +++ b/tests/aims/test_flows/test_phonon_workflow.py @@ -1,6 +1,5 @@ """Test various makers""" -import json from pathlib import Path import pytest @@ -11,6 +10,34 @@ si_structure_file = Path(__file__).parents[2] / "test_data/structures/Si_diamond.cif" +phonopy_settings_schema = { + "description": "Collection to store computational settings for " + "the phonon computation.", + "properties": { + "npoints_band": { + "description": "number of points for band structure computation", + "title": "Npoints Band", + "anyOf": [{"type": "integer"}, {"type": "null"}], + "default": None, + }, + "kpath_scheme": { + "anyOf": [{"type": "string"}, {"type": "null"}], + "default": None, + "description": "indicates the kpath scheme", + "title": "Kpath Scheme", + }, + "kpoint_density_dos": { + "anyOf": [{"type": "integer"}, {"type": "null"}], + "default": None, + "description": "number of points for computation of free energies " + "and densities of states", + "title": "Kpoint Density Dos", + }, + }, + "title": "PhononComputationalSettings", + "type": "object", +} + def test_phonon_flow(clean_dir, mock_aims, species_dir): import numpy as np @@ -67,39 +94,14 @@ def test_phonon_flow(clean_dir, mock_aims, species_dir): # validation the outputs of the job output = responses[flow.job_uuids[-1]][1].output - phonopy_settings_schema = { - "description": "Collection to store computational settings for the " - "phonon computation.", - "properties": { - "npoints_band": { - "default": "number of points for band structure computation", - "title": "Npoints Band", - "type": "integer", - }, - "kpath_scheme": { - "default": "indicates the kpath scheme", - "title": "Kpath Scheme", - "type": "string", - }, - "kpoint_density_dos": { - "default": "number of points for computation of free energies and" - " densities of states", - "title": "Kpoint Density Dos", - "type": "integer", - }, - }, - "title": "PhononComputationalSettings", - "type": "object", - } assert output.code == "aims" assert output.born is None assert not output.has_imaginary_modes - assert output.temperatures == list(range(0, 500, 100)) - assert output.heat_capacities[0] == 0.0 - assert np.round(output.heat_capacities[-1], 2) == 21.95 - assert output.phonopy_settings.schema_json() == json.dumps(phonopy_settings_schema) - assert np.round(output.phonon_bandstructure.bands[-1, 0], 2) == 15.1 + assert output.heat_capacity(0.0) == 0.0 + assert output.heat_capacity(400.0) == pytest.approx(21.95, abs=1e-2) + assert output.post_process_settings.schema() == phonopy_settings_schema + assert np.round(output.phonon_bandstructure.frequencies[-1][0], 2) == 15.1 @pytest.mark.skip(reason="Currently not mocked and needs FHI-aims binary") @@ -157,39 +159,14 @@ def test_phonon_socket_flow(si, clean_dir, mock_aims, species_dir): # validation the outputs of the job output = responses[flow.job_uuids[-1]][1].output - phonopy_settings_schema = { - "description": "Collection to store computational settings for the " - "phonon computation.", - "properties": { - "npoints_band": { - "default": "number of points for band structure computation", - "title": "Npoints Band", - "type": "integer", - }, - "kpath_scheme": { - "default": "indicates the kpath scheme", - "title": "Kpath Scheme", - "type": "string", - }, - "kpoint_density_dos": { - "default": "number of points for computation of free energies and" - " densities of states", - "title": "Kpoint Density Dos", - "type": "integer", - }, - }, - "title": "PhononComputationalSettings", - "type": "object", - } assert output.code == "aims" assert output.born is None assert not output.has_imaginary_modes - assert output.temperatures == list(range(0, 500, 10)) - assert output.heat_capacities[0] == 0.0 - assert np.round(output.heat_capacities[-1], 2) == 23.06 - assert output.phonopy_settings.schema_json() == json.dumps(phonopy_settings_schema) - assert np.round(output.phonon_bandstructure.bands[-1, 0], 2) == 14.41 + assert output.heat_capacity(0.0) == 0.0 + assert output.heat_capacity(400.0) == pytest.approx(23.06, abs=1e-2) + assert output.post_process_settings.schema() == phonopy_settings_schema + assert np.round(output.phonon_bandstructure.frequencies[-1][0], 2) == 14.41 def test_phonon_default_flow(si, clean_dir, mock_aims, species_dir): @@ -225,39 +202,15 @@ def test_phonon_default_flow(si, clean_dir, mock_aims, species_dir): # validation the outputs of the job output = responses[flow.job_uuids[-1]][1].output - phonopy_settings_schema = { - "description": "Collection to store computational settings for " - "the phonon computation.", - "properties": { - "npoints_band": { - "default": "number of points for band structure computation", - "title": "Npoints Band", - "type": "integer", - }, - "kpath_scheme": { - "default": "indicates the kpath scheme", - "title": "Kpath Scheme", - "type": "string", - }, - "kpoint_density_dos": { - "default": "number of points for computation of free energies " - "and densities of states", - "title": "Kpoint Density Dos", - "type": "integer", - }, - }, - "title": "PhononComputationalSettings", - "type": "object", - } + assert output.code == "aims" assert output.born is None assert not output.has_imaginary_modes - assert output.temperatures == list(range(0, 500, 10)) - assert output.heat_capacities[0] == 0.0 - assert np.round(output.heat_capacities[-1], 2) == 22.85 - assert output.phonopy_settings.schema_json() == json.dumps(phonopy_settings_schema) - assert np.round(output.phonon_bandstructure.bands[-1, 0], 2) == 15.02 + assert output.heat_capacity(0.0) == 0.0 + assert output.heat_capacity(490.0) == pytest.approx(22.85, abs=1e-2) + assert output.post_process_settings.schema() == phonopy_settings_schema + assert np.round(output.phonon_bandstructure.frequencies[-1][0], 2) == 15.02 if aims_sd is not None: SETTINGS["AIMS_SPECIES_DIR"] = aims_sd @@ -299,39 +252,14 @@ def test_phonon_default_socket_flow(si, clean_dir, mock_aims, species_dir): # validation the outputs of the job output = responses[flow.job_uuids[-1]][1].output - phonopy_settings_schema = { - "description": "Collection to store computational settings for " - "the phonon computation.", - "properties": { - "npoints_band": { - "default": "number of points for band structure computation", - "title": "Npoints Band", - "type": "integer", - }, - "kpath_scheme": { - "default": "indicates the kpath scheme", - "title": "Kpath Scheme", - "type": "string", - }, - "kpoint_density_dos": { - "default": "number of points for computation of free energies " - "and densities of states", - "title": "Kpoint Density Dos", - "type": "integer", - }, - }, - "title": "PhononComputationalSettings", - "type": "object", - } assert output.code == "aims" assert output.born is None assert not output.has_imaginary_modes - assert output.temperatures == list(range(0, 500, 10)) - assert output.heat_capacities[0] == 0.0 - assert np.round(output.heat_capacities[-1], 2) == 22.85 - assert output.phonopy_settings.schema_json() == json.dumps(phonopy_settings_schema) - assert np.round(output.phonon_bandstructure.bands[-1, 0], 2) == 15.02 + assert output.heat_capacity(0.0) == 0.0 + assert output.heat_capacity(490.0) == pytest.approx(22.85, abs=1e-2) + assert output.post_process_settings.schema() == phonopy_settings_schema + assert np.round(output.phonon_bandstructure.frequencies[-1][0], 2) == 15.02 if aims_sd is not None: SETTINGS["AIMS_SPECIES_DIR"] = aims_sd diff --git a/tests/common/schemas/test_phonons.py b/tests/common/schemas/test_phonons.py deleted file mode 100644 index f0ddd9315f..0000000000 --- a/tests/common/schemas/test_phonons.py +++ /dev/null @@ -1,60 +0,0 @@ -import json - -import numpy as np -import pytest -from monty.json import MontyEncoder -from pydantic import ValidationError - -from atomate2.common.schemas.phonons import ( - PhononBSDOSDoc, - PhononComputationalSettings, - PhononJobDirs, - PhononUUIDs, - ThermalDisplacementData, -) - - -def test_thermal_displacement_data(): - doc = ThermalDisplacementData(freq_min_thermal_displacements=0.0) - validated = ThermalDisplacementData.model_validate_json( - json.dumps(doc, cls=MontyEncoder) - ) - assert isinstance(validated, ThermalDisplacementData) - - -def test_phonon_bs_dos_doc(): - kwargs = { - "total_dft_energy": None, - "supercell_matrix": np.eye(3), - "primitive_matrix": np.eye(3), - "code": "test", - "phonopy_settings": PhononComputationalSettings( - npoints_band=1, kpath_scheme="test", kpoint_density_dos=1 - ), - "thermal_displacement_data": None, - "jobdirs": None, - "uuids": None, - } - doc = PhononBSDOSDoc(**kwargs) - # check validation raises no errors - validated = PhononBSDOSDoc.model_validate_json(json.dumps(doc, cls=MontyEncoder)) - assert isinstance(validated, PhononBSDOSDoc) - - # test invalid supercell_matrix type fails - with pytest.raises(ValidationError): - doc = PhononBSDOSDoc(**kwargs | {"supercell_matrix": (1, 1, 1)}) - - # test optional material_id - doc = PhononBSDOSDoc(**kwargs | {"material_id": 1234}) - assert doc.material_id == 1234 - - # test extra="allow" option - doc = PhononBSDOSDoc(**kwargs | {"extra_field": "test"}) - assert doc.extra_field == "test" - - -# schemas where all fields have default values -@pytest.mark.parametrize("model_cls", [PhononJobDirs, PhononUUIDs]) -def test_model_validate(model_cls): - validated = model_cls.model_validate_json(json.dumps(model_cls(), cls=MontyEncoder)) - assert isinstance(validated, model_cls) diff --git a/tests/common/schemas/test_qha.py b/tests/common/schemas/test_qha.py index a3da77eecd..9a6643abf8 100644 --- a/tests/common/schemas/test_qha.py +++ b/tests/common/schemas/test_qha.py @@ -3,8 +3,8 @@ from pymatgen.core.structure import Structure from ruamel.yaml import YAML -from atomate2.common.jobs.phonons import PhononBSDOSDoc from atomate2.common.jobs.qha import PhononQHADoc, analyze_free_energy +from atomate2.common.schemas.qha import PhononSummaryData def test_analyze_free_energy(tmp_dir, test_dir): @@ -66,7 +66,7 @@ def test_analyze_free_energy(tmp_dir, test_dir): fe = [v["free_energy"] * 1000.0 for v in thermal_properties] phonon_docs.append( - PhononBSDOSDoc( + PhononSummaryData( free_energies=fe, heat_capacities=cv, entropies=entropy, @@ -143,7 +143,7 @@ def test_analyze_free_energy_small(tmp_dir, test_dir): fe = [v["free_energy"] * 1000.0 for v in thermal_properties] phonon_docs.append( - PhononBSDOSDoc( + PhononSummaryData( free_energies=fe, heat_capacities=cv, entropies=entropy, diff --git a/tests/forcefields/flows/test_phonon.py b/tests/forcefields/flows/test_phonon.py index d92b8a14a2..c7476786be 100644 --- a/tests/forcefields/flows/test_phonon.py +++ b/tests/forcefields/flows/test_phonon.py @@ -2,18 +2,17 @@ from pathlib import Path import pytest +from emmet.core.phonon import ( + CalcMeta, + PhononBS, + PhononBSDOSDoc, + PhononComputationalSettings, + PhononDOS, +) 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.schemas.phonons import ( - PhononBSDOSDoc, - PhononComputationalSettings, - PhononJobDirs, - PhononUUIDs, -) from atomate2.forcefields.flows.phonons import PhononMaker from atomate2.forcefields.jobs import ForceFieldRelaxMaker @@ -63,23 +62,15 @@ def test_phonon_wf_force_field( ph_bs_dos_doc = responses[flow[-1].uuid][1].output assert isinstance(ph_bs_dos_doc, PhononBSDOSDoc) - assert_allclose( - ph_bs_dos_doc.free_energies, - [5058.4521752, 4907.4957516, 3966.5493299, 2157.8178928, -357.5054580], - atol=1000, - ) - ph_band_struct = ph_bs_dos_doc.phonon_bandstructure - assert isinstance(ph_band_struct, PhononBandStructureSymmLine) + assert isinstance(ph_band_struct, PhononBS) ph_dos = ph_bs_dos_doc.phonon_dos - assert isinstance(ph_dos, PhononDos) + assert isinstance(ph_dos, PhononDOS) assert ph_bs_dos_doc.thermal_displacement_data is None assert isinstance(ph_bs_dos_doc.structure, Structure) - assert_allclose(ph_bs_dos_doc.temperatures, [0, 100, 200, 300, 400]) assert ph_bs_dos_doc.force_constants is None - assert isinstance(ph_bs_dos_doc.jobdirs, PhononJobDirs) - assert isinstance(ph_bs_dos_doc.uuids, PhononUUIDs) + assert all(isinstance(cm, CalcMeta) for cm in ph_bs_dos_doc.calc_meta) assert_allclose(ph_bs_dos_doc.total_dft_energy, -5.37245798, 4) assert ph_bs_dos_doc.born is None assert ph_bs_dos_doc.epsilon_static is None @@ -93,24 +84,37 @@ def test_phonon_wf_force_field( atol=1e-8, ) assert ph_bs_dos_doc.code == "forcefields" - assert isinstance(ph_bs_dos_doc.phonopy_settings, PhononComputationalSettings) - assert ph_bs_dos_doc.phonopy_settings.npoints_band == 101 - assert ph_bs_dos_doc.phonopy_settings.kpath_scheme == "seekpath" - assert ph_bs_dos_doc.phonopy_settings.kpoint_density_dos == 7_000 - assert_allclose( - ph_bs_dos_doc.entropies, - [0.0, 4.78393981, 13.99318695, 21.88641334, 28.19110667], - atol=2, - ) - assert_allclose( - ph_bs_dos_doc.heat_capacities, - [0.0, 8.86060586, 17.55758943, 21.08903916, 22.62587271], - atol=2, + assert isinstance(ph_bs_dos_doc.post_process_settings, PhononComputationalSettings) + assert ph_bs_dos_doc.post_process_settings.npoints_band == 101 + assert ph_bs_dos_doc.post_process_settings.kpath_scheme == "seekpath" + assert ph_bs_dos_doc.post_process_settings.kpoint_density_dos == 7_000 + + ref_vals = { + "entropy": [0.0, 7.45806197, 24.99582177, 40.53981354, 53.0450785], + "heat_capacity": [0.0, 15.9212379, 34.32542093, 41.73809612, 44.95600976], + "internal_energy": [ + 10510.17946131, + 11038.76862405, + 13676.21828021, + 17534.72238986, + 21889.29538244, + ], + "free_energy": [ + 10510.17946131, + 10292.96242722, + 8677.0539271, + 5372.77832663, + 671.26398379, + ], + } + thermo_props = ph_bs_dos_doc.compute_thermo_quantities( + [0, 100, 200, 300, 400], normalization=None ) - assert_allclose( - ph_bs_dos_doc.internal_energies, - [5058.44158791, 5385.88058579, 6765.19854165, 8723.78588089, 10919.0199409], - atol=1000, + + assert all( + thermo_props[k][i] == pytest.approx(val, rel=0.1) + for k, vals in ref_vals.items() + for i, val in enumerate(vals) ) # check phonon plots exist diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/inputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/INCAR.gz new file mode 100644 index 0000000000..df8dbd7760 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..9f1c1ad5e8 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POSCAR.gz new file mode 100644 index 0000000000..a393c90ee2 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..9b2f5f8175 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/CONTCAR.gz new file mode 100644 index 0000000000..c84d13a711 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/outputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/INCAR.gz new file mode 100644 index 0000000000..df8dbd7760 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..9f1c1ad5e8 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/OUTCAR.gz new file mode 100644 index 0000000000..8037e1d4fd Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POSCAR.gz new file mode 100644 index 0000000000..a393c90ee2 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..9b2f5f8175 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/dielectric/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..652ee17466 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/dielectric/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/INCAR.gz new file mode 100644 index 0000000000..933bcccb96 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..98ee85e2d8 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/KPOINTS.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/KPOINTS.gz new file mode 100644 index 0000000000..6024584b1d Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..2ef5ae4de7 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POSCAR.gz new file mode 100644 index 0000000000..bd6627e23d Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..0d45f1dab8 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/CONTCAR.gz new file mode 100644 index 0000000000..3c05a8c698 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/INCAR.gz new file mode 100644 index 0000000000..933bcccb96 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..98ee85e2d8 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/KPOINTS.gz new file mode 100644 index 0000000000..6024584b1d Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..2ef5ae4de7 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/OUTCAR.gz new file mode 100644 index 0000000000..e030514420 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POSCAR.gz new file mode 100644 index 0000000000..bd6627e23d Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..0d45f1dab8 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..2ce85dbb84 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_1_2/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/INCAR.gz new file mode 100644 index 0000000000..16e6f2ca36 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..24b6a95d02 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/KPOINTS.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/KPOINTS.gz new file mode 100644 index 0000000000..d6bf6427a4 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..dc147b5962 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POSCAR.gz new file mode 100644 index 0000000000..053d239170 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..4923c8cf94 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/CONTCAR.gz new file mode 100644 index 0000000000..87d01fccb9 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/INCAR.gz new file mode 100644 index 0000000000..16e6f2ca36 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..24b6a95d02 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/KPOINTS.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/KPOINTS.gz new file mode 100644 index 0000000000..d6bf6427a4 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/KPOINTS.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/KPOINTS.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/KPOINTS.orig.gz new file mode 100644 index 0000000000..dc147b5962 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/KPOINTS.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/OUTCAR.gz new file mode 100644 index 0000000000..07b205a9df Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POSCAR.gz new file mode 100644 index 0000000000..053d239170 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..4923c8cf94 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..f44a08d39b Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/phonon_static_2_2/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/inputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/static/inputs/INCAR.gz new file mode 100644 index 0000000000..de79338e0e Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/static/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..18a00cc740 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/inputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/static/inputs/POSCAR.gz new file mode 100644 index 0000000000..a5aea21035 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/static/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..8209f6e0aa Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/static/inputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_pheasy/static/outputs/CONTCAR.gz new file mode 100644 index 0000000000..23797aae36 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/outputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/static/outputs/INCAR.gz new file mode 100644 index 0000000000..de79338e0e Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/static/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..18a00cc740 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_pheasy/static/outputs/OUTCAR.gz new file mode 100644 index 0000000000..c6e4600a7c Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/outputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/static/outputs/POSCAR.gz new file mode 100644 index 0000000000..a5aea21035 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/static/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..8209f6e0aa Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/static/outputs/POTCAR.spec.gz new file mode 100644 index 0000000000..03d35c79c0 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/static/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_pheasy/static/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..c4ddad2006 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/static/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/INCAR.gz new file mode 100644 index 0000000000..2699cdc1e4 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..504bf3725d Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/POSCAR.gz new file mode 100644 index 0000000000..a0d31abdfe Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..cf967cea39 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/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_pheasy/tight_relax_1/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/CONTCAR.gz new file mode 100644 index 0000000000..f9b41a941d Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/INCAR.gz new file mode 100644 index 0000000000..2699cdc1e4 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..504bf3725d Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/OUTCAR.gz new file mode 100644 index 0000000000..47230006c9 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/POSCAR.gz new file mode 100644 index 0000000000..a0d31abdfe Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..cf967cea39 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/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_pheasy/tight_relax_1/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..4c181d8ee7 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_1/outputs/vasprun.xml.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/INCAR.gz new file mode 100644 index 0000000000..48def910ba Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/INCAR.orig.gz new file mode 100644 index 0000000000..ac7a6dbe93 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/POSCAR.gz new file mode 100644 index 0000000000..8b9450aac9 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/POSCAR.orig.gz new file mode 100644 index 0000000000..0a4823de33 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/inputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/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_pheasy/tight_relax_2/inputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/CONTCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/CONTCAR.gz new file mode 100644 index 0000000000..fc60000106 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/CONTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/INCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/INCAR.gz new file mode 100644 index 0000000000..48def910ba Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/INCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/INCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/INCAR.orig.gz new file mode 100644 index 0000000000..ac7a6dbe93 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/INCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/OUTCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/OUTCAR.gz new file mode 100644 index 0000000000..6877560f70 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/OUTCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/POSCAR.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/POSCAR.gz new file mode 100644 index 0000000000..8b9450aac9 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/POSCAR.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/POSCAR.orig.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/POSCAR.orig.gz new file mode 100644 index 0000000000..0a4823de33 Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/POSCAR.orig.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/POTCAR.spec.gz b/tests/test_data/vasp/Si_pheasy/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_pheasy/tight_relax_2/outputs/POTCAR.spec.gz differ diff --git a/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/vasprun.xml.gz b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/vasprun.xml.gz new file mode 100644 index 0000000000..b7a814eaac Binary files /dev/null and b/tests/test_data/vasp/Si_pheasy/tight_relax_2/outputs/vasprun.xml.gz differ diff --git a/tests/vasp/flows/test_pheasy.py b/tests/vasp/flows/test_pheasy.py new file mode 100644 index 0000000000..4a323a4dcf --- /dev/null +++ b/tests/vasp/flows/test_pheasy.py @@ -0,0 +1,129 @@ +from emmet.core.phonon import ( + CalcMeta, + PhononBS, + PhononBSDOSDoc, + PhononComputationalSettings, + PhononDOS, + ThermalDisplacementData, +) +from jobflow import run_locally +from numpy.testing import assert_allclose +from pymatgen.core.structure import Structure + +from atomate2.common.flows.pheasy import BasePhononMaker +from atomate2.common.powerups import add_metadata_to_flow +from atomate2.vasp.flows.pheasy import PhononMaker +from atomate2.vasp.jobs.base import BaseVaspMaker +from atomate2.vasp.powerups import update_user_incar_settings + + +def test_pheasy_wf_vasp(mock_vasp, clean_dir, si_structure: Structure, test_dir): + # mapping from job name to directory containing test files + ref_paths = { + "tight relax 1": "Si_pheasy/tight_relax_1", + "tight relax 2": "Si_pheasy/tight_relax_2", + "phonon static 1/2": "Si_pheasy/phonon_static_1_2", + "phonon static 2/2": "Si_pheasy/phonon_static_2_2", + "static": "Si_pheasy/static", + "dielectric": "Si_pheasy/dielectric", + } + + # 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", "KSPACING"]}, + "tight relax 2": {"incar_settings": ["NSW", "ISMEAR", "KSPACING"]}, + "phonon static 1/2": {"incar_settings": ["NSW", "ISMEAR"]}, + "phonon static 2/2": {"incar_settings": ["NSW", "ISMEAR"]}, + "static": {"incar_settings": ["NSW", "ISMEAR"]}, + "dielectric": {"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_pheasy/tight_relax_1/inputs/POSCAR.gz" + ) + + job = PhononMaker( + force_diagonal=True, + min_length=12, + cal_anhar_fcs=False, + # use_symmetrized_structure="primitive" + ).make(structure=si_struct) + + job = update_user_incar_settings( + job, + { + "ENCUT": 600, + "ISMEAR": 0, + "SIGMA": 0.05, + "KSPACING": 0.15, + "ISPIN": 1, + "EDIFFG": -1e-04, + "EDIFF": 1e-07, + }, + ) + job = add_metadata_to_flow( + flow=job, + additional_fields={"mp_id": "mp-149", "unit_testing": "yes"}, + class_filter=(BaseVaspMaker, BasePhononMaker, PhononMaker), + ) + + # run the flow or job and ensure that it finished running successfully + responses = run_locally(job, create_folders=True, ensure_success=True) + ph_doc = responses[job.jobs[-1].uuid][1].output + + # validate the outputs + assert isinstance(ph_doc, PhononBSDOSDoc) + + assert isinstance( + ph_doc.phonon_bandstructure, + PhononBS, + ) + assert isinstance(ph_doc.phonon_dos, PhononDOS) + assert isinstance( + ph_doc.thermal_displacement_data, + ThermalDisplacementData, + ) + assert isinstance(ph_doc.structure, Structure) + assert ph_doc.has_imaginary_modes is False + assert isinstance(ph_doc.force_constants, tuple) + assert all(isinstance(cm, CalcMeta) for cm in ph_doc.calc_meta) + assert_allclose(ph_doc.total_dft_energy, -5.7466748) + assert_allclose( + ph_doc.born, + [ + ((0.0, 0.0, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0)), + ((0.0, 0.0, 0.0), (0.0, 0.0, 0.0), (0.0, 0.0, 0.0)), + ], + ) + assert_allclose( + ph_doc.epsilon_static, + ( + (13.31020238, 0.0, -0.0), + (0.0, 13.31020238, 0.0), + (0.0, -0.0, 13.31020238), + ), + atol=1e-8, + ) + assert_allclose( + ph_doc.supercell_matrix, + [[4.0, 0.0, 0.0], [0.0, 4.0, 0.0], [0.0, 0.0, 4.0]], + ) + assert_allclose( + ph_doc.primitive_matrix, + [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], + rtol=1e-5, + atol=1e-10, + ) + assert ph_doc.code == "vasp" + assert isinstance( + ph_doc.post_process_settings, + PhononComputationalSettings, + ) + assert ph_doc.post_process_settings.npoints_band == 101 + assert ph_doc.post_process_settings.kpath_scheme == "seekpath" + assert ph_doc.post_process_settings.kpoint_density_dos == 7000 + + assert ph_doc.chemsys == "Si" diff --git a/tests/vasp/flows/test_phonons.py b/tests/vasp/flows/test_phonons.py index 057a803988..ad5c7d78ad 100644 --- a/tests/vasp/flows/test_phonons.py +++ b/tests/vasp/flows/test_phonons.py @@ -1,20 +1,21 @@ import numpy as np import pytest -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.schemas.phonons import ( +from emmet.core.base import CalcMeta +from emmet.core.phonon import ( + PhononBS, PhononBSDOSDoc, PhononComputationalSettings, - PhononJobDirs, - PhononUUIDs, + PhononDOS, ThermalDisplacementData, ) +from jobflow import run_locally +from numpy.testing import assert_allclose +from pymatgen.core.structure import Structure + from atomate2.vasp.flows.phonons import PhononMaker +# TODO: re-enable `total_dft_energy` checks once patched on emmet-core side + def test_phonon_wf_vasp_only_displacements3( mock_vasp, clean_dir, si_structure: Structure @@ -49,72 +50,67 @@ def test_phonon_wf_vasp_only_displacements3( 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, - [6115.980051, 6059.749756, 5490.929122, 4173.234384, 2194.164562], - ) + ph_doc = responses[job.jobs[-1].uuid][1].output + assert isinstance(ph_doc, PhononBSDOSDoc) 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 responses[job.jobs[-1].uuid][1].output.thermal_displacement_data is None - assert isinstance(responses[job.jobs[-1].uuid][1].output.structure, Structure) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.temperatures, [0, 100, 200, 300, 400] + ph_doc.phonon_bandstructure, + PhononBS, ) - assert responses[job.jobs[-1].uuid][1].output.has_imaginary_modes is False - assert responses[job.jobs[-1].uuid][1].output.force_constants is None - 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.total_dft_energy, -5.74555232 - ) - assert responses[job.jobs[-1].uuid][1].output.born is None - assert responses[job.jobs[-1].uuid][1].output.epsilon_static is None + assert isinstance(ph_doc.phonon_dos, PhononDOS) + assert ph_doc.thermal_displacement_data is None + assert isinstance(ph_doc.structure, Structure) + + assert ph_doc.has_imaginary_modes is False + assert ph_doc.force_constants is None + assert all(isinstance(cm, CalcMeta) for cm in ph_doc.calc_meta) + + # assert_allclose( + # ph_doc.total_dft_energy, -5.74555232 + # ) + + assert ph_doc.born is None + assert ph_doc.epsilon_static is None assert_allclose( - responses[job.jobs[-1].uuid][1].output.supercell_matrix, + ph_doc.supercell_matrix, np.eye(3), ) assert_allclose( - responses[job.jobs[-1].uuid][1].output.primitive_matrix, + ph_doc.primitive_matrix, (np.ones((3, 3)) - np.eye(3)) / 2, atol=1e-8, ) - assert responses[job.jobs[-1].uuid][1].output.code == "vasp" + assert ph_doc.code == "vasp" assert isinstance( - responses[job.jobs[-1].uuid][1].output.phonopy_settings, + ph_doc.post_process_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 - == 7_000 + assert ph_doc.post_process_settings.npoints_band == 101 + assert ph_doc.post_process_settings.kpath_scheme == "seekpath" + assert ph_doc.post_process_settings.kpoint_density_dos == 7_000 + + ref_vals = { + "entropy": [0.0, 0.548554, 2.369651, 4.17177, 5.675544], + "heat_capacity": [0.0, 1.437528, 3.852217, 4.958031, 5.460526], + "internal_energy": [ + 1528.995013, + 1569.79283, + 1846.662406, + 2294.839547, + 2818.758881, + ], + "free_energy": [1528.995013, 1514.937439, 1372.73228, 1043.308596, 548.54114], + } + thermo_props = ph_doc.compute_thermo_quantities( + [0, 100, 200, 300, 400], normalization="atoms" ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.entropies, - [0.0, 2.194216, 9.478603, 16.687079, 22.702177], - atol=1e-6, - ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.heat_capacities, - [0.0, 5.750113, 15.408866, 19.832123, 21.842104], - atol=1e-6, + assert all( + thermo_props[k][i] == pytest.approx(val) + for k, vals in ref_vals.items() + for i, val in enumerate(vals) ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.internal_energies, - [6115.980051, 6279.17132, 7386.649622, 9179.358187, 11275.035523], - atol=1e-6, - ) - assert responses[job.jobs[-1].uuid][1].output.chemsys == "Si" + assert ph_doc.chemsys == "Si" # structure will be kept in the format that was transferred @@ -151,64 +147,61 @@ def test_phonon_wf_vasp_only_displacements_no_structural_transformation( responses = run_locally(job, create_folders=True, ensure_success=True) # validate settings - assert responses[job.jobs[-1].uuid][1].output.code == "vasp" + ph_doc = responses[job.jobs[-1].uuid][1].output + assert ph_doc.code == "vasp" assert isinstance( - responses[job.jobs[-1].uuid][1].output.phonopy_settings, + ph_doc.post_process_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" - ) - phonopy_settings = responses[job.jobs[-1].uuid][1].output.phonopy_settings - assert phonopy_settings.kpoint_density_dos == 7_000 + assert ph_doc.post_process_settings.npoints_band == 101 + assert ph_doc.post_process_settings.kpath_scheme == "seekpath" + post_process_settings = ph_doc.post_process_settings + assert post_process_settings.kpoint_density_dos == 7_000 # validate the outputs - assert isinstance(responses[job.jobs[-1].uuid][1].output, PhononBSDOSDoc) + assert isinstance(ph_doc, PhononBSDOSDoc) 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 responses[job.jobs[-1].uuid][1].output.thermal_displacement_data is None - assert isinstance(responses[job.jobs[-1].uuid][1].output.structure, Structure) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.temperatures, [0, 100, 200, 300, 400] - ) - assert responses[job.jobs[-1].uuid][1].output.has_imaginary_modes is False - assert responses[job.jobs[-1].uuid][1].output.force_constants is None - 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.total_dft_energy, -5.74555232 - ) - assert responses[job.jobs[-1].uuid][1].output.born is None - assert responses[job.jobs[-1].uuid][1].output.epsilon_static is None - assert responses[job.jobs[-1].uuid][1].output.supercell_matrix == tuple( - map(tuple, np.eye(3)) + ph_doc.phonon_bandstructure, + PhononBS, ) + assert isinstance(ph_doc.phonon_dos, PhononDOS) + assert ph_doc.thermal_displacement_data is None + assert isinstance(ph_doc.structure, Structure) + + assert ph_doc.has_imaginary_modes is False + assert ph_doc.force_constants is None + assert all(isinstance(cm, CalcMeta) for cm in ph_doc.calc_meta) + # assert_allclose( + # ph_doc.total_dft_energy, -5.74555232 + # ) + assert ph_doc.born is None + assert ph_doc.epsilon_static is None + assert ph_doc.supercell_matrix == tuple(map(tuple, np.eye(3))) assert_allclose( - responses[job.jobs[-1].uuid][1].output.primitive_matrix, + ph_doc.primitive_matrix, ((0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0)), atol=1e-8, ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.free_energies, - [6115.980051, 6059.749756, 5490.929122, 4173.234384, 2194.164562], - ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.entropies, - [0.0, 2.194216, 9.478603, 16.687079, 22.702177], - atol=1e-6, - ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.heat_capacities, - [0.0, 5.750113, 15.408866, 19.832123, 21.842104], + + ref_vals = { + "entropy": [0.0, 0.548554, 2.369651, 4.17177, 5.675544], + "heat_capacity": [0.0, 1.437528, 3.852217, 4.958031, 5.460526], + "internal_energy": [ + 1528.995013, + 1569.79283, + 1846.662406, + 2294.839547, + 2818.758881, + ], + "free_energy": [1528.995013, 1514.937439, 1372.73228, 1043.308596, 548.54114], + } + thermo_props = ph_doc.compute_thermo_quantities( + [0, 100, 200, 300, 400], normalization="atoms" ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.internal_energies, - [6115.980051, 6279.17132, 7386.649622, 9179.358187, 11275.035523], + assert all( + thermo_props[k][i] == pytest.approx(val) + for k, vals in ref_vals.items() + for i, val in enumerate(vals) ) @@ -247,35 +240,46 @@ def test_phonon_wf_vasp_only_displacements_kpath( ph_doc = responses[job.jobs[-1].uuid][1].output assert isinstance(ph_doc, PhononBSDOSDoc) - assert_allclose( - ph_doc.free_energies, - [5776.14995034, 5617.74737777, 4725.50269363, 3043.81827626, 694.49078355], - atol=1e-3, + ref_vals = { + "temperature": [0, 100, 200, 300, 400], + "free_energy": [ + 5776.14995034, + 5617.74737777, + 4725.50269363, + 3043.81827626, + 694.49078355, + ], + } + assert all( + ph_doc.free_energy(t) == pytest.approx(ref_vals["free_energy"][i]) + for i, t in enumerate(ref_vals["temperature"]) ) - assert isinstance(ph_doc.phonon_bandstructure, PhononBandStructureSymmLine) - assert isinstance(ph_doc.phonon_dos, PhononDos) + assert isinstance(ph_doc.phonon_bandstructure, PhononBS) + assert isinstance(ph_doc.phonon_dos, PhononDOS) assert isinstance(ph_doc.thermal_displacement_data, ThermalDisplacementData) assert isinstance(ph_doc.structure, Structure) - assert_allclose(ph_doc.temperatures, [0, 100, 200, 300, 400]) assert ph_doc.has_imaginary_modes is False - force_const = ph_doc.force_constants.force_constants[0][0][0][0] + force_const = ph_doc.force_constants[0][0][0][0] assert force_const == pytest.approx(13.032324) - assert isinstance(ph_doc.jobdirs, PhononJobDirs) - assert isinstance(ph_doc.uuids, PhononUUIDs) + assert all(isinstance(cm, CalcMeta) for cm in ph_doc.calc_meta) assert ph_doc.total_dft_energy is None assert ph_doc.born is None assert ph_doc.epsilon_static is None - assert ph_doc.supercell_matrix == ((-1, 1, 1), (1, -1, 1), (1, 1, -1)) - assert ph_doc.primitive_matrix == ((1, 0, 0), (0, 1, 0), (0, 0, 1)) + assert_allclose( + ph_doc.supercell_matrix, ((-1, 1, 1), (1, -1, 1), (1, 1, -1)), atol=1e-6 + ) + assert_allclose( + ph_doc.primitive_matrix, ((1, 0, 0), (0, 1, 0), (0, 0, 1)), atol=1e-6 + ) assert ph_doc.code == "vasp" assert isinstance( - ph_doc.phonopy_settings, + ph_doc.post_process_settings, PhononComputationalSettings, ) - assert ph_doc.phonopy_settings.npoints_band == 101 - assert ph_doc.phonopy_settings.kpath_scheme == kpath_scheme - assert ph_doc.phonopy_settings.kpoint_density_dos == 7_000 + assert ph_doc.post_process_settings.npoints_band == 101 + assert ph_doc.post_process_settings.kpath_scheme == kpath_scheme + assert ph_doc.post_process_settings.kpoint_density_dos == 7_000 # test supply of Born charges, epsilon, DFT energy, supercell @@ -358,68 +362,64 @@ def test_phonon_wf_vasp_only_displacements_add_inputs( responses = run_locally(job, create_folders=True, ensure_success=True) # validate the outputs - # print(type(responses)) - assert isinstance(responses[job.jobs[-1].uuid][1].output, PhononBSDOSDoc) + ph_doc = responses[job.jobs[-1].uuid][1].output + assert isinstance(ph_doc, PhononBSDOSDoc) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.free_energies, - [5776.14995034, 5617.74737777, 4725.50269363, 3043.81827626, 694.49078355], - atol=1e-3, + ref_vals = { + "temperature": [0, 100, 200, 300, 400], + "free_energy": [ + 5776.14995034, + 5617.74737777, + 4725.50269363, + 3043.81827626, + 694.49078355, + ], + } + assert all( + ph_doc.free_energy(t) == pytest.approx(ref_vals["free_energy"][i]) + for i, t in enumerate(ref_vals["temperature"]) ) assert isinstance( - responses[job.jobs[-1].uuid][1].output.phonon_bandstructure, - PhononBandStructureSymmLine, + ph_doc.phonon_bandstructure, + PhononBS, ) - assert isinstance(responses[job.jobs[-1].uuid][1].output.phonon_dos, PhononDos) + assert isinstance(ph_doc.phonon_dos, PhononDOS) assert isinstance( - responses[job.jobs[-1].uuid][1].output.thermal_displacement_data, + ph_doc.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, 100, 200, 300, 400] - ) - assert responses[job.jobs[-1].uuid][1].output.has_imaginary_modes is False + assert isinstance(ph_doc.structure, Structure) + assert ph_doc.has_imaginary_modes is False assert_allclose( - responses[job.jobs[-1].uuid][1].output.force_constants.force_constants[0][0][0][ - 0 - ], + ph_doc.force_constants[0][0][0][0], 13.032324, ) - 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.born, born) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.total_dft_energy, - total_dft_energy_per_formula_unit, - ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.epsilon_static, epsilon_static, atol=1e-8 - ) + + assert all(isinstance(cm, CalcMeta) for cm in ph_doc.calc_meta) + assert_allclose(ph_doc.born, born) + # assert_allclose( + # ph_doc.total_dft_energy, + # total_dft_energy_per_formula_unit, + # ) + assert_allclose(ph_doc.epsilon_static, epsilon_static, atol=1e-8) assert_allclose( - responses[job.jobs[-1].uuid][1].output.supercell_matrix, + ph_doc.supercell_matrix, [[-1, 1, 1], [1, -1, 1], [1, 1, -1]], ) assert_allclose( - responses[job.jobs[-1].uuid][1].output.primitive_matrix, + ph_doc.primitive_matrix, [[1, 0, 0], [0, 1, 0], [0, 0, 1]], atol=1e-8, ) - assert responses[job.jobs[-1].uuid][1].output.code == "vasp" + assert ph_doc.code == "vasp" assert isinstance( - responses[job.jobs[-1].uuid][1].output.phonopy_settings, + ph_doc.post_process_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 - == 7_000 - ) + assert ph_doc.post_process_settings.npoints_band == 101 + assert ph_doc.post_process_settings.kpath_scheme == "seekpath" + assert ph_doc.post_process_settings.kpoint_density_dos == 7_000 # test optional parameters @@ -451,70 +451,66 @@ def test_phonon_wf_vasp_only_displacements_optional_settings( responses = run_locally(job, create_folders=True, ensure_success=True) # validate the outputs - assert isinstance(responses[job.jobs[-1].uuid][1].output, PhononBSDOSDoc) + ph_doc = responses[job.jobs[-1].uuid][1].output + assert isinstance(ph_doc, PhononBSDOSDoc) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.free_energies, - [5776.14995034, 5617.74737777, 4725.50269363, 3043.81827626, 694.49078355], - atol=1e-3, - ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.entropies, - [0, 4.79066818, 13.03470621, 20.37400284, 26.41425489], - atol=1e-8, - ) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.heat_capacities, - [0.0, 8.05373626, 15.98005669, 19.98031234, 21.88513476], - atol=1e-8, + ref_vals = { + "free_energy": [ + 5776.14995034, + 5617.74737777, + 4725.50269363, + 3043.81827626, + 694.49078355, + ], + "entropy": [0, 4.79066818, 13.03470621, 20.37400284, 26.41425489], + "heat_capacity": [0.0, 8.05373626, 15.98005669, 19.98031234, 21.88513476], + "internal_energy": [ + 5776.14995034, + 6096.81419519, + 7332.44393529, + 9156.01912756, + 11260.1927412, + ], + } + thermo_props = ph_doc.compute_thermo_quantities( + [0, 100, 200, 300, 400], normalization="atoms" ) - - assert_allclose( - responses[job.jobs[-1].uuid][1].output.internal_energies, - [5776.14995034, 6096.81419519, 7332.44393529, 9156.01912756, 11260.1927412], - atol=1e-8, + assert all( + thermo_props[k][i] == pytest.approx(val) + for k, vals in ref_vals.items() + for i, val in enumerate(vals) ) assert isinstance( - responses[job.jobs[-1].uuid][1].output.phonon_bandstructure, - PhononBandStructureSymmLine, + ph_doc.phonon_bandstructure, + PhononBS, ) - assert isinstance(responses[job.jobs[-1].uuid][1].output.phonon_dos, PhononDos) - assert responses[job.jobs[-1].uuid][1].output.thermal_displacement_data is None - assert isinstance(responses[job.jobs[-1].uuid][1].output.structure, Structure) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.temperatures, [0, 100, 200, 300, 400] - ) - assert responses[job.jobs[-1].uuid][1].output.has_imaginary_modes is False - assert responses[job.jobs[-1].uuid][1].output.force_constants is None - assert isinstance(responses[job.jobs[-1].uuid][1].output.jobdirs, PhononJobDirs) - assert isinstance(responses[job.jobs[-1].uuid][1].output.uuids, PhononUUIDs) - assert responses[job.jobs[-1].uuid][1].output.total_dft_energy is None - assert responses[job.jobs[-1].uuid][1].output.born is None - assert responses[job.jobs[-1].uuid][1].output.epsilon_static is None + assert isinstance(ph_doc.phonon_dos, PhononDOS) + assert ph_doc.thermal_displacement_data is None + assert isinstance(ph_doc.structure, Structure) + assert ph_doc.has_imaginary_modes is False + assert ph_doc.force_constants is None + assert all(isinstance(cm, CalcMeta) for cm in ph_doc.calc_meta) + assert ph_doc.total_dft_energy is None + assert ph_doc.born is None + assert ph_doc.epsilon_static is None assert_allclose( - responses[job.jobs[-1].uuid][1].output.supercell_matrix, + ph_doc.supercell_matrix, [[-1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, -1.0]], ) assert_allclose( - responses[job.jobs[-1].uuid][1].output.primitive_matrix, + ph_doc.primitive_matrix, [[1, 0, 0], [0, 1, 0], [0, 0, 1]], atol=1e-8, ) - assert responses[job.jobs[-1].uuid][1].output.code == "vasp" + assert ph_doc.code == "vasp" assert isinstance( - responses[job.jobs[-1].uuid][1].output.phonopy_settings, + ph_doc.post_process_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 - == 7_000 - ) + assert ph_doc.post_process_settings.npoints_band == 101 + assert ph_doc.post_process_settings.kpath_scheme == "seekpath" + assert ph_doc.post_process_settings.kpoint_density_dos == 7_000 # test run including all steps of the computation for Si @@ -552,66 +548,66 @@ def test_phonon_wf_vasp_all_steps(mock_vasp, clean_dir, si_structure: Structure) responses = run_locally(job, create_folders=True, ensure_success=True) # validate the outputs - assert isinstance(responses[job.jobs[-1].uuid][1].output, PhononBSDOSDoc) + ph_doc = responses[job.jobs[-1].uuid][1].output + assert isinstance(ph_doc, PhononBSDOSDoc) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.free_energies, - [5853.74150399, 5692.29089555, 4798.67784919, 3122.48296003, 782.17345333], + ref_vals = { + "temperature": [0, 100, 200, 300, 400], + "free_energy": [ + 5853.74150399, + 5692.29089555, + 4798.67784919, + 3122.48296003, + 782.17345333, + ], + } + assert all( + ph_doc.free_energy(t) == pytest.approx(ref_vals["free_energy"][i]) + for i, t in enumerate(ref_vals["temperature"]) ) assert isinstance( - responses[job.jobs[-1].uuid][1].output.phonon_bandstructure, - PhononBandStructureSymmLine, + ph_doc.phonon_bandstructure, + PhononBS, ) - assert isinstance(responses[job.jobs[-1].uuid][1].output.phonon_dos, PhononDos) + assert isinstance(ph_doc.phonon_dos, PhononDOS) assert isinstance( - responses[job.jobs[-1].uuid][1].output.thermal_displacement_data, + ph_doc.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, 100, 200, 300, 400] - ) - assert responses[job.jobs[-1].uuid][1].output.has_imaginary_modes is False + assert isinstance(ph_doc.structure, Structure) + assert ph_doc.has_imaginary_modes is False assert_allclose( - responses[job.jobs[-1].uuid][1].output.force_constants.force_constants[0][0][0][ - 0 - ], + ph_doc.force_constants[0][0][0][0], 13.41185599, ) - 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.born, np.zeros((2, 3, 3))) - assert_allclose( - responses[job.jobs[-1].uuid][1].output.total_dft_energy, -5.74629058 - ) + + assert all(isinstance(cm, CalcMeta) for cm in ph_doc.calc_meta) + assert_allclose(ph_doc.born, np.zeros((2, 3, 3))) + # assert_allclose( + # ph_doc.total_dft_energy, -5.74629058 + # ) assert_allclose( - responses[job.jobs[-1].uuid][1].output.epsilon_static, + ph_doc.epsilon_static, ((13.19242034, -0.0, 0.0), (-0.0, 13.19242034, 0.0), (0.0, 0.0, 13.19242034)), atol=1e-8, ) assert_allclose( - responses[job.jobs[-1].uuid][1].output.supercell_matrix, + ph_doc.supercell_matrix, [[-1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [1.0, 1.0, -1.0]], ) assert_allclose( - responses[job.jobs[-1].uuid][1].output.primitive_matrix, + ph_doc.primitive_matrix, [[1, 0, 0], [0, 1, 0], [0, 0, 1]], ) - assert responses[job.jobs[-1].uuid][1].output.code == "vasp" + assert ph_doc.code == "vasp" assert isinstance( - responses[job.jobs[-1].uuid][1].output.phonopy_settings, + ph_doc.post_process_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 - == 7_000 - ) + assert ph_doc.post_process_settings.npoints_band == 101 + assert ph_doc.post_process_settings.kpath_scheme == "seekpath" + assert ph_doc.post_process_settings.kpoint_density_dos == 7_000 # use a structure where Born charges are actually useful for the computation and change diff --git a/tests/vasp/flows/test_qha.py b/tests/vasp/flows/test_qha.py index 9c8b4f6670..680a74e0ed 100644 --- a/tests/vasp/flows/test_qha.py +++ b/tests/vasp/flows/test_qha.py @@ -35,4 +35,4 @@ def test_qha(mock_vasp, clean_dir, si_diamond: Structure): # run the flow or job and ensure that it finished running successfully responses = run_locally(qha_maker, create_folders=True, ensure_success=True) - assert len(responses) == 9 + assert len(responses) == 10 diff --git a/tutorials/grueneisen_workflow.ipynb b/tutorials/grueneisen_workflow.ipynb index 77d4549b11..7721f530a4 100644 --- a/tutorials/grueneisen_workflow.ipynb +++ b/tutorials/grueneisen_workflow.ipynb @@ -325,7 +325,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.16" + "version": "3.12.0" } }, "nbformat": 4, diff --git a/tutorials/phonon_workflow.ipynb b/tutorials/phonon_workflow.ipynb index 3715885b1d..0668b762fe 100644 --- a/tutorials/phonon_workflow.ipynb +++ b/tutorials/phonon_workflow.ipynb @@ -158,8 +158,6 @@ "metadata": {}, "outputs": [], "source": [ - "from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine\n", - "from pymatgen.phonon.dos import PhononDos\n", "from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter\n", "\n", "job_store.connect()\n", @@ -182,12 +180,14 @@ "metadata": {}, "outputs": [], "source": [ - "ph_bs = PhononBandStructureSymmLine.from_dict(\n", - " result[\"output\"][\"phonon_bandstructure\"]\n", - ") # get pymatgen bandstructure object\n", - "ph_dos = PhononDos.from_dict(\n", - " result[\"output\"][\"phonon_dos\"]\n", - ") # get pymatgen phonon dos object\n", + "from emmet.core.phonon import PhononBS, PhononDOS\n", + "\n", + "ph_bs = PhononBS(\n", + " **result[\"output\"][\"phonon_bandstructure\"]\n", + ").to_pmg # get pymatgen bandstructure object\n", + "ph_dos = PhononDOS(\n", + " **result[\"output\"][\"phonon_dos\"]\n", + ").to_pmg # get pymatgen phonon dos object\n", "\n", "# initialize dos plotter and visualize dos plot\n", "dos_plot = PhononDosPlotter()\n", diff --git a/tutorials/phonon_workflow_aims.ipynb b/tutorials/phonon_workflow_aims.ipynb index 18bc13445f..bd2b639460 100644 --- a/tutorials/phonon_workflow_aims.ipynb +++ b/tutorials/phonon_workflow_aims.ipynb @@ -159,8 +159,6 @@ "metadata": {}, "outputs": [], "source": [ - "from pymatgen.phonon.bandstructure import PhononBandStructureSymmLine\n", - "from pymatgen.phonon.dos import PhononDos\n", "from pymatgen.phonon.plotter import PhononBSPlotter, PhononDosPlotter\n", "\n", "job_store.connect()\n", @@ -179,12 +177,14 @@ "metadata": {}, "outputs": [], "source": [ - "ph_bs = PhononBandStructureSymmLine.from_dict(\n", - " result[\"output\"][\"phonon_bandstructure\"]\n", - ") # get pymatgen bandstructure object\n", - "ph_dos = PhononDos.from_dict(\n", - " result[\"output\"][\"phonon_dos\"]\n", - ") # get pymatgen phonon dos object\n", + "from emmet.core.phonon import PhononBS, PhononDOS\n", + "\n", + "ph_bs = PhononBS(\n", + " **result[\"output\"][\"phonon_bandstructure\"]\n", + ").to_pmg # get pymatgen bandstructure object\n", + "ph_dos = PhononDOS(\n", + " **result[\"output\"][\"phonon_dos\"]\n", + ").to_pmg # get pymatgen phonon dos object\n", "\n", "# initialize dos plotter and visualize dos plot\n", "dos_plot = PhononDosPlotter()\n", @@ -344,7 +344,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.16" + "version": "3.12.0" } }, "nbformat": 4, diff --git a/tutorials/qha_workflow.ipynb b/tutorials/qha_workflow.ipynb index 97f4915daa..b15ee73021 100644 --- a/tutorials/qha_workflow.ipynb +++ b/tutorials/qha_workflow.ipynb @@ -375,7 +375,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.16" + "version": "3.12.0" } }, "nbformat": 4,