diff --git a/docs/about/contributors.md b/docs/about/contributors.md index b5dbf78fd3..3dc99f7083 100644 --- a/docs/about/contributors.md +++ b/docs/about/contributors.md @@ -28,6 +28,7 @@ Additional contributions were made by the individuals listed [here](https://gith - [@yw-fang](https://github.com/yw-fang): VASP Non-SCF recipe - [@espottesmith](http://github.com/espottesmith): Some ORCA and Q-Chem workflows, mostly re: IRC - [@honghuikim](https://github.com/honghuikim): Dynamic changing of quacc settings with workflow engines +- [@benshi97](https://github.com/benshi97): VASP frequency recipe ## Inspiration diff --git a/docs/user/recipes/recipes_list.md b/docs/user/recipes/recipes_list.md index d4f287b8e6..4530df72db 100644 --- a/docs/user/recipes/recipes_list.md +++ b/docs/user/recipes/recipes_list.md @@ -251,6 +251,7 @@ The list of available quacc recipes is shown below. The "Req'd Extras" column sp | VASP Non-SCF | `#!Python @job` | [quacc.recipes.vasp.core.non_scf_job][] | | | VASP Slab Static | `#!Python @job` | [quacc.recipes.vasp.slabs.static_job][] | | | VASP Slab Relax | `#!Python @job` | [quacc.recipes.vasp.slabs.relax_job][] | | +| VASP Frequency | `#!Python @job` | [quacc.recipes.vasp.core.freq_job][] | | | VASP Bulk to Slabs | `#!Python @flow` | [quacc.recipes.vasp.slabs.bulk_to_slabs_flow][] | | | VASP Slab to Adsorbates | `#!Python @flow` | [quacc.recipes.vasp.slabs.slab_to_ads_flow][] | | | VASP MP GGA Relax | `#!Python @job` | [quacc.recipes.vasp.mp.mp_gga_relax_job][] | `quacc[mp]` | diff --git a/src/quacc/recipes/vasp/_base.py b/src/quacc/recipes/vasp/_base.py index 0b88d583cf..59ea4809a7 100644 --- a/src/quacc/recipes/vasp/_base.py +++ b/src/quacc/recipes/vasp/_base.py @@ -6,11 +6,13 @@ from quacc.calculators.vasp import Vasp from quacc.runners.ase import Runner +from quacc.runners.thermo import ThermoRunner +from quacc.schemas.ase import summarize_vib_and_thermo from quacc.schemas.vasp import summarize_vasp_opt_run, vasp_summarize_run from quacc.utils.dicts import recursive_dict_merge if TYPE_CHECKING: - from typing import Any + from typing import Any, Literal from ase.atoms import Atoms @@ -20,6 +22,7 @@ SourceDirectory, VaspASEOptSchema, VaspSchema, + VibThermoSchema, ) @@ -124,3 +127,76 @@ def run_and_summarize_opt( report_mp_corrections=report_mp_corrections, additional_fields=additional_fields, ) + + +def run_and_summarize_vib_and_thermo( + atoms: Atoms, + energy: float = 0.0, + temperature: float = 298.15, + pressure: float = 1.0, + thermo_method: Literal["ideal_gas", "harmonic"] = "ideal_gas", + preset: str | None = None, + calc_defaults: dict[str, Any] | None = None, + calc_swaps: dict[str, Any] | None = None, + vib_kwargs: dict[str, Any] | None = None, + additional_fields: dict[str, Any] | None = None, + copy_files: SourceDirectory | dict[SourceDirectory, Filenames] | None = None, +) -> VibThermoSchema: + """ + Base job function for VASP recipes with ASE vibrational analysis. + + Parameters + ---------- + atoms + Atoms object + energy + Energy of the system + temperature + Temperature of the system + pressure + Pressure of the system + thermo_method + Method to use for thermochemistry. Options are "harmonic" or "ideal_gas". + preset + Preset to use from `quacc.calculators.vasp.presets`. + calc_defaults + Default parameters for the recipe. + calc_swaps + Dictionary of custom kwargs for the Vasp calculator. Set a value to + `None` to remove a pre-existing key entirely. For a list of available + keys, refer to [quacc.calculators.vasp.vasp.Vasp][]. + vib_kwargs + Dictionary of custom kwargs for [quacc.runners.ase.Runner.run_vib][] + additional_fields + Additional fields to supply to the summarizer. + copy_files + Files to copy (and decompress) from source to the runtime directory. + + Returns + ------- + VibThermoSchema + Dictionary of results, specified in [quacc.schemas.ase.summarize_vib_and_thermo][]. + See the type-hint for the data structure. + """ + + # Set defaults + calc_flags = recursive_dict_merge(calc_defaults, calc_swaps) + + calc = Vasp(atoms, preset=preset, **calc_flags) + vibrations = Runner(atoms, calc, copy_files=copy_files).run_vib( + vib_kwargs=vib_kwargs + ) + thermo_runner = ThermoRunner(atoms, vibrations.get_frequencies(), energy=energy) + if thermo_method == "harmonic": + thermo_object = thermo_runner.run_harmonic_thermo() + elif thermo_method == "ideal_gas": + thermo_object = thermo_runner.run_ideal_gas() + + return summarize_vib_and_thermo( + vibrations, + thermo_object, + atoms=atoms, + temperature=temperature, + pressure=pressure, + additional_fields=additional_fields, + ) diff --git a/src/quacc/recipes/vasp/core.py b/src/quacc/recipes/vasp/core.py index d1fb2fc946..52797c73c5 100644 --- a/src/quacc/recipes/vasp/core.py +++ b/src/quacc/recipes/vasp/core.py @@ -11,7 +11,11 @@ from pymatgen.io.vasp import Vasprun from quacc import flow, job -from quacc.recipes.vasp._base import run_and_summarize, run_and_summarize_opt +from quacc.recipes.vasp._base import ( + run_and_summarize, + run_and_summarize_opt, + run_and_summarize_vib_and_thermo, +) if TYPE_CHECKING: from typing import Any @@ -25,6 +29,8 @@ SourceDirectory, VaspASEOptSchema, VaspSchema, + VibKwargs, + VibThermoSchema, ) @@ -335,3 +341,65 @@ def non_scf_job( additional_fields={"name": "VASP Non-SCF"}, copy_files={prev_dir: ["CHGCAR*", "WAVECAR*"]}, ) + + +@job +def freq_job( + atoms: Atoms, + preset: str | None = "BulkSet", + energy: float = 0.0, + temperature: float = 298.15, + pressure: float = 1.0, + thermo_method: Literal["harmonic", "ideal_gas"] = "ideal_gas", + vib_kwargs: VibKwargs | None = None, + copy_files: SourceDirectory | dict[SourceDirectory, Filenames] | None = None, + **calc_kwargs, +) -> VibThermoSchema: + """ + Run a frequency job and calculate thermochemistry. + + Parameters + ---------- + atoms + Atoms object + preset + Preset to use from `quacc.calculators.vasp.presets`. + energy + Potential energy in eV. If 0, then the output is just the correction. + temperature + Temperature in Kelvins. + pressure + Pressure in bar. + thermo_method + Method to use for thermochemistry. Options are "harmonic" or "ideal_gas". + vib_kwargs + Dictionary of kwargs for the [ase.vibrations.Vibrations][] class. + copy_files + Files to copy (and decompress) from source to the runtime directory. + **calc_kwargs + Custom kwargs for the Vasp calculator. Set a value to + `None` to remove a pre-existing key entirely. For a list of available + keys, refer to [quacc.calculators.vasp.vasp.Vasp][]. + + Returns + ------- + VibThermoSchema + Dictionary of results, specified in [quacc.schemas.ase.summarize_vib_and_thermo][]. + See the type-hint for the data structure. + """ + calc_defaults = {"ediff": 1e-7, "isym": 0, "lcharg": False, "lwave": True, "nsw": 0} + vib_kwargs = vib_kwargs or {} + + return run_and_summarize_vib_and_thermo( + atoms, + energy=energy, + temperature=temperature, + pressure=pressure, + thermo_method=thermo_method, + preset=preset, + calc_defaults=calc_defaults, + calc_swaps=calc_kwargs, + vib_kwargs=vib_kwargs, + copy_files=copy_files, + additional_fields={"name": "VASP Frequency and Thermo"}, + ) diff --git a/src/quacc/runners/thermo.py b/src/quacc/runners/thermo.py index 3a9c2c9034..b8fe8680a4 100644 --- a/src/quacc/runners/thermo.py +++ b/src/quacc/runners/thermo.py @@ -7,7 +7,7 @@ import numpy as np from ase import units -from ase.thermochemistry import IdealGasThermo +from ase.thermochemistry import HarmonicThermo, IdealGasThermo from emmet.core.symmetry import PointGroupData from pymatgen.io.ase import AseAtomsAdaptor @@ -53,6 +53,7 @@ def run_ideal_gas(self, spin_multiplicity: int | None = None) -> IdealGasThermo: ------- IdealGasThermo object """ + # Ensure all negative modes are made complex for i, f in enumerate(self.vib_freqs): if not isinstance(f, complex) and f < 0: @@ -82,11 +83,11 @@ def run_ideal_gas(self, spin_multiplicity: int | None = None) -> IdealGasThermo: spin = 0 # Get symmetry for later use - natoms = len(self.atoms) mol = AseAtomsAdaptor().get_molecule(self.atoms, charge_spin_check=False) point_group_data = PointGroupData().from_molecule(mol) # Get the geometry + natoms = len(self.atoms) if natoms == 1: geometry = "monatomic" elif point_group_data.linear: @@ -103,3 +104,28 @@ def run_ideal_gas(self, spin_multiplicity: int | None = None) -> IdealGasThermo: spin=spin, ignore_imag_modes=True, ) + + def run_harmonic_thermo(self) -> HarmonicThermo: + """ + Create a HarmonicThermo object for a molecule from a given vibrational analysis. + This works for free gases, solids and adsorbates. + + Returns + ------- + HarmonicThermo object + The ASE `HarmonicThermo` object. + """ + + # Ensure all negative modes are made complex + for i, f in enumerate(self.vib_freqs): + if not isinstance(f, complex) and f < 0: + self.vib_freqs[i] = complex(0 - f * 1j) + + # Convert vibrational frequencies to energies + vib_energies = [f * units.invcm for f in self.vib_freqs] + + return HarmonicThermo( + vib_energies=vib_energies, + potentialenergy=self.energy, + ignore_imag_modes=True, + ) diff --git a/src/quacc/schemas/ase.py b/src/quacc/schemas/ase.py index eb6b3663cb..37e97e6537 100644 --- a/src/quacc/schemas/ase.py +++ b/src/quacc/schemas/ase.py @@ -7,6 +7,7 @@ import numpy as np from ase import units from ase.io import read +from ase.thermochemistry import HarmonicThermo, IdealGasThermo from ase.vibrations import Vibrations from ase.vibrations.data import VibrationsData @@ -24,7 +25,6 @@ from ase.io import Trajectory from ase.md.md import MolecularDynamics from ase.optimize.optimize import Optimizer - from ase.thermochemistry import IdealGasThermo from maggma.core import Store from quacc.types import ( @@ -288,7 +288,8 @@ def summarize_md_run( def summarize_vib_and_thermo( vib: Vibrations, - igt: IdealGasThermo, + thermo_object: IdealGasThermo | HarmonicThermo, + atoms: Atoms | None = None, temperature: float = 298.15, pressure: float = 1.0, charge_and_multiplicity: tuple[int, int] | None = None, @@ -303,8 +304,10 @@ def summarize_vib_and_thermo( ---------- vib ASE Vibrations object. - igt - ASE IdealGasThermo object. + thermo_object + ASE IdealGasThermo or HarmonicThermo object. + atoms + ASE Atoms object following a calculation. temperature Temperature in Kelvins. pressure @@ -328,12 +331,18 @@ def summarize_vib_and_thermo( vib_task_doc = _summarize_vib_run( vib, charge_and_multiplicity=charge_and_multiplicity ) - thermo_task_doc = _summarize_ideal_gas_thermo( - igt, - temperature=temperature, - pressure=pressure, - charge_and_multiplicity=charge_and_multiplicity, - ) + + if isinstance(thermo_object, HarmonicThermo): + thermo_task_doc = _summarize_harmonic_thermo( + atoms, thermo_object, temperature=temperature, pressure=pressure + ) + elif isinstance(thermo_object, IdealGasThermo): + thermo_task_doc = _summarize_ideal_gas_thermo( + thermo_object, + temperature=temperature, + pressure=pressure, + charge_and_multiplicity=charge_and_multiplicity, + ) unsorted_task_doc = recursive_dict_merge( vib_task_doc, thermo_task_doc, additional_fields @@ -510,3 +519,59 @@ def _summarize_ideal_gas_thermo( ) return atoms_metadata | inputs | results + + +def _summarize_harmonic_thermo( + atoms: Atoms, + harmonic_thermo: HarmonicThermo, + temperature: float = 298.15, + pressure: float = 1.0, +) -> ThermoSchema: + """ + Get tabulated results from an ASE HarmonicThermo object and store them in a + database-friendly format. + + Parameters + ---------- + atoms + ASE Atoms object used for the vibrational frequency calculation. + harmonic_thermo + ASE HarmonicThermo object. + temperature + Temperature in Kelvins. + pressure + Pressure in bar. + + Returns + ------- + ThermoSchema + Dictionary representation of the task document + """ + settings = get_settings() + + inputs = { + "parameters_thermo": { + "temperature": temperature, + "pressure": pressure, + "vib_freqs": [e / units.invcm for e in harmonic_thermo.vib_energies], + "vib_energies": harmonic_thermo.vib_energies.tolist(), + "n_imag": harmonic_thermo.n_imag, + } + } + + results = { + "results": { + "energy": harmonic_thermo.potentialenergy, + "helmholtz_energy": harmonic_thermo.get_helmholtz_energy( + temperature, verbose=settings.DEBUG + ), + "internal_energy": harmonic_thermo.get_internal_energy( + temperature, verbose=settings.DEBUG + ), + "entropy": harmonic_thermo.get_entropy(temperature, verbose=settings.DEBUG), + "zpe": harmonic_thermo.get_ZPE_correction(), + } + } + + atoms_metadata = atoms_to_metadata(atoms) + return atoms_metadata | inputs | results diff --git a/tests/core/recipes/vasp_recipes/mocked/test_vasp_recipes.py b/tests/core/recipes/vasp_recipes/mocked/test_vasp_recipes.py index 1a65211ed8..443d1f08a0 100644 --- a/tests/core/recipes/vasp_recipes/mocked/test_vasp_recipes.py +++ b/tests/core/recipes/vasp_recipes/mocked/test_vasp_recipes.py @@ -21,6 +21,7 @@ from quacc.recipes.vasp.core import ( ase_relax_job, double_relax_flow, + freq_job, non_scf_job, relax_job, static_job, @@ -849,3 +850,21 @@ def test_mp_relax_flow_custom(tmp_path, patch_nonmetallic_taskdoc): job_params={"mp_metagga_relax_job": {"nsw": 0}}, ) assert output["relax2"]["parameters"]["nsw"] == 0 + + +def test_freq_job(): + atoms = molecule("H2") + atoms.pbc = True + atoms.center(vacuum=1) + + output = freq_job(atoms, kpts=(1, 1, 1), thermo_method="harmonic") + assert output["parameters"]["ediff"] == 1e-07 + # Check that "sigma" (only used in ideal_gas) isn't a key in parameters_thermo + assert "sigma" not in output["parameters_thermo"] + assert len(output["results"]["vib_freqs_raw"]) == 3 * len(atoms) + + output = freq_job(atoms, kpts=(1, 1, 1)) + assert output["parameters"]["ediff"] == 1e-07 + assert output["parameters_thermo"]["sigma"] == 2.0 + assert len(output["results"]["vib_freqs_raw"]) == 3 * len(atoms) + assert len(output["results"]["vib_freqs"]) == 3 * len(atoms) - 5 diff --git a/tests/core/runners/test_thermo.py b/tests/core/runners/test_thermo.py index b81d02b235..722cd36888 100644 --- a/tests/core/runners/test_thermo.py +++ b/tests/core/runners/test_thermo.py @@ -26,3 +26,13 @@ def test_run_ideal_gas(): co2.calc.results["magmom"] = 1.0 igt = ThermoRunner(co2, [-12, 526, 526, 1480, 2565]).run_ideal_gas() assert igt.get_ZPE_correction() == pytest.approx(2548.5 * invcm) + + +def test_run_harmonic_thermo(): + co2 = molecule("CO2") + ht = ThermoRunner(co2, [526, 526, 1480, 2565]).run_harmonic_thermo() + assert ht.get_ZPE_correction() == pytest.approx(2548.5 * invcm) + + co2 = molecule("CO2") + ht = ThermoRunner(co2, [-12, 526, 526, 1480, 2565]).run_harmonic_thermo() + assert ht.get_ZPE_correction() == pytest.approx(2548.5 * invcm) diff --git a/tests/core/schemas/test_ase.py b/tests/core/schemas/test_ase.py index 34b98fd284..80568c9e36 100644 --- a/tests/core/schemas/test_ase.py +++ b/tests/core/schemas/test_ase.py @@ -9,18 +9,21 @@ from ase.calculators.emt import EMT from ase.io import read from ase.optimize import BFGS -from ase.thermochemistry import IdealGasThermo +from ase.thermochemistry import HarmonicThermo, IdealGasThermo from ase.units import invcm from ase.vibrations import Vibrations from maggma.stores import MemoryStore from monty.json import MontyDecoder, jsanitize from monty.serialization import loadfn +from numpy.testing import assert_allclose from quacc.schemas.ase import ( + _summarize_harmonic_thermo, _summarize_ideal_gas_thermo, _summarize_vib_run, summarize_opt_run, summarize_run, + summarize_vib_and_thermo, ) FILE_DIR = Path(__file__).parent @@ -345,6 +348,70 @@ def test_summarize_ideal_gas_thermo(tmp_path, monkeypatch): _summarize_ideal_gas_thermo(igt, charge_and_multiplicity=[0, 1]) +def test_summarize_harmonic_thermo(tmp_path, monkeypatch): + monkeypatch.chdir(tmp_path) + atoms = molecule("H2") + + # Make sure metadata is made + ht = HarmonicThermo([0.34]) + results = _summarize_harmonic_thermo(atoms, ht) + assert results["parameters_thermo"]["vib_energies"] == [0.34] + assert results["parameters_thermo"]["vib_freqs"] == [0.34 / invcm] + assert results["results"]["energy"] == 0 + assert_allclose( + results["results"]["helmholtz_energy"], 0.16999995401497991, rtol=1e-5 + ) + assert_allclose( + results["results"]["internal_energy"], 0.1700006085385999, rtol=1e-5 + ) + assert_allclose(results["results"]["entropy"], 2.1952829783392438e-09, rtol=1e-5) + assert_allclose(results["results"]["zpe"], 0.17, rtol=1e-5) + + # test document can be jsanitized and decoded + d = jsanitize(results, strict=True, enum_values=True) + MontyDecoder().process_decoded(d) + + +def test_summarize_vib_and_thermo(tmp_path, monkeypatch): + monkeypatch.chdir(tmp_path) + + # Test harmonic thermo + atoms = molecule("H2") + atoms.calc = EMT() + + ht = HarmonicThermo([0.34]) + vib = Vibrations(atoms) + vib.run() + results = summarize_vib_and_thermo(vib, ht, atoms=atoms) + + assert results["parameters_thermo"]["vib_energies"] == [0.34] + assert results["parameters_thermo"]["vib_freqs"] == [0.34 / invcm] + assert results["results"]["energy"] == 0 + assert_allclose( + results["results"]["helmholtz_energy"], 0.16999995401497991, rtol=1e-5 + ) + assert_allclose( + results["results"]["internal_energy"], 0.1700006085385999, rtol=1e-5 + ) + assert_allclose(results["results"]["entropy"], 2.1952829783392438e-09, rtol=1e-5) + assert_allclose(results["results"]["zpe"], 0.17, rtol=1e-5) + + # Test ideal gas thermo + atoms = molecule("N2") + atoms.calc = EMT() + igt = IdealGasThermo([0.34], "linear", atoms=atoms, spin=0, symmetrynumber=2) + vib = Vibrations(atoms) + vib.run() + results = summarize_vib_and_thermo(vib, igt, atoms=atoms) + + assert results["natoms"] == len(atoms) + assert results["atoms"] == atoms + assert results["parameters_thermo"]["vib_energies"] == [0.34] + assert results["parameters_thermo"]["vib_freqs"] == [0.34 / invcm] + assert results["results"]["energy"] == 0 + assert "pymatgen_version" in results["builder_meta"] + + def test_errors(tmp_path, monkeypatch): monkeypatch.chdir(tmp_path)