As a test case we computed here the simple oxides: NiO and CoO, both containing TM-cation in + 2 oxidation state. The Hubbard U parameters derived for TMs and the resulting band gaps in these two compounds are provided in Table 1. CoO and NiO are wide band gap charge transfer Mott insulators (Zhang et al., 2021). For such materials, the band gap should be equivalent to half of the Hubbard U parameter (expected shift of the unoccupied d states). The Hubbard U parameters computed here for these compounds are large and indicate wide band gaps (Table 1). The band gaps computed with the DFT + U method are consistent with the measurements. This is, however, somewhat misleading, as it results from incorrect occupations of d states. Because the atomic orbitals are applied as projectors, these result in significant, fractional occupations of empty orbitals (Table 2) and overestimated total occupancy of the d states. This shortcoming is corrected with the Wannier-functions as projectors. The resulting occupation matrix and number of d electrons are also reported in Table 2. These reflect the expected values, which has also significant impact on the computed band gaps. In set of previous studies we obtained similar improvement for occupations of d and f orbitals of various other elements [e.g., Murphy et al. (2021); Kowalski et al. (2021)]. The band gaps computed with the DFT + U(WF) method are larger than the ones predicted with the standard DFT + U approach. We note that for this case, the hybrid functionals such as HSE06 also predict correct band gaps of the considered oxide materials (Seo et al., 2015). This is opposite to the case of lithiated transition metal oxides, which show significant level of d elements delocalization and overestimation of the band gaps of these materials by the hybrid functionals approaches (Seo et al., 2015).

Metal phases consisting of transition metal elements play a key role in the electrochemical devices. The strong correlation character of d electrons these systems contain has been illustrated by the computation of large Hubbard U parameters in a series of computational studies (Schnell et al., 2002; Nakamura et al., 2006; Şaşıoğlu et al., 2011; Tesch and Kowalski, 2022). This would imply the importance of Hubbard model corrections also for these systems. On the other hand, it has been demonstrated that the standard DFT method results in very good prediction of metal properties, including good description of the x-ray photoelectron spectroscopy (XPS) spectra by the resulting density of states (DOSes) (Hofmann et al., 2012). This apparent contradiction could be explained by the delocalized character of d electrons in metals and inapplicability of the fully localized limit (FLL) of standard DFT to metals. An “around mean field” (AMF) version of DFT has been proposed instead for computation of metals (Czyżyk and Sawatzky, 1994; Petukhov et al., 2003; Himmetoglu et al., 2014; Ryee and Han, 2018). In this approach, the occupations are forced to reproduce the mean values, instead of the preference of fully occupied or unoccupied orbitals, as in the standard DFT + U approach. This causes no overall change in the DOSes when the AMF version of the DFT + U method is applied. Nevertheless, it has been also proposed that the correct computation of metals should be done with a combination of the FLL and AMF methods (Petukhov et al., 2003). Having wide popularity of the FLL version of the DFT + U approach, it is thus of importance to understand the performance of this method when applied to metallic systems.

Tesch and Kowalski (2022) derived the Hubbard U parameters for all the 3d, 4d and 5d transition metals applying the linear response method (Cococcioni and de Gironcoli, 2005). The values (up to 11 eV) and trends of increasing the U parameter values with increasing number of valence electrons, observed in previous studies (Şaşıoğlu et al., 2011), have been reproduced. These show that the 3d metals exhibit the largest Hubbard U parameter values, followed by 4d and 5d metals. Having these results, Tesch and Kowalski (2022) performed a thoughtful analysis of the performance of the DFT + U method for metals. When computing transition metals with the DFT + U approach, they observed a shift of the d-band to lower energies, comparing to the experimental XPS spectra. Interestingly, they found a very similar behavior when hybrid functionals (like PBE0) were applied. As mentioned, this behavior could be prevented when the AMF approach is applied. However, Tesch and Kowalski (2022) also demonstrated the importance of the projectors of d orbitals for the overall results. When the Wannier functions-based projectors of the occupancy of d states were applied, the spurious shift of the d-band was prevented. A correct description of metal d-bands is extremely important, since they are often used as descriptors for catalytic activity, like in the famous d-band model (Hammer and Norskov, 1995; Ruban et al., 1997).

Layered oxides are the most widely used electrode materials in rechargeable lithium-ion batteries. Among them, LiCoO₂ (LCO) is one of the successfully commercialized cathodes (Blomgren, 2016). However, the increased demand of cobalt resulted in a high cost of such materials, triggering a widespread research effort to find alternative cobalt-free cathode materials, applicable especially for automotive applications (Olivetti et al., 2017). In the past decades, LiNiO₂ (LNO) has been considered as an alternative cathode material for rechargeable lithium batteries. This is because of its high availability and associated low cost, as well as safety and energy-density when compared to the commercial cobalt counterpart (Kalyani and Kalaiselvi, 2005; Mukai et al., 2010). If accurately described by first principles methods, a valuable, molecular-level insight into the properties that determine the cycle stability and charge capacity, when these materials are used in electrochemical compounds, could be obtained (Chakraborty et al., 2018). There exist several first-principle studies on LNO [e.g., Yoshida et al. (2019); Zhu et al. (2011); Chen et al. (2011); Mock et al. (2021))]. However, DFT fails to correctly predict the electronic structure of these systems: the band gap is severely underestimated and oxidation states of cations are incorrectly described (Pavarini et al.,2012).

The DFT + U method improves the materials description in the outlined aspects (Seo et al., 2015; Kowalski et al., 2021). However, even this correction can fail for the TM elements with more than five d-electrons, such as iron or nickel (Kowalski et al., 2021). The description of d-electrons in nickel can be further improved by representing the orbitals with strongly correlated electrons with Wannier-type functions (Marzari et al., 2012; Gu et al., 2020; Kowalski et al., 2021; Murphy et al., 2021).

Here, we computed the Hubbard U parameter for Ni in the layered rhombohedral ( R 3 ̄ m ) structure and obtained the value of 7.96 eV. It is well known that LCO has a stronger covalency compared to CoO, which induces less charge localization on the TM (van Elp et al., 1991; Galakhov et al., 1996). Seo et al. (2015) studied the LCO and LNO compounds with first-principle calculation using the HSE06 hybrid functional and showed that the optimal exact exchange admixing parameter α for LCO and LNO is substantially lower than the default one (α = 0.25), namely 0.17 for LCO and 0.18 for LNO. We assume here that the strong covalency in LCO and LNO will also affect the derived Hubbard U parameter and that the more realistic value should be these rescaled to U = 0.18/0.25 ⋅ 7.96 eV = 5.73 eV. As shown in Figure 1, comparing to the derived U parameter, the rescaled parameter describes more accurately the density of states of LiNiO₂ (Chakraborty et al., 2018).

In order to track the effect of delithiation, series of Li _x NiO₂ compounds was computed with different Li concentration, namely with x = 1.00, 0.75, 0.50, 0.25 and 0.00. The geometries were relaxed and the Hubbard U parameter was computed and rescaled for all the considered structures. Beyond the expected oxidation of Ni³⁺ into Ni⁴⁺ (three atoms at a time with increasing Li content), no disproportionation of Ni³⁺ into Ni²⁺ and Ni⁴⁺ species was observed during delithiation process, after wannierization scheme was applied. This suggests that the instability of Ni⁴⁺ reported in literature is due to other chemical forces during the synthesis of the material rather than to an inherent instability of the final structure itself (Li et al., 2022). The DFT + U(WF) method resulted in overall more stable structures, lower in energy by at least 10 kJ/mol, and the correct occupations of d orbitals, as in the case of NiO and CoO.

The intercalation voltage was estimated using standard approach (Aydinol et al., 1997; Dixit et al., 2016): V = − E Li x + d x NiO 2 − E Li x NiO 2 + d x ⋅ E Li bcc d x , (2) with dx = 0.25 as a step change of Li content, E (Li _x NiO ₂) the total energy per formula unit at Li content x, and E (Li_bcc) the total energy per formula unit of Li in the bulk material. As shown in Figure 2, the average calculated value of 4.17 eV is in good agreement with the measured value (Ohzuku et al., 1993).

In order to test the performance of the applied method for predicting the formation enthalpy of LNO, we also computed the Li₂O and NiO compounds, and the molecular O₂. Assuming the following LNO formation reaction: 0.5 L i 2 O + N i O + 0.25 O 2 → L i N i O 2 , (3) we estimate the formation energy of LNO from oxide to be 65.12, kJ/mol. When the entropy of gas phase (oxygen) is considered (S = 205 J/mol/K at 293 K), the resulting free energy of formation is 50.31, kJ/mol, is in good agreement with the measured value of 53.35, kJ/mol (Wang and Navrotsky, 2004).

During the delithiation process, as a result of oxidation of Ni atoms, removal of Li causes a linear increase of the Hubbard parameter U (Figure 3). Such a dependence of the Hubbard U parameter on the oxidation state is consistent with our previous studies [e.g., Beridze et al. (2016); Beridze and Kowalski (2014); Kowalski et al. (2021)]. The change in the oxidation state of Ni from + 3 to + 4 results in change in the DOS and the band gap (Figure 4), and in shortening of the Ni-Ni, Ni-O and O-O bond lengths, and increase of Ni-Li distances, as shown in Table 3. The cell parameters of the optimized supercell show excellent agreement with the crystallographic data.

It is worth noting that, while the applied Wannierization scheme contributes to the correct description of the oxidation states of Ni cations, it also replicates the small experimental band gap of 0.4 eV in LiNiO₂ (Laubach et al., 2009; Shishkin and Sato, 2016), as indicated in Figure 4. On the other hand, in most of the computational studies, perfect-layered LNO was computed as half-metal even with the computationally intensive hybrid functionals (Laubach et al., 2009; Shishkin and Sato, 2016; Chakraborty et al., 2018). In our studies, the accurate computation of occupation of d orbitals was prioritized, as the changes in the oxidation state of Ni are crucial to correctly describe the delithiation process.

Similar computations were performed for LiCoO₂ in the layered rhombohedral ( R 3 ̄ m ) structure. We derived a Hubbard U parameter of 6.8 eV, which for the same reason as in the case of LNO, we assume to be overestimated. The computed Hubbard U parameter value was rescaled, taking the optimal mixing parameter α for the exact Hartree-Fock exchange as a reference (Seo et al., 2015). This procedure gave a value of 4.6 eV. The computed band gaps with different U parameter values and with the DFT + U and DFT + U(WF) methods are shown in Table 4. A good agreement with the experimental value of 2.7 eV (Galakhov et al., 1996) is obtained with the rescaled value of Hubbard U parameter.

Applying the same computational approach as for LNO, we also computed the formation energy of LCO from oxides, following the reaction: 0.5 L i 2 O + C o O + 0.25 O 2 → L i C o O 2 . (4)

We obtained the value of 128.50 kJ/mol. Though less accurate than the one obtained for LNO, the computed LCO formation enthalpy from oxides shows reasonable agreement with the experimental value of 140.18 kJ/mol (Wang and Navrotsky, 2004; Takahashi et al., 2007).

Spent nuclear fuel contains significant amount of fission products or minor actinides (Np, Am, Cm), elements which can alter the material structure and performance (Bosbach et al., 2020). Doping of a phase with an element of other oxidation state can trigger a phase transition and stabilization of the new phase. This is for instance the case of yttria-stabilized zirconia, where doping of zirconium dioxide phase with ∼8% of tri-valent element like yttrium, stabilizes the cubic phase of that compound, which is one of the fastest known ionic conductors (Kowalski et al., 2021). Here we investigate the phase transition upon doping U₃O₈ phase (realized as a phase with the space group C2mm) with small amount of Am. This system has been investigated experimentally by Caisso et al. (2016). They found stabilization of a phase with the space group P62m phase upon doping with ∼ 10 % of Am. Here we confirm this finding with the aid of atomistic simulations.

First, we computed the Hubbard U parameters for U and Am in both structures. The resulting values are reported in Table 5. These are consistent with our previous studies (Beridze and Kowalski, 2014; Beridze et al., 2016). We notice, however, that the values derived for U species in U₃ O ₈ are larger than the ones obtained in Kvashnina et al. (2018). This is because in those studies another structure of U₃ O ₈ was investigated.

The computed formation enthalpy difference between the phases with space groups P62m and C2mm as a function of Am content are shown in Figure 5. Following studies of Caisso et al. (2016) we assumed that Am is incorporated as +3 species on U +5 site and the charge is compensated by the conversion of the two U +5 species to the oxidation state of +6. We predict the phase transition at very low content of Am of ∼ 1.2 % . It is thus not surprising that Caisso et al. (2016) measured the P62m phase at higher Am content of ∼ 10 % . The obtained volumes of pure and Am doped U₃O₈ are also well consistent with the measurements of Caisso et al. (2016) (see Table 6). The computed values are within 1% of those measured, and the computed increase in volume upon doping with Am is also well consistent with the measured values.

To show the importance of solvation for computation of surface electrochemical pathways, as a test case we investigated the effect of explicit and implicit solvation on the CO₂ reduction reaction towards CO and HCOOH on the Pb(100) surface. We considered the following reaction pathways for the formation of CO and HCOOH species: C O 2 g + H + + e − → H C O O * (5) H C O O * + H + + e − → H C O O H g (6) C O 2 g + H + + e − → C O O H * (7) C O O H * + H + + e − → C O g + H 2 O (8)

We assumed that the CO₂ forms either a HCOO* intermediate, which reacts further to formic acid (HCOOH), or COOH* which further reacts to CO and water.

Figure 6 shows the energy diagram for both pathways computed in three ways: (M1) assuming no solvation, (M2) in the presence of explicit solvation and (M3) by applying the continuum solvation model. In the approach M2, the aqueous phase was modeled by explicitly adding 12 water molecules. The case M3 was computed using the implicit solvation model implemented in VASPsol method (Mathew et al., 2019). The adsorption energies for HCOO* and COOH* intermediates computed with M1 are -0.45 eV and 0.73 eV, respectively. In the presence of implicit solvation, the HCOO* and COOH* adsorption energy increases to −0.61 eV and 0.52 eV. These indicated bonded HCOO* and unbonded COOH*. Applying method M2, the adsorption energy of the two species increases by 0.27 eV and 1.00 eV, respectively. These results show the importance of including solvent effects, when computing chemical pathways on metal surfaces.

The method M3 is a fast way to include solvent effects in the DFT calculations. In the considered cases, it results in increase in the adsorption energies, but significantly smaller than the explicit solvation schemes. The adsorption energies are changing by 0.1 eV only.

In Figure 6, we also show the results of hybrid solvation approach, where on top of the 12 explicit water molecules we have applied the implicit solvation, to preserve continuity of the solvent medium. While for the HCOO* intermediate, the effect of implicit solvation leads to an insignificant change, COOH* is stabilized by more than 0.2 eV. This shows that for the adsorbed species, the effect of first hydration shell dominates the solvent effect on the bonding energy between the adsorbate and the metal surface. In case of the weak bonding, the long-range interaction between the species and the solvent plays an important role.

Simulation of electrochemical solid/liquid interfaces at an applied electrode potential is a challenge in computational electrochemistry (Schwarz and Sundararaman, 2020). It requires an approach that can consistently describe charging of the electrode as well as effects of an electrolyte solution. Both conditions are realized by the DFT-based effective screening medium reference interaction site method (ESM-RISM) (Nishihara and Otani, 2017). The RISM (Chandler and Andersen, 1972; Hansen and McDonald, 2013) is an implicit solvent model that relies on the classical theory of liquids. It computes in a statistical way correlations between electrolyte species and thus takes into account the electrolyte structure. Electrolyte–electrolyte as well as electrode–electrolyte interactions are described by the parameterized interparticle interaction potentials. The method requires a choice of classical water model for description of the aqueous electrolyte solution, a set of Lennard–Jones parameters for the electrode–electrolyte interaction and the partitioning of the system into quantum-mechanically and classically treated parts (i.e., treatment of near-surface water layers at the level of DFT) (Nishihara and Otani, 2017; Tesch et al.,2021).

We applied this method in our previous study for computation of the partially oxidized Pt(111)/electrolyte interface (Tesch et al., 2021). Taking into account the potential-dependent oxygen coverage of Pt(111) surface, we were able to accurately describe the chemisorption and charging state of the interface, correctly reproducing the chemisorption-induced non-monotonic charging relation of this specific interface [see Figure 8 of Tesch et al. (2021)]. Here, we applied the ESM-RISM approach to model the interface between the Au(111) electrode and an 0.1 M aqueous HCl solution. Figure 7 shows the obtained interface structure of the electrolyte solution as represented by the density distribution functions of water and electrolyte ions species. The near-surface structure of water solvent shows alternating shells of water, with the main peaks positions that agree well with the AIMD simulations (Goldsmith et al., 2021). The effects of the applied potential and the resulting surface charge on excess or depletion of ions near the interface, as expected from the electrostatic arguments, are clearly seen. This indicates that the electrode charge is correctly balanced by the arrangement of electrolyte ions. The relation between the surface charge and the applied potential is linear (see Figure 8). This is expected because the Au(111) surface under these conditions (the considered potential range) is characterized by the formation of double layer and is not covered by adsorbed species, as indicated by the cyclic voltammetry (Hamelin, 1996). Other computational studies have derived very similar charging relation (Letchworth-Weaver and Arias (2012); Goldsmith et al. (2021), see in Figure 8). The trends are consistent with the experimental findings of linearly decreasing (i.e. opposite) charge at the outer Helmholtz plane (Saha and Zenyuk, 2021). The double layer capacitance, as derived from the slope of linear fit, is 23.4 μF/cm² and compares well to the experimental (Vasiljevic et al., 2004; Garlyyev et al., 2018) and other computational (Kastlunger et al., 2018; Hörmann et al., 2019; Goldsmith et al., 2021) results.

The reliable computation of surface chemical reactions at electrochemical conditions requires correct estimates of thermodynamic parameters. In particular, entropy effects are usually omitted in atomistic modeling studies [e.g., Klinkova et al. (2016))]. However, most of the electrocatalytic reactions on a solid catalyst surface happen in the presence of a dense electrolyte. The entropic contributions from the solvation effects can have a non-negligible effect on the reactions free energy. There are ways to compute the entropy of species in a solvent phase. These are based either on the theoretical consideration [e.g., Garza (2019)] or simulation based approaches, such as thermodynamic integration (Baranyai, 2018) or analysis of the velocity autocorrelation function (Pascal et al., 2011). Garza (2019) proposed a simplified approach to the computation of entropy of liquid phases and obtained agreements with the measured values within 12 J/mol/K. Regarding entropy at the interfaces, Jung et al. (2021) applied the two-phase thermodynamics (2PT) approach (Lin et al., 2010; 2003; Pascal et al., 2011) and computed entropy of aqueous phase as a function of distance from the Pt(111) interface. Here we applied the 2PT approach to evaluate different entropy contributions for HCOOH molecules submerged in aqueous phase, in the bulk and at the Pb(100) and Ag (100) surfaces. The so estimated total entropy is given in Table 7.

The computed molar entropy of HCOOH in bulk water is in good agreement with the experimental values. For HCOOH on Pb(100) surface and Ag (100) surface, we observe a reduction in entropy by ∼ 30 J / m o l / K . However, we observe a slightly different magnitude in entropy reduction at the Ag (100) and Pb(100) surfaces, with the Ag (100) surface causing larger reduction in entropy. This can be explained by more hydrophilic character of Ag (100) surface than that of Pb(100), and thus easier movement of HCOOH when adsorbed in the latter case. The reduction of ∼ 30 J / m o l / K , moving a species from the bulk water to the metal surface, contributes to the reaction free energy by ∼ 0.1 e V , which is a non-negligible effect. This indicates necessity of accounting for the solvation entropy effects when computing the electrochemical reactions under realistic reaction conditions.