Equation of State

Set up environment (optional)

These steps are required for Google Colab, but may work on other systems too:

# import locale
# locale.getpreferredencoding = lambda: "UTF-8"

# ! pip uninstall torch torchaudio torchvision numpy -y
# ! uv pip install janus-core[all] data-tutorials torch==2.5.1 --system
# get_ipython().kernel.do_shutdown(restart=True)

from ase import Atoms
from ase.io import read
from weas_widget import WeasWidget

from data_tutorials.data import get_data

from janus_core.calculations.eos import EoS

Use data_tutorials to get the data required for this tutorial:

get_data(
    url="https://raw.githubusercontent.com/stfc/janus-tutorials/main/data/",
    filename=["beta_quartz.cif"],
    folder="data",
)

Equation of state for α-quartz

Build the structure:

α_quartz = Atoms(
    symbols=(*["Si"] * 3, *["O"] * 6),
    scaled_positions=[
        [0.469700, 0.000000, 0.000000],
        [0.000000, 0.469700, 0.666667],
        [0.530300, 0.530300, 0.333333],
        [0.413500, 0.266900, 0.119100],
        [0.266900, 0.413500, 0.547567],
        [0.733100, 0.146600, 0.785767],
        [0.586500, 0.853400, 0.214233],
        [0.853400, 0.586500, 0.452433],
        [0.146600, 0.733100, 0.880900],
    ],
    cell=[
        [4.916000, 0.000000, 0.000000],
        [-2.45800, 4.257381, 0.000000],
        [0.000000, 0.000000, 5.405400],
    ],
    pbc=True,
)

v=WeasWidget()
v.from_ase(α_quartz)
v.avr.model_style = 1
v.avr.show_hydrogen_bonds = True
v

Calculate the equation of state using the MACE-MP potential:

mace_eos = EoS(
    struct=α_quartz.copy(),
    arch="mace_mp",
    device="cpu",
    model_path="small",
    calc_kwargs={"default_dtype": "float64"},
    minimize_kwargs={"filter_func": None},
    min_volume=0.75,
    max_volume=1.25,
    n_volumes=20,
).run()

mace_eos["eos"].plot(show=True)

Equation of state for β-quartz

Perform the same calculation for β-quartz:

β_quartz = read("data/beta_quartz.cif")

v=WeasWidget()
v.from_ase(β_quartz)
v.avr.model_style = 1
v.avr.show_hydrogen_bonds = True
v

mace_eos_beta = EoS(
    struct=β_quartz.copy(),
    arch="mace_mp",
    device="cpu",
    model_path="small",
    calc_kwargs={"default_dtype": "float64"},
    minimize_kwargs={"filter_func": None},
    min_volume=0.75,
    max_volume=1.25,
    n_volumes=20,
).run()

mace_eos_beta["eos"].plot(show=True)

Combining plots for α-quartz and β-quartz:

import matplotlib.pyplot as plt

ax = plt.gca()

data_alpha = mace_eos["eos"].getplotdata()
data_beta = mace_eos_beta["eos"].getplotdata()

ax.plot(data_alpha[4], data_alpha[5], ls="-", color="C3", label="α-quartz")
ax.plot(data_alpha[6], data_alpha[7], ls="", marker="x", color="C4", mfc="C4")

ax.plot(data_beta[4], data_beta[5], ls="-", color="C0", label="β-quartz")
ax.plot(data_beta[6], data_beta[7], ls="", marker="x", color="C2", mfc="C2")

ax.set_xlabel("volume [Å$^3$]")
ax.set_ylabel("energy [eV]")
ax.legend()

plt.show()

Comparing MACE to CHGNeT and SevenNet

sevennet_eos = EoS(
    struct=α_quartz.copy(),
    arch="sevennet",
    device="cpu",
    minimize_kwargs={"filter_func": None},
).run()

sevennet_eos["eos"].plot(show=True)

chgnet_eos = EoS(
    struct=α_quartz.copy(),
    arch="chgnet",
    device="cpu",
    minimize_kwargs={"filter_func": None},
).run()

chgnet_eos["eos"].plot(show=True)

print(f"MACE energy [eV]: {mace_eos['e_0']}")
print(f"SevenNet energy [eV]: {sevennet_eos['e_0']}")
print(f"CHGNeT energy [eV]: {chgnet_eos['e_0']}")

print()

print(f"MACE volume [Å^3]: {mace_eos['v_0']}")
print(f"SevenNet volume [Å^3]: {sevennet_eos['v_0']}")
print(f"CHGNeT volume [Å^3]: {chgnet_eos['v_0']}")

print()

print(f"MACE bulk_modulus [GPa]: {mace_eos['bulk_modulus']}")
print(f"SevenNet bulk_modulus [GPa]: {sevennet_eos['bulk_modulus']}")
print(f"CHGNeT bulk_modulus [GPa]: {chgnet_eos['bulk_modulus']}")