Nudged Elastic Band

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)

You can toggle the following to investigate different models:

model_params = {"arch": "mace_mp", "model_path": "medium-0b3"}
# model_params = {"arch": "mace_mp", "model_path": "medium-mpa-0"}
# model_params = {"arch": "mace_mp", "model_path": "medium-omat-0"}
# model_params = {"arch": "chgnet"}
# model_params = {"arch": "sevennet"}

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 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:

get_data(
    # url="https://raw.githubusercontent.com/stfc/janus-core/main/tests/data/",
    url="https://raw.githubusercontent.com/stfc/janus-tutorials/main/neb/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