Monte-Carlo simulation combined with density functional theory to investigate the equilibrium thermodynamics of electrode materials: lithium titanates as model compounds
This post is a detailed review and discussion of the paper “Monte-Carlo simulation combined with density functional theory to investigate the equilibrium thermodynamics of electrode materials: lithium titanates as model compounds” by Ozaki, Tada, and Kiyobayashi, published in 2019.
Lithium-titanate-oxide is an anode material whose open-circuit potential stays almost flat over the entire range of cycleable stoichiometry bounds:
The paper answers the question: what is the physical explanation of LTO’s open-circuit voltage plateau?
This page is also a sequel of the page about electrode lithiation stages which focused on graphite anode.
LTO lattice structure
There are two flavours of anode lithium-titanate-oxide: Li_{1-2}Ti_2O_4 and Li_{4/3-7/3}Ti_{5/3}O_4. Compared to the first flavour, in the second flavour one-sixth of Titanium atoms in the lattice are replaced with Lithium atoms. LiTi_2O_4 and Li_{4/3}Ti_{5/3}O_4 correspond to 0% stoichiometry, i. e. there is always some Lithium in LTO anode, unlike in graphite which could be fully delithiated.
Lithium ions can occupy two or three distinct types of sites in the lattice: 8a and 16c, and 16d (the latter only in Li_{4/3-7/3}Ti_{5/3}O_4).
"8a", "16c", "16d", and "32e" are Wyckoff positions, but could be considered more as mnemonics here. The numeric part of the position label indicates the relative proportionality of the positions in the structure: in the case of LTO, it's important to note that there are two times as many 16c sites as 8a sites.
Geometrically, each 8a site shares a face with four 16c sites (and since 8a sites are tetrahedra that have four faces, this means that 8a sites don't share faces with any other types of sites in LTO apart from 16c), but each 16c site shares a face with only two 8a sites. 16c sites don't share faces among themselves. Every distance between the 8a and 16c sites are identical.
LTO almost doesn't expand or contract upon lithiation or delithiation, unlike graphite, so Van der Waals forces don't affect the chemical potential of Lithium intercalation and therefore the shape of LTO's open-circuit voltage function. (However, using Lattice Gas Model to find the OCV of graphite produces acceptable results, too, despite Lattice Gas Model also ignores expansion and contraction of anode particles.)
Energetic model of Lithium intercalated in LTO lattice
The authors assume ε_8a and ε_16c to be the Lithium intercalation energies on 8a and 16c sites, respectively, and J to be the interaction energy between Lithium ions occupying neighbouring 8a and 16c sites, i. e. two sites that share a face. They model the potential energy of a particle with 3L lattice sites for Lithium intercalation (L sites of 8a type and 2L sites of 16c type) H_C(N) where N is the number of Lithium atoms intercalated the sum of energies of all 8a and 16c sites occupied (ε_8a and ε_16c, respectively), plus the number of pairs of neighbouring occupied 8a and 16c sites multiplied by J.
The chemical potential of an Li+ ion intercalation at a certain stoichiometry is calculated as the differential of the potential energy H_C(N) averaged over dozens of stochastic simulation runs
The model decides whether there is a Lithium atom in every individual site in the lattice, unlike Lattice Gas Model. This demands a different method for obtaining Lithium atom configuration that minimises the potential energy: the authors use Markov chain Monte Carlo stochastic simulation (using the Metropolis algorithm) of a particle with 1200 sites (i. e., L = 400) whereas for Lattice Gas Model that is numerical optimisation of the expression of the Gibbs free energy (modelled as potential energy plus configurational entropy) of a system with less than 10 layers, with aggregate stoichiometries of layers as optimisation variables.
Note that the stochasticity of Metropolis simulation (where Boltzmann factor is used as the probability of accepting a randomly generated update) plays the same role as configuration entropy term in the analytical expression of free Gibbs energy in Lattice Gas Model, so we shouldn't additionally include entropy into the H_C(N) energy expression that guides the simulation.
The chemical potential of a Lithium atom intercalation μ (i. e., the open-circuit voltage divided by the charge of Li+ ion: 1e) can be calculated as avg(H_C(N+1)) – avg(H_C(N)), where avg(H_C(N)) is the final H_C(N) of a Monte Carlo simulation run averaged over several dozens of runs (the authors used 36–72), after performing about a million steps in each run, with the temperature initially set to 3000K (which allows almost any Lithium ion "jumps" in the lattice) and then gradually decreased towards normal temperatures (e. g., 300K). This trick is called simulated annealing, it helps the simulation to converge on the globally minimal energy state. A real cell, of course, doesn't need to be heated to 3000K to achieve open-circuit potential: it can do this at a normal temperature when rested for a few hours.
The authors didn't publish the variance within the outcomes of the sets of 36–72 Monte Carlo simulation runs that they averaged to obtain the open-circuit potential function of LTO.
This Monte Carlo approach to calculating open-circuit voltage takes vastly more computation than the approach used with Lattice Gas Model. The latter is not applicable to LTO because there are no separate layers in the LTO structure like in graphite. The positive side is that the fine-grained model helps to better understand the dynamics of the material when visualised:
Density functional theory calculations show that ε_8a, ε_16c, and J energies are roughly the same regardless of the Lithium population of 8a and 16c sites
To validate the hypothesis that the energetic model with just three parameters: ε_8a, ε_16c, and J can be used, the authors used DFT to "more truthfully" compute the potential energies of about a hundred configurations of Li_{1-2}Ti_2O_4 particles with several hundred lattice sites for Lithium atoms and various populations of 8a and 16c sites.
Personally, I don't understand what is said in this sentence, but quote it here for a reference:
The Vienna Ab initio Simulation Package (VASP) was used for the DFT calculations, in which we adopted the Perdew, Burke, and Ernzerhof formula as the exchange-correlation functional and the plane-wave basis with a projector-augmented wave scheme as the basis set.
The authors found values for ε_8a, ε_16c, and J energy parameters (using multiple linear regression) so that the values of the potential energy expression H_C(N) lay within approximately ±5% of the potential energy calculated using DFT for all tested particle configurations. Qualitatively, this means that the potential energy of a Lithium ion in the LTO particle depends primarily on whether the adjacent 8a and 16c sites are occupied, whereas the contribution of second-nearest neighbours, and the mean-field depending on the local stoichiometry is much smaller.
The authors also found that this does not hold for Li_{4/3-7/3}Ti_{5/3}O_4, evidently because the potential energy of a Lithium ion at a 16c site depends significantly on whether it is surrounded by 16d sites occupied by Lithium or Titanium.
Another thing to note about the above conclusion is that although the total potential energy of all Lithium intercalated into a particle predicted by the energetic model described on this page and predicted by DFT differ by at most 5%, the variance among individual Lithium ions could have been higher: for instance, the energetic model could have significantly underestimated the energy of Lithium ions at 16c sites whose vertex-adjacent 16c sites are also occupied, and underestimated the energy of Lithium in other situations. However, since at 0% and 100% almost all Lithium atoms occupy 8a and 16c sites (as explained below), respectively, and the DFT and the energetic model predictions are still quite similar at these extremes, there seems to be no room for significant "hidden" variance.
Completely flat LTO's open-circuit voltage is one big equipotential phase transition: Lithium ions move between 8a to 16c sites
The energy values ε_8a, ε_16c, and J, defined and calculated as described above, have three important properties:
ε_8a is lower than ε_16c.
J is positive, i. e. Lithium ions occupying adjacent 8a and 16c sites repulse from each other.
Δε = ε_16c – ε_8a is about two times smaller than J.
It intuitively follows from this that to minimise the potential energy, in LiTi_2O_4 (0% stoichiometry), almost all Lithium ions stay at 8a sites; in Li_2Ti_2O_4 (100% stoichiometry), almost all Lithium ions stay at 16c sites (any shift from a 16c to an 8a site means three new interactions with neighbouring ions at 16c, so the energy balance of such shift, 3J - Δε, is negative); in Li_xTi_2O_4 where x is between 1 and 2 (stoichiometry between 0% and 100%), the system "prefers" to maintain mostly homogeneous areas of a particle with 8a and 16c sites occupied to minimise the overall number of neighbouring Lithium ions at 8a and 16c sites.
As the stoichiometry rises from 0% to 100%, the areas of the particle bulk (or a single area, e. g. the surface of a particle) where 16c sites are occupied become progressively larger, and the areas (or a single area, e. g. the interior of a particle) where 8a sites are occupied become progressively smaller. During this whole phase transition, the open-circuit voltage of a particle stays at 2ε_16c – ε_8a = ε_8a + 2Δε. This page explains this arithmetic in the case of phase transitions from Stage 1 (LiC_6) to Stage 2 (LiC_12) of graphite which has a 2:1 proportion between stage populations as well LTO with its 8a and 16c sites.
Discussion
The proposed energetic model with three parameters: (with intercalation energies ε_8a and ε_16c into 8a and 16c lattice site types of LTO, respectively, and the repulsion energy J between adjacent 8a and 16c sites both occupied by Lithium atoms) qualitatively explains the LTO's open-circuit potential plateau as a phase transition in which Lithium atoms move from 8a sites (at 0% stoichiometry) to 16c sites (at 100% stoichiometry). This explanation is robust because the energetic model's predictions seem to be accurate within 5%, whereas the difference between Δε = ε_16c - ε_8a and J is much larger: 50%, and the system "should" have two distinct phases as long as J is larger than Δε and the temperature of the particle is not crazily high (e. g., less than 400K). Besides, the phase transition dynamic was experimentally verified, too [2].
The above applies to Li_{1-2}Ti_2O_4 flavour of Lithium-titanate-oxide. For Li_{4/3-7/3}Ti_{5/3}O_4, the model developed in the paper is not sufficient for making a confident assertion that particles undergo simple phase transition between 8a and 16c sites during lithiation or delithiation, yet this is apparently what happens in reality anyway [2]. The actually measured open-circuit potential of Li_{4/3-7/3}Ti_{5/3}O_4 also has a single big plateau as well as Li_{1-2}Ti_2O_4.
Quantitatively, the proposed energetic model predicts the position of LTO's open-circuit voltage plateau within about 3% (e. g., the prediction might be 1.45V whereas the real plateau might be 1.4V, depending on the specific LTO anode: particle size, the preparation process, binders, etc.). This precision might be sufficient for some applications and not sufficient for others, yet the main purpose of the model is a qualitative description of the material dynamics, not accurate fitting of the actual OCV relationship. This is in line with the conjecture made on the page about physics-based vs. equivalent circuit cell models: physics-based models should guide scientific understanding, whereas the models used in a production battery management system (BMS) should not necessarily be physical if learned models are more accurate.