Workflow for Quasi-harmonic approximation (forcefields and VASP)

src/atomate2/common/flows/eos.py

@@ -81,7 +81,9 @@ def make(self, structure: Structure, prev_dir: str | Path = None) -> Flow:
             ("relax", "static") if self.static_maker else ("relax",)
         )
         flow_output: dict[str, dict] = {
-            key: {quantity: [] for quantity in ("energy", "volume", "stress")}
+            key: {
+                quantity: [] for quantity in ("energy", "volume", "stress", "structure")
+            }
             for key in job_types
         }
 
@@ -166,6 +168,8 @@ def make(self, structure: Structure, prev_dir: str | Path = None) -> Flow:
                 flow_output[key]["energy"] += [output.energy]
                 flow_output[key]["volume"] += [output.structure.volume]
                 flow_output[key]["stress"] += [output.stress]
+                # TODO: make this potentially optional?
+                flow_output[key]["structure"] += [output.structure]
 
         if self.postprocessor is not None:
             min_points = self.postprocessor.min_data_points

src/atomate2/common/flows/qha.py

@@ -0,0 +1,179 @@
+"""Define common QHA flow agnostic to electronic-structure code."""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from dataclasses import dataclass, field
+from typing import TYPE_CHECKING, Literal
+
+from jobflow import Flow, Maker
+
+from atomate2.common.flows.eos import CommonEosMaker
+from atomate2.common.jobs.qha import analyze_free_energy, get_phonon_jobs
+
+if TYPE_CHECKING:
+    from pathlib import Path
+
+    from pymatgen.core import Structure
+
+    from atomate2.common.flows.phonons import BasePhononMaker
+    from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker
+    from atomate2.vasp.jobs.base import BaseVaspMaker
+
+supported_eos = frozenset(("vinet", "birch_murnaghan", "murnaghan"))
+
+
+@dataclass
+class CommonQhaMaker(Maker, ABC):
+    """
+    Use the quasi-harmonic approximation.
+
+    First relax a structure.
+    Then we scale the relaxed structure, and
+    then compute harmonic phonons for each scaled
+    structure with Phonopy.
+    Finally, we compute the Gibb's free energy and
+    other thermodynamic properties available from
+    the quasi-harmonic approximation.
+
+    Note: We do not consider electronic free energies so far.
+    This might be problematic for metals (see e.g.,
+    Wolverton and Zunger, Phys. Rev. B, 52, 8813 (1994).)
+
+    Note: Magnetic Materials have never been computed with
+    this workflow.
+
+    Parameters
+    ----------
+    name: str
+        Name of the flows produced by this maker.
+    initial_relax_maker: .ForceFieldRelaxMaker | .BaseVaspMaker | None
+        Maker to relax the input structure.
+    eos_relax_maker: .ForceFieldRelaxMaker | .BaseVaspMaker | None
+        Maker to relax deformed structures for the EOS fit.
+        The volume has to be fixed!
+    phonon_displacement_maker: .ForceFieldStaticMaker | .BaseVaspMaker | None
+    phonon_static_maker: .ForceFieldStaticMaker | .BaseVaspMaker | None
+    phonon_maker_kwargs: dict
+    linear_strain: tuple[float, float]
+        Percentage linear strain to apply as a deformation, default = -5% to 5%.
+    number_of_frames: int
+        Number of strain calculations to do for EOS fit, default = 6.
+    t_max: float | None
+        Maximum temperature until which the QHA will be performed
+    pressure: float | None
+        Pressure at which the QHA will be performed (default None, no pressure)
+    ignore_imaginary_modes: bool
+        By default, volumes where the harmonic phonon approximation shows imaginary
+        will be ignored
+    eos_type: str
+        Equation of State type used for the fitting. Defaults to vinet.
+    """
+
+    name: str = "QHA Maker"
+    initial_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None
+    eos_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None
+    phonon_displacement_maker: ForceFieldStaticMaker | BaseVaspMaker | None = None
+    phonon_static_maker: ForceFieldStaticMaker | BaseVaspMaker | None = None
+    phonon_maker_kwargs: dict = field(default_factory=dict)
+    linear_strain: tuple[float, float] = (-0.05, 0.05)
+    number_of_frames: int = 6
+    t_max: float | None = None
+    pressure: float | None = None
+    ignore_imaginary_modes: bool = False
+    eos_type: Literal["vinet", "birch_murnaghan", "murnaghan"] = "vinet"
+    analyze_free_energy_kwargs: dict = field(default_factory=dict)
+    # TODO: implement advanced handling of
+    #  imaginary modes in phonon runs (i.e., fitting procedures)
+
+    def make(self, structure: Structure, prev_dir: str | Path = None) -> Flow:
+        """Run an EOS flow.
+
+        Parameters
+        ----------
+        structure : Structure
+            A pymatgen structure object.
+        prev_dir : str or Path or None
+            A previous calculation directory to copy output files from.
+
+        Returns
+        -------
+        .Flow, a QHA flow
+        """
+        if self.eos_type not in supported_eos:
+            raise ValueError(
+                "EOS not supported.",
+                "Please choose 'vinet', 'birch_murnaghan', 'murnaghan'",
+            )
+
+        # In this way, one can easily exchange makers and enforce postprocessor None
+        self.eos = CommonEosMaker(
+            initial_relax_maker=self.initial_relax_maker,
+            eos_relax_maker=self.eos_relax_maker,
+            static_maker=None,
+            postprocessor=None,
+            number_of_frames=self.number_of_frames,
+        )
+        self.phonon_maker = self.initialize_phonon_maker(
+            phonon_displacement_maker=self.phonon_displacement_maker,
+            phonon_static_maker=self.phonon_static_maker,
+            bulk_relax_maker=None,
+            phonon_maker_kwargs=self.phonon_maker_kwargs,
+        )
+        eos_job = self.eos.make(structure)
+        phonon_jobs = get_phonon_jobs(
+            phonon_maker=self.phonon_maker, eos_output=eos_job.output
+        )
+        analysis = analyze_free_energy(
+            phonon_jobs.output,
+            structure=structure,
+            t_max=self.t_max,
+            pressure=self.pressure,
+            ignore_imaginary_modes=self.ignore_imaginary_modes,
+            eos_type=self.eos_type,
+            **self.analyze_free_energy_kwargs,
+        )
+
+        return Flow([eos_job, phonon_jobs, analysis])
+
+    @abstractmethod
+    def initialize_phonon_maker(
+        self,
+        phonon_displacement_maker: ForceFieldStaticMaker | BaseVaspMaker | None,
+        phonon_static_maker: ForceFieldStaticMaker | BaseVaspMaker | None,
+        bulk_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | None,
+        phonon_maker_kwargs: dict,
+    ) -> BasePhononMaker | None:
+        """Initialize phonon maker.
+
+        This implementation will be different for
+        any newly implemented QHAMaker.
+
+        Parameters
+        ----------
+        phonon_displacement_maker: ForceFieldStaticMaker|BaseVaspMaker|None
+            Maker for displacement calculations.
+        phonon_static_maker: ForceFieldStaticMaker|BaseVaspMaker|None
+            Maker for additional static calculations.
+        bulk_relax_maker: : ForceFieldRelaxMaker|BaseVaspMaker|None
+            Maker for optimization. Here: None.
+        phonon_maker_kwargs: dict
+            Additional keyword arguments for phonon maker.
+
+        Returns
+        -------
+        .BasePhononMaker
+
+        """
+
+    @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)
+        """

src/atomate2/common/jobs/qha.py

@@ -0,0 +1,122 @@
+"""Jobs for running qha calculations."""
+
+from __future__ import annotations
+
+import logging
+from typing import TYPE_CHECKING
+
+from jobflow import Flow, Response, job
+
+from atomate2.common.schemas.qha import PhononQHADoc
+
+if TYPE_CHECKING:
+    from pymatgen.core.structure import Structure
+
+    from atomate2.common.flows.phonons import BasePhononMaker
+    from atomate2.common.schemas.phonons import PhononBSDOSDoc
+
+logger = logging.getLogger(__name__)
+
+
+@job
+def get_phonon_jobs(phonon_maker: BasePhononMaker, eos_output: dict) -> Flow:
+    """
+    Start all relevant phonon jobs.
+
+    Parameters
+    ----------
+    phonon_maker: .BasePhononMaker
+        Maker to start harmonic phonon runs.
+    eos_output: dict
+        Output from EOSMaker
+
+    """
+    phonon_jobs = []
+    outputs = []
+    for structure in eos_output["relax"]["structure"]:
+        phonon_job = phonon_maker.make(structure)
+        phonon_jobs.append(phonon_job)
+        outputs.append(phonon_job.output)
+
+    return Response(replace=phonon_jobs, output=outputs)
+
+
+@job(
+    output_schema=PhononQHADoc,
+)
+def analyze_free_energy(
+    phonon_outputs: list[PhononBSDOSDoc],
+    structure: Structure,
+    t_max: float = None,
+    pressure: float = None,
+    ignore_imaginary_modes: bool = False,
+    eos_type: str = "vinet",
+    **kwargs,
+) -> Flow:
+    """Analyze the free energy from all phonon runs.
+
+    Parameters
+    ----------
+    phonon_outputs: list[PhononBSDOSDoc]
+        list of PhononBSDOSDoc objects
+    structure: Structure object
+        Corresponding structure object.
+    t_max: float
+        Max temperature for QHA in Kelvin.
+    pressure: float
+        Pressure for QHA in GPa.
+    ignore_imaginary_modes: bool
+        If True, all free energies will be used
+        for EOS fit
+    kwargs: dict
+        Additional keywords to pass to this job
+    """
+    # only add free energies if there are no imaginary modes
+    # tolerance has to be tested
+    electronic_energies: list[list[float]] = []
+    free_energies: list[list[float]] = []
+    heat_capacities: list[list[float]] = []
+    entropies: list[list[float]] = []
+    temperatures: list[float] = []
+    formula_units: list[int] = []
+    volume: list[float] = [
+        output.volume_per_formula_unit * output.formula_units
+        for output in phonon_outputs
+    ]
+
+    for itemp, temp in enumerate(phonon_outputs[0].temperatures):
+        temperatures.append(float(temp))
+        sorted_volume = []
+        electronic_energies.append([])
+        free_energies.append([])
+        heat_capacities.append([])
+        entropies.append([])
+
+        for _, output in sorted(zip(volume, phonon_outputs)):
+            # check if imaginary modes
+            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)
+                heat_capacities[itemp].append(output.heat_capacities[itemp])
+                entropies[itemp].append(output.entropies[itemp])
+                sorted_volume.append(output.volume_per_formula_unit)
+                formula_units.append(output.formula_units)
+
+    if len(set(formula_units)) != 1:
+        raise ValueError("There should be only one formula unit.")
+
+    return PhononQHADoc.from_phonon_runs(
+        volumes=sorted_volume,
+        free_energies=free_energies,
+        electronic_energies=electronic_energies,
+        entropies=entropies,
+        heat_capacities=heat_capacities,
+        temperatures=temperatures,
+        structure=structure,
+        t_max=t_max,
+        pressure=pressure,
+        formula_units=next(iter(set(formula_units))),
+        eos_type=eos_type,
+        **kwargs,
+    )