diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index 07e81ea6b8..7405467740 100644 --- a/.github/workflows/testing.yml +++ b/.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,7 +93,7 @@ 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' @@ -93,7 +101,7 @@ jobs: 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: diff --git a/docs/user/codes/vasp.md b/docs/user/codes/vasp.md index 60b8f27c3f..f283cd436d 100644 --- a/docs/user/codes/vasp.md +++ b/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. diff --git a/pyproject.toml b/pyproject.toml index 8cfc833508..8d9b27680c 100644 --- a/pyproject.toml +++ b/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"] diff --git a/src/atomate2/common/flows/cte.py b/src/atomate2/common/flows/cte.py new file mode 100644 index 0000000000..9868a76f33 --- /dev/null +++ b/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. + """ diff --git a/src/atomate2/common/jobs/cte.py b/src/atomate2/common/jobs/cte.py new file mode 100644 index 0000000000..982d398513 --- /dev/null +++ b/src/atomate2/common/jobs/cte.py @@ -0,0 +1,102 @@ +"""Jobs for thermal expansion from mode Grueneisen tensors.""" + +from __future__ import annotations + +from pathlib import Path +from typing import TYPE_CHECKING + +from jobflow import job + +from atomate2.common.jobs.gruneisen import PhononDoc, _get_taskdoc_run_dir +from atomate2.common.jobs.pheasy import _ANHARMONIC_FIT_METHODS, _DEFAULT_FILE_PATHS +from atomate2.common.schemas.cte import CTEDocument + +if TYPE_CHECKING: + from collections.abc import Sequence + + from emmet.core.math import MatrixVoigt + from pymatgen.core import Structure + + +@job(output_schema=CTEDocument) +def compute_cte( + phonon_output: PhononDoc, + elastic_tensor: MatrixVoigt, + elastic_structure: Structure, + anhar_fit_methods: Sequence[str] = ("one-shot",), + temperatures: Sequence[float] = tuple(range(0, 1001, 10)), + mesh: tuple[int, int, int] | float = (12, 12, 12), + tol_imaginary_modes: float = 0.1, + min_frequency: float = 1e-3, + symprec: float = 1e-5, +) -> CTEDocument: + """ + Compute the thermal expansion from the pheasy force constants. + + The second- and third-order force constants of each fit method are read + from the folder of the pheasy fit, so this job must run where that folder + can be read. The cocktail fit uses the second-order force constants of the + harmonic fit. The one-shot fit uses its own. + + Parameters + ---------- + phonon_output: PhononDoc + Output document of the pheasy phonon flow, run with cal_anhar_fcs=True. + elastic_tensor: MatrixVoigt + Elastic tensor in GPa and Voigt notation, as ElasticDocument's + elastic_tensor.raw. + elastic_structure: Structure + Structure of the elastic calculation. Its lattice must match the unit + cell of the phonon run, so that both tensors are in the same frame. + anhar_fit_methods: Sequence[str] + Fit methods whose force constants are used, "cocktail" and/or "one-shot". + temperatures: Sequence[float] + Temperatures in K, not negative. + mesh: tuple[int, int, int] | float + q-point mesh, 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. + symprec: float + Symmetry precision passed to phono3py. + + Returns + ------- + CTEDocument + """ + unknown = set(anhar_fit_methods) - set(_ANHARMONIC_FIT_METHODS) + if unknown or not anhar_fit_methods: + raise ValueError( + f"anhar_fit_methods must be a non-empty subset of " + f"{_ANHARMONIC_FIT_METHODS}, not {list(anhar_fit_methods)}." + ) + job_dir_name = _get_taskdoc_run_dir(phonon_output) + if job_dir_name is None: + raise ValueError("The phonon output does not record the pheasy job folder.") + job_dir = Path(job_dir_name) + + # the files written by the pheasy fits + one_shot_dir = job_dir / _DEFAULT_FILE_PATHS["one_shot_dir"] + files = { + "cocktail": ( + job_dir / _DEFAULT_FILE_PATHS["force_constants"], + job_dir / "fc3.hdf5", + ), + "one-shot": (one_shot_dir / "fc2.hdf5", one_shot_dir / "fc3.hdf5"), + } + force_constant_files = {method: files[method] for method in anhar_fit_methods} + + return CTEDocument.from_force_constants( + phonopy_yaml=job_dir / _DEFAULT_FILE_PATHS["phonopy"], + force_constant_files=force_constant_files, + elastic_tensor=elastic_tensor, + structure=elastic_structure, + temperatures=temperatures, + mesh=mesh, + tol_imaginary_modes=tol_imaginary_modes, + min_frequency=min_frequency, + symprec=symprec, + ) diff --git a/src/atomate2/common/schemas/cte.py b/src/atomate2/common/schemas/cte.py new file mode 100644 index 0000000000..dcb0a1c080 --- /dev/null +++ b/src/atomate2/common/schemas/cte.py @@ -0,0 +1,372 @@ +"""Schemas for the thermal expansion workflow outputs.""" + +from __future__ import annotations + +import warnings +from pathlib import Path +from typing import TYPE_CHECKING + +import numpy as np +from emmet.core.math import Matrix3D, MatrixVoigt +from emmet.core.structure import StructureMetadata +from pydantic import BaseModel, Field +from pymatgen.core import Structure +from pymatgen.io.phonopy import get_pmg_structure +from pymatgen.io.vasp import Kpoints +from scipy.constants import Boltzmann, Planck + +if TYPE_CHECKING: + from collections.abc import Mapping, Sequence + + from phonopy import Phonopy + from typing_extensions import Self + +# Voigt order of the symmetric 3x3 tensor components, as in pymatgen +_VOIGT_INDICES = ((0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)) + + +def get_cte( + frequencies: np.ndarray, + gruneisen_tensors: np.ndarray, + weights: np.ndarray, + elastic_tensor: np.ndarray, + volume: float, + temperatures: Sequence[float], + min_frequency: float = 1e-3, +) -> tuple[np.ndarray, list[np.ndarray | None]]: + """ + Get the thermal expansion tensor from mode Grueneisen tensors. + + The thermal stress of each mode is its heat capacity times its Grueneisen + tensor. The strain that relaxes the summed thermal stress is the thermal + expansion, alpha = S sum(c * gamma) / (N_q * V). Here S is the elastic + compliance, c are the modal heat capacities, N_q is the sum of the q-point + weights and V is the cell volume. The mode Grueneisen tensors are + symmetrized first, since the strain is symmetric. Modes below + min_frequency, including the acoustic modes at Gamma, are left out. + + Parameters + ---------- + frequencies: np.ndarray + Phonon frequencies in THz, with shape (n_qpoints, n_bands). + gruneisen_tensors: np.ndarray + Mode Grueneisen tensors, with shape (n_qpoints, n_bands, 3, 3). + weights: np.ndarray + Weight of each q-point, with shape (n_qpoints,). + elastic_tensor: np.ndarray + Elastic tensor in GPa and Voigt notation, in the same Cartesian frame as + the Grueneisen tensors. + volume: float + Volume of the cell used for the phonons, in Angstrom^3. + temperatures: Sequence[float] + Temperatures in K. + min_frequency: float + Modes below this frequency in THz are left out. + + Returns + ------- + tuple[np.ndarray, list[np.ndarray | None]] + The thermal expansion tensors in 1/K, with shape (n_temperatures, 3, 3), + and the heat-capacity weighted mean Grueneisen tensor at each + temperature. The mean is None where the heat capacity is zero. + """ + frequencies = np.asarray(frequencies, dtype=float) + gruneisen_tensors = np.asarray(gruneisen_tensors, dtype=float) + weights = np.asarray(weights, dtype=float) + + # S in 1/Pa and V in m^3 + compliance = np.linalg.inv(np.asarray(elastic_tensor, dtype=float) * 1e9) + volume_m3 = volume * 1e-30 + + # the Grueneisen tensors of the left-out modes can be large or undefined, + # so zero them + kept = frequencies > min_frequency + gruneisen_tensors = np.where(kept[..., None, None], gruneisen_tensors, 0.0) + gruneisen_tensors = (gruneisen_tensors + np.swapaxes(gruneisen_tensors, -1, -2)) / 2 + gruneisen_voigt = np.stack( + [gruneisen_tensors[..., i, j] for i, j in _VOIGT_INDICES], axis=-1 + ) + energies = Planck * frequencies * 1e12 # J + + alphas, mean_gruneisen = [], [] + for temperature in temperatures: + if temperature <= 0: + heat_capacities = np.zeros_like(frequencies) + else: + x = np.where(kept, energies / (Boltzmann * temperature), 1.0) + # x^2 e^x / (e^x - 1)^2, written with e^-x so that it cannot overflow + heat_capacities = np.where( + kept, Boltzmann * x**2 * np.exp(-x) / np.expm1(-x) ** 2, 0.0 + ) # J/K + weighted = heat_capacities * weights[:, None] + thermal_stress = np.einsum("qb,qbv->v", weighted, gruneisen_voigt) + thermal_stress /= weights.sum() * volume_m3 # Pa/K + + # the compliance gives engineering shear strains, twice the tensor ones + alpha_voigt = compliance @ thermal_stress + alpha = np.empty((3, 3)) + for k, (i, j) in enumerate(_VOIGT_INDICES): + alpha[i, j] = alpha[j, i] = alpha_voigt[k] if k < 3 else alpha_voigt[k] / 2 + alphas.append(alpha) + + total_heat_capacity = weighted.sum() + if total_heat_capacity > 0: + mean = np.einsum("qb,qbij->ij", weighted, gruneisen_tensors) + mean_gruneisen.append(mean / total_heat_capacity) + else: + mean_gruneisen.append(None) + + return np.array(alphas), mean_gruneisen + + +def _expand_born_to_unitcell(phonon: Phonopy) -> np.ndarray: + """Map the Born charges of the phonopy primitive cell onto the unit cell atoms.""" + primitive = phonon.primitive + born = np.asarray(phonon.nac_params["born"]) + if len(primitive) == len(phonon.unitcell): + return born + primitive_indices = [ + primitive.p2p_map[primitive.s2p_map[s]] for s in phonon.supercell.u2s_map + ] + return born[primitive_indices] + + +class CTEResult(BaseModel): + """Thermal expansion from one set of second- and third-order force constants.""" + + fit_method: str = Field( + description='Anharmonic fit that gave the force constants, "cocktail" or ' + '"one-shot".' + ) + lowest_frequency: float = Field( + description="Lowest phonon frequency on the sampling mesh in THz. Imaginary " + "frequencies are negative." + ) + has_imaginary_modes: bool = Field( + description="Whether a frequency on the sampling mesh lies below " + "-tol_imaginary_modes. The thermal expansion is not computed in that case." + ) + thermal_expansion_tensor: list[Matrix3D] | None = Field( + None, + description="Thermal expansion tensor in 1/K at each temperature, in the " + "Cartesian frame of the structure.", + ) + thermal_expansion: list[float] | None = Field( + None, + description="Volumetric thermal expansion in 1/K at each temperature, the " + "trace of the thermal expansion tensor.", + ) + average_gruneisen: list[Matrix3D | None] | None = Field( + None, + description="Mode Grueneisen tensor averaged with the mode heat capacities " + "at each temperature. None where the heat capacity is zero.", + ) + + +class CTEDocument(StructureMetadata): + """Thermal expansion from mode Grueneisen tensors and the elastic tensor.""" + + structure: Structure | None = Field( + None, description="Structure used for the phonon and elastic calculations." + ) + supercell_matrix: Matrix3D | None = Field( + None, description="Supercell matrix of the phonon calculation." + ) + temperatures: list[float] | None = Field(None, description="Temperatures in K.") + mesh: tuple[int, int, int] | None = Field( + None, description="q-point mesh used for the mode Grueneisen tensors." + ) + elastic_tensor: MatrixVoigt | None = Field( + None, + description="Elastic tensor in GPa and Voigt notation, in the Cartesian " + "frame of the structure.", + ) + min_frequency: float | None = Field( + None, + description="Modes below this frequency in THz are left out of the thermal " + "expansion.", + ) + tol_imaginary_modes: float | None = Field( + None, + description="The thermal expansion of a fit is not computed if a frequency " + "on the mesh is below -tol_imaginary_modes in THz.", + ) + phonon_job_dir: str | None = Field( + None, description="Directory of the pheasy fit that wrote the force constants." + ) + results: list[CTEResult] | None = Field( + None, description="Thermal expansion for each anharmonic fit method." + ) + + @classmethod + def from_force_constants( + cls, + phonopy_yaml: str | Path, + force_constant_files: Mapping[str, tuple[str | Path, str | Path]], + elastic_tensor: MatrixVoigt, + structure: Structure, + temperatures: Sequence[float], + mesh: tuple[int, int, int] | float, + tol_imaginary_modes: float, + min_frequency: float, + symprec: float, + ) -> Self: + """ + Compute the thermal expansion from second- and third-order force constants. + + phono3py gives the frequencies and mode Grueneisen tensors on the q-point + mesh. Only time reversal symmetry is used to reduce the mesh. The + non-analytical term correction is applied to the dynamical matrix when + the phonopy yaml file holds the Born charges and the dielectric tensor. + Only the third-order force constants enter the strain derivative of the + dynamical matrix. The results of each fit are written to + gruneisen_.hdf5 in the current directory. + + Parameters + ---------- + phonopy_yaml: str or Path + phonopy.yaml of the pheasy fit, with the unit cell, the supercell + matrix and the Born charges. + force_constant_files: Mapping + For each fit method, "cocktail" or "one-shot", the files with the + second- and third-order force constants. + elastic_tensor: MatrixVoigt + Elastic tensor in GPa and Voigt notation. + structure: Structure + Structure of the elastic calculation. Its lattice must match the + unit cell in phonopy_yaml, so that both tensors are in the same frame. + temperatures: Sequence[float] + Temperatures in K, not negative. + mesh: tuple[int, int, int] or float + q-point mesh, 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. + symprec: float + Symmetry precision passed to phono3py. + + Returns + ------- + CTEDocument + """ + import h5py + import phonopy + from phono3py import Phono3py + from phono3py.file_IO import read_fc2_from_hdf5, read_fc3_from_hdf5 + from phono3py.phonon3.gruneisen import Gruneisen + from phonopy.file_IO import parse_FORCE_CONSTANTS + + if not force_constant_files: + raise ValueError("No force constant files were given.") + if min(temperatures) < 0: + raise ValueError("The temperatures must not be negative.") + + phonon = phonopy.load(phonopy_yaml, produce_fc=False, log_level=0) + if not np.allclose(phonon.unitcell.cell, structure.lattice.matrix, atol=1e-5): + raise ValueError( + "The lattice of the elastic calculation differs from the unit cell " + "of the phonon calculation, so the two tensors are not in the same " + "frame." + ) + if list(phonon.unitcell.symbols) != [site.specie.symbol for site in structure]: + raise ValueError( + "The atoms of the elastic calculation differ from the unit cell of " + "the phonon calculation." + ) + + # pheasy writes compact force constants for the unit cell in POSCAR, so + # the unit cell is also the primitive cell here. "P" keeps it as it is. + ph3 = Phono3py( + phonon.unitcell, + supercell_matrix=phonon.supercell_matrix, + primitive_matrix="P", + symprec=symprec, + ) + nac_params = None + if phonon.nac_params is not None: + nac_params = { + **phonon.nac_params, + "born": _expand_born_to_unitcell(phonon), + } + + if isinstance(mesh, int | float | np.number): + kpoints = Kpoints.automatic_density( + structure=get_pmg_structure(ph3.primitive), + kppa=float(mesh), + force_gamma=True, + ) + mesh_numbers = tuple(int(m) for m in kpoints.kpts[0]) + else: + mesh_numbers = tuple(int(m) for m in mesh) + + results = [] + for method, (fc2_file, fc3_file) in force_constant_files.items(): + if Path(fc2_file).suffix == ".hdf5": + fc2 = read_fc2_from_hdf5(fc2_file) + else: + fc2 = parse_FORCE_CONSTANTS(filename=fc2_file) + fc3 = read_fc3_from_hdf5(fc3_file) + + gruneisen = Gruneisen( + fc2, fc3, ph3.supercell, ph3.primitive, nac_params=nac_params + ) + gruneisen.set_sampling_mesh(mesh_numbers, primitive_symmetry=None) + gruneisen.run() + filename = f"gruneisen_{method.replace('-', '_')}" + gruneisen.write(filename=filename) + # the full third-order force constants can take several GB + del gruneisen, fc3 + with h5py.File(f"{filename}.hdf5") as file: + frequencies = file["frequency"][:] + gruneisen_tensors = file["gruneisen_tensor"][:] + weights = file["weight"][:] + + lowest = float(frequencies.min()) + has_imaginary_modes = lowest < -tol_imaginary_modes + result = { + "fit_method": method, + "lowest_frequency": lowest, + "has_imaginary_modes": has_imaginary_modes, + } + if has_imaginary_modes: + warnings.warn( + f"The {method} force constants give a frequency of " + f"{lowest:.3f} THz, below -{tol_imaginary_modes} THz. The " + "thermal expansion is not computed for this fit.", + stacklevel=2, + ) + else: + alphas, mean_gruneisen = get_cte( + frequencies, + gruneisen_tensors, + weights, + np.asarray(elastic_tensor), + ph3.primitive.volume, + temperatures, + min_frequency=min_frequency, + ) + result["thermal_expansion_tensor"] = alphas.tolist() + result["thermal_expansion"] = np.trace( + alphas, axis1=1, axis2=2 + ).tolist() + result["average_gruneisen"] = [ + None if mean is None else mean.tolist() for mean in mean_gruneisen + ] + results.append(CTEResult(**result)) + + return cls.from_structure( + meta_structure=structure, + structure=structure, + supercell_matrix=phonon.supercell_matrix.tolist(), + temperatures=list(temperatures), + mesh=mesh_numbers, + elastic_tensor=np.asarray(elastic_tensor).tolist(), + min_frequency=min_frequency, + tol_imaginary_modes=tol_imaginary_modes, + phonon_job_dir=str(Path(phonopy_yaml).parent), + results=results, + ) diff --git a/src/atomate2/forcefields/flows/cte.py b/src/atomate2/forcefields/flows/cte.py new file mode 100644 index 0000000000..d083ddf463 --- /dev/null +++ b/src/atomate2/forcefields/flows/cte.py @@ -0,0 +1,149 @@ +"""Define the force field thermal expansion maker.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any + +from atomate2.common.flows.cte import BaseCTEMaker +from atomate2.forcefields.flows.elastic import ElasticMaker +from atomate2.forcefields.flows.pheasy import PhononMaker +from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker + +if TYPE_CHECKING: + from typing_extensions import Self + + from atomate2.forcefields import MLFF + +_DEFAULT_FORCE_FIELD = "MACE-MP-0" + + +def _get_makers( + force_field_name: str | MLFF | dict, calculator_kwargs: dict | None = None +) -> dict: + """Get the relaxation, phonon and elastic makers for one force field.""" + calculator: dict[str, Any] = { + "force_field_name": force_field_name, + "calculator_kwargs": dict(calculator_kwargs or {}), + } + # the relaxation settings are those of the force field elastic flow + return { + "bulk_relax_maker": ForceFieldRelaxMaker( + relax_cell=True, + relax_kwargs={"fmax": 0.00001}, + fix_symmetry=True, + **calculator, + ), + "phonon_maker": PhononMaker( + min_length=12.0, + bulk_relax_maker=None, + static_energy_maker=None, + phonon_displacement_maker=ForceFieldStaticMaker(**calculator), + cal_anhar_fcs=True, + displacement_anhar=0.03, + anhar_fit_methods=("one-shot",), + ), + "elastic_maker": ElasticMaker( + bulk_relax_maker=None, + elastic_relax_maker=ForceFieldRelaxMaker( + relax_cell=False, + relax_kwargs={"fmax": 0.00001}, + fix_symmetry=True, + **calculator, + ), + ), + } + + +@dataclass +class CTEMaker(BaseCTEMaker): + """ + Maker to calculate the thermal expansion with a force field and pheasy. + + By default, the structure is converted to the standard primitive cell and + relaxed tightly. 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. Neither flow relaxes the structure + again. Finally, phono3py gives the mode Grueneisen tensors, and the thermal + expansion tensor follows from them and the elastic tensor. The frequencies + are not renormalized with temperature. + + By default, all steps use MACE-MP-0. Use :obj:`from_force_field_name` to + run every step with another force field. The phonon flow fits the force + constants with the one-shot method from randomly displaced supercells with + 0.03 A displacements, and builds the supercells with min_length=12.0. + + 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 pheasy fits can fail + for cells that are not in a standard setting. + bulk_relax_maker: .ForceFieldRelaxMaker or None + A maker to perform a tight relaxation on the bulk. Set to None to skip + the relaxation. + phonon_maker: .PhononMaker + 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: .ElasticMaker + 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" + bulk_relax_maker: ForceFieldRelaxMaker | None = field( + default_factory=lambda: _get_makers(_DEFAULT_FORCE_FIELD)["bulk_relax_maker"] + ) + phonon_maker: PhononMaker = field( + default_factory=lambda: _get_makers(_DEFAULT_FORCE_FIELD)["phonon_maker"] + ) + elastic_maker: ElasticMaker = field( + default_factory=lambda: _get_makers(_DEFAULT_FORCE_FIELD)["elastic_maker"] + ) + + @classmethod + def from_force_field_name( + cls, + force_field_name: str | MLFF | dict, + calculator_kwargs: dict | None = None, + **kwargs, + ) -> Self: + """ + Create a thermal expansion maker that uses one force field for all steps. + + Parameters + ---------- + force_field_name: str or .MLFF or dict + The name of the force field. + calculator_kwargs: dict or None + Keyword arguments passed to the force field calculator. + **kwargs + Further keyword arguments passed to CTEMaker. A maker given here + replaces the one built for the force field. + + Returns + ------- + CTEMaker + """ + return cls(**{**_get_makers(force_field_name, calculator_kwargs), **kwargs}) + + @property + def prev_calc_dir_argname(self) -> None: + """Name of the argument that passes the previous calculation directory.""" + return diff --git a/src/atomate2/forcefields/utils.py b/src/atomate2/forcefields/utils.py index d6388d845a..2cd486e43e 100644 --- a/src/atomate2/forcefields/utils.py +++ b/src/atomate2/forcefields/utils.py @@ -496,7 +496,12 @@ def revert_default_dtype() -> Generator[None]: Originally added for use with MACE(Relax|Static)Maker. https://github.com/ACEsuit/mace/issues/328 """ - import torch + try: + import torch + except ImportError: + # force fields that do not use torch, such as EMT, have no dtype to revert + yield + return orig = torch.get_default_dtype() yield diff --git a/src/atomate2/vasp/flows/cte.py b/src/atomate2/vasp/flows/cte.py new file mode 100644 index 0000000000..42be710e18 --- /dev/null +++ b/src/atomate2/vasp/flows/cte.py @@ -0,0 +1,93 @@ +"""Define the VASP thermal expansion maker.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +from atomate2.common.flows.cte import BaseCTEMaker +from atomate2.vasp.flows.core import DoubleRelaxMaker +from atomate2.vasp.flows.elastic import ElasticMaker +from atomate2.vasp.flows.pheasy import PhononMaker +from atomate2.vasp.jobs.core import TightRelaxMaker + +if TYPE_CHECKING: + from atomate2.vasp.jobs.base import BaseVaspMaker + + +@dataclass +class CTEMaker(BaseCTEMaker): + """ + Maker to calculate the thermal expansion with VASP, pheasy and phono3py. + + By default, the structure is converted to the standard primitive cell, and + a tight double 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. Neither flow relaxes the + structure again. Finally, phono3py gives the mode Grueneisen tensors, and + the thermal expansion tensor follows from them and the elastic tensor. The + frequencies are not renormalized with temperature. + + By default, the phonon flow fits the force constants with the one-shot + method from randomly displaced supercells with 0.03 A displacements, and + builds the supercells with min_length=12.0. The phonon flow skips the + static energy calculation, which the thermal expansion does not need. The + elastic flow is the atomate2 elastic flow without its own + relaxation. The stress is more sensitive to ENCUT than the forces are, so + check the ENCUT convergence of the elastic tensor for your material. + + 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 pheasy fits can fail + for cells that are not in a standard setting. + bulk_relax_maker: .BaseVaspMaker or None + A maker to perform a tight relaxation on the bulk. Set to None to skip + the relaxation. + phonon_maker: .PhononMaker + 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: .ElasticMaker + 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" + bulk_relax_maker: BaseVaspMaker | None = field( + default_factory=lambda: DoubleRelaxMaker.from_relax_maker(TightRelaxMaker()) + ) + phonon_maker: PhononMaker = field( + default_factory=lambda: PhononMaker( + min_length=12.0, + bulk_relax_maker=None, + static_energy_maker=None, + cal_anhar_fcs=True, + displacement_anhar=0.03, + anhar_fit_methods=("one-shot",), + ) + ) + elastic_maker: ElasticMaker = field( + default_factory=lambda: ElasticMaker(bulk_relax_maker=None) + ) + + @property + def prev_calc_dir_argname(self) -> str: + """Name of the argument that passes the previous calculation directory.""" + return "prev_dir" diff --git a/tests/common/jobs/test_cte.py b/tests/common/jobs/test_cte.py new file mode 100644 index 0000000000..f2173216b5 --- /dev/null +++ b/tests/common/jobs/test_cte.py @@ -0,0 +1,132 @@ +"""Tests for get_cte and the Born charge mapping of the CTE document. + +The force field workflow is tested in tests/forcefields/flows/test_cte.py. +""" + +import numpy as np +import pytest +from ase.build import bulk +from phonopy import Phonopy +from phonopy.structure.atoms import PhonopyAtoms +from pymatgen.analysis.elasticity import ElasticTensor +from scipy.constants import Boltzmann, Planck +from scipy.spatial.transform import Rotation + +from atomate2.common.schemas.cte import _expand_born_to_unitcell, get_cte + + +def _heat_capacity(frequency: float, temperature: float) -> float: + x = Planck * frequency * 1e12 / (Boltzmann * temperature) + return Boltzmann * x**2 * np.exp(x) / (np.exp(x) - 1) ** 2 + + +def test_get_cte_cubic(): + """For a cubic crystal with gamma = g * I, alpha = g * Cv / (3 * B * V).""" + c11, c12, c44 = 170.0, 120.0, 75.0 + elastic = np.zeros((6, 6)) + elastic[:3, :3] = c12 + np.fill_diagonal(elastic[:3, :3], c11) + elastic[3:, 3:] = np.eye(3) * c44 + bulk_modulus = (c11 + 2 * c12) / 3 * 1e9 # Pa + + frequencies = np.array([[2.0, 5.0, 7.0], [3.0, 4.0, 6.0]]) + weights = np.array([1, 3]) + gruneisen = 1.7 * np.broadcast_to(np.eye(3), (2, 3, 3, 3)) + volume, temperature = 45.0, 300.0 + + alpha, mean_gruneisen = get_cte( + frequencies, gruneisen, weights, elastic, volume, [temperature] + ) + + heat_capacity = ( + sum( + w * _heat_capacity(f, temperature) + for w, row in zip(weights, frequencies, strict=True) + for f in row + ) + / weights.sum() + ) + expected = 1.7 * heat_capacity / (3 * bulk_modulus * volume * 1e-30) + assert alpha[0] == pytest.approx(expected * np.eye(3), rel=1e-10, abs=1e-20) + assert mean_gruneisen[0] == pytest.approx(1.7 * np.eye(3)) + + +def test_get_cte_rotation(): + """Rotating the inputs must rotate alpha, which checks the shear terms. + + The mode Grueneisen tensors from phono3py are not symmetric, so the input + tensors here are not symmetric either. + """ + # hexagonal elastic tensor, with C66 = (C11 - C12) / 2 + c11, c12, c13, c33, c44 = 350.0, 120.0, 90.0, 400.0, 110.0 + elastic = np.array( + [ + [c11, c12, c13, 0, 0, 0], + [c12, c11, c13, 0, 0, 0], + [c13, c13, c33, 0, 0, 0], + [0, 0, 0, c44, 0, 0], + [0, 0, 0, 0, c44, 0], + [0, 0, 0, 0, 0, (c11 - c12) / 2], + ] + ) + rng = np.random.default_rng(7) + frequencies = rng.uniform(1.0, 10.0, size=(4, 6)) + gruneisen = rng.normal(1.0, 0.5, size=(4, 6, 3, 3)) + weights = np.ones(4) + temperatures = [100.0, 500.0] + + alpha, _ = get_cte(frequencies, gruneisen, weights, elastic, 30.0, temperatures) + rotation = Rotation.from_euler("zxz", [30, 50, 70], degrees=True).as_matrix() + rotated_elastic = ElasticTensor.from_voigt(elastic).rotate(rotation).voigt + rotated_gruneisen = np.einsum("ik,qbkl,jl->qbij", rotation, gruneisen, rotation) + rotated_alpha, _ = get_cte( + frequencies, rotated_gruneisen, weights, rotated_elastic, 30.0, temperatures + ) + + expected = np.einsum("ik,tkl,jl->tij", rotation, alpha, rotation) + assert np.abs(expected[:, 0, 1]).max() > 1e-7 # the shear terms are tested + assert rotated_alpha == pytest.approx(expected, rel=1e-8, abs=1e-15) + + +def test_get_cte_left_out_modes(): + """Modes below min_frequency are left out, and T = 0 gives zero.""" + elastic = np.diag([200.0, 200.0, 200.0, 80.0, 80.0, 80.0]) + frequencies = np.array([[0.0, 0.0, 0.0, 4.0], [-0.5, 3.0, 5.0, 6.0]]) + gruneisen = np.ones((2, 4, 3, 3)) + gruneisen[0, :3] = np.nan # NaN checks that the left-out modes are ignored + + alpha, mean_gruneisen = get_cte( + frequencies, gruneisen, np.ones(2), elastic, 40.0, [0.0, 300.0] + ) + assert np.all(alpha[0] == 0.0) + assert mean_gruneisen[0] is None + assert mean_gruneisen[1] == pytest.approx(np.ones((3, 3))) + + # only the modes at 4, 3, 5 and 6 THz count, over two q-points + heat_capacity = sum(_heat_capacity(f, 300.0) for f in (4.0, 3.0, 5.0, 6.0)) / 2 + thermal_stress = heat_capacity / (40.0 * 1e-30) + expected = np.full((3, 3), thermal_stress / 80e9 / 2) + np.fill_diagonal(expected, thermal_stress / 200e9) + assert alpha[1] == pytest.approx(expected, rel=1e-10) + + +def test_expand_born_to_unitcell(): + """Born charges of the primitive cell are mapped onto the conventional cell.""" + atoms = bulk("MgO", "rocksalt", a=4.2, cubic=True) + unitcell = PhonopyAtoms( + symbols=atoms.get_chemical_symbols(), + cell=atoms.cell[:], + scaled_positions=atoms.get_scaled_positions(), + ) + phonon = Phonopy(unitcell, supercell_matrix=np.eye(3), primitive_matrix="auto") + assert len(phonon.primitive) == 2 + born_primitive = {"Mg": 1.9, "O": -1.9} + phonon.nac_params = { + "born": [np.eye(3) * born_primitive[s] for s in phonon.primitive.symbols], + "dielectric": np.eye(3) * 3.0, + } + + born = _expand_born_to_unitcell(phonon) + assert born.shape == (8, 3, 3) + for symbol, charges in zip(unitcell.symbols, born, strict=True): + assert charges == pytest.approx(np.eye(3) * born_primitive[symbol]) diff --git a/tests/common/jobs/test_pheasy.py b/tests/common/jobs/test_pheasy.py index 58b72dc7a0..fc1d816d29 100644 --- a/tests/common/jobs/test_pheasy.py +++ b/tests/common/jobs/test_pheasy.py @@ -29,7 +29,6 @@ _get_num_irreducible_fcs, _run_band_structure_and_plot, ) -from atomate2.forcefields.flows.pheasy import PhononMaker # fcs_cutoff_radius in Bohr. 8 Bohr (4.2 A) covers the first two neighbour # shells of fcc Cu (2.55 and 3.61 A) and stays inside the 10.8 A supercell. @@ -261,28 +260,6 @@ def test_get_num_anharmonic_supercells(monkeypatch): _get_num_anharmonic_supercells(num_disp_anhar=0, **kwargs) -def test_get_supercell_size_kwargs(monkeypatch): - received = {} - transformation = pheasy_jobs.CubicSupercellTransformation - - def record_kwargs(**kwargs): - received.update(kwargs) - return transformation(**kwargs) - - monkeypatch.setattr(pheasy_jobs, "CubicSupercellTransformation", record_kwargs) - - # the maker passes get_supercell_size_kwargs on to the job - maker = PhononMaker(get_supercell_size_kwargs={"angle_tolerance": 0.1}) - job = maker.get_supercell_matrix(_cu_structure()) - assert job.function_kwargs == {"angle_tolerance": 0.1} - - # the job passes them to CubicSupercellTransformation. The other default - # is kept. - job.function(*job.function_args, **job.function_kwargs) - assert received["angle_tolerance"] == 0.1 - assert received["allow_orthorhombic"] is False - - def test_check_lasso_alpha(tmp_dir): log_file = Path("pheasy_anharmonic_fit.log") diff --git a/tests/forcefields/flows/test_cte.py b/tests/forcefields/flows/test_cte.py new file mode 100644 index 0000000000..8d4a493df3 --- /dev/null +++ b/tests/forcefields/flows/test_cte.py @@ -0,0 +1,174 @@ +"""Tests for the force field thermal expansion workflow. + +The tests run the force field CTEMaker with ASE's EMT potential, so they need no +DFT reference data and no machine-learned force field. pheasy, ALM and phono3py +are only installed in the numpy-limited forcefield CI job, so the tests are +skipped in the other forcefield jobs. +""" +# ruff: noqa: E402 + +from pathlib import Path + +import pytest + +pytest.importorskip("pheasy") +gruneisen_module = pytest.importorskip("phono3py.phonon3.gruneisen") + +import numpy as np +import phonopy +from ase.build import bulk +from jobflow import run_locally +from pymatgen.io.ase import AseAtomsAdaptor + +from atomate2.common.jobs.cte import compute_cte +from atomate2.common.schemas.cte import CTEDocument +from atomate2.forcefields.flows.cte import CTEMaker + +EMT = {"@module": "ase.calculators.emt", "@callable": "EMT"} + + +def test_cte_maker_emt(clean_dir, monkeypatch): + """Run the whole force field workflow with EMT forces on fcc Cu.""" + structure = AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + # the conventional cell keeps the 4-atom cubic unit cell used below + maker = CTEMaker.from_force_field_name( + EMT, + use_symmetrized_structure="conventional", + temperatures=[0, 100, 300], + mesh=(8, 8, 8), + ) + # a 2x2x2 supercell of the cubic cell, 32 atoms, and an fc3 cutoff of 6 Bohr + # (3.2 A), which covers the nearest neighbours at 2.55 A + maker.phonon_maker.min_length = 7.0 + maker.phonon_maker.fcs_cutoff_radius = [-1, 6, 6] + maker.phonon_maker.anhar_fit_methods = ["cocktail", "one-shot"] + + flow = maker.make(structure) + # the uuid of the phonon flow output that compute_cte reads + phonon_uuid = flow.jobs[-1].function_kwargs["phonon_output"].uuid + responses = run_locally(flow, create_folders=True, ensure_success=True) + doc = responses[flow.output.uuid][1].output + assert isinstance(doc, CTEDocument) + assert doc.mesh == (8, 8, 8) + assert np.array(doc.supercell_matrix) == pytest.approx(2 * np.eye(3)) + assert [result.fit_method for result in doc.results] == ["cocktail", "one-shot"] + + # EMT elastic constants of Cu in GPa + assert doc.elastic_tensor[0][0] == pytest.approx(172.6, rel=0.01) + assert doc.elastic_tensor[0][1] == pytest.approx(115.4, rel=0.01) + assert doc.elastic_tensor[3][3] == pytest.approx(89.9, rel=0.01) + + for result in doc.results: + assert not result.has_imaginary_modes + alpha = np.array(result.thermal_expansion_tensor) + assert np.all(alpha[0] == 0.0) + # cubic, so alpha is isotropic in every frame + assert alpha[2] == pytest.approx(alpha[2, 0, 0] * np.eye(3), abs=1e-12) + assert result.thermal_expansion[2] == pytest.approx(3 * alpha[2, 0, 0]) + # linear thermal expansion at 300 K in 1/K, and the mean Grueneisen parameter. + # The one-shot fit has few supercells in this small cell, so its value is only + # compared with the cocktail value. + cocktail, one_shot = doc.results + assert cocktail.thermal_expansion_tensor[2][0][0] == pytest.approx( + 1.859e-5, rel=0.02 + ) + assert np.trace(cocktail.average_gruneisen[2]) / 3 == pytest.approx(2.237, rel=0.02) + assert one_shot.thermal_expansion_tensor[2][0][0] == pytest.approx( + cocktail.thermal_expansion_tensor[2][0][0], rel=0.2 + ) + assert Path(doc.phonon_job_dir, "one_shot", "fc3.hdf5").exists() + + # the same force constants, with a q-point density instead of a mesh. The + # negative tolerance flags every frequency, to test the imaginary mode check. + phonon_output = responses[phonon_uuid][1].output + job = compute_cte( + phonon_output=phonon_output, + elastic_tensor=doc.elastic_tensor, + elastic_structure=doc.structure, + anhar_fit_methods=["one-shot"], + temperatures=[300], + mesh=100.0, + tol_imaginary_modes=-10.0, + ) + with pytest.warns(UserWarning, match="thermal expansion is not"): + responses = run_locally(job, create_folders=True, ensure_success=True) + flagged_doc = responses[job.uuid][1].output + assert flagged_doc.mesh == (2, 2, 2) + (flagged,) = flagged_doc.results + assert flagged.has_imaginary_modes + assert flagged.thermal_expansion_tensor is None + + # the non-analytical term correction with zero Born charges leaves alpha + # unchanged. pheasy only stores Born charges for VASP, so they are added here. + # The phonopy primitive cell has one atom and the unit cell four, so the + # charges must be expanded to four atoms before they reach phono3py. + phonon_job_dir = Path(doc.phonon_job_dir) + phonon = phonopy.load(phonon_job_dir / "phonopy.yaml", produce_fc=False) + assert len(phonon.primitive) == 1 + phonon.nac_params = { + "born": np.zeros((1, 3, 3)), + "dielectric": np.eye(3) * 10.0, + "factor": 14.399652, + } + phonon.save("phonopy_nac.yaml") + + nac_params_used = [] + original_gruneisen = gruneisen_module.Gruneisen + + def _gruneisen(*args, **kwargs): + nac_params_used.append(kwargs["nac_params"]) + return original_gruneisen(*args, **kwargs) + + monkeypatch.setattr(gruneisen_module, "Gruneisen", _gruneisen) + cocktail_files = { + "cocktail": (phonon_job_dir / "FORCE_CONSTANTS", phonon_job_dir / "fc3.hdf5") + } + settings = { + "force_constant_files": cocktail_files, + "elastic_tensor": doc.elastic_tensor, + "temperatures": [300], + "mesh": (8, 8, 8), + "tol_imaginary_modes": 0.1, + "min_frequency": 1e-3, + "symprec": 1e-5, + } + nac_doc = CTEDocument.from_force_constants( + phonopy_yaml="phonopy_nac.yaml", structure=doc.structure, **settings + ) + (nac_params,) = nac_params_used + assert nac_params["born"].shape == (4, 3, 3) + assert nac_params["dielectric"] == pytest.approx(np.eye(3) * 10.0) + assert np.array(nac_doc.results[0].thermal_expansion_tensor[0]) == pytest.approx( + np.array(cocktail.thermal_expansion_tensor[2]), rel=1e-6, abs=1e-15 + ) + + # a structure with another lattice or other atoms is refused + strained = doc.structure.copy() + strained.apply_strain(0.01) + with pytest.raises(ValueError, match="not in the same frame"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", structure=strained, **settings + ) + other_atoms = doc.structure.copy() + other_atoms.replace_species({"Cu": "Au"}) + with pytest.raises(ValueError, match="atoms of the elastic calculation"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", + structure=other_atoms, + **settings, + ) + with pytest.raises(ValueError, match="must not be negative"): + CTEDocument.from_force_constants( + phonopy_yaml=phonon_job_dir / "phonopy.yaml", + structure=doc.structure, + **{**settings, "temperatures": [-10, 300]}, + ) + + +def test_cte_maker_force_field_defaults(): + """The default phonon maker uses min_length=12.0 for the supercells.""" + maker = CTEMaker() + assert maker.use_symmetrized_structure == "primitive" + assert maker.phonon_maker.min_length == 12.0 + assert maker.phonon_maker.cal_anhar_fcs + assert maker.phonon_maker.displacement_anhar == 0.03 diff --git a/tests/forcefields/flows/test_pheasy.py b/tests/forcefields/flows/test_pheasy.py new file mode 100644 index 0000000000..3315963a1b --- /dev/null +++ b/tests/forcefields/flows/test_pheasy.py @@ -0,0 +1,78 @@ +"""Tests for the force field pheasy workflow. + +pheasy and ALM are only installed in the numpy-limited forcefield CI job, so the +tests are skipped in the other forcefield jobs. +""" + +import pytest + +pytest.importorskip("pheasy") + +from ase.build import bulk +from emmet.core.phonon import ( + PhononBS, + PhononBSDOSDoc, + PhononDOS, + ThermalDisplacementData, +) +from jobflow import run_locally +from pymatgen.core import Structure +from pymatgen.io.ase import AseAtomsAdaptor + +import atomate2.common.jobs.pheasy as pheasy_jobs +from atomate2.forcefields.flows.pheasy import PhononMaker +from atomate2.forcefields.jobs import ForceFieldRelaxMaker, ForceFieldStaticMaker + +EMT = {"@module": "ase.calculators.emt", "@callable": "EMT"} + + +def _cu_structure() -> Structure: + return AseAtomsAdaptor.get_structure(bulk("Cu", "fcc", a=3.61, cubic=True)) + + +def test_get_supercell_size_kwargs(monkeypatch): + received = {} + transformation = pheasy_jobs.CubicSupercellTransformation + + def record_kwargs(**kwargs): + received.update(kwargs) + return transformation(**kwargs) + + monkeypatch.setattr(pheasy_jobs, "CubicSupercellTransformation", record_kwargs) + + # the maker passes get_supercell_size_kwargs on to the job + maker = PhononMaker(get_supercell_size_kwargs={"angle_tolerance": 0.1}) + job = maker.get_supercell_matrix(_cu_structure()) + assert job.function_kwargs == {"angle_tolerance": 0.1} + + # the job passes them to CubicSupercellTransformation. The other default + # is kept. + job.function(*job.function_args, **job.function_kwargs) + assert received["angle_tolerance"] == 0.1 + assert received["allow_orthorhombic"] is False + + +def test_pheasy_wf_force_field(clean_dir): + """Run the harmonic force field pheasy workflow with EMT forces on fcc Cu.""" + # a 2x2x2 supercell of the 4-atom cubic cell + maker = PhononMaker( + min_length=7.0, + use_symmetrized_structure="conventional", + create_thermal_displacements=True, + bulk_relax_maker=ForceFieldRelaxMaker( + force_field_name=EMT, relax_kwargs={"fmax": 0.00001} + ), + static_energy_maker=ForceFieldStaticMaker(force_field_name=EMT), + phonon_displacement_maker=ForceFieldStaticMaker(force_field_name=EMT), + ) + flow = maker.make(_cu_structure()) + responses = run_locally(flow, create_folders=True, ensure_success=True) + ph_doc = responses[flow.jobs[-1].uuid][1].output + + assert isinstance(ph_doc, PhononBSDOSDoc) + assert isinstance(ph_doc.phonon_bandstructure, PhononBS) + assert isinstance(ph_doc.phonon_dos, PhononDOS) + assert isinstance(ph_doc.thermal_displacement_data, ThermalDisplacementData) + assert isinstance(ph_doc.structure, Structure) + assert ph_doc.has_imaginary_modes is False + assert isinstance(ph_doc.force_constants, list) diff --git a/tests/vasp/flows/test_cte.py b/tests/vasp/flows/test_cte.py new file mode 100644 index 0000000000..98c392dde6 --- /dev/null +++ b/tests/vasp/flows/test_cte.py @@ -0,0 +1,73 @@ +import pytest +from jobflow import Flow +from pymatgen.core.structure import Structure + +from atomate2.vasp.flows.cte import CTEMaker +from atomate2.vasp.flows.elastic import ElasticMaker +from atomate2.vasp.flows.pheasy import PhononMaker + + +def test_cte_maker_vasp_flow(si_structure: Structure): + """The relaxed structure goes to both flows, and their outputs to compute_cte.""" + maker = CTEMaker(temperatures=[100, 300], mesh=(10, 10, 10)) + flow = maker.make(si_structure) + prim, relax, phonon_flow, elastic_flow, cte = flow.jobs + assert prim.name == "structure_to_primitive" + assert relax.jobs[0].function_args[0].uuid == prim.output.uuid + assert isinstance(phonon_flow, Flow) + assert isinstance(elastic_flow, Flow) + + # both flows start from the relaxed structure + relaxed = relax.output.structure + for sub_flow in (phonon_flow, elastic_flow): + structure = sub_flow.jobs[0].function_args[0] + assert structure.uuid == relaxed.uuid + assert structure.attributes == relaxed.attributes + + kwargs = cte.function_kwargs + assert kwargs["phonon_output"].uuid == phonon_flow.output.uuid + assert kwargs["elastic_tensor"].uuid == elastic_flow.output.uuid + assert kwargs["elastic_structure"].uuid == elastic_flow.output.uuid + assert kwargs["anhar_fit_methods"] == ("one-shot",) + assert kwargs["temperatures"] == [100, 300] + assert kwargs["mesh"] == (10, 10, 10) + assert kwargs["symprec"] == maker.phonon_maker.symprec + assert kwargs["min_frequency"] == maker.min_frequency + assert maker.phonon_maker.min_length == 12.0 + + assert kwargs["elastic_tensor"].attributes == ( + ("a", "elastic_tensor"), + ("a", "raw"), + ) + + # the elastic fit gets the stress of the relaxation, as in the elastic flow + fit = next(job for job in elastic_flow.jobs if job.name == "fit_elastic_tensor") + stress = fit.function_kwargs["equilibrium_stress"] + assert stress.uuid == relax.output.uuid + assert stress.attributes == (("a", "output"), ("a", "stress")) + + +@pytest.mark.parametrize( + ("phonon_kwargs", "match"), + [ + ({"bulk_relax_maker": None, "cal_anhar_fcs": False}, "cal_anhar_fcs=True"), + ( + { + "bulk_relax_maker": None, + "cal_anhar_fcs": True, + "use_symmetrized_structure": "primitive", + }, + "use_symmetrized_structure=None", + ), + ({"cal_anhar_fcs": True}, "The phonon maker needs bulk_relax_maker=None"), + ], +) +def test_cte_maker_checks_phonon_maker(phonon_kwargs, match): + phonon_maker = PhononMaker(**phonon_kwargs) + with pytest.raises(ValueError, match=match): + CTEMaker(phonon_maker=phonon_maker) + + +def test_cte_maker_checks_elastic_maker(): + with pytest.raises(ValueError, match="The elastic maker needs bulk_relax_maker"): + CTEMaker(elastic_maker=ElasticMaker()) diff --git a/tutorials/cte_workflow.ipynb b/tutorials/cte_workflow.ipynb new file mode 100644 index 0000000000..92572b3b9f --- /dev/null +++ b/tutorials/cte_workflow.ipynb @@ -0,0 +1,771 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Thermal expansion workflow with pheasy and machine-learned potentials\n", + "\n", + "This notebook computes the thermal expansion with the `CTEMaker` workflow. The\n", + "forces come from three MACE potentials, so no DFT is needed.\n", + "\n", + "Section 1 runs five materials with the default settings. MgO is worked through\n", + "step by step. Section 2 converges the third-order cutoff and the supercell size\n", + "for Si, where the default settings are not enough." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Background\n", + "\n", + "The workflow computes the thermal expansion tensor from mode Grüneisen tensors.\n", + "The thermal stress of each phonon mode is its heat capacity times its Grüneisen\n", + "tensor. The strain that relaxes the summed thermal stress is the thermal\n", + "expansion:\n", + "\n", + "$$\\alpha = \\frac{S}{N_q V} \\sum_{q\\nu} c_{q\\nu} \\gamma_{q\\nu}$$\n", + "\n", + "Here $S$ is the elastic compliance, $c_{q\\nu}$ and $\\gamma_{q\\nu}$ are the heat\n", + "capacity and Grüneisen tensor of each mode, $N_q$ is the number of q-points and\n", + "$V$ is the cell volume.\n", + "\n", + "The flow has four steps:\n", + "\n", + "1. The structure is converted to the standard primitive cell and relaxed tightly.\n", + "2. The pheasy phonon workflow fits the second- and third-order force constants\n", + " with LASSO. It uses randomly displaced supercells with 0.03 Å displacements.\n", + "3. The elastic workflow fits the elastic tensor on the same relaxed structure.\n", + "4. phono3py computes the mode Grüneisen tensors on a 12x12x12 q-point mesh, and\n", + " the thermal expansion follows from the formula above.\n", + "\n", + "Steps 2 and 3 do not depend on each other, so a workflow manager can run them\n", + "at the same time. The same workflow runs with VASP as\n", + "`atomate2.vasp.flows.cte.CTEMaker`." + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Installation\n", + "\n", + "The `pheasy`, `phono3py` and `ase` extras are needed, plus MACE and the\n", + "Materials Project client:\n", + "\n", + "```\n", + "pip install 'atomate2[pheasy,phono3py,ase]'\n", + "pip install 'mace-torch>=0.3.16' mp-api\n", + "```\n", + "\n", + "The `pheasy` extra pulls in pheasy and ALM, and the `phono3py` extra pulls in\n", + "phono3py. ALM is compiled from source. The hiPhive tutorial and the atomate2\n", + "VASP documentation list what the build needs." + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "## The potentials\n", + "\n", + "The notebook uses three MACE foundation models. All three are released under the\n", + "Academic Software License.\n", + "\n", + "- **MACE-OMAT-0-medium.** `mace_mp(model=\"medium-omat-0\")` downloads it once and\n", + " caches it under `~/.cache/mace`.\n", + "- **MACE-MATPES-PBE-0** and **MACE-MATPES-r2SCAN-0.** These are MACE-OMAT-0\n", + " fine-tuned on the MatPES data set at the PBE and r2SCAN levels. They need\n", + " `mace-torch>=0.3.10`. Download the two files from the `mace_matpes_0` release of\n", + " [mace-foundations](https://github.com/ACEsuit/mace-foundations/releases/tag/mace_matpes_0):\n", + " - [MACE-matpes-pbe-omat-ft.model](https://github.com/ACEsuit/mace-foundations/releases/download/mace_matpes_0/MACE-matpes-pbe-omat-ft.model)\n", + " - [MACE-matpes-r2scan-omat-ft.model](https://github.com/ACEsuit/mace-foundations/releases/download/mace_matpes_0/MACE-matpes-r2scan-omat-ft.model)\n", + "\n", + "Set the paths to the two downloaded files in the next cell. The MgO example\n", + "below only needs MACE-OMAT-0-medium. Change `device` to `\"cuda\"` to run MACE on a\n", + "GPU.\n", + "\n", + "`float64` is worth the cost here. The third-order force constants come from\n", + "small force differences between displaced supercells, where `float32` noise is\n", + "not negligible." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import warnings\n", + "\n", + "# macOS: conda's llvm-openmp and torch's bundled libomp both load and the\n", + "# duplicate aborts the process. This must be set before torch is imported.\n", + "os.environ.setdefault(\"KMP_DUPLICATE_LIB_OK\", \"TRUE\")\n", + "warnings.filterwarnings(\"ignore\")\n", + "\n", + "from mace.calculators import mace_mp # noqa: E402\n", + "\n", + "# Replace the two paths with the files you downloaded.\n", + "MODEL_FILES = {\n", + " \"OMAT-0-medium\": \"medium-omat-0\",\n", + " \"MATPES-PBE-0\": \"/path/to/MACE-matpes-pbe-omat-ft.model\",\n", + " \"MATPES-r2SCAN-0\": \"/path/to/MACE-matpes-r2scan-omat-ft.model\",\n", + "}\n", + "\n", + "\n", + "def calc_kwargs(model: str) -> dict:\n", + " \"\"\"Return the calculator settings for one of the models above.\"\"\"\n", + " return {\"model\": MODEL_FILES[model], \"device\": \"cpu\", \"default_dtype\": \"float64\"}\n", + "\n", + "\n", + "_ = mace_mp(**calc_kwargs(\"OMAT-0-medium\")) # downloads and caches on first use" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "## One calculator for all jobs\n", + "\n", + "Every force field job builds its own calculator. `run_locally` runs all of\n", + "them in this one Python process, and the workflow has a few hundred jobs.\n", + "Loading MACE for each of them grows the memory by about 150 MB per job. This\n", + "cell makes the jobs share one calculator. It is only needed when many jobs run\n", + "in one process, not with a workflow manager." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "from ase.calculators.calculator import Calculator\n", + "\n", + "import atomate2.forcefields.utils as ff_utils\n", + "\n", + "_ase_calculator = ff_utils.ase_calculator\n", + "_calculators: dict[str, Calculator] = {}\n", + "\n", + "\n", + "def shared_ase_calculator(calculator_meta: object, **kwargs: object) -> Calculator:\n", + " \"\"\"Build each calculator once and reuse it for every job.\"\"\"\n", + " key = repr((calculator_meta, sorted(kwargs.items())))\n", + " if key not in _calculators:\n", + " _calculators[key] = _ase_calculator(calculator_meta, **kwargs)\n", + " _calculators[key].reset()\n", + " return _calculators[key]\n", + "\n", + "\n", + "ff_utils.ase_calculator = shared_ase_calculator" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "## 1. Five materials with the default settings\n", + "\n", + "The first part runs MgO with MACE-OMAT-0-medium step by step. The other four\n", + "materials and the other two potentials follow with the same code." + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## The structure\n", + "\n", + "The primitive cell of MgO (mp-1265) comes from the Materials Project. Set your\n", + "key first with `export MP_API_KEY=...`.\n", + "\n", + "The Materials Project serves this cell rotated away from the cubic axes, and\n", + "the pheasy fits fail for such a cell. `CTEMaker` therefore converts the input to\n", + "the standard primitive cell before the relaxation. This is the default,\n", + "`use_symmetrized_structure=\"primitive\"`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "from mp_api.client import MPRester\n", + "\n", + "if not os.environ.get(\"MP_API_KEY\"):\n", + " raise OSError(\"export MP_API_KEY before running this cell\")\n", + "\n", + "with MPRester() as mpr:\n", + " structure = mpr.get_structure_by_material_id(\"mp-1265\")\n", + "structure.lattice" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## Building the workflow\n", + "\n", + "`CTEMaker.from_force_field_name` sets one force field for the relaxation, the\n", + "phonons and the elastic tensor. Any MACE model routes through `mace_mp`, and\n", + "the model name in `calculator_kwargs` selects the weights.\n", + "\n", + "Everything else stays at the workflow defaults. The supercells have\n", + "`min_length=12` Å, which gives 128 atoms for MgO. The third-order force\n", + "constants use the one-shot fit, and the temperatures run from 0 K to 1000 K in\n", + "steps of 10 K." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "from atomate2.forcefields.flows.cte import CTEMaker\n", + "\n", + "maker = CTEMaker.from_force_field_name(\n", + " \"MACE-MP-0\", calculator_kwargs=calc_kwargs(\"OMAT-0-medium\")\n", + ")\n", + "flow = maker.make(structure)\n", + "flow.draw_graph().show()" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "## Running the workflow\n", + "\n", + "The flow has a little over 200 jobs, most of them force calculations on the\n", + "displaced supercells. On 32 CPU cores of a Perlmutter node the whole run took\n", + "about 8 minutes.\n", + "\n", + "`create_folders=True` is needed, because the thermal expansion job reads the\n", + "force constants from the folder of the pheasy fit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "from jobflow import JobStore, run_locally\n", + "from maggma.stores import MemoryStore\n", + "\n", + "job_store = JobStore(MemoryStore(), additional_stores={\"data\": MemoryStore()})\n", + "responses = run_locally(flow, store=job_store, create_folders=True, ensure_success=True)\n", + "doc = responses[flow.output.uuid][1].output" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "## Results\n", + "\n", + "The output is a `CTEDocument`. It holds the thermal expansion tensor at each\n", + "temperature for each fit method. The value at 0 K is zero, since the heat\n", + "capacity is zero there. MgO is cubic, so the three diagonal components are\n", + "equal." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "(result,) = doc.results\n", + "temperatures = np.array(doc.temperatures)\n", + "alpha = np.array(result.thermal_expansion_tensor)\n", + "i300 = int(np.argmin(np.abs(temperatures - 300)))\n", + "\n", + "fig, ax = plt.subplots()\n", + "ax.plot(temperatures, alpha[:, 0, 0] * 1e6)\n", + "ax.set_xlabel(\"Temperature (K)\")\n", + "ax.set_ylabel(\"Linear thermal expansion (10$^{-6}$ K$^{-1}$)\")\n", + "plt.show()\n", + "\n", + "{\n", + " \"alpha(300 K) in 1/K\": float(alpha[i300, 0, 0]),\n", + " \"mean Grueneisen parameter at 300 K\": float(\n", + " np.trace(result.average_gruneisen[i300]) / 3\n", + " ),\n", + " \"has imaginary modes\": result.has_imaginary_modes,\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "With these settings we get $\\alpha$(300 K) = 11.2e-6 K⁻¹. Experiment gives\n", + "about 1.0e-5 K⁻¹ for MgO at room temperature." + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Checking the fit\n", + "\n", + "If the thermal expansion comes out as zero or very small, check the\n", + "third-order force constants first. All zeros means the LASSO fit dropped them.\n", + "The penalty chosen by cross validation should also lie inside the search\n", + "range, not on one of its bounds." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "import re\n", + "from pathlib import Path\n", + "\n", + "import h5py\n", + "\n", + "fit_dir = Path(doc.phonon_job_dir) / \"one_shot\"\n", + "log = (fit_dir / \"pheasy_anharmonic_fit.log\").read_text()\n", + "with h5py.File(fit_dir / \"fc3.hdf5\") as file:\n", + " fc3 = file[\"fc3\"][:]\n", + "\n", + "{\n", + " \"LASSO penalty\": re.findall(r\"alpha_(?:min|max|opt):\\s*\\S+\", log),\n", + " \"max |fc3| in eV/A^3\": float(np.abs(fc3).max()),\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "## The other materials and potentials\n", + "\n", + "The same workflow runs for NaCl, KCl, CaO and GaAs, and with all three\n", + "potentials. The function below builds and runs one workflow in its own folder.\n", + "Section 2 uses it as well.\n", + "\n", + "The loop runs 15 workflows, so it is off by default. Set `RUN_ALL = True` to\n", + "run it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "from pymatgen.core import Structure\n", + "\n", + "from atomate2.common.schemas.cte import CTEDocument\n", + "\n", + "MATERIALS = {\n", + " \"NaCl\": \"mp-22862\",\n", + " \"KCl\": \"mp-23193\",\n", + " \"MgO\": \"mp-1265\",\n", + " \"CaO\": \"mp-2605\",\n", + " \"GaAs\": \"mp-2534\",\n", + "}\n", + "\n", + "\n", + "def run_cte(\n", + " structure: Structure,\n", + " model: str,\n", + " root_dir: str,\n", + " min_length: float | None = None,\n", + " c3: float | None = None,\n", + ") -> CTEDocument:\n", + " \"\"\"Run the thermal expansion workflow and return the CTEDocument.\n", + "\n", + " min_length sets the supercell size in Angstrom. c3 sets the third-order\n", + " cutoff in Bohr. None keeps the workflow default.\n", + " \"\"\"\n", + " maker = CTEMaker.from_force_field_name(\n", + " \"MACE-MP-0\", calculator_kwargs=calc_kwargs(model)\n", + " )\n", + " if min_length is not None:\n", + " maker.phonon_maker.min_length = min_length\n", + " if c3 is not None:\n", + " maker.phonon_maker.fcs_cutoff_radius = [-1, c3, 10]\n", + " flow = maker.make(structure)\n", + " os.makedirs(root_dir, exist_ok=True)\n", + " responses = run_locally(\n", + " flow,\n", + " store=JobStore(MemoryStore(), additional_stores={\"data\": MemoryStore()}),\n", + " create_folders=True,\n", + " root_dir=root_dir,\n", + " ensure_success=True,\n", + " )\n", + " return responses[flow.output.uuid][1].output\n", + "\n", + "\n", + "def alpha_300k(doc: CTEDocument) -> float:\n", + " \"\"\"Linear thermal expansion at 300 K in 1/K, from the xx component.\"\"\"\n", + " (result,) = doc.results\n", + " i300 = int(np.argmin(np.abs(np.array(doc.temperatures) - 300)))\n", + " return float(np.array(result.thermal_expansion_tensor)[i300, 0, 0])\n", + "\n", + "\n", + "RUN_ALL = False\n", + "if RUN_ALL:\n", + " with MPRester() as mpr:\n", + " structures = {\n", + " name: mpr.get_structure_by_material_id(mp_id)\n", + " for name, mp_id in MATERIALS.items()\n", + " }\n", + " alpha_all = {\n", + " (name, model): alpha_300k(run_cte(structure, model, f\"runs/{name}_{model}\"))\n", + " for name, structure in structures.items()\n", + " for model in MODEL_FILES\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "# alpha(300 K) in 1e-6/K from our runs with the default settings\n", + "ALPHA_DEFAULT = {\n", + " \"OMAT-0-medium\": {\"NaCl\": 39.4, \"KCl\": 41.6, \"MgO\": 11.2, \"CaO\": 13.9, \"GaAs\": 5.9},\n", + " \"MATPES-PBE-0\": {\"NaCl\": 44.4, \"KCl\": 47.5, \"MgO\": 31.4, \"CaO\": 14.3, \"GaAs\": 4.4},\n", + " \"MATPES-r2SCAN-0\": {\n", + " \"NaCl\": 32.8,\n", + " \"KCl\": 17.7,\n", + " \"MgO\": 11.3,\n", + " \"CaO\": 11.1,\n", + " \"GaAs\": 2.5,\n", + " },\n", + "}\n", + "# finite displacements in the largest supercell we ran, 686 to 2662 atoms\n", + "ALPHA_FD_LARGE = {\n", + " \"OMAT-0-medium\": {\"NaCl\": 39.3, \"KCl\": 41.8, \"MgO\": 11.0, \"CaO\": 14.0, \"GaAs\": 6.7},\n", + " \"MATPES-PBE-0\": {\"NaCl\": 40.0, \"KCl\": 44.4, \"MgO\": 39.2, \"CaO\": 15.8, \"GaAs\": 6.7},\n", + " \"MATPES-r2SCAN-0\": {\n", + " \"NaCl\": 33.7,\n", + " \"KCl\": 19.4,\n", + " \"MgO\": 12.3,\n", + " \"CaO\": 11.8,\n", + " \"GaAs\": 4.1,\n", + " },\n", + "}\n", + "\n", + "pd.concat(\n", + " {\n", + " \"default settings\": pd.DataFrame(ALPHA_DEFAULT),\n", + " \"finite displacements\": pd.DataFrame(ALPHA_FD_LARGE),\n", + " },\n", + " axis=1,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "The first set of columns is what the workflow gives with the default settings.\n", + "To check it, we computed the force constants of each run by finite\n", + "displacements with phono3py, with no cutoff and no fit. In the supercell of\n", + "the run the fits agree with finite displacements within 6%. The one exception\n", + "is GaAs with MATPES-r2SCAN-0, at 10%. So the fit is right for the supercell it\n", + "uses.\n", + "\n", + "We then repeated the finite displacements in larger supercells. The second set\n", + "of columns gives the value in the largest one, with 686 to 2662 atoms. Between\n", + "the two largest supercells each value changes by about 3% or less.\n", + "\n", + "- With MACE-OMAT-0-medium the default settings are within 2% of the\n", + " large-supercell value for the four rock-salt materials.\n", + "- With the two MatPES potentials the default settings are off by up to 11% for\n", + " the rock-salt materials, and by 20% for MgO with MATPES-PBE-0.\n", + "- For GaAs the default settings are 11% to 38% too low with all three\n", + " potentials. GaAs has the zinc-blende structure of Si. Section 2 shows the\n", + " same effect for Si in more detail.\n", + "\n", + "The potentials differ much more than these errors. Two results stand out.\n", + "\n", + "- KCl with MATPES-r2SCAN-0 gives 19.4e-6 K⁻¹ in the large supercell. This is\n", + " less than half the value of the other two potentials.\n", + "- MgO with MATPES-PBE-0 gives 39.2e-6 K⁻¹. This is more than three times the\n", + " value of the other two potentials and of experiment. MATPES-PBE-0 also gives\n", + " a much softer MgO. Its bulk modulus is 100 GPa, against 154 GPa with\n", + " MACE-OMAT-0-medium and 163 GPa with MATPES-r2SCAN-0. A softer lattice\n", + " expands more." + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "## 2. Converging the cutoff and the supercell for Si\n", + "\n", + "In Si the transverse acoustic modes near the zone boundary have negative\n", + "Grüneisen parameters. Their contribution to the thermal stress partly cancels\n", + "that of the other modes. The thermal expansion is the small difference of two\n", + "larger terms. Small errors in the third-order force constants therefore change\n", + "it a lot. This is why the default settings are not enough for Si.\n", + "\n", + "Two settings control the third-order force constants:\n", + "\n", + "- `phonon_maker.min_length` sets the supercell. For the Si primitive cell, 12, 16\n", + " and 21 Å give the 4x4x4, 5x5x5 and 6x6x6 supercells with 128, 250 and 432\n", + " atoms. The default is 12 Å.\n", + "- `phonon_maker.fcs_cutoff_radius` sets the cutoff radius of each order in Bohr.\n", + " The second entry is the third-order cutoff. The default is\n", + " `[-1, 12, 10]`, so 12 Bohr.\n", + "\n", + "A cutoff should stay below half the shortest distance between periodic images\n", + "of the supercell. Above it an atom starts to interact with its own images. For\n", + "Si this limit is 14.6, 18.3 and 21.9 Bohr for the 4x4x4, 5x5x5 and 6x6x6\n", + "supercells.\n", + "\n", + "A larger cutoff adds third-order force constants. The workflow picks the number\n", + "of displaced supercells so that the fit has about 100 equations per free force\n", + "constant, up to 600 supercells. The cost therefore grows quickly with the\n", + "cutoff. The 6x6x6 run at 18 Bohr used 201 displaced supercells. At 20 Bohr it used 527,\n", + "and the LASSO fit needed more than the 220 GB of memory we gave it. We stopped\n", + "at 18 Bohr.\n", + "\n", + "The scan below runs 11 settings for each potential, so it is off by default." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "SI_SCAN = { # min_length in Angstrom: third-order cutoffs in Bohr\n", + " 12: (12, 14, 16),\n", + " 16: (12, 14, 16, 18),\n", + " 21: (12, 14, 16, 18),\n", + "}\n", + "\n", + "RUN_SI_SCAN = False\n", + "if RUN_SI_SCAN:\n", + " with MPRester() as mpr:\n", + " si = mpr.get_structure_by_material_id(\"mp-149\")\n", + " alpha_si = {\n", + " (model, min_length, c3): alpha_300k(\n", + " run_cte(\n", + " si, model, f\"runs/Si_{model}_{min_length}A_{c3}bohr\", min_length, c3\n", + " )\n", + " )\n", + " for model in MODEL_FILES\n", + " for min_length, cutoffs in SI_SCAN.items()\n", + " for c3 in cutoffs\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "# alpha(300 K) in 1e-6/K from our runs, by potential, supercell and cutoff in Bohr\n", + "ALPHA_SI = {\n", + " \"OMAT-0-medium\": {\n", + " \"4x4x4\": {12: 1.42, 14: 1.90, 16: 1.87},\n", + " \"5x5x5\": {12: 1.53, 14: 2.04, 16: 2.27, 18: 2.06},\n", + " \"6x6x6\": {12: 1.64, 14: 1.99, 16: 2.30, 18: 2.07},\n", + " },\n", + " \"MATPES-PBE-0\": {\n", + " \"4x4x4\": {12: 3.25, 14: 4.13, 16: 4.06},\n", + " \"5x5x5\": {12: 3.35, 14: 4.30, 16: 3.91, 18: 3.53},\n", + " \"6x6x6\": {12: 3.43, 14: 4.32, 16: 3.85, 18: 3.48},\n", + " },\n", + " \"MATPES-r2SCAN-0\": {\n", + " \"4x4x4\": {12: -0.83, 14: -0.46, 16: -0.49},\n", + " \"5x5x5\": {12: -0.85, 14: -0.32, 16: -0.61, 18: 0.03},\n", + " \"6x6x6\": {12: -0.52, 14: -0.37, 16: -0.64, 18: 0.08},\n", + " },\n", + "}\n", + "IMAGE_LIMIT = {\"4x4x4\": 14.6, \"5x5x5\": 18.3, \"6x6x6\": 21.9} # Bohr\n", + "# finite displacements in the 10x10x10 supercell, from the next section\n", + "ALPHA_FD_CONVERGED = {\n", + " \"OMAT-0-medium\": 2.90,\n", + " \"MATPES-PBE-0\": 3.71,\n", + " \"MATPES-r2SCAN-0\": 0.58,\n", + "}\n", + "\n", + "fig, axes = plt.subplots(1, 3, figsize=(12, 3.6), sharey=True)\n", + "for ax, (model, by_cell) in zip(axes, ALPHA_SI.items(), strict=True):\n", + " for cell, by_cutoff in by_cell.items():\n", + " cutoffs = np.array(list(by_cutoff))\n", + " values = np.array(list(by_cutoff.values()))\n", + " (line,) = ax.plot(cutoffs, values, \"o-\", label=cell)\n", + " beyond = cutoffs > IMAGE_LIMIT[cell]\n", + " ax.plot(\n", + " cutoffs[beyond], values[beyond], \"o\", color=line.get_color(), mfc=\"white\"\n", + " )\n", + " ax.axhline(\n", + " ALPHA_FD_CONVERGED[model], color=\"black\", ls=\":\", label=\"finite displacements\"\n", + " )\n", + " ax.axhline(2.6, color=\"gray\", ls=\"--\", label=\"experiment\")\n", + " ax.set_title(model)\n", + " ax.set_xlabel(\"Third-order cutoff (Bohr)\")\n", + "axes[0].set_ylabel(\"$\\\\alpha$(300 K) (10$^{-6}$ K$^{-1}$)\")\n", + "axes[0].legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "Open markers are cutoffs above the image limit of the supercell. The dotted\n", + "line is the converged finite-displacement value from the next section. The\n", + "dashed line is the experimental value of about 2.6e-6 K⁻¹ (Okada and Tokumaru,\n", + "J. Appl. Phys. 56, 314 (1984)).\n", + "\n", + "- The default 12 Bohr cutoff is far from the converged value. With\n", + " MACE-OMAT-0-medium it gives about half of it. With MATPES-r2SCAN-0 even the\n", + " sign is wrong.\n", + "- Between 12 and 18 Bohr the value goes up and down as each new shell of\n", + " neighbors enters the fit. It does not settle within the cutoffs we could\n", + " afford.\n", + "- From 14 Bohr on, the value at a fixed cutoff changes little between the 5x5x5\n", + " and 6x6x6 supercells. The cutoff limits the result, more than the supercell." + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "## Finite displacements as a reference\n", + "\n", + "To find the converged value, we computed the force constants of Si by finite\n", + "displacements with phono3py, in supercells from 4x4x4 to 10x10x10. There is no\n", + "cutoff and no fit. The 10x10x10 supercell has 2000 atoms and needed 6201\n", + "displaced supercells for each potential." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "# alpha(300 K) in 1e-6/K from finite displacements, by potential and supercell\n", + "ALPHA_FD = {\n", + " \"OMAT-0-medium\": {4: 1.74, 5: 1.96, 6: 2.32, 7: 2.64, 8: 2.83, 9: 2.89, 10: 2.90},\n", + " \"MATPES-PBE-0\": {4: 3.75, 5: 3.41, 6: 3.45, 7: 3.58, 8: 3.67, 9: 3.70, 10: 3.71},\n", + " \"MATPES-r2SCAN-0\": {\n", + " 4: -0.53,\n", + " 5: -0.10,\n", + " 6: 0.14,\n", + " 7: 0.44,\n", + " 8: 0.56,\n", + " 9: 0.57,\n", + " 10: 0.58,\n", + " },\n", + "}\n", + "\n", + "fig, ax = plt.subplots(figsize=(5, 3.6))\n", + "for model, by_size in ALPHA_FD.items():\n", + " ax.plot(list(by_size), list(by_size.values()), \"o-\", label=model)\n", + "ax.axhline(2.6, color=\"gray\", ls=\"--\", label=\"experiment\")\n", + "ax.set_xlabel(\"Supercell size n (n x n x n)\")\n", + "ax.set_ylabel(\"$\\\\alpha$(300 K) (10$^{-6}$ K$^{-1}$)\")\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, + "source": [ + "- From 9x9x9 to 10x10x10 the value changes by less than 0.02e-6 K⁻¹. It is\n", + " converged at 2.90e-6 K⁻¹ with MACE-OMAT-0-medium, 3.71e-6 K⁻¹ with\n", + " MATPES-PBE-0 and 0.58e-6 K⁻¹ with MATPES-r2SCAN-0.\n", + "- We also cut the 8x8x8 finite-displacement force constants at the pheasy\n", + " cutoffs and recomputed the thermal expansion. At 16 and 18 Bohr this agrees\n", + " with the 6x6x6 fits within 0.05e-6 K⁻¹. The fits are right for their cutoff.\n", + "- The value only settles in the 9x9x9 and 10x10x10 supercells. They hold\n", + " interactions up to 17 and 19 Å. Each of these MACE potentials sees the\n", + " neighbors within 12 Å of an atom, through two layers with a 6 Å cutoff. So\n", + " the force constants can couple atoms up to 24 Å apart.\n", + "- A pheasy fit with a cutoff of about 19 Å, or 36 Bohr, is far beyond what the\n", + " workflow can do. The fit at 20 Bohr already ran out of memory.\n", + "- MACE-OMAT-0-medium gives 12% more than experiment and MATPES-PBE-0 gives 43%\n", + " more. MATPES-r2SCAN-0 gives about a fifth of the experimental value.\n", + "\n", + "Si is a hard case, because its thermal expansion is a small difference of two\n", + "larger terms. For such a material, compare the result with finite displacements\n", + "in growing supercells. In section 1 the default supercell is within 2% for the\n", + "rock-salt materials with MACE-OMAT-0-medium. With the other potentials, and for\n", + "GaAs, a larger supercell changes the result by up to 38%." + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, + "source": [ + "## Known limitations\n", + "\n", + "The frequencies and Grüneisen tensors are those of the relaxed structure. They\n", + "are not renormalized with temperature, so the result is least reliable at high\n", + "temperature.\n", + "\n", + "A machine-learned potential gives no Born charges, so the non-analytical term\n", + "correction is not applied here. With VASP, the Born charges of the phonon run\n", + "are used.\n", + "\n", + "If a frequency on the q-point mesh is below -0.1 THz, a warning is raised and\n", + "the thermal expansion of that fit is not computed." + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tutorials/tutorials.md b/tutorials/tutorials.md index e30858a4b9..487c6da5e9 100644 --- a/tutorials/tutorials.md +++ b/tutorials/tutorials.md @@ -18,6 +18,7 @@ pheasy_workflow hiphive_workflow force_fields/phonon_workflow grueneisen_workflow +cte_workflow qha_workflow torchsim_tutorial ```