Single Point

janus-core contains various machine learnt interatomic potentials (MLIPs), including MACE based models (MACE-MP, MACE-OFF), CHGNet, SevenNet and more, full list on https://github.com/stfc/janus-core.

Other will be added as their utility is proven beyond a specific material.

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)

Use data_tutorials to get the data required for this tutorial:

from data_tutorials.data import get_data

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

Single point calculations for periodic system

from ase.build import bulk

from janus_core.calculations.single_point import SinglePoint

# Change to cuda if you have a gpu or mps for apple silicon
device = "cpu"

NaCl = bulk("NaCl", "rocksalt", a=5.63, cubic=True)
sp = SinglePoint(
    struct=NaCl,
    arch="mace_mp",
    device=device,
    calc_kwargs={"model_paths": "small", "default_dtype": "float64"},
)
res_mace = sp.run()

NaCl = bulk("NaCl", "rocksalt", a=5.63, cubic=True)
sp = SinglePoint(
    struct=NaCl,
    arch="sevennet",
    device=device,
)
res_sevennet = sp.run()

NaCl = bulk("NaCl", "rocksalt", a=5.63, cubic=True)
sp = SinglePoint(
    struct=NaCl,
    arch="chgnet",
    device=device,
)
res_chgnet = sp.run()

print(f"    MACE[eV]: {res_mace['energy']}")
print(f"SevenNet[eV]: {res_sevennet['energy']}")
print(f"  CHGNeT[eV]: {res_chgnet['energy']}")

Simple Molecules

from ase.build import molecule
from weas_widget import WeasWidget

from janus_core.calculations.single_point import SinglePoint

sp = SinglePoint(
    struct=molecule("H2O"),
    arch="mace_off",
    device=device,
    calc_kwargs={"model_paths": "medium"},
)
res = sp.run()
print(res)
v=WeasWidget()
v.from_ase(sp.struct)
v.avr.model_style = 1
v.avr.show_hydrogen_bonds = True
v

Sugar on salt

from ase.build import add_adsorbate, bulk
from ase.io import read, write

a = 5.63
NaCl = bulk(
    "NaCl", crystalstructure="rocksalt", cubic=True, orthorhombic=True, a=5.63
) * (6, 6, 3)

NaCl.center(vacuum=20.0, axis=2)
sugar = read("data/sucrose.xyz")
add_adsorbate(slab=NaCl, adsorbate=sugar, height=4.5, position=(10, 10))
write("slab.xyz", NaCl)

sp = SinglePoint(
    struct="slab.xyz",
    arch="mace_mp",
    device=device,
    calc_kwargs={"model_paths": "small"},
)
res = sp.run()
print(res)

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

Calculate an entire collection of data frames

from ase.io import read
import matplotlib.pyplot as plt
import numpy as np

frames = read("data/NaCl-set.xyz", format="extxyz", index=":")
dft_energy = np.array([s.info["dft_energy"] / len(s) for s in frames])

sp = SinglePoint(
    struct="data/NaCl-set.xyz",
    arch="mace_mp",
    device=device,
    calc_kwargs={"model_paths": "medium-0b3"},
)
sp.run()

mace_mp_energy = np.array([s.info["mace_mp_energy"] / len(s) for s in sp.struct])
rmse_mace = np.linalg.norm(mace_mp_energy - dft_energy) / np.sqrt(len(dft_energy))

sp = SinglePoint(struct="data/NaCl-set.xyz", arch="chgnet", device=device)
sp.run()

chgnet_energy = np.array([s.info["chgnet_energy"] / len(s) for s in sp.struct])
rmse_chgnet = np.linalg.norm(chgnet_energy - dft_energy) / np.sqrt(len(dft_energy))

sp = SinglePoint(struct="data/NaCl-set.xyz", arch="sevennet", device=device)
sp.run()

sevennet_energy = np.array([s.info["sevennet_energy"] / len(s) for s in sp.struct])
rmse_sevennet = np.linalg.norm(sevennet_energy - dft_energy) / np.sqrt(len(dft_energy))

print(
    f"rmse: mace_mp = {rmse_mace}, chgnet = {rmse_chgnet}, "
    f"sevennet = {rmse_sevennet} eV/atom"
)

fig, ax = plt.subplots()

ax.scatter(dft_energy, mace_mp_energy, marker="o", label="mace-mp-0")
ax.scatter(dft_energy, sevennet_energy, marker="x", label="sevennet")
ax.scatter(dft_energy, chgnet_energy, marker="+", label="chgnet")
ax.legend()
plt.xlabel("MLIP [eV/atom]")
plt.ylabel("DFT [eV/atom]")

plt.show()