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)
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,orb,chgnet,visualise] data-tutorials --system # Install janus-core with MACE, Orb, 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.5
Prepare data, modules, and model parameters¶
You can toggle the following to investigate different models:
[3]:
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": "orb"}
[4]:
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
/home/runner/work/janus-core/janus-core/.venv/lib/python3.12/site-packages/codecarbon/core/gpu.py:4: FutureWarning: The pynvml package is deprecated. Please install nvidia-ml-py instead. If you did not install pynvml directly, please report this to the maintainers of the package that installed pynvml for you.
import pynvml
Use data_tutorials to get the data required for this tutorial:
[5]:
get_data(
url="https://raw.githubusercontent.com/stfc/janus-core/main/docs/source/tutorials/data/",
filename="LiFePO4_supercell.cif",
folder="data",
)
try to download LiFePO4_supercell.cif from https://raw.githubusercontent.com/stfc/janus-core/main/docs/source/tutorials/data/ and save it in data/LiFePO4_supercell.cif
saved in data/LiFePO4_supercell.cif
Preparing end structures¶
The initial structure can be downloaded from the Materials Project (mp-19017):
[6]:
LFPO = read("data/LiFePO4_supercell.cif")
v=WeasWidget()
v.from_ase(LFPO)
v.avr.model_style = 1
v.avr.show_hydrogen_bonds = True
v
[6]:
First, we will relax the supercell:
[7]:
GeomOpt(struct=LFPO, **model_params).run()
v1=WeasWidget()
v1.from_ase(LFPO)
v1.avr.model_style = 1
v1.avr.show_hydrogen_bonds = True
v1
/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.
Downloading MACE model from 'https://github.com/ACEsuit/mace-mp/releases/download/mace_mp_0b3/mace-mp-0b3-medium.model'
Cached MACE model to /home/runner/.cache/mace/macemp0b3mediummodel
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/macemp0b3mediummodel
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:197: 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 16:15:01 -762.842666 0.654854
LBFGS: 1 16:15:03 -763.038742 0.437508
LBFGS: 2 16:15:05 -763.090096 0.392867
LBFGS: 3 16:15:07 -763.119488 0.374220
LBFGS: 4 16:15:09 -763.172830 0.346620
LBFGS: 5 16:15:10 -763.213488 0.334927
LBFGS: 6 16:15:12 -763.256842 0.327520
LBFGS: 7 16:15:14 -763.297792 0.328964
LBFGS: 8 16:15:16 -763.342419 0.332009
LBFGS: 9 16:15:18 -763.385108 0.328513
LBFGS: 10 16:15:20 -763.427187 0.308144
LBFGS: 11 16:15:21 -763.474267 0.274918
LBFGS: 12 16:15:23 -763.528185 0.234548
LBFGS: 13 16:15:25 -763.579927 0.246328
LBFGS: 14 16:15:27 -763.627605 0.245408
LBFGS: 15 16:15:29 -763.684752 0.231582
LBFGS: 16 16:15:31 -763.766360 0.283626
LBFGS: 17 16:15:32 -763.855346 0.296599
LBFGS: 18 16:15:34 -763.931149 0.188692
LBFGS: 19 16:15:36 -763.965778 0.154808
LBFGS: 20 16:15:38 -763.996599 0.213159
LBFGS: 21 16:15:40 -764.039520 0.255107
LBFGS: 22 16:15:41 -764.108589 0.255570
LBFGS: 23 16:15:43 -764.175608 0.200210
LBFGS: 24 16:15:45 -764.208816 0.096659
[7]:
Next, we will create the start and end structures:
[8]:
# 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.
[9]:
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,
)
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/macemp0b3mediummodel
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:197: 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)
[10]:
results = neb_b.run()
Step Time Energy fmax
LBFGS: 0 16:15:47 -758.975432 1.930990
LBFGS: 1 16:15:49 -759.140669 0.711957
LBFGS: 2 16:15:50 -759.209571 0.514600
LBFGS: 3 16:15:52 -759.298948 0.404890
LBFGS: 4 16:15:54 -759.316304 0.265548
LBFGS: 5 16:15:55 -759.336518 0.252446
LBFGS: 6 16:15:57 -759.351338 0.330410
LBFGS: 7 16:15:59 -759.368546 0.309163
LBFGS: 8 16:16:01 -759.378014 0.229741
LBFGS: 9 16:16:02 -759.386280 0.221552
LBFGS: 10 16:16:04 -759.395290 0.278156
LBFGS: 11 16:16:06 -759.404981 0.291370
LBFGS: 12 16:16:07 -759.411852 0.192426
LBFGS: 13 16:16:09 -759.415355 0.075502
Step Time Energy fmax
LBFGS: 0 16:16:11 -758.975436 1.930961
LBFGS: 1 16:16:13 -759.140671 0.711930
LBFGS: 2 16:16:14 -759.209572 0.514582
LBFGS: 3 16:16:16 -759.298949 0.404887
LBFGS: 4 16:16:18 -759.316305 0.265544
LBFGS: 5 16:16:19 -759.336519 0.252446
LBFGS: 6 16:16:21 -759.351338 0.330395
LBFGS: 7 16:16:23 -759.368546 0.309158
LBFGS: 8 16:16:24 -759.378015 0.229745
LBFGS: 9 16:16:26 -759.386281 0.221556
LBFGS: 10 16:16:28 -759.395291 0.278147
LBFGS: 11 16:16:30 -759.404982 0.291359
LBFGS: 12 16:16:31 -759.411852 0.192423
LBFGS: 13 16:16:33 -759.415355 0.075502
Step Time fmax
NEBOptimizer[ode]: 0 16:17:09 1.5090
NEBOptimizer[ode]: 1 16:17:21 1.1032
NEBOptimizer[ode]: 2 16:17:33 0.9983
NEBOptimizer[ode]: 3 16:17:45 0.8704
NEBOptimizer[ode]: 4 16:17:57 0.8013
NEBOptimizer[ode]: 5 16:18:09 0.7712
NEBOptimizer[ode]: 6 16:18:21 0.7441
NEBOptimizer[ode]: 7 16:18:33 0.6738
NEBOptimizer[ode]: 8 16:18:45 0.4346
NEBOptimizer[ode]: 9 16:19:09 0.4299
NEBOptimizer[ode]: 10 16:19:21 0.4226
NEBOptimizer[ode]: 11 16:19:33 0.3946
NEBOptimizer[ode]: 12 16:19:45 0.2843
NEBOptimizer[ode]: 13 16:20:09 0.2726
NEBOptimizer[ode]: 14 16:20:21 0.2652
NEBOptimizer[ode]: 15 16:20:33 0.2612
NEBOptimizer[ode]: 16 16:20:45 0.2552
NEBOptimizer[ode]: 17 16:20:57 0.2323
NEBOptimizer[ode]: 18 16:21:21 0.2219
NEBOptimizer[ode]: 19 16:21:33 0.2141
NEBOptimizer[ode]: 20 16:21:45 0.2087
NEBOptimizer[ode]: 21 16:21:57 0.2055
NEBOptimizer[ode]: 22 16:22:09 0.2005
NEBOptimizer[ode]: 23 16:22:21 0.1809
NEBOptimizer[ode]: 24 16:22:45 0.1729
NEBOptimizer[ode]: 25 16:22:57 0.1664
NEBOptimizer[ode]: 26 16:23:09 0.1624
NEBOptimizer[ode]: 27 16:23:22 0.1598
NEBOptimizer[ode]: 28 16:23:34 0.1559
NEBOptimizer[ode]: 29 16:23:46 0.1409
NEBOptimizer[ode]: 30 16:23:58 0.2612
NEBOptimizer[ode]: 31 16:24:22 0.0826
The results include the barrier (without any interpolation between highest images) and maximum force at the point in the simulation
[11]:
print(results)
{'barrier': np.float64(0.2847173696816294), 'delta_E': np.float64(-1.6499575394846033e-07), 'max_force': np.float64(0.09009922957482004)}
We can also plot the band:
[12]:
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
[12]:
Path c¶
We can calculate the barrier height along path c similarly
[13]:
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,
)
Using Materials Project MACE for MACECalculator with /home/runner/.cache/mace/macemp0b3mediummodel
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:197: 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)
[14]:
results = neb_c.run()
Step Time Energy fmax
LBFGS: 0 16:24:27 -759.415355 0.075502
Step Time Energy fmax
LBFGS: 0 16:24:29 -758.975431 1.930990
LBFGS: 1 16:24:31 -759.140669 0.711958
LBFGS: 2 16:24:33 -759.209571 0.514601
LBFGS: 3 16:24:34 -759.298948 0.404890
LBFGS: 4 16:24:36 -759.316304 0.265548
LBFGS: 5 16:24:38 -759.336518 0.252446
LBFGS: 6 16:24:40 -759.351338 0.330410
LBFGS: 7 16:24:41 -759.368546 0.309163
LBFGS: 8 16:24:43 -759.378014 0.229741
LBFGS: 9 16:24:45 -759.386280 0.221553
LBFGS: 10 16:24:47 -759.395290 0.278156
LBFGS: 11 16:24:48 -759.404981 0.291370
LBFGS: 12 16:24:50 -759.411852 0.192426
LBFGS: 13 16:24:52 -759.415355 0.075502
Step Time fmax
NEBOptimizer[ode]: 0 16:25:28 2.3204
NEBOptimizer[ode]: 1 16:25:40 1.6499
NEBOptimizer[ode]: 2 16:25:52 0.9977
NEBOptimizer[ode]: 3 16:26:04 0.6537
NEBOptimizer[ode]: 4 16:26:28 0.4186
NEBOptimizer[ode]: 5 16:26:40 0.3969
NEBOptimizer[ode]: 6 16:26:53 0.3176
NEBOptimizer[ode]: 7 16:27:17 0.2998
NEBOptimizer[ode]: 8 16:27:41 0.2690
NEBOptimizer[ode]: 9 16:27:53 0.2608
NEBOptimizer[ode]: 10 16:28:05 0.2293
NEBOptimizer[ode]: 11 16:28:29 0.2186
NEBOptimizer[ode]: 12 16:28:41 0.2622
NEBOptimizer[ode]: 13 16:28:52 0.2009
NEBOptimizer[ode]: 14 16:29:04 0.1946
NEBOptimizer[ode]: 15 16:29:16 0.1883
NEBOptimizer[ode]: 16 16:29:28 0.1653
NEBOptimizer[ode]: 17 16:29:52 0.1600
NEBOptimizer[ode]: 18 16:30:04 0.1538
NEBOptimizer[ode]: 19 16:30:16 0.1475
NEBOptimizer[ode]: 20 16:30:29 0.1379
NEBOptimizer[ode]: 21 16:30:41 0.1346
NEBOptimizer[ode]: 22 16:30:53 0.1288
NEBOptimizer[ode]: 23 16:31:05 0.1263
NEBOptimizer[ode]: 24 16:31:17 0.1169
NEBOptimizer[ode]: 25 16:31:29 0.1806
NEBOptimizer[ode]: 26 16:31:53 0.0821
[15]:
print(results)
{'barrier': np.float64(1.787575440772347), 'delta_E': np.float64(4.249386620358564e-09), 'max_force': np.float64(0.09762768267248513)}
[16]:
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
[16]:
