Voltage profile from Grand-Canonical Monte Carlo
************************************************

ReaxFF Batteries Discharge voltage profile Chemical potential GCMC


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.

Note:

  The GCMC calculation we are going to run later in this tutorial
  takes approximately an hour to finish on a typical desktop computer.


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..

[image]

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_8 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: [image: ADFPanel][image] **→** [image: ReaxFFPanel][image]
   **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:

   E_\mathrm{S} \approx -8540~\textrm{kcal/mol}


Setting 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:

   E_\mathrm{Li,bulk} \approx -19304~\textrm{kcal/mol}

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

   \mu_\mathrm{Li} = -37.7~\textrm{kcal/mol}

Note:

  The exact value depends to some minor extent on the chosen force
  field and for consistency one should always use the value predicted
  by the force field at use. Throughout this tutorial we use the
  "LiS.ff" forcefield.


Setting up the GCMC calculation
===============================

Our previously optimized S_8 structure will serve as the system:

[image]

   **1.** Open the sulfur optimization job in AMSinput
   **2.** Click the [image: UpdateMolecule][image] button to load the optimized geometry.
   **3.** **Task → GCMC**
   **4.** Switch to the GCMC panel **Model → GCMC** (or click on [image: MoreBtn][image])
   **5.** Select **MC ensemble → Mu-PT**
   **6.** 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 is irrelevant.
   **5.** Go to the **Model → Lattice** panel for the *Li atom system* and set **Periodicity → None**. Double-check that the periodicity for the *Sulfur system* is still **Bulk**.
   **5.** *Optional:* Rename the tabs in the viewport (by clicking into the label): **Mol-1 → S_cathode** and **Mol-1 → Li_atom**
   **7.** 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.)
 [image]

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.

[image]

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 in AMSmovie by clicking **Properties → Lattice Info → 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

   V(x) =  - F \frac{G_{Li_x} - x \cdot G_{Li} - G_S}{x}

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:

   #!/usr/bin/env amspython
   import sys

   from scm.plams import AMSJob, Units


   def main():
       if len(sys.argv) != 2:
           print(f"\nUsage: {sys.argv[0]} <path to .results folder>\n")
           sys.exit(1)

       resultsdir = sys.argv[1]

       # 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
       ha_to_v = 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", f"Energy({iHistEntry:d})")
           voltage = -ha_to_v * (totalEnergy - initEnergy - numLiAtoms * muLi) / numLiAtoms
           if numLiAtoms not in voltages:
               voltages[numLiAtoms] = []
           voltages[numLiAtoms].append(voltage)

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

       # Calculate 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
               )
           )


   if __name__ == "__main__":
       main()

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

"$AMSBIN/amspython LiVoltageProfile.py 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:

[image]

[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).
