Electron affinity of liquid water

Understanding redox and photochemical reactions in aqueous environments requires a precise knowledge of the ionization potential and electron affinity of liquid water. The former has been measured, but not the latter. We predict the electron affinity of liquid water and of its surface from first principles, coupling path-integral molecular dynamics with ab initio potentials, and many-body perturbation theory. Our results for the surface (0.8 eV) agree well with recent pump-probe spectroscopy measurements on amorphous ice. Those for the bulk (0.1–0.3 eV) differ from several estimates adopted in the literature, which we critically revisit. We show that the ionization potential of the bulk and surface are almost identical; instead their electron affinities differ substantially, with the conduction band edge of the surface much deeper in energy than that of the bulk. We also discuss the significant impact of nuclear quantum effects on the fundamental gap and band edges of the liquid.

In Eqs. (1) and (2), ε slab i and ε bulk i are energy levels referred to an arbitrary zero of energy; ∆V slab = V vacuum slab − V water slab is the difference between the plane-average electrostatic potential computed in the vacuum and water regions of the slab; V bulk is the average electrostatic potential computed in bulk water. The values of all these quantities for quantum and classical water determined and utilized in this work are summarized in Supplementary Table 1.
We explored two ways of computing the valence and conduction band edges: (i) we equated the position of the valence band maximum (VBM) and conduction band minimum (CBM) to the highest occupied and lowest unoccupied energy levels obtained in our calculations; (ii) we extrapolated the peaks in the EDOS to zero energy. To do so, we computed the second derivative of the EDOS as a function of the single particle energy to find inflection points, followed by calculation of the tangent line, as shown in Supplementary Figure 1.
As pointed out by Ambrosio et al. [12], different definitions of the band edges have different rate of convergence with respect to the system size. To understand size effects, we first explored the convergence of EDOS at the DFT level with respect to the number of k points included in the electronic structure calculations. Supplementary Figure 2 shows that the EDOS of occupied states is insensitive to the number of k points, while the shape of the EDOS of the empty states substantially depends on the number of k points included in the calculation, as first noted in Supplementary Ref. 13. However the energies of the highest occupied (HOMO) and lowest unoccupied (LUMO) orbitals change by less than 0.01 eV when using 64 k-points compared to the Γ-point only calculation.
A comparison of the positions of HOMO and LUMO energies computed with cells of different size, as well as the extrapolated valence band edge are summarized in Supplementary Table 2. The Table shows that the HOMO energy exhibits slow convergence with respect to the system size [12]. The extrapolated position of the 1b 1 band edge has instead a weaker system-size dependence.
On the other hand, extrapolation of the band edge for the conduction band minimum does not exhibit a satisfactory convergence trend with respect to the system size. This is shown in Supplementary Figure 3, which compares EDOS of unoccupied states of bulk water averaged along snapshots of the 64-and 256-molecule trajectories. Straightforward extrapolation of the conduction band minimum using the first, LUMO peak of the EDOS yields −2.54 eV for the 64-molecule system and −2.60 eV for the 256-molecule system, respectively. These values are substantially lower in energy ( by ∼ 0.3 eV) than those obtained using the energies of the LUMO states, −2.22 and −2.29 eV respectively (Supplementary Table 2). Since the converged dependence of the unoccupied EDOS as a function of energy is close to linear, extrapolations should be performed over a larger energy range, and not just by considering the first peak. We verified that when considering a larger energy range between −2.4 and −2.8 eV (see Supplementary Figure 2), the extrapolation of the EDOS computed at Γ and with 64 k-points yields a value for the band edge which is within 0.05 eV and 0.005 eV of that obtained considering the LUMO energy, respectively.
The HOMO energy computed for slabs behaves similarly to that computed for bulk water, as a function of the supercell size (see Supplementary Table 2). However we note that the fluctuations of the energy of the valence and conduction band edges as a function of the simulation time are much larger for slabs (0.4-0.6 eV) than in bulk samples (0.1-0.4 eV), leading to a larger error estimate of the valence and conduction band edges near the water surface, compared to the bulk. The LUMO energy of the slab converges slower with respect to the system size than that of the bulk, and it appears to be converged for the 216-water-molecule sample, within statistical error bars of 0.1 eV. Since performing hybrid DFT and G 0 W 0 calculations for 216-molecule slabs was prohibitively expensive, we computed band edges for the 108-molecule slabs at PBE, PBE0, RSH, and sc-hybrid levels of theory, and corrected the position of the CBM by the difference between the LUMO of the 108-and 216-molecule slabs (−0.22 eV), computed with the PBE functional. No size corrections were necessary for the extrapolated 1b 1 edge position for the 108-molecule supercell.
Based on the results reported in this section, we used a linearly extrapolated edge of the 1b 1 peak as a measure of the valence band maximum of both the bulk liquid and the surface. We used instead the LUMO energy as a definition of the conduction band minimum, corrected by −0.22 eV in the case of the water slab to take into account size effects, as discussed above.
We performed G 0 W 0 calculations starting with converged PBE and hybrid functional wavefunctions for classical and quantum bulk water. G 0 W 0 /PBE calculations were carried out for 128 snapshots along each trajectory, while the G 0 W 0 calculations starting with hybrid functionals were performed for 4 representative snapshots selected from the quantum and classical trajectories. Supplementary  Tables 3 and 4 show that the statistical error for LUMO estimated over a set of 4 snapshots is up to 0.1 eV, while the statistical error for the energy differences ∆G 0 W 0 is smaller, up to 0.03 eV. Therefore we computed the G 0 W 0 /hybrid conduction band edge by correcting the DFT LUMO energy with ∆G 0 W 0 values from Supplementary Tables 3 and 4. The statistical error bars on the resulting CBM were estimated as the sum of the error of the DFT energy and the error of the correction ∆G 0 W 0 . The valence band edge was always determined using the extrapolated 1b 1 peak of the G 0 W 0 /hybrid EDOS, without applying any corrections.
We also investigated the convergence of the G 0 W 0 energy of the LUMO orbital with respect to the number of eigenpotentials, N PDEP . We computed the LUMO energy of a single quantum water snapshot using G 0 W 0 /PBE as a function of the number of eigenpotentials (Supplementary Table 5). We then fit the resulting LUMO energies to the function a + b/N PDEP where a is the LUMO energy at infinite number of eigenpotentials. We find that the extrapolation procedure lowers the LUMO energy by 0.32 eV compared to calculations with 1600 eigenpotentials performed in this work. This correction of −0.32 eV determined at G 0 W 0 /PBE level of theory was applied to all G 0 W 0 LUMO energies computed in this work.
G 0 W 0 calculations for the water/vacuum interface were prohibitive, from a computational standpoint, if conduced at the same statistical level as for the bulk water. Hence, we obtained G 0 W 0 energies for the slab in the same way as we did for the G 0 W 0 /hybrid results-DFT energies computed for the slab, plus the ∆G 0 W 0 correction from Supplementary Table 3 and 4, extrapolated to the infinite number of eigenpotentials. The analysis provided above showed that this scheme adds a modest additional statistical error of ∼ 0.03 eV. To check that this protocol still applies to water in contact with vacuum, we performed G 0 W 0 /PBE calculations for one snapshot of the quantum water/vacuum interface at several numbers of eigenpotentials, and then extrapolated the LUMO energy to the infinite N PDEP . The results, shown in Supplementary Table 6, indicate that ε LUMO computed as a correction to the DFT energy differs from the extrapolated G 0 W 0 value by 0.01 eV, thus justifying our computational procedure.

Supplementary Note
This Section summarizes the values of the electronic properties of classical and quantum bulk water (Supplementary Tables 7 and 8), as well as water/vacuum interface (Supplementary Tables 9 and 10).
Nuclear quantum effects (NQE) on water structure have been discussed in Supplementary Ref. 14. Here we report on NQE influence on the electronic structure. Supplementary Figure 4 compares the experimental photoelectron spectrum to the density of states, computed following the protocol and the best level of theory (G 0 W 0 quasiparticle energies obtained using RSH wavefunctions) presented in Supplementary Refs. 9 and 15. The Supplementary Figure 4 shows that NQE have a softening effect on the shape of the spectrum. The width of the bands increases and their magnitude decreases, consistent with the recent study of Supplementary Ref. 16. Overall, the spectrum computed over the PIMD trajectory appears to be closer to experimental results, although photoionization cross-sections would need to be taken into account for a full detailed comparison [9].
Supplementary Figures 5 and 6 compare electronic properties (effective molecular polarizabilities, dipole moments, distribution of Wannier centers) of water molecules in bulk water (classical and quantum) and at the interface of quantum water with vacuum. ) portions of the slab simulated with (Quantum water) and without (Classical water) nuclear quantum effects, computed with different density functionals (PBE [17], PBE0 [18], RSH [19], and sc-hybrid [20]). The sum of ∆V slab = V vacuum slab − V water slab and V bulk (0.12 eV) used in Eq. (2) is also reported. G 0 W 0 /DFT eigenvalues were aligned with vacuum using respective DFT values. All values are in eV.

Quantum water
Classical water     Table 3: Density-functional theory (DFT) and G 0 W 0 /DFT lowest unoccupied molecular orbital (LUMO) energies averaged for select snapshots extracted from the classical MB-pol trajectory for bulk water. ∆G 0 W 0 refers to differences between the G 0 W 0 quasiparticle energies and the respective DFT eigenvalue, and represents a correction provided by the G 0 W 0 method to the LUMO energy. "Average (4 samples)" denotes the average computed using 4 snapshots taken from each trajectory, while "Average (trajectory)" refers to the average computed along the entire trajectory (128 snapshots). The latter is available for all DFT methods and for G 0 W 0 /PBE, but is not available for G 0 W 0 /hybrid calculations. All errors are computed as three standard deviations of the mean on the respective quantity. The energies (in eV) are given on an absolute scale, referred to vacuum using Eq. (2) Table 3 but for the quantum MB-pol trajectory for bulk water. ∆G 0 W 0 refers to differences between the G 0 W 0 quasiparticle energies and the respective DFT eigenvalue, and represents a correction provided by the G 0 W 0 method. "Average (4 samples)" denotes the average computed using 4 snapshots taken from each trajectory, while "Average (trajectory)" refers to the average computed along the entire trajectory (128 snapshots). The latter is available for all DFT methods and for G 0 W 0 /PBE, but not for G 0 W 0 /hybrid calculations. All errors were computed as three standard deviations of the mean on the respective quantity. The energies (in eV) are given on an absolute scale, referred to vacuum using Eq. (2) Table 7: Positions of peaks in the electronic densities of states (EDOS) computed along the classical MB-pol trajectory of bulk liquid water using a series of density functionals (PBE [17], PBE0 [18], RSH [19], and sc-hybrid [20]) and many-body perturbation theory (G 0 W 0 ). The mean absolute error (MAE) of peak positions reflects the average deviation of computed peaks from the experimental ones. The valence band maximum (VBM) was determined by the linear extrapolation of the 1b 1 peak of water to zero. The conduction band minimum (CBM) is the LUMO energy, and in the case of G 0 W 0 calculations were extrapolated to infinite number of eigenpotentials. The energies were aligned with respect to vacuum using Eq. (2)  Supplementary Table 8: Positions of peaks in the electronic densities of states (EDOS) computed along the quantum MB-pol trajectory of bulk liquid water using a series of density functionals (PBE [17], PBE0 [18], RSH [19], and sc-hybrid [20]) and many-body perturbation theory (G 0 W 0 ). The mean absolute error (MAE) of peak positions reflects the average deviation of computed peaks from the experimental ones. The valence band maximum (VBM) was determined by the linear extrapolation of the 1b 1 peak of water to zero. The conduction band minimum (CBM) is the LUMO energy, and in the case of G 0 W 0 calculations extrapolated to infinite number of eigenpotentials. The energies were aligned with respect to vacuum using Eq. (2) and the values from Supplementary Supplementary Table 9: Positions of peaks in the electronic densities of states (EDOS) computed along the classical MB-pol trajectory of water/vacuum interface using a series of density functionals (PBE [17], PBE0 [18], RSH [19], and sc-hybrid [20]) and many-body perturbation theory (G 0 W 0 ). The mean absolute error   Figure 1: Graphical representation of our calculation of the valence band edge by extrapolating the 1b 1 peak of the electronic density of states (EDOS) as a function of the single-particle energy (ε). EDOS(ε) was obtained using the PBE functional along the quantum MB-pol trajectory of bulk water in the 64-molecule supercell.