docs/source/tutorials/python/neb.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/stfc/janus-core/blob/main/docs/source/tutorials/python/neb.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nudged Elastic Band"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "DFT references energies are:\n",
    "\n",
    "- Barrier heights:\n",
    "  - path b = 0.27 eV\n",
    "  - path c = 2.5 eV\n",
    "  \n",
    "(see table 1 in https://doi.org/10.1039/C5TA05062F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can toggle the following to investigate different models:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_params = {\"arch\": \"mace_mp\", \"model\": \"medium-0b3\"}\n",
    "# model_params = {\"arch\": \"mace_mp\", \"model\": \"medium-mpa-0\"}\n",
    "# model_params = {\"arch\": \"mace_mp\", \"model\": \"medium-omat-0\"}\n",
    "# model_params = {\"arch\": \"chgnet\"}\n",
    "# model_params = {\"arch\": \"sevennet\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up environment (optional)\n",
    "\n",
    "These steps are required for Google Colab, but may work on other systems too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TF-EiWxyuMc7"
   },
   "outputs": [],
   "source": [
    "# import locale\n",
    "# locale.getpreferredencoding = lambda: \"UTF-8\"\n",
    "\n",
    "# ! pip uninstall torch torchaudio torchvision numpy -y\n",
    "# ! uv pip install janus-core[all] data-tutorials torch==2.5.1 --system\n",
    "# get_ipython().kernel.do_shutdown(restart=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from weas_widget import WeasWidget\n",
    "from ase.io import read\n",
    "from data_tutorials.data import get_data\n",
    "\n",
    "from janus_core.calculations.geom_opt import GeomOpt\n",
    "from janus_core.calculations.neb import NEB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `data_tutorials` to get the data required for this tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_data(\n",
    "    # url=\"https://raw.githubusercontent.com/stfc/janus-core/main/tests/data/\",\n",
    "    url=\"https://raw.githubusercontent.com/stfc/janus-tutorials/main/neb/data/\",\n",
    "    filename=\"LiFePO4_supercell.cif\",\n",
    "    folder=\"data\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "857K7R9Cenca"
   },
   "source": [
    "## Preparing end structures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial structure can be downloaded from the Materials Project (mp-19017):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LFPO = read(\"data/LiFePO4_supercell.cif\")\n",
    "v=WeasWidget()\n",
    "v.from_ase(LFPO)\n",
    "v.avr.model_style = 1\n",
    "v.avr.show_hydrogen_bonds = True\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will relax the supercell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cgch8VQ--AES"
   },
   "outputs": [],
   "source": [
    "GeomOpt(struct=LFPO, **model_params).run()\n",
    "\n",
    "v1=WeasWidget()\n",
    "v1.from_ase(LFPO)\n",
    "v1.avr.model_style = 1\n",
    "v1.avr.show_hydrogen_bonds = True\n",
    "v1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will create the start and end structures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NEB path along b and c directions have the same starting image.\n",
    "# For start bc remove site 5\n",
    "LFPO_start_bc = LFPO.copy()\n",
    "del LFPO_start_bc[5]\n",
    "\n",
    "# For end b remove site 11 \n",
    "LFPO_end_b = LFPO.copy()\n",
    "del LFPO_end_b[11]\n",
    "\n",
    "# For end c remove site 4\n",
    "LFPO_end_c = LFPO.copy()\n",
    "del LFPO_end_c[4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Path b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the barrier height along path b.\n",
    "\n",
    "This also includes running geometry optimization on the end points of this path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_images = 7\n",
    "interpolator=\"pymatgen\" # ASE interpolation performs poorly in this case\n",
    "\n",
    "neb_b = NEB(\n",
    "    init_struct=LFPO_start_bc,\n",
    "    final_struct=LFPO_end_b,\n",
    "    n_images=n_images,\n",
    "    interpolator=interpolator,\n",
    "    minimize=True,\n",
    "    fmax=0.1,\n",
    "    **model_params,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = neb_b.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results include the barrier (without any interpolation between highest images) and maximum force at the point in the simulation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the band:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = neb_b.nebtools.plot_band()\n",
    "v1=WeasWidget()\n",
    "v1.from_ase(neb_b.nebtools.images)\n",
    "v1.avr.model_style = 1\n",
    "v1.avr.show_hydrogen_bonds = True\n",
    "v1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Path c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate the barrier height along path c similarly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_images = 7\n",
    "interpolator=\"pymatgen\"\n",
    "\n",
    "neb_c = NEB(\n",
    "    init_struct=LFPO_start_bc,\n",
    "    final_struct=LFPO_end_c,\n",
    "    n_images=n_images,\n",
    "    interpolator=interpolator,\n",
    "    minimize=True,\n",
    "    fmax=0.1,\n",
    "    **model_params,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = neb_c.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = neb_c.nebtools.plot_band()\n",
    "v2=WeasWidget()\n",
    "v2.from_ase(neb_c.nebtools.images)\n",
    "v2.avr.model_style = 1\n",
    "v2.avr.show_hydrogen_bonds = True\n",
    "v2"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "authorship_tag": "ABX9TyNvtIsPHVgkt0NvUv51T6ZG",
   "gpuType": "T4",
   "include_colab_link": true,
   "private_outputs": true,
   "provenance": []
  },
  "kernelspec": {
   "display_name": "janus",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}