# LSTM-based approach to detect cyber attacks on market-based congestion
management methods

Omniyah Gul M Khan1, Amr Youssef2, Ehab El-Saadany34, and Magdy Salama5
1Electrical and Computer Engineering, University of Waterloo, Canada, Email:
<EMAIL_ADDRESS>2Information Systems Engineering, Concordia University,
Canada, Email<EMAIL_ADDRESS>3Advanced Power and Energy Center,
EECS Department, Khalifa University, Email<EMAIL_ADDRESS>4Electrical and Computer Engineering, University of Waterloo, Email:
<EMAIL_ADDRESS>5Electrical and Computer Engineering, University of
Waterloo, Email<EMAIL_ADDRESS>
###### Abstract

Market-based congestion management methods adopt Demand Side Management (DSM)
techniques to alleviate congestion in the day-ahead market. Reliance of these
methods on the communication layer makes it prone to cyber attacks affecting
the security, reliability, and economic operation of the distribution network.
In this paper, we focus on Load Altering Attacks that would compromise the
operation of market-based congestion management methods. A detection technique
is proposed using Long Short-term Memory (LSTM) Recurrent Neural Networks
(RNN). IEEE 33 bus system is used as a case study to demonstrate the
effectiveness of the proposed technique. An accuracy of 97% was obtained
proving the capability of using LSTM-RNN to detect a load altering cyber
attack compromising aggregators in the network.

###### Index Terms:

Cyber attack, Load Altering Attack, Long Short-Term Memory (LSTM), Recurrent
Neural Network (RNN).

##  I Introduction

Electrification of the transportation system and space heating, and the
increased penetration of Distributed Energy Resources have changed the
operating conditions of the distribution system. Distribution System Operators
(DSO) are facing congestion problems due to power flow exceeding the network’s
transfer capability. To manage congestion that could cause voltage violations
and/or thermal overloading, possibly damaging devices such as distribution
transformers and feeders, market-based congestion management methods have been
developed. Such methods utilize existing flexibility in the market through
Demand Side Management (DSM), avoiding incurring huge costs to reinforce
network assets [1].

Market-based congestion management methods [2] [3] employ the communication
layer to relay consumer preferences to their respective aggregators, price
tariffs from DSO to aggregators, and load schedules from aggregators to the
DSO to mitigate congestion. The success of such congestion management methods,
hence, relies on the cyber security of the entities involved and their
communication links. Ensuring the confidentiality, integrity, and availability
of the DSO, aggregators, consumers’ meters, and the information exchanged is
vital to resolve congestion using DSM.

Cyber attacks compromising consumers’ smart meters and their communication
lines have been studied in the literature. [4] studied the impact of closed-
loop load altering attacks compromising the stability of the power system. In
[5], three types of attackers were studied, with varying cyber attacking
capabilities, that aim to steal electricity using false data injection (FDI)
by tampering with the smart meters. On the other hand, in [6], the impact of a
cyber attack manipulating the readings of consumers’ meters responsible for
monitoring renewable DG units to overcharge the utility.

Attacking aggregators or their communication links with the DSO is another
aspect that can be utilized by a cyber attacker to affect the performance of
market-based congestion management methods. Aggregators hold a unique position
in the distribution network in terms of their connection to their respective
DER equipment and the DSO [7]. Compromising the cyber security of aggregators
could risk the overall security of the grid. Financial losses can be caused as
a result of an adversary causing congestions not to be detected in the day-
ahead market by the DSO. This would cause unnecessary load shedding that could
be avoided if the attack was mitigated. Moreover, cyber attacks can be
executed by DSO competitors to fake congestion resulting in higher congestion
tariffs to be imposed on the consumers. This would encourage consumers to
change their utility provider. Hence, it is critical to study the impact of
cyber attacks on aggregators and how to mitigate or detect such attacks.
However, no literature exists on studying attacks compromising aggregators
explicitly except our work in [8]. In [8], the impact of compromising
aggregators, as a result of Load Altering Attacks (LAA), on congestion and
congestion tariff has been evaluated. But, up to the author’s knowledge, no
literature exists on detection of such attacks that compromise aggregators
load profiles communicated to the DSO affecting congestion in the network.

Motivated by the research gap, the main contribution of this paper is focused
on developing a detection scheme that would determine the presence or absence
of a load altering cyber attack compromising aggregators in the distribution
network. The detection scheme is executed by the DSO in the day-ahead market
on receiving aggregators’ load profiles prior to computing the congestion
tariff.

##  II Related Work

Load Altering Attacks (LAAs) has been previously studied in the literature. In
[9], static LAAs compromising Internet connected vulnerable loads were
studied. A cost-efficient load protection was then proposed identifying the
vulnerable loads that need to be secured without overloading the distribution
system. In [4], closed-loop dynamic LAA were studied in which non-secured load
was cyber-attacked over time in a pre-programmed trajectory based on system
frequency. A non-convex pole placement optimization problem is then solved to
determine the minimum load required to be secured without impacting power
system stability. It should be noted, however, that [9]-[4] were studied based
on power system dynamics.

Detection of FDI attacks utilizing data-based methods has been also studied.
Data-based methods’ ability to scale to larger systems with low computational
cost makes it an attractive technique [10]. [11] tested and compared the
performance of linear and Gaussian Support Vector Machines (SVM), K-Nearest
Neighbours (KNN), and single layer perceptrons as FDI detection models.
Similarly, [12] utilized SVM, KNN, and Extended Nearest Neighbours to detect
the most accurate method to detect falsification of measurement and state
data. A SVM-based classifier was adopted to detect energy theft attacks as a
result of hacking smart meters in [13].

##  III Case Study System Details

To simulate load altering attacks on market-based congestion management
techniques and the proposed detection model, the IEEE 33 bus system [14],
illustrated in Figure 1, is adopted. Having 32 load buses, residential
consumers are represented using both flexible and base loads. Residential base
load profiles are designed to have a maximum load equivalent to the default
loading of the IEEE 33 test bus system. On the other hand, residential
flexible loads consist of one Heat Pump (HP) and one Electric Vehicle (EV)
with a battery size of 36kWh, per residence. Consumers’ preferences in terms
of their comfortable indoor temperature range and initial State of Charge
(SOC) of EV are randomized to be uniformly distributed between
$18-24\degree\text{C}$ and $20-30\%$ respectively. The residential charging
power of EVs was chosen as 11kW. All the load buses of the IEEE 33 bus system,
are assumed to be residential, having both flexible and non-flexible loads,
except bus 23 and 24. Having the largest base-load, load bus 23 and 24 were
utilized as commercial buses consisting of Fast Charging Direct Current (FCDC)
EV charging parking lots of 50kW. Arrival and departure of EVs from the FCDC
parking lot were simulated using data obtained from Toronto parking authority
[15]. The absence of an EV from the commercial parking lot was assumed to
indicate its presence at its residence. All residential and commercial buses
were represented by 9 aggregators having comparable loads, as illustrated in
Figure 1. Aggregators 6 and 7 are used to represent commercial buses. Day-
ahead electricity prices are obtained from the PJM market via the dataminer
tool [16]. GAMS [17] was used for executing the optimizations in this paper.
MATLAB was used for determining the congestion tariffs imposed by the DNO to
relieve congestion and for implementing the proposed detection scheme.

Figure 1: IEEE 33 bus system involving nine aggregators [14]

##  IV Impact of LAAs compromising aggregators

The general methodology adopted in market-based congestion management methods
is illustrated in Figure 2. Three entities are involved: DSO, aggregators, and
consumers. DSO is the central system operator that is responsible for ensuring
the absence of network congestion. It utilizes the Wide Area Network (WAN) to
receive optimized load schedules from the aggregators, performs load flow
analysis, and then determines the congestion tariff to be imposed to alleviate
congestion. These tariffs are then sent to the aggregators to re-optimize
their load schedules based on their customers’ preferences. Aggregators are
responsible for receiving their customers’ preferences using the Neighborhood
Area Network (NAN) and optimizing their load schedules to minimize their
costs. Consumers are the owners of flexible and non-flexible loads and they
hire aggregators to represent their needs in the electricity market. Utilizing
the Home Area Network (HAN), consumers transmit their loads and their
preferences to their respective aggregators.

Figure 2: General representation of DSM-based congestion management

DSM-based congestion management methods’ reliance on the communication layer
to relieve congestion makes them prone to cyber attacks. A Load Altering cyber
Attack (LAA) can be implemented by compromising aggregators and altering the
load schedules that are transmitted to the DSO. The optimal load schedule of
aggregators can be defined as follows,

$\begin{split}P_{i,t}^{a}=E(P_{c,t}^{a,nf}+P_{c,t}^{a,f}),\forall a\in
N_{a},\forall i\in N_{d},\\\ \forall c\in N_{m},\forall t\in N_{T}\\\
\end{split}$ (1)

where $P_{i,t}^{a}\in\mathbb{R}^{N_{d}}$ is the load of bus $i$ of aggregator
$a$ at time $t$. $N_{a}$, $N_{d}$ and $N_{m}$ are the set of aggregators,
buses, and customers in the network. $N_{T}$ is the total set of time slots
considered in the day-ahead market. $E$ is the customer-to-bus mapping matrix
where $E\in\mathbb{R}^{{N_{d}}\times{N_{m}}}$. $P_{c,t}^{a,nf}$ and
$P_{c,t}^{a,f}$ are the non-flexible and flexible loads of customer, $c$, of
aggregator $a$ at time $t$, where
$\\{P_{c,t}^{a,nf},P_{c,t}^{a,f}\\}\in\mathbb{R}^{N_{m}}$. The load schedule
of compromised aggregator can be defined as follows,

$\begin{split}P_{i,t}^{a^{c}}=P_{i,t}^{a^{c}}+\Delta_{i}^{a^{c}}\text{, for
}0\leq\Delta_{i}^{a^{c}}\leq\gamma\end{split}$ (2)

where $a^{c}$ is the aggregator that is compromised and $\Delta_{i}^{a^{c}}$
represents the alteration of the compromised aggregator’s load. To not raise
any alarms with the DSO, $\Delta_{i}^{a^{c}}$ needs to be as low as possible.
$\gamma$ is the upper limit for the change in the compromised aggregators’
load schedule. In this paper, $\gamma$ is chosen as 0.1 to ensure the
stealthiness of the attack while achieving the objective of the attacker.

Adopting the case study explained in Section II, and utilizing the DDT
congestion management method [18] to mitigate congestion, the load schedules
of aggregators 3 and 7 are compromised to study the impact of the cyber
attack. Figure 3 illustrates the effect of LAA on power flow in line 0-1 of
the IEEE 33 bus system as a comparison between the attacked and not attacked
scenarios, prior to imposing any congestion tariff. At peak time, $t=12$, the
line 0-1 was already congested. However, the LAA caused an increase in the
congestion level. Moreover, line 0-1 was not congested at times 10 and 11.
However, on scaling the load schedules of the two aggregators, congestion was
created. This has an impact on electricity price as illustrated in Figure 4.
In the absence of a LAA, the imposed congestion tariff increased the price of
electricity at t = 12 by 15%. However, due to the load altering attack, the
imposed tariff increased the price by a not needed 10% at t = 12 which
customers have to pay for a higher congestion level that is not real. Also, as
a result of fake congestion being created at times 10 and 11, a congestion
tariff is imposed by the DNO to solve the fake congestion. Customers now have
to pay 4.4% more at t = 10 and 6.4% more at t = 11.

Figure 3: Effect of a LAA on power flow levels Figure 4: Effect of LAA on
electricity price

##  V Modeling of LAA Detection Scheme

To detect LAAs, Recurrent Neural Networks (RNN), which is a class of
Artificial Neural Networks (ANN), is proposed in this paper. Traditional
backpropagation NNs, though capable of learning nonlinear functions, cannot
capture sequential information in the input data which is required for dealing
with sequence data. Neglecting data’s temporal correlation impacts the model’s
accuracy [19]. RNNs, on the other hand, have a recurrent connection in their
hidden states making them capable of capturing sequential information present
in input data. The network, hence, develops a memory of past events and
encodes it in its hidden state. RNNs are ideal for applications in which there
is a temporal dependency in the data [20]. Hence, an RNN’s output at $t$ could
have an effect on its decision at $t+1$. There are many types of RNNs based on
their input and output data: sequence input (also known as many-to-one),
sequence output (also known as one-to-many), sequence input and output (also
known as many-to-many), or synced sequence inputs and outputs (also known as
bidirectional many-to-many). RNNs are trained using backpropagation. However,
for long sequences, the gradients, which carry information utilized for RNN
parameter update, measured in the last layers becomes smaller. Thus, as a
consequence parameter updates become insignificant. LSTM RNNs is an
improvement over RNNs since it reduces RNNs vanishing gradient problem. LSTM’s
insensitivity to the gap between the relevant information in the past and the
time when that information is needed in the present gives it an advantage over
alternative sequential models, such as Hidden Markov Models [21].

Being in the day-ahead market, and receiving the aggregators’ load schedules
as a full day sequence makes utilization of LSTM as a cyber attack detection
scheme attractive for the DSO prior to computation of the congestion tariff.
To model the detection scheme, the steps illustrated in Figure 5, and
explained in the following subsections, needs to be executed.

Figure 5: LSTM-based LAA Detection Model

###  V-A Data set, Feature Extraction and Preprocessing

To model the LSTM network, hourly customer real power load, electricity price,
and weather data for Michigan, USA was utilized. Data was collected for random
days in November to March of 2015 to 2020. These winter months were chosen as
the loads are at their peak, and there is a high chance of the network being
congested. The real power load per bus in the distribution network is
determined using the weather and day-ahead electricity prices to optimize
aggregators’ consumers’ load. Hence, the four types of sequences that were
chosen as input features are as follows:

  * •

Historical daily weather data.

  * •

Historical day-ahead electricity prices [16].

  * •

Three days historical real power load per bus.

  * •

Power load per bus.

  * •

Power flow in the main feeder of the distribution network.

All the data were normalized using min-max normalization. For every feature,
the minimum value was transformed to a 0, the maximum value was changed to a
1, and every other value gets represented by a decimal between 0 and 1. 1550
days of data were used to simulate attacked and not attacked scenarios. An
attacked day is labeled as 1, and 0 otherwise.

###  V-B LSTM Modelling

To model LAA LSTM-based detection scheme, a many-to-one LSTM model is adopted.
The input to the model is in the form of sequential data representing the
24-hour weather forecast, day-ahead electricity prices, real power load per
bus, and the power flow in the feeder of concern. The output to the model, on
the other hand, is a classification decision. A decision equal to $1$
represents an attacked scenario, while a decision equal to $0$ represents a
not-attacked scenario. The LSTM-RNN model includes multiple LSTM-RNN memory
units that replaces each node in the hidden layer of a neural network model.
Each unit consists of an input gate, $i(t)$, an output gate, $o(t)$, a forget
gate, $f(t)$, and a memory cell state, $c(t)$, as illustrated in Figure 6
[22]. The three gates are responsible for selecting and rejecting the
information that passes through the network in its cell state. Initially, the
LSTM needs to determine which information it will remove from the cell state,
$c(t-1)$. This is done by the sigmoid layer forget gate, $f(t)$, which takes
its previous cell output, $h(t-1)$, and the new cell input, $x(t)$, to output
a number between 0 and 1 for $c(t-1)$. A 1 indicates retaining the contents of
$c(t-1)$ while a 0 means deleting the contents.

Figure 6: Structure of the LSTM-RNN memory-unit

$f(t)=\delta(W_{f}[x(t),h(t-1)]+b_{f})\\\ $ (3)

A sigmoid layer input gate then takes the current input, $x(t)$, and previous
cell output, $h(t-1)$, to determine which cell states to update. A tanh layer
is responsible for creating an array of new updated values that would be added
to the state.

$i(t)=\delta(W_{i}[x(t),h(t-1)]+b_{i})\\\ $ (4)
$\tilde{c}(t)=tanh(W_{c}[x(t)+h(t-1)+b_{c})\\\ $ (5)

The previous cell state, $c(t-1)$, is then updated by multiplying it with
$f(t)$ (3) to forget the contents that was decided to be forgotten. The result
is added to the product of $i(t)$ (4) and $\tilde{c}(t)$ (5) to update the
cell state values.

$\tilde{c}(t)=f(t)\times{c}(t-1)+i(t)\times\tilde{c}(t)\\\ $ (6)

The output gate, $o(t)$, utilizes the sigmoid activation function to decide
which part of the current cell makes it to the output. These values are then
filtered using a tanh function to be between -1 and 1.

$o(t)=\delta(W_{o}[x(t),h(t-1)]+b_{o})\\\ $ (7) $h_{t}=o_{t}\tanh{c(t)}\\\ $
(8)

$W_{f}$, $W_{i}$, $W_{c}$, and $W_{o}$ represent the weight matrices for the
inputs of the network activation functions. Also, $b_{f}$, $b_{i}$, $b_{c}$,
and $b_{o}$ represent the biases for the different gates. The above process is
repeated for all time steps, as illustrated in Figure 7. Weights and biases
are updated by minimizing the differences between the LSTM outputs and the
training samples.

Figure 7: Overall Structure of the LSTM-RNN Model

###  V-C Results and Discussion

$70\%$ of the data explained in the section earlier, corresponding to 1200
days of data representing attacked and not attacked scenarios, were used for
training and validating the LAA detection scheme. The remaining $30\%$ of the
data, corresponding to 350 days, were used for testing the performance of the
model. the performance of the LSTM detection scheme was compared with the
traditional back propagation neural network (ANN) model. To evaluate its
performance, confusion matrix and accuracy were used as performance metrics. A
confusion matrix assesses the performance of the model in terms of true
positives, true negatives, false positives, and false negatives. False
positives refer to the number of cases that were classified as attacked when
in reality they were not attacked. On the other hand, a false negative refers
to the number of cases that were not classified as attacked by the model. True
positives and true negatives refer to case scenarios representing attacked and
not attacked days respectively and the model classified it correctly. On the
other hand, accuracy represents the correctness of the classifier and is
calculated as follows:

$\begin{split}\text{Accuracy }=\frac{\text{ (True Positive + True
Negative)}}{\text{Number of Cases}}\end{split}$ (9)

Table I illustrates the confusion matrix of the results obtained using the
LSTM and the ANN models. As observed, in both the models the number of false
positive is higher compared to that of the false negative cases. The ANN
resulted had nearly $9.4\%$ of false positive cases compared to $1.7\%$ of
false negatives. Similarly, LSTM resulted in nearly $2.9\%$ of false positives
compared to $0.6\%$ of false negatives. Using more data to train the model on
more attack scenarios can decrease this number so that the model could clearly
distinguish between an attacked and not attacked scenario. Compared to the ANN
model, the LSTM model performed better. Summing up the number of true positive
and true negative cases, as illustrated in (9), results in an accuracy of
$96.57\%$ for the LSTM model and accuracy of $88.86\%$ for the ANN model. The
results signify the advantage of using LSTM sequential classifier to model the
LAA detection scheme.

In the day-ahead market, on receipt of the aggregators’ load schedules, the
DSO executes the proposed LSTM model to detect if an attack has occurred. If
the LSTM classifier detects no attack, then the DSO would go ahead to compute
the congestion tariffs. However, if an attack is detected the DSO would inform
the aggregators that a cyber attack has occurred in the network so that the
aggregators could check their cyber security. Until the cyber attack is
eliminated, the DSO would rely on its historical data to determine if there is
congestion in the network. Future work should incorporate the detection of the
aggregator that is compromised in the network.

Table I: Confusion Matrix of LAA Detection Models | ANN | LSTM
---|---|---
Output  Actual | Attack |  No Attack | Attack |  No Attack
Attack | 137 | 33 | 141 | 10
No Attack | 6 | 174 | 2 | 197

##  VI Conclusions

Cyber vulnerabilities of market-based congestion management techniques were
studied in this paper. The impact of load altering attacks that creates fake
congestion was evaluated in terms of congestion level and congestion tariffs.
A detection technique using LSTM recurrent neural networks was proposed. On
receiving aggregators’ load schedules in the day-ahead market, the DSO would
simulate the LSTM detection model to determine if any of the data received for
each bus has been compromised. Distinctive sequences were chosen as features
for the detection model to classify the presence or absence of an attack in
the day-ahead market. Although the attacks simulated were made to be stealthy
and hence hard to be detected, an accuracy of $96.57\%$ was obtained
reflecting the capability of the proposed model in detecting load altering
cyber attacks that compromises aggregators in the network.

## References

  * [1] N. Haque, A. Tomar, P. Nguyen, and G. Pemen, “Dynamic tariff for day-ahead congestion management in agent-based LV distribution networks,” _Energies_ , vol. 13, 2020.
  * [2] S. Huang, Q. Wu, L. Cheng, and Z. Liu, “Optimal reconfiguration-based dynamic tariff for congestion management and line loss reduction in distribution networks,” _IEEE Transactions on Smart Grid_ , vol. 7, pp. 1295–1303, May 2016.
  * [3] Y. Gu, J. Xie, X. Chen, K. Yu, Z. Chen, and Z. Li, “Dynamic tariff method for congestion management in distribution networks,” in _IEEE Conf. on Energy Internet and Energy System Integration_ , 2017, pp. 1–6.
  * [4] S. Amini, F. Pasqualetti, and H. Mohsenian-Rad, “Dynamic load altering attacks against power system stability: Attack models and protection schemes,” _IEEE Transactions on Smart Grid_ , vol. 9, pp. 2862–2872, July 2018.
  * [5] X. Liu, P. Zhu, Y. Zhang, and K. Chen, “A collaborative intrusion detection mechanism against false data injection attack in advanced metering infrastructure,” _IEEE Transactions on Smart Grid_ , vol. 6, pp. 2435–2443, 2015.
  * [6] M. Ismail, M. F. Shaaban, M. Naidu, and E. Serpedin, “Deep learning detection of electricity theft cyber-attacks in renewable distributed generation,” _IEEE Trans. on Smart Grid_ , vol. 11, pp. 3428–3437, 2020\.
  * [7] C. Lai, N. Jacobs, S. Hossain-Mckenzie, C. Carter, P. Cordeiro, I. Onunkwo, and J. Johnson, “Cyber security primer for der vendors, aggregators, and grid operators,” 2017.
  * [8] O. G. M. Khan, E. El-Saadany, K. Saleh, M. Shaaban, and A. Youssef, “Cyber attacks on distributed congestion management methods,” in _2019 IEEE Power Energy Society General Meeting (PESGM)_ , 2019, pp. 1–5.
  * [9] A. Mohsenian-Rad and A. Leon-Garcia, “Distributed internet-based load altering attacks against smart power grids,” _IEEE Transactions on Smart Grid_ , vol. 2, no. 4, pp. 667–674, 2011.
  * [10] J. Sakhnini, H. Karimipour, and A. Dehghantanha, “Smart grid cyber attacks detection using supervised learning and heuristic feature selection,” in _IEEE Conf. on Smart Energy Grid Engineering_ , 2019, pp. 108–112.
  * [11] M. Ozay, I. Esnaola, F. T. Yarman Vural, S. R. Kulkarni, and H. V. Poor, “Machine learning methods for attack detection in the smart grid,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 27, pp. 1773–1786, 2016.
  * [12] J. Yan, B. Tang, and H. He, “Detection of false data attacks in smart grid with supervised learning,” in _2016 International Joint Conference on Neural Networks (IJCNN)_ , 2016, pp. 1395–1402.
  * [13] P. Jokar, N. Arianpoo, and V. C. M. Leung, “Electricity theft detection in ami using customers’ consumption patterns,” _IEEE Transactions on Smart Grid_ , vol. 7, no. 1, pp. 216–226, 2016.
  * [14] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” _IEEE Transactions on Power Delivery_ , vol. 4, no. 2, pp. 1401–1407, April 1989.
  * [15] E. A. Rezai, M. Shaaban, E. F. El-Saadany, and F. Karray, “Online intelligent demand management of plug-in electric vehicles in future smart parking lots,” _IEEE Systems Journal_ , vol. 10, pp. 483–494, 2016.
  * [16] “Data miner.” [Online]. Available: https://www.pjm.com/markets-and-operations/etools/data-miner-2.aspx
  * [17] G. D. Corporation, “General algebraic modeling system (gams) release 24.2.1,” Washington, DC, USA, 2013.
  * [18] S. Huang, Q. Wu, H. Zhao, and C. Li, “Distributed optimization based dynamic tariff for congestion management in distribution networks,” _IEEE Transactions on Smart Grid_ , 2018.
  * [19] J. Li, D. Deng, J. Zhao, D. Cai, W. Hu, M. Zhang, and Q. Huang, “A novel hybrid short-term load forecasting method of smart grid using MLR and LSTM NN,” _IEEE Trans. on Industrial Informatics_ , pp. 2443–2452, 2021.
  * [20] F. M. Bianchi, E. Maiorino, M. C. Kampffmeyer, A. Rizzi, and R. Jenssen, _Recurrent Neural Networks for Short-Term Load Forecasting: An Overview and Comparative Analysis_. Springer, 2017\.
  * [21] P. Lea, _IoT and Edge Computing for Architects: Implementing edge and IoT systems from sensors to clouds with communication systems, analytics, and security_. Computers, 2020.
  * [22] R. S. F. Gers and F. Cummins, “Learning to forget: continual prediction with lstm,” _Neural Comput._ , vol. 12, pp. 2451–2471, 2000.


††thanks: These two authors contributed equally.††thanks: These two authors
contributed equally.

# Modeling intercalation chemistry with multi-redox reactions by sparse
lattice models in disordered rocksalt cathodes

Peichen Zhong Department of Materials Science and Engineering, University of
California, Berkeley, California 94720, United States  Materials Sciences
Division, Lawrence Berkeley National Laboratory, California 94720, United
States  Fengyu Xie Department of Materials Science and Engineering, University
of California, Berkeley, California 94720, United States  Materials Sciences
Division, Lawrence Berkeley National Laboratory, California 94720, United
States  Luis Barroso-Luque Department of Materials Science and Engineering,
University of California, Berkeley, California 94720, United States  Materials
Sciences Division, Lawrence Berkeley National Laboratory, California 94720,
United States  Liliang Huang Materials Sciences Division, Lawrence Berkeley
National Laboratory, California 94720, United States  Gerbrand Ceder
<EMAIL_ADDRESS>Department of Materials Science and Engineering,
University of California, Berkeley, California 94720, United States  Materials
Sciences Division, Lawrence Berkeley National Laboratory, California 94720,
United States

###### Abstract

Modern battery materials can contain many elements with substantial site
disorder, and their configurational state has been shown to be critical for
their performance. The intercalation voltage profile is a critical parameter
to evaluate the performance of energy storage. The application of commonly
used cluster expansion techniques to model the intercalation thermodynamics of
such systems from ab-initio is challenged by the combinatorial increase in
configurational degrees of freedom as the number of species grows. Such
challenges necessitate efficient generation of lattice models without over-
fitting and proper sampling of the configurational space under charge balance
in ionic systems. In this work, we introduce a combined approach that
addresses these challenges by (1) constructing a robust cluster-expansion
Hamiltonian using the sparse regression technique, including
$\ell_{0}\ell_{2}$-norm regularization and structural hierarchy; and (2)
implementing semigrand-canonical Monte Carlo to sample charge-balanced ionic
configurations using the table-exchange method and an ensemble-average
approach. These techniques are applied to a disordered rocksalt oxyfluoride
Li1.3-xMn0.4Nb0.3O1.6F0.4 (LMNOF) which is part of a family of promising
earth-abundant cathode materials. The simulated voltage profile is found to be
in good agreement with experimental data and particularly provides a clear
demonstration of the Mn and oxygen contribution to the redox potential as a
function of Li content.

###### pacs:

##  I Introduction

The market for electric-vehicle high-energy-density Li-ion batteries has
witnessed tremendous growth over the past decade [1, 2]. The increasing demand
for electrical energy storage requires further development of the high-energy-
density-based cathode materials in rechargeable Li-ion batteries. Current
cathode materials are mostly limited to layered Li(Ni,Mn,Co)O2 (NMC) variants.
The high cost of Ni and Co limits the large-scale expansion of Li-ion
batteries with NMC-type cathodes. Recently developed disordered rocksalt with
Li-excess (DRX) are promising earth-abundant cathode materials, which can
enable scaling of Li-ion energy storage to several TWh/year production.
Although these compounds have no long-range cation order, the interactions
between species generate short-range ordering (SRO), which critically affects
the electrochemical performance [3]. The broad chemical flexibility and the
wide variety of chemical environments that can be created for Li and the
transition metals (TMs) that can be used by SRO and chemistry provide new
opportunities to improve the cathode performance. Examples include improving
cyclability via fluorine/vacancy doping of the anion sublattice [4, 5],
enhancing the rate capability by engineering the cation SRO [3], and achieving
zero-strain cathodes for solid-state batteries [6, 7, 8].

Over the last two decades, many properties of Li-ion battery materials have
been successfully predicted by first-principles calculations [9, 10, 11]. The
equilibrium voltage is one of the fundamental quantities that help
characterize the electrochemical performance of a particular material and is
defined by the difference in Li chemical potentials between the cathode and
anode:

$V=-\frac{\mu_{\text{Li}}^{\text{cathode}}-\mu_{\text{Li}}^{\text{anode}}}{zF}.$
(1)

In Eq. (1), $z$ is the charge transferred per ion, $F$ is Faraday’s constant,
and $\mu_{\text{Li}}$ is the chemical potential of Li. For example,
considering a Li transition-metal oxide LixTMO2 ($x_{1}\leq x\leq x_{2}$) as
the cathode and Li metal as the anode with the cell reaction

$\text{Li}_{x_{1}}\text{TMO}_{2}\longrightarrow\text{Li}_{x_{2}}\text{TMO}_{2}+(x_{1}-x_{2})\text{Li},$
(2)

the approximated equilibrium voltage can be computed as [12]

$\bar{V}(x_{1},x_{2})\approx-\frac{E_{\text{Li}_{x_{1}}\text{TM}\text{O}_{2}}-E_{\text{Li}_{x_{2}}\text{TM}\text{O}_{2}}-(x_{1}-x_{2})E_{\text{Li}}}{F(x_{1}-x_{2})}.$
(3)

The internal energy of the bcc Li metal $E_{\text{Li}}$, the lithiated
structure $E_{\text{Li}_{x_{1}}\text{M}\text{O}_{2}}$, and the delithiated
structure ($E_{\text{Li}_{x_{2}}\text{M}\text{O}_{2}}$) can be obtained from
first-principles density functional theory (DFT). In this approach, the
entropic effect is assumed to be small at low temperatures, and the change in
internal energy is used to approximate the chemical potential change. By
computing the formation energy of LixTMO2 structures with varied Li
concentrations $x$, a convex hull can be constructed from the energy of
ground-state structures at each concentration. A piece-wise voltage profile
can then be built from Eq. (3) via the ground states on the convex hull by
using relevant constructive values of $x$ on the hull.

Figure 1: An overview of reported methods for computing voltage profiles in
(a) ordered and (b) disordered electrode materials. The green arrows represent
the configurational samplings that can be accelerated by using cluster
expansion as an effective Hamiltonian fitted from DFT calculations.

Unlike simple Li transition-metal oxides (LiTMO2), a typical DRX cathode
(Li1+xM’aM”bO2-yFy) is composed of three major components: (1) the redox-
active species M’, which provides electron redox; (2) the inert high-valent TM
M”, which charge compensates for the Li excess and stabilizes disordered
structures [13]; (3) fluorine, which can improve the cyclability and allows
more Li excess to be accommodated without losing TM redox by lowering the
overall anion valence [14]. When the cathode is charged/discharged, Li is
removed/inserted into the cathode structure accompanied by
oxidation/reduction, resulting in various oxidation states among the redox-
active metal M’ and oxygen atoms. Different oxidation states of a TM can
exhibit very different local chemistry preferences (e.g., Mn3+ has a
substantial Jahn–Teller effect compared to Mn2+/4+). To capture these chemical
differences in simulations, different valence states of the same elements must
be treated as distinct species. This treatment is called _charge decoration_ ,
which has been demonstrated to be essential in capturing the electronic
entropy effect to construct the correct phase diagram in some compounds [15].
Charge decoration intrinsically increases the number of components, and
therefore the complexity of modeling the intercalation voltage profiles of
DRXs.

###  I.1 Intractability of composition enumeration

To obtain the voltage profile of a DRX material, most previous studies have
used the convex-hull construction approach by finding the ground states (GS),
which we refer to as the _GS-algo_. In this approach, one tries to find the
low-energy structures at varied Li content $x$ using a variety of algorithms
[16, 17]. These low-energy configurations are calculated by DFT to construct
the piece-wise voltage profile following Eq. (3). The GS-algo has shown
reasonable predictions for voltages and redox mechanisms [18, 19, 20, 17, 21,
16]. When a high number of components and valence states are present, the GS-
algo can become impractical as all possible valence combinations at each stage
of delithiation must be enumerated. For example, when evaluating a delithiated
supercell of composition Li${}_{21}\square_{18}$Mn12Nb9O48F12, the combination
of valence in Mn and O can take Mn${}^{3+}_{6}$Mn${}^{4+}_{6}$,
Mn${}^{2+}_{1}$Mn${}^{3+}_{4}$Mn${}^{4+}_{7}$,
Mn${}^{2+}_{2}$Mn${}^{3+}_{2}$Mn${}^{4+}_{8}$, and even
Mn${}^{3+}_{7}$Mn${}^{4+}_{5}$O${}^{-}_{1}$, etc. Enumerating all the possible
compositions and searching for the possible ground states under each charge-
decorated composition are NP-hard problems and become intractable, in
particular when the supercell grows. To resolve the enumeration problem, Monte
Carlo (MC) sampling is a better choice for studying configurational energetics
in a high-dimensional space.

###  I.2 Fast growth of cluster basis caused by charge decoration

To bridge the gap between 0K ground states and sampling at finite
temperatures, MC simulation with a cluster expansion (CE) as an effective
Hamiltonian is typically used for intercalation chemistry in ordered cathode
materials [22, 23, 24]. The CE casts the energy as a function of the occupancy
of atoms on a set of predefined sites. For example, in a multicomponent
system, the energy is expanded as

$E(\boldsymbol{\boldsymbol{\sigma}})=\sum_{\beta}m_{\beta}J_{\beta}\left\langle\Phi_{\boldsymbol{\alpha}\in\beta}\right\rangle_{\beta}+\frac{E_{0}}{\varepsilon_{r}},~{}\Phi_{\boldsymbol{\alpha}}=\prod_{i=1}^{N}\phi_{\alpha_{i}}(\sigma_{i}).$
(4)

A configuration $\boldsymbol{\sigma}$ represents a specific occupancy state of
species on all the system sites, where $\sigma_{i}$ describes which species
sits on the $i$-th site of the lattice. The cluster basis function
$\Phi_{\boldsymbol{\alpha}}=\prod_{i=1}^{N}\phi_{\alpha_{i}}(\sigma_{i})$ are
the product of site basis functions $\phi_{\alpha_{i}}(\sigma_{i})$ across a
collection $\alpha$ of multiple sites. They are taken as the average over the
crystal symmetry orbits $\beta$, forming a complete basis to expand the scalar
energy function on the configuration space. The expansion coefficients
$J_{\beta}$ are called effective cluster interactions (ECIs). The
electrostatic energy (Ewald energy $E_{0}/\varepsilon_{r}$) is also included
to capture long-range electrostatic interactions [25] ($E_{0}$ is the
unscreened electrostatic energy, and $1/\varepsilon_{r}$ is fitted as one of
the ECIs ($1/\varepsilon_{r}\geq 0$). We refer readers to Ref. [26] for a more
comprehensive description of the CE formalism in ionic systems.

Figure 2: An illustration of cluster basis growth: number of correlation
functions vs. number of components included in the CE with
pair/triplet/quadruplet cutoff radius of 7/4/4 Å based on a rocksalt primitive
cell with lattice parameter $a=3$ Å.

Since the cluster site basis functions are defined by the number of components
[27], the charge decoration can significantly increase the model complexity of
the CE, resulting in a rapid growth in the number of clusters. Figure 2
illustrates how the number of cluster basis functions grows as the number of
cation species included in the DRX increases, where the cutoff radius is fixed
with pair/triplet/quadruplet interactions up to 7/4/4 Å. For the full
description of Mn and O redox, more than 800 ECIs are predefined and need to
be fitted, whereas the number of DFT calculations is typically much smaller
than that of predefined ECIs. When building a CE, the ECIs
$\boldsymbol{J}^{*}=\\{J_{\beta}\\}$ can be obtained by fitting the DFT energy
$\boldsymbol{E}_{\text{DFT},S}$ of training structures $S$ to their
correlation functions with (regularized) linear regression:

$\boldsymbol{J}^{*}=\arg\min_{\boldsymbol{J}}~{}||\boldsymbol{E}_{\text{DFT},S}-\boldsymbol{\Pi}_{S}\boldsymbol{J}||^{2}_{2}+\rho(\boldsymbol{J}),$
(5)

where $\boldsymbol{\Pi}_{S}$ is the feature matrix formed by the correlation
functions and $\rho(\boldsymbol{J})$ is a regularization term. High-component
CE can easily be overfitted as the rank of the feature matrix
($\boldsymbol{\Pi}_{S}$) is typically smaller than the dimension of ECIs
dim($\boldsymbol{J}$). The low-rank structure (referred to as _rank
deficiency_) of the feature matrix $\boldsymbol{\Pi}_{S}$ in Eq. (5) requires
selecting the most physically informative ECIs and avoiding overfitting in CE
fitting [28, 29]. The rank deficiency can also arise from other physical
constraints in ionic systems. For example, charge balance creates a linear
dependency between the number of charge-decorated species and the
corresponding correlation functions. And another challenge is that the
training structures calculated with DFT are predominantly low in energy, and
such low-energy structures often narrow the configurational space that can be
represented in the training structures, which in principle, high-energy
configurations could be included. DFT tends not to cover such configurations
or relax them to lower-energy configurations by moving ions and electrons. The
number of non-zero ECIs must be constrained
($||\boldsymbol{J}||_{0}\leq\text{rank}(\boldsymbol{\Pi}_{S})$) to prevent ECI
fitting from being an underdetermined problem. This constraint can be achieved
by properly introducing $\mu||\boldsymbol{J}||_{0}$ as a regularization term
in Eq. (5) to penalize the number of non-zero ECIs and impose sparsity [30].

###  I.3 Charge-neutrality constraint in MC sampling

After fitting the ECIs, MC simulations can be used to sample the energy of
configurations under finite temperatures. Applying semigrand-canonical Monte
Carlo (sGCMC) sampling is most suited for calculating voltage profiles [31,
32]. The relation between the Li content $x$ and a Li chemical potential can
be obtained by applying the Metropolis–Hastings algorithm with the Boltzmann
distribution

$f(E(\boldsymbol{\sigma}),\mu)\propto\exp\left(-\frac{E(\boldsymbol{\sigma})-\mu_{\text{Li}}\cdot
x_{\text{Li}}N}{k_{B}T}\right).$ (6)

In Eq. (6), $E$ is the energy of the configuration $\boldsymbol{\sigma}$ given
by the CE Hamiltonian, $\mu_{\text{Li}}$ is the Li chemical potential,
$x_{\text{Li}}$ is the Li content in the configuration $\boldsymbol{\sigma}$,
$N$ is the total number of Li and vacancy sites, $k_{B}$ is the Boltzmann
constant, and $T$ is the simulation temperature. Computing voltage profiles
using sGCMC has been successful in the study of several simple binary
electrode materials, such as LiCoO2/LiNiO2 [32, 33, 34], MgTiS2 [35], and
Li3V2O5 [36].

When charge-decorated CEs are used, the requirement of charge neutrality must
be enforced in MC sampling. Because the training set only includes charge-
balanced structures with no information about the charge-unbalanced
structures, the energy predicted by the CE will be unphysical if
configurations with non-zero net charge are assessed during sGCMC. Techniques
for enforcing strict charge neutrality in sGCMC have been applied to several
electrolyte systems [37, 38]; however, few have been demonstrated in a system
with complex redox reactions such as DRX [39].

To overcome all the abovementioned issues, we propose a voltage-calculation
framework that combines several state-of-the-art methods in CE-MC. With this
framework, we demonstrate how to correctly model the intercalation voltage
profile in DRX and, more generally, any complex ionic systems with redox-
active ions and configurational disorder. In the following methodology
sections, we will introduce the construction of a robust and predictive
cluster-expansion Hamiltonian with $\ell_{0}\ell_{2}$-norm regularized sparse
regression [30], demonstrate an effective sampling strategy of the
intercalation stages with sGCMC under charge balance using the table-exchange
(TE) method [39], and illustrate an ensemble average method over
representative structures to handle various chemical environments. In the
results section, the equilibrium voltage profile of Li1.3-xMn0.4Nb0.3O1.6F0.4
(LMNOF) is presented and compared with the experiments. To explain the redox
mechanism, we analyze the proportion of multiple redox-active species at
varied Li content. We find that the calculated voltage profile and redox
mechanism agree well with experiments and argue that the ability of our method
to describe oxygen redox in the Li-excess cathode accurately is particularly
noticeable.

##  II Methodology

###  II.1 Training structure generation for DFT

To describe DRX materials well, two configurational degrees of freedom need to
be accurately represented: the Li/vacancy interactions and the different local
chemical environments (i.e., the SRO of TMs and anions). We propose the
following two-step procedure for generating the training set:

(1) _Pristine states_ : Generate several fully lithiated structures with
different transition-metal and anion configurations in relatively small
supercells (e.g., supercell structures with $10\times$ or $20\times$ the
formula unit of Li1.3Mn0.4Nb0.3O1.6F0.4). These structures can be generated
from canonical MC samplings using a pre-fitted cluster expansion as the
effective Hamiltonian or using solely the electrostatic energy for simplicity.

(2) _Delithiated states_ : Starting from the structures generated in Step 1,
fix the TM and O/F orderings and enumerate different Li/vacancy configurations
at varied Li contents (e.g., $x=0.3/0.5/0.7$ in Li1.3-x). As the total number
of enumerated structures can be large, one can further sort the structures at
each Li content by their electrostatic energy and only keep the low-energy
ones.

All the sampled structures are in the chemical space of
Li1.3-xMn0.4Nb0.3O1.6F0.4, and this two-step procedure covers different
Li/vacancy orderings in varied local chemical environments formed by TM and
anion SRO to be calculated with DFT.

###  II.2 Sparse regression for charge-decorated CE

Figure 3: (a) An illustration of the rank deficiency in the ECI fitting. The
left empty bar represents the dimension of predefined ECIs with
dim($\boldsymbol{J}$) = 858. The middle empty bar indicates the number of
structures sampled in the feature matrix with
dim$(\boldsymbol{\Pi}_{S})=(463\times 858)$. The blue shaded area overlapping
the middle bar indicates a feature matrix of rank$(\boldsymbol{\Pi}_{S})=287$.
The right solid blue bar represents the number of non-zero ECIs after an
$\ell_{0}\ell_{2}$-norm regularized sparse regression, giving
$||\boldsymbol{J}||_{0}=169$. (b) An illustration of the TEs used in charge-
balanced sGCMC. (c)–(d) The distribution of total on-site magnetizations of Mn
(c) and O (d) atoms among all DFT-r2SCAN calculated structures in the chemical
space of Li1.3-xMn0.4Nb0.3O1.6F0.4. The valence of each Mn and O atom is
determined by the on-site magnetization. From the histogram, the
classification boundary between Mn4+/3+ and Mn3+/2+ is estimated to be
$3.25\mu_{B}$ and $4.1\mu_{B}$, and the O- classification is estimated to be
$0.5\mu_{B}$. (As the percentage of O2- with a low magnetization is too large
compared to the percentage of O-, the panel (d) is truncated on the y-axis.)

To obtain effective valence states of redox-active Mn and O species from a DFT
calculated configuration, the on-site magnetization can be used [26, 40]. For
example, Figure 3(c) and (d) show the distribution of magnetic moments
representing Mn2+/3+/4+ and O2-/- in our set of 463 DFT-calculated structures.
The valence of each Mn and O atom is classified by the site magnetization
using 3.25$\mu_{B}$ for distinguishing Mn4+/3+, 4.1$\mu_{B}$ for separating
Mn3+/2+, and 0.5$\mu_{B}$ to indicate O- [41] ($\mu_{B}$ is the Bohr
magneton). For the sparse regression of ECIs, we apply the
$\ell_{0}\ell_{2}$-norm regularization with hierarchy constraints [30]. The
ECIs are optimized in the following mixed-integer quadratic programming
problem:

$\displaystyle\min_{\boldsymbol{J}}~{}\boldsymbol{J}^{T}\boldsymbol{\Pi}_{S}^{T}\boldsymbol{\Pi}_{S}\boldsymbol{J}^{T}-2\boldsymbol{E}_{\text{DFT}}^{T}\boldsymbol{\Pi}_{S}\boldsymbol{J}+\mu_{0}\sum_{c\in\boldsymbol{C}}z_{0,c}+\mu_{2}||\boldsymbol{J}||_{2}^{2}$
(7) $\displaystyle\begin{array}[]{r@{~}l<EMAIL_ADDRESS>Mz_{0,c}~{}&\geq J_{c},~{}\forall c\in\boldsymbol{C}\\\ Mz_{0,c}~{}&\geq-
J_{c},~{}\forall c\in\boldsymbol{C}\\\ z_{0,b}~{}&\leq z_{0,a},~{}\forall
a\subset b,~{}\\{a,b\\}\in\boldsymbol{C}\\\ z_{0,c}~{}&\in\\{0,1\\},~{}\forall
c\in\boldsymbol{C},\end{array}$ (12)

where $\boldsymbol{\Pi}_{S}$ is the feature matrix, $\boldsymbol{J}$ are the
ECIs, $z_{0,c}$ is the slack variable representing $J_{c}=0,z_{0,c}=0$, and
$J_{c}\neq 0,z_{0,c}=1$. $M=100$ is set to constrain the optimization
boundaries, and $||\boldsymbol{J}||_{2}^{2}=\boldsymbol{J}^{T}\boldsymbol{J}$
is a ridge regression term ($\ell_{2}$-norm). We refer the readers to Ref.
[30] for a detailed description of this approach which we found to lead to
relatively sparse but accurate CE.

The CE Hamiltonian was constructed with pair interactions up to 7 Å, triplet
interactions up to 4 Å, and quadruplet interactions up to 4 Å based on a
rocksalt primitive cell with lattice parameter $a=3$ Å leading to a possible
858 ECIs (including the constant term $J_{0}$). The ECIs were fitted using 463
training structures, forming a feature matrix of
rank$(\boldsymbol{\Pi}_{S})=287$. The resulting ECIs using the sparse
regression in Eq. (7) contain 169 non-zero elements
($||\boldsymbol{J}||_{0}=169$). The relationship between dimension, rank and
the number of non-zero elements is illustrated in Fig. 3(a).

###  II.3 Charge-balanced Monte Carlo sampling

Figure 4: An illustration of charge-balanced sGCMC on an ensemble of
disordered structures (Green: Li, white: vacancy, other colors: different TM
species, anions are not displayed). An ensemble of fully lithiated disordered
structures is first generated and referred to as the pristine state. The sGCMC
is performed to topotactically delithiate each pristine structure at
decreasing Li chemical potentials $\mu_{\text{Li}}$ (i.e., increasing voltages
$V$). The content of species is averaged over each sGCMC sampling and all the
pristine structures.

We applied sGCMC simulations on the Li/vacancy occupancy and the charge
decoration degrees of freedom to obtain the voltage curve of
Li1.3-xMn0.4Nb0.3O1.6F0.4 composition. The delithiation is achieved by step-
wise removal of Li atoms. In each MC step, the Li+ is removed/inserted,
accompanied by the oxidation/reduction of an Mn or O atom. The net charge of
each configuration is maintained at zero by only executing a combination of
site occupancy changes that are charge neutral. This type of MC step is
referred to as table exchange (TE) [39, 37]. In our calculations, we used the
following three elemental classes of perturbations:

  1. 1.

Li+ \+ Mn${}^{2+}\rightarrow$ Mn3+ \+ Vac.

  2. 2.

Li+ \+ Mn${}^{3+}\rightarrow$ Mn4+ \+ Vac.

  3. 3.

Li+ \+ O${}^{2-}\rightarrow$ O- \+ Vac.

Any other charge-conserving MC step can be expressed as a linear combination
of these three classes and their inverses. For example, charge transfer such
as $2\text{Mn}^{3+}\rightarrow\text{Mn}^{2+}+\text{Mn}^{4+}$ and
$\text{Mn}^{4+}+\text{O}^{2-}\rightarrow\text{Mn}^{3+}+\text{O}^{-}$ can be
achieved by a combination of elementary perturbations. As the intercalation is
assumed to be topotactic, Mn, O, and F ions do not change sites. The
acceptance probabilities of each MC step are scaled to ensure detailed balance
(see Appendix).

In a series of sGCMC simulations, the chemical potentials are scanned between
two limiting values $\mu_{\text{Li}}\in[\mu_{\min},\mu_{\max}]$ at finite
temperature. For each sGCMC simulation with a given Li chemical potential
$\mu$, the content of each charge-decorated species is averaged over the
sampled MC structures from the equilibrium. The Li chemical potential can be
converted into the cathode voltage using $V=-(\mu-E_{\text{Li}})/e$.

In this work, the numerical simulations were implemented with the open-source
software package smol [42].

###  II.4 The ensemble average method

Unlike ordered intercalation compounds which are characterized by a small
number of local environments, the configurational disorder in DRX creates an
abundance of local chemical environments, leading to a significant variation
of Li extraction energy from different sites [43, 44]. This multitude of
environments is difficult to capture in a single small unit cell. Instead, we
sample the delithiation from multiple distinct structures and obtain the true
lithiation curve as the ensemble average of them.

We use a canonical MC simulation to generate an ensemble of fully lithiated
(i.e., pristine) configurations $\\{\boldsymbol{\sigma}_{i}\\}$, from which
multiple representative structures are drawn. The CE model used to generate
the canonical MC structures is reported in Ref. [45]. To recover the actual
SRO in DRX, we require that the average of physical quantities in the
collection of representative structures be approximately equal to the actual
ensemble average, assuming infinitely many structures are drawn from the
ensemble
($\langle\boldsymbol{\Pi}\rangle_{\text{ens}}\approx\langle\boldsymbol{\Pi}\rangle_{\infty}=\bar{\boldsymbol{\Pi}}_{\text{SRO}}$).
For simplicity, in this work, we verify that the average energy $\langle
E\rangle_{\text{ens}}=\langle\boldsymbol{\Pi}\rangle_{\text{ens}}\cdot\boldsymbol{J}\approx\bar{\boldsymbol{\Pi}}_{\text{SRO}}\cdot\boldsymbol{J}$
converges to the ensemble average with an increasing number of selected
structures.

Each simulated fully lithiated structure contains 120 atoms
(Li39Mn12Nb9O48F12) with a supercell lattice constant of $\sim$ 10 Å. We
rationalize the choice of such a supercell size because $\sim$ 10 Å is a good
cut-off to maintain enough distance between each atom and its periodic images
and to encapsulate all cluster-interaction distances. Subsequently, sGCMC
simulations are performed for every structure in the ensemble, and for a given
chemical potential, the species (e.g., Li) content is computed using the
following average

$\langle
x_{\text{Li}}^{\mu}\rangle=\frac{1}{N_{\text{ens}}}(\bar{x}_{1}+\bar{x}_{2}+...),$
(13)

where $\\{\bar{x}_{1},\bar{x}_{2}...\\}$ are the averaged Li content given by
the thermally equilibrated sGCMC sample starting from each structure. In this
work, we used 30 disordered structures for the ensemble average.

##  III Simulation results

Figure 5: (a) Calculated (circle) and experimental (solid line) voltage
profiles of (Li1.3-xMn0.4Nb0.3O1.6F0.4). (b) Calculated content of Mn2+, Mn3+,
Mn4+ and O- per f.u. as a function of Li content ($x$). (c) Comparison of
voltage profiles generated with a charge-decorated CE (green dots) and a CE
without charge decoration (red dots). The shaded area represents the
distribution of the derivative quantity $|\frac{dx}{dV}|$, which reflects the
Li-site energy distribution during the intercalation process.

We demonstrate the procedures presented above to model the intercalation
thermodynamics of Li1.3-xMn0.4Nb0.3O1.6F0.4 (LMNOF). Figure 5(a) presents the
simulated and the experimental voltage profiles [45]. The solid green line is
the experimental charging profile under a low current density (20 mA/g)
between 1.5 and 5.0 V. The bottom of the discharge state (1.5 V) is aligned to
be the fully lithiated state (Li1.3Mn0.4Nb0.3O1.6F0.4). The green dots
represent the computed voltage profile at T = 300 K. The slope and turning
point in the slope are in good agreement with the experimental profile for
$0.4\leq x\leq 1.3$, indicating that the Li/vacancy interaction and redox
potentials of Mn and O are well described in our model. In the highly charged
region ($x<0.7$) specifically, the simulation shows remarkably good agreement
with the experiment reproducing the fact that the slope of the profile becomes
flatter at $x\sim 0.7$ (marked by the gray dashed line in Fig. 5(a)). The
computed voltage profile is systematically lower than the experimental one,
which is well-known for most DFT functionals that are not augmented with a
Hubbard $U$ correction [46]. Even though the SCAN functional removes more
self-interaction than previous LDA and GGA [47], it does not fully remove it
and underestimates the intercalation voltage.

The fraction of each redox-active species during intercalation is presented in
Fig. 5(b). Mn2+/3+/4+ are represented by orange/green/blue dots, and O- are
represented by red dots. Comparing Fig. 5(a) and (b) makes it apparent that
the Li content where the voltage profile flattens ($x\sim 0.7$) corresponds to
the start of the oxidation of O2- to O-. The lowering of the voltage slope as
the system changes from TM redox to O redox is consistent with the higher
dilution of the oxygen charge compared to TM redox centers. Figure 5(b) also
reveals several key points about the redox mechanism in LMNOF. The oxidation
of Mn does not appear consecutively from Mn2+ to Mn3+ to Mn4+. The amount of
Mn3+ and Mn4+ increase simultaneously as Mn2+ begins to be oxidized. The fact
that different oxidation states of Mn co-exist over a wide range of Li content
is likely due to the variety of local chemical environments induced by cation
disorder. The co-existence is consistent with the marginal stability of Mn3+
and its propensity to disproportionate into Mn2 and Mn4+ when it cannot exist
in an environment where it can lower its energy significantly through a
Jahn–Teller distortion [48]. O-redox occurs after all the Mn2+ has been
consumed but before all the Mn atoms are fully oxidized to +4 valence. The
hybridized redox mechanism between Mn and O has been confirmed by previous
synchrotron characterization experiments [45] and on related materials
Li1.2Mn0.4Ti0.4O2.0 by mass spectroscopy [49].

##  IV Discussion

Typical modeling of the intercalation energetics only includes the Li/vacancy
degree of freedom, assuming that the electronic degrees of freedom can be
integrated out. Integrating out degrees of freedom is based on two key
assumptions [50]: (1) One assumes that the degree of freedom that is
variationally removed is always optimized in the DFT calculation. More
specifically, this would require that for a given Li/vacancy configuration,
the DFT calculation can easily find the charge decoration with the lowest
energy. While this is a reasonable assumption for systems with a highly
delocalized charge such as LixCoO2 [51], it is unlikely to be the case for
disordered systems where the variation of local environments and the highly
localized charge easily lead to charge metastability in DFT. (2) The second
assumption when degrees of freedom are variationally removed is that their
entropy contribution can be neglected, as their variational optimization is
supposed to find their ground state. The role of electronic entropy in
intercalation systems is yet to be fully understood, but at least in LixFePO4
it has been shown to be critical to reproduce the correct intercalation
behavior [15]. Both assumptions are unlikely to be valid in DRX cathodes,
given their multitude of possible redox couples and large variations of
energetic environments.

To illustrate the significance of accurately capturing the electronic degree
of freedom in DRX materials, we fitted another CE parameterizing only the
Li/vacancy interactions. Mn and O atoms are considered single charge-less
species regardless of their oxidation states in such a CE. The computed
voltage profiles with the charge-decorated CE (in green) and the undecorated
CE (in red) are compared in Fig. 5(c). The undecorated CE yields a featureless
voltage profile and virtually no change of slope as a function of Li content.
The mechanism-related details are poorly portrayed in the profile. The shaded
areas in Fig. 5(c) are the relative derivative of capacity with respect to
voltage ($dx/dV$). Two distinct peaks are observed in the charge-decorated
calculation (green), indicating the contribution from the Mn-redox and
O-redox, whereas they cannot be adequately distinguished in the undecorated
version (red).

Figure 6: (a) Average energy (per atom) among pristine structures as a
function of the number of pristine structures selected in an ensemble. The
red-shaded region indicates a variation of $\pm 1$ meV/atom of $\langle
E\rangle_{\infty}$. The orange/green lines represent the results of the
supercell structure with 120/960 atoms, respectively. (b) Comparison of
computational efficiency (defined in Eq. (C.1)) between supercells with 120
and 960 atoms, computed and averaged over all the pristine structures in the
ensembles.  Figure 7: Contents of each redox species per formula unit as a
function of the Li content $x$ in LixMn0.4Nb0.3O1.6F0.4 as computed from sGCMC
simulations at T = 300 K. The statistics are computed based on an ensemble of
disordered structures of in supercells of Li30xMn12Nb9O48F12. The average
values over different disordered structures are marked with colored dots, the
median values are marked with orange lines, and the variances are marked with
error bars around the medians.

We also highlight the importance of using the ensemble average method to
capture the multitude of environments in disordered systems. In principle, one
can approach the representation of disordered systems by using a single
structure in a very large supercell (e.g., special quasi-random structure
(SQS) approach [52]) or by making an ensemble of multiple smaller-sized
supercell structures and taking an average of the computed quantities. These
two approaches are statistically equivalent, given that enough structures have
been used with the smaller-sized supercells. Figure 6(a) shows the number of
the structures required for the average energy per atom over structures to
convergence within $\pm 1$ meV/atom of the equilibrium average $\langle
E\rangle_{\infty}$ when choosing from an ensemble of representative structures
in canonical MC [45]. The green line represents the results in a supercell of
960 atoms, whereas the orange line represents the results in a supercell of
120 atoms. The smaller-sized supercell requires drawing $\sim$ ten times more
structures to converge, with the number of structures required roughly scaling
with the supercell size ($960=120\times 8)$. A large supercell approaches the
actual distribution of SRO with fewer structures. However, a much longer MC
sampling time is required in a large supercell, whereas in a smaller
supercell, one can easily benefit from parallelizing multiple pristine
structures to require much less total time consumption. To illustrate this
factor, we introduce computational efficiency as a quantitative benchmark for
sGCMC. Figure 6(b) shows the computational efficiency in sGCMC at varied
voltages in supercells with 120 (orange) and 960 (green) atoms. Higher
efficiency indicates that less computational time is required (see details
about the definition in the Appendix). Approximating the computational
efficiency, the sampling time required in each pristine structure scales $\sim
O(N^{1})$ with the supercell size. Therefore, with the help of
parallelization, using an ensemble of smaller structures is statistically as
accurate but practically more efficient than using fewer large structures. An
ensemble of smaller structures is also more tractable in DFT for computing
other properties when necessary, such as electronic structure, whereas DFT is
computationally prohibitive in large supercells even when the SRO can be well
presented.

Finally, we like to discuss the necessity of using enough structure samples
when smaller-sized supercells are used to model disordered systems. Figure 7
shows the variance of redox-active species with respect to the averaged Li
content. As a result of the disorder, the variance in the content of each
species among pristine structures is not negligible. This further discourages
the methods based on convex-hull and direct DFT evaluations (e.g., GS-algo) to
accurately determine voltage profiles for materials with significant disorder
such as DRX. When the size of the supercell in which one enumerates possible
Li/vacancy and charge configurations is too small or the number of structure
samples is too few, as is often the case in GS-algo, the variance resulting
from different local chemical environments cannot be adequately captured.

##  V Summary & Outlook

In DRX materials, the abundance of available redox-active transition metals
creates new opportunities for cathode design, such as using the redox
reactions between $\text{Cr}^{3+}/\text{Cr}^{4+}/\text{Cr}^{6+}$ [53, 54],
$\text{V}^{3+}/\text{V}^{4+}/\text{V}^{5+}$ [55], or even
$\text{Fe}^{2+}/\text{Fe}^{3+}/\text{Fe}^{4+}$ [56]. However, the high number
of chemical components, distinguished by their valence states in DRX materials
renders the “curse of dimensionality” (CoD) the main simulation obstacle [57].
As the complexity of the energy model grows with the number of components and
charge decoration, the computational cost grows exponentially fast. For
example, CoD has been an essential impediment for the computational design of
high-entropy cathodes [6, 54].

In this work, we propose the following procedure to obtain accurate
intercalation voltage profiles in DRX with multi-redox reactions: (1) build a
training set containing different fully lithiated structures and enumerate the
Li/vacancy orderings at varied delithiation levels calculated by DFT; (2)
construct a charge-decorated CE Hamiltonian and fit the ECIs using the sparse
regression technique to address the fast growth of the cluster basis and rank-
deficiency issues (e.g., $\ell_{0}\ell_{2}$-norm regularized regression); and
(3) run sGCMC in an ensemble of disordered structures under charge balance to
obtain physically rational energetics and compute the voltage profiles.

As demonstrated in this work, the workflow above provides an effective way to
study the intercalation chemistry of oxides with Mn and O redox activity.
Looking forward, the recent development of machine-learning force-fields
(MLFF) may provide new opportunities to accelerate the training structure
generation process by accurately approximating DFT calculations [58, 59, 60,
61]. The CE Hamiltonian can be further obtained by coarse-graining the MLFF
predicted configurational energy, especially using the charge-informed MLFF to
include the heterovalent states of transition metals [60]. We believe such an
approach has the potential to bridge up first-principles calculations, force
fields, and cluster expansions and give a higher accuracy sampling at a lower
computational cost to study intercalation chemistries in energy storage
materials.

##  VI Acknowledgements

This work was funded by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences, Materials Sciences and Engineering Division
under Contract No. DE-AC0205CH11231 (Materials Project program KC23MP). The
work was also supported by the computational resources provided by the Extreme
Science and Engineering Discovery Environment (XSEDE), supported by National
Science Foundation grant number ACI1053575; the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility located at Lawrence Berkeley National Laboratory; and
the Lawrencium Computational Cluster resource provided by the IT Division at
the Lawrence Berkeley National Laboratory.

##  Appendix A DFT Calculations

DFT calculations were performed with the Vienna ab initio simulation package
(VASP) using the projector-augmented wave method [62, 63], a plane-wave basis
set with an energy cutoff of 520 eV, and a reciprocal space discretization of
25 k-points per Å-1. All the calculations were converged to $10^{-6}$ eV in
total energy for electronic loops and 0.02 eV/Å in interatomic forces for
ionic loops. To model the Li-Mn-Nb-O-F oxyfluoride system, we relied on the
regularized strongly constrained and appropriately normed meta-GGA exchange-
correlation functional (r2SCAN) [47, 64], which is believed to better capture
the cation–anion hybridization and Li coordination preference [65]. r2SCAN
provides better computational efficiency than the earlier version of SCAN
[66].

##  Appendix B Table Exchange in sGCMC

Because the TE Metropolis steps are proposed by randomly selecting and
exchanging species in a configuration, the forward and backward proposal
probability for a TE step can be imbalanced. To ensure detailed balance in the
simulation, the acceptance probability of a TE step is adjusted as follows:

$\displaystyle p_{\bm{\sigma\sigma^{\prime}}}$
$\displaystyle=\min\left\\{1,r\right\\},$ (B.1) $\displaystyle r$
$\displaystyle=\frac{\prod_{u_{s}\neq 0}n_{s}!}{\prod_{u_{s}\neq
0}(n_{s}+u_{s})!}\exp\left[-\frac{\left(\Delta
E_{\boldsymbol{\sigma\sigma^{\prime}}}-\sum_{s}\mu_{s}u_{s}\right)}{k_{B}T}\right],$

where $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma^{\prime}}$ are the
configurations before and after applying the TE step, respectively; $n_{s}$ is
the number of species $s$ in configuration $\bm{\sigma}$ before the step;
$u_{s}$ represents the change of species $s$ amount after the step; and the
vector value $\boldsymbol{u}=(u_{1},u_{2},...,u_{s})$ must fall into one of
the three aforementioned TE categories or their inverses. All TE categories
and their inverses were selected with equal probability.

Additionally, in the TE sGCMC, a percentage ($w=20\%$) of canonical swaps
between the Li/vacancy and different oxidation states of the same elements are
mixed with TE steps to improve the efficiency of exploring different Li and
electronic configurations within the same compositions. For each chemical
potential in sGCMC, we applied simulated annealing from T = 5000 K to 100 K to
approach the Li/vacancy and oxidation states ground-state configuration. Then,
starting from such a ground state, the sGCMC was run at T = 300 K with 500,000
steps to thermally equilibrate the system, and another 500,000 steps were used
to generate a sample of configurations for analysis [39].

##  Appendix C Computational Efficiency

At each Metropolis step $p$ in semigrand-canonical MC, the Li content
($x_{\text{Li}}$) at the current configuration is recorded as
$x_{\text{Li},p}$. We denote $\overline{x_{\text{Li},[p,q]}}$ as the average
of $x_{\text{Li}}$ in a block from simulation step $p$ to step $q$. After
thermal equilibration, we define the _block mean variance_ at block length
$L$, $\text{Var}(\overline{x_{\text{Li}}}_{L})$, as the variance of the _block
means_
$\overline{x_{\text{Li},[p,p+L]}},\overline{x_{\text{Li},[p+L,p+2L]}},...$ for
each block containing $L$ samples. The block mean variance can be used as a
measure of uncertainty when estimating $\overline{\theta}$ using a block mean.
With the block mean variance, we further define the _computational efficiency_
,

$\text{eff}_{t}(x_{\text{Li}})=\frac{\tau^{2}}{\overline{T}_{L}\text{Var}(\overline{x_{\text{Li}}}_{L})},$
(C.1)

where $\tau^{2}$ is the ensemble variance of Li content $x_{\text{Li}}$
approximated by the variance of all MC steps after thermal equilibration,
$\overline{T}_{L}$ is the average CPU time spent in each $L$-steps block.
Given the same set of hardware used in MC simulations, the higher the
computational efficiency, the less sampling time required for the uncertainty
of the average Li content (i.e., the block mean variance) to be decreased
below the same threshold. In brief, higher sampling efficiency means less
sampling time.

## References

  * Olivetti  _et al._ [2017] E. A. Olivetti, G. Ceder, G. G. Gaustad, and X. Fu, Lithium-Ion Battery Supply Chain Considerations: Analysis of Potential Bottlenecks in Critical Metals, Joule 1, 229 (2017).
  * Goodenough and Kim [2010] J. B. Goodenough and Y. Kim, Challenges for Rechargeable Li Batteries, Chemistry of Materials 22, 587 (2010).
  * Ji  _et al._ [2019a] H. Ji, A. Urban, D. A. Kitchaev, D.-H. Kwon, N. Artrith, C. Ophus, W. Huang, Z. Cai, T. Shi, J. C. Kim, H. Kim, and G. Ceder, Hidden structural and chemical order controls lithium transport in cation-disordered oxides for rechargeable batteries, Nature Communications 10, 592 (2019a).
  * Richards  _et al._ [2018] W. D. Richards, S. T. Dacek, D. A. Kitchaev, and G. Ceder, Fluorination of Lithium‐Excess Transition Metal Oxide Cathode Materials, Advanced Energy Materials 8, 1701533 (2018).
  * Huang  _et al._ [2022] J. Huang, P. Zhong, Y. Ha, Z. Lun, Y. Tian, M. Balasubramanian, W. Yang, and G. Ceder, Oxygen Vacancy Introduction to Increase the Capacity and Voltage Retention in Li‐Excess Cathode Materials, Small Structures , 2200343 (2022).
  * Zhang  _et al._ [2022] R. Zhang, C. Wang, P. Zou, R. Lin, L. Ma, L. Yin, T. Li, W. Xu, H. Jia, Q. Li, S. Sainio, K. Kisslinger, S. E. Trask, S. N. Ehrlich, Y. Yang, A. M. Kiss, M. Ge, B. J. Polzin, S. J. Lee, W. Xu, Y. Ren, and H. L. Xin, Compositionally complex doping for zero-strain zero-cobalt layered cathodes, Nature 610, 67 (2022).
  * Zhou  _et al._ [2022] K. Zhou, Y. Li, Y. Ha, M. Zhang, W. Dachraoui, H. Liu, C. Zhang, X. Liu, F. Liu, C. Battaglia, W. Yang, J. Liu, and Y. Yang, A Nearly Zero-Strain Li-Rich Rock-Salt Oxide with Multielectron Redox Reactions as a Cathode for Li-Ion Batteries, Chemistry of Materials 34, 9711 (2022).
  * Zhao  _et al._ [2022] X. Zhao, Y. Tian, Z. Lun, Z. Cai, T. Chen, B. Ouyang, and G. Ceder, Design principles for zero-strain Li-ion cathodes, Joule 6, 1654 (2022).
  * Meng and Arroyo-De Dompablo [2009] Y. S. Meng and M. E. Arroyo-De Dompablo, First principles computational materials design for energy storage materials in lithium ion batteries, Energy and Environmental Science 2, 589 (2009).
  * Urban  _et al._ [2016] A. Urban, D.-H. Seo, and G. Ceder, Computational understanding of Li-ion batteries, npj Computational Materials 2, 16002 (2016).
  * Van Der Ven  _et al._ [2020] A. Van Der Ven, Z. Deng, S. Banerjee, and S. P. Ong, Rechargeable Alkali-Ion Battery Materials: Theory and Computation, Chemical Reviews 120, 6977 (2020).
  * Aydinol  _et al._ [1997] 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, Physical Review B 56, 1354 (1997).
  * Urban  _et al._ [2017] A. Urban, A. Abdellahi, S. Dacek, N. Artrith, and G. Ceder, Electronic-Structure Origin of Cation Disorder in Transition-Metal Oxides, Physical Review Letters 119, 1 (2017).
  * Lun  _et al._ [2019] Z. Lun, B. Ouyang, D. A. Kitchaev, R. J. Clément, J. K. Papp, M. Balasubramanian, Y. Tian, T. Lei, T. Shi, B. D. McCloskey, J. Lee, and G. Ceder, Improved Cycling Performance of Li-Excess Cation-Disordered Cathode Materials upon Fluorine Substitution, Advanced Energy Materials 9, 1802959 (2019).
  * Zhou  _et al._ [2006] F. Zhou, T. Maxisch, and G. Ceder, Configurational Electronic Entropy and the Phase Diagram of Mixed-Valence Oxides: The Case of LixFePO4, Physical Review Letters 97, 155704 (2006).
  * Lee  _et al._ [2021] J. Lee, C. Wang, R. Malik, Y. Dong, Y. Huang, D. Seo, and J. Li, Determining the Criticality of Li‐Excess for Disordered‐Rocksalt Li‐Ion Battery Cathodes, Advanced Energy Materials 11, 2100204 (2021).
  * Li  _et al._ [2021] B. Li, M. T. Sougrati, G. Rousse, A. V. Morozov, R. Dedryvère, A. Iadecola, A. Senyshyn, L. Zhang, A. M. Abakumov, M.-L. Doublet, and J.-M. Tarascon, Correlating ligand-to-metal charge transfer with voltage hysteresis in a Li-rich rock-salt compound exhibiting anionic redox, Nature Chemistry 13, 1070 (2021).
  * Yao  _et al._ [2018] Z. Yao, S. Kim, J. He, V. I. Hegde, and C. Wolverton, Interplay of cation and anion redox in Li4Mn2O5 cathode material and prediction of improved Li4(Mn,M)2O5 electrodes for Li-ion batteries, Science Advances 4, eaao6754 (2018).
  * Li  _et al._ [2018] Q. Li, J. Wu, Z. Yao, Y. Xu, M. M. Thackeray, C. Wolverton, and V. P. Dravid, Dynamic imaging of metastable reaction pathways in lithiated cobalt oxide electrodes, Nano Energy 44, 15 (2018).
  * Ji  _et al._ [2019b] H. Ji, D. A. Kitchaev, Z. Lun, H. Kim, E. Foley, D. H. Kwon, Y. Tian, M. Balasubramanian, M. Bianchini, Z. Cai, R. J. Clément, J. C. Kim, and G. Ceder, Computational Investigation and Experimental Realization of Disordered High-Capacity Li-Ion Cathodes Based on Ni Redox, Chemistry of Materials 31, 2431 (2019b).
  * McColl  _et al._ [2022] K. McColl, R. A. House, G. J. Rees, A. G. Squires, S. W. Coles, P. G. Bruce, B. J. Morgan, and M. S. Islam, Transition metal migration and O2 formation underpin voltage hysteresis in oxygen-redox disordered rocksalt cathodes, Nature Communications 13, 5275 (2022).
  * Ångqvist  _et al._ [2019] M. Ångqvist, W. A. Muñoz, J. M. Rahm, E. Fransson, C. Durniak, P. Rozyczko, T. H. Rod, and P. Erhart, Icet–a python library for constructing and sampling alloy cluster expansions, Advanced Theory and Simulations 2, 1900015 (2019).
  * Wu  _et al._ [2016] Q. Wu, B. He, T. Song, J. Gao, and S. Shi, Cluster expansion method and its application in computational materials science, Computational Materials Science 125, 243 (2016).
  * van de Walle  _et al._ [2002] A. van de Walle, M. Asta, and G. Ceder, The alloy theoretic automated toolkit: A user guide, Calphad 26, 539 (2002).
  * Seko and Tanaka [2014] A. Seko and I. Tanaka, Cluster expansion of multicomponent ionic systems with controlled accuracy: importance of long-range interactions in heterovalent ionic systems, Journal of Physics: Condensed Matter 26, 115403 (2014).
  * Barroso-Luque  _et al._ [2022a] L. Barroso-Luque, P. Zhong, J. H. Yang, F. Xie, T. Chen, B. Ouyang, and G. Ceder, Cluster expansions of multicomponent ionic materials: Formalism and methodology, Physical Review B 106, 144202 (2022a).
  * van de Walle [2009] A. van de Walle, Multicomponent multisublattice alloys, nonconfigurational entropy and other additions to the Alloy Theoretic Automated Toolkit, Calphad 33, 266 (2009).
  * Nelson  _et al._ [2013] L. J. Nelson, G. L. W. Hart, F. Zhou, and V. Ozoliņš, Compressive sensing as a paradigm for building physics models, Physical Review B 87, 035125 (2013).
  * Seko  _et al._ [2014] A. Seko, A. Takahashi, and I. Tanaka, Sparse representation for a potential energy surface, Phys. Rev. B 90, 024101 (2014).
  * Zhong  _et al._ [2022] P. Zhong, T. Chen, L. Barroso-Luque, F. Xie, and G. Ceder, An $\ell_{0}\ell_{2}$-norm regularized regression model for construction of robust cluster expansions in multicomponent systems, Physical Review B 106, 024203 (2022).
  * Puchala  _et al._ [2023] B. Puchala, J. C. Thomas, A. R. Natarajan, J. G. Goiri, S. S. Behara, J. L. Kaufman, and A. Van der Ven, CASM — A software package for first-principles based study of multicomponent crystalline solids, Computational Materials Science 217, 111897 (2023).
  * Van der Ven  _et al._ [1998] A. Van der Ven, M. K. Aydinol, G. Ceder, G. Kresse, and J. Hafner, First-principles investigation of phase stability in LixCoO2, Physical Review B 58, 2975 (1998).
  * Wolverton and Zunger [1998] C. Wolverton and A. Zunger, First-principles prediction of vacancy order-disorder and intercalation battery voltages in LixCoO2, Physical Review Letters 81, 606 (1998).
  * Houchins and Viswanathan [2020] G. Houchins and V. Viswanathan, An accurate machine-learning calculator for optimization of Li-ion battery cathodes, The Journal of Chemical Physics 153, 054124 (2020).
  * Kolli and Van Der Ven [2018] S. K. Kolli and A. Van Der Ven, First-Principles Study of Spinel MgTiS2 as a Cathode Material, Chemistry of Materials 30, 2436 (2018).
  * Guo  _et al._ [2023] X. Guo, C. Chen, and S. P. Ong, Intercalation Chemistry of the Disordered Rocksalt Li3V2O5 Anode from Cluster Expansions and Machine Learning Interatomic Potentials, Chemistry of Materials 35, 1537 (2023).
  * Valleau and Cohen [1980] J. P. Valleau and L. K. Cohen, Primitive model electrolytes. I. Grand canonical Monte Carlo computations, The Journal of Chemical Physics 72, 5935 (1980).
  * Deng  _et al._ [2020] Z. Deng, G. Sai Gautam, S. K. Kolli, J. N. Chotard, A. K. Cheetham, C. Masquelier, and P. Canepa, Phase Behavior in Rhombohedral NaSiCON Electrolytes and Electrodes, Chemistry of Materials 32, 7908 (2020).
  * Xie  _et al._ [2023] F. Xie, P. Zhong, L. Barroso-Luque, B. Ouyang, and G. Ceder, Semigrand-canonical Monte-Carlo simulation methods for charge-decorated cluster expansions, Computational Materials Science 218, 112000 (2023).
  * Yang  _et al._ [2022] J. H. Yang, T. Chen, L. Barroso-Luque, Z. Jadidi, and G. Ceder, Approaches for handling high-dimensional cluster expansions of ionic systems, npj Computational Materials 8, 133 (2022).
  * Seo  _et al._ [2016] D.-H. Seo, J. Lee, A. Urban, R. Malik, S. Kang, and G. Ceder, The structural and chemical origin of the oxygen redox activity in layered and cation-disordered Li-excess cathode materials, Nature Chemistry 8, 692 (2016).
  * Barroso-Luque  _et al._ [2022b] L. Barroso-Luque, J. H. Yang, F. Xie, T. Chen, R. L. Kam, Z. Jadidi, P. Zhong, and G. Ceder, smol: A Python package for cluster expansions and beyond, Journal of Open Source Software 7, 4504 (2022b).
  * Abdellahi  _et al._ [2016] A. Abdellahi, A. Urban, S. Dacek, and G. Ceder, The Effect of Cation Disorder on the Average Li Intercalation Voltage of Transition-Metal Oxides, Chemistry of Materials 28, 3659 (2016).
  * Squires and Scanlon [2023] A. G. Squires and D. O. Scanlon, Understanding the limits to short-range order suppression in many-component disordered rock salt lithium-ion cathode materials, Journal of Materials Chemistry A 11, 13765 (2023).
  * Huang  _et al._ [2023a] L. Huang, P. Zhong, Y. Ha, Z. Cai, Y. Byeon, T. Huang, Y. Sun, F. Xie, H. Hau, H. Kim, M. Balasubramanian, B. D. McCloskey, W. Yang, and G. Ceder, Optimizing Li‐Excess Cation‐Disordered Rocksalt Cathode Design Through Partial Li Deficiency, Advanced Energy Materials 13, 2202345 (2023a).
  * Zhou  _et al._ [2004] F. Zhou, M. Cococcioni, C. A. Marianetti, D. Morgan, and G. Ceder, First-principles prediction of redox potentials in transition-metal compounds with LDA$+U$, Physical Review B 70, 235121 (2004), 0406382 .
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional, Physical Review Letters 115, 036402 (2015).
  * Reed and Ceder [2004] J. Reed and G. Ceder, Role of Electronic Structure in the Susceptibility of Metastable Transition-Metal Oxide Structures to Transformation, Chemical Reviews 104, 4513 (2004).
  * Huang  _et al._ [2023b] T. Y. Huang, Z. Cai, M. J. Crafton, L. A. Kaufman, Z. M. Konz, H. K. Bergstrom, E. A. Kedzie, H. M. Hao, G. Ceder, and B. D. McCloskey, Quantitative Decoupling of Oxygen-Redox and Manganese-Redox Voltage Hysteresis in a Cation-Disordered Rock Salt Cathode, Advanced Energy Materials 2300241, 1 (2023b).
  * Ceder [1993] G. Ceder, A derivation of the Ising model for the computation of phase diagrams, Computational Materials Science 1, 144 (1993).
  * Ménétrier  _et al._ [1999] 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, 1135 (1999).
  * Zunger  _et al._ [1990] A. Zunger, S.-H. Wei, L. G. Ferreira, and J. E. Bernard, Special quasirandom structures, Physical Review Letters 65, 353 (1990).
  * Huang  _et al._ [2021] J. Huang, P. Zhong, Y. Ha, D.-h. Kwon, M. J. Crafton, Y. Tian, M. Balasubramanian, B. D. McCloskey, W. Yang, and G. Ceder, Non-topotactic reactions enable high rate capability in Li-rich cathode materials, Nature Energy 6, 706 (2021).
  * Lun  _et al._ [2021] Z. Lun, B. Ouyang, D.-h. Kwon, Y. Ha, E. E. Foley, T.-Y. Huang, Z. Cai, H. Kim, M. Balasubramanian, Y. Sun, J. Huang, Y. Tian, H. Kim, B. D. McCloskey, W. Yang, R. J. Clément, H. Ji, and G. Ceder, Cation-disordered rocksalt-type high-entropy cathodes for Li-ion batteries, Nature Materials 20, 214 (2021).
  * Nakajima and Yabuuchi [2017] M. Nakajima and N. Yabuuchi, Lithium-Excess Cation-Disordered Rocksalt-Type Oxide with Nanoscale Phase Segregation: Li1.25Nb0.25V0.5O2, Chemistry of Materials 29, 6927 (2017).
  * Lebens-Higgins  _et al._ [2021] Z. Lebens-Higgins, H. Chung, I. Temprano, M. Zuba, J. Wu, J. Rana, C. Mejia, M. A. Jones, L. Wang, C. P. Grey, _et al._ , Electrochemical Utilization of Iron IV in the Li1.3Fe0.4Nb0.3O2 Disordered Rocksalt Cathode, Batteries & Supercaps 4, 771 (2021).
  * Weinan [2020] E. Weinan, Machine learning and computational mathematics, arXiv preprint arXiv:2009.14596 (2020).
  * Xie  _et al._ [2022] J. Z. Xie, X. Y. Zhou, D. Luan, and H. Jiang, Machine Learning Force Field Aided Cluster Expansion Approach to Configurationally Disordered Materials: Critical Assessment of Training Set Selection and Size Convergence, Journal of Chemical Theory and Computation 18, 3795 (2022).
  * Chen and Ong [2022] C. Chen and S. P. Ong, A universal graph deep learning interatomic potential for the periodic table, Nature Computational Science 2, 718 (2022).
  * Deng  _et al._ [2023] B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel, and G. Ceder, CHGNet: Pretrained universal neural network potential for charge-informed atomistic modeling, arXiv preprint arXiv:2302.14231 (2023).
  * Takamoto  _et al._ [2022] S. Takamoto, C. Shinagawa, D. Motoki, K. Nakago, W. Li, I. Kurata, T. Watanabe, Y. Yayama, H. Iriguchi, Y. Asano, T. Onodera, T. Ishii, T. Kudo, H. Ono, R. Sawada, R. Ishitani, M. Ong, T. Yamaguchi, T. Kataoka, A. Hayashi, N. Charoenphakdee, and T. Ibuka, Towards universal neural network potential for material discovery applicable to arbitrary combination of 45 elements, Nature Communications 13, 1 (2022).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Physical Review B 59, 1758 (1999).
  * Furness  _et al._ [2020] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation, The Journal of Physical Chemistry Letters 11, 8208 (2020), 2008.03374 .
  * Zhang  _et al._ [2018] Y. Zhang, D. A. Kitchaev, J. Yang, T. Chen, S. T. Dacek, R. A. Sarmiento-Pérez, M. A. Marques, H. Peng, G. Ceder, J. P. Perdew, _et al._ , Efficient first-principles prediction of solid stability: Towards chemical accuracy, npj Computational Materials 4, 1 (2018).
  * Kingsbury  _et al._ [2022] R. Kingsbury, A. S. Gupta, C. J. Bartel, J. M. Munro, S. Dwaraknath, M. Horton, and K. A. Persson, Performance comparison of r2SCAN and SCAN metaGGA density functionals for solid materials via an automated, high-throughput computational workflow, Physical Review Materials 6, 013801 (2022).