Add thermal expansion workflow (CTEMaker)

.github/workflows/testing.yml

@@ -70,6 +70,14 @@ jobs:
           cp -r tests/test_data/abinit/pseudos/ONCVPSP-PBE-SR-PDv0.4 ~/.abinit/pseudos
           uv pip install .[strict,strict-forcefields-${{ matrix.dep-group }},abinit,approxneb,aims] --group tests
 
+      # the thermal expansion tests use EMT, so one group is enough to run them
+      - name: Install pheasy, ALM and phono3py for the thermal expansion tests
+        if: matrix.dep-group == 'numpy-limited'
+        run: |
+          micromamba install -n a2 -c conda-forge compilers boost eigen=3.3 cmake spglib --yes
+          micromamba activate a2
+          uv pip install .[strict,strict-forcefields-${{ matrix.dep-group }},abinit,approxneb,aims,pheasy,phono3py] --group tests
+
       - name: Install pymatgen from master if triggered by pymatgen repo dispatch
         if: github.event_name == 'repository_dispatch' && github.event.action == 'pymatgen-ci-trigger'
         run: |
@@ -85,15 +93,15 @@ jobs:
         # xdist spawns multiple workers each holding a copy of the models.
         run: |
           micromamba activate a2
-          pytest --durations=5 -n ${{ matrix.dep-group == 'torch-limited' && '1' || 'auto' }} --cov=atomate2 --cov-report=xml tests/forcefields
+          pytest -v --durations=5 -n ${{ matrix.dep-group == 'torch-limited' && '1' || 'auto' }} --cov=atomate2 --cov-report=xml tests/forcefields
 
       - name: Forcefield tutorial
         if: matrix.dep-group == 'torch-limited'
         env:
           MP_API_KEY: ${{ secrets.MP_API_KEY }}
         run: |
           micromamba activate a2
-          pytest --durations=5 -n auto --nbmake ./tutorials/force_fields
+          pytest -v --durations=5 -n auto --nbmake ./tutorials/force_fields
 
       - uses: codecov/codecov-action@v1
         if: matrix.python-version == '3.11' && github.repository == 'materialsproject/atomate2'
@@ -147,7 +155,7 @@ jobs:
           python -m pip install --upgrade pip
           mkdir -p ~/.abinit/pseudos
           cp -r tests/test_data/abinit/pseudos/ONCVPSP-PBE-SR-PDv0.4 ~/.abinit/pseudos
-          uv pip install .[strict,abinit,approxneb,aims,pheasy,alamode,hiphive] --group tests
+          uv pip install .[strict,abinit,approxneb,aims,pheasy,alamode,hiphive,phono3py] --group tests
           uv pip install torch-runstats torch_dftd
           uv pip install --no-deps nequip==0.5.6
 
@@ -168,7 +176,7 @@ jobs:
         # However this `splitting-algorithm` means that tests cannot depend sensitively on the order they're executed in.
         run: |
           micromamba activate a2
-          pytest --durations=5 -n auto --splits 3 --group ${{ matrix.split }} --durations-path tests/.pytest-split-durations --splitting-algorithm least_duration --ignore=tests/ase --ignore=tests/openff_md --ignore=tests/openmm_md --ignore=tests/forcefields --ignore=tests/torchsim --cov=atomate2 --cov-report=xml
+          pytest -v --durations=5 -n auto --splits 3 --group ${{ matrix.split }} --durations-path tests/.pytest-split-durations --splitting-algorithm least_duration --ignore=tests/ase --ignore=tests/openff_md --ignore=tests/openmm_md --ignore=tests/forcefields --ignore=tests/torchsim --cov=atomate2 --cov-report=xml
 
 
       - uses: codecov/codecov-action@v1
@@ -295,7 +303,7 @@ jobs:
           MP_API_KEY: ${{ secrets.MP_API_KEY }}
         run: |
           micromamba activate a2
-          pytest --durations=5 -n auto --nbmake ./tutorials --ignore=./tutorials/openmm_tutorial.ipynb --ignore=./tutorials/force_fields --ignore=./tutorials/torchsim_tutorial.ipynb --ignore=./tutorials/lammps_workflow.ipynb --ignore=./tutorials/pheasy_workflow.ipynb --ignore=./tutorials/hiphive_workflow.ipynb 
+          pytest --durations=5 -n auto --nbmake ./tutorials --ignore=./tutorials/openmm_tutorial.ipynb --ignore=./tutorials/force_fields --ignore=./tutorials/torchsim_tutorial.ipynb --ignore=./tutorials/lammps_workflow.ipynb --ignore=./tutorials/pheasy_workflow.ipynb --ignore=./tutorials/hiphive_workflow.ipynb --ignore=./tutorials/cte_workflow.ipynb 
 
       - name: Test ASE
         env:

docs/user/codes/vasp.md

@@ -437,6 +437,68 @@ gruneisen_flow = GruneisenMaker(
 ).make(structure=structure)
 ```
 
+### Thermal expansion workflow
+
+`CTEMaker` calculates the thermal expansion tensor from third-order force constants, with the help of [Pheasy](https://doi.org/10.48550/arXiv.2508.01020) and [phono3py](https://doi.org/10.1088/1361-648X/acd831).
+It needs the `pheasy` and `phono3py` extras, `pip install "atomate2[pheasy,phono3py]"`.
+For the pheasy extra, see the Pheasy section above.
+
+```{warning}
+This workflow is new and has not been tested widely.
+It might still change in future versions.
+```
+
+First, the structure is converted to the standard primitive cell, and a tight structural relaxation is performed.
+The pheasy fits can fail for cells that are not in a standard setting.
+Set `use_symmetrized_structure="conventional"` to use the standard conventional cell instead.
+The relaxed structure is then passed to the pheasy phonon workflow and to the elastic constant workflow.
+The two do not depend on each other, so a workflow manager can run them at the same time.
+Neither of them relaxes the structure again, so both use the same structure.
+The phonon workflow fits the third-order force constants with LASSO, on randomly displaced supercells with 0.03 Å displacements.
+The phonon maker builds its supercells with `min_length=12.0`.
+By default, it uses the one-shot fit, which fits the second-order force constants together with them.
+The cocktail fit keeps the second-order force constants of the harmonic fit.
+The thermal expansion of each fit uses the second-order force constants of that fit.
+phono3py then gives the mode Grüneisen tensors on a 12x12x12 q-point mesh.
+The thermal expansion tensor follows from the mode heat capacities, the Grüneisen tensors and the elastic compliance.
+It is computed for each fit in `anhar_fit_methods` of the phonon maker, from 0 K to 1000 K in steps of 10 K by default.
+If a frequency on the mesh is below -0.1 THz, a warning is raised and the thermal expansion of that fit is not computed.
+
+The mode Grüneisen tensors come from the third-order force constants at the relaxed structure.
+The phonon frequencies are not renormalized with temperature.
+The workflow uses PBEsol by default.
+The stress is more sensitive to ENCUT than the forces are, so check the ENCUT convergence of the elastic tensor for your material.
+For metals, set `born_maker=None` in the phonon maker to skip the Born charge calculation.
+The `compute_cte` job reads the force constant files from the folder of the pheasy fit, so it must run where that folder can be read.
+
+A thermal expansion workflow for VASP can be started as follows:
+```python
+from atomate2.vasp.flows.cte import CTEMaker
+from pymatgen.core.structure import Structure
+
+structure = Structure(
+    lattice=[[0, 2.13, 2.13], [2.13, 0, 2.13], [2.13, 2.13, 0]],
+    species=["Mg", "O"],
+    coords=[[0, 0, 0], [0.5, 0.5, 0.5]],
+)
+
+cte_flow = CTEMaker().make(structure=structure)
+```
+
+`update_user_incar_settings` changes the INCAR of every VASP job in the flow.
+For example, this switches all of them to r2SCAN:
+```python
+from atomate2.vasp.powerups import update_user_incar_settings
+
+cte_flow = update_user_incar_settings(cte_flow, {"GGA": None, "METAGGA": "R2SCAN"})
+```
+The Born charge job uses DFPT with `LEPSILON = True`.
+Check that your VASP version runs DFPT with r2SCAN before you use this setting.
+
+The same workflow runs with a force field via `from atomate2.forcefields.flows.cte import CTEMaker`.
+`CTEMaker.from_force_field_name` sets one force field for the relaxation, the phonons and the elastic tensor.
+The notebook `tutorials/cte_workflow.ipynb` runs it for MgO with MACE-OMAT-0-medium.
+
 ### Quasi-harmonic Workflow
 
 Uses the quasi-harmonic approximation with the help of Phonopy to compute thermodynamic properties.

pyproject.toml

@@ -55,6 +55,9 @@ phonons = ["phonopy>=2.43,<5", "seekpath>=2.0.0"]
 # the pheasy fork adds the --disp_matrix_file and --force_matrix_file options
 # and fixes the --scell option
 pheasy = ["atomate2[phonons,alamode]", "hiphive==1.5", "numpy<=2.2", "pheasy @ git+https://gitlab.com/hpsahasrabuddhe/pheasy.git@9f24162a4ed0f0ab8911d382fd2617944aade55d"]
+# phono3py computes the mode Grueneisen tensors in the thermal expansion workflow
+# phonors 0.5 breaks the Grueneisen q-point mesh of phonopy and phono3py 4.5
+phono3py = ["atomate2[phonons]", "phono3py==4.5.0", "phonors<0.5"]
 alamode = ["alm @ git+https://github.com/ttadano/ALM.git@f1d668fdee66e7e7218a04c88daf19d0e14fce0c#subdirectory=python"]
 hiphive = ["hiphive==1.5", "trainstation>=1.0", "atomate2[phonons,alamode]"]
 lobster = ["ijson>=3.2.2", "lobsterpy>=0.6.0"]

src/atomate2/common/flows/cte.py

@@ -0,0 +1,178 @@
+"""Flow for the thermal expansion from third-order force constants."""
+
+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 pymatgen.util.due import Doi, due
+
+from atomate2.common.jobs.cte import compute_cte
+from atomate2.common.jobs.utils import structure_to_conventional, structure_to_primitive
+
+if TYPE_CHECKING:
+    from pathlib import Path
+
+    from pymatgen.core.structure import Structure
+
+    from atomate2.common.flows.elastic import BaseElasticMaker
+    from atomate2.common.flows.pheasy import BasePhononMaker
+    from atomate2.forcefields.jobs import ForceFieldRelaxMaker
+    from atomate2.vasp.jobs.base import BaseVaspMaker
+
+
+@due.dcite(
+    Doi("10.1088/1361-648X/acd831"),
+    description="Implementation strategies in phonopy and phono3py.",
+)
+@due.dcite(
+    Doi("10.7566/JPSJ.92.012001"),
+    description="Phonopy and phono3py.",
+)
+@dataclass
+class BaseCTEMaker(Maker, ABC):
+    """
+    Maker to calculate the thermal expansion from third-order force constants.
+
+    By default, the structure is first converted to the standard primitive
+    cell, and a tight structural relaxation follows. The relaxed structure is
+    then passed to the pheasy phonon flow, which fits the second- and
+    third-order force constants, and to the elastic flow. The two flows do not
+    depend on each other, so a workflow manager can run them at the same time.
+    Neither flow may relax the structure again, so that both use the same
+    structure in the same frame. Finally, phono3py gives the mode
+    Grueneisen tensors on a q-point mesh, and the thermal expansion tensor
+    follows from the heat-capacity weighted Grueneisen tensors and the elastic
+    compliance. The mode Grueneisen tensors come from the third-order force
+    constants at the relaxed structure. The frequencies are not renormalized
+    with temperature.
+
+    This workflow is new and has not been tested widely. It might still change
+    in future versions.
+
+    Parameters
+    ----------
+    name: str
+        Name of the flows produced by this maker.
+    use_symmetrized_structure: str or None
+        Convert the input structure to the standard "primitive" or
+        "conventional" cell before the relaxation. The phonon and elastic flows
+        both use the converted structure. The pheasy fits can fail for cells
+        that are not in a standard setting, so only set this to None for an
+        input structure that already is.
+    bulk_relax_maker: .ForceFieldRelaxMaker, .BaseVaspMaker, or None
+        A maker to perform a tight relaxation on the bulk. Set to None to skip
+        the relaxation.
+    phonon_maker: .BasePhononMaker
+        The pheasy phonon maker. It must have cal_anhar_fcs=True,
+        use_symmetrized_structure=None and bulk_relax_maker=None. Its
+        anhar_fit_methods set which force constants are used for the thermal
+        expansion.
+    elastic_maker: .BaseElasticMaker
+        Maker for the elastic tensor. It must have bulk_relax_maker=None.
+    temperatures: list[float]
+        Temperatures in K.
+    mesh: tuple[int, int, int] | float
+        q-point mesh for the mode Grueneisen tensors, or a q-point density used
+        as kppa in pymatgen's Kpoints.automatic_density for the unit cell.
+    tol_imaginary_modes: float
+        If a frequency on the mesh is below -tol_imaginary_modes in THz, a
+        warning is raised and the thermal expansion of that fit is not computed.
+    min_frequency: float
+        Modes below this frequency in THz are left out of the thermal expansion.
+    """
+
+    name: str = "cte"
+    use_symmetrized_structure: Literal["primitive", "conventional"] | None = "primitive"
+    bulk_relax_maker: ForceFieldRelaxMaker | BaseVaspMaker | None = None
+    phonon_maker: BasePhononMaker = None
+    elastic_maker: BaseElasticMaker = None
+    temperatures: list[float] = field(default_factory=lambda: list(range(0, 1001, 10)))
+    mesh: tuple[int, int, int] | float = (12, 12, 12)
+    tol_imaginary_modes: float = 0.1
+    min_frequency: float = 1e-3
+
+    def __post_init__(self) -> None:
+        """Check that the phonon and elastic makers fit this workflow."""
+        if not self.phonon_maker.cal_anhar_fcs:
+            raise ValueError(
+                "The phonon maker needs cal_anhar_fcs=True, since the thermal "
+                "expansion needs the third-order force constants."
+            )
+        if self.phonon_maker.use_symmetrized_structure is not None:
+            raise ValueError(
+                "The phonon maker needs use_symmetrized_structure=None, so that the "
+                "phonon and elastic calculations use the same frame."
+            )
+        for label, maker in (
+            ("phonon", self.phonon_maker),
+            ("elastic", self.elastic_maker),
+        ):
+            if maker.bulk_relax_maker is not None:
+                raise ValueError(
+                    f"The {label} maker needs bulk_relax_maker=None. Otherwise the "
+                    "structure is relaxed again, and the phonon and elastic "
+                    "calculations do not use the same structure."
+                )
+
+    def make(self, structure: Structure, prev_dir: str | Path | None = None) -> Flow:
+        """
+        Make a flow to calculate the thermal expansion.
+
+        Parameters
+        ----------
+        structure: Structure
+            A pymatgen structure. Start with a structure that is nearly fully
+            optimized, as the relaxation settings are strict.
+        prev_dir: str or Path or None
+            A previous calculation directory to use for copying outputs.
+        """
+        jobs = []
+        if self.use_symmetrized_structure == "primitive":
+            prim_job = structure_to_primitive(structure, self.phonon_maker.symprec)
+            jobs.append(prim_job)
+            structure = prim_job.output
+        elif self.use_symmetrized_structure == "conventional":
+            conv_job = structure_to_conventional(structure, self.phonon_maker.symprec)
+            jobs.append(conv_job)
+            structure = conv_job.output
+
+        equilibrium_stress = None
+        if self.bulk_relax_maker is not None:
+            bulk_kwargs = {}
+            if self.prev_calc_dir_argname is not None:
+                bulk_kwargs[self.prev_calc_dir_argname] = prev_dir
+            bulk = self.bulk_relax_maker.make(structure, **bulk_kwargs)
+            jobs.append(bulk)
+            structure = bulk.output.structure
+            prev_dir = bulk.output.dir_name
+            # as in the elastic flow when it runs its own relaxation
+            equilibrium_stress = bulk.output.output.stress
+
+        phonon_flow = self.phonon_maker.make(structure, prev_dir=prev_dir)
+        elastic_flow = self.elastic_maker.make(
+            structure, prev_dir=prev_dir, equilibrium_stress=equilibrium_stress
+        )
+        cte = compute_cte(
+            phonon_output=phonon_flow.output,
+            elastic_tensor=elastic_flow.output.elastic_tensor.raw,
+            elastic_structure=elastic_flow.output.structure,
+            anhar_fit_methods=self.phonon_maker.anhar_fit_methods,
+            temperatures=self.temperatures,
+            mesh=self.mesh,
+            tol_imaginary_modes=self.tol_imaginary_modes,
+            min_frequency=self.min_frequency,
+            symprec=self.phonon_maker.symprec,
+        )
+        jobs += [phonon_flow, elastic_flow, cte]
+        return Flow(jobs, output=cte.output, name=self.name)
+
+    @property
+    @abstractmethod
+    def prev_calc_dir_argname(self) -> str | None:
+        """Name of the argument that passes the previous calculation directory.
+
+        It differs between codes, so each subclass sets it.
+        """