Battery discharge voltage profiles using Grand Canonical Monte Carlo¶

Overview¶

In this tutorial we use the ReaxFF engine and the GCMC method in the AMS driver to calculate the discharge voltage profile of a LiS battery.

Contents of this tutorial:

• Importing a CIF file from an external database and equilibrating the structure.
• Calculating the chemical potential for Li.
• Setting up a Grand Canonical Monte Carlo (GCMC) simulation with AMS.
• Analysis: Writing a Python script that calculates the discharge voltage profile.

This tutorial assumes a basic familiarity with the GUI. If you are not familiar with using the GUI yet, you might take a look at the GUI introduction tutorials before continuing.

The System¶

This tutorial uses a small alpha-sulfur system consisting of 128 atoms. Both the system and the workflow presented here are originally described in the publication “ReaxFF molecular dynamics simulations on lithiated sulfur cathode materials” by M.M. Islam et al..

The discharge process is simulated using ReaxFF in a Grand Canonical Monte Carlo scheme. Volume changes upon lithiation are accounted for by using an N_SPT-μ_Lischeme. The discharge voltage can be calculated from the total energies of the lithiated compounds.

Importing and optimizing the Sulfur(α) crystal structure¶

The crystal structure can be directly imported from a CIF file. There are several resources for crystallographic data available online and you can choose according to your liking. Here we use a structure of S₈ alpha sulfur from the American Mineralogist Crystal Structure Database.

1. Download the CIF file by clicking here

2. In AMSinput: File → Import Coordinates

3. Select the CIF file you just downloaded using the file dialogue window

Before adding any Li-ions to the system, we need to relax the structure using a geometry and cell optimization, i.e. including the optimization of lattice vectors.

1. Select the ReaxFF panel: →

2. Select the LiS.ff force field

3. Select Task → Geometry Optimization

4. Switch to Details → Geometry Optimization

5. Tick the Optimize lattice box

Save and run the calculation to yield a total energy of:

Calculating the chemical potential for Li¶

Following the literature approach, we fix the external chemical potential of Li at the total energy of a single lithium atom in body-centered cubic lithium. The structure of bcc Lithium can be created easily via the crystal builder from within AMSinput:

1. Create a new job by clicking SCM → New Input

2. Edit → Crystal → Cubic → bcc

3. Select Li from the Presets and click OK

4. Edit → Crystal → Generate Supercell

5. Enter 8, 8 , 8 on the diagonal to create an 8x8x8 supercell (512 atoms)

Optimize the resulting structure and lattice using the same settings as with the sulfur, yielding a total energy of:

The chemical potential is calculated as the total energy / number of atoms, which for our supercell with 512 atoms yields:

Setting up the GCMC calculation¶

Our previously optimized S₈ structure will serve as the system:

1. Task → GCMC

2. Switch to the GCMC panel Model → GCMC (or click on )

3. Select MC ensemble → Mu-PT

4. Set Number of GCMC iterations to 1000

This will be enough steps to fully lithiate the cathode in the current calculation. However, this value is not known beforehand and one might want to choose a larger number of iterations for larger systems.

1. Add molecules within 6 Å but not closer than 1.2 Å

These fields correspond to the Rmax and Rmin values in the GCMC input, and refer to the distance to other atoms of the system. As a completely empirical rule of thumb: setting Rmin to roughly ½ of the shortest expected bond and Rmax approx. half of the shortest lattice vector seems a good starting point.

Now we need to specify that we want to add Lithium atoms to the cathode.

1. Press the button labeled + next to the Molecules label

2. In the appearing Molecule drop down menu, select New Molecule

3. A second tab in the viewport on the left will appear and be selected.

4. Add a single Li atom to the empty unit cell using the periodic table tool from the toolbar below. The position of the Li atom within the cell is irrelevant.

5. Optional: Rename the tabs in the viewport (by clicking into the label): Mol-1 → S_cathode and Mol-2 → Li_atom

6. Set the Chemical potential of the Li atom to -37.7 kcal/mol. (The default unit there is Hartree, but you can click on the unit and select kcal/mol.)

Save your changes and run the calculation. The calculation should take around an hour on a typical desktop computer, but you can follow its progress in AMSmovie.

You should see the battery expand as it discharges. There is about a factor of two in volume between the charged and discharged states!

1. Plot the cell volume by clicking Graph → Cell Volume.

2. Enable showing the unit cell as a box by clicking View → Periodic → Show Unit Cell.

You can download the movie here if it does not play in your browser.

GCMC Troubleshoot¶

No MC-moves are accepted:

• Check if you set the correct chemical potential and that you use the right units. Remember to use the one calculated with the method you are also using for the GCMC calculation!
• Change Rmax and Rmin settings.
• Try to optimize the system (here: the sulfur cathode) with tighter optimization settings.

The calculation takes a lot of time:

• Loosen the convergence criteria, e.g. less steps and a lower convergence criterion.
• The obvious: Try a smaller system.

Analyzing the results¶

The discharge voltage profiles can be calculated as a function of Li intercalation content from

where \(G\) denotes the Gibbs free energy and \(x\) the concentration of Li-ions. \(F\) is Faraday’s constant. The enthalpic (pV) and entropic (TS) contributions can be neglected and thus the Gibbs free energy replaced by the ground state energy.

In this case and many other cases of non-standard trajectory analysis, writing a short Python script using the PLAMS library is the most efficient way to obtain results. Remember there is no need to take any further action than writing the script: Both Python and PLAMS are shipped with every copy of the Amsterdam Modeling Suite. Our script might look like this:

# Make sure this script was called correctly
if not 'resultsdir' in globals():
    print("\n USAGE: $AMSBIN/plams LiVoltageProfile.py -v resultsdir=[path to .results folder]\n")
    quit()

# Load the job associated with that .results folder
gcmc = AMSJob.load_external(resultsdir)

# Get energy and number of atoms of the initial state (before adding any Li)
initEnergy = gcmc.results.readrkf('GCMC', 'InitialEnergy')
initNumAtoms = len(gcmc.results.get_input_molecule())

# Get chemical potential of Li
muLi = gcmc.results.readrkf('GCMC','Mol1.chemPot')
# Get Faraday's constant in atomic units
Ha2V = Units.convert(1, 'Hartree','kcal/mol')/23.06
# (Note: 23.06 is Faraday's constant in kcal per volt gram equivalent)

# Loop over entries in the History and make dictionaries for the voltages and volumes,
# using the number of Li atoms as the keys to the dictionaries.
# (Note: We start at history step 2, because the first one is the initial system without any Li.)
voltages = {}
volumes  = {}
for iHistEntry in range(2,gcmc.results.readrkf('History', 'nEntries')+1):

    # Get the geometry for the accepted MC step
    iSys = gcmc.results.get_history_molecule(iHistEntry)

    # Find out how many Li atoms it has
    numLiAtoms = sum(atom.symbol == 'Li' for atom in iSys.atoms)
    if numLiAtoms == 0: continue

    totalEnergy = gcmc.results.readrkf('History', 'Energy({:d})'.format(iHistEntry))
    voltage = - Ha2V * (totalEnergy - initEnergy - numLiAtoms*muLi) / numLiAtoms
    if not numLiAtoms in voltages: voltages[numLiAtoms] = []
    voltages[numLiAtoms].append(voltage)

    volume = iSys.unit_cell_volume()
    if not numLiAtoms in volumes: volumes[numLiAtoms] = []
    volumes[numLiAtoms].append(volume)

# Calulate and print averages for voltages and volumes.
print("# Nr. Li atoms, fraction of Li-atoms, Voltage in V, Volume in Angstrom**3")
for numLiAtoms in voltages:
    voltage = sum(voltages[numLiAtoms]) / len(voltages[numLiAtoms])
    volume  = sum(volumes[numLiAtoms])  / len(volumes[numLiAtoms])
    print("{:3d} {:2.3f} {:2.3f} {:2.2f}" \
          .format(numLiAtoms, numLiAtoms/initNumAtoms, voltage, volume))

The script called LiVoltageProfile.py is available here and can be run as follows from the command line with

$AMSBIN/plams LiVoltageProfile.py -v resultsdir=GCMC.results

assuming that your job was called GCMC and your terminal is in the directory where it was run. The results are printed to the terminal after a few seconds, but you can of course redirect them into a file, which you can plot with standard tools:

[1] M. M. Islam, A. Ostadhossein, O. Borodin, A. Todd Yeates, W. W. Tipton, R. G. Hennig, N. Kumar and A. C. T. van Duin ReaxFF molecular dynamics simulations on lithiated sulfur cathode materials, Phys. Chem. Chem. Phys. 17, 3383-3393 (2015).