Nudged Elastic Band¶

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.

In this tutorial, we will determine the activation energies of Li diffusion along the [010] and [001] directions (referred to as paths b and c respectively) in lithium iron phosphate (LiFePO_4), a cathode material for lithium ion batteries.

DFT references energies are:

Barrier heights:
- path b = 0.27 eV
- path c = 2.5 eV

(see table 1 in https://doi.org/10.1039/C5TA05062F)

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.

[ ]:

import sys

if "google.colab" in sys.modules:
    ! uv pip install janus-core[mace,orb,chgnet,visualise] data-tutorials --system # Install janus-core with MACE, Orb, CHGNet, and WeasWidget, and data-tutorials

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/

[ ]:

from janus_core import __version__

print(__version__)

Prepare data, modules, and model parameters¶

You can toggle the following to investigate different models:

[ ]:

model_params = {"arch": "mace_mp", "model": "medium-0b3", "device": "cuda"}
# model_params = {"arch": "mace_mp", "model": "medium-mpa-0", "device": "cuda"}
# model_params = {"arch": "mace_mp", "model": "medium-omat-0", "device": "cuda"}
# model_params = {"arch": "chgnet"}
# model_params = {"arch": "orb"}

[ ]:

from weas_widget import WeasWidget
from ase.io import read
from data_tutorials.data import get_data

from janus_core.calculations.geom_opt import GeomOpt
from janus_core.calculations.neb import NEB

Use data_tutorials to get the data required for this tutorial – not needed if you are in the provided cloud:

[ ]:

if "google.colab" in sys.modules:
    get_data(
        url="https://raw.githubusercontent.com/ccp5UK/summerschool/main/data/",
        filename="LiFePO4_supercell.cif",
        folder="data",
    )

Preparing end structures¶

The initial structure can be downloaded from the Materials Project (mp-19017):

[ ]:

LFPO = read("data/LiFePO4_supercell.cif")
v=WeasWidget()
v.from_ase(LFPO)
v.avr.model_style = 1
v.avr.show_hydrogen_bonds = True
v

First, we will relax the supercell:

[ ]:

GeomOpt(struct=LFPO, **model_params).run()

v1=WeasWidget()
v1.from_ase(LFPO)
v1.avr.model_style = 1
v1.avr.show_hydrogen_bonds = True
v1

Next, we will create the start and end structures:

[ ]:

# NEB path along b and c directions have the same starting image.
# For start bc remove site 5
LFPO_start_bc = LFPO.copy()
del LFPO_start_bc[5]

# For end b remove site 11
LFPO_end_b = LFPO.copy()
del LFPO_end_b[11]

# For end c remove site 4
LFPO_end_c = LFPO.copy()
del LFPO_end_c[4]

Path b¶

We can now calculate the barrier height along path b.

This also includes running geometry optimization on the end points of this path.

[ ]:

n_images = 7
interpolator="pymatgen" # ASE interpolation performs poorly in this case

neb_b = NEB(
    init_struct=LFPO_start_bc,
    final_struct=LFPO_end_b,
    n_images=n_images,
    interpolator=interpolator,
    minimize=True,
    fmax=0.1,
    **model_params,
)

[ ]:

results = neb_b.run()

The results include the barrier (without any interpolation between highest images) and maximum force at the point in the simulation

[ ]:

print(results)

We can also plot the band:

[ ]:

fig = neb_b.nebtools.plot_band()
v1=WeasWidget()
v1.from_ase(neb_b.nebtools.images)
v1.avr.model_style = 1
v1.avr.show_hydrogen_bonds = True
v1

Path c¶

We can calculate the barrier height along path c similarly

[ ]:

n_images = 7
interpolator="pymatgen"

neb_c = NEB(
    init_struct=LFPO_start_bc,
    final_struct=LFPO_end_c,
    n_images=n_images,
    interpolator=interpolator,
    minimize=True,
    fmax=0.1,
    **model_params,
)

[ ]:

results = neb_c.run()

[ ]:

print(results)

[ ]:

fig = neb_c.nebtools.plot_band()
v2=WeasWidget()
v2.from_ase(neb_c.nebtools.images)
v2.avr.model_style = 1
v2.avr.show_hydrogen_bonds = True
v2

Continuing NEBs¶

You can also continue a NEB to a tighter fmax:

[ ]:

neb_b.optimize(fmax=0.05)
neb_b.run_nebtools()
print(neb_b.results)

fig = neb_b.nebtools.plot_band()
v1=WeasWidget()
v1.from_ase(neb_b.nebtools.images)
v1.avr.model_style = 1
v1.avr.show_hydrogen_bonds = True
v1

It is also possible to change NEB parameters, accessed through the attributes of NEB.neb, such as to continue or retry a NEB using different methods.

Note that although the “success” of the NEB optimisation can be accessed through NEB.converged, this is typically only set to False if a OptimizerConvergenceError, not when the maximum number of steps is reached. We therefore recommend checking NEB.opt.nsteps or NEB.results["max_force"].

[ ]:

neb_c.neb.climb = True
neb_c.run(fmax=0.05)
print(neb_c.results)
print(neb_c.opt.nsteps)

fig = neb_c.nebtools.plot_band()
v1=WeasWidget()
v1.from_ase(neb_c.nebtools.images)
v1.avr.model_style = 1
v1.avr.show_hydrogen_bonds = True
v1