Vibrational analysis recipe for VASP

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
 

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]` |

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,
+ )

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"},
+ )

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,
+ )

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