Thermal expansion workflow with pheasy and machine-learned potentials¶
This notebook computes the thermal expansion with the CTEMaker workflow. The
forces come from three MACE potentials, so no DFT is needed.
Section 1 runs five materials with the default settings. MgO is worked through step by step. Section 2 converges the third-order cutoff and the supercell size for Si, where the default settings are not enough.
Background¶
The workflow computes the thermal expansion tensor from mode Grüneisen tensors. The thermal stress of each phonon mode is its heat capacity times its Grüneisen tensor. The strain that relaxes the summed thermal stress is the thermal expansion:
Here \(S\) is the elastic compliance, \(c_{q\nu}\) and \(\gamma_{q\nu}\) are the heat capacity and Grüneisen tensor of each mode, \(N_q\) is the number of q-points and \(V\) is the cell volume.
The flow has four steps:
The structure is converted to the standard primitive cell and relaxed tightly.
The pheasy phonon workflow fits the second- and third-order force constants with LASSO. It uses randomly displaced supercells with 0.03 Å displacements.
The elastic workflow fits the elastic tensor on the same relaxed structure.
phono3py computes the mode Grüneisen tensors on a 12x12x12 q-point mesh, and the thermal expansion follows from the formula above.
Steps 2 and 3 do not depend on each other, so a workflow manager can run them
at the same time. The same workflow runs with VASP as
atomate2.vasp.flows.cte.CTEMaker.
Installation¶
The pheasy, phono3py and ase extras are needed, plus MACE and the
Materials Project client:
pip install 'atomate2[pheasy,phono3py,ase]'
pip install 'mace-torch>=0.3.16' mp-api
The pheasy extra pulls in pheasy and ALM, and the phono3py extra pulls in
phono3py. ALM is compiled from source. The hiPhive tutorial and the atomate2
VASP documentation list what the build needs.
The potentials¶
The notebook uses three MACE foundation models. All three are released under the Academic Software License.
MACE-OMAT-0-medium.
mace_mp(model="medium-omat-0")downloads it once and caches it under~/.cache/mace.MACE-MATPES-PBE-0 and MACE-MATPES-r2SCAN-0. These are MACE-OMAT-0 fine-tuned on the MatPES data set at the PBE and r2SCAN levels. They need
mace-torch>=0.3.10. Download the two files from themace_matpes_0release of mace-foundations:
Set the paths to the two downloaded files in the next cell. The MgO example
below only needs MACE-OMAT-0-medium. Change device to "cuda" to run MACE on a
GPU.
float64 is worth the cost here. The third-order force constants come from
small force differences between displaced supercells, where float32 noise is
not negligible.
import os
import warnings
# macOS: conda's llvm-openmp and torch's bundled libomp both load and the
# duplicate aborts the process. This must be set before torch is imported.
os.environ.setdefault("KMP_DUPLICATE_LIB_OK", "TRUE")
warnings.filterwarnings("ignore")
from mace.calculators import mace_mp # noqa: E402
# Replace the two paths with the files you downloaded.
MODEL_FILES = {
"OMAT-0-medium": "medium-omat-0",
"MATPES-PBE-0": "/path/to/MACE-matpes-pbe-omat-ft.model",
"MATPES-r2SCAN-0": "/path/to/MACE-matpes-r2scan-omat-ft.model",
}
def calc_kwargs(model: str) -> dict:
"""Return the calculator settings for one of the models above."""
return {"model": MODEL_FILES[model], "device": "cpu", "default_dtype": "float64"}
_ = mace_mp(**calc_kwargs("OMAT-0-medium")) # downloads and caches on first use
One calculator for all jobs¶
Every force field job builds its own calculator. run_locally runs all of
them in this one Python process, and the workflow has a few hundred jobs.
Loading MACE for each of them grows the memory by about 150 MB per job. This
cell makes the jobs share one calculator. It is only needed when many jobs run
in one process, not with a workflow manager.
from ase.calculators.calculator import Calculator
import atomate2.forcefields.utils as ff_utils
_ase_calculator = ff_utils.ase_calculator
_calculators: dict[str, Calculator] = {}
def shared_ase_calculator(calculator_meta: object, **kwargs: object) -> Calculator:
"""Build each calculator once and reuse it for every job."""
key = repr((calculator_meta, sorted(kwargs.items())))
if key not in _calculators:
_calculators[key] = _ase_calculator(calculator_meta, **kwargs)
_calculators[key].reset()
return _calculators[key]
ff_utils.ase_calculator = shared_ase_calculator
1. Five materials with the default settings¶
The first part runs MgO with MACE-OMAT-0-medium step by step. The other four materials and the other two potentials follow with the same code.
The structure¶
The primitive cell of MgO (mp-1265) comes from the Materials Project. Set your
key first with export MP_API_KEY=....
The Materials Project serves this cell rotated away from the cubic axes, and
the pheasy fits fail for such a cell. CTEMaker therefore converts the input to
the standard primitive cell before the relaxation. This is the default,
use_symmetrized_structure="primitive".
from mp_api.client import MPRester
if not os.environ.get("MP_API_KEY"):
raise OSError("export MP_API_KEY before running this cell")
with MPRester() as mpr:
structure = mpr.get_structure_by_material_id("mp-1265")
structure.lattice
Building the workflow¶
CTEMaker.from_force_field_name sets one force field for the relaxation, the
phonons and the elastic tensor. Any MACE model routes through mace_mp, and
the model name in calculator_kwargs selects the weights.
Everything else stays at the workflow defaults. The supercells have
min_length=12 Å, which gives 128 atoms for MgO. The third-order force
constants use the one-shot fit, and the temperatures run from 0 K to 1000 K in
steps of 10 K.
from atomate2.forcefields.flows.cte import CTEMaker
maker = CTEMaker.from_force_field_name(
"MACE-MP-0", calculator_kwargs=calc_kwargs("OMAT-0-medium")
)
flow = maker.make(structure)
flow.draw_graph().show()
Running the workflow¶
The flow has a little over 200 jobs, most of them force calculations on the displaced supercells. On 32 CPU cores of a Perlmutter node the whole run took about 8 minutes.
create_folders=True is needed, because the thermal expansion job reads the
force constants from the folder of the pheasy fit.
from jobflow import JobStore, run_locally
from maggma.stores import MemoryStore
job_store = JobStore(MemoryStore(), additional_stores={"data": MemoryStore()})
responses = run_locally(flow, store=job_store, create_folders=True, ensure_success=True)
doc = responses[flow.output.uuid][1].output
Results¶
The output is a CTEDocument. It holds the thermal expansion tensor at each
temperature for each fit method. The value at 0 K is zero, since the heat
capacity is zero there. MgO is cubic, so the three diagonal components are
equal.
import matplotlib.pyplot as plt
import numpy as np
(result,) = doc.results
temperatures = np.array(doc.temperatures)
alpha = np.array(result.thermal_expansion_tensor)
i300 = int(np.argmin(np.abs(temperatures - 300)))
fig, ax = plt.subplots()
ax.plot(temperatures, alpha[:, 0, 0] * 1e6)
ax.set_xlabel("Temperature (K)")
ax.set_ylabel("Linear thermal expansion (10$^{-6}$ K$^{-1}$)")
plt.show()
{
"alpha(300 K) in 1/K": float(alpha[i300, 0, 0]),
"mean Grueneisen parameter at 300 K": float(
np.trace(result.average_gruneisen[i300]) / 3
),
"has imaginary modes": result.has_imaginary_modes,
}
With these settings we get \(\alpha\)(300 K) = 11.2e-6 K⁻¹. Experiment gives about 1.0e-5 K⁻¹ for MgO at room temperature.
Checking the fit¶
If the thermal expansion comes out as zero or very small, check the third-order force constants first. All zeros means the LASSO fit dropped them. The penalty chosen by cross validation should also lie inside the search range, not on one of its bounds.
import re
from pathlib import Path
import h5py
fit_dir = Path(doc.phonon_job_dir) / "one_shot"
log = (fit_dir / "pheasy_anharmonic_fit.log").read_text()
with h5py.File(fit_dir / "fc3.hdf5") as file:
fc3 = file["fc3"][:]
{
"LASSO penalty": re.findall(r"alpha_(?:min|max|opt):\s*\S+", log),
"max |fc3| in eV/A^3": float(np.abs(fc3).max()),
}
The other materials and potentials¶
The same workflow runs for NaCl, KCl, CaO and GaAs, and with all three potentials. The function below builds and runs one workflow in its own folder. Section 2 uses it as well.
The loop runs 15 workflows, so it is off by default. Set RUN_ALL = True to
run it.
from pymatgen.core import Structure
from atomate2.common.schemas.cte import CTEDocument
MATERIALS = {
"NaCl": "mp-22862",
"KCl": "mp-23193",
"MgO": "mp-1265",
"CaO": "mp-2605",
"GaAs": "mp-2534",
}
def run_cte(
structure: Structure,
model: str,
root_dir: str,
min_length: float | None = None,
c3: float | None = None,
) -> CTEDocument:
"""Run the thermal expansion workflow and return the CTEDocument.
min_length sets the supercell size in Angstrom. c3 sets the third-order
cutoff in Bohr. None keeps the workflow default.
"""
maker = CTEMaker.from_force_field_name(
"MACE-MP-0", calculator_kwargs=calc_kwargs(model)
)
if min_length is not None:
maker.phonon_maker.min_length = min_length
if c3 is not None:
maker.phonon_maker.fcs_cutoff_radius = [-1, c3, 10]
flow = maker.make(structure)
os.makedirs(root_dir, exist_ok=True)
responses = run_locally(
flow,
store=JobStore(MemoryStore(), additional_stores={"data": MemoryStore()}),
create_folders=True,
root_dir=root_dir,
ensure_success=True,
)
return responses[flow.output.uuid][1].output
def alpha_300k(doc: CTEDocument) -> float:
"""Linear thermal expansion at 300 K in 1/K, from the xx component."""
(result,) = doc.results
i300 = int(np.argmin(np.abs(np.array(doc.temperatures) - 300)))
return float(np.array(result.thermal_expansion_tensor)[i300, 0, 0])
RUN_ALL = False
if RUN_ALL:
with MPRester() as mpr:
structures = {
name: mpr.get_structure_by_material_id(mp_id)
for name, mp_id in MATERIALS.items()
}
alpha_all = {
(name, model): alpha_300k(run_cte(structure, model, f"runs/{name}_{model}"))
for name, structure in structures.items()
for model in MODEL_FILES
}
import pandas as pd
# alpha(300 K) in 1e-6/K from our runs with the default settings
ALPHA_DEFAULT = {
"OMAT-0-medium": {"NaCl": 39.4, "KCl": 41.6, "MgO": 11.2, "CaO": 13.9, "GaAs": 5.9},
"MATPES-PBE-0": {"NaCl": 44.4, "KCl": 47.5, "MgO": 31.4, "CaO": 14.3, "GaAs": 4.4},
"MATPES-r2SCAN-0": {
"NaCl": 32.8,
"KCl": 17.7,
"MgO": 11.3,
"CaO": 11.1,
"GaAs": 2.5,
},
}
# finite displacements in the largest supercell we ran, 686 to 2662 atoms
ALPHA_FD_LARGE = {
"OMAT-0-medium": {"NaCl": 39.3, "KCl": 41.8, "MgO": 11.0, "CaO": 14.0, "GaAs": 6.7},
"MATPES-PBE-0": {"NaCl": 40.0, "KCl": 44.4, "MgO": 39.2, "CaO": 15.8, "GaAs": 6.7},
"MATPES-r2SCAN-0": {
"NaCl": 33.7,
"KCl": 19.4,
"MgO": 12.3,
"CaO": 11.8,
"GaAs": 4.1,
},
}
pd.concat(
{
"default settings": pd.DataFrame(ALPHA_DEFAULT),
"finite displacements": pd.DataFrame(ALPHA_FD_LARGE),
},
axis=1,
)
The first set of columns is what the workflow gives with the default settings. To check it, we computed the force constants of each run by finite displacements with phono3py, with no cutoff and no fit. In the supercell of the run the fits agree with finite displacements within 6%. The one exception is GaAs with MATPES-r2SCAN-0, at 10%. So the fit is right for the supercell it uses.
We then repeated the finite displacements in larger supercells. The second set of columns gives the value in the largest one, with 686 to 2662 atoms. Between the two largest supercells each value changes by about 3% or less.
With MACE-OMAT-0-medium the default settings are within 2% of the large-supercell value for the four rock-salt materials.
With the two MatPES potentials the default settings are off by up to 11% for the rock-salt materials, and by 20% for MgO with MATPES-PBE-0.
For GaAs the default settings are 11% to 38% too low with all three potentials. GaAs has the zinc-blende structure of Si. Section 2 shows the same effect for Si in more detail.
The potentials differ much more than these errors. Two results stand out.
KCl with MATPES-r2SCAN-0 gives 19.4e-6 K⁻¹ in the large supercell. This is less than half the value of the other two potentials.
MgO with MATPES-PBE-0 gives 39.2e-6 K⁻¹. This is more than three times the value of the other two potentials and of experiment. MATPES-PBE-0 also gives a much softer MgO. Its bulk modulus is 100 GPa, against 154 GPa with MACE-OMAT-0-medium and 163 GPa with MATPES-r2SCAN-0. A softer lattice expands more.
2. Converging the cutoff and the supercell for Si¶
In Si the transverse acoustic modes near the zone boundary have negative Grüneisen parameters. Their contribution to the thermal stress partly cancels that of the other modes. The thermal expansion is the small difference of two larger terms. Small errors in the third-order force constants therefore change it a lot. This is why the default settings are not enough for Si.
Two settings control the third-order force constants:
phonon_maker.min_lengthsets the supercell. For the Si primitive cell, 12, 16 and 21 Å give the 4x4x4, 5x5x5 and 6x6x6 supercells with 128, 250 and 432 atoms. The default is 12 Å.phonon_maker.fcs_cutoff_radiussets the cutoff radius of each order in Bohr. The second entry is the third-order cutoff. The default is[-1, 12, 10], so 12 Bohr.
A cutoff should stay below half the shortest distance between periodic images of the supercell. Above it an atom starts to interact with its own images. For Si this limit is 14.6, 18.3 and 21.9 Bohr for the 4x4x4, 5x5x5 and 6x6x6 supercells.
A larger cutoff adds third-order force constants. The workflow picks the number of displaced supercells so that the fit has about 100 equations per free force constant, up to 600 supercells. The cost therefore grows quickly with the cutoff. The 6x6x6 run at 18 Bohr used 201 displaced supercells. At 20 Bohr it used 527, and the LASSO fit needed more than the 220 GB of memory we gave it. We stopped at 18 Bohr.
The scan below runs 11 settings for each potential, so it is off by default.
SI_SCAN = { # min_length in Angstrom: third-order cutoffs in Bohr
12: (12, 14, 16),
16: (12, 14, 16, 18),
21: (12, 14, 16, 18),
}
RUN_SI_SCAN = False
if RUN_SI_SCAN:
with MPRester() as mpr:
si = mpr.get_structure_by_material_id("mp-149")
alpha_si = {
(model, min_length, c3): alpha_300k(
run_cte(
si, model, f"runs/Si_{model}_{min_length}A_{c3}bohr", min_length, c3
)
)
for model in MODEL_FILES
for min_length, cutoffs in SI_SCAN.items()
for c3 in cutoffs
}
# alpha(300 K) in 1e-6/K from our runs, by potential, supercell and cutoff in Bohr
ALPHA_SI = {
"OMAT-0-medium": {
"4x4x4": {12: 1.42, 14: 1.90, 16: 1.87},
"5x5x5": {12: 1.53, 14: 2.04, 16: 2.27, 18: 2.06},
"6x6x6": {12: 1.64, 14: 1.99, 16: 2.30, 18: 2.07},
},
"MATPES-PBE-0": {
"4x4x4": {12: 3.25, 14: 4.13, 16: 4.06},
"5x5x5": {12: 3.35, 14: 4.30, 16: 3.91, 18: 3.53},
"6x6x6": {12: 3.43, 14: 4.32, 16: 3.85, 18: 3.48},
},
"MATPES-r2SCAN-0": {
"4x4x4": {12: -0.83, 14: -0.46, 16: -0.49},
"5x5x5": {12: -0.85, 14: -0.32, 16: -0.61, 18: 0.03},
"6x6x6": {12: -0.52, 14: -0.37, 16: -0.64, 18: 0.08},
},
}
IMAGE_LIMIT = {"4x4x4": 14.6, "5x5x5": 18.3, "6x6x6": 21.9} # Bohr
# finite displacements in the 10x10x10 supercell, from the next section
ALPHA_FD_CONVERGED = {
"OMAT-0-medium": 2.90,
"MATPES-PBE-0": 3.71,
"MATPES-r2SCAN-0": 0.58,
}
fig, axes = plt.subplots(1, 3, figsize=(12, 3.6), sharey=True)
for ax, (model, by_cell) in zip(axes, ALPHA_SI.items(), strict=True):
for cell, by_cutoff in by_cell.items():
cutoffs = np.array(list(by_cutoff))
values = np.array(list(by_cutoff.values()))
(line,) = ax.plot(cutoffs, values, "o-", label=cell)
beyond = cutoffs > IMAGE_LIMIT[cell]
ax.plot(
cutoffs[beyond], values[beyond], "o", color=line.get_color(), mfc="white"
)
ax.axhline(
ALPHA_FD_CONVERGED[model], color="black", ls=":", label="finite displacements"
)
ax.axhline(2.6, color="gray", ls="--", label="experiment")
ax.set_title(model)
ax.set_xlabel("Third-order cutoff (Bohr)")
axes[0].set_ylabel("$\\alpha$(300 K) (10$^{-6}$ K$^{-1}$)")
axes[0].legend()
plt.show()
Open markers are cutoffs above the image limit of the supercell. The dotted line is the converged finite-displacement value from the next section. The dashed line is the experimental value of about 2.6e-6 K⁻¹ (Okada and Tokumaru, J. Appl. Phys. 56, 314 (1984)).
The default 12 Bohr cutoff is far from the converged value. With MACE-OMAT-0-medium it gives about half of it. With MATPES-r2SCAN-0 even the sign is wrong.
Between 12 and 18 Bohr the value goes up and down as each new shell of neighbors enters the fit. It does not settle within the cutoffs we could afford.
From 14 Bohr on, the value at a fixed cutoff changes little between the 5x5x5 and 6x6x6 supercells. The cutoff limits the result, more than the supercell.
Finite displacements as a reference¶
To find the converged value, we computed the force constants of Si by finite displacements with phono3py, in supercells from 4x4x4 to 10x10x10. There is no cutoff and no fit. The 10x10x10 supercell has 2000 atoms and needed 6201 displaced supercells for each potential.
# alpha(300 K) in 1e-6/K from finite displacements, by potential and supercell
ALPHA_FD = {
"OMAT-0-medium": {4: 1.74, 5: 1.96, 6: 2.32, 7: 2.64, 8: 2.83, 9: 2.89, 10: 2.90},
"MATPES-PBE-0": {4: 3.75, 5: 3.41, 6: 3.45, 7: 3.58, 8: 3.67, 9: 3.70, 10: 3.71},
"MATPES-r2SCAN-0": {
4: -0.53,
5: -0.10,
6: 0.14,
7: 0.44,
8: 0.56,
9: 0.57,
10: 0.58,
},
}
fig, ax = plt.subplots(figsize=(5, 3.6))
for model, by_size in ALPHA_FD.items():
ax.plot(list(by_size), list(by_size.values()), "o-", label=model)
ax.axhline(2.6, color="gray", ls="--", label="experiment")
ax.set_xlabel("Supercell size n (n x n x n)")
ax.set_ylabel("$\\alpha$(300 K) (10$^{-6}$ K$^{-1}$)")
ax.legend()
plt.show()
From 9x9x9 to 10x10x10 the value changes by less than 0.02e-6 K⁻¹. It is converged at 2.90e-6 K⁻¹ with MACE-OMAT-0-medium, 3.71e-6 K⁻¹ with MATPES-PBE-0 and 0.58e-6 K⁻¹ with MATPES-r2SCAN-0.
We also cut the 8x8x8 finite-displacement force constants at the pheasy cutoffs and recomputed the thermal expansion. At 16 and 18 Bohr this agrees with the 6x6x6 fits within 0.05e-6 K⁻¹. The fits are right for their cutoff.
The value only settles in the 9x9x9 and 10x10x10 supercells. They hold interactions up to 17 and 19 Å. Each of these MACE potentials sees the neighbors within 12 Å of an atom, through two layers with a 6 Å cutoff. So the force constants can couple atoms up to 24 Å apart.
A pheasy fit with a cutoff of about 19 Å, or 36 Bohr, is far beyond what the workflow can do. The fit at 20 Bohr already ran out of memory.
MACE-OMAT-0-medium gives 12% more than experiment and MATPES-PBE-0 gives 43% more. MATPES-r2SCAN-0 gives about a fifth of the experimental value.
Si is a hard case, because its thermal expansion is a small difference of two larger terms. For such a material, compare the result with finite displacements in growing supercells. In section 1 the default supercell is within 2% for the rock-salt materials with MACE-OMAT-0-medium. With the other potentials, and for GaAs, a larger supercell changes the result by up to 38%.
Known limitations¶
The frequencies and Grüneisen tensors are those of the relaxed structure. They are not renormalized with temperature, so the result is least reliable at high temperature.
A machine-learned potential gives no Born charges, so the non-analytical term correction is not applied here. With VASP, the Born charges of the phonon run are used.
If a frequency on the q-point mesh is below -0.1 THz, a warning is raised and the thermal expansion of that fit is not computed.