Teaching oxidation states to neural networks

Abstract

While the accurate description of redox reactions remains a challenge for first-principles calculations, it has been shown that extended Hubbard functionals (DFT+U+V) can provide a reliable approach, mitigating self-interaction errors, in materials with strongly localized d or f electrons. Here, we first show that DFT+U+V molecular dynamics is capable of following the adiabatic evolution of oxidation states over time, using representative Li-ion cathode materials. In turn, this allows to develop redox-aware machine-learning potentials. We show that considering atoms with different oxidation states (as accurately predicted by DFT+U+V) as distinct species in the training leads to potentials that are able to identify the correct ground state and pattern of oxidation states for redox elements present. This can be achieved, e.g., through a systematic combinatorial search for the lowest-energy configuration or with stochastic methods. This brings the advantages of machine-learning potentials to key technological applications (e.g., rechargeable batteries), which require an accurate description of the evolution of redox states.

Introduction

The concept of oxidation state (OS) is ubiquitous in physics and chemistry. The evolution of oxidation states on the atoms of a system follows redox reactions, electrolysis, and many other crucial electrochemical processes. These processes are essential for understanding numerous phenomena underlying contemporary technologies^{1, 2–3}. Redox reactions, for example, are at the heart of the functioning of rechargeable batteries, enabling the conversion and storage of energy^{4, 5, 6–7}. Modeling these processes is fundamental for understanding the underlying physical mechanism, allowing the identification and optimization of key factors influencing the behavior of materials, and laying the groundwork for the development of possible innovative solutions.

In this domain, density-functional theory (DFT)^8,9 plays a fundamental role, enabling first-principles modeling of atomic systems by drawing on the fundamental tenets of quantum mechanics. This approach has led to significant advances in the analysis, characterization, and discovery of materials. Moreover, the combination of DFT with molecular dynamics (i.e., first-principles molecular dynamics, FPMD) has further extended first-principles studies from static properties (at zero temperature) to dynamic properties (at finite temperature), providing a comprehensive and general investigative tool¹⁰.

However, within DFT, the accurate description of redox reactions remains a challenge. One difficulty stems from the fact that a formal definition of oxidation state is lacking. The generally accepted IUPAC definition states¹¹: “The oxidation state of an atom is the charge of this atom after the ionic approximation of its heteronuclear bonds.” While the charge is often used as a descriptor for the oxidation state in atomistic simulations (e.g., in Bader’s¹² and Voronoi’s¹³ methods), it has been shown that this can sometimes lead to misleading descriptions^14,15. The challenge of accurately and clearly defining oxidation states has driven the development of novel methods for calculating them from first principles. For instance, the oxidation state can be defined using topological considerations^16,17 based on the modern theory of polarization^18,19, or it can be determined from the electronic population derived from the eigenvalues of the atomic occupation matrix obtained by projecting the electronic wavefunctions onto specific electronic manifolds (e.g., d orbitals of transition-metal elements)²⁰. Whereas the OS as defined in ref. ¹⁶ has proven to be effective for transport processes²¹, the method of ref. ²⁰ is particularly well suited for the battery cathode materials, as shown in ref. ²².

Alongside this general aspect, it is important to note that standard DFT, which uses approximated exchange-correlation functionals, is affected by self-interaction errors (SIEs), leading to unphysical delocalizations of electrons. This limitation prevents accurate modeling of processes where changes in atomic oxidation states are crucial, particularly in systems with strongly localized d or f electrons. Such systems, including transition-metal oxides, are often of considerable interest for practical applications, such as in rechargeable battery technology^23,24. Various methods have been proposed to correct or alleviate SIEs: DFT+U^{25, 26, 27–28} along with the DFT+U+V extension^{29, 30–31}, meta-GGA functionals, such as SCAN and other derivatives^{32, 33–34} (as well as SCAN+U^{35, 36, 37–38}), or hybrid functionals (e.g., PBE0³⁹ and HSE06^40,41), to name a few. In this work, we focus on DFT+U+V, which we will analyze in detail throughout the article, due to its accuracy and affordable computational cost. In this regard, we note how it has recently been shown to improve the prediction of electrochemical and thermodynamic properties when localization is accompanied by significant hybridization, as in ref. ⁴². Following the methodology suggested in ref. ²⁰ to compute oxidation states, DFT+U+V has been shown to deliver sharp ("digital”) changes in oxidation states for transition metals within intercalation cathode materials—such as Mn and Fe phospho-olivines (LMPO, LFPO)—as lithium concentration varies, accurately mirroring the shifts in oxidation states of the transition-metal atoms²².

Thus, applying DFT+U+V to phospho-olivine cathodes provides an ideal starting point to investigate the landscape of oxidation states into material and its evolution. In this study, we specifically explore finite-temperature effects in LMPO using DFT+U+V FPMD, where, during the dynamics, we observe a clear evolution of OSs and can track the adiabatic ground state across the redox-active elements. While electron-transfer reactions are most often non-adiabatic processes⁴³ and can, for example, be described by the general framework of the Marcus theory^{44, 45–46}, it is essential to be able to capture faithfully first the adiabatic ground state driven by the new atomic arrangements resulting from the FPMD evolution.

FPMD simulations with DFT+U+V offer significant advantages, such as reduced computational cost and enhanced accuracy, compared to hybrid functionals, especially for transition-metal oxides²². A key reason is that the screened exchange interaction, captured by the Hubbard parameters, is inherently material-specific in the linear-response formulation of DFT+U+V, while hybrid functionals, in addition to cost and complexity, adopt a general recipe that is not always optimal depending on the material or system of interest.

We showcase here the significant potential of DFT+U+V FPMD; however, the cost of first-principles simulations remains a limiting factor, potentially restricting studies to length and timescales that may be insufficient to observe new phenomena such as ion migration, phase transformations, and chemical reactions. Nowadays, it is well established that machine-learning interatomic potentials can help with these challenges. These potentials use a limited set of first-principles calculations for training and are capable of delivering predictions with first-principles accuracy at a cost comparable to that of classical force field simulations⁴⁷. A wide variety of methods exist for building machine learning interatomic potentials; recently, graph neural networks that incorporate equivariant symmetry constraints have been shown to achieve state-of-the-art performance, and notable examples include NequIP⁴⁸ and MACE⁴⁹. The equivariant nature of these models allows them to outperform some invariant models (e.g., DeepMD⁵⁰) using a fraction of the data for training. However, despite their success, there is still a limiting factor, i.e., lack of control over the atomic OSs. The importance of the atomic OSs stems from the fact that ions with different OS behave differently. This is generally true, but let us consider the example of Mn we will also encounter later: high spin Mn⁴⁺ is a non-bonding spherical ion that almost always adopts octahedral coordination with oxygen atoms, Mn³⁺ is a Jahn-Teller active ion that radically distorts its environment, while Mn²⁺ preferably exists in tetrahedral coordination⁵¹.

Several schemes have been proposed to describe systems with multiple charge states at the same geometry. Models based on charge equilibration approaches^{52, 53, 54, 55, 56–57} are expected to incorporate information about oxidation states; however, whether and to what extent they actually capture this information and effectively describe redox chemistry remains an open question, despite recent rationalizations that have offered some insight⁵⁸. Most of these models adopt a dual-learning approach, requiring additional neural networks to learn atomic charges alongside energies and forces. Other models, such as that in ref. ⁵⁹, function as universal charge-informed potentials but depend on additional training variables, such as magnetic moments. Furthermore, some models incorporate the concept of oxidation state through geometrical considerations⁶⁰. While the atomic environment and geometry may at times be sufficient to determine the OS of an atom, this is not universally the case. In liquid environments, for example, the mobility of solvent molecules and other chemical species prevents the establishment of a fixed, well-defined atomic environment that could reliably serve as a label for defining the oxidation state. Therefore, a simple model oxidation state informed that could leverage recent state-of-the-art tools (e.g., NequIP⁴⁸ or MACE⁴⁹), requiring an uncomplicated training procedure without additional dual-learning schemes or training variables, and capable of going beyond information deduced solely from geometry, would represent a valuable option in many cases of study.

In this work, building on the fact that different oxidation states of transition-metal ions behave differently from each other as they were different elements, and emphasizing the accuracy over oxidation states achieved through DFT+U+V, we propose an equivariant neural network potential that treats atoms of the same element with different oxidation states as different species. Once the potential has been generated, we demonstrate not only its accuracy but also that the correct arrangement of oxidation states can be determined through a combinatorial search for the lowest-energy configuration, reproducing correctly the adiabatic rearrangement of OSs observed in the FPMD.

Results

Static DFT+U+V: oxidation states and voltages

DFT+U+V is briefly outlined in the “Methods” section. Here, we discuss the main results where we investigate the phospho-olivine cathode material Li_xMnPO₄, comparing standard DFT and DFT+U+V. The crystal structure is orthorhombic at x = 0 and x = 1 with a Pnma space group. In the present simulations, the unit cell contains four formula units, i.e., 24 atoms for x = 0 and 28 atoms for x = 1. The four Mn atoms are each coordinated by six oxygen atoms, forming a MnO₆ octahedra with the Mn atom centrally positioned. In the DFT+U+V framework, the U correction is applied to the 3d orbitals of Mn, while the V is set between the 3d orbitals of Mn and the 2p orbitals of the surrounding O atoms. The phospho-olivines are known to be antiferromagnetic: we use the magnetic configuration that minimizes the total energy (labeled AF₁ in ref. ⁴²). We consider all possible concentrations of Li x = 0, 1/4, 1/2, 3/4, 1.

The OS of Mn atoms within the material can be determined by analyzing the electronic population of their 3d shells. Following the approach of ref. ²⁰, occupation numbers of the 3d shells are derived from the eigenvalues of the site-diagonal atomic occupation matrix (i.e., I = J in Eq. (3)), which has a 5 × 5 size in the respective spin-up and spin-down channels. Hence, this approach provides the 10 electronic occupation numbers for the 3d shells of Mn. Next, we count how many d states are “fully occupied” (i.e., states with corresponding occupation number approximately 1) and compare to the valence electronic configuration of the Mn atom, which has 7 valence electrons: 3d⁵4s². For example, in the fully de-lithiated case with x = 0, the DFT+U+V calculation provides the following occupations for all four Mn atoms: $n_i^{\uparrow }=0.50$, 0.99,0.99,1.00,1.00 for the spin-up channel, and $n_i^{\downarrow }=0.05$, 0.06, 0.08, 0.09, 0.22 for the spin-down channel. In this case, occupations equal or smaller than ~0.5 indicate an original unoccupied 3d state of Mn and, ~1 a fully occupied state. The fractional occupations are a result of the mixing between Mn-3d and O-2p orbitals; even though the original Mn-3d orbital is empty, contribution to the occupations might arise from originally occupied O-2p orbitals that contribute to the projections. Furthermore, greater deviations from 0 in the eigenvalues signal a stronger mixing, as illustrated by the first spin-up eigenvalue. Thus, the fully occupied states (indicated in bold) are 4; this indicates that 4 electrons are assigned to each Mn atom in the compound. Since Mn has 7 valence electrons, this implies that 7 − 4 = 3 electrons are involved in bonding with the oxygen environment, yielding an oxidation state of 3+. Now, let us consider the fully lithiated case with x = 1, where the system is fully lithiated (with 4 additional Li⁺ ions and 4 extra electrons). For each Mn atom, our DFT+U+V calculation provides the following occupations: $n_i^{\uparrow }=\mathbf{0.99}$, 0.99,1.00,1.00,1.00 and $n_i^{\downarrow }=0.02$, 0.02, 0.03, 0.07, 0.08. The fully occupied states (in bold) indicate that each Mn atom now accommodates 5 electrons, resulting in an oxidation state of 2+: each Mn atom has gained an electron and has been reduced.

We now move on to examining the system’s behavior at intermediate Li concentrations. For simplicity, we can use as a descriptor of oxidation state the sum of the electronic occupations $n={\sum }_i(n_i^{\uparrow }+n_i^{\downarrow })$. This quantity, often referred to as Löwdin occupation, is particularly useful for bookkeeping¹⁴ and especially for describing the (de-)lithiation process, as will be discussed below. For example, applying this definition to the DFT+U+V 3d-shells occupations of Mn presented above, we find n = 4.98 and n = 5.21 for the fully de-lithiated and fully lithiated systems, respectively. In Fig. 1a–c, we present the Löwdin occupations of the 3d shells of Mn atoms in Li_xMnPO₄ for intermediate Li concentrations (i.e., at x = 0, 1/4, 1/2, 3/4, 1), computed using three different approaches: (1) standard DFT, (2) DFT+U+V with self-consistent U and V, and (3) DFT+U+V with U and V parameters averaged over all Mn sites and all possible Li concentrations. The four Mn atoms within the unit cell are represented by bars in different shades. In going from x = 0 to x = 1/4, one Li atom is added to the system, which introduces one Li⁺ cation and one electron. Then, additional Li atoms are incrementally introduced until x = 1 is reached with four Li⁺ cations and four additional electrons in the system. Figure 1a shows that standard DFT delocalizes every electron added across all Mn sites, resulting in no significant differences in the Löwdin occupations across sites as the Li concentration changes; no individual Mn ion undergoes a redox reaction. Figure 1b presents the same study using DFT+U+V, and calculating self-consistent U and V parameters. To distinguish more clearly the oxidation states of Mn, we indicate with a blue dashed horizontal line the occupation level corresponding to the concentration x = 0, where all Mn atoms are 3+ (n = 4.98), and with a red dashed line the corresponding occupation level for the concentration x = 1, where all Mn atoms are 2+ (n = 5.21). DFT+U+V shows a clear “digital” change in the Löwdin occupations: during the lithiation process, the addition of one Li⁺ ion and one electron to the cathode alters the occupation of a single Mn ion from 4.98 to 5.21, corresponding to a change in oxidation state from 3+ to 2+, while all other Mn ions remain unaffected. This process continues with further Li intercalations, ultimately resulting in the reduction of all Mn ions from 3+ to 2+. Consequently, DFT+U+V with self-consistent parameters accurately captures the mixed-valence character of the Li_xMnPO₄ compound, which includes two distinct Mn ion states, Mn³⁺ and Mn²⁺, at x = 1/4, 1/2, 3/4. In contrast, standard DFT fails to localize the additional electrons on a specific Mn ion; instead, the charge density is delocalized, spreading almost equally across all Mn ions, resulting in approximately equal occupations, as shown in Fig. 1a. Thus, in DFT, at x = 1/4, 1/2, 3/4, there is effectively a single type of Mn ion with intermediate occupation values that progressively shift with Li content, corresponding to “unchemical” oxidation states of Mn^2.25+, Mn^2.5+, Mn^2.75+.

[See PDF for image]

Fig. 1

Löwdin occupations and voltage predictions for Li_xMnPO₄.

a–c Löwdin occupations n of the Mn 3d shells in Li_xMnPO₄ at Li concentrations x = 0, 1/4, 1/2, 3/4, and 1, computed using: a standard DFT, b DFT+U+V with self-consistent Hubbard parameters, and c DFT+U+V with averaged parameters (U = 5.1 eV, V = 0.7 eV). Each bar corresponds to one of the four Mn atoms in the unit cell and is distinguished by a different shade of red. In (b) and (c), the red dashed horizontal line indicates the Löwdin occupation corresponding to Mn²⁺ (n = 5.21), and the blue dashed line to Mn³⁺ (n = 4.98). d Voltages (in V) for Li_xMnPO₄ computed using DFT+U+V with averaged self-consistent U and V, compared with standard DFT, HSE06, DFT+U, and DFT+U+V with self-consistent U and V²². Experimental data are from refs. ^62,63.

We now want to explore the use of the DFT+U+V at finite temperatures by means of FPMD. During MD simulations, the positions of the atoms change, which can lead to variations in the Hubbard parameters. Ideally, U and V should be re-evaluated from first principles for each MD frame (self-consistently at fixed geometry). However, it is currently impractical to add the computationally intensive calculation of Hubbard parameters to each frame of an already demanding FPMD simulation. To address this challenge, machine learning approaches employed to predict the Hubbard parameters will greatly help⁶¹. In any case, we expect the use of average self-consistent U and V parameters to be enough in the majority of cases, with the added advantage that U and V do not change across configurations. We now verify this for our system. First, we examine whether the use of average U and V parameters, kept fixed during the MD runs, can still describe accurately the “digital” change of OSs of Mn atoms shown in Fig. 1b. These U and V parameters are calculated as the average of the self-consistent values (see Supplementary Table 1) across all Li concentrations, resulting in U = 5.1 eV and V = 0.7 eV. The results are shown in Fig. 1c. In order to quantify the effect of using average self-consistent U and V as opposed to the fully self-consistent ones, we display the same Löwdin occupation treshold for Mn²⁺ (n = 4.98, red dashed line) and Mn³⁺ (n = 5.21, blue dashed line) as obtained with the self-consistent parameters shown in Fig. 1b. Clearly, even when using average parameters the digital variations in the OSs of Mn atoms are preserved. The only effect is a slight reduction in the contrast between the 2+ and 3+ levels. As a practical note, we emphasize that the cells at different concentrations were generated from the fully lithiated cell and relaxed independently in parallel, consistently with the considered level of theory (i.e., DFT, self-consistent DFT+U+V, or averaged DFT+U+V).

Finally, in Fig. 1d, we report the voltages Φ in comparison with different methods²² and experiments^62,63. The voltages are computed following the relation⁴² −eΦ = E_x=1 − E_x=0 − E_Li, where −e is the electron charge, E_x=1 and E_x=0 are the total energies per formula unit of Li_xMnPO₄ at x = 1 and x = 0, respectively. E_Li is the total energy of bulk Li (representing the anode). The voltage computed with average self-consistent U and V is 3.96 V, slightly underestimating the one computed with self-consistent U and V, 4.20 V. Nevertheless, it is still closer to the experiment than standard DFT, hybrid functionals (HSE06), and DFT+U. Note also that we do not consider here meta-GGA functionals (SCAN or SCAN+U) because they do not seem to achieve the same level of accuracy of DFT+U+V or hybrids functionals, in particular when dealing with systems with strong hybridization. The reason is that they describe the screened exchange less accurately, resulting in a poorer description of the energetics⁶⁴.

In summary, the study of average U and V parameters across oxidation states and voltages supports their applicability as fixed parameters throughout the DFT+U+V FPMD simulations, with details and further investigations presented in the following section.

DFT+U+V molecular dynamics

We performed FPMD on the Li_xMnPO₄ system, where energies and forces are calculated using DFT+U+V (Eq. (1) and Eq. (4)) with average U and V parameters (U = 5.1 eV and V = 0.7 eV). Simulation were carried out for all possible Li concentrations: x = 0, 1/4, 1/3, 3/4, 1. The simulations were performed with a timestep δt = 4 fs within the NVT ensemble by using the stochastic-velocity rescaling algorithm⁶⁵. The change in oxidation states is mostly associated with local volume variations that are well described by NVT sampling.

For each Li concentration, we performed two simulations at different temperatures: T = 400 K and T = 900 K. In passing, we note that a higher temperature is necessary to allow Mn atoms to explore all possible configurations of their OSs (i.e., all possible combinations) in the timescale of the simulations.

In Fig. 2, we show an example of the evolution of the Löwdin occupations over time for the four Mn atoms, for the system with Li concentration x = 1/2, consisting of two Li⁺ ions and two additional electrons that can localize on any Mn atom depending on the actual atomic configuration. To facilitate the analysis, we include the occupation levels previously discussed in the static case to distinguish the two OSs of Mn: the 3+ levels (n = 4.98) indicated by the blue line and the 2+ levels (n = 5.21) indicated by the red line. In the simulation at T = 400 K (Fig. 2a), the occupations fluctuate within a small range of values but remain close to the initial level associated with the specific OS. Specifically, Mn1 and Mn4 remain in the 3+ state, while Mn2 and Mn3 remain in the 2+ state. The situation changes when we consider the temperature T = 900 K (Fig. 2b). The evolution of the occupations reveals distinct jumps between the two OSs. Electrons are exchanged adiabatically between Mn atoms, so there are always two Mn²⁺ and two Mn³⁺. For example, we see that approximately in the first 2 picoseconds, Mn1 is 3+, and correspondingly, Mn3 is 2+. Subsequently, we observe a sharp shift in OSs: Mn1 switches to the 2+ state, while simultaneously, Mn3 shifts to the 3+ state, indicating that the electron has moved adiabatically from Mn3 to Mn1. As can be observed, such transitions in OSs continue over time and also involve the other two Mn atoms. The adiabatic transitions occur within a time interval spanning from a few frames (about 2 or 3) to several frames (about 10). The oxygen environment readjusts, and a larger average volume of the oxygen cage is associated with the 2+ state. In reality, electron transfer will be driven by tunneling and happen almost instantaneously on these timescales, but not necessarily at the same time as the adiabatic crossing. We reiterate that this is not intended to reproduce the physical dynamics of the polaron hopping and tunneling but just the exploration of the electronic ground states. In the Supplementary Material, we further discuss the adiabatic transitions, also examining the individual electronic occupations.

[See PDF for image]

Fig. 2

Dynamics of oxidation states in DFT+U+V first-principles molecular dynamics of Li_xMnPO₄.

a–c Time evolution of Löwdin occupations for the four Mn atoms at x = 1/2. The horizontal red and blue dashed lines indicate reference values for Mn²⁺ (n = 5.21) and Mn³⁺ (n = 4.98), respectively. Shaded areas represent thermal fluctuation ranges. aT = 400 K, bT = 900 K, c filtered T = 900 K trajectory excluding intermediate n values. d Atomic configurations at t = A (top) and t = B (bottom) showing oxidation state changes. Mn²⁺ atoms are shown in red, Mn³⁺ in blue; Li in green, O in pink, and P in orange.

The finite temperature causes small fluctuations in the occupations, during which the atoms reasonably keep their OSs. However, when the new atomic configuration in the dynamics favors a different OS pattern, sharp transitions occur, which are distinct and well separated from the thermal fluctuations. In Fig. 2a–c, the extent of the fluctuations is represented by light red and light blue bands, within which Mn²⁺ and Mn³⁺ atoms preserve their own OS.

In Fig. 2b, c, two specific configurations A and B are highlighted (with vertical violet dashed lines) where there is a change in the distribution of OSs. In A, the four Mn atoms—Mn1, Mn2, Mn3, and Mn4—are 2+, 3+, 3+, and 2+, respectively. In B, they are 3+, 2+, 2+, and 3+, respectively. As an example, in Fig. 2d, the corresponding atomic configurations are shown. Mn²⁺ atoms are red, Mn³⁺ blue, and they are surrounded by other LMPO atoms, depicted in different colors. The upper cell shows configuration A, while the lower one shows configuration B, visually highlighting how, in these particular examples, the Mn atoms make a transition with exactly opposite OSs.

In the following section, we will explain how we selected data from the FPMD simulations to train the machine learning potential. Then, we test the potential and investigate its predictive power in determining the atomic OSs of Mn atoms on specific configurations.

Neural network potential and oxidation state identification

The FPMD trajectories are utilized to train an equivariant NN potential. For this purpose, we employ NequIP, which has demonstrated state-of-the-art performance⁴⁸. The training and validation datasets are obtained following the procedure described below (a schematic representation of the workflow is also provided in Supplementary Fig. 1). DFT+U+V FPMD simulations are conducted for each x_Li concentration of Li at temperatures T = 400 K and T = 900 K. Each trajectory spans approximately 9 ps, yielding a total combined trajectory length of around 92 ps. The trajectories are filtered as shown in Fig. 2c, ensuring that all Mn atoms have occupations within the red and blue bands —that is, excluding configurations where at least one Mn atom has Löwdin occupations n 5.05 < n < 5.15 out of the average fluctuations associated with each OS. The red and blue stripes have been chosen to encompass the amplitude of fluctuations observed in the T = 400 K simulations, where no jumps occur, and it is clear which fluctuation range corresponds to each oxidation state. Physically, the fluctuation in the occupations arises from the motion of Mn atoms and their surrounding environment. During these movements (in the absence of transitions), the Mn maintains its OS, and we heuristically chose these thresholds. After the filtering procedure, we obtain training segments summing up to 80 ps. Then, we proceed to identify all possible OSs patterns; those depend on the Li-ion concentration, which determines the number of additional electrons and, in turn, the number of Mn²⁺ atoms. Furthermore, for a fixed Li-ion concentration (i.e., at fixed number of Mn²⁺ and Mn³⁺), the OSs patterns of the Mn atoms can vary, as the different atomic configurations in the molecular dynamics simulations may promote the redistribution of localized electrons across Mn sites. Let us consider the example of x_Li = 1/4 with one Li⁺ cation and one additional electron in the system (with respect to the fully de-lithiated structure where all Mn atoms are 3+). The electron localizes on one of the Mn atoms, shifting its OS from Mn³⁺ to Mn²⁺, while the others remain Mn³⁺. Since any of the four Mn³⁺ atoms can receive the electron and become Mn²⁺, there are four possible OSs patterns. For simplicity, if we represent Mn²⁺ as 0 and Mn³⁺ as 1, the four patterns correspond to the permutations of 4 elements with 3 repeated (i.e., 4!/3! = 4): 0111, 1011, 1101, and 1110, where each number represents, in order, the OS of Mn1, Mn2, Mn3, and Mn4. Similarly, the concentration x_Li = 3/4 includes 4 possible OSs patterns, having now three Mn²⁺ and one Mn³⁺: 1000, 0100, 0010, 0001. For the concentration x_Li = 1/2, there are (4!/2!/2! = 6) 6 possible OSs patterns: 0110, 0011, 0101, 1001, 1010, 1100. Lastly, at x_Li = 0, all Mn atoms are in the 3+ state (1111), while at x_Li = 1, all Mn atoms are in the 2+ state (0000): each contributing with a single OS pattern. Thus, by adding up all possible OSs patterns for all Li-ion concentrations, we get a total number of 16 possible OSs patterns. In our simulations, we find that the FPMD at T = 400 K spans 9 possible of them across all concentrations, while the FPMD at T = 900 K includes the whole set of 16 patterns, although some correspond to a very small number of configurations. For example, the patterns 1010, 0011, and 0101 appear in only 31, 44, and 13 snapshots, respectively.

Then, each atomic configuration is assigned to a particular pattern of OSs, and this information is used to build the ensemble of data for the training and validation of the neural network potential. We proceed as follows (see also Supplementary Fig. 1): for each temperature and OS pattern, we select up to a maximum of 100 snapshots (if available), randomly chosen in order to minimize correlations as much as possible in the available trajectories. As a result of this procedure, we automatically include all possible Li concentrations. This selection ends up with a total number of 2288 snapshots that are divided into a training (80 %) and a validation set (20 %), to monitor training over time.

As mentioned, we train an equivariant neural network potential by using NequIP⁴⁸ and using atomic positions and DFT+U+V forces; Mn in different OSs are defined as different atomic types. Then, for testing, we select a maximum number of 50 random frames (if available) for each OS pattern and temperature, not included in the training and validation set. This results in an ensemble of approximately 1000 frames with which we perform error analysis and evaluate the performance of the model. The mean absolute error (MAE) on energies is 16.6 meV, the MAE on energies per atom is 0.64 meV, and the MAE on forces is 41 meV/Å (calculated on the mentioned test data). Then, the same testing dataset is used for parity plots in Fig. 3, where we compare the machine learning predictions of energies and forces vs. the DFT+U+V corresponding results reported in Fig. 3. The total energies are naturally clustered, each cluster associated with a different Li concentration. Both studies on energies and forces show excellent linear correlation (linear correlation coefficient r ~ 0.995 in all cases).

[See PDF for image]

Fig. 3

Parity plots of ML predictions vs. DFT+U+V for Li_xMnPO₄.

a–e Machine learning predictions of total energies vs. DFT+U+V total energies across all considered Li concentrations. f Machine learning predictions of forces vs. DFT+U+V forces (Cartesian components in color). All predictions show high agreement with DFT+U+V (linear correlation coefficient ~ 0.995), demonstrating excellent accuracy of the trained NequIP model.

As described above, during the FPMD simulations, Mn atoms change their OSs in response to the rearrangement of the atomic positions in the system. We now conduct an in-depth analysis of the ability of the neural network potential to predict such OSs transitions. The study is presented in Fig. 4 and described in the following. First, only configurations outside the training and validation datasets are chosen. In Fig. 4a, we consider a segment of the FPMD trajectory of Li_xMnPO₄ with x = 1/4 at T = 900 K. At this concentration, we have a single Li⁺ ion and one additional electron that localizes on one Mn atom. Then, we consider two atomic configurations, A and B, which are separated by an adiabatic “jump” of the electron from one Mn to another, indicated in the figure by dashed black vertical lines. For configuration A, the electron is on the Mn2 atom, which is 2+, while all the other Mn atoms are 3+. For configuration B, the electron has switched to the Mn1 atom, which now is the one with OS 2+. We refer to this rearrangement of OSs as “single jump,” meaning that only one couple of Mn atoms and one electron are involved, and we denote it as Mn2 → e⁻ → Mn1. The OSs patterns of Mn atoms in configuration A is 1011, and in configuration B is 0111. On the right side of Fig. 4a, we report block diagrams comparing DFT+U+V total energies of these configurations A and B (green) with ML predictions (orange) for all possible 4 OSs patterns. The study reveals that the DFT+U+V FPMD configuration is accurately identified as the one with the lowest energy by the machine-learning combinatorial search. This occurs for both the initial and final states, A and B. Thus, remarkably, the methodology precisely predicts the correct pattern of OSs.

[See PDF for image]

Fig. 4

Oxidation state transitions in Li_xMnPO₄ identified with the ML combinatorial search.

a Single-jump oxidation state transition at x = 1/4 and T = 900 K: one electron is transferred from Mn2 (initially Mn²⁺) to Mn1 (final Mn²⁺), changing the oxidation state pattern from 1011 to 0111. Here, 0 indicates Mn²⁺ and 1 indicates Mn³⁺. b Similar single-jump event at x = 3/4. c Double-jump transition at x = 1/2: Mn3 → e⁻ → Mn1 and Mn2 → e⁻ → Mn4, modifying the oxidation state pattern accordingly. d Independent DFT+U+V FPMD trajectory at x = 1/2 and T = 800 K confirms the same double-jump behavior. In all cases, the ML model identifies the DFT+U+V oxidation state pattern as the lowest-energy configuration, demonstrating its ability to accurately recover oxidation state rearrangements across different concentrations.

In Fig. 4b, we select instead a FPMD trajectory at 900 K with a Li concentration 3/4, meaning there are three additional electrons in the system and, as before, still 4 OSs patterns. The selected configurations A and B still lead to a “single jump” process where the electron moves according to Mn2 → e⁻ → Mn1. Evaluating the various OSs patterns with the neural network potential, once again, the DFT+U+V one is identified as the one having the lowest energy among all possible combinations.

In Fig. 4c, a trajectory at a concentration 1/2 is considered, where there are two additional electrons and 6 possible OSs patterns. The process is now a “double jump,” for which, between configurations A and B, electrons move as follows: Mn3 → e⁻ → Mn1 and Mn2 → e⁻ → Mn4. Once again, the machine-learning potential identifies the correct OS pattern for both the initial state A and the final state B, accurately reproducing the rearrangement of OSs.

Last, we performed a new, independent DFT+U+V FPMD simulation at a Li concentration of x = 1/2 and a different temperature of T = 800 K. In Fig. 4d, we show a portion of this new trajectory, which is completely independent of the previous datasets. We observe a “dual jump” process similar to that reported in Fig. 4c. Remarkably, the exploration and evaluation of various OSs patterns by the machine-learning potential once again proves highly effective in accurately describing the rearrangement of OSs.

In Fig. 5, we present the same investigation for an even larger supercell. The supercell is built by doubling the fully lithiated unit cell along the z direction and randomly excluding 5 Li atoms (among the available 8) to obtain the unseen concentration x = 3/8. The supercell thus includes 51 atoms. A portion of the FPMD trajectory at 800 K is presented in Fig. 5a, where we report the time evolution of the Löwdin occupations for the eight Mn atoms present in the supercell. Between the configurations A and B (indicated with vertical dashed lines in the figure), four Mn atoms shift their OS, namely Mn1, Mn3, Mn4 and Mn5. In this case, a larger number of OSs patterns are observed, specifically 8!/5!/3! = 56. The energy of all of them is evaluated with the neural network potential for both A and B configurations, and the results are shown in Fig. 5b, c. The pattern of OSs given by DFT+U+V is successfully predicted by the machine learning as the minimum energy configuration for both A and B states, reproducing the OSs transitions observed in the dynamics. Of course, in much larger supercells, it would not be intended to perform exponentially exploding combinatorial searches, substituting instead local explorations or Monte Carlo moves. The reader can find more transferability tests in the Supplementary Material and a comparison with CHGNet⁵⁹.

[See PDF for image]

Fig. 5

Oxidation state transitions in a larger supercell of Li_xMnPO₄ at x = 3/8.

a Time evolution of Löwdin occupations for eight Mn atoms in a 51-atom supercell at T = 800 K. Configurations A and B are marked by vertical violet dashed lines. Mn1, Mn3, Mn4, and Mn5 undergo oxidation state transitions. b Machine learning predicted energies for all 56 possible oxidation state patterns at configuration A. c Same as (b) for configuration B. In both cases, the DFT+U+V oxidation state pattern is identified by the ML model as the lowest-energy configuration.

We conclude this section with a comment on the application of the methodology to more complex systems. As proof of principle, we used the Mn phospho-olivine, LMPO. However, we expect this approach to be applicable in more complex cases as well, involving more challenging or multiple TM elements. For instance, it has been shown that within the mixed phospho-olivine LFMPO, both Fe and Mn exhibit in DFT+U+V a “digital” variation in the oxidation states²², so we expect the machine-learning potential to faithfully capture this. More complex cases may involve systems where the same OS is associated with different local environments, or even liquid systems that lack a well-defined local environment. In these cases, the occupation matrix method, not relying on particular geometric information, is still expected to distinguish between the different oxidation states. As long as representative configurations of interest (e.g., configurations with atoms in the same OS but different environments) are included in the training data, we expect the methodology to remain valid. Furthermore, as a note on the filtering procedure (for solids but also for liquid systems), it is always recommended to first analyze electronic occupations during the dynamics, preferably starting at moderate temperatures (in order to keep the atomic environments as stable as possible) and small systems (in order to have a complete control over all possible OS patterns). In general, the assignment of the appropriate filters results naturally from this analysis and from the fluctuations of the electronic occupations at fixed OS where no jumps are involved (e.g., Fig. 2a). Finally, it is always worth verifying the consistency of the OSs distribution after filtering: e.g., are electrons conserved? In a system like LMPO, inconsistencies might arise when only one occupation threshold is used to discriminate OSs. This could include transition states and lead to ambiguous assignments of OSs, resulting in an incorrect electron count for certain configurations. For example, in our system at x = 0.5 (with 2 Li cations, 2 additional electrons and 4 Mn), exactly 2 Mn²⁺ should be present.

Possible challenges are expected not in teaching oxidation states to the machine learning potential, but rather in having an electronic structure method that is accurate: for systems with low-spin states that might exhibit a multi-reference character, even the present advanced functionals might still be not sufficient. However, DFT+U+V is promising in its accuracy for a large class of materials.

Discussion

In the following, we schematically outline the main findings of this work:

• The DFT+U+V molecular dynamics approach has proven effective in exploring the ground-state landscape of oxidation states for redox-active elements; here, with the specific example of Mn in the LMPO cathode system. To achieve this, we used fixed U and V parameters derived from the ab initio self-consistent ones, averaging over all Mn atoms and Li concentrations. This was justified by the fact that average U and V still provide electrochemical properties—e.g., voltages—in better agreement with experiments than hybrid functionals, as well as a sharp description of oxidation state transitions. In passing, we add that DFT+U+V molecular dynamics is much less computationally intensive than hybrid functionals molecular dynamics. So, it is both more convenient and more accurate. If on-the-fly re-optimization of U and V parameters where necessary during the dynamics, machine learning models predicting Hubbard parameters⁶¹ may be of great help. For all these reasons, we believe DFT+U+V FPMD represents the state-of-the-art in descriptions of systems where localization and hybridization of d and f electrons play a central role.

• We showed how a machine learning potential can be trained to identify the correct oxidation states and switch those whenever necessary. The training is done on atomic positions and forces, but treating redox-active elements with different oxidation states as distinct types, where the oxidation states are determined by means of DFT+U+V. Notably, this procedure is expected to be agnostic to the method used for determining oxidation states (DFT+U+V) as well as to the specific machine learning machinery employed (NequIP). Regarding the potential’s training, the inclusion of various Li concentrations plays a crucial role in capturing the range of oxidation state patterns for subsequent predictions. However, it is reasonable to assume that not all concentrations are necessary. Recent studies suggest that including only two non-trivial concentrations may suffice to yield low errors for both energies and forces⁶⁶. Nevertheless, further investigations are needed to determine how reducing the number of concentrations might still yield accurate predictions of the oxidation state patterns. Future developments include automating the process, and ideally coupling it with Monte Carlo algorithms to simplify the combinatorial search of OSs patterns (especially required for supercells where the number of OSs patterns greatly increases). This process could be further streamlined by an initial guess based on the atomic environment geometry of redox-active elements.

  We believe this work could greatly enhance predictive, atomistic electrochemistry across multiple fields of science, due to its simplicity and generality. For example, we expect it to be easily extendable to scenarios where geometric information alone is clearly inadequate for determining oxidation states, such as in liquids.

Methods

The main ingredients

The DFT+U+V method was originally introduced in ref. ²⁹ as a generalization of DFT+U (in its simplest, rotationally invariant formulation introduced by Dudarev et al.²⁸) and it is based on an extended Hubbard model that contains both on-site U and intersite V electronic interactions. The physical rationale for the Hubbard U and V corrections lies in their ability to address spurious deviations from the piecewise linearity (PWL) of the DFT total energy with respect to the fractional addition or removal of charge^{67, 68, 69, 70, 71, 72–73}, which are associated with the self-interaction errors (SIEs). In DFT+U+V, such deviations from PWL are handled by adding an extended Hubbard corrective term E_U+V to the standard DFT energy E_DFT:

1

$E_{DFT+U+V}=E_{DFT}+E_{U+V}.$The extended Hubbard correction energy for a manifold with angular momentum l is written as:

2

$E_{U+V}=\sum _{I,\sigma }\frac{U^I}{2}Tr\left[{\mathbf{n}}_l^{II\sigma }\right(\mathbf{1}-{\mathbf{n}}_l^{II\sigma }\left)\right]-\sum _{I,J,\sigma }^{*}\frac{V^{IJ}}{2}Tr\left[{\mathbf{n}}_l^{IJ\sigma }{\mathbf{n}}_l^{JI\sigma }\right],$where I and J are atomic site indices, σ labels the spin of electrons. U^J and V^IJ are effective on-site and inter-site Hubbard parameters, respectively. The asterisk in the sum signifies that, for each atom I, the index J includes all its neighbors up to a given distance. The generalized occupation matrices ${\mathbf{n}}_l^{IJ\sigma }$ are derived from the projection of Kohn-Sham (KS) states onto localized atom-centered orbitals ${\varphi }_m^J$ (Hubbard projector functions):

3

${\mathbf{n}}_l^{IJ\sigma }\equiv n_{m,m^{\prime }}^{IJ\sigma }=\sum _{\mathbf{k},v}{f}_{\mathbf{k}v}^{\sigma }⟨{\psi }_{\mathbf{k}v}^{\sigma }\mid {\varphi }_{m^{\prime }}^J⟩⟨{\varphi }_m^I\mid {\psi }_{\mathbf{k}v}^{\sigma }⟩,$where v and σ represent the band and spin labels of the KS wavefuncitions ${\psi }_{\mathbf{k}v}^{\sigma }$, respectively, k denotes points in the first Brillouin zone, and $f_{\mathbf{k}v}^{\sigma }$ are the occupations of the KS states. We also point out that the choice of the projector function exerts a significant influence on the numerical values of the calculated Hubbard parameters and, more generally, on the prediction of material properties⁷⁴. The on-site U^I and intersite V^IJ terms work in opposition to each other. The on-site term encourages localization at atomic sites, reducing hybridization with neighboring atoms, while the inter-site term favors the formation of hybridized states involving neighboring atoms. As a result, the values of U^I and V^IJ play a crucial role in balancing localization and hybridization in Hubbard-corrected DFT. However, since they are not known a priori, they must be determined through an appropriate method.

While in the simplest cases these parameters could be obtained through semi-empirical tuning (but then negating the predictive power of the approach, and the capability to deal with complex and very diverse local environments, that require atom-specific U and V), unbiased predictions identify Hubbard parameters self-consistently through linear-response calculations^42,69,75, particularly efficient when density-functional perturbation theory (DFPT) is deployed^76,77. Moreover, it has been shown that jointly optimizing Hubbard parameters and the crystal structure, rather than relying on the equilibrium geometry obtained using (semi-)local functionals, can significantly improve the accuracy of the final properties of interest⁷⁴. To do so, a self-consistent procedure combining DFPT and structural optimizations can be used^76,78.

Finally, we mention that, from the total energy expressed in Eqs. (1) and (2), it is possible to calculate the extended Hubbard forces. Under the Born-Oppenheimer approximation, the force acting on atom I situated at position R_I is defined as:

4

${\mathbf{F}}_I=-\frac{d{E}_{DFT+U+V}}{d{\mathbf{R}}_I}.$Formally, the forces are computed by means of the Hellman-Feynman theorem. For details on the calculation, we refer to ref. ⁷⁴, and we note that fixed average U and V do not require the calculation of $\frac{dU}{dR}$ or $\frac{dV}{dR}$. The forces of Eq. (4) are then used to update the nuclear positions according to the laws of classical dynamics in the Born-Oppenheimer molecular dynamics.

In summary, we have outlined the essential ingredients of this work. Eqs. (1) and (4) provide the energies and forces necessary for training the machine learning potential, while Eq. (3), the occupation matrix, yields—once diagonalized—the electronic occupations of the Mn 3d shells, allowing us to infer its OS.

Computational parameters

All the calculations were performed using the plane-wave pseudopotential implementation of DFT contained in the Quantum ESPRESSO distribution^79,80. We use the PBEsol exchange-correlation functional⁸¹ and pseudopotentials from the SSSP library v1.1⁸², which are either ultrasoft (US) or projector-augmented-wave (PAW). For manganese, we have used mn_pbesol_v1.5.uspp.F.UPF from the GBRV v1.5 library⁸³, for oxygen O.pbesol-n-kjpaw_psl.0.1.UPF from the Pslibrary v0.3.1⁸⁴, for phosphorus P.pbesol-n-rrkjus_psl.1.0.0.UPF from the Pslibrary v1.0.0⁸⁴, and for lithium li_pbesol_v1.4.uspp.F.UPF from the GBRV v1.4 library⁸³. All DFT+U+V calculations are performed using orthogonalized atomic orbitals as described in detail in ref. ⁸⁵.

DFT+U+V FPMD is performed using the pw. x code of Quantum ESPRESSO. We used averaged self-consistent parameters U = 5.1 eV and V = 0.7 eV. We use the uniform Γ-centered k-points grid of size 3 × 4 × 5. The KS wavefunctions and potentials are expanded in PWs up to a kinetic-energy cutoff of 65 and 780 Ry, respectively. We used a Gaussian smearing with a broadening parameter of 0.005 Ry. We conducted the FPMD simulations on the previous structurally optimized cells with atomic positions relaxed at T = 0 K. The Verlet algorithm is used to integrate the classical equations of motion with a timestep of 4 fs. The dynamics was performed in the NVT ensemble (at fixed number of particles N, volume V and temperature T) by using the stochastic-velocity rescaling algorithm⁶⁵.

For the structural optimizations, the Brillouin zone was sampled using the uniform Γ-centered k point mesh of size 5 × 8 × 9. KS wavefunctions and potentials are expanded in PWs up to a kinetic-energy cutoff of 90 and 1080 Ry, respectively. The crystal structure was optimized using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm⁸⁶, with a convergence threshold for the total energy of 10⁻⁶ Ry, for forces of 10⁻⁵ Ry/Bohr, and for pressure of 0.5 Kbar.

The DFPT calculations of Hubbard parameters are performed using the hp. x code⁸⁵ of Quantum ESPRESSO using the uniform Γ-centered k and q point meshes of size 3 × 4 × 5 and 1 × 2 × 3, respectively. The KS wavefunctions and potentials are expanded in PWs up to a kinetic-energy cutoff of 65 and 780 Ry, respectively, for the calculation of Hubbard parameters.

ML architecture

The model architecture features a convolution filter with a cutoff radius of 4 Å and consists of four interaction blocks. It employs a maximum rotation order of 2, with feature multiplicities set to 32, including features of odd mirror parity. The radial neural network comprises eight basis radial functions, two radial invariant layers, and 64 hidden neurons. Training was performed over 4200 epochs using the Adam optimizer.

Acknowledgements

We thank S. Muy and I. Timrov for useful discussions. We acknowledge support by the European Commission through the MaX Centre of Excellence for supercomputing applications (grant number 101093374). We gratefully acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) under Germany’s Excellence Strategy (EXC 2077, No. 390741603, University Allowance, University of Bremen) and Lucio Colombi Ciacchi, the host of the “U Bremen Excellence Chair Program.” We acknowledge support by the NCCR MARVEL, a National Centre of Competence in Research, funded by the Swiss National Science Foundation (grant number 205602).

Author contributions

C.M.: conceptualization, methodology, software, formal analysis, investigation, data curation, writing—original draft, visualization; N.M.: supervision, conceptualization, methodology, project administration, funding acquisition. All authors: writing—review & editing.

Funding

Open Access funding enabled and organized by Projekt DEAL.

Data availability

The data used to produce the results of this work are available in the Materials Cloud Archive⁸⁷.

Code availability

DFT+U+V calculations are performed using QUANTUM ESPRESSO v7.2, which is open-source and can be freely downloaded from https://www.quantum-espresso.org. The machine learning framework used is NequIP, which is open-source and available at https://github.com/mir-group/nequip.

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41524-025-01709-z.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Walsh, A; Sokol, AA; Buckeridge, J; Scanlon, DO; Catlow, CRA. Oxidation states and ionicity. Nat. Mater.; 2018; 17, pp. 958-964.1:CAS:528:DC%2BC1cXhvVCqtrjJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30275565][DOI: https://dx.doi.org/10.1038/s41563-018-0165-7]

2. Barton, JL. Electrification of the chemical industry. Science; 2020; 368, pp. 1181-1182.1:CAS:528:DC%2BB3cXhtFGitLbM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32527815][DOI: https://dx.doi.org/10.1126/science.abb8061]

3. Jablonka, KM; Ongari, D; Moosavi, SM; Smit, B. Using collective knowledge to assign oxidation states of metal cations in metal–organic frameworks. Nat. Chem.; 2021; 13, pp. 771-777.1:CAS:528:DC%2BB3MXhsV2ksrjN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34226703][DOI: https://dx.doi.org/10.1038/s41557-021-00717-y]

4. Nykvist, B; Nilsson, M. Rapidly falling costs of battery packs for electric vehicles. Nat. Clim. Change; 2015; 5, pp. 329-332. [DOI: https://dx.doi.org/10.1038/nclimate2564]

5. Sarma, DD; Shukla, AK. Building better batteries: a travel back in time. ACS Energy Lett.; 2018; 3, pp. 2841-2845.1:CAS:528:DC%2BC1cXitVKkt7bL [DOI: https://dx.doi.org/10.1021/acsenergylett.8b01966]

6. Ponrouch, A; Palacín, MR. Post-Li batteries: promises and challenges. Philos. Trans. R. Soc. A; 2019; 377, 37720180297. [DOI: https://dx.doi.org/10.1098/rsta.2018.0297]

7. Liu, J et al. Pathways for practical high-energy long-cycling lithium metal batteries. Nat. Energy; 2019; 4, pp. 180-186.1:CAS:528:DC%2BC1MXmslKjsr0%3D [DOI: https://dx.doi.org/10.1038/s41560-019-0338-x]

8. Hohenberg, P; Kohn, W. Inhomogeneous electron gas. Phys. Rev.; 1964; 136, B864. [DOI: https://dx.doi.org/10.1103/PhysRev.136.B864]

9. Kohn, W; Sham, L. Self-consistent equations including exchange and correlation effects. Phys. Rev.; 1965; 140, A1133. [DOI: https://dx.doi.org/10.1103/PhysRev.140.A1133]

10. Car, R; Parrinello, M. Unified approach for molecular dynamics and density-functional theory. Phys. Rev. Lett.; 1985; 55, pp. 2471-2474.1:CAS:528:DyaL28XhtVOlug%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10032153][DOI: https://dx.doi.org/10.1103/PhysRevLett.55.2471]

11. McNaught, A. & Wilkinson, A. Compendium of Chemical Terminology—The “Gold Book” 2nd edn (Blackwell, 1997).

12. Bader, R. Atoms in Molecules. A Quantum Theory (Oxford University Press, 1990).

13. Bickelhaupt, FM; van Eikema Hommes, NJR; Fonseca Guerra, C; Baerends, EJ. The carbon-lithium electron pair bond in (CH₃Li)_n (n = 1, 2, 4). Organometallics; 1996; 15, pp. 2923-2931.1:CAS:528:DyaK28XjtlKntbg%3D [DOI: https://dx.doi.org/10.1021/om950966x]

14. Raebiger, H; Lany, S; Zunger, A. Charge self-regulation upon changing the oxidation state of transition metals in insulators. Nature; 2008; 453, pp. 763-766.1:CAS:528:DC%2BD1cXmvVGmsb8%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18528391][DOI: https://dx.doi.org/10.1038/nature07009]

15. Resta, R. Charge states in transition. Nature; 2008; 453, pp. 735-735.1:CAS:528:DC%2BD1cXmvVGms70%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18528387][DOI: https://dx.doi.org/10.1038/453735a]

16. Jiang, L; Levchenko, SV; Rappe, AM. Rigorous definition of oxidation states of ions in solids. Phys. Rev. Lett.; 2012; 108, 166403. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22680739][DOI: https://dx.doi.org/10.1103/PhysRevLett.108.166403]

17. Grasselli, F; Baroni, S. Topological quantization and gauge invariance of charge transport in liquid insulators. Nat. Phys.; 2019; 15, pp. 967-972.1:CAS:528:DC%2BC1MXhtlamtrvK [DOI: https://dx.doi.org/10.1038/s41567-019-0562-0]

18. Resta, R. & Vanderbilt, D. Theory of Polarization: A Modern Approach, 31–68 (Springer, 2007).

19. Vanderbilt, D. Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators (Cambridge University Press, 2018).

20. Sit, PH-L; Car, R; Cohen, MH; Selloni, A. Simple, unambiguous theoretical approach to oxidation state determination via first-principles calculations. Inorg. Chem.; 2011; 50, 10259.1:CAS:528:DC%2BC3MXhtFOitLzE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21923086][DOI: https://dx.doi.org/10.1021/ic2013107]

21. Pegolo, P; Grasselli, F; Baroni, S. Oxidation states, Thouless’ pumps, and nontrivial ionic transport in nonstoichiometric electrolytes. Phys. Rev. X; 2020; 10, 041031.1:CAS:528:DC%2BB3MXjtFSksbo%3D

22. Timrov, I; Aquilante, F; Cococcioni, M; Marzari, N. Accurate electronic properties and intercalation voltages of olivine-type li-ion cathode materials from extended hubbard functionals. Phys. Rev. X Energy; 2022; 1, 033003.

23. Hautier, G et al. Phosphates as lithium-ion battery cathodes: an evaluation based on high-throughput ab initio calculations. Chem. Mat.; 2011; 23, pp. 3495-3508.1:CAS:528:DC%2BC3MXovVKrur0%3D [DOI: https://dx.doi.org/10.1021/cm200949v]

24. Hautier, G et al. Novel mixed polyanions lithium-ion battery cathode materials predicted by high-throughput ab initio computations. J. Mater. Chem.; 2011; 21, 17147.1:CAS:528:DC%2BC3MXhtlKgtLrI [DOI: https://dx.doi.org/10.1039/c1jm12216a]

25. Anisimov, V; Zaanen, J; Andersen, O. Band theory and Mott insulators: Hubbard U instead of Stoner I. Phys. Rev. B; 1991; 44, 943.1:CAS:528:DyaK3MXkvVaqsLs%3D [DOI: https://dx.doi.org/10.1103/PhysRevB.44.943]

26. Liechtenstein, A; Anisimov, V; Zaanen, J. Density-functional theory and strong interactions: orbital ordering in Mott-Hubbard insulators. Phys. Rev. B; 1995; 52, R5467.1:CAS:528:DyaK2MXnslOisbw%3D [DOI: https://dx.doi.org/10.1103/PhysRevB.52.R5467]

27. Anisimov, V; Aryasetiawan, F; Lichtenstein, A. First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA+U method. J. Phys. Condens. Matter; 1997; 9, 767.1:CAS:528:DyaK2sXhtFalt78%3D [DOI: https://dx.doi.org/10.1088/0953-8984/9/4/002]

28. Dudarev, S; Botton, G; Savrasov, S; Humphreys, C; Sutton, A. Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+U study. Phys. Rev. B; 1998; 57, 1505.1:CAS:528:DyaK1cXlsVarsQ%3D%3D [DOI: https://dx.doi.org/10.1103/PhysRevB.57.1505]

29. Campo Jr, VL; Cococcioni, M. Extended DFT+U+V method with on-site and inter-site electronic interactions. J. Phys. Condens. Matter; 2010; 22, 055602. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21386347][DOI: https://dx.doi.org/10.1088/0953-8984/22/5/055602]

30. Tancogne-Dejean, N; Rubio, A. Parameter-free hybridlike functional based on an extended Hubbard model: DFT+U+V. Phys. Rev. B; 2020; 102, 155117.1:CAS:528:DC%2BB3cXitlCht7%2FP [DOI: https://dx.doi.org/10.1103/PhysRevB.102.155117]

31. Lee, S-H; Son, Y-W. First-principles approach with a pseudohybrid density functional for extended Hubbard interactions. Phys. Rev. Res.; 2020; 2, 043410.1:CAS:528:DC%2BB3MXjtVygtb4%3D [DOI: https://dx.doi.org/10.1103/PhysRevResearch.2.043410]

32. Sun, J; Ruzsinszky, A; Perdew, J. Strongly constrained and appropriately normed semilocal density functional. Phys. Rev. Lett.; 2015; 115, 036402. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26230809][DOI: https://dx.doi.org/10.1103/PhysRevLett.115.036402]

33. Bartók, A; Yates, J. Regularized SCAN functional. J. Chem. Phys.; 2019; 150, 161101. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31042928][DOI: https://dx.doi.org/10.1063/1.5094646]

34. Furness, J; Kaplan, A; Ning, J; Perdew, J; Sun, J. Accurate and numerically efficient r²SCAN meta-generalized gradient approximation. J. Phys. Chem. Lett.; 2020; 11, 8208.1:CAS:528:DC%2BB3cXhslequ77N [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32876454][DOI: https://dx.doi.org/10.1021/acs.jpclett.0c02405]

35. Gautam, G; Carter, E. Evaluating transition metal oxides within DFT-SCAN and SCAN+U frameworks for solar thermochemical applications. Phys. Rev. Mater.; 2018; 2, 095401. [DOI: https://dx.doi.org/10.1103/PhysRevMaterials.2.095401]

36. Long, O; Gautam, G; Carter, E. Evaluating optimal U for 3d transition-metal oxides within the SCAN+U framework. Phys. Rev. Mater.; 2020; 4, 045401.1:CAS:528:DC%2BB3cXhtVeksbzJ [DOI: https://dx.doi.org/10.1103/PhysRevMaterials.4.045401]

37. Kaczkowski, J; Pugaczowa-Michalska, M; Pĺowaś-Korus, I. Comparative density functional studies of pristine and doped bismuth ferrite polymorphs by GGA+U and meta-GGA SCAN+U. Phys. Chem. Chem. Phys.; 2021; 23, 8571.1:CAS:528:DC%2BB3MXltlegu7k%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33876019][DOI: https://dx.doi.org/10.1039/D0CP06157C]

38. Artrith, N; Torres, J; Urban, A; Hybertsen, M. Data-driven approach to parameterize SCAN+U for an accurate description of 3d transition metal oxide thermochemistry. Phys. Rev. Mater.; 2022; 6, 035003.1:CAS:528:DC%2BB38XpsFKmtrk%3D [DOI: https://dx.doi.org/10.1103/PhysRevMaterials.6.035003]

39. Adamo, C; Barone, V. Toward reliable density functional methods without adjustable parameters: the PBE0 model. J. Chem. Phys.; 1999; 110, 6158.1:CAS:528:DyaK1MXitVCmt7Y%3D [DOI: https://dx.doi.org/10.1063/1.478522]

40. Heyd, J; Scuseria, G; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys.; 2003; 118, 8207.1:CAS:528:DC%2BD3sXjtlSisLw%3D [DOI: https://dx.doi.org/10.1063/1.1564060]

41. Heyd, J; Scuseria, G; Ernzerhof, M. Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)]. J. Chem. Phys.; 2006; 124, 219906. [DOI: https://dx.doi.org/10.1063/1.2204597]

42. Cococcioni, M; Marzari, N. Energetics and cathode voltages of LiMPO₄ olivines (M = Fe, Mn) from extended Hubbard functionals. Phys. Rev. Mater.; 2019; 3, 033801. [DOI: https://dx.doi.org/10.1103/PhysRevMaterials.3.033801]

43. Sit, PH-L; Cococcioni, M; Marzari, N. Realistic quantitative descriptions of electron transfer reactions: diabatic free-energy surfaces from first-principles molecular dynamics. Phys. Rev. Lett.; 2006; 97, 028303. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16907484][DOI: https://dx.doi.org/10.1103/PhysRevLett.97.028303]

44. Newton, MD; Sutin, N. Electron transfer reactions in condensed phases. Annu. Rev. Phys. Chem.; 1984; 35, pp. 437-480.1:CAS:528:DyaL2cXmtFKmt7Y%3D [DOI: https://dx.doi.org/10.1146/annurev.pc.35.100184.002253]

45. Marcus, R; Sutin, N. Electron transfers in chemistry and biology. Biochim. Biophys. Acta; 1985; 811, pp. 265-322.1:CAS:528:DyaL2MXltFygs78%3D [DOI: https://dx.doi.org/10.1016/0304-4173(85)90014-X]

46. Warshel, A; Parson, WW. Computer simulations of electron-transfer reactions in solution and in photosynthetic reaction centers. Annu. Rev. Phys. Chem.; 1991; 42, pp. 279-309.1:CAS:528:DyaK38XhsF2ksrk%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/1747189][DOI: https://dx.doi.org/10.1146/annurev.pc.42.100191.001431]

47. Behler, J; Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys. Rev. Lett.; 2007; 98, 146401. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17501293][DOI: https://dx.doi.org/10.1103/PhysRevLett.98.146401]

48. Batzner, S et al. E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nat. Commun.; 2022; 13, 1:CAS:528:DC%2BB38Xht1ShtLrE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35508450][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9068614][DOI: https://dx.doi.org/10.1038/s41467-022-29939-5] 2453.

49. Batatia, I., Kovacs, D. P., Simm, G., Ortner, C. & Csanyi, G. Mace: higher order equivariant message passing neural networks for fast and accurate force fields. Preprint at arXivhttps://doi.org/10.48550/arXiv.2206.07697 (2022).

50. Zhang, L; Han, J; Wang, H; Car, R; E, W. Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics. Phys. Rev. Lett.; 2018; 120, 143001.1:CAS:528:DC%2BC1MXltFSksrg%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29694129][DOI: https://dx.doi.org/10.1103/PhysRevLett.120.143001]

51. Reed, J; Ceder, G. Role of electronic structure in the susceptibility of metastable transition-metal oxide structures to transformation. Chem. Rev.; 2004; 104, pp. 4513-4534.1:CAS:528:DC%2BD2cXnvVKltrk%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15669161][DOI: https://dx.doi.org/10.1021/cr020733x]

52. Rappe, AK; Goddard III, WA. Charge equilibration for molecular dynamics simulations. J. Phys. Chem.; 1991; 95, pp. 3358-3363.1:CAS:528:DyaK3MXhvFaktbo%3D [DOI: https://dx.doi.org/10.1021/j100161a070]

53. Ghasemi, SA; Hofstetter, A; Saha, S; Goedecker, S. Interatomic potentials for ionic systems with density functional accuracy based on charge densities obtained by a neural network. Phys. Rev. B; 2015; 92, 045131. [DOI: https://dx.doi.org/10.1103/PhysRevB.92.045131]

54. Ko, TW; Finkler, JA; Goedecker, S; Behler, J. A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer. Nat. Commun.; 2021; 12, 1:CAS:528:DC%2BB3MXhsl2qtrg%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33452239][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7811002][DOI: https://dx.doi.org/10.1038/s41467-020-20427-2] 398.

55. Staacke, CG et al. Kernel charge equilibration: efficient and accurate prediction of molecular dipole moments with a machine-learning enhanced electron density model. Mach. Learn. Sci. Technol.; 2022; 3, 015032. [DOI: https://dx.doi.org/10.1088/2632-2153/ac568d]

56. Shaidu, Y; Pellegrini, F; Küçükbenli, E; Lot, R; de Gironcoli, S. Incorporating long-range electrostatics in neural network potentials via variational charge equilibration from shortsighted ingredients. NPJ Comput. Mater.; 2024; 10, 47.1:CAS:528:DC%2BB2cXltVGkurg%3D [DOI: https://dx.doi.org/10.1038/s41524-024-01225-6]

57. Chen, X. et al. Atomistic learning in the electronically grand-canonical ensemble. NPJ Comput. Mater. 9, 73 (2023).

58. Kocer, E., Haouari, R. E., Dellago, C. & Behler, J. Machine learning potentials for redox chemistry in solution. Preprint at https://arxiv.org/abs/2410.03299 (2024).

59. Deng, B., Ahong, P., Jun, K., Riebesell, J. & Han, K. CHGNet as a pretrained universal neural network potential for charge-informed atomistic modelling. Nat. Mach. Intell. 5, 1031–1041 (2023).

60. Eckhoff, M; Lausch, KN; Blöchl, PE; Behler, J. Predicting oxidation and spin states by high-dimensional neural networks: applications to lithium manganese oxide spinels. J. Chem. Phys.; 2020; 153, 164107.1:CAS:528:DC%2BB3cXitFyit7rP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33138439][DOI: https://dx.doi.org/10.1063/5.0021452]

61. Uhrin, M., Zadoks, A., Binci, L., Marzari, N. & Timrov, I. Machine learning Hubbard parameters with equivariant neural networks. NPJ Comput. Mater.11, 19 (2024).

62. Muraliganth, T; Manthiram, A. Understanding the shifts in the redox potentials of olivines LiM_1−yM_yPO₄ (M = Fe, Mn, Co, and Mg) solid solution cathodes. J. Phys. Chem. C; 2010; 114, 15530.1:CAS:528:DC%2BC3cXhtVKrtrrI [DOI: https://dx.doi.org/10.1021/jp1055107]

63. Kobayashi, G et al. Shift of redox potential and kinetics in Li_x(Mn_yFe_1-y)PO₄. J. Power Sources; 2009; 189, pp. 397-401.1:CAS:528:DC%2BD1MXjsVegs7g%3D [DOI: https://dx.doi.org/10.1016/j.jpowsour.2008.07.085]

64. Kam, R. L. et al. Interplay between electron localization, magnetic order, and Jahn-Teller distortion dictates LiMnO₂ phase stability. Phys. Rev. B111, 245132 (2025).

65. Bussi, G; Donadio, D; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys.; 2007; 126, 014101. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17212484][DOI: https://dx.doi.org/10.1063/1.2408420]

66. Goodwin, ZAH et al. Transferability and accuracy of ionic liquid simulations with equivariant machine learning interatomic potentials. J. Phys. Chem. Lett.; 2024; 15, pp. 7539-7547.1:CAS:528:DC%2BB2cXhsFKktrnN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39023916][DOI: https://dx.doi.org/10.1021/acs.jpclett.4c01942]

67. Perdew, JP; Parr, RG; Levy, M; Balduz, JL. Density-functional theory for fractional particle number: derivative discontinuities of the energy. Phys. Rev. Lett.; 1982; 49, pp. 1691-1694.1:CAS:528:DyaL3sXit1CisQ%3D%3D [DOI: https://dx.doi.org/10.1103/PhysRevLett.49.1691]

68. Cococcioni, M. A LDA+U Study of Selected Iron Compounds. PhD thesis, Scuola Internazionale Superiore di Studi Avanzati (2002).

69. Cococcioni, M; de Gironcoli, S. Linear response approach to the calculation of the effective interaction parameters in the LDA+U method. Phys. Rev. B; 2005; 71, 035105. [DOI: https://dx.doi.org/10.1103/PhysRevB.71.035105]

70. Kulik, H; Cococcioni, M; Scherlis, D; Marzari, N. Density functional theory in transition-metal chemistry: a self-consistent Hubbard U approach. Phys. Rev. Lett.; 2006; 97, 103001. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17025809][DOI: https://dx.doi.org/10.1103/PhysRevLett.97.103001]

71. Mori-Sánchez, P; Cohen, A; Yang, W. Many-electron self-interaction error in approximate density functionals. J. Chem. Phys.; 2006; 125, 201102. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17144681][DOI: https://dx.doi.org/10.1063/1.2403848]

72. Mori-Sánchez, P; Cohen, AJ; Yang, W. Discontinuous nature of the exchange-correlation functional in strongly correlated systems. Phys. Rev. Lett.; 2009; 102, 066403. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19257614][DOI: https://dx.doi.org/10.1103/PhysRevLett.102.066403]

73. Zhao, Q; Ioannidis, EI; Kulik, HJ. Global and local curvature in density functional theory. J. Chem. Phys.; 2016; 145, 054109. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27497541][DOI: https://dx.doi.org/10.1063/1.4959882]

74. Timrov, I; Aquilante, F; Binci, L; Cococcioni, M; Marzari, N. Pulay forces in density-functional theory with extended Hubbard functionals: from nonorthogonalized to orthogonalized manifolds. Phys. Rev. B; 2020; 102, 235159.1:CAS:528:DC%2BB3MXivVOlsbY%3D [DOI: https://dx.doi.org/10.1103/PhysRevB.102.235159]

75. Shishkin, M; Sato, H. Self-consistent parametrization of DFT+U framework using linear response approach: application to evaluation of redox potentials of battery cathodes. Phys. Rev. B; 2016; 93, 085135. [DOI: https://dx.doi.org/10.1103/PhysRevB.93.085135]

76. Timrov, I; Marzari, N; Cococcioni, M. Self-consistent Hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formulations. Phys. Rev. B; 2021; 103, 045141.1:CAS:528:DC%2BB3MXkt1Grt7c%3D [DOI: https://dx.doi.org/10.1103/PhysRevB.103.045141]

77. Timrov, I; Marzari, N; Cococcioni, M. Hubbard parameters from density-functional perturbation theory. Phys. Rev. B; 2018; 98, 085127. [DOI: https://dx.doi.org/10.1103/PhysRevB.98.085127]

78. Hsu, H; Umemoto, K; Cococcioni, M; Wentzcovitch, R. First-principles study for low-spin LaCoO₃ with a structurally consistent Hubbard U. Phys. Rev. B; 2009; 79, 125124. [DOI: https://dx.doi.org/10.1103/PhysRevB.79.125124]

79. Giannozzi, P et al. Quantum ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys. Condens. Matter; 2009; 21, 395502. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21832390][DOI: https://dx.doi.org/10.1088/0953-8984/21/39/395502]

80. Giannozzi, P et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys. Condens. Matter; 2017; 29, 465901.1:STN:280:DC%2BC1M7jvFOjuw%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29064822][DOI: https://dx.doi.org/10.1088/1361-648X/aa8f79]

81. Perdew, J. et al. Restoring the density-gradient expansion for exchange in solids and surfaces. Phys. Rev. Lett. 100, 136406 (2008).

82. Bosoni, E. er al. How to verify the precision of density-functional-theory implementations via reproducible and universal workflows. Nat. Rev. Phys. 6, 45–58 (2023).

83. Garrity, K., Bennett, J., Rabe, K. & Vanderbilt, D. Pseudopotentials for high-throughput DFT calculations. Comput. Mater. Sci. 81, 446 (2014).

84. Dal Corso, A. Pseudopotentials periodic table: from H to Pu. Comput. Mater. Sci.; 2014; 95, 337.1:CAS:528:DC%2BC2cXhtlKht7vP [DOI: https://dx.doi.org/10.1016/j.commatsci.2014.07.043]

85. Timrov, I; Marzari, N; Cococcioni, M. Comput. Phys. Commun.; 2022; 279, 108455.1:CAS:528:DC%2BB38Xhs1Ckur%2FE [DOI: https://dx.doi.org/10.1016/j.cpc.2022.108455]

86. Fletcher, R. Practical Methods of Optimization 2nd edn (Wiley, 1987).

87. Malica, C. & Marzari, N. Teaching oxidation states to neural networks. Materials Cloud Archive https://doi.org/10.24435/materialscloud:w7-k1 (2024).