Transport barriers and Voltage profile

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.

Create initial and final structures for a NEB calculation, that will determine the transition state
- If time permits you can study the influence of changing the interlayer graphite distance on the energy barrier.
Create structures for a Li vacancy in LiFePO₄ and a single Li in FePO₄
Calculate the Li vacancy/insertion energies and compare them to the equilibrium potential
- What can they tell you about the charge/discharge potential curves?

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.

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.

Transport barrier of Li in graphite

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.

You can use your work from Day 2, but for simplicity you are advised to load in the initial atomic configuration from file.

from ase import Atoms
from ase.mep import NEB
from ase.optimize import BFGS
from gpaw import GPAW, PW
from ase.constraints import FixAtoms

ccdist = 1.39
a = ccdist * 3**0.5
c = 3.7

initial = Atoms(
    'C6Li',
    positions=[(0, 0, 0),
               (0, ccdist, 0),
               (a, 0, 0),
               (-a, 0, 0),
               (-a / 2, -ccdist / 2, 0),
               (a / 2, -ccdist / 2, 0),
               (0, -ccdist, c / 2)],
    cell=([1.5 * a, -1.5 * ccdist, 0],
          [1.5 * a, 1.5 * ccdist, 0],
          [0, 0, c]),
    pbc=(1, 1, 1))

You will now make a final structure, where the Li atom has been moved to a neighbouring equivalent site. The get_positions(), set_positions() and get_cell() methods are highly useful for such a task. HINT: Displace the Li atom \(\frac{1}{n} (\vec{a}+\vec{b})\)

final = initial.copy()
...
...

Visualize that you have made the final strcuture correctly.

from ase.visualize import view
view(final)

Make a band consisting of 7 images including the initial and final.

# These will become the minimum energy path images:
images = [initial]
images += [initial.copy() for i in range(5)]
images += [final]

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.

neb = NEB(images, method='improvedtangent')
neb.interpolate()

Add a view(images) line to your script and run it so that you can visualize the NEB images.

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 the positions of the carbon atoms to save computational time.

Each image in the NEB requires a unique calculator.

This very simple case is highly symmetric. To better illustrate how the NEB method works, the symmetry is broken using the ase.Atoms.rattle() method.

# Constrain C-atoms:
mask = [atom.symbol == 'C' for atom in initial]
constraint = FixAtoms(mask=mask)
for image in images[0:7]:
    calc = GPAW(
        mode=PW(500),
        kpts=(5, 5, 6),
        xc='LDA',
        txt=None,
        symmetry={'point_group': False})
    image.calc = calc
    image.set_constraint(constraint)

images[3].rattle(stdev=0.05, seed=42)

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.

Note, that this can take a while if you opt to do it inside the notebook.

images[0].get_potential_energy()
images[0].get_forces()
images[6].get_potential_energy()
images[6].get_forces()

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.

optimizer = BFGS(neb, trajectory='neb.traj', logfile='neb.log')
optimizer.run(fmax=0.10)

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.

$ mq submit NEB.py -R 8:1h  # submits the calculation to 8 cores, 1 hour
$ ...
$ mq ls
$ ...
$ tail neb.log

You can move on while you wait for the calculation to finish.

Once the maximum force (fmax) in the log is below 0.1, the calculation is finished. Load in the full trajectory.

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.

$ ase gui neb.traj@-7:

In the GUI use Tools ‣ NEB.

Now inspect how the TS image has developed.

$ ase gui neb.traj@3::7

For more complicated MEP’s, use the Climbing image method to determine the transition state. Why is it not required here?

Bonus

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.

from ase.io import read
from gpaw import GPAW, PW

images = read('neb.traj@-7:')
IS_image = images[0]
TS_image = images[3]

Now calculate the energy of the initial state (IS) image and the transition state (TS) image using the get_potential_energy() method.

epot_IS = IS_image.get_potential_energy()
epot_TS = TS_image.get_potential_energy()
barrier = epot_TS - epot_IS
print('Energy barrier:', barrier)

New change the graphite layer distance by changing the the size of the unit cell in the z direction by ±3 %. and use the same calculator object as you did above and calculate the potential energy of the compressed initial and final state.

for strain in [-0.03, 0.03]:
    cell = IS_image.get_cell()
    cell[2] *= 1 + strain
    IS_image97 = IS_image.copy()
    IS_image97.set_cell(cell, scale_atoms=True)
    TS_image97 = TS_image.copy()
    TS_image97.set_cell(cell, scale_atoms=True)
    # Create GPAW calculator and calculate energies:

How does the energy barrier change?

FePO₄ with one Li

You will now calculate the energy gain of adding a single Li atom into the FePO₄ cell you made on Day 3. This corresponds to a charge of 25 %. You can compare this energy to the equilibrium potential.

Start preparing a new Python script (say, fepo4_1li.py) and load in the FePO₄ structure you wrote to file on in a previous exercise and add Li. Assume that the cell dimension remain unchanged.

from ase.io import read, write
from gpaw import GPAW, PW, FermiDirac, Mixer
from ase.dft.bee import BEEFEnsemble
from ase.parallel import paropen

fepo4 = read('fepo4_out.traj')
fepo4_1li = fepo4.copy()
print(fepo4_1li.get_initial_magnetic_moments())

Add a Li atoms at \((x,y,z)=(0,0,0)\):

fepo4_1li.append('Li')
fepo4_1li.positions[-1] = (0, 0, 0)  # not really needed

Now finish the script:

calculate the energy with the BEEF functional
write result to a .traj file: from ase.io import write and write('fepo4_1li_out.traj', fepo4_1li)
calculate a BEEF-ensemble

and submit the job to the queue.

Once the calculation is finished, you are 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 new DFT calculations.

li_metal = read('li_metal.traj')
fepo4 = read('fepo4_out.traj')
fepo4_1li = read('fepo4_1li_out.traj')
epot_li_metal = li_metal.get_potential_energy() / len(li_metal)
epot_fepo4 = ...
epot_fepo4_1li = ...

Calculate the energy of intercalting a single Li in the FePO₄ cell. How does this energy compare with the equilibirum potential? What can it tell you about the charge/discharge potential curves?

Bonus: LiFePO₄ with one vacancy

If time permits, you will now do a similar calculation but this time with LiFePO₄ contraining one vacancy. Once again you should assume that the cell dimension remain unchanged compaired to LiFePO₄.

There are numerous ways to obtain this structure. You can get inspiration from the way LiFePO₄ was made on Exercise day 3. Use del atoms[index], the ase.Atoms.pop() method or even the GUI to delete an atom and save the structure afterwards.

When you have made your script (say lifepo4_vac.py), submit it to the HPC cluster. Once the calculation has finished you are ready to calculate the energy cost of creating a Li vacancy in the fully lithiated LiFePO₄. Start by loading in the relevant reference structures and obtain the potential energies. This should not require any calculations.

# Loading in files from exercise day 3.
li_metal = read('li_metal.traj')   # you should have already read this in above
lifepo4 = read('lifepo4_out.traj')
epot_li_metal = li_metal.get_potential_energy() / len(li_metal)
epot_lifepo4 = lifepo4.get_potential_energy()
vac_cost = ...
print(vac_cost)

How does this energy compare with the equilibirum potential? What can it tell you about the charge/discharge potential curves?

Bonus

Calculate the error estimates of the energy for the added Li atom and vacancy formation using the ensembles.

Transport barriers and Voltage profile

Transport barrier of Li in graphite

Bonus

FePO4 with one Li

Bonus: LiFePO4 with one vacancy

Bonus

FePO₄ with one Li

Bonus: LiFePO₄ with one vacancy