{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# Battery project\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Day 3 - Equilibrium potential\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Today you will study the LiFePO$_4$ cathode. You will calculate the equilibrium potential and use Bayesian error estimation to quantify how sensitive the calculated equilibrium potential is towards choice of functional. After today you should be able to discuss:\n",
        "\n",
        "-  The volume change during charge/discharge.\n",
        "\n",
        "-  The maximum gravimetric and volumetric energy density of a FePO$_4$/C battery assuming the majority of weight and volume will be given by the electrodes.\n",
        "\n",
        "-  Uncertainty in the calculations.\n",
        "\n",
        "Some of calculations you will perform today will be tedious to be run in this notebook. You will automatically submit some calculations to the HPC cluster directly from this notebook. When you have to wait for calculations to finish you can get started on addressing the bullet points above.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Initialize\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%matplotlib notebook\n",
        "import numpy as np\n",
        "from ase.visualize import view\n",
        "import matplotlib.pyplot as plt\n",
        "from ase.io import read, write, Trajectory\n",
        "from ase.parallel import paropen\n",
        "from gpaw import GPAW, FermiDirac, Mixer, PW\n",
        "from ase.dft.bee import BEEFEnsemble\n",
        "from ase import Atoms"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## FePO$_4$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "First we will construct an atoms object for FePO$_4$. ASE can read files from in a large number of different [formats](https://wiki.fysik.dtu.dk/ase/ase/io/io.html?highlight=read%20formats#file-input-and-output). However, in this case you will build it from scratch using the below information:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Positions:\n",
        "# Fe      2.73015081       1.46880951       4.56541172\n",
        "# Fe      2.23941067       4.40642872       2.14957739\n",
        "# Fe      7.20997230       4.40642925       0.26615813\n",
        "# Fe      7.70070740       1.46880983       2.68199421\n",
        "# O       1.16033403       1.46881052       3.40240205\n",
        "# O       3.80867172       4.40642951       0.98654342\n",
        "# O       8.77981469       4.40642875       1.42923946\n",
        "# O       6.13142032       1.46881092       3.84509827\n",
        "# O       4.37288562       1.46880982       0.81812712\n",
        "# O       0.59764596       4.40643021       3.23442747\n",
        "# O       5.56702590       4.40642886       4.01346264\n",
        "# O       9.34268360       1.46880929       1.59716233\n",
        "# O       1.64001691       0.26061277       1.17298291\n",
        "# O       3.32931769       5.61463705       3.58882629\n",
        "# O       8.30013707       3.19826250       3.65857000\n",
        "# O       6.61076951       2.67698811       1.24272700\n",
        "# O       8.30013642       5.61459688       3.65856912\n",
        "# O       6.61076982       0.26063178       1.24272567\n",
        "# O       1.64001666       2.67700652       1.17298270\n",
        "# O       3.32931675       3.19822249       3.58882660\n",
        "# P       0.90585688       1.46880966       1.89272372\n",
        "# P       4.06363530       4.40642949       4.30853266\n",
        "# P       9.03398503       4.40642957       2.93877879\n",
        "# P       5.87676435       1.46881009       0.52297232\n",
        "\n",
        "\n",
        "# Unit cell:\n",
        "#            periodic     x          y          z\n",
        "#   1. axis:    yes    9.94012    0.00000    0.00000\n",
        "#   2. axis:    yes    0.00000    5.87524    0.00000\n",
        "#   3. axis:    yes    0.00000    0.00000    4.83157\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You *can* use the cell below as a starting point.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# fepo4 = Atoms('Fe4O...',\n",
        "#               positions=[[x0, y0, z0],[x1, y1, z1]...],\n",
        "#               cell=[x, y, z],\n",
        "#               pbc=[True, True, True])\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Visualize the structure you have made. Explore the different functions in the visualizer and determine the volume of the cell (`View -> Quick Info`).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(fepo4)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "For better convergence of calculations you should specify initial magnetic moments to iron. The iron will in this structure be Fe$^{3+}$ as it donates two *4s* electrons and one *3d* electron to PO$_4$$^{3-}$. What is the magnetic moment of iron? For simplicity you should assume that FePO$_4$ is ferromagnetic.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Now examine the initial magnetic moments of the system using an [appropriate method](https://wiki.fysik.dtu.dk/ase/ase/atoms.html?highlight=get_initial#list-of-all-methods).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "magmoms = fepo4.xxx()\n",
        "print(magmoms)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Write your atoms object to file.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "write('fepo4.traj', fepo4)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "For this calculation you will use the BEEF-vdW functional developed by [Wellendorff et al.](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.85.235149) Although there are better alternatives for calculating the energy of bulk systems, the BEEF-vdW has a build-in ensemble for error estimation of calculated energies. In the set-up of this calculator you will append relevant keyword values into a dictionary, which is inputted the calculator object.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "params_GPAW = {}"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "To save computational time while keeping the calculations physically sound, the following should be used:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "params_GPAW['mode']        = PW(500)                      # The used plane wave energy cutoff\n",
        "params_GPAW['nbands']      = -40                          # The number on empty bands had the system been spin-paired\n",
        "params_GPAW['kpts']        = {'size': (2, 4, 5),          # The k-point mesh\n",
        "                              'gamma': True}\n",
        "params_GPAW['spinpol']     = True                         # Performing spin polarized calculations\n",
        "params_GPAW['xc']          = 'BEEF-vdW'                   # The used exchange-correlation functional\n",
        "params_GPAW['occupations'] = FermiDirac(width=0.1,        # The smearing\n",
        "                                        fixmagmom=True)   # Total magnetic moment fixed to the initial value\n",
        "params_GPAW['convergence'] = {'eigenstates': 1.0e-4,      # eV^2 / electron\n",
        "                              'energy':      2.0e-4,      # eV / electron\n",
        "                              'density':     1.0e-3}\n",
        "params_GPAW['mixer']       = Mixer(0.1, 5, weight=100.0)  # The mixer used during SCF optimization\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "DFT suffers from a so-called self-interaction error. An electron interacts with the system electron density, to which it contributes itself. The error is most pronounced for highly localized orbitals. [Hubbard U correction](https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/energetics/hubbardu/hubbardu.html) is used to mitigate the self-interaction error of the highly localized *3d*-electrons of Fe. This is done in GPAW using the `setups` keyword.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "params_GPAW['setups']      = {'Fe': ':d,4.3'}             # U=4.3 applied to d orbitals\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Make a GPAW calculator and attach it to the atoms object. Here you will use [get_potential_energy](https://wiki.fysik.dtu.dk/ase/ase/atoms.html#ase.Atoms.get_potential_energy) to start the calculation.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "calc = GPAW(**params_GPAW)\n",
        "fepo4.calc = calc\n",
        "epot_fepo4_cell = fepo4.get_potential_energy()\n",
        "print(epot_fepo4_cell)\n",
        "write('fepo4_out.traj', fepo4)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You will use the ensemble capability of the BEEF-vdW functional. You will need this later so you should write it to file so you do not have to start all over again later. Start by obtaining the required data from the calculator, i.e., the individual energy of each term in the BEEF-vdW functional expansion. Get the energy difference compared to BEEF-vdW for 2000 ensemble functionals.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.dft.bee import BEEFEnsemble\n",
        "\n",
        "ens = BEEFEnsemble(calc)\n",
        "dE = ens.get_ensemble_energies(2000)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Print the energy differences to file. This is not the most efficient way of printing to file but can allow easier subsequent data treatment.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "with paropen('ensemble_fepo4.dat', 'a') as result:\n",
        "    for e in dE:\n",
        "        print(e, file=result)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You now have what you need to make a full script. Make it in the cell below and execute it to make sure the script runs. Once you have made sure the calculation is able to run, stop it by `interupt the kernel`.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# %%writefile 'fepo4.py'\n",
        "from ase.parallel import paropen\n",
        "from ase.io import read, write\n",
        "from ase.dft.bee import BEEFEnsemble\n",
        "from gpaw import GPAW, FermiDirac, Mixer, PW\n",
        "\n",
        "# Read in the structure you made and wrote to file above\n",
        "fepo4 = read('fepo4.traj')\n",
        "\n",
        "params_GPAW = {...}\n",
        "\n",
        "# do calculation ...\n",
        "# BEEF ...\n",
        "# write ensemble_fepo4.dat file ...\n",
        "\n",
        "write('fepo4_out.traj', fepo4)\n",
        "\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Uncomment the `%%writefile` line and execute the cell again and submit the calculation to the HPC cluster. The calculation should take around 10 minutes.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq submit fepo4.py -R 8:1h  # submits the calculation to 8 cores, 1 hour"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Run the below cell to examine the status of your calculation. If no output is returned, the calculation has either finished or failed.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq ls"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Once the calculation begins, you can run the cells below to open the error log and output of the calculation in a new window. This can be done while the calculation is running.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Error log\n",
        "!cat \"$(ls -t fepo4.py.*err | head -1)\""
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Output\n",
        "!cat \"$(ls -t fepo4.py.*out | head -1)\""
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Once the calculation has finished, load in the result. You can skip past this cell and return later.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "try:\n",
        "    fepo4 = read('fepo4_out.traj')\n",
        "    print('Calculation finished')\n",
        "except FileNotFoundError:\n",
        "    print('Calculation has not yet finished')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## LiFePO$_4$\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You will now do similar for LiFePO$_4$. In this case you will load in a template structure called `lifepo4_wo_li.traj` missing only the Li atoms. It is located in the resources folder.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "lifepo4_wo_li = read('lifepo4_wo_li.traj')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Visualize the structure.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(lifepo4_wo_li)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You should now add Li into the structure using the fractional coordinates below:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Li  0    0    0\n",
        "# Li  0    0.5  0\n",
        "# Li  0.5  0.5  0.5\n",
        "# Li  0.5  0    0.5\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Add Li atoms into the structure, e.g., by following the example in [this ASE tutorial](https://wiki.fysik.dtu.dk/ase/gettingstarted/manipulating_atoms/manipulating_atoms.html?highlight=set_cell#manipulating-atoms).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from numpy import identity\n",
        "from ase import Atom\n",
        "\n",
        "cell = lifepo4_wo_li.get_cell()\n",
        "\n",
        "# ...\n",
        "\n",
        "# lifepo4 = lifepo4_wo_li.copy()\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Visualize the structure with added Li.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(lifepo4)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Ensure that the magnetic moments are as they should be, once again assuming ferromagnetism for simplicity.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "At this point you should save your structure by writing it to a trajectory file.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "write('lifepo4.traj', lifepo4)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You should now calculate the potential energy of this sytem using the method and same calculational parameters as for FePO$_4$ above. Make a full script in the cell below similar to what you did above for FePO$_4$ and make sure that it runs.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# %% writefile 'lifepo4.py'\n",
        "from ase.parallel import paropen\n",
        "from ase.io import read, write\n",
        "from ase.dft.bee import BEEFEnsemble\n",
        "from gpaw import GPAW, FermiDirac, Mixer, PW\n",
        "\n",
        "# Read in the structure you made and wrote to file above\n",
        "lifepo4 = read('lifepo4.traj')\n",
        "\n",
        "params_GPAW = {...}\n",
        "\n",
        "# ...\n",
        "# ...\n",
        "# ...\n",
        "\n",
        "# write('lifepo4_out.traj', lifepo4)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "If the code runs, submit to the HPC cluster as you did above. The calculation takes approximately 10 minutes.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq submit lifepo4.py -R 8:1h  # submits the calculation to 8 cores, 1 hour"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Run the below cell to examine the status of your calculation. If no output is returned, the calculation has either finished or failed.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq ls"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Once the calculation begins, you can run the cells below to open the error log and output of the calculation in a new window.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Error log\n",
        "!cat \"$(ls -t lifepo4.py.*err | head -1)\""
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Output\n",
        "!cat \"$(ls -t lifepo4.py.*out | head -1)\""
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "When calculation has finished, load in the result. You can skip past this cell and return later.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "try:\n",
        "    lifepo4=read('lifepo4_out.traj')\n",
        "    print('Calculation finished')\n",
        "except FileNotFoundError:\n",
        "    print('Calculation has not yet finished')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "### Li metal\n",
        "\n",
        "We use a Li metal reference to calculate the equilibrium potential. On exercise day 2 you used a Li metal reference to calculate the intercalation energy in the graphite anode. The approach is similar here. Just read in the result of the calculation with DFTD3 and attach a new calulator to it.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase import Atoms\n",
        "from gpaw import GPAW, FermiDirac, PW\n",
        "from ase.io import read\n",
        "\n",
        "li_metal = read('Li-metal-DFTD3.traj')  # Change file name accordingly\n",
        "\n",
        "calc = GPAW(mode=PW(500),\n",
        "            kpts=(8, 8, 8),\n",
        "            occupations=FermiDirac(0.15),\n",
        "            nbands=-10,\n",
        "            txt=None,\n",
        "            xc='BEEF-vdW')\n",
        "\n",
        "li_metal.calc = calc\n",
        "li_metal.get_potential_energy()\n",
        "li_metal.write('li_metal.traj')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Now calculate the ensemble in the same way as for FePO$_4$ and LiFePO$_4$.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "ens = BEEFEnsemble(calc)\n",
        "li_metal_ens_cell= ens.get_ensemble_energies(2000)\n",
        "with paropen('ensemble_li_metal.dat', 'a') as result:\n",
        "    for e in li_metal_ens_cell:\n",
        "        print(e, file=result)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Calculate equilibrium potential and uncertainty\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can now calculate the equilibrium potential for the case of a FePO$_4$/Li metal battery from the intercallation energy of Li in FePO$_4$. For simplicity, use that assumption that all vibrational energies and entropic terms cancel each other. You should now have completed all submitted calculations before you proceed.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "The calculated energies are for the full cells. Convert them to the energy per formula unit. The [len(...)](https://docs.python.org/3.6/library/functions.html#len) function can be quite helpful for this.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "epot_fepo4_cell=fepo4.get_potential_energy()\n",
        "epot_lifepo4_cell=lifepo4.get_potential_energy()\n",
        "epot_li_metal_cell=li_metal.get_potential_energy()\n",
        "print('epot_fepo4_cell =', epot_fepo4_cell)\n",
        "print('epot_lifepo4_cell =', epot_lifepo4_cell)\n",
        "print('epot_li_metal_cell =', epot_li_metal_cell)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "epot_fepo4 = ...\n",
        "epot_lifepo4 = ...\n",
        "epot_li_metal = ...\n",
        "# print(epot_fepo4, ...)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "No calculate the equilibrium potential under the assumption that it is given by $V_{eq} = \\Delta U /e $, where $U$ is the electronic potential energy of the system and $e$ is the number of electrons transfered.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# V_eq = ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You will now calculate the error estimate for the Li intercallation energy in FePO$_4$ using the BEEF ensemble results. Start by loading in the files. Wait a few minutes and rerun the cell if the number is not 2000 for all of them.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "fepo4_ens_cell = np.genfromtxt('ensemble_fepo4.dat', max_rows=2000)\n",
        "lifepo4_ens_cell = np.genfromtxt('ensemble_lifepo4.dat', max_rows=2000)\n",
        "\n",
        "print('number of functionals in ensemble=', len(fepo4_ens_cell))\n",
        "print('number of functionals in ensemble=', len(lifepo4_ens_cell))\n",
        "print('number of functionals in ensemble=', len(li_metal_ens_cell))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Note that these are energies per cell and not per formula unit. Convert them as you did the potential energies above. Note that you are now performing the operation on a list of length 2000 and not a single float value as before.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# fepo4_ens = fepo4_ens_cell / ...\n",
        "# ...\n",
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Make a list of equilibrium potentials.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# V_eq_ens = lifepo4_ens - ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Use the plot command below to visualize the distribution.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "plt.hist(V_eq_ens, 50)\n",
        "plt.grid(True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Use the [NumPy function standard deviation function](https://docs.scipy.org/doc/numpy/reference/generated/numpy.std.html#numpy-std) to obtain the standard deviation of the ensemble.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# error = ...\n",
        "# print(error)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "The equilibrium potential for a FePO$_4$/Li battery is thus as a good estimate:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "print(f'{V_eq:.2f} V +- {error:.2f} V')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can get the equilibrium potential for the FePO$_4$/C battery using the intercallation energy of Li in graphite, that you calculated on Day 2. What equilibrium potential do you find? How does that compare to the cell voltage you can obtain from FePO$_4$/C batteries?\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# You can use this cell for FePO4/C potential calculation"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Make sure you are able to discuss the bullet points at the top of this notebook. You can use the cell below for calculations.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Bonus\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "How does the predicted error estimate change if you consider the full reaction from Li in graphite + FePO4  to empty graphite + LiFePO4.\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "CAMD2022",
      "language": "python",
      "name": "camd2022"
    },
    "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.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}