{
"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 2 - Li intercalation energy\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Today we will calculate the energy cost/gain associated with intercalating a lithium atom into graphite using approaches at different levels of theory. After today you should be able to:\n",
"\n",
"- Setup structures and do relaxations using ASE and GPAW.\n",
"- Discuss which level of theory is required to predict the Li intercalation energy and why?\n",
"\n",
"The Li intercalation reaction is:\n",
"$$Li(s) + C_x^{graphite} \\rightarrow LiC_x$$\n",
"and the intercalation energy will then be\n",
"$$E_{Li@graphite} = E_{LiC_x} - (E_{Li(s)} + E_{C_x^{graphite}})$$\n",
"Where $x$ is the number of Carbon atoms in your computational cell.\n",
"\n",
"We will calculate these energies using Density Functional Theory (DFT) with different exchange-correlation functionals. Graphite is characterised by graphene layers with strongly bound carbon atoms in hexagonal rings and weak long-range van der Waals interactions between the layers.\n",
"\n",
"In the sections below we will calculate the lowest energy structures of each compound. This is needed to determine the intercalation energy.\n",
"\n",
"- [Graphite](#graphite)\n",
"- [Li metal](#limetal)\n",
"- [Li intercalation](#liintercalation)\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Graphite\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"To try some of the methods that we are going to use we will start out by utilizing the interatomic potential [EMT](https://wiki.fysik.dtu.dk/ase/ase/calculators/emt.html) (Effective Medium Theory) calculator. This will allow to quickly test some of the optimization procedures in ASE, and ensure that the scripts do what they are supposed to, before we move on to the more precise and time consuming DFT calculations. Initially we will calculate the C-C distance and interlayer spacing in graphite in a two step procedure.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# The graphite structure is set up using a tool in ASE\n",
"from ase.lattice.hexagonal import Graphite\n",
"\n",
"# One has only to provide the lattice constants\n",
"structure = Graphite('C', latticeconstant={'a': 1.5, 'c': 4.0})"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"To verify that the structures is as expected we can check it visually.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from ase.visualize import view\n",
"\n",
"view(structure) # This will pop up the ase gui window"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Next we will use the EMT calculator to get the energy of graphite. Remember absolute energies are not meaningful, we will only use energy differences.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from ase.calculators.emt import EMT\n",
"\n",
"calc = EMT()\n",
"structure.calc = calc\n",
"energy = structure.get_potential_energy()\n",
"print(f'Energy of graphite: {energy:.2f} eV')\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The cell below requires some input. We set up a loop that should calculate the energy of graphite for a series of C-C distances.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"import numpy as np\n",
"\n",
"# Provide some initial guesses on the distances\n",
"ccdist = 1.3\n",
"layerdist = 3.7\n",
"\n",
"# We will explore distances around the initial guess\n",
"dists = np.linspace(ccdist * .8, ccdist * 1.2, 10)\n",
"# The resulting energies are put in the following list\n",
"energies = []\n",
"for cc in dists:\n",
" # Calculate the lattice parameters a and c from the C-C distance\n",
" # and interlayer distance, respectively\n",
" a = cc * np.sqrt(3)\n",
" c = 2 * layerdist\n",
"\n",
" gra = Graphite('C', latticeconstant={'a': a, 'c': c})\n",
"\n",
" # Set the calculator and attach it to the Atoms object with\n",
" # the graphite structure\n",
" calc = EMT()\n",
" gra.calc = calc\n",
"\n",
" energy = gra.get_potential_energy()\n",
" # Append results to the list\n",
" energies.append(energy)\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Determine the equilibrium lattice constant. We use a [np.polynomial.Polynomial](https://www.numpy.org/devdocs/reference/generated/numpy.polynomial.polynomial.Polynomial.html#numpy.polynomial.polynomial.Polynomial)\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# Fit a polynomial:\n",
"poly = np.polynomial.Polynomial.fit(dists, energies, 3)\n",
"# and plot it:\n",
"%matplotlib notebook\n",
"import matplotlib.pyplot as plt\n",
"fig, ax = plt.subplots(1, 1)\n",
"ax.plot(dists, energies, '*r')\n",
"x = np.linspace(1, 2, 100)\n",
"ax.plot(x, poly(x))\n",
"ax.set_xlabel('C-C distance [Ang]')\n",
"ax.set_ylabel('energy [eV]')"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "code",
"metadata": {},
"source": [
"poly1 = poly.deriv()\n",
"poly1.roots() # two extrema"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# Find the minimum:\n",
"emin, ccdist = min((poly(d), d) for d in poly1.roots())\n",
"ccdist"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# alternatively:\n",
"poly2 = poly1.deriv()\n",
"for ccdist in poly1.roots():\n",
" if poly2(ccdist) > 0:\n",
" break\n",
"ccdist"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Make a script that calculates the interlayer distance with EMT in the cell below\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# This script will calculate the energy of graphite for a series of inter-layer distances.\n",
"from ase.calculators.emt import EMT\n",
"from ase.lattice.hexagonal import Graphite\n",
"from ase.eos import EquationOfState\n",
"import numpy as np\n",
"\n",
"# ccdist is already defined in the previous cell\n",
"# Start from a small guess\n",
"layerdist = 2.0\n",
"\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Now both the optimal C-C distance and the interlayer distance is evaluated with EMT. Try to compare with the experimental numbers provided below.\n",
"\n",
"| | Experimental values | EMT | LDA | PBE | PBE+DFTD3 |\n",
"|-------------------------|---------------------|-----|-----|-----|-----------|\n",
"| C-C distance / \u00c5 | 1.420 | | | | |\n",
"| Interlayer distance / \u00c5 | 3.35 | | | | |\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Not surprisingly we need to use more sophisticated theory to model this issue. Below we will use DFT as implemented in GPAW to determine the same parameters.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"First we set up an initial guess of the structure as before.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"import numpy as np\n",
"from ase.lattice.hexagonal import Graphite\n",
"\n",
"ccdist = 1.40\n",
"layerdist = 3.33\n",
"a = ccdist * np.sqrt(3)\n",
"c = 2 * layerdist\n",
"gra = Graphite('C', latticeconstant={'a': a, 'c': c})\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Then we create the GPAW calculator object. The parameters are explained [here](https://wiki.fysik.dtu.dk/gpaw/documentation/basic.html#parameters), see especially the section regarding [mode](https://wiki.fysik.dtu.dk/gpaw/documentation/basic.html#manual-mode). This graphite structure has a small unit cell, thus plane wave mode, `mode=PW()`, will be faster than the grid mode, `mode='fd'`. Plane wave mode also lets us calculate the strain on the unit cell - useful for optimizing the lattice parameters.\n",
"\n",
"We will start by using the LDA exchange-correlation functional. Later you will try other functionals.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from gpaw import GPAW, PW\n",
"\n",
"xc = 'LDA'\n",
"calcname = f'graphite-{xc}'\n",
"calc = GPAW(mode=PW(500), kpts=(10, 10, 6), xc=xc,\n",
" txt=calcname + '.log')\n",
"\n",
"gra.calc = calc # Connect system and calculator\n",
"print(gra.get_potential_energy())\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Check out the contents of the output file, all relevant information about the scf cycle are printed therein:\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"n = calcname + '.log'\n",
"!cat $n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Then we optimize the unit cell of the structure. We will take advantage of the [StrainFilter](https://wiki.fysik.dtu.dk/ase/ase/constraints.html#the-strainfilter-class) class. This allows us to simultaneously optimize both C-C distance and interlayer distance. We employ the [BFGS](http://aria42.com/blog/2014/12/understanding-lbfgs) algorithm to minimize the strain on the unit cell.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from ase.filters import StrainFilter\n",
"from ase.optimize.bfgs import BFGS\n",
"from ase.io import Trajectory\n",
"\n",
"sf = StrainFilter(gra, mask=[1, 1, 1, 0, 0, 0])\n",
"opt = BFGS(sf, trajectory=calcname + '.traj')\n",
"# traj = Trajectory(calcname + '.traj', 'w', gra)\n",
"# opt.attach(traj)\n",
"opt.run(fmax=0.01)\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Read in the result of the relaxation and determine the C-C and interlayer distances.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from ase.io import read\n",
"import numpy as np\n",
"\n",
"atoms = read(calcname + '.traj')\n",
"a = atoms.cell[0, 0]\n",
"h = atoms.cell[2, 2]\n",
"# Determine the rest from here\n",
"print(a / np.sqrt(3))\n",
"print(h / 2)"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Now we need to try a GGA type functional (e.g. PBE) and also try to add van der Waals forces on top of PBE (i.e. PBE+DFTD3). These functionals will require more computational time, thus the following might be beneficial to read.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"If the relaxation takes too long time, we can submit it to be run in parallel on the computer cluster. Remember we can then run different calculations simultaneously. There are two ways to achieve this:\n",
"\n",
"A. Log in to the gbar terminal, save a full script in a file (e.g. `calc.py`) and submit that file to be calculated in parallel (e.g. `mq submit calc.py -R 8:5h` (5 hours on 8 cores)).\n",
"\n",
"or\n",
"\n",
"B. Stay in the notebook and submit the calculations using some commands we will give you, these are explained in the following cells.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"You need to make a full working script in a cell that does not depend on variables defined in any other cells. Do that in the cell below. Once you are done, run it for a few seconds to make sure there is no error, then stop it by `interrupt the kernel`. This will give you an error message, ignore that and move on, we will submit the calculation in the next cell.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# %%writefile graphite_LDA.py\n",
"\n",
"# Full script\n",
"# from ase...\n",
"\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"The following cell submits a calculation using [`myqueue`](https://myqueue.readthedocs.io/en/latest/).\n",
"\n",
"Note the line `#%%writefile graphite_LDA.py` in the previous cell. Remove the `#` symbol and execute the cell again. This will write the cell to a file that you can use submit to the queue:\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"!mq submit graphite_LDA.py -R 8:1h # submits the calculation to 8 cores, 1 hour\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Check the status of the calculation by running the cell below. The column with `S` heading gives you a character code of the status: `Q`: still in queue, `R`: running, `C`: finished.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"!mq ls\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"When the calculation finishes the result can be interpreted by [reading](https://wiki.fysik.dtu.dk/ase/ase/io/io.html#ase.io.read) in the Atoms object. This assumes that the trajectory file is called: `graphite-LDA.traj`, if not change accordingly.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from ase.io import read\n",
"\n",
"# The full relaxation\n",
"relaxation = read('graphite-LDA.traj', index=':')\n",
"# The final image of the relaxation containing the relevant energy\n",
"atoms = read('graphite-LDA.traj')"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# Extract the relevant information from the calculation\n",
"\n",
"# Energy\n",
"print(atoms.get_potential_energy())\n",
"\n",
"# Unit cell\n",
"print(atoms.get_cell())\n",
"\n",
"# See the steps of the optimization\n",
"from ase.visualize import view\n",
"view(relaxation)"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Li metal\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"Now we need to calculate the energy of Li metal. We will use the same strategy as for graphite, i.e. first determine the lattice constant, then use the energy of that structure. This time though you will have to do most of the work.\n",
"\n",
"Some hints:\n",
"1. The crystal structure of Li metal is shown in the image above. That structure is easily created with one of the functions in the [ase.build](https://wiki.fysik.dtu.dk/ase/ase/build/build.html) module\n",
"2. A k-point density of approximately 2.4 points / Angstrom will be sufficient\n",
"3. The DFTD3 correction is done _a posteriori_, this means the calculator should be created a little differently, see [the second example here](https://wiki.fysik.dtu.dk/ase/ase/calculators/dftd3.html#examples)\n",
"4. See also the [equation of state module](https://wiki.fysik.dtu.dk/ase/ase/eos.html)\n",
"\n",
"In the end try to compare the different functionals with experimental values:\n",
"\n",
"| | Experimental values | LDA | PBE | PBE+DFTD3 |\n",
"|-------|---------------------|-----|-----|-----------|\n",
"| a / \u00c5 | 3.51 | | | |\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Get the lattice information and compare with experimental values\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"a =\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Li intercalation in graphite\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"Now we will calculate the intercalation of Li in graphite. For simplicity we will represent the graphite with only one layer. Also try and compare the C-C and interlayer distances to experimental values.\n",
"\n",
"| | Experimental values | LDA | PBE | PBE+DFTD3 |\n",
"|-------------------------|---------------------|-----|-----|-----------|\n",
"| C-C distance / \u00c5 | 1.441 | | | |\n",
"| Interlayer distance / \u00c5 | 3.706 | | | |\n",
"\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"# A little help to get you started with the combined structure\n",
"from ase.lattice.hexagonal import Graphene\n",
"from ase import Atom\n",
"\n",
"# We will have to optimize these distances\n",
"ccdist = 1.40\n",
"layerdist = 3.7\n",
"\n",
"a = ccdist * np.sqrt(3)\n",
"c = layerdist\n",
"\n",
"# We will require a larger cell, to accomodate the Li\n",
"Li_gra = Graphene('C', size=(2, 2, 1), latticeconstant={'a': a, 'c': c})\n",
"# Append an Li atom on top of the graphene layer\n",
"Li_gra.append(Atom('Li', (a / 2, ccdist / 2, layerdist / 2)))\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "code",
"metadata": {},
"source": [],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Now calculate the intercalation energy of Li in graphite with the following formula:\n",
"$$E_{Li@graphite} = E_{LiC_x} - (E_{Li(s)} + x * E_{C^{graphite}})$$\n",
"where $x$ is the number of Carbon atoms in your Li intercalated graphene cell. Remember to adjust the energies so that the correct number of atoms is taken into account.\n",
"\n",
"These are the experimental values to compare with\n",
"In the end try to compare the different functionals with experimental values:\n",
"\n",
"| | Experimental values | LDA | PBE | BEEF-vdW |\n",
"|---------------------------|---------------------|-----|-----|----------|\n",
"| Intercalation energy / eV | -0.124 | | | |\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [],
"outputs": [],
"execution_count": null
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Great job! When you have made it this far it is time to turn your attention to the cathode. If you have lots of time left though see if you can handle the bonus exercise.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Bonus\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"In the calculation of the intercalated Li we used a graphene layer with 8 Carbon atoms per unit cell. We can actually use only 6 Carbon by rotating the x and y cell vectors. This structure will be faster to calculate and still have neglible Li-Li interaction.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"source": [
"from ase import Atoms\n",
"from ase.visualize import view\n",
"\n",
"ccdist = 1.40\n",
"layerdist = 3.7\n",
"\n",
"a = ccdist * np.sqrt(3)\n",
"c = layerdist\n",
"\n",
"Li_gra = Atoms('CCCCCCLi') # Fill out the positions and cell vectors\n",
"\n"
],
"outputs": [],
"execution_count": null
},
{
"cell_type": "code",
"metadata": {},
"source": [],
"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
}