Open In Colab

Equation of State

Set up environment (optional)

These steps are required to run this tutorial with Google Colab. To do so, uncomment and run the cell below.

This will replace pre-installed versions of numpy and torch in Colab with versions that are known to be compatible with janus-core.

It may be possible to skip the steps that uninstall and reinstall torch, which will save a considerable amount of time.

These instructions but may work for other systems too, but it is typically preferable to prepare a virtual environment separately before running this notebook if possible.

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

# ! pip uninstall numpy -y # Uninstall pre-installed numpy

# ! pip uninstall torch torchaudio torchvision transformers -y # Uninstall pre-installed torch
# ! uv pip install torch==2.5.1 # Install pinned version of torch

# ! uv pip install janus-core[mace,sevennet,chgnet,visualise] data-tutorials --system # Install janus-core with MACE, SevenNet, CHGNet, and WeasWidget, and data-tutorials

# get_ipython().kernel.do_shutdown(restart=True) # Restart kernel to update libraries. This may warn that your session has crashed.

To ensure you have the latest version of janus-core installed, compare the output of the following cell to the latest version available at https://pypi.org/project/janus-core/

[2]:
from janus_core import __version__

print(__version__)
0.8.3

Prepare data and modules

[3]:
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:

[4]:
get_data(
    url="https://raw.githubusercontent.com/stfc/janus-core/main/docs/source/tutorials/data/",
    filename=["beta_quartz.cif"],
    folder="data",
)
try to download beta_quartz.cif from https://raw.githubusercontent.com/stfc/janus-core/main/docs/source/tutorials/data/ and save it in data/beta_quartz.cif
saved in data/beta_quartz.cif

Equation of state for α-quartz

Build the structure:

[5]:
α_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
[5]:

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

[6]:
mace_eos = EoS(
    struct=α_quartz.copy(),
    arch="mace_mp",
    device="cpu",
    model="small",
    calc_kwargs={"default_dtype": "float64"},
    minimize_kwargs={"filter_class": None},
    min_volume=0.75,
    max_volume=1.25,
    n_volumes=20,
).run()
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/e3nn/o3/_wigner.py:10: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
  _Jd, _W3j_flat, _W3j_indices = torch.load(os.path.join(os.path.dirname(__file__), 'constants.pt'))
cuequivariance or cuequivariance_torch is not available. Cuequivariance acceleration will be disabled.
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/20231210mace128L0_energy_epoch249model
Using float64 for MACECalculator, which is slower but more accurate. Recommended for geometry optimization.
Using head Default out of ['Default']
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/mace/calculators/mace.py:143: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
  torch.load(f=model_path, map_location=device)
       Step     Time          Energy          fmax
LBFGS:    0 17:00:44      -71.020250        0.406340
LBFGS:    1 17:00:44      -71.032943        0.324324
LBFGS:    2 17:00:44      -71.057296        0.251886
LBFGS:    3 17:00:44      -71.064274        0.188059
LBFGS:    4 17:00:45      -71.066728        0.112945
LBFGS:    5 17:00:45      -71.070293        0.166075
LBFGS:    6 17:00:45      -71.074647        0.187184
LBFGS:    7 17:00:45      -71.078118        0.141462
LBFGS:    8 17:00:45      -71.079346        0.074487
[7]:
mace_eos["eos"].plot(show=True)
../../_images/tutorials_python_eos_15_0.png
[7]:
<Axes: title={'center': 'birchmurnaghan: E: -71.091 eV, V: 114.353 Å$^3$, B: 185.324 GPa'}, xlabel='volume [Å$^3$]', ylabel='energy [eV]'>

Equation of state for β-quartz

Perform the same calculation for β-quartz:

[8]:
β_quartz = read("data/beta_quartz.cif")

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

[8]:
[9]:
mace_eos_beta = EoS(
    struct=β_quartz.copy(),
    arch="mace_mp",
    device="cpu",
    model="small",
    calc_kwargs={"default_dtype": "float64"},
    minimize_kwargs={"filter_class": None},
    min_volume=0.75,
    max_volume=1.25,
    n_volumes=20,
).run()
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/20231210mace128L0_energy_epoch249model
Using float64 for MACECalculator, which is slower but more accurate. Recommended for geometry optimization.
Using head Default out of ['Default']
       Step     Time          Energy          fmax
LBFGS:    0 17:00:46      -70.779216        0.490905
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/mace/calculators/mace.py:143: UserWarning: Environment variable TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD detected, since the`weights_only` argument was not explicitly passed to `torch.load`, forcing weights_only=False.
  torch.load(f=model_path, map_location=device)
LBFGS:    1 17:00:46      -70.798248        0.413527
LBFGS:    2 17:00:46      -70.843094        0.020448
[10]:
mace_eos_beta["eos"].plot(show=True)
../../_images/tutorials_python_eos_20_0.png
[10]:
<Axes: title={'center': 'birchmurnaghan: E: -71.024 eV, V: 124.056 Å$^3$, B: 172.905 GPa'}, xlabel='volume [Å$^3$]', ylabel='energy [eV]'>

Combining plots for α-quartz and β-quartz:

[11]:
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()
../../_images/tutorials_python_eos_22_0.png

Comparing MACE to CHGNeT and SevenNet

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

sevennet_eos["eos"].plot(show=True)
       Step     Time          Energy          fmax
LBFGS:    0 17:00:50      -71.035088        0.385548
LBFGS:    1 17:00:50      -71.046516        0.309889
LBFGS:    2 17:00:50      -71.069382        0.239566
LBFGS:    3 17:00:50      -71.075714        0.174657
LBFGS:    4 17:00:50      -71.077728        0.100454
LBFGS:    5 17:00:50      -71.080536        0.149302
LBFGS:    6 17:00:51      -71.084152        0.173768
LBFGS:    7 17:00:51      -71.087196        0.135815
LBFGS:    8 17:00:51      -71.088348        0.074208
../../_images/tutorials_python_eos_24_1.png
[12]:
<Axes: title={'center': 'birchmurnaghan: E: -71.098 eV, V: 114.510 Å$^3$, B: 183.490 GPa'}, xlabel='volume [Å$^3$]', ylabel='energy [eV]'>
[13]:
chgnet_eos = EoS(
    struct=α_quartz.copy(),
    arch="chgnet",
    device="cpu",
    minimize_kwargs={"filter_class": None},
).run()

chgnet_eos["eos"].plot(show=True)
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/nvidia_smi.py:810: SyntaxWarning: invalid escape sequence '\A'
  mem = 'N\A'
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/nvidia_smi.py:831: SyntaxWarning: invalid escape sequence '\A'
  maxMemoryUsage = 'N\A'
CHGNet v0.3.0 initialized with 412,525 parameters
CHGNet will run on cpu
       Step     Time          Energy          fmax
LBFGS:    0 17:00:52      -74.947891        0.359846
LBFGS:    1 17:00:53      -74.958481        0.298901
LBFGS:    2 17:00:53      -74.981049        0.196990
LBFGS:    3 17:00:53      -74.985588        0.135577
LBFGS:    4 17:00:53      -74.986984        0.078268
../../_images/tutorials_python_eos_25_2.png
[13]:
<Axes: title={'center': 'birchmurnaghan: E: -75.009 eV, V: 115.249 Å$^3$, B: 176.635 GPa'}, xlabel='volume [Å$^3$]', ylabel='energy [eV]'>
[14]:
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']}")
MACE energy [eV]: -71.09060046602642
SevenNet energy [eV]: -71.09814494621779
CHGNeT energy [eV]: -75.00901087194896

MACE volume [Å^3]: 114.35325408752064
SevenNet volume [Å^3]: 114.51022735844657
CHGNeT volume [Å^3]: 115.24904049894886

MACE bulk_modulus [GPa]: 185.3243187627527
SevenNet bulk_modulus [GPa]: 183.4904268046179
CHGNeT bulk_modulus [GPa]: 176.63511301008066