# LixCoO2 phase stability studied by machine learning-enabled scale bridging
between electronic structure, statistical mechanics and phase field theories

G.H. Teichert  S. Das  M. Aykol  C. Gopal  V. Gavini  K. Garikipati
<EMAIL_ADDRESS>Department of Mechanical Engineering, University of Michigan
Toyota Research Institute, Los Altos, CA  Department of Materials Science &
Engineering, University of Michigan  Department of Mathematics, University of
Michigan  Michigan Institute for Computational Discovery & Engineering,
University of Michigan

###### Abstract

LixTMO2 (TM=Ni, Co, Mn) are promising cathodes for Li-ion batteries, whose
electrochemical cycling performance is strongly governed by crystal structure
and phase stability as a function of Li content at the atomistic scale. Here,
we use LixCoO2 (LCO) as a model system to benchmark a scale-bridging framework
that combines density functional theory (DFT) calculations at the atomistic
scale with phase field modeling at the continuum scale to understand the
impact of phase stability on microstructure evolution. This scale bridging is
accomplished by incorporating traditional statistical mechanics methods with
integrable deep neural networks, which allows formation energies for specific
atomic configurations to be coarse-grained and incorporated in a neural
network description of the free energy of the material. The resulting
realistic free energy functions enable atomistically informed phase-field
simulations. These computational results allow us to make connections to
experimental work on LCO cathode degradation as a function of temperature,
morphology and particle size.

††journal:

##  1 Introduction

Layered oxides of type LixTMO2 (TM=Ni, Co, Mn) are commonly used cathodes in
Li-ion batteries. Key to their versatility is the ability to tune the
chemistry specific to an application. While LiCoO2 (LCO) has traditionally
been the cathode of choice for consumer electronics, Ni-rich compositions are
attractive candidates for high energy density batteries, which are
increasingly in demand with electrification of transportation. However, the
cycling performance of Ni-rich NMC cathodes is known to suffer from surface
phase transitions [1, 2], mechanical failures [3, 4], and undesired
degradation reactions with electrolytes [5, 6]. A common strategy to improve
the cycling and structural stability of layered oxides is extrinsic doping
with trace amounts of foreign elements [7, 8, 9]. Understanding the defect-
thermodynamics of dopants and how that impacts the crystal structure stability
and long-term electrochemical properties of NMCs from first principles
requires a modeling framework that span lengthscales ranging from the
atomistic to the continuum. However, scale-bridging is inherently a hard
problem to address and prior studies, as we outline below, have looked at the
different length-scales in isolation. Here, we benchmark a recently published
scale-bridging framework [10, 11] against LixCoO2 (LCO), a well-characterized
model cathode, by computing its thermodynamic phase behavior and
microstructure evolution under Li intercalation.

The rich phase-behavior of LCO has been extensively studied by experiments. A
first order phase transition has been identified for compositions $0.75\leq
x\leq 0.94$, caused by a metal-insulator transition [12, 13]. Reimers and Dahn
[12] reported on an ordering that appears at $x=0.5$, confirmed to be a row
ordering [14, 15]. Shao-Horn et al. [14] found evidence of ordering at
$x=0.33$ at low temperatures ($\sim$100 K). Additionally, charge ordering Co3+
and Co4+ atoms at compositions of $x=0.5$ and $x=0.67$ was observed at a
temperature of 175 K, while Co3+ and Co4+ occupation was random at room
temperature [16, 15]. Multiple studies have also identified a transition from
the O3 host structure to the H1-3 structure at compositions below $x=0.33$
[17, 18].

Alongside the experimental work have been the numerous computational studies
that have also explored LCO. _Ab initio_ thermodynamic approaches combining
density functional theory (DFT) calculations and statistical mechanics [19,
20] have been used to systematically predict phase diagrams, free energy data,
and voltage curves. The first-principles methods also have been combined with
experimental results and CALPHAD models to define functions for the Gibbs
energy [21, 18]. However, quantities computed from the atomistic and lattice
models have not informed the more coarse-grained methods, such as phase field
models to study the thermodynamics and kinetics of phase transitions in these
systems. Recently, a continuum level phase field approach has combined free
energy and misfit strain information from experiments with the predicted
diffusivity from first-principles calculations to model the metal-insulator
phase transition that has been seen experimentally [22]. In effect, while the
link between _ab initio_ and statistical mechanics computations is well
established, there remains an absence, so far, of comprehensive scale bridging
frameworks.

The current work seeks to demonstrate a systematic framework combining some of
the computational approaches described above with advances in machine learning
to bridge electronic structure calculations, statistical mechanics approaches
and continuum scale simulations, Figure 1. Our use of O3 layered LCO to
benchmark the approach is motivated by the observation that, to the best of
our knowledge, this path to scale bridging has not been presented for this
cathode material.

We perform DFT$+U$ calculations using the well-known Hubbard correction [23]
to capture electronic correlation in transition metal oxides with localized
$d$ orbitals. Combined with van der Waals functionals, this ensures closer
agreement with experimental voltages in layered LCO [24]. A cluster expansion
Hamiltonian is parameterized by the DFT and used in Monte Carlo sampling to
compute the point defect thermodynamics as a function of temperature and
chemical potential, and to construct temperature - composition phase diagrams
[25, 26].

For bridging statistical mechanics and continuum representations of the free
energy density function (and its derivatives), we turn to integrable deep
neural networks (IDNNs) introduced previously [10, 11]. IDNNs can be trained
to derivative data, e.g. chemical potential data generated through Monte Carlo
sampling, and analytically integrated to recover the antiderivative function,
e.g. the free energy function needed in phase field models. The chemical
potentials and free energies thus obtained are used to confirm a match with
extensively reported phase diagrams for LCO, as well as voltage curves over
the composition range of $[0.33,1]$. Finally, using the obtained free
energies, we carry out phase field studies that take a step beyond the
construction of equilibrium phase diagrams by accessing dynamics of phase
formations and particle morphologies. Specifically, we study spinodal
decomposition as well as charge-discharge cycles on single LCO particles at
temperatures of practical interest, and for different particle sizes. These
computational results allow us to make connections to experimental work on LCO
cathode degradation as a function of temperature, morphology and particle
size. By this, they allow us to address experimental and technologically
relevant scenarios. They also serve as a culminating demonstration of
predictive computations enabled by our scale bridging framework.

Figure 1: Flowchart outlining the data and methods used to bridge from the
atomic scale up to the continuum scale.

The paper is laid out as follows: Section 2 sets forth the methods used in the
paper consisting of the DFT calculations, statistical mechanics sampling, IDNN
setup and training, and the phase field equations. Section 3 describes and
discusses the results of each step of the framework while making comparisons
with experimental observations. Discussion and conclusions are presented in
Section 4.

##  2 Methods

In this section we detail the methods that were previewed in the introduction,
starting with the DFT calculations and moving on to the statistical mechanics
simulations that they inform. This is followed by the machine learning
techniques, specifically IDNNs [10, 11] that enable bridging of the atomic to
continuum scales via free energy density representations. This provides a
principled basis for continuum simulations of the dynamics of LCO particle
morphology via phase field models whose computational treatment is the
culmination of the scale bridging framework (Figure 1).

###  2.1 Density functional theory methods

The formation energy $E$ of a configuration with Li composition $x$ can be
calculated using total energy values computed using density functional theory
(DFT) according to the following equation [19]:

$\displaystyle E$
$\displaystyle=E^{\mathrm{tot}}-xE^{\mathrm{tot}}_{\mathrm{LiCoO}_{2}}-(1-x)E^{\mathrm{tot}}_{\mathrm{CoO}_{2}}$
(1)

where $E^{\mathrm{tot}}$, $E^{\mathrm{tot}}_{\mathrm{LiCoO}_{2}}$, and
$E^{\mathrm{tot}}_{\mathrm{CoO}_{2}}$ are the total energies for the given
configuration, LiCoO2, and CoO2, respectively. In performing DFT calculations
for LixCoO2, Aykol and co-authors found that it is important to include not
only the Hubbard correction (DFT+$U$), which reduces self-interaction errors
[27, 28, 29]), but to also include the van der Waals (vdW) interactions [24].
The effect of the vdW interactions allows the voltage predicted by DFT to
match the experimental voltage when using an appropriately tuned value of $U$.
While both vdW-corrections and vdW-density functionals (vdW-DF) improve the
DFT results, predictions are most consistent and accurate with vdW-DF.

In this workflow, we use a simplified rotational-invariant formulation of
DFT+$U$ [23] and a vdW-DF exchange correlation functional, namely the optB88
exchange correlation functional [30, 31, 32, 33, 24, 34, 35], to calculate the
formation energy on a chosen subset of LCO configurations. We perform the
ground-state DFT calculations with geometry optimization in Quantum Espresso
[36, 37], using projector augmented-wave (PAW) pseudopotentials calculated
with the Perdew-Burke-Ernzerhof (PBE) functional from PSlibrary 1.0.0 [38,
39]. The values for the wave function and charge density cutoffs are chosen in
two steps. Starting with the cutoff values suggested in the pseudopotential
file for Co, we first increase only the charge density cutoff until the total
energy converges to $<$1 meV/atom. Next, we increase both the wave function
and the charge density cutoffs, maintaining the ratio of the two values, until
the total energy again converges to $<$1 meV/atom, giving a wave function
cutoff of 55 Ha and a charge density cutoff of 301.5 Ha. A k-point grid of
$6\times 6\times 3$ is also used to ensure total energy convergence within
$<$1 meV/atom. Structural optimization is performed until cell stress and
ionic forces are under 0.5 kbar and 0.00005 Ha/Bohr, respectively. The crystal
structure for LiCoO2 is shown in Figure 2.

(a) Ball-and-stick

(b) Polyhedral

Figure 2: Two representations of the crystal structure of layered O3 LiCoO2,
where the Li atoms are green, Co are blue, and O are red (visuals generated
using VESTA [40]).

In this work, we calibrate the Hubbard $U$ parameter to match experimental
average lithiation voltages over various ranges of Li composition. To do this,
we compute the voltage at increasing $U$ for given composition values of
$x_{1}$ and $x_{2}$ [41], using crystal structures previously reported from
experiments. The average voltage from $x_{1}$ to $x_{2}$ is calculated using
the following equation [42, 24]:

$\displaystyle
V=-\frac{E[\mathrm{Li}_{x_{2}}\mathrm{CoO}_{2}]-E[\mathrm{Li}_{x_{1}}\mathrm{CoO}_{2}]-(x_{2}-x_{1})E[\mathrm{Li}]}{(x_{2}-x_{1})e}$
(2)

where $E$ is the calculated total energy from DFT and $e$ is the elementary
charge. We compute the average voltage over the following ranges of Li
composition: $x\in[0,1]$, $x\in[0,1/2]$, and $x\in[1/2,1]$. By comparing with
experimental voltages (Figure 3), we select an appropriate value of $U$.

Figure 3: Calculated voltage across various compositions as a function of U,
compared with experimental voltage [24].

###  2.2 Statistical mechanics

In this work, we follow the statistical mechanics approach outlined by Van der
Ven, et al. [19] in their first-principles study of LixCoO2 (LCO), and of
Natarajan, et al. [43], among others.

####  2.2.1 Cluster expansion for formation energy

Given the expense of DFT calculations, we adopt cluster expansions to access
the large numbers of configurations needed in the statistical mechanics
studies. The formation energy is written as $E(\boldsymbol{\sigma})$, where
$\boldsymbol{\sigma}$ is the configuration vector with the occupancy variable
$\sigma_{i}=1$ if Li occupies the site $i$ and $\sigma_{i}=0$ if the site
contains a vacancy. As is common, we define $E(\boldsymbol{\sigma})$ by
computing the formation energy with DFT for a subset of configurations and use
these values to parameterize a cluster expansion [44, 45] as a rapidly
queryable surrogate for the formation energy. We use the CASM software[46],
which facilitates the construction and parameterization of cluster expansion
Hamiltonians and their use in Monte Carlo simulations, to select
configurations for the DFT computations and perform the statistical mechanics
calculations in this work [47, 48, 49]. Candidate configurations are chosen
for O3-LixCoO2 with Li compositions from $x=0$ to $x=1$.

A cluster is a collection of potential Li sites. Given cluster $\alpha$ of
sites $\alpha=\\{i,j,\dots,k\\}$, a polynomial $\phi_{\alpha}$ can be defined
as the product of occupancy variables of those sites, i.e.
$\phi_{\alpha}=\sigma_{i}\sigma_{j}\cdots\sigma_{k}$. A cluster expansion is a
linear combination of the polynomials $\phi_{\alpha}$, leading to the
following form for the formation energy:

$\displaystyle E(\boldsymbol{\sigma})$
$\displaystyle=V_{0}+\sum_{\alpha}V_{\alpha}\phi_{\alpha}(\boldsymbol{\sigma}),$
(3)

where the coefficients $V_{0}$ and $V_{\alpha}$ are optimized for a “best
fit”, and are called effective cluster interactions (ECI).

We fit a cluster expansion to the formation energy calculated using DFT for
333 configurations for the O3 host structure using the CASM code. A sparse
regression technique that combines a genetic algorithm with weighted linear
regression is used to perform the fit. In this approach, a number of candidate
basis functions are created using pairs, triplets, etc. of lattice sites that
are within a given distance. The genetic algorithm is used to select a subset
of these basis functions to include in the cluster expansion. The coefficients
of each subset of basis functions are calculated using linear regression. We
use cluster types up to quadruplets, with maximum lattice site distances shown
in Table 1, for a total of 221 candidate basis functions. During the linear
regression, we apply weights to bias the fit towards greater accuracy for
configurations on or near the convex hull. Additionally, the cluster expansion
is selected such that the convex hull constructed from the cluster expansion
predictions for all 333 configurations consists of the same configurations as
in the DFT convex hull.

Table 1: Basis function specifications used for cluster expansion, resulting in a total of 221 candidate basis functions. cluster type | max. distance (Å)
---|---
pairs | 24
triplets | 8
quadruplets | 6

####  2.2.2 Monte Carlo sampling

Given $E(\boldsymbol{\sigma})$, we sample within the semi-grand canonical
ensemble, in which the chemical potential is specified and the corresponding
composition and/or order parameters are determined through ensemble averaging.
The partition function for the semi-grand canonical ensemble is the following:

$\displaystyle\Theta$
$\displaystyle=\sum_{\boldsymbol{\sigma}}\exp{\left(-\frac{E(\boldsymbol{\sigma})-M\left[x(\boldsymbol{\sigma})\mu_{0}+\boldsymbol{\eta(\boldsymbol{\sigma})}\cdot\boldsymbol{\mu}\right]}{k_{B}T}\right)}$
(4)

where $\boldsymbol{\eta}$ represents the vector of order parameters explained
below, $\mu_{0}$ is the chemical potential associated with composition,
$\boldsymbol{\mu}$ is the vector of chemical potentials associated with the
order parameters, $M$ is the number of reference supercells that tile the
configuration, $k_{B}$ is the Boltzmann constant, and $T$ is the temperature.
Each chemical potential is the derivative of the free energy with respect to
its corresponding composition or order parameter.

The chemical potential associated with Li composition is related to the
voltage as a function of composition, $V(x)$ as

$V(x)=-\frac{\mu_{\mathrm{Li}}^{\mathrm{cathode}}(x)-\mu_{\mathrm{Li}}^{\mathrm{anode}}}{e},$
(5)

where $\mu_{\mathrm{Li}}^{\mathrm{cathode}}$ is the chemical potential with
respect to Li for the cathode (LCO) and $\mu_{\mathrm{Li}}^{\mathrm{anode}}$
is the chemical potential of the anode. Taking the anode to be Li metal,
$\mu_{\mathrm{Li}}^{\mathrm{anode}}$ is a constant equal to the Gibbs free
energy of Li [42]. It is necessary for the two chemical potentials to have a
consistent reference, however. We approximate the Gibbs free energy of Li with
the total energy calculated using DFT, but the chemical potential for LCO
reported from Monte Carlo sampling is based on formation energies that were
re-referenced to the end members, according to the linear relation in Eq. (1).
Since the chemical potential is the derivative of the free energy with respect
to composition, this re-referencing shifted the chemical potential of LCO by a
constant $k$ equal to the derivative of Eq. (1) with respect to $x$, i.e.
$k:=-E^{\mathrm{tot}}_{\mathrm{LiCoO}_{2}}+E^{\mathrm{tot}}_{\mathrm{CoO}_{2}}$.
Therefore, we reverse the shift in the chemical potential by subtracting $k$
from $\mu_{\mathrm{Li}}^{\mathrm{cathode}}(x)$ before computing the voltage
according to Eq. (5). This imposes a consistent reference for both
$\mu_{\mathrm{Li}}^{\mathrm{cathode}}$ and
$\mu_{\mathrm{Li}}^{\mathrm{anode}}$.

A limitation with Monte Carlo sampling within the semi-grand canonical
ensemble is that it does not produce data within the unstable regions
associated with phase separation. While this is sufficient for producing phase
diagrams and predicting voltage, phase field simulations require free energy
information within these two-phase regions in order to consistently resolve
phase interfaces. For this reason, we also employ umbrella sampling methods,
in which bias potentials are applied to sample within the unstable regions
that would otherwise be missed [43, 50, 51, 52]. The partition function used
in this case, then, is the following:

$\displaystyle\Theta$
$\displaystyle=\sum_{\boldsymbol{\sigma}}\exp{\left(-\frac{E(\boldsymbol{\sigma})+M\left[\phi_{0}(x(\boldsymbol{\sigma})-\kappa_{0})^{2}+\sum_{i=1}^{n}\phi_{i}(\eta_{i}(\boldsymbol{\sigma})-\kappa_{i})^{2}\right]}{k_{B}T}\right)}$
(6)

where $\phi_{i}$ and $\kappa_{i}$ determine the curvature and center of the
bias potential, respectively, and $n$ is the number of order parameters.
Usually, $\phi_{i}$ can be fixed at an appropriate value while $\kappa_{i}$ is
varied to sample across the desired composition and order parameter space
within the Monte Carlo sampling routine. The ensemble average of the
composition $\langle x\rangle$ and each order parameter
$\langle\eta_{i}\rangle$ is related to its corresponding chemical potential
through the bias parameters:

$\displaystyle\frac{1}{M}\mu_{0}\Big{|}_{\langle
x\rangle,\langle\boldsymbol{\eta}\rangle}$ $\displaystyle=-2\phi_{0}(\langle
x\rangle-\kappa_{0})$ (7) $\displaystyle\frac{1}{M}\mu_{i}\Big{|}_{\langle
x\rangle,\langle\boldsymbol{\eta}\rangle}$
$\displaystyle=-2\phi_{i}(\langle\eta_{i}\rangle-\kappa_{i}),\qquad
i=1,\ldots,n$ (8)

####  2.2.3 Symmetry-adapted order parameters

In addition to composition, symmetry-adapted order parameters are useful for
observing and tracking order-disorder transitions, as well as for identifying
the various translational and rotational variants of a given ordering. The
method for identifying symmetry-adapted order parameters is laid out
elsewhere; e.g., see Natarajan, et al. [43]. In this process, a supercell
capable of representing all variants of the ordering is found, called a
mutually commensurate supercell. The sublattice compositions of this supercell
form a basis that can describe each variant as a vector $\boldsymbol{x}$,
where each component is 1 if the corresponding sublattice sites are fully
occupied by Li and 0 otherwise. A matrix $\boldsymbol{Q}$ that transforms the
sublattice compositions $\boldsymbol{x}$ to a vector of order parameters
$\boldsymbol{\eta}$ is constructed following an algorithm described by Thomas
and Van der Ven [53].

The algorithm involves enumerating the group $P$ of matrices that transform
between all variants of the given ordering. A matrix is formed, defined by its
invariance to all transformation matrices in $P$. The eigenvalue decomposition
of this $P$-invariant matrix is used to construct $\boldsymbol{Q}$. While
Thomas and Van der Ven deal with a general system, in simpler cases each
nonzero eigenvalue is associated with either the average composition or the
set of variants for a particular ordering. The eigenvectors form the rows of
the matrix $\boldsymbol{Q}$, defining the transformation from sublattice
compositions to order parameters. In some cases, $\boldsymbol{Q}$ maps to more
order parameters than are necessary to describe a single ordering. When that
occurs, a subset of rows can be extracted from $\boldsymbol{Q}$ to create a
reduced set of relevant order parameters. This algorithmic process is easily
understood via an example, which is described in section 3.2.3.

###  2.3 Integrable deep neural networks for free energy representations

The IDNN representation is obtained for the free energy density function by
training on derivative data in the form of pairs of chemical potential–the
label–and corresponding composition or order parameter–the features [10, 11].
The integrability of the IDNN is built into its structure by constructing it
as the derivative form of a standard fully-connected deep neural network
(DNN). While it is straightforward to differentiate the equations describing
the standard DNN to derive the equations for the IDNN, modern deep learning
libraries make this step unnecessary. The user can simply construct a standard
DNN and apply a gradient operation to create the IDNN, which is then used for
training. Mathematically, a DNN can be denoted by a function
$\boldsymbol{Y}(\boldsymbol{X},\boldsymbol{W},\boldsymbol{b})$ representing
the chemical potentials $\langle\mu_{i}\rangle$, with arguments or inputs
$\boldsymbol{X}$ representing the composition $\langle x\rangle$, order
parameters $\langle\eta_{i}\rangle$, weights $\boldsymbol{W}$, and biases
$\boldsymbol{b}$. Training the DNN involves an optimization problem for the
weights and biases, given the dataset
$\\{(\boldsymbol{\widehat{X}}_{\theta},\widehat{Y}_{\theta_{k}})\\}$:

$\displaystyle\boldsymbol{\widehat{W}},\boldsymbol{\widehat{b}}=\underset{\boldsymbol{W},\boldsymbol{b}}{\mathrm{arg\,min}}\,\mathrm{MSE}\left(\boldsymbol{Y}(\boldsymbol{X},\boldsymbol{W},\boldsymbol{b})\Big{|}_{\boldsymbol{\widehat{X}}_{\theta}},\widehat{Y}_{\theta}\right)$
(9)

The case of interest here is when the dataset
$\\{(\boldsymbol{\widehat{X}}_{\theta},\widehat{Y}_{\theta_{k}})\\}$ is not
available. Instead, we have the derivative dataset
$\\{(\boldsymbol{\widehat{X}}_{\theta},\widehat{y}_{\theta_{k}})\\}$, where
$\widehat{y}_{\theta_{k}}$ corresponds to the $k$-th partial derivative of
$\widehat{Y}_{\theta}$. To use these data, the IDNN is defined as the gradient
of $Y$ with respect to its inputs $\boldsymbol{X}$, i.e.
$\partial\boldsymbol{Y}(\boldsymbol{X},\boldsymbol{W},\boldsymbol{b})/\partial
X_{k}$. The training is defined as follows:

$\displaystyle\boldsymbol{\widehat{W}},\boldsymbol{\widehat{b}}=\underset{\boldsymbol{W},\boldsymbol{b}}{\mathrm{arg\,min}}\,\sum_{k=1}^{n}\mathrm{MSE}\left(\frac{\partial\boldsymbol{Y}(\boldsymbol{X},\boldsymbol{W},\boldsymbol{b})}{\partial
X_{k}}\Big{|}_{\boldsymbol{\widehat{X}}_{\theta}},\widehat{y}_{\theta_{k}}\right)$
(10)

The resulting optimized weights $\boldsymbol{\widehat{W}}$ and biases
$\boldsymbol{\widehat{b}}$ can be used with the function
$\partial\boldsymbol{Y}(\boldsymbol{X},\boldsymbol{\widehat{W}},\boldsymbol{\widehat{b}})/\partial
X_{k}$ to return a prediction of the chemical potential. Its antiderivative is
exactly represented by using the same weights and biases in the function
$\boldsymbol{Y}(\boldsymbol{X},\boldsymbol{\widehat{W}},\boldsymbol{\widehat{b}})$.
For the current work,
$\partial\boldsymbol{Y}(\boldsymbol{X},\boldsymbol{\widehat{W}},\boldsymbol{\widehat{b}})/\partial
X_{k}$ gives the IDNN representation of the chemical potentials, and
$\boldsymbol{Y}(\boldsymbol{X},\boldsymbol{\widehat{W}},\boldsymbol{\widehat{b}})$
is the DNN representation of the free energy.

We use the Keras and Tensorflow libraries to train IDNNs using the chemical
potential and composition data for temperatures of 260 K, 300 K, and 340 K.
The IDNN used consists of two hidden layers. For our initial studies we have
used one input (Li composition) and one output (chemical potential with
respect to composition). The IDNNs for 260 K, 300 K, and 340 K are constructed
with 100, 60, and 20 neurons per layer, respectively, arrived at by
hyperparameter optimization, due to differing complexity in the data at
different temperatures. The hyberbolic tangent is used as the activation
function, since we found it to better capture the small fluctuations in
chemical potential around $x=0.5$ than the softplus function used by us in
previous work [10, 11]. The IDNNs are trained for 50,000 epochs using the
RMSprop optimizer and a batch size of 20. An initial learning rate of 0.002 is
used, with periodic reductions in the learning rate applied during training.
The training sets consist of 90, 85, and 69 data points for 260 K, 300 K, and
340 K, respectively.

###  2.4 Phase field theory and associated computational framework

The evolution of microstructure and phase changes can be modeled using the
phase field equations. The Cahn-Hilliard equation [54] models the dynamics of
conserved quantities, such as composition, while nonconserved order parameter
fields are modeled using the Allen-Cahn equation [55]. When neglecting elastic
effects, the total free energy of the system with $n$ order parameters can be
described as follows:

$\displaystyle\Pi[x,\boldsymbol{\eta}]=\int\limits_{\Omega}\left(f(x,\boldsymbol{\eta})+\frac{1}{2}\chi_{0}|\nabla
x|^{2}+\sum_{i=1}^{n}\frac{1}{2}\chi_{i}|\nabla\eta_{i}|^{2}\right)\,\mathrm{d}V$
(11)

where $\chi_{i}$ are the gradient parameters, and $f(x,\boldsymbol{\eta})$ is
the free energy density, to be represented by an analytically integrated DNN
in this work.

The chemical potentials $\widetilde{\mu}_{i}$ used in the phase field
equations are given by the variational derivatives of the total free energy,
such that $\widetilde{\mu}_{0}:=\delta\Pi/\delta x$ and
$\widetilde{\mu}_{i}:=\delta\Pi/\delta\eta_{i}$, $i=1,\ldots,n$:

$\displaystyle\widetilde{\mu}_{0}$ $\displaystyle=\frac{\partial f}{\partial
x}-\chi_{0}\nabla^{2}x$ (12) $\displaystyle\widetilde{\mu}_{i}$
$\displaystyle=\frac{\partial
f}{\partial\eta_{i}}-\chi_{i}\nabla^{2}\eta_{i},\qquad i=1,\ldots,n$ (13)

The Cahn-Hilliard and Allen-Cahn equations, respectively, are the following:

$\displaystyle\frac{\partial x}{\partial t}$
$\displaystyle=\nabla\cdot(\widetilde{M}\nabla\widetilde{\mu}_{0})$ (14)
$\displaystyle\frac{\partial\eta_{i}}{\partial t}$
$\displaystyle=-L\widetilde{\mu}_{i},\qquad i=1,\ldots,n$ (15)

where $\widetilde{M}$ is the mobility and $L$ is a kinetic coefficient. We
substitute in the equations for the chemical potentials and write the
governing equations in weak form to be solved using a mixed finite element
method:

$\displaystyle 0$ $\displaystyle=\int_{\Omega}\left(w_{x}\frac{\partial
x}{\partial t}+\widetilde{M}\nabla
w_{x}\cdot\nabla\tilde{\mu}_{0}\right)\mathrm{d}V-\int_{\partial\Omega}w_{x}j_{n}\mathrm{d}S$
(16) $\displaystyle 0$
$\displaystyle=\int_{\Omega}\left[w_{\tilde{\mu}_{0}}\left(\tilde{\mu}_{0}-\frac{\partial
f}{\partial x}\right)-\chi_{0}\nabla w_{\tilde{\mu}_{0}}\cdot\nabla
x\right]\mathrm{d}V$ (17) $\displaystyle 0$
$\displaystyle=\int_{\Omega}\left[w_{i}\frac{\partial\eta_{i}}{\partial
t}+L\left(w_{i}\frac{\partial f}{\partial\eta_{i}}+\chi_{i}\nabla
w_{i}\cdot\nabla\eta_{i}\right)\right]\mathrm{d}V,\qquad i=1,\ldots,n$ (18)

where $w_{x}$, $w_{\tilde{\mu}_{0}}$, and $w_{i}$ are weighting functions. For
the equations written in this mixed formulation, the following Neumann
boundary conditions have been applied to $x$, $\tilde{\mu}_{0}$, and
$\eta_{i}$, $i=1,\ldots,n$, on $\partial\Omega$, where $\boldsymbol{n}$ is the
outward unit normal and $j_{n}$ is an influx:

$\displaystyle\nabla x\cdot\boldsymbol{n}$ $\displaystyle=0$ (19)
$\displaystyle\widetilde{M}\nabla\tilde{\mu}_{0}\cdot\boldsymbol{n}$
$\displaystyle=j_{n}$ (20) $\displaystyle\nabla\eta_{i}\cdot\boldsymbol{n}$
$\displaystyle=0,\qquad i=1,\ldots,n$ (21)

Figure 4: Surrogate function for the diffusivity of LCO in the Li composition
range $x\in[0.5,1]$, plotted against first-principles diffusivity data of LCO
at 300 K from Ref. [56].

To compare the behavior of LCO at different temperatures, we perform phase
field simulations at 260 K, 300 K, and 340 K, informed by data from the Monte
Carlo computations at these temperatures. The analytically integrated free
energy DNNs for each of the three temperatures are used. We define the
following composition-dependent surrogate function for the diffusivity $D$ at
300 K, using the predicted values from Ref. [56]:

$D=0.01\exp(-274(1.05-x)(0.47-x)(1.-x))$ (22)

The diffusivity surrogate function and the predicted data are plotted in
Figure 4, where the effective vibrational frequency $\nu^{*}$ is reported to
be on the order of $10^{13}$ s-1 [56]. The mobility $\widetilde{M}$ is related
to $D$ by the equation $\widetilde{M}=x/(k_{B}T)$ [57]. The value of $D$ is
multiplied by 4 to approximate the diffusivity at 340 K and is divided by 4
for 260 K. The factor of 4 for every 40 K is based on the predicted
diffusivity at 400 K being approximately 30 times the diffusivity at 300 K
[56]. Since the free energy for these initial simulations is a function of Li
composition only, i.e., a conserved order parameter, the Cahn-Hilliard
equation from Eq. (14) is sufficient to model the system. The non-conservative
Allen-Cahn equation in Eq. (15) is not relevant without inclusion of the non-
conserved order parameters, which are obtained from rows 2-7 of the
$\widehat{\boldsymbol{Q}}$ matrix presented in Section 3.2.3. The simulations
are performed using the finite element method with the mechanoChemFEM
code111Code available at github.com/mechanoChem/mechanoChemFEM, which is based
on the deal.II [58] library, and run on the ConFlux HPC cluster at the
University of Michigan. Adaptive meshing with hanging nodes and adaptive time
stepping are used.

##  3 Results

###  3.1 Formation energies and configurations from density functional theory

The predicted average voltage for Li composition ranges [0,1], [0,0.5], and
[0.5,1] and increasing values of $U$ are plotted in Figure 3. By comparing
with experimental voltages from Ref. [24], we selected a value of $U=2.5$ eV
to use in the calculation of formation energies for LCO. This provides a good
match with the experimental voltage for $x\in[0,1]$ and $x\in[0.5,1]$, with a
small error of 0.03 eV for $x\in[0,0.5]$.

CASM was used to identify a set of configurations with the O3 crystal
structure for use in parametrizing the cluster expansion. Using the tuned
value for $U$, fully-relaxed DFT calculations were completed for 333 of these
configurations, with CASM automatically adjusting the size of the k-points
grid according to atomic configuration. Due to difficulty in convergence for
some configurations, it was necessary in those cases to incrementally increase
$U$ up to 2.5 to achieve a fully converged result. The calculated formation
energies and associated convex hull with the predicted ground states are
plotted in Figure 5. The unit cell volume for LiCoO2 was calculated to be
32.502 Å3.

Figure 5: Formation energies and convex hull from the DFT calculations.

(a) Configuration A

(b) Configuration B

(c) Configuration C

Figure 6: Three configurations with Li concentration $c=0.5$ with the
following calculated formation energies: (a) -0.150, (b) -0.164, and (c)
-0.170 eV/primitive cell. The filled circles represent Li atoms within a Li
layer, and the empty circles represent Li atoms in the adjacent Li layer.

Three configurations at $x=0.5$ are shown in Figure 6 that have been discussed
in the literature, all with low formation energies. Note that configurations A
and B both exhibit a row ordering and are distinguished only by the relative
stacking of the Li layers, whereas configuration C shows a zig-zag ordering.
Experimental work has found evidence for the row ordering [12, 14, 15], where
Takahashi, et al. [15] present configuration A as the crystal structure and
Shao-Horn, et al. [14] suggest configuration B might be possible. Charge
ordering of Co3+ and Co4+ atoms has been experimentally observed at low
temperatures and Li compositions of $x=1/2$ and $x=2/3$ [16]. The development
of charge ordering can be monitored through the magnetic moments in the DFT
calculations [24]. Interestingly, evidence of charge splitting appears in the
DFT results in this work for configuration A, but not in configurations B or
C. However, a more accurate investigation of charge ordering would require use
of advanced DFT functionals such as hybrid exchange functionals [59, 60, 61]
or extended Hubbard models (DFT+$U$+$V$) [62, 63].

In their first-principles study, Van der Ven, et al. [19] found configuration
B (as labeled in Figure 6) to have the lowest formation energy, followed
closely by configuration C. They make the observation, however, that the
energy difference is small enough that it is difficult to determine with
certainty the true ground state solely from the first-principles calculations.
In our work, we take the zig-zag ordering of configuration C to be the ground
state based on the calculated formation energies, but we repeat the caveat
that its calculated energy is close to that of the configurations with row
ordering. The differences in energy between our results and those of earlier
work are likely due to our using DFT+$U$, vdW-DF, and generalized gradient
approximations (GGA), as opposed to using local-density approximations (LDA).

###  3.2 Statistical mechanics

####  3.2.1 Cluster expansion

Figure 7: ECI values associated with each basis function in the cluster
expansion. Figure 8: Formation energies and convex hull predicted by the
cluster expansion (CLEX) are compared with those from the DFT calculations.
The configurations predicted by the cluster expansion to lie on the convex
hull are the same as those on the DFT convex hull.

In parameterizing the cluster expansion, the genetic algorithm selected 35
basis functions. The RMS error for the cluster expansion prediction of the
formation energy was 2.49 meV per primitive cell. A plot of the effective
cluster interaction (ECI) for each cluster in the cluster expansion is shown
in Figure 7. The formation energies predicted by the cluster expansion are
also plotted with the DFT calculated formation energies in Figure 8, along
with the their respective convex hulls. It can be seen the convex hull
constructed from the cluster expansion predictions matches the convex hull of
the DFT data.

####  3.2.2 Chemical potential data and phase diagram from Monte Carlo
simulations

Figure 9: Chemical potential values at 300 K, predicted using Monte Carlo
simulations without bias potentials.

Li content as a function of chemical potential and temperature was computed
through Monte Carlo statistical mechanics calculations for temperatures
ranging from 200 K to 400 K, initially without using the bias potentials.
Since the bias potentials were not used in these calculations the two-phase
regions are easily seen as gaps in the composition values of the data (see,
for example, Figure 9). Regions with orderings are recognized by the steepness
and positive slope of the curve traced out by connecting the chemical
potential data points. This indicates a strongly convex free energy well and
thermodynamic preference for these structures over the corresponding
composition intervals. Two-phase regions appear on either side of these
intervals. These data were used to construct an approximate phase diagram for
O3 LCO, for Li compositions greater than $x\geq 1/3$. The calculated phase
diagram of the LixCoO2 system is shown in Figure 10.

Figure 10: Phase diagram for the O3 structure for LCO with orderings at
$x=1/3$ and $x=1/2$, based on Monte Carlo results. Figure 11: Comparison of
the predicted voltage at 300 K (blue) with the experimental voltage (red) from
Ref. [64].

From the Monte Carlo results, the order-disorder transition temperature for
the ordering at $x=0.5$ is slightly greater than 330 K, which is consistent
with experimental results [12]. An ordering is also predicted by the Monte
Carlo results at $x=0.33$, appearing below approximately 235 K. This ordering
has not been seen experimentally at room temperature, but some reported
evidence shows that it may exist at low temperatures [14], which is consistent
with our calculations. The improvement in accuracy of the predicted order-
disorder temperatures over previous first-principles work by Van der Ven, et
al. [19] is likely due to advancements in the DFT methods over the last two
decades. A fairly good agreement between our computed voltage profile at 300 K
and the experimental measurements [64] is also observed in Figure 11. Similar
to past first-principles studies on LCO, however, because of the lack of the
metal-insulator transition that leads to the two-phase region observed between
$\sim 0.7$-$0.9$ in calculations, the voltage plateau in the same region is
not captured. We also note that the predicted voltage, as with the other Monte
Carlo results in this work, is only for the O3 structure and, therefore, only
valid for compositions greater than $x\approx 1/3$.

Additional Monte Carlo calculations were performed for temperatures of 260 K,
300 K, and 340 K, this time using bias potentials to sample within the
unstable regions of order-disorder transition. Since the fluctuations in the
chemical potential around $x=0.5$ are critical to representing the ordering in
the phase field models at temperatures below $\sim$330 K, dense sampling of
chemical potential data was performed in that region for 260 K and 300 K.
Additional sampling was also performed near the extreme compositions of $x=0$
and $x=1$ to help capture the divergent nature of the chemical potential (see
Figure 12(a)).

(a) Chemical potential, $\mathrm{d}f/\mathrm{d}c$

(b) $\mathrm{d}^{2}f/\mathrm{d}c^{2}$

(c) Free energy, $f(c)$

Figure 12: Comparison of IDNN fits with the data for 260 K, 300 K, and 340 K:
(a) the chemical potential data used for training, as sampled with the bias
potentials, and the IDNN prediction, (b) numerically differentiated data and
derivative of the chemical potential IDNN, and (c) the analytically integrated
free energy DNN.

####  3.2.3 Symmetry-adapted order parameters

Symmetry-adapted order parameters are useful for tracking the order-disorder
transition and the different ordering variants that can form. Since the Monte
Carlo and DFT results show a zig-zag ordering forming at $x=0.5$, we developed
the corresponding order parameters here. The zig-zag ordering of Li atoms on
the triangular lattice has 12 variants (3 rotations $\times$ 4 translations),
as shown in Figure 13. Following the algorithm to construct the mutually
commensurate supercell [43] yielded a supercell with 16 unique sublattice
sites. Using the 16 sublattice sites as a basis, the twelve variants can be
written as a set of twelve vectors
$\\{\boldsymbol{x}^{(1)},\ldots,\boldsymbol{x}^{(12)}\\}$.

Figure 13: The zig-zag ordering has 12 variants, resulting from combinations
of 3 rotations and 4 translations.

Additionally, 48 symmetry operations between these twelve variants were
enumerated in the group $P$, each represented by a $16\times 16$ matrix.
Following the algorithm from Thomas and Van der Ven [53], we constructed a
$P$-invariant matrix and performed its eigenvalue decomposition. It resulted
in four nonzero eigenvalues (one distinct and three repeated) and four sets of
eigenvectors, with a total of 16 eigenvectors. The eigenvectors formed the
rows of a $16\times 16$ transformation matrix $\boldsymbol{Q}$, to map a
sublattice composition vector $\boldsymbol{x}$ to a vector $\boldsymbol{\eta}$
with 16 order parameters.

The subset of order parameters relevant to a given ordering can be identified
by operating with $\boldsymbol{Q}$ on the vectors $\boldsymbol{x}^{(i)}$
describing the variants of that ordering. The superscript in
$\boldsymbol{x}^{(i)}$ indexes the variants. Components of $\boldsymbol{\eta}$
that are zero for all variants are irrelevant for describing that ordering,
and the corresponding rows can be removed from $\boldsymbol{Q}$. Multiplying
$\boldsymbol{Q}$ by the vectors describing the twelve zig-zag variants, it was
determined that seven order parameters were relevant. One was associated with
the distinct eigenvalue and corresponded to the composition averaged over all
16 sublattice sites. The other six order parameters were associated with one
of the degenerate eigenvalues which had six corresponding eigenvectors. Thus,
a reduced $7\times 16$ matrix $\widehat{\boldsymbol{Q}}$ can be used if only
the zig-zag ordering is of interest, where the first row gives the average
composition, such that $\eta_{0}=x$:

$\displaystyle\widehat{\boldsymbol{Q}}$
$\displaystyle=\frac{1}{16}\setcounter{MaxMatrixCols}{16}\begin{bmatrix}[r]1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1\\\
1&-1&1&-1&-1&1&-1&1&-1&1&-1&1&1&-1&1&-1\\\
1&1&-1&-1&-1&1&1&-1&-1&-1&1&1&1&-1&-1&1\\\
-1&1&-1&1&-1&1&-1&1&1&-1&1&-1&1&-1&1&-1\\\
1&-1&-1&1&-1&1&1&-1&1&-1&-1&1&-1&1&1&-1\\\
1&-1&-1&1&1&1&-1&-1&-1&1&1&-1&-1&-1&1&1\\\
1&1&-1&-1&-1&-1&1&1&1&1&-1&-1&-1&-1&1&1\end{bmatrix}$ (23)

Of the other two degenerate eigenvalues, one had three associated eigenvectors
and described the row ordering, and the other had six eigenvectors and
described a double-row ordering.

###  3.3 Trained integrable deep neural networks

The chemical potential IDNNs for 260 K, 300 K, and 340 K, trained as functions
of Li composition, are plotted on top of the chemical potential data in Figure
12(a), showing a fairly accurate match. Since the regions of phase separation
occur where the second derivative of the free energy is negative, the chemical
potential data were numerically differentiated and plotted with the analytical
derivatives of the chemical potential IDNNs, as shown in Figure 12(b). Phase
separation regions do indeed appear in both the Monte Carlo data, where they
are resolved with the biased potential, and in the neural network fit on
either side of the ordering at $x=0.5$ at 260 K and 300 K, although the two-
phase region at $x\approx 0.52$ for 260 K appears to occur at a slightly lower
composition in the IDNN than in the data. The analytically integrated free
energy neural networks are plotted in Figure 12(c). We also emphasize that
these data and functions represent the O3 structure and are, therefore, most
relevant for $x\geq 1/3$, since previous experimental and computational work
show that LCO adopts the H1-3 structure below $x=1/3$.

###  3.4 Phase field simulations

Here, we demonstrate our integrated scale bridging framework via three case
studies of phase field simulations that help gain some insight to the phase
stability of LCO cathodes in the composition neighborhood of the order-
disorder transition around $x=0.5$. The thermodynamic description entering the
phase field model is quantitatively accurate, since it is the result of the
scale bridging framework. Furthermore, we use appropriate kinetic data.
However, some physics remains to be included in our model, such as strain
effects, which are known to cause LCO particle degradation [65, 66, 67],
interphase interface energies, explicit representations of the different
orderings, and spatio-temporal heat generation/transport. Therefore, we do not
claim to present fully quantitative predictions, but demonstrate at least how
morphology evolves under various microstructure-specific, environmental,
kinetic and cycling conditions at a resolution not accessible in experiments
and at a computational fidelity not previously accessed with regard to the
length scales. From such simulations we make connections to observations in
the experimental literature on degradation of LCO [65, 66, 67, 68].

Figure 14: 2D phase field simulation results at 260, 300, and 340 K showing
the Li composition in a 1 $\mu$m diameter particle, with initial Li
composition randomly perturbed about $x=0.53$ and no boundary flux. The
zoomed-in phase diagram on the right identifies the phases that are observed
in the simulations.

In the first set of phase field simulations, the Li composition was
initialized as a randomly perturbed field about $x=0.53$ for each of the three
temperatures, within a circular domain with a diameter of 1 $\mu$m. No fluxes
were applied, and the composition field was allowed to evolve until a steady
state was achieved. The simulation at 340 K quickly reached a uniform
composition of $x=0.53$, consistent with the completely convex free energy at
that temperature. At 300 K and 260 K, however, spinodal decomposition took
place, followed by Ostwald ripening, as shown in Figure 14. In both
simulations, there was separation into an ordered phase and a disordered
phase. Slight differences in the compositions of these phases can be seen,
however, at the different temperatures. At 260 K, the ordered phase appears at
$x=0.512$ and the disordered phase at $x=0.551$, whereas the ordered and
disordered phases occur at $x=0.509$ and $x=0.539$, respectively, at 300 K.
This shift toward higher compositions at the lower temperature for this
particular two-phase region is qualitatively consistent with the phase diagram
in Figure 10, although slight mismatches between the nonconvex region of the
IDNN and the data, particularly at 260 K, has shifted the predicted
compositions slightly lower than the values shown in the phase diagram. The
fact that the phase boundaries occur at different compositions for different
temperatures is also seen in the ratio of ordered to disordered phase in the
transient microstructures observed at 300 and 260 K. Since the ratio of
ordered to disordered phase is smaller at 300 K, circular precipitates (in 2D)
of the ordered phase appear, whereas the ordered phase regions are elongated
and more tortuous at 260 K.

(a) 340 K

(b) 300 K

(c) 260 K

Figure 15: Phase field results showing the Li composition field resulting from
applying a cycling current density for (a) 340, (b) 300, and (c) 260 K, for a
1 $\mu$m diameter particle. The applied C-rate is set to 0.625, -1, 0.625,
-0.3125, and 0 C at 0, 500, 800, 1300, and 2000 s, respectively, where a
positive sign denotes discharging and negative is charging. The blue line
shows the applied C-rate, and the red line shows the corresponding average Li
composition of the particle.

As seen in Figure 14, notable variations can arise in the morphology of the
charged cathode particle in a temperature range that overlaps the window of
relevance for actual batteries. For example, based on the figure, battery
designers should anticipate morphology changes in charged lithium ion
batteries in an electric vehicle parked outdoors in the summer during day
_versus_ night in many geographical locations. As observed by Pender et al.
[67], lattice dimension changes in the neighborhood of $x=0.5$ cause
mechanical damage. Order-disorder transitions are pronounced over the
$260-300$ K range in Figure 14 and as the associated lattice dimension vary,
the mismatches across the order-disorder interphase interfaces will cause
stresses that may play a role in driving the damage [65]. We also note that
the higher temperature disordered phase has a more homogeneous microstructure,
and the simulation in Figure 14 predicts that it persists across the
composition range studied. This could motivate doping with elements that
stabilize the disordered structure. Such a mechanism has relevance to the Al
and La doped stabilization of the order-disorder transition as observed
recently [67]. The phase field treatment also predicts higher driving forces
for Li flux in regions of rapidly varying interface curvature, which are seen
at the lower temperatures. The lower mobility at $260$ K relative to $300$ K
slows down transport despite the correspondingly greater curvature-induced
driving force.

Additional phase field simulations show the effect of cycling of a LixCoO2
particle at 260, 300, and 340 K, as presented in Figure 15 [57]. The Li
composition was initialized to 0.5, following the recommended lower limit of
composition for LCO, below which large c-axis variations of the lattice cause
enhanced degradation [67]. This corresponds to an initial predicted voltage of
4.23 V. An initial current density of 0.1335 A m-2 was applied. The applied
current density was adjusted to -0.2137, 0.1335, -0.0668, and 0 A m-2 at 500,
800, 1300, and 2000 s, respectively. This cycling corresponds to C-rates of
0.625, -1, 0.625, -0.3125, and 0 C, where a positive sign denotes discharging
and negative is charging. This cycle results in a net discharging, as
indicated by the voltage-composition plot in Figure 11, with a final predicted
voltage of 4.14, 4.17, or 4.19 V for 260, 300, or 340 K, respectively.
Equilibrium was reached in all simulations by 2013 s. Since 340 K is above the
order-disorder transition at $x=0.5$, no ordered phase forms. However, the
presence of the order-disorder interfaces at lower temperatures suggests that
such regimes of charging/discharging of LCO particles can induce transient
misfit stresses. Such cycling between the ordered and disordered states at the
lower temperatures can promote degradation as discussed in the experimental
literature [66, 68, 67].

Figure 16: 2D phase field simulation results at 260, 300, and 340 K showing
the Li composition in a 50 nm diameter particle, with initial Li composition
randomly perturbed about $x=0.53$ and no boundary flux. Compare with the
results using 1 $\mu$m particles in Figure 14. Figure 17: Phase field results
showing the Li composition field resulting from applying a cycling current
density for a 50 nm diameter particle. Results are visually indistinguishable
for 260, 300, and 340 K, and so only one set of results is included here. The
applied current density is reduced by a factor of 20 compared with the 1
$\mu$m particles to maintain the same C-rates in both particle sizes (see
Figure 15). The blue line shows the applied C-rate, and the red line shows the
corresponding average Li composition of the particle.

While the phase field simulations in this work up to this point have used 1
$\mu$m diameter particles, experimental work by Choi et al. [66] has compared
performance using micron sized LCO particles with that of 30-70 nm particles.
Among other findings, they report that the cycle performance with these
nanometer-sized particles is consistent, without deterioration, across their
selected range of operating temperatures (258 to 333 K), and that the smaller
sized particles lead to increase reliability. With this as motivation,
spinodal decomposition and cycling phase field simulations were performed
using a 50 nm diameter particle, again with the temperatures 260, 300, and 340
K. The current density applied in the cycling simulations was scaled down by a
factor of 20 to maintain the same C-rates in the smaller particles as in the
larger, micron sized particles. The same mobility functions and phase field
gradient coefficient were used as with the larger particles.

Results from the phase field simulations initialized with Li composition
randomly perturbed about $x=0.53$, with no boundary flux, are shown in Figure
16 for the 50 nm particle. As with the larger particles, spinodal
decomposition occurs at 260 and 300 K, but not 340 K. However, equilibrium is
reached between 100 and 10,000 times more quickly with the smaller particles,
depending on the temperature. Additionally, while the equilibrium
microstructures are comparable for both sizes, the transient microstructure is
much simpler for the smaller particle. This is because the length scale of the
transient microstructure in the simulations is defined by the gradient
coefficient, which is related to the interfacial energy and is independent of
the particle size. Since the length scale of the microstructure is maintained
but the size of the particle is decreased, the result is a less complex
microstructure. This reduced complexity could lead to diminished misfit stress
effects and increased reliability.

The results from the cycling simulation for the 50 nm particles presented in
Figure 17 show that, unlike with the larger particles, the Li composition
field is essentially the same for all three temperatures. This is because the
particles are small enough that the relative values of the mobility, particle
size, and microstructure length scale lead the Li atoms to diffuse uniformly
and extremely rapidly in comparison with the time scale of the simulation such
that any temperature effects are not seen. This, again, matches the
experimental findings of Choi et al. [66] which show consistent performance
across the temperature range of 258 to 333 K.

##  4 Discussion and conclusions

In this work, we have presented a computational framework for modeling
materials physics that bridges between atomic and continuum scales, as applied
to the widely used battery material LixCoO2. The framework begins with DFT
calculations of the formation energy of potential configurations and
compositions of LixCoO2, bridging scales up through statistical mechanics
Monte Carlo simulations and integrable deep neural network training,
culminating in continuum phase field modeling. In particular, through the DFT
and Monte Carlo simulations, the framework identified the order-disorder
transition around composition $x=0.5$ at room temperature. This phase
transition was subsequently captured within the neural network representation
of the free energy and the formation and evolution of this transition was
modeled in the phase field simulations.

Our study focused on phase stability has been developed, thus far, for the O3
structure of LCO. However, since H1-3 form of LCO is the stable structure
below $x=1/3$, our results do not hold in that range. The corresponding DFT
studies are underway, and the results from them will be subject to our scale-
bridging framework to ultimately present a complete study that will have
rigorous predictive value over the entire composition range. We note that,
since LCO has two distinct free energy curves in the composition ranges
$0<x\leq 1/3$ and $1/3<x<1$, a phase field treatment would not be
thermodynamically rigorous over the entire composition range. Sharp interface
methods such as levelsets present one appropriate framework in this setting.

With regard to charge splitting we note that our DFT+U calculations used the
Hubbard correction with a van der Waals functional, and we highlighted three
configurations at $x=1/2$, which are reported in Figure 6. Of these, the
lowest energy configuration, which we adopted as the ground state, was found
to be the zigzag ordering (configuration C in Figure 6). Charge splitting was
not observed for this ground state configuration, although we noted that a
more extensive investigation would require more advanced hybrid exchange
functionals [59, 60, 61] or extended Hubbard models (DFT $+U+V$) [62, 63]. If
such detailed studies do predict charge splitting configurations, this
information could be propagated up the scales by incorporating them in the
cluster Hamiltonians and defining additional order parameters that would be
identified from the Monte Carlo sampling. The active learning approach would
play an important role in directing this sampling, which would take place in a
much higher dimensional space to generate training data for the IDNN.

The choice of the $U$ parameter was guided by matching with experimental
voltage measurements in the $x\in[0,1/2]$ and $x\in[0,1]$ composition ranges
(see Figure 3). This choice is finally reflected in the voltage reported
across the entire $x\in[0,1]$ composition range in Figure 11. We note that a
different approach to calibration could be to choose $U$ for an optimal fit
across the composition range in Figure 11, although this would be a much more
expensive and complex process that would have to fundamentally leverage the
scale bridging between DFT and statistical mechanics. We further note that
empirical approaches to determine the $U$ parameter can be completely
circumvented by using the recently developed self-consistent Hubbard method
[69] where the site specific values of $U$ are determined using density-
functional perturbation theory and iterated self-consistently until
convergence of the DFT $+U$ ground-state. This approach seamlessly ties the
$U$ values to the changing composition through its dependence on the local
chemical environment.

We re-emphasize that while the LCO material system has been the subject of a
vast body of computational study at the atomic as well as at the mesoscopic
continuum scale and thermodynamic investigations, all supported by
experimental campaigns, the extent of scale bridging has not extended as
widely as in this communication, to the best of our knowledge. Our choice of
LCO has been motivated by this deep and broad literature, which provides a
robust basis for comparison as well as the interest in this material system.
More broadly, this study also serves a demonstration of method development
that will be applied to other interesting materials systems, with systematic
incorporation of additional effects. On this last note, our first extension
will be to the non-conserved symmetry-adapted order parameters,
$\boldsymbol{\eta}$ discussed in Section 3.2.3. Other extensions are also
possible to coupling the phase field description with a thermodynamically
rigorous treatment of heat generation and transport that will more completely
reflect the phase diagram. Scale-bridged treatments of strain effects, misfit
elasticity, charge splitting as well as the dependence of coefficients on all
these fields are also accessible. An extension to the equations of
electrochemistry at the active particle-electrolyte scale will allow the study
of phenomena such as the mosaic instability observed in nanoparticulate
batteries, wherein phase segregation occurs between distinct particles rather
than within a single particle [70] as studied here. Even without such a level
of rigorously obtained detailed physics, we note that the phase field
treatment exposes trends that both align with battery material particle and
cell scale experiments, and present a platform for making technologically
relevant projections as discussed in Section 3.4.

## Author contributions

The study’s scale bridging methodology was developed by KG and GHT; electronic
structure studies were planned by all the authors; electronic structure
calculations were performed by GHT and SD; statistical mechanics and phase
field studies were carried out by GHT; the paper was written by all the
authors.

## Acknowledgements

We thank Brian Puchala for his insight regarding the CASM software and related
methods. We also gratefully acknowledge the support of Toyota Research
Institute, Award #849910, “Computational framework for data-driven,
predictive, multi-scale and multi-physics modeling of battery materials”. We
acknowledge the support of the U.S. Army Research Office through the DURIP
grant W911NF1810242, which provided the computational resources for this work.
This work has also been supported in part by National Science Foundation DMREF
grant #1729166, “Integrated Framework for Design of Alloy-Oxide Structures”.
Additional support was provided by Defense Advanced Research Projects Agency
(DARPA) under Agreement No. HR0011199002, “Artificial Intelligence guided
multi-scale multi-physics framework for discovering complex emergent materials
phenomena”. Simulations in this work were also performed in part using
resources were provided by the Extreme Science and Engineering Discovery
Environment (XSEDE) Comet at the San Diego Supercomputer Center and Stampede2
at the Texas Advanced Computing Center through allocations TG-DMR180072 and
TG-MCH200011. The phase field simulations were performed using resources
provided by the NSF via grant 1531752 MRI: Acquisition of Conflux, A Novel
Platform for Data-Driven Computational Physics (Tech. Monitor: Ed Walker),
with past and current support by the University of Michigan. XSEDE is
supported by National Science Foundation grant number ACI-1548562.

## References

  * [1] Shiyao Zheng, Chaoyu Hong, Xiaoyun Guan, Yuxuan Xiang, Xiangsi Liu, Gui-Liang Xu, Rui Liu, Guiming Zhong, Feng Zheng, Yixiao Li, Xiaoyi Zhang, Yang Ren, Zonghai Chen, Khalil Amine, and Yong Yang.  Correlation between long range and local structural changes in ni-rich layered materials during charge and discharge process.  Journal of Power Sources, 412:336–343, 2019.
  * [2] Linqin Mu, Qingxi Yuan, Chixia Tian, Chenxi Wei, Kai Zhang, Jin Liu, Piero Pianetta, Marca M. Doeff, Yijin Liu, and Feng Lin.  Propagation topography of redox phase transformations in heterogeneous layered oxide cathode materials.  Nature Communications, 9(1):2810, Jul 2018.
  * [3] Yang Yang, Rong Xu, Kai Zhang, Sang-Jun Lee, Linqin Mu, Pengfei Liu, Crystal K. Waters, Stephanie Spence, Zhengrui Xu, Chenxi Wei, David J. Kautz, Qingxi Yuan, Yuhui Dong, Young-Sang Yu, Xianghui Xiao, Han-Koo Lee, Piero Pianetta, Peter Cloetens, Jun-Sik Lee, Kejie Zhao, Feng Lin, and Yijin Liu.  Quantification of heterogeneous degradation in li-ion batteries.  Advanced Energy Materials, 9(25):1900674, 2019.
  * [4] Yuwei Mao, Xuelong Wang, Sihao Xia, Kai Zhang, Chenxi Wei, Seongmin Bak, Zulipiya Shadike, Xuejun Liu, Yang Yang, Rong Xu, Piero Pianetta, Stefano Ermon, Eli Stavitski, Kejie Zhao, Zhengrui Xu, Feng Lin, Xiao-Qing Yang, Enyuan Hu, and Yijin Liu.  High-voltage charging-induced strain, heterogeneity, and micro-cracks in secondary particles of a nickel-rich layered cathode material.  Advanced Functional Materials, 29(18):1900247, 2019.
  * [5] Wen Liu, Pilgun Oh, Xien Liu, Min-Joon Lee, Woongrae Cho, Sujong Chae, Youngsik Kim, and Jaephil Cho.  Nickel-rich layered lithium transition-metal oxide for high-energy lithium-ion batteries.  Angewandte Chemie International Edition, 54(15):4440–4457, 2015\.
  * [6] Peiyu Hou, Jiangmei Yin, Meng Ding, Jinzhao Huang, and Xijin Xu.  Surface/interfacial structure and chemistry of high-energy nickel-rich layered oxide cathodes: Advances and perspectives.  Small, 13(45):1701802, 2017.
  * [7] Un-Hyuck Kim, Seung-Taek Myung, Chong S. Yoon, and Yang-Kook Sun.  Extending the battery life using an al-doped li[ni0.76co0.09mn0.15]o2 cathode with concentration gradients for lithium ion batteries.  ACS Energy Letters, 2(8):1848–1854, 2017.
  * [8] Defei Kong, Jiangtao Hu, Zhefeng Chen, Kepeng Song, Cheng Li, Mouyi Weng, Maofan Li, Rui Wang, Tongchao Liu, Jiajie Liu, Mingjian Zhang, Yinguo Xiao, and Feng Pan.  Ti-gradient doping to stabilize layered surface structure for high performance high-ni oxide cathode of li-ion battery.  Advanced Energy Materials, 9(41):1901756, 2019.
  * [9] Tina Weigel, Florian Schipper, Evan M. Erickson, Francis Amalraj Susai, Boris Markovsky, and Doron Aurbach.  Structural and electrochemical aspects of lini0.8co0.1mn0.1o2 cathode materials doped by various cations.  ACS Energy Letters, 4(2):508–516, 2019.
  * [10] G.H. Teichert, A.R. Natarajan, A. Van der Ven, and K. Garikipati.  Machine learning materials physics: Integrable deep neural networks enable scale bridging by learning free energy functions.  Computer Methods in Applied Mechanics and Engineering, 353:201 – 216, 2019.
  * [11] G.H. Teichert, A.R. Natarajan, A. Van der Ven, and K. Garikipati.  Scale bridging materials physics: Active learning workflows and integrable deep neural networks for free energy function representations in alloys.  Computer Methods in Applied Mechanics and Engineering, 371:113281, 2020.
  * [12] Jan N. Reimers and J. R. Dahn.  Electrochemical and in situ x-ray diffraction studies of lithium intercalation in LixCoO2.  Journal of The Electrochemical Society, 139(8):2091–2097, Aug 1992\.
  * [13] M. Ménétrier, I. Saadoune, S. Levasseur, and C. Delmas.  The insulator-metal transition upon lithium deintercalation from LiCoO2: Electronic properties and 7Li NMR study.  Journal of Materials Chemistry, 9(5):1135–1140, 1999.  cited By 337.
  * [14] Y. Shao-Horn, S. Levasseur, F. Weill, and C. Delmas.  Probing lithium and vacancy ordering in O3 layered LixCoO2 (x$\approx$0.5).  Journal of The Electrochemical Society, 150(3):A366, 2003.
  * [15] Y Takahashi, N Kijima, K Tokiwa, T Watanabe, and J Akimoto.  Single-crystal synthesis, structure refinement and electrical properties of Li0.5CoO2.  Journal of Physics: Condensed Matter, 19(43):436202, Sep 2007.
  * [16] T. Motohashi, T. Ono, Y. Sugimoto, Y. Masubuchi, S. Kikkawa, R. Kanno, M. Karppinen, and H. Yamauchi.  Electronic phase diagram of the layered cobalt oxide system LixCoO2 $(0.0\leq x\leq 1.0)$.  Phys. Rev. B, 80:165114, Oct 2009.
  * [17] Zhaohui Chen and J.R. Dahn.  Methods to obtain excellent capacity retention in LiCoO2 cycled to 4.5 v.  Electrochimica Acta, 49(7):1079 – 1090, 2004.
  * [18] Keke Chang, Bengt Hallstedt, Denis Music, Julian Fischer, Carlos Ziebert, Sven Ulrich, and Hans J. Seifert.  Thermodynamic description of the layered O3 and O2 structural LiCoO2–CoO2 pseudo-binary systems.  Calphad, 41:6 – 15, 2013.
  * [19] A. Van der Ven, M. K.. Aydinol, G. Ceder, G. Kresse, and J. Hafner.  First-principles investigation of phase stability in LixCoO2.  Phys. Rev. B, 58:2975–2987, 1998.
  * [20] C. Wolverton and Alex Zunger.  First-principles prediction of vacancy order-disorder and intercalation battery voltages in LixCoO2.  Phys. Rev. Lett., 81:606–609, Jul 1998.
  * [21] Taichi Abe and Toshiyuki Koyama.  Thermodynamic modeling of the LiCoO2–CoO2 pseudo-binary system.  Calphad, 35(2):209 – 218, 2011.
  * [22] Neel Nadkarni, Tingtao Zhou, Dimitrios Fraggedakis, Tao Gao, and Martin Z. Bazant.  Modeling the metal–insulator phase transition in LixCoO2 for energy and information storage.  Advanced Functional Materials, 29(40):1902821, 2019.
  * [23] Matteo Cococcioni and Stefano de Gironcoli.  Linear response approach to the calculation of the effective interaction parameters in the $\mathrm{LDA}+\mathrm{U}$ method.  Phys. Rev. B, 71:035105, Jan 2005.
  * [24] Muratahan Aykol, Soo Kim, and C. Wolverton.  van der waals interactions in layered lithium cobalt oxides.  The Journal of Physical Chemistry C, 119(33):19053–19058, 2015\.
  * [25] Chirranjeevi Balaji Gopal and Axel van de Walle.  Ab initio thermodynamics of intrinsic oxygen vacancies in ceria.  Phys. Rev. B, 86:134117, Oct 2012.
  * [26] ME Arroyo y de Dompablo, A Van der Ven, and G Ceder.  First-principles calculations of lithium ordering and phase stability on li x nio 2.  Physical Review B, 66(6):064112, 2002.
  * [27] Vladimir I. Anisimov, Jan Zaanen, and Ole K. Andersen.  Band theory and Mott insulators: Hubbard U instead of Stoner I.  Phys. Rev. B, 44:943–954, Jul 1991.
  * [28] F. Zhou, M. Cococcioni, C. A. Marianetti, D. Morgan, and G. Ceder.  First-principles prediction of redox potentials in transition-metal compounds with $\mathrm{LDA}+U$.  Phys. Rev. B, 70:235121, Dec 2004.
  * [29] Muratahan Aykol and C. Wolverton.  Local environment dependent $\text{GGA}+u$ method for accurate thermochemistry of transition metal compounds.  Phys. Rev. B, 90:115105, Sep 2014.
  * [30] T. Thonhauser, Valentino R. Cooper, Shen Li, Aaron Puzder, Per Hyldgaard, and David C. Langreth.  Van der Waals density functional: Self-consistent potential and the nature of the van der Waals bond.  Phys. Rev. B, 76:125112, Sep 2007.
  * [31] Jiří Klimeš, David R Bowler, and Angelos Michaelides.  Chemical accuracy for the van der Waals density functional.  Journal of Physics: Condensed Matter, 22(2):022201, Dec 2009.
  * [32] D C Langreth, B I Lundqvist, S D Chakarova-Käck, V R Cooper, M Dion, P Hyldgaard, A Kelkkanen, J Kleis, Lingzhu Kong, Shen Li, P G Moses, E Murray, A Puzder, H Rydberg, E Schröder, and T Thonhauser.  A density functional for sparse matter.  Journal of Physics: Condensed Matter, 21(8):084203, Jan 2009.
  * [33] Riccardo Sabatini, Emine Küçükbenli, Brian Kolb, T Thonhauser, and Stefano de Gironcoli.  Structural evolution of amino acid crystals under stress from a non-empirical density functional.  Journal of Physics: Condensed Matter, 24(42):424209, Oct 2012.
  * [34] T. Thonhauser, S. Zuluaga, C. A. Arter, K. Berland, E. Schröder, and P. Hyldgaard.  Spin signature of nonlocal correlation binding in metal-organic frameworks.  Phys. Rev. Lett., 115:136402, Sep 2015.
  * [35] Kristian Berland, Valentino R Cooper, Kyuho Lee, Elsebeth Schröder, T Thonhauser, Per Hyldgaard, and Bengt I Lundqvist.  van der waals forces in density functional theory: a review of the vdW-DF method.  Reports on Progress in Physics, 78(6):066501, May 2015.
  * [36] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, Andrea Dal Corso, Stefano de Gironcoli, Stefano Fabris, Guido Fratesi, Ralph Gebauer, Uwe Gerstmann, Christos Gougoussis, Anton Kokalj, Michele Lazzeri, Layla Martin-Samos, Nicola Marzari, Francesco Mauri, Riccardo Mazzarello, Stefano Paolini, Alfredo Pasquarello, Lorenzo Paulatto, Carlo Sbraccia, Sandro Scandolo, Gabriele Sclauzero, Ari P Seitsonen, Alexander Smogunov, Paolo Umari, and Renata M Wentzcovitch.  Quantum espresso: a modular and open-source software project for quantum simulations of materials.  Journal of Physics: Condensed Matter, 21(39):395502 (19pp), 2009\.
  * [37] P Giannozzi, O Andreussi, T Brumme, O Bunau, M Buongiorno Nardelli, M Calandra, R Car, C Cavazzoni, D Ceresoli, M Cococcioni, N Colonna, I Carnimeo, A Dal Corso, S de Gironcoli, P Delugas, R A DiStasio Jr, A Ferretti, A Floris, G Fratesi, G Fugallo, R Gebauer, U Gerstmann, F Giustino, T Gorni, J Jia, M Kawamura, H-Y Ko, A Kokalj, E Küçükbenli, M Lazzeri, M Marsili, N Marzari, F Mauri, N L Nguyen, H-V Nguyen, A Otero de-la Roza, L Paulatto, S Poncé, D Rocca, R Sabatini, B Santra, M Schlipf, A P Seitsonen, A Smogunov, I Timrov, T Thonhauser, P Umari, N Vast, X Wu, and S Baroni.  Advanced capabilities for materials modelling with quantum espresso.  Journal of Physics: Condensed Matter, 29(46):465901, 2017.
  * [38] Andrea Dal Corso.  Pseudopotentials periodic table: From H to Pu.  Computational Materials Science, 95:337 – 350, 2014.
  * [39] Pslibrary 1.0.0.  https://github.com/dalcorso/pslibrary.
  * [40] K. Momma and F. Izumi.  VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data.  J. Appl. Crystallogr., 44:1272–1276, 2011.
  * [41] B. Meredig, A. Thompson, H. A. Hansen, C. Wolverton, and A. van de Walle.  Method for locating low-energy solutions within $\text{DFT}+u$.  Phys. Rev. B, 82:195128, Nov 2010.
  * [42] M. K. Aydinol, A. F. Kohan, G. Ceder, K. Cho, and J. Joannopoulos.  Ab initio study of lithium intercalation in metal oxides and metal dichalcogenides.  Phys. Rev. B, 56:1354–1365, Jul 1997.
  * [43] Anirudh Raju Natarajan, John C. Thomas, Brian Puchala, and Anton Van der Ven.  Symmetry-adapted order parameters and free energies for solids undergoing order-disorder phase transitions.  Phys. Rev. B, 96:134204, Oct 2017.
  * [44] J. M. Sanchez, F. Ducastelle, and D. Gratias.  Generalized cluster description of multicomponent systems.  Physica A, 128:334–350, 1984.
  * [45] D. De Fontaine.  Cluster approach to order-disorder transformations in alloys.  volume 47 of Solid State Physics, pages 33 – 176. Academic Press, 1994.
  * [46] https://github.com/prisms-center/CASMcode.  CASM: A Clusters Approach to Statistical Mechanics, v0.3.dev, 2018.
  * [47] A. Van der Ven, J. C. Thomas, Q. Xu, and J. Bhattacharya.  Linking the electronic structure of solids to their thermodynamic and kinetic properties.  Mathematics and Computers in Simulation, 80(7):1393–1410, 2010\.
  * [48] John C. Thomas and Anton Van der Ven.  Finite-temperature properties of strongly anharmonic and mechanically unstable crystal phases from first principles.  Physical Review B, 88(21):214111–214111, December 2013.
  * [49] B. Puchala and A. Van der Ven.  Thermodynamics of the Zr-O system from first-principles calculations.  Phys. Rev. B, 88:094108, Sep 2013.
  * [50] Y. Mishin.  Atomistic modeling of the $\gamma$ and $\gamma^{\prime}$-phases of the Ni–Al system.  Acta Materialia, 52(6):1451 – 1467, 2004.
  * [51] Babak Sadigh, Paul Erhart, Alexander Stukowski, Alfredo Caro, Enrique Martinez, and Luis Zepeda-Ruiz.  Scalable parallel monte carlo algorithm for atomistic simulations of precipitation in alloys.  Phys. Rev. B, 85:184203, May 2012.
  * [52] Babak Sadigh and Paul Erhart.  Calculation of excess free energies of precipitates via direct thermodynamic integration across phase boundaries.  Phys. Rev. B, 86:134204, Oct 2012.
  * [53] John C. Thomas and Anton Van der Ven.  The exploration of nonlinear elasticity and its efficient parameterization for crystalline materials.  Journal of the Mechanics and Physics of Solids, 107:76–95, October 2017.
  * [54] J. W. Cahn and J. E. Hilliard.  Free energy of a nonuniform system. I Interfacial energy.  The Journal of Chemical Physics, 28:258–267, 1958.
  * [55] S. M. Allen and J. W. Cahn.  A microscopic theory for antiphase boundary motion and its application to antiphase boundary coarsening.  Acta Metallurgica, 27:1085–1091, 1979.
  * [56] A. Van der Ven.  Lithium diffusion in layered LixCoO2.  Electrochemical and Solid-State Letters, 3(7):301, 2000.
  * [57] Tonghu Jiang, Shiva Rudraraju, Anindya Roy, Anton Van der Ven, Krishna Garikipati, and Michael L. Falk.  Multiphysics simulations of lithiation-induced stress in Li1+xTi2O4 electrode particles.  The Journal of Physical Chemistry C, 120(49):27871–27881, 2016\.
  * [58] Daniel Arndt, Wolfgang Bangerth, Denis Davydov, Timo Heister, Luca Heltai, Martin Kronbichler, Matthias Maier, Jean-Paul Pelteret, Bruno Turcksin, and David Wells.  The deal.II library, version 8.5.  Journal of Numerical Mathematics, 25(3), April 2017.
  * [59] Axel D. Becke.  Density‐functional thermochemistry. III. The role of exact exchange.  The Journal of Chemical Physics, 98(7):5648–5652, 1993.
  * [60] John P. Perdew, Matthias Ernzerhof, and Kieron Burke.  Rationale for mixing exact exchange with density functional approximations.  The Journal of Chemical Physics, 105(22):9982–9985, 1996.
  * [61] Jochen Heyd, Gustavo E. Scuseria, and Matthias Ernzerhof.  Hybrid functionals based on a screened coulomb potential.  The Journal of Chemical Physics, 118(18):8207–8215, 2003.
  * [62] A. Amaricci, A. Camjayi, K. Haule, G. Kotliar, D. Tanasković, and V. Dobrosavljević.  Extended hubbard model: Charge ordering and wigner-mott transition.  Phys. Rev. B, 82:155102, Oct 2010.
  * [63] Vivaldo Leiria Campo Jr and Matteo Cococcioni.  Extended DFT+ U+ V method with on-site and inter-site electronic interactions.  Journal of Physics: Condensed Matter, 22(5):055602, 2010.
  * [64] G. G. Amatucci, J. M. Tarascon, and L. C. Klein.  CoO2, the end member of the LixCoO2 solid solution.  Journal of The Electrochemical Society, 143(3):1114–1123, Mar 1996\.
  * [65] Haifeng Wang, Young-Il Jang, Biying Huang, Donald R Sadoway, and Yet-Ming Chiang.  Tem study of electrochemical cycling-induced damage and disorder in licoo2 cathodes for rechargeable lithium batteries.  Journal of the Electrochemical Society, 146(2):473, 1999.
  * [66] Sun Hee Choi, Ji-Won Son, Young Soo Yoon, and Joosun Kim.  Particle size effects on temperature-dependent performance of licoo2 in lithium batteries.  Journal of power sources, 158(2):1419–1424, 2006.
  * [67] Joshua P Pender, Gaurav Jha, Duck Hyun Youn, Joshua M Ziegler, Ilektra Andoni, Eric J Choi, Adam Heller, Bruce S Dunn, Paul S Weiss, Reginald M Penner, et al.  Electrode degradation in lithium-ion batteries.  ACS nano, 14(2):1243–1295, 2020.
  * [68] Feng Leng, Cher Ming Tan, and Michael Pecht.  Effect of temperature on the aging rate of li ion battery operating above room temperature.  Scientific reports, 5(1):1–12, 2015.
  * [69] Iurii Timrov, Nicola Marzari, and Matteo Cococcioni.  Self-consistent hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formulations.  Phys. Rev. B, 103:045141, Jan 2021.
  * [70] Bernardo Orvananos, Todd R Ferguson, Hui-Chia Yu, Martin Z Bazant, and Katsuyo Thornton.  Particle-level modeling of the charge-discharge behavior of nanoparticulate phase-separating li-ion battery electrodes.  Journal of The Electrochemical Society, 161(4):A535, 2014.