Machine Learning-Aided First-Principles Calculations of Redox Potentials

Redox potentials of electron transfer reactions are of fundamental importance for the performance and description of electrochemical devices. Despite decades of research, accurate computational predictions for the redox potential of even simple metals remain very challenging. Here we use a combination of first principles calculations and machine learning to predict the redox potentials of three redox couples, $\mathrm{Fe}^{2+}$/$\mathrm{Fe}^{3+}$, $\mathrm{Cu}^{+}$/$\mathrm{Cu}^{2+}$ and $\mathrm{Ag}^{+}$/$\mathrm{Ag}^{2+}$. Using a hybrid functional with a fraction of 25\% exact exchange (PBE0) the predicted values are 0.92, 0.26 and 1.99 V in good agreement with the best experimental estimates (0.77, 0.15, 1.98 V). We explain in detail, how we combine machine learning, thermodynamic integration from machine learning to semi-local functionals, as well as a combination of thermodynamic perturbation theory and $\Delta$-machine learning to determine the redox potentials for computationally expensive hybrid functionals. The combination of these approaches allows one to obtain statistically accurate results.


INTRODUCTION
Green energy and a circular economy are one of the key paradigms that our human society needs to realize in the next few decades.This implies that we need to give up on the combustion of fossil fuels.A key element to achieve this paradigm shift is the use of electrochemistry, be it for batteries and fuel cells, to convert electrical energy to hydrogen or other valuable chemicals, or to convert hydrogen back to energy without direct combustion in air.
The redox potential of electron transfer (ET), Ox + ne − → Red in liquids, is an essential property for a variety of electrochemical energy conversion devices, such as batteries, fuel cells, and electrochemical fuel synthesis.It determines the alignment of redox levels relative to the Fermi level of a metal, or valence band maximum (VBM) and conduction band minimum (CBM) of semiconductor and insulator electrodes.It also determines the stability windows of ions and molecules in solutions, that is the range of voltages within which a specific ion or molecule can undergo electrochemical reactions.This information is vital to design redox species and solvent molecules, such as redox couples for redox-flow batteries [1], solvents and additives for Li-ion batteries [2][3][4], radical scavengers for fuel cells [5] and electrocatalysts for fuel synthesis [6,7].
and applied this method to several redox reactions in aqueous solutions [9,15].They discovered that the use of a semi-local functional leads to errors exceeding 0.5 V.This discrepancy arises because the functional inaccurately yields the shallow valence band edge and the deep conduction band edge, resulting in incorrect hybridization with the redox levels.Similar magnitudes of errors have also been observed in other FP calculations that employ semi-local approximations [16,17].As a result, Sprik and co-workers opted for a hybrid functional.Nonetheless, they observed a significant spread of values for two metal ion couples, with the Cu 2+ /Cu + couple ranging from −1.13 to −0.20 V (experimental value 0.16 V) and the Ag 2+ /Ag + couple ranging from 0.90 to 1.72 V (experimental value 1.98 V) [9].These variations were attributed to differences in the pseudopotential and the computational code base (CMPD versus CP2K).While the "best" values obtained using the hybrid functional and highly accurate pseudopotentials are relatively close to experimental values (−0.20 V for Cu, and 1.72 V for Ag), the agreement is still far from being quantitative.Due to the high computational cost of hybrid functionals, most calculations have been performed using approximated methods, such as continuum solvation models [18][19][20][21] and QM/MM models [22,23].Although these models can reproduce the experimental redox potentials of ions and molecules with convenient accuracy, the computational results heavily rely on many approximations, making it unclear which predictions are strictly correct.Here, we briefly note that these FP and approximated methods have been extended to electrochemistry at liquid-solid interfaces [24][25][26][27][28][29][30][31][32].Nowadays, these methods have become indispensable for elucidating electrochemical interfacial phenomena and designing advanced materials [33][34][35][36][37][38].However, even in the calculation of redox reactions at interfaces, approximations are made in most of applications, such as representing the motion of atomic nuclei with simple statistical models like the harmonic oscillator model [24,25,33,34,36,39,40], modeling solvents by molecular mechanics [29], or modeling by continuum mediums [26-28, 30, 32].A rigorous FP method that eliminates these approximations is also desired in the field of interfacial electrochemistry.
The main goals of the present work are three-fold: First, we want to accurately calculate the redox potential of metal ions in water for three prototypical cases: Ag, Cu, and Fe.Ag 2+ ions are among the most aggressive oxidants with a large redox potential, whereas the redox potential of Cu 2+ ions is fairly shallow, and the Fe 3+ /Fe 2+ reaction lies in between.The first two redox reactions involve large changes in the ion water coordination, which makes the calculation challenging, whereas the redox reaction of Fe is a so called simple outer sphere ET reaction and has been the subject of numerous experimental and theoretical studies [42].The Fe ions are conceived to be particularly challenging for density functional theory.Second, we want to establish a computationally feasible pathway that yields statistically accurate results.Last but not least, we want to systematically explore different density functionals to set a guideline for future studies.

OVERVIEW OF THEORY AND MODELLING
We begin with an overview of the used theory and modeling.Further details can be found in the Methods section.The reactions evaluated in this study are electron transfer reactions in water: Fe 3+ + e − ↔ Fe 2+ , Cu 2+ + e − ↔ Cu + , and Ag 2+ + e − ↔ Ag + .Here, we assume that other side reactions do not occur, and only the valency of redox species changes due to the reaction similarly to the previous study [9].
The redox potential U redox is determined by the Gibbs free energy difference ∆A between the reduced and oxidized states as where e is the elementary charge and n is the number of electrons involved in the reaction.The free energy difference ∆A can be exactly determined by thermodynamic integration (TI) [43,44]: Here, ⟨X⟩ λ denotes the expectation value of X for an ensemble created by the Hamiltonian at coupling λ .The ensemble required for computing the Gibbs free energy difference is the isothermal-isobaric ensemble.However, in practice, due to the minor impact of volumetric change caused by the reaction, we utilized the canonical ensemble with the empirical density of liquid water, similar to previous publications [9,13].The integral seamlessly connects the oxidized state (λ = 0) to the reduced state (λ = 1) along a coupling path [45,46].The potential energy surface upon which atoms move is described by the grand potential Ω of the system opened for electrons [47].
Consequently, the Hamiltonian of the system is described as follows: where N a is the number of atoms, m i and p i are the mass and momentum vector of i-th atom, and µ and N are the chemical potential and the number of electrons.The chemical potential µ is fixed at the reservoir level, whereas N varies by n along the coupling path.U represents the potential energy surface at λ , equating to the sum of the Helmholtz free energy of the electronic subsystem and the electrostatic interactions among nuclei.Following previous studies [45,46], U can be described as where U 0 and U 1 are the potential energies of the oxidized and reduced states, respectively.Hence, the free energy difference ∆A is written as If the structural changes are significant from the oxidized to the reduced species -recall that this is the case for Ag and Cu -many integration steps are required to accurately determine the energy difference.The application of this approach entails two difficulties.(i) Clearly, it implies huge computational cost if applied directly to hybrid functionals; if 100.000 timesteps using a complete plane wave basis set are required to obtain good statistical accuracy, several 10 mio core hours are necessary.(ii) Second, during the reaction one electron needs to be transferred from the reservoir, characterizing the chemical potential of the electrons.The vacuum level is the best suited reference chemical potential that allows one to align the redox levels and band edges of the electrode in the absolute potential scale.However, in FP calculations of bulk systems under periodic boundary conditions, the vacuum level is a quantity that cannot be directly accessed during simulations.Chemical potential of electrons: We will address the second point (ii) first.Jiao and co-workers [48] suggested to use the average electrostatic potential as suitable reference point, and Leung [49] calculated the position of the average electrostatic potential with respect to the vacuum level in a second independent calculation involving a water slab.We refine this approach in a conceptually easy to understand way that simultaneously reduces finite size errors.As a reference, instead of using the vacuum level, we employ the O 1s level of water, which is fixed relative to the vacuum level and can be conveniently calculated with the FP code used in this study.Our approach is schematically illustrated in Fig. 1.In FP calcu-lations of a solution system under a periodic boundary condition, the energy contribution ⟨U 1 −U 0 ⟩ λ in Eq.( 6) is equal to the negative electron affinity of the oxidized species scaled to the average local potential of the system.The same calculation can also determine the O 1s level ε 1s,bulk λ of water, sufficiently far from the redox species and unaffected by the reactant, scaled to the average local potential.Therefore, measuring the redox level using the O 1s level as a reference results in ⟨U 1 −U 0 ⟩ λ /n − ε 1s,bulk , as highlighted in orange letters in Fig. 1 (c).In practice, ε 1s,bulk λ may slightly vary along the coupling path due to finite size effects [refer to Table S4 in the Supporting Information (SI)].By aligning the potentials between the 'defect' and the 'host' within the same supercell in this manner, the finite size effects can be mitigated [50,51].The vacuum level referenced to the O 1s level can be calculated using a slab model.As depicted in Fig. 1 (b), when referencing the O 1s level of water molecules located in the middle layer of the water slab, the vacuum level can be expressed as µ − ε 1s,slab , as indicated in blue letters in Fig. 1  (c).The difference between the redox level and vacuum level scaled to the O 1s level results in the redox level scaled to the vacuum level, as shown in red letters in Fig. 1 (c).Consequently, the free energy difference ∆A on an absolute scale is written as where the vacuum level µ is set to zero.As illustrated by the green letters in Fig. 1 (c), ∆ φ accounts for the difference between the local potential at the vacuum in the slab model and the one in the bulk solution model.Equation ( 7) is similar to the one used for the computational standard hydrogen electrode (CSHE) method proposed in the previous studies [13,14].However, there are three differences in our approach compared to the previous method.Firstly, our approach computes the absolute vacuum reference instead of referencing to the standard hydrogen electrode (SHE), thereby enabling the evaluation of absolute potentials of half-cell reactions.Second, aligning using the O 1s level in the same supercell is expected to mitigate the finite size effects.Third, machine learned (ML) force fields (FFs) can create many statistically independent configurations for the water slab.We do this by on-the-fly learning an H 2 O force field for the bulk and then for the surface, and performing finally extensive million step (total 1.5 ns) ML molecular dynamics for the surface.From this simulation we draw 3000 statistically independent snapshots.Only for these snapshots, FP calculations are performed to determine the average O 1s level with respect to the vacuum level.This substantially reduces the required computational time from 1 mio core hours for brute force runs using the semi-local functional to only 2200 core hours for the FP calculations on 3000 structures, including the ML simula- Thermodynamic integration: To address the problem to compute the free energy difference, i.e. point (i), we propose the ML-aided scheme as depicted in Fig. 2. Here, we use the abbreviations FP nl (Ox/Red), FP sl (Ox/Red), and ML(Ox/Red) to denote calculations using a non-local hybrid functional, a semi-local functional and machine-learned force field for the oxidized and reduced cases, respectively.Naively, one could just perform the required TI using ML surrogate models.As we will show in Results and Discussion, this yields only acceptable accuracy.Errors in ∆A resulting from inaccuracies in the trajectory and the energy predictions by the ML potential can be corrected by performing TI from the ML potential to the FP potential for both the oxidized and reduced states.We will adopt this strategy for the FP sl method.So this involves two calculations: TI from the oxidized to the reduced species using ML surrogate models via Eq.( 11) in Methods section, ML (Ox) → ML (Red), and then for each oxidation state, TI from MLFF to the FP sl Hamiltonian via Eq.( 13), FIG. 2. Schematic of ML-aided TI and TPT to compute the free energy difference ∆A.The notations ML, FP sl and FP nl mean the machine-learning force field, FP method with semi-local functional, FP method with non-local hybrid functional, respectively.See details of equations in the Methods section.ML(Ox)→FP sl (Ox) and ML(Red)→FP sl (Red).This twostep integration has three advantages as summarized below: • The integration ML (Ox) → ML (Red) using the MLFF takes into account most of the non-linear components of the integrand in the TI (see Figs. S8 and S9 in SI).Excellent statistical accuracy can be reached for this step.
• The MLFFs also provide well-equilibrated initial structures required for other calculational steps.
• The integrands in ML(Ox)→FP sl (Ox) and ML(Red)→FP sl (Red) are small and almost linear in the coupling parameter (see Fig. S10) owing to the accurate reproduction of the FP sl structures by the MLFF (see Results and Discussion).Hence, these demanding integrals (evaluation of FP sl calculation in every MD step) converge using a few tens of pico second MD simulations.
There is one final obstacle though: performing TI to a potential calculated by a hybrid functional that includes nonlocal exchange (FP nl ) is still exceedingly demanding and challenging.So in this specific case, as depicted in Fig. 2, we have decided to apply the ∆-machine learning (∆-ML) [52][53][54][55][56][57][58][59] which learns the difference ∆U between the FP sl potential and the FP nl potential.Due to the very smooth energy difference between the FP sl functional and the FP nl functional, it is possible to learn an extremely accurate ML representation of ∆U with just a few tens of FP nl calculations, with errors significantly smaller by an order of magnitude or more compared to those associated with MLFF models (see details in Figs.S2 to S4 and Fig. S17).In the current implementation, the TI integration has been replaced with thermodynamic perturbation theory (TPT), where β is the inverse temperature, and the symbol ∆U denotes the potential energy difference between the two end points.Although Eq. ( 9) is in principle exact, the potential energy difference might need to be evaluated for thousands or even many ten thousand configurations if the ensembles generated by the two potentials are too distinct.This implies the significantly expensive FP nl calculations.The ∆-ML scheme allows for the circumvention of this issue, enabling the reduction of the required FP nl calculations to merely tens.Thanks to the remarkable accuracy of the ∆-ML models, it is possible to obtain exceedingly accurate free energy differences between different FP methods without further correction (see Fig. S12 in SI).This is one of the key advances of the present work.The computational cost is ultimately only limited by generating sufficient configurations using the FP sl .Thus, the required compute time for direct TI using the FP nl method is reduced from 20 mio core hours to 16800 core hours for the FPMD simulations that generate configurations using the FP sl method (see details of the estimation in Section S2 in SI).

RESULTS AND DISCUSSION
We now detail our results, and will show that the adopted procedure yields statistically highly accurate results.The calculations were performed using VASP [60,61] and the projector augmented wave (PAW) method [62,63].For the ML force fields (MLFFs) the implementation detailed in previous publications is used [64][65][66].Similar to the pioneering ML approaches [54,67,68], the potential energy in our MLFF method is approximated as a summation of local energies [see Eq. ( 20)].The local energy is approximated as a weighted sum of kernel basis functions [see Eq. ( 21)].A Bayesian formulation allows to efficiently predict energies, forces and stress tensor components as well as their uncertainties.The predicted uncertainty enables the reliable sampling of the reference structures on-the-fly during the FPMD simulation.Details of the equations, parameters and training conditions are summarized in the Methods section and Section S1 in SI.As in the previous studies [64,[69][70][71], the MLFFs trained on a semi-local functional with dispersion corrections achieve root mean square errors (RMSEs) of 1-5 meV atom −1 and 0.04-0.11eV Å −1 for energies and forces (see error distributions in Figs.S1 to S4 in SI).The three ET reactions are examined in water by using a semi-local functional [72] with a dispersion correction [73,74] (RPBE+D3) and hybrid functionals [75,76] with and without a dispersion correction (PBE0 and PBE0+D3).Systematic comparisons of different functionals help us to study the effects of the exact exchange as well as dispersion corrections.As shown in Table 1 [see lines of PBE0 (0.25) and PBE0+D3 (0.25)] good agreement with experiment is achieved using the hybrid functional with 1/4 exact exchange, regardless of whether dispersion corrections are used or not.Water surface calculations: For RPBE+D3, the present MLFF provides a surface tension of 79±5 mN m −1 for the 128 molecular system and 84±5 mN m −1 for the 1024 molecular system at 298 K. Here, the surface tension was computed as [77] where x and y define the directions parallel to the macroscopic interface, z defines the direction normal to the interface, L z is the length of the unit cell in the z-direction, and p i j are the pressure tensor.The results are slightly larger than the value of 68±2 mN m −1 calculated by a neural network potential [77] and experimental value of 72 mN m −1 [78] while it is within the range (50-90 mN m −1 ) of previous MD results by FP [79] and classical force fields [80,81].Distributions of interfacial The small letters 'w/ ML' under FP sl mean that the ML result was corrected by the scheme shown in Fig. 2. The letters 'w/o ML' mean the result calculated by FP sl using Eqs.( 18) and ( 19) without MLFF.Here, FP nl means PBE0+D3 (0.25) result obtained by the scheme shown in Fig. 2. Experimental values are taken from Ref. [12].
water dipole moments for both, 128 and 1024 molecular systems, are shown in Fig. S5.They consistently indicate that the orientation of interfacial water molecules is bimodal as reported in previous MD studies employing the classical SPC/E force field [80].The distributions are also consistent with the results of sum frequency generation (SFG) analyses [82].
Metal water coordination: Figure 3 shows metal-oxygen radial distribution functions (RDFs) and running integration numbers (RINs) at the reduced and oxidized states calculated by the MLFF and FP sl methods.The MLFFs well reproduce RDFs and RINs of the FP sl method.Both methods show that the coordination number of Fe ions is 6 independent of the charge state.In contrast, the value for Cu changes from 5-6 in the oxidized state (Cu 2+ ) to 2-3 in the reduced state (Cu + ).The coordination number of Ag also changes from 5-6 in the oxidized state (Ag 2+ ) to 4-5 in the reduced state (Ag + ).These hydration structures agree with the ones reported in previous MD studies using FPMD methods [83][84][85] and empirical force fields [83].Although there are slight deviations in the Fe-O distance and shoulders for Cu-O and Ag-O in the RDFs likely related to the short FPMD simulation time and errors in the MLFFs, overall, our MLFFs reproduce the first-principles energies and structures of the hydrated metal cations with good accuracy.
Redox potentials: After verifying the size effect on the redox potentials obtained on the FP sl level with using unit cells containing 32, 64 and 96 water molecules (see U redox in Fig. S13), calculations were conducted on the bulk solutions containing 64 water molecules in the unit cell.The computed redox potentials are compared with the experimental ones in Fig. 4. All relevant data (⟨U 1 −U 0 ⟩ λ and ∆ φ ), as well as results of other functionals with error bars, are summarized in Sections S3, S4 and S5.The MLFFs trained on FP sl (RPBE+D3) (see ML in Fig. 4) lead to non-negligible deviations of 30-250 mV from the values of a full FP sl calculations without any MLFF (see FP sl w/o ML) depending on the training data size (see Section S9 in SI).The deviations can be corrected by two TI integrations [ML(Ox) → FP sl (Ox) and ML(Red) → FP sl (Red) in Fig. 2] as shown by FP sl w/ ML.However, the semi-local functional result in fairly large and non-systematic errors.For Ag the redox potential is underestimated, whereas for Cu it is significantly overestimated compared to experiment.
The errors can be significantly decreased to 0.11 V on average using hybrid functionals with one quarter exact exchange.As tabulated in Table 1, U redox for the Cu 2+ /Cu + couple decreases with increasing fraction of the exact exchange, whereas the redox potential for the Ag 2+ /Ag + couple increases with increasing the fraction.For Fe 3+ /Fe 2+ , the trend is not so obvious (first increase then slight decrease).Overall the present trends agree with the results obtained using semi-local and hybrid functionals as reported by Liu and co-workers [9].Finally, the effects of Grimme's dispersion correction is small for all redox couples.This implies that changes of the electronic properties (such as water valence band maximum and conduction band minimum) are most relevant, whereas all the functionals give a similar and good account of the solvation structure.It remains unclear, however, why one quarter of exact exchange results in balanced accuracy.The functional form of PBE0 was rationalized by the adiabatic connection from the uncorrelated exact exchange to the fully interacting energy, which is approximated by the PBE functional [75,86].Nonetheless, the ratio of exact exchange continues to be a parameter.One quarter of exact exchange is known to achieve balanced accuracy for the geometries, thermochemistry, and spectroscopic properties of molecules.However, as reported in previous studies [87], this functional underestimates the band gap of liquid water, even though it provides a more accurate prediction than the HSE06 functional.The mechanism behind this remains an open question.
Another key observation in this study lies in the relationship between the error of the ML surrogate model and the error in the redox potential.Our MLFF models achieve an RMSE of a few meV atom −1 for energy and tens of meV Å −1 for force.These accuracies can be considered standard level compared to ML models generated in past research [54,64,67,68,71,[88][89][90], yet they yield non-negligible deviations in the redox potential from the FP method.In comparison, ∆-ML models, which attained a RMSE substantially lower by more than an order of magnitude, markedly diminished the deviation in the redox potential to below 10 mV (refer to Fig. S12 in the SI).These results suggest that, in aiming for an accuracy of 10 mV in reproducing the redox potential of the FP method, an RMSE at least an order of magnitude smaller than that shown by standard MLFFs is required.Achieving this level of accuracy is highly challenging for MLFFs, even if they are trained on larger training datasets, as demonstrated in the previous study on liquid water [91].While the accuracy of emerging MLFFs continues to improve [92,93], there is always a risk that machine learning models may produce errors concerning the structure of extrapolation regions outside the training data.Even in a future where machine learning models have further advanced, our ML correction schemes will serve as a powerful method for quantifying errors and providing results from accurate FP calculations.
In summary, our approach enables efficient statistical sampling that is indispensable for accurate computations of the free energies of aqueous systems.The TI and TPT calculations allow to improve the accuracy from the ML model to the semi-local functional and from the semi-local functional to the hybrid functional step-by-step.Combining TPT and ∆machine learning are particularly promising, since this allows to obtain statistically highly accurate results even for expensive functionals in very little compute time.For instance, it is well conceivable that one could also use methods beyond density functional theory for the final step.Our final results reproduce the redox potentials of the three transition metal cations with excellent accuracy using a standard hybrid functional.The integration pathways chosen here are generalizable to a wide variety of electron transfer reactions.We believe that the scheme will pave the way to first-principles electrochemistry predicting the key property of redox reactions in energy conversion devices.

TI and TPT
The TI and TPT shown in Fig. 2 in the main text is conducted by using the equations listed below.
• ML(Ox) → ML(Red) • ML(Ox) → FP sl (Ox) and ML(Red) → FP sl (Red) • FP sl (Ox) → FP nl (Ox) and FP sl (Red) → FP nl (Red) The symbols U FP nl κ , U between the non-local and semi-local functionals.In Eq. ( 15), the second-order cumulant expansion is employed.The expansion is exact if the probability distribution of ∆U ∆ML κ is Gaussian (see derivation in Section S8 in SI).The condition is reasonably satisfiled as shown in Fig. S11.Preliminary TI and TPT simulations using the MLFFs also indicate that the TPT calculation reproduces TI results as shown in Section S6.
The free energy differences of the FP sl and FP nl methods are obtained as To validate the MLFF-aided computations of the free energy difference ∆A FP sl , the same property was also computed by the TI without using the ML method: The TI calculation in Eq. ( 2) can be decomposed into the two terms on the right-hand side of Eq. (7).The integrand in the first term nonlinearly varies along the coupling path (see Figs. S8 and S9 in SI) while the integrand in Eq. ( 8), which is relevant to the second term in Eq. ( 7), varies only slightly (see Table S4 in SI).To perform the integration of the first term in Eq. ( 7), the Simpson's rule with equidistant five points was used following the previous FP study by Blumberger and coworkers [45].For the integration in Eq. ( 8), the average of the O 1s levels in the fully reduced and oxidized states was used based on the trapezoidal rule.For each point, the ensemble average was taken over an 80-ps-NVT-ensemble MD simulation at 298 K after a 20 ps equilibration.Similar to the MLFF calculations, the Simpson's rule with equidistant five points was used for the TI calculation in Eq. (18).For each grid, the ensemble average was taken over a 20-ps-MD simulation starting from the final structure of the TI calculation using the MLFF at the same grid point.Each initial structure of the MD simulations was prepared by annealing the system from 400 K to 298 K by a 100-ps-NVT-ensemble MD simulation using the MLFF after annealing the same system from 1000 K to 400 K by a 1-ns-NVT ensemble MD simulation using the polymer consistent force field (PCFF) [94] implemented in a homemade MD program [95].Figures S8 and S9 in SI show the integrands of Eqs.(11) and (18), respectively, as functions of the coupling parameter λ .In the same figures, probability distributions of ∆U ML = U ML 1 −U ML 0 and ∆U FP sl = U Hence, the second cumulant expansion Eq. ( S12) is not applicable to the whole integration from the oxidized state to the reduced state.The TI calculations in Eq. ( 13) were conducted using the trapezoidal rule with three equidistant points.At each point, a 10-ps-NVT-ensemble MD simulation at 298 K was performed.The integrands shown in Fig. S10 are smaller than the ones shown in Figs.S8 and S9, respectively.They are also nearly proportional to the coupling parameter η.
In the TPT calculations using the ∆-ML models, the ensemble average in Eq. ( 15) was taken over 1400 configurations selected randomly from 70-ps-NVT-ensemble FPMD simulations using the FP sl method.Although these FPMD simulations are expensive, the overall computational time is still much smaller than full FP simulations.To ensure the applications of the second-order cumulant expansion, we show the probability distributions of the energy difference ∆U ∆ML κ in Fig. S11.The distribution are well fitted by Gaussian functions, indicating that Eq. ( 15) is reasonable approximations.
The MD simulations were performed using a Langevin thermostat [96].For efficient sampling, the mass of hydrogen and time step were set to 2 amu and 1 fs.

MLFF and ∆-ML
Similar to previous machine-learning approaches [67,68], the potential energy U of a structure with N a atoms in our MLFF method is approximated as a summation of local energies U i written as Following the Gaussian approximation potential pioneered by Bártok and co-workers [68], the local energy U i is approximated as a weighted sum of functions K ( The coefficients {w i B |i B = 1, ..., N B } are optimized to best reproduce the FP energies, forces, and stress tensor components as obtained by the FPMD simulations.The descriptor x i used in this study is a vector containing two and three body contributions [71]: Here, β (2) and β (3) (= 1 − β (2) ) are the weights on the two and three body descriptors, x i and x i , respectively.The vectors x (2) i and x (3) i collect the expansion coefficients of two and three body distribution functions with respect to the orthonormal radial and angular basis sets [64,71]: The two and three body distribution functions ρ (2) i and ρ (3) i are defined as: The function g is the smoothed δ -function, and f cut is a cutoff function that smoothly eliminates the contribution from atoms outside a given cutoff radius R cut .For χ nl and P l , normalized spherical Bessel functions χ nl = j l (q n r) and Legendre polynomials of order l are used in this work, respectively.For the kernel basis functions, the smooth overlap of atomic positions (SOAP) kernel [54] is employed The hat symbol xi denotes a normalized vector of x i .The normalization and exponentiation in Eq. ( 28) introduce nonlinear terms that mix two-and three-body contributions.The same formulation is used for the ∆-ML method.In the ∆-ML method, differences of potential energies and forces between two FP methods, semi-local and non-local functionals in this study, are used as the training data.
Parameter sets of the descriptors and kernel basis functions used in previous publications were employed in this study [64,66,71].The parameters are tabulated in Table S1.
Bulk solutions containing the redox species were modeled by systems as shown in Fig. 1.The number of water molecules was set to 32, 64, and 96.Three different model sizes were examined to clarify the system size effect.The sizes of the unit cells were set to obtain a water density of 0.99 g cm −3 .The size of the unit cell for the 32 water molecules is same as the one used in previous FPMD studies [9,15,45,84].For each of the reduced and oxidized states, MLFF and ∆-ML models were constructed.All MLFF models were generated on the fly during a 100-ps-NVT-MD simulation at 400 K by using the active-learning algorithm developed in our previous study [64].The temperature for the training runs was set to a value higher than the target one of 298 K for production runs, to ensure that training data and kernel basis functions were provided in a wider phase space.A Langevin thermostat [96] was used to control the temperature.Exchange-correlation interactions between electrons were described by the semi-local RPBE functional [72] with Grimme's D3 dispersion corrections [73,74].Probability distributions of the errors of the constructed MLFFs for energies and forces on test data are shown in Figs.S1 to S4 in SI.The RMSEs are similar to those of MLFFs used in previous studies [64,[69][70][71].
After examining the system size effect using the semi-local functional (see results in Fig. S13), ∆-ML models were constructed on systems containing 64 water molecules.Each ∆-ML model was trained on FP energies and forces of 40 structures selected randomly from a trajectory of a 20 ps NVTensemble FPMD simulation at 298K.The FPMD simulation was performed using the RPBE+D3 functional.Differences in energies and forces between the non-local and semi-local functionals for these 40 structures were used as training data.PBE0 [76] with and without the Grimme's D3 dispersion correction [73,74] was employed as the non-local functional because the functional is known to accurately predict properties of water [97].The fraction of the exact exchange was set to 0.25 and 0.50 to determine how this influences the redox potentials.Error distributions of the ∆-ML models on test structures are shown in Figs.S2 to S4.The RMSEs are one to two orders of magnitude smaller than the errors of the RPBE+D3 MLFFs.
The vacuum-water interface for the production run was modeled by a pure water slab without the redox species composed of 128 water molecules per unit cell.Following the previous study [49], a rectangular cell of 12.5 × 12.5 × 50 Å 3 was employed.Similar to the MLFFs for the bulk solution systems containing the redox species, the MLFF for the interface was also generated by using the active-learning scheme.The systems used for the training were a pure water bulk composed of 64 water molecules in a 12.4 × 12.4 × 12.4 Å 3 cubic cell and pure water slab composed of 64 water molecules in an 8.8 × 8.8 × 40.8 Å 3 rectangular cell.Training simulations for both the bulk and slab were performed by NVT-ensemble MD simulations at 300, 400, 600 and 800 K.As shown in Fig. S1, the constructed MLFF realizes small errors on test data taken from a 100-ps-MD simulation of a water slab composed of 128 water molecules at 298 K.
The annealing procesure used for the production runs explained in the previous subsection was also used to prepare for the initial structures for the training runs.All FP calculations were performed using VASP [60,61].A 2×2×2 k-point mesh was used for the bulk systems containing 32 water molecules.For other systems, Γ-point was used.Plane-wave cutoff energy was set to 520 eV.The PAW [62,63] distributed in VASP 5.4 was used in all FP calculations.The PAW atomic reference configuration was 1s 1 for H, 2s 2 2p 4 for O, 3d 7 4s 1 for Fe, and 4d 10 5s 1 for Ag.The comparison of two atomic configurations for Cu, specifically 3d 10 4p 1 and 3p 6 3d 10 4p 1 , was conducted to examine the impact of semi-core electron relaxations on the redox potential.Upon verification that these effects are minimal within the PAW framework in VASP, as detailed in Section S7 in SI, we employed the less computationally demanding 3d 10 4p 1 electronic configuration.The parameters for the MD simulations are same as the ones described in the previous subsection.

FIG. 1 .
FIG. 1. Aligning energy levels based on the O 1s level of water molecules: (a) aligning the redox level based on the O 1s level of water molecules far from the redox species in the bulk solution model under a periodic boundary condition, (b) aligning the O 1s level of water molecules at the middle of the slab based on the local potential at the middle of the vacuum layer in the slab model, and (c) schematic of the alignment.The figure inset in (b) shows the snapshot of the water slab and computed local potential profile across the water slab.The graphics showing bulk and interfacial models are made by VESTA [41].
tions and training runs, while retaining statistical accuracy, as demonstrated by the local potential profile shown in the inset of Fig. 1(b) [see details of the estimation of compute time in Section S2 in SI].

FIG. 3 .
FIG. 3. Metal-oxygen radial distribution functions (g X−O ) and running integration numbers (n X−O ) provided by 100 ps MLFF-MD and 10 ps FPMD simulations.Gray and black lines are for the reduced and oxidized states, respectively.Solid and dashed lines are results obtained by the MLFF and FP sl , respectively.Graphics in the insets show first solvation structures at the reduced and oxidized states.

FIG. 4 .
FIG.4.Computed and experimental redox potentials.ML means the results obtained by the MLFF trained on FP sl (RPBE+D3) without any correction.The small letters 'w/ ML' under FP sl mean that the ML result was corrected by the scheme shown in Fig.2.The letters 'w/o ML' mean the result calculated by FP sl using Eqs.(18) and (19) without MLFF.Here, FP nl means PBE0+D3 (0.25) result obtained by the scheme shown in Fig.2.Experimental values are taken from Ref.[12].
energies for the oxidized (κ = 0) and reduced (κ = 1) states calculated by the non-local functional, semi-local functional and MLFF trained on the semi-local functional, respectively.The symbol ∆U ∆ML κ denotes the potential energy difference calculated by the ∆-ML model trained on the potential energy difference U FP nl κ −U FP sl κ each λ are also shown.For all redox couples, the variance of the distribution varies with changing λ , and thus, the integrand is non-linear with respect to λ [see Eq. (S14) in SI].