{
  "cells": [
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Today you will calculate the energy barriers for transport of Li intercalated in the graphite anode. You will examine how sensitive this barrier is to the interlayer distance in graphite.  You will also examine the energy of intermediate states during the charge/discharge process. This will allow some basic discussion of the voltage profile of the battery.\n",
        "\n",
        "You will in general be provided less code than yesterday, especially towards the end of this notebook. You will have to use what you have already seen and learned so far.\n",
        "\n",
        "There will be some natural pauses while you wait for calculations to finish. If you do not finish this entire notebook today, do not despair.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Initialize\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%matplotlib notebook\n",
        "from ase import Atom\n",
        "from ase.visualize import view\n",
        "import matplotlib.pyplot as plt\n",
        "from ase.io import read, write\n",
        "from ase.mep import NEB\n",
        "from ase.optimize import BFGS\n",
        "from ase.parallel import paropen\n",
        "from gpaw import GPAW, FermiDirac, Mixer, PW\n",
        "from ase.constraints import FixAtoms"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Transport barrier of Li in graphite\n",
        "\"\"\"\n",
        "# %%\n",
        "\"\"\"\n",
        "You will now calculate the energy barrier for Li diffusion in the graphite anode. You will do this using the [Nudged Elastic Band (NEB) method](https://wiki.fysik.dtu.dk/ase/ase/neb.html#module-ase.neb)\n",
        "\n",
        "You can use your work from Day 2, but for simplicity you are advised to load in the initial atomic configuration from file.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "initial = read('NEB_init.traj')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Visualize the structure.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(initial)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You will now make a final structure, where the Li atom has been moved to a neighbouring equivalent site. The [`get_positions`](https://wiki.fysik.dtu.dk/ase/ase/atoms.html?highlight=get_positions#ase.Atoms.get_positions), [`set_positions`](https://wiki.fysik.dtu.dk/ase/ase/atoms.html?highlight=get_positions#ase.Atoms.set_positions) and [`get_cell`](https://wiki.fysik.dtu.dk/ase/ase/atoms.html?highlight=get_positions#ase.Atoms.get_cell) functions are highly useful for such a task. HINT: Displace the Li atom $\\frac{1}{n} (\\vec{a}+\\vec{b})$\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "final = initial.copy()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# ...\n",
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Visualize that you have made the final strcuture correctly.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(final)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Make a band consisting of 7 images including the initial and final.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "images = [initial]\n",
        "images += [initial.copy() for i in range(5)]  # These will become the minimum energy path images.\n",
        "images += [final]"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "It this point `images` consist of 6 copies of `initial` and one entry of `final`. Use the `NEB` method to create an initial guess for the minimum energy path (MEP). In the cell below a simple interpolation between the `initial` and `final` image is used as initial guess.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "neb = NEB(images)\n",
        "neb.interpolate()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Visualize the NEB images.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(images)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "It turns out, that while running the NEB calculation, the largest amount of resources will be spend translating the carbon layer without any noticeable buckling. You will thus [constrain](https://wiki.fysik.dtu.dk/ase/ase/constraints.html#constraints) the positions of the carbon atoms to save computational time.\n",
        "\n",
        "Each image in the NEB requires a unique calculator.\n",
        "\n",
        "This very simple case is highly symmetric. To better illustrate how the NEB method works, the symmetry is broken using the [rattle](https://wiki.fysik.dtu.dk/ase/ase/atoms.html#ase.Atoms.rattle) function.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "for image in images[0:7]:\n",
        "    calc = GPAW(mode=PW(500), kpts=(5, 5, 6), xc='LDA', txt=None, symmetry={'point_group': False})\n",
        "    image.calc = calc\n",
        "    image.set_constraint(FixAtoms(mask=[atom.symbol == 'C' for atom in image]))\n",
        "\n",
        "images[3].rattle(stdev=0.05, seed=42)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Start by calculating the energy and forces of the first (`initial`) and last (`final`) images as this is not done during the actual NEB calculation.\n",
        "\n",
        "Note, that this can take a while if you opt to do it inside the notebook.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "images[0].get_potential_energy()\n",
        "images[0].get_forces()\n",
        "images[6].get_potential_energy()\n",
        "images[6].get_forces()\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can run the NEB calculation by running an optimization on the NEB object the same way you would on an atoms object. Note the `fmax` is larger for this tutorial example than you would normally use.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "optimizer = BFGS(neb, trajectory='neb.traj', logfile='neb.log' )\n",
        "optimizer.run(fmax=0.10)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Submit the calculation to the HPC cluster. Do this by first building a complete script in the cell below using the cells above (minus the `view()` commands). Make sure the cell runs and then interrupt the kernel.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#from ase.io import read, write\n",
        "#from ase.mep import NEB\n",
        "#from ase.optimize import BFGS\n",
        "#from ase.parallel import paropen\n",
        "#from gpaw import GPAW, FermiDirac, Mixer, PW\n",
        "#from ase.constraints import FixAtoms\n",
        "\n",
        "# initial = read('NEB_init.traj')\n",
        "\n",
        "# final = ...\n",
        "\n",
        "# ...\n",
        "# ...\n",
        "\n",
        "# optimizer.run(fmax=0.10)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can use the cell below to submit the calculation in the same way as on earlier days.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq submit NEB.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.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq ls"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "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": [
        "!cat \"$(ls -t NEB.py.*err | head -1)\""
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!cat \"$(ls -t NEB.py.*out | head -1)\""
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "The optimiziation progress can be seen by running the below cell.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!cat neb.log"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can move on while you wait for the calculation to finish.\n",
        "\n",
        "Once the maximum force (`fmax`) in the log is below 0.1, the calculation is finished.\n",
        "Load in the full trajectory.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "full_neb = read('neb.traj@:')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You will use the `ase gui` to inspect the result. The below line reads in the last 7 images in the file. In this case the MEP images.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!ase gui neb.traj@-7:"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "In the GUI use `Tools` $\\rightarrow$ `NEB`.\n",
        "\n",
        "Now inspect how the TS image has developed.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!ase gui neb.traj@3::7"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "For more complicated MEP's, use the [climbing image method](https://wiki.fysik.dtu.dk/ase/ase/neb.html?highlight=neb#climbing-image) to determine the transition state. Why is it not required here?\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Bonus\n",
        "\n",
        "You will now study the influence of changing the interlayer graphite distance on the energy barrier. Due to the high degree of symmetry, this can be done easily in this case. Load in the initial state (IS) and transition state (TS) images from the converged MEP.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "IS_image = images[0]\n",
        "TS_image = images[3]"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Now calculate the energy of the initial state (IS) image and the transition state (TS) image using [`get_potential_energy()`](https://wiki.fysik.dtu.dk/ase/ase/atoms.html?highlight=get_potential_energy#ase.Atoms.get_potential_energy)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "epot_IS = IS_image.get_potential_energy()\n",
        "#epot_TS= ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "barrier = epot_TS - epot_IS\n",
        "print('Energy barrier:', barrier)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Why does this not fully align with what you found before? New reduce the graphite layer distance by change the the size of the unit cell in the *z* direction by 3 %.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "cell = IS_image.get_cell()\n",
        "IS_image97 = IS_image.copy()\n",
        "IS_image97.set_cell([cell[0], cell[1], cell[2] * 0.97], scale_atoms=True)\n",
        "TS_image97 = TS_image.copy()\n",
        "TS_image97.set_cell([cell[0], cell[1], cell[2] * 0.97], scale_atoms=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Use the same calculator object as you did above and calculate the potential energy of the compressed initial and final state.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# calc = ...\n",
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Now calculate the energy of the compressed IS and TS.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# epot_TS97 = ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "What is the energy barrier now?\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# barrier97 = ...\n",
        "# print('Energy barrier:', barrier97)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Now repeat the procedure but expanding the intergraphite distance by 3 %.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# IS_image103 = IS_image.copy()\n",
        "# IS_image103.set_cell(...\n",
        "\n",
        "# calc ...\n",
        "\n",
        "\n",
        "# epot_TS103 = ...\n",
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "What is the energy barrier now?\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# barrier103 = ...\n",
        "# print('Energy barrier:', barrier103)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## FePO$_4$ with one Li\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You will now calculate the energy gain of adding a single Li atom into the FePO$_4$ cell you made on Day 3. This corresponds to a charge of 25 %. You can compare this energy to the equilibrium potential.\n",
        "\n",
        "Load in the FePO$_4$ structure you wrote to file on in a previous exercise and add Li. Assume that the cell dimension remain unchanged.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "#fepo4=read('fepo4.traj')\n",
        "#fepo4_1li=fepo4.copy()\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# fepo4_1li.append(...)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Visualize the structure you made.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(fepo4_1li)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Adjust the total magnetic moment of the cell such that it is 19.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "for atom in fepo4_1li:\n",
        "    if atom.symbol == 'Fe':\n",
        "        atom.magmom = 4.75\n",
        "\n",
        "print(sum(fepo4_1li.get_initial_magnetic_moments()))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Write your atoms object to file giving it the name `fepo4_1li.traj`.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "write('fepo4_1li.traj', fepo4_1li)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Make a full script in the cell below similar to those you made yesterday. Make sure the cell runs before interupting the notebook kernel.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# %%writefile 'fepo4_1li.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_1li = read('fepo4_1li.traj')\n",
        "\n",
        "#...\n",
        "#...\n",
        "\n",
        "# write('fepo4_1li_out.traj', fepo4_1li)\n",
        "\n",
        "# ens = BEEFEnsemble(calc)\n",
        "# with paropen('ensemble_fepo4_1li.dat', 'a') as result:\n",
        "#     for e in dE:\n",
        "#         print(e, file=result)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Submit this calculation to the HPC cluster as you did on exercise day 3.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq submit fepo4_1li.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.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq ls"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "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": [
        "!cat \"$(ls -t fepo4_1li.py.*err | head -1)\""
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!cat \"$(ls -t fepo4_1li.py.*out | head -1)\""
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "You can move on while you wait for the calculation to finish. Once the calculation is finished load in the structure by running the cell below.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "try:\n",
        "    fepo4_1li=read('fepo4_1li_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",
        "You are now ready to calculate the energy gained by intercalating a single Li ion into the cathode. Start by loading in the relevant reference structures and obtain the potential energies. This should not require any DFT calculations.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Loading in files from exercise day 3.\n",
        "li_metal = read('li_metal.traj')\n",
        "fepo4 = read('fepo4_out.traj')\n",
        "\n",
        "epot_li_metal = li_metal.get_potential_energy() / len(li_metal)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# epot_fepo4 = ...\n",
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Calculate the energy of intercalting a single Li in the FePO$_4$ cell. How does this energy compare with the equilibirum potential? What can it tell you about the charge/discharge potential curves?\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# ...\n",
        "# print(...)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Bonus: LiFePO$_4$ with one vacancy\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "If time permits, you will now do a similar calculation but this time with LiFePO$_4$ contraining one vacancy. Once again you should assume that the cell dimension remain unchanged compaired to LiFePO$_4$.\n",
        "\n",
        "There are numerous ways to obtain this structure. You can get inspiration from the way LiFePO$_4$ was made on Exercise day 3, use the [`del` or `pop()` methods](https://wiki.fysik.dtu.dk/ase/ase/atoms.html?highlight=pop#list-methods), or even use the GUI to delete an atom and save the structure afterwards.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# In this cell you create the vacancy in LiFePO4\n",
        "\n",
        "# lifepo4_vac = ...\n",
        "\n",
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Visualize the structure\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "view(lifepo4_vac)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Now ensure that the total magnetic moment is equal to 17.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "for atom in fepo4_1li:\n",
        "    if atom.symbol == 'Fe':\n",
        "        atom.magmom = 4.25\n",
        "\n",
        "print(sum(fepo4_1li.get_initial_magnetic_moments()))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Write your atoms object to file giving it the name `lifepo4_vac.traj`.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Make a full script in the cell below similar to that you made above. Make sure the cell runs before interupting the notebook kernel.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# %%writefile 'lifepo4_vac.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_vac = read('lifepo4_vac.traj')\n",
        "\n",
        "\n",
        "# ...\n",
        "\n",
        "# write('lifepo4_vac_out.traj', lifepo4_vac)\n",
        "\n",
        "# ens = BEEFEnsemble(calc)\n",
        "# dE = ens.get_ensemble_energies(2000)\n",
        "# with paropen('ensemble_lifepo4_vac.dat','a') as results:\n",
        "#     for e in dE:\n",
        "#         print(e, file=result)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Once you have made sure the cell runs, submit it to the HPC cluster.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "!mq submit lifepo4_vac.py -R 8:1h  # submits the calculation to 8 cores, 1 hour"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Once the calculation has finished, load in the trajectory.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "try:\n",
        "    lifepo4_vac=read('lifepo4_vac_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",
        "Once the calculation has finished you are ready to calculate the energy cost of creating a li vacancy in the fully lithiated LiFePO$_4$. Start by loading in the relevant reference structures and obtain the potential energies. This should not require any calculations.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Loading in files from exercise day 3.\n",
        "li_metal = read('li_metal.traj')   # you should have already read this in above\n",
        "lifepo4 = read('lifepo4_out.traj')\n",
        "\n",
        "epot_li_metal = li_metal.get_potential_energy() / len(li_metal)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# epot_lifepo4 = ...\n",
        "# ...\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# vac_cost = ...\n",
        "# print(vac_cost)\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "How does this energy compare with the equilibirum potential? What can it tell you about the charge/discharge potential curves?\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Bonus\n",
        "Calculate the error estimates of the energy for the added Li atom and vacancy formation using the ensembles.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Cell for bonus question"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Cell for bonus question"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Cell for bonus question\n"
      ],
      "outputs": [],
      "execution_count": null
    }
  ],
  "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
}