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 using various machine learnt interatomic potentials(MLIP).

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)

The foundation potentials (called universal by some) provide an easy way to use MLIPs. In this tutorial we will see how well some of these models perform against each other and against the published DFT literature.

A partial top chart of main MLIPs is available

https://matbench-discovery.materialsproject.org/

References:

Set up environment (optional)¶

These steps are required to run this tutorial with Google Colab. To do so, run the cell below, it may work on other environments but is not tested

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

NOTE 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 locale
 locale.getpreferredencoding = lambda: "UTF-8"

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

! 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/

[ ]:

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"}
# model_params = {"arch": "mace_mp", "model": "medium-mpa-0"}
# model_params = {"arch": "mace_mp", "model": "medium-omat-0"}
# model_params = {"arch": "chgnet"}
# model_params = {"arch": "sevennet"}

[ ]:

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/ccp5UK/summerschool/main/Day_11ML/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