Workflow for Quasi-harmonic approximation (forcefields and VASP) (#903)

* add structures to eos output

* fix defaults in phonon schema

* add qha flows and one full workflow test

* really add qha flows

* fix description of qha schema file

* remove tests

* add schema test instead

* new test files

* remove future import

* Fix tests by changing the warning treatment back to default after the QHA api

* fix linting

* add switch to eos types

* change how warnings are treated

* add safer code for literal

* add more tests

* fix tests add kwargs

* Fix tests

* Fix tests

* reformat

* linting

* update prev calc

* Apply suggestions from code review

* fix test tolerance

* draft a first vasp workflow

* fix return

* fix return

* make analysis optional and and names for jobs

* add vasp test

* fix document

* fix dc

* update

* fix phonon tests

* fix linting

* fix schema

* change to phonon static

* change phonon maker

* fix tests

* fix test

* update qha_fresh

* fix comments

* fix comments

* connect flows

* fix linting

* correct documentation

* document workflows, fix imports

Author: JaGeo

docs/user/codes/vasp.md

@@ -260,6 +260,21 @@ With the help of phonopy, these forces are then converted into a dynamical matri
 The dynamical matrices of three structures are then used as an input to the phonopy Grueneisen api
 to compute mode-dependent Grueneisen parameters.
 
+
+### Quasi-harmonic Workflow
+Uses the quasi-harmonic approximation with the help of Phonopy to compute thermodynamic properties.
+First, a tight relaxation is performed. Subsequently, several optimizations at different constant
+volumes are performed. At each of the volumes, an additional phonon run is performed as well.
+Afterwards, equation of state fits are performed with phonopy.
+
+
+
+### Equation of State Workflow
+An equation of state workflow is implemented. First, a tight relaxation is performed. Subsequently, several optimizations at different constant
+volumes are performed. Additional static calculations might be performed afterwards to arrive at more
+accurate energies. Then, an equation of state fit is performed with pymatgen.
+
+
 ### LOBSTER
 
 Perform bonding analysis with [LOBSTER](http://cohp.de/) and [LobsterPy](https://github.com/jageo/lobsterpy)

src/atomate2/common/flows/eos.py

@@ -81,7 +81,10 @@ 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", "dir_name")
+            }
             for key in job_types
         }
 
@@ -163,9 +166,12 @@ def make(self, structure: Structure, prev_dir: str | Path = None) -> Flow:
         for key in job_types:
             for idx in range(len(jobs[key])):
                 output = jobs[key][idx].output.output
+                dir_name = jobs[key][idx].output.dir_name
                 flow_output[key]["energy"] += [output.energy]
                 flow_output[key]["volume"] += [output.structure.volume]
                 flow_output[key]["stress"] += [output.stress]
+                flow_output[key]["structure"] += [output.structure]
+                flow_output[key]["dir_name"] += [dir_name]
 
         if self.postprocessor is not None:
             min_points = self.postprocessor.min_data_points

src/atomate2/common/flows/qha.py

@@ -0,0 +1,178 @@
+"""Define common QHA flow agnostic to electronic-structure code."""
+
+from __future__ import annotations
+
+import warnings
+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
+    from atomate2.vasp.jobs.core 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_maker: .BasePhononMaker | None
+        Maker to compute phonons. The volume has to be fixed!
+        The beforehand relaxation could be switched off.
+    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)
+    skip_analysis: bool
+        Skips the analysis step and only performs EOS and phonon computations.
+    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_maker: BasePhononMaker | None = None
+    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
+    skip_analysis: 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'",
+            )
+
+        qha_jobs = []
+
+        # 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,
+        )
+
+        eos_job = self.eos.make(structure)
+        qha_jobs.append(eos_job)
+
+        phonon_jobs = get_phonon_jobs(
+            phonon_maker=self.phonon_maker, eos_output=eos_job.output
+        )
+        qha_jobs.append(phonon_jobs)
+        if not self.skip_analysis:
+            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,
+            )
+            qha_jobs.append(analysis)
+
+        return Flow(qha_jobs)
+
+    def __post_init__(self) -> None:
+        """Test settings during the initialisation."""
+        if self.phonon_maker.bulk_relax_maker is not None:
+            warnings.warn(
+                "An additional bulk_relax_maker has been added "
+                "to the phonon workflow. Please be aware "
+                "that the volume needs to be kept fixed.",
+                stacklevel=2,
+            )
+        # if self.phonon_maker.symprec != self.symprec:
+        #     warnings.warn(
+        #         "You are using different symmetry precisions "
+        #         "in the phonon makers and other parts of the "
+        #         "QHA workflow.",
+        #         stacklevel=2,
+        #     )
+        if self.phonon_maker.static_energy_maker is None:
+            warnings.warn(
+                "A static energy maker "
+                "is needed for "
+                "this workflow."
+                " Please add the static_energy_maker.",
+                stacklevel=2,
+            )
+
+    @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/eos.py

@@ -82,6 +82,8 @@ def fit(self, eos_flow_output: dict[str, Any]) -> None:
                         "energy": list, <required>
                         "volume": list, <required>
                         "stress": list <optional>
+                        "structure": list <not needed for the fit>
+                        "dir_name": list <optional for the fit>
                     },
                     "initial_<key>": {"E0": float, "V0": float} <optional>,
                         for <key> in ("relax", "static")

src/atomate2/common/jobs/qha.py

@@ -0,0 +1,133 @@
+"""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.phonons import PhononBSDOSDoc
+from atomate2.common.schemas.qha import PhononQHADoc
+
+if TYPE_CHECKING:
+    from pymatgen.core.structure import Structure
+
+    from atomate2.common.flows.phonons import BasePhononMaker
+
+logger = logging.getLogger(__name__)
+
+
+@job(
+    data=[PhononBSDOSDoc],
+)
+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 istructure, structure in enumerate(eos_output["relax"]["structure"]):
+        if eos_output["relax"]["dir_name"][istructure] is not None:
+            phonon_job = phonon_maker.make(
+                structure, prev_dir=eos_output["relax"]["dir_name"][istructure]
+            )
+        else:
+            phonon_job = phonon_maker.make(structure)
+        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)
+    return Response(replace=replace_flow)
+
+
+@job(
+    output_schema=PhononQHADoc,
+    data=["free_energies", "heat_capacities", "entropies", "helmholtz_volume"],
+)
+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)
+
+    # potentially implement a space group check in the future
+
+    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,
+    )