[WIP] Implementation of quasi-harmonic approximation

inset_write/Input_Dict

@@ -0,0 +1,27 @@
+$molecule
+ 0 1
+ O      0.0000000000      0.0000000000      0.1210000000
+ H     -0.7830000000      0.0000000000     -0.4850000000
+ H      0.7830000000      0.0000000000     -0.4850000000
+$end
+
+$rem
+   job_type = sp
+   basis = def2-tzvppd
+   max_scf_cycles = 100
+   gen_scfman = true
+   xc_grid = 3
+   thresh = 14
+   s2thresh = 16
+   scf_algorithm = diis
+   resp_charges = true
+   symmetry = false
+   sym_ignore = true
+   method = wb97mv
+$end
+
+$nbo
+$end
+
+$geom_opt
+$end

output

@@ -0,0 +1 @@
+No hook with id `../` in stage `pre-commit`

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/phonons.py

@@ -273,6 +273,7 @@ def make(
             )
             jobs.append(static_job)
             total_dft_energy = static_job.output.output.energy
+            stress = static_job.output.output.stress
             static_run_job_dir = static_job.output.dir_name
             static_run_uuid = static_job.output.uuid
             prev_dir = static_job.output.dir_name
@@ -283,6 +284,7 @@ def make(
             )
             jobs.append(compute_total_energy_job)
             total_dft_energy = compute_total_energy_job.output
+            stress = None
 
         # get a phonon object from phonopy
         displacements = generate_phonon_displacements(
@@ -347,6 +349,7 @@ def make(
             born_run_uuid=born_run_uuid,
             optimization_run_job_dir=optimization_run_job_dir,
             optimization_run_uuid=optimization_run_uuid,
+            stress=stress,
             create_thermal_displacements=self.create_thermal_displacements,
             store_force_constants=self.store_force_constants,
             **self.generate_frequencies_eigenvectors_kwargs,

src/atomate2/common/flows/qha.py

@@ -0,0 +1,158 @@
+"""Define common EOS 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
+
+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.base import BaseVaspMaker
+
+
+@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
+    initial_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker
+    eos_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker
+    phonon_displacement_maker: ForceFieldRelaxMaker | BaseVaspMaker
+    phonon_static_maker: ForceFieldRelaxMaker | BaseVaspMaker
+    phonon_maker_kwargs: dict
+    linear_strain: tuple[float, float]
+    number_of_frames: int
+    t_max: float | None
+    pressure: float | None
+    ignore_imaginary_modes: bool
+
+    name : str
+        Name of the flows produced by this maker.
+    initial_relax_maker : .Maker | None
+        Maker to relax the input structure, defaults to None (no initial relaxation).
+    eos_relax_maker : .Maker
+        Maker to relax deformed structures for the EOS fit.
+    phonon_static_maker : .Maker | None
+        Maker to generate statics after each relaxation, defaults to None.
+    strain : tuple[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.
+    #postprocessor : .atomate2.common.jobs.EOSPostProcessor
+    #    Optional postprocessing step, defaults to
+    #    `atomate2.common.jobs.PostProcessEosEnergy`.
+    #_store_transformation_information : .bool = False
+    #    Whether to store the information about transformations. Unfortunately
+    #    needed at present to handle issues with emmet and pydantic validation
+    #    TODO: remove this when clash is fixed
+    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
+    """
+
+    name: str = "QHA Maker"
+    initial_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None
+    eos_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None
+    phonon_displacement_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None
+    phonon_static_maker: ForceFieldRelaxMaker | 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
+    # 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
+        """
+        # 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)
+        # Todo: think about whether to keep the tight relax here
+        phonon_jobs = get_phonon_jobs(self.phonon_maker, eos_job.output)
+
+        # Todo: reuse postprocessor from equation of state to make fits of free energy curves
+        # get free energy fits and perform qha
+        analysis = analyze_free_energy(
+            phonon_jobs.output,
+            structure=structure,
+            t_max=self.t_max,
+            pressure=self.pressure,
+            ignore_imaginary_modes=self.ignore_imaginary_modes,
+        )
+
+        return Flow([eos_job, phonon_jobs, analysis])
+
+    @abstractmethod
+    def initialize_phonon_maker(
+        self,
+        phonon_displacement_maker,
+        phonon_static_maker,
+        bulk_relax_maker,
+        phonon_maker_kwargs,
+    ) -> BasePhononMaker | None:
+        """
+
+        :param phonon_displacement_maker:
+        :param phonon_static_maker:
+        :param bulk_relax_maker:
+        :param phonon_maker_kwargs:
+        :return:
+        """

src/atomate2/common/jobs/phonons.py

@@ -198,6 +198,7 @@ def generate_frequencies_eigenvectors(
     total_dft_energy: float,
     epsilon_static: Matrix3D = None,
     born: Matrix3D = None,
+    stress: Matrix3D = None,
     **kwargs,
 ) -> PhononBSDOSDoc:
     """
@@ -229,6 +230,8 @@ def generate_frequencies_eigenvectors(
         The high-frequency dielectric constant
     born: Matrix3D
         Born charges
+    stress: Optional[Matrix3D]
+        stress tensor
     kwargs: dict
         Additional parameters that are passed to PhononBSDOSDoc.from_forces_born
     """
@@ -247,6 +250,7 @@ def generate_frequencies_eigenvectors(
         total_dft_energy=total_dft_energy,
         epsilon_static=epsilon_static,
         born=born,
+        stress=stress,
         **kwargs,
     )