Solvation eﬀects on the band edge positions of photocatalysts from first principles †

The band edge positions of photocatalysts relative to the redox potentials of water play an important role in determining the eﬃciency of photoelectrochemical cells. These band positions depend on the structure of the solid–liquid interface, but direct ab initio molecular dynamics calculations of these interfaces, while expected to be accurate, are too computationally demanding for high-throughput materials screening. Thus rapid theoretical screening of new photocatalyst materials requires simplified continuum solvation models that are suitable for treating solid–liquid interfaces. In this paper, we evaluate the accuracy of the recently developed CANDLE and SaLSA continuum solvation models for predicting solvation effects on the band positions of several well-studied surfaces [Si(111), TiO 2 (110), IrO 2 (110) and WO 3 (001)] in water. We find that the solvation effects vary considerably, ranging from o 0.5 eV for hydrophobic surfaces, 0.5–1 eV for many hydrophilic oxide surfaces, to B 2 eV for oxygen-deficient surfaces. The solvation model predictions are in excellent agreement (within B 0.1 eV) with ab initio molecular dynamics results where available, and in good agreement (within B 0.2–0.3 eV) with experimental measurements. We also predict the energetics for surface oxygen vacancies and their effect on the band positions of the hydrated WO 3 (001) surface, leading to an explanation for why the solvation shift observed experimentally is substantially larger than predicted for the ideal surface. Based on the correlation between solvation shift and the type of surface and solvent, we suggest approaches to engineer the band positions of surfaces in aqueous and non-aqueous solutions.


Introduction
Artificial photosynthesis, the reduction of water to H 2 or CO 2 to carbon-based fuels using the energy from sunlight in a photoelectrochemical cell (PEC), provides a promising path towards clean renewable energy while reducing CO 2 emissions. 1,2In addition to high catalytic activity for the reaction of interest, the reliability and efficiency of photocatalysts in a PEC depend critically on the alignment of their band edges, for example, with the redox potentials of water for water splitting solar cells. 3][6][7][8] However, most theoretical predictions are based on electronic structure calculations in vacuum 5,9,10 which neglect the substantial solvation effects on the electronic states of the photocatalysts.Depending on the interaction between water and the photocatalyst surfaces, the shift of electronic states due to solvation can be as large as 1-2 eV. 11Consequently, neglecting solvation effects can qualitatively change predictions for certain materials.
In fact, ab initio molecular dynamics (AIMD) calculations of explicit liquid water/photocatalyst interfaces have been validated to provide accurate estimates of the band positions. 11Unfortunately, these calculations are highly computationally expensive, requiring long simulation times to obtain statistically meaningful results.Evaluating reaction mechanisms at the solid/liquid interface introduces additional challenges and computational cost in such AIMD calculations, making it impractical today as a systematic tool for rapid theoretical screening of new photocatalysts.
Continuum solvation models directly abstract the thermodynamically-averaged effect of the liquid in an electronic structure calculation of the solute alone, and potentially provide an accurate but computationally affordable way of including solvent effects in the rapid theoretical screening of photocatalyst materials.However, these models are parameterized primarily to describe the solvated free energies of organic molecules, [12][13][14] and do not extrapolate well to strongly ionic surfaces. 15The absence of unambiguous thermodynamic data for transition metal oxide surfaces, analogous to the solvation energies available for molecules, precludes the direct parameterization of empirical solvation models for these systems.Consequently, band offsets in solvated ionic surfaces have not yet been validated, requiring the development of reliable non-empirical methods.
Here, we evaluate the accuracy of the recently-developed non-empirical SaLSA solvation model, 16 and the empirical CANDLE solvation model 17 derived from SaLSA, for the prediction of band positions in solution.We find good agreement for both models with available experimental measurements and with AIMD simulations of well-studied photocatalysts, such as Si, TiO 2 and WO 3 .We also explore the effect of surface oxygen vacancies on the solvated band edge positions of WO 3 , which is known to exhibit highly oxygen-deficient surfaces. 18,19Our study compares solvation effects on various surfaces, including hydrophobic, hydrophilic, ionic and non-ionic in order to provide a general understanding of the magnitude of solvation shifts for different types of surfaces and to guide discovery of new materials.
We start with a brief description of solvation models and our protocol for calculating the solvation shift on the electronic states in the Methods section.The Results section begins with a comparison of solvation models with AIMD predictions for functionalized Si surfaces, and presents comparisons between theoretical predictions and experimental measurements on the solvation shift of the band positions of the stable TiO 2 , IrO 2 and WO 3 surfaces.We also discuss the effect of oxygen vacancies on the WO 3 band positions and conclude the Results section with the trends of solvation shifts for surfaces in different solvents.The final two sections summarize these results, point out the general trends in solvation shifts that are relevant for ab initio photocatalyst design, and suggest future approaches to engineer band positions in solution.

Methods
We performed first principles calculations of Si, TiO 2 , IrO 2 and WO 3 slabs in vacuum and solution using the open-source plane-wave density functional theory software, JDFTx. 20This software is designed specifically for electronic structure calculations of systems in solution within the framework of joint density functional theory, 21 and enables the rapid development of solvation models (such as SaLSA 16 and CANDLE 17 ) using the algebraic formulations for electronic and classical density functional theories. 22,23We used the PBE 24 generalized gradient approximation (GGA) along with DFT-D2 pair-potential dispersion corrections, 25 and GBRV ultrasoft pseudopotentials. 26See the ESI, † for further computational details of the DFT calculations and detailed structures of all calculations presented below.
Density functional predictions for the band gaps and absolute positions of the valence band maximum (VBM) and conduction band minimum (CBM) are inaccurate in comparison to experiment, 27 especially for semi-local exchange-correlation functionals such as the PBE GGA.Consequently, we calculate the shift of the electronic states due to solvation which is mainly an electrostatic effect; and compare it to the experimental solvation shift deduced from vacuum and electrochemical measurements, as detailed below.Specifically, we calculate the theoretical solvation shift by subtracting the planar-averaged electrostatic potential obtained from solvated and vacuum calculations.
This electrostatic potential difference is insensitive to the choice of exchange-correlation approximation with differences only B0.1 eV between PBE, hybrid functionals and even many-body perturbation theory (GW approximation) calculations, [28][29][30][31] and can be reliably compared to experiment.
We constructed inversion-symmetric slab models for each material with sufficient number of layers (4 to 6) to converge the electrostatic potential in vacuum as well as solution, and optimized all ionic positions self-consistently both in vacuum and solvent.Inversion symmetry ensures identical top and bottom surface configurations of the slab with zero net dipole in the unit cell in all directions.We employ truncated Coulomb potentials [32][33][34] to eliminate spurious interactions between periodic images and accelerate the convergence with respect to the thickness of the vacuum or solvent region.The electrostatic potential decays to zero away from the slab in all cases, because of the truncation in vacuum calculations and the Debye screening due to the electrolyte in the solvated ones. 35This establishes a well-defined vacuum reference (zero potential at infinity) in solvent as well as vacuum, thereby enabling reliable electrostatic potential difference calculations.
Solvation models replace the thermodynamically-averaged effect of the liquid by the electrostatic response of a continuum dielectric cavity along with corrections for cavity formation and dispersion energies. 12The distance of the cavity surface from the solute atoms determines the strength of the electrostatic response and hence the solvation shift in the electronic states.This distance is conventionally fit to reproduce the solvation energies of organic solutes but such models are often unreliable for application to the ionic inorganic surfaces of interest here. 15he SaLSA solvation model 16 directly captures the nonlocal response of the solvent and accurately captures its electrostatic response without any fit parameters, § and is therefore ideal for applications to photocatalyst surfaces which are far removed from the fit sets of conventional solvation models.
However, SaLSA systematically underestimates the solvation energies of organic anions which, although unimportant for the present study, might be relevant for future studies of reaction mechanisms on these catalyst surfaces.The CANDLE solvation model 17 is an empirical solvation model derived from SaLSA which adds two empirical parameters: one to capture the nonlocal response with an effective local response and another parameter to correct for the asymmetric error in anion solvation.CANDLE is systematically more accurate for solvation free energies with a mean absolute error of 1.8 kcal mol À1 on an extensive set of 240 neutral molecules, 51 cations and 55 anions, and 3.5 kcal mol À1 for the anions alone, 17 compared to 4.5 kcal mol À1 and 20 kcal mol À1 respectively for SaLSA on the same set. 16owever, the accuracy of neither SaLSA nor CANDLE has been tested previously for ionic surfaces, so we evaluate the accuracy of both models for predicting the solvation shift on the band positions of several photocatalyst surfaces.Continuum solvation models can describe electrostatic solvent-solute interactions, but ignore stronger chemical bonding effects of the solvent molecules with the surface.In order to investigate these effects, we apply the solvation model both to the bare surface in solution and to the case in which we include a single layer of explicit water molecules adsorbed or chemisorbed on the surface.The solvent monolayer describes chemical bonding well, but without AIMD, it may not accurately describe the solvation shift due to liquid water because fluctuations between numerous low-energy configurations may introduce fluctuations in the interfacial dipoles.To estimate this error, we also consider multiple low-energy configurations of the explicit solvent layer.
Experimentally, the band edge positions in liquid are estimated from the flat-band potential U fb , the applied potential at which the band bending of the semiconductor at the solidliquid interface disappears.The flat-band potential measured relative to a reference electrode, such as the normal hydrogen electrode (NHE), satisfies where A is the band edge positions of the majority carrier referenced to vacuum; which equals to the electron affinity (conduction band minimum E c ) for an n-type semiconductor and the ionization potential (valence band maximum E v ) for a p-type semiconductor.DE F is the difference between the Fermi energy and A, and E NHE 0 = 4.44 eV 37 is the absolute position of the reference NHE potential below vacuum.¶ V pH captures the variation of the flat-band potential with pH and is defined to be zero at the pH of zero charge (pH PZC ).Our calculations do not account for explicit adsorption of ions on the surface; thus they correspond to the pH PZC condition where the net charge of adsorbed ions at the surface is zero.V dip includes only the potential drop across the solid-liquid interface due to the effect of interfacial dipoles.This is exactly what the solvation model captures, and so we can compare V dip directly with its theoretical counterpart: the difference between the electrostatic potential in vacuum and solvated calculations.Note that some studies assume V H V dip + V pH = 0 at pH PZC i.e.V dip = 0, 40 which completely neglects the effect of dipoles at the solid-liquid interface in direct contrast to our findings.
Experimentally, A in eqn ( 1) is obtained from photoemission spectra in vacuum.4][45] We could therefore make theoretical predictions for the absolute band positions in liquid based on eqn (1), but as discussed above, we focus on the solvation shifts (V dip ), since they are insensitive to the choice of the electronic structure method.
Finally, we need to account for the effect of pH on the band edge positions in order to extract the interfacial dipole contribution V dip from experimental data.Fortunately, the band edge positions in an electrolyte follow the Nernst equation for many oxides 3 including TiO 2 46 and WO 3 , 47 because H + and OH À are the only charge-determining ions on these surfaces.Consequently, since we defined V pH to be zero at pH PZC .][3] In summary, we calculate the solvation shift from the electrostatic potential difference between vacuum and solvated calculations, and compare it to the experimental difference in band edge positions from photoemission and electrochemical flat-band potential measurements referenced to pH PZC via the Nernst equation.

Functionalized Si(111) surfaces
Crystalline silicon is an efficient and popular semiconductor for PEC applications because of its suitable band gap and band edge positions relative to water redox potentials.9][50][51][52] Of particular interest is the recent AIMD study of band edge positions and electrostatic potentials in functionalized Si(111)/liquid H 2 O interfaces. 11As discussed previously, although computationally expensive, AIMD describes the dynamical interaction between H 2 O and solid surfaces fairly accurately and serves as an ideal benchmark for our solvation approach before comparison to experiment for more complicated surfaces.
Table 1 compares the solvation shift in the band edge positions, V dip , predicted using the SaLSA and CANDLE solvation models with the AIMD results from ref. 11, for the functionalized Si(111) surfaces shown in Fig. 1.AIMD predicts Table 1 Predicted solvation shifts (V dip ) of functionalized Si(111) surfaces, compared to ab initio molecular dynamics (AIMD) results. 11Positive V dip shifts up the bands towards the vacuum level (lower electron affinity).The final column, -COOH x , presents solvation model results for a partially deprotonated (x = 0.92) surface using snapshots from the AIMD trajectory.The agreement with AIMD for H and CH 3 terminated surfaces validates the solvation models, while the COOH and COOH x cases underscore the importance of identifying the correct surface composition a small solvation shift B0.2-0.3 eV towards the vacuum level for the hydrophobic H and CH 3 terminated surfaces which interact weakly with liquid water.The SaLSA model predictions agree perfectly with AIMD for these surfaces and the CANDLE predictions agree reasonably well with AIMD with a small systematic error BÀ0.08 eV (i.e.further from the vacuum level) compared with the SaLSA results.
The hydrophilic COOH-terminated Si(111) surface, on the other hand, exhibits a large solvation shift 41.5 eV in AIMD and a much smaller shift in the solvation model predictions.However, the AIMD calculations for this surface find significant deprotonation of the COOH group, which implies that the surface is charged on average.Therefore, the AIMD simulation does not correspond to pH PZC , while the solvation model predictions are for the neutral surface at pH PZC , causing the discrepancy.
To verify this hypothesis for the discrepancy between AIMD and the solvation model results for the COOH-terminated surface, we applied the solvation model to partially deprotonated surface configurations extracted from the AIMD trajectory. 11In the 20 configurations including explicit water molecules used in the band position calculations of the hydrated COOH surface in ref. 11, we separate the surface by removing all the water molecules and any H atoms not bound to the surface carboxylate groups.8For each H atom that we remove, we add an electron to the calculation, so that the surface nominally consists of COOH and COO-anions.We find that the COOH terminations are deprotonated approximately 8% on average, and we refer to this surface as COOH x (with x E 0.92).**We then calculate band positions using the solvation models on these charged configurations, and find results in much better agreement with the AIMD results: SaLSA predicts a shift 0.15 eV more positive than AIMD while CANDLE predicts a 0.2 eV smaller shift (see Table 1).
We note that our solvation model is capable of explicitly dealing with charged systems in solution, e.g.Si-COOH x (x E 0.92) surface, because Debye screening in the electrolyte automatically produces a physically-meaningful counter-charge in the liquid that keeps the unit cell neutral and makes absolute potentials well-defined. 35(We do not need to deal with charged surfaces in vacuum because the Si-COOH surface does not deprotonate in vacuum and remains neutral.) The above analysis shows that our solvation model treats the effect of liquid water with remarkable accuracy once we determine the composition of the surface.We show below with comparisons to experimental data that this remains true even for strongly hydrophilic oxide surfaces.

Rutile TiO 2 (110) surface
Having benchmarked our solvation approach against AIMD for functionalized Si surfaces, we proceed with comparisons to experimental data for ionic metal oxide surfaces to demonstrate the generality of our method.Rutile TiO 2 is one of the most well-studied and widely used photoelectrode materials because of its in operando stability, the proper alignment of its band edge positions with the H 2 O redox potentials, and its capability to act as a conductive protective layer for higher efficiency photoelectrodes such as Si. 2,53,54][57][58] We apply the solvation models to the TiO 2 (110) surface both directly and with a single layer of explicit water molecules in order to account for any strong bonding effects.We find two distinct classes of stable H 2 O configurations with comparable energy: the 'flat' configuration with intact water molecules and the 'dissociated' configuration with hydroxides and protons as shown in Fig. 2. In both cases, the additional O binds to the 5-coordinated surface Ti and the H atoms in such positions form a hydrogen bonding network parallel to the surface.† † Table 2 shows that water binds strongly to the surface, with a binding energy of 1.39 eV and 1.24 eV per water molecule in vacuum respectively for the flat and dissociated configurations.The binding energies in solution are smaller because solvation stabilizes isolated water molecules more than adsorbed ones, but note that the binding energy difference between the two configurations is identical in the vacuum, SaLSA and CANDLE calculations.(See ESI, † for details of the binding energy calculations.)The higher binding energy in the flat configuration suggests that the first layer of liquid water would be dominated by the flat motif.Indeed the calculated V dip for this configuration best matches the experimental estimates below.
We calculate the experimental solvation shift as the difference between band positions obtained from photoemission and from electrochemical measurements, as discussed earlier.Ultraviolet photoemission (UPS) measurements 58 find the work function to be 5.2 eV, which is the position of the CBM below the vacuum level since TiO 2 is naturally n-type.Flat-band potential measurements 60 find the CBM at 0.0-0.1 eV relative to NHE, which corresponds to an absolute position of 4.44-4.54eV below the vacuum level.Therefore, using eqn (1) and ( 2) and the experimental pH PZC of 4.5-6.0, 59we find that the experimental estimate of V dip ranges from 5.2-4.54+ 0.059 Â 4.5 = 0.93 eV to 5.2-4.44 + 0.059 Â 6 = 1.11 eV.Additionally, two-photon photoemission measurements 57 find that the work function shift of the TiO 2 (110) surface due to adsorbed water molecules saturates to 0.9-1.0eV after five layers of water molecules, in agreement with the above electrochemical estimates.
Table 2 compares these experimental estimates to the solvation shifts predicted using the solvation model with and without explicit water molecules.The flat configuration of water molecules is the most stable one, and the corresponding calculated shifts lie in the experimental range for the CANDLE model and is overestimated by 0.1 eV for the SaLSA.The calculated shifts for the less stable dissociated configuration disagree with the experimental estimates suggesting that the first water layer in experiment is indeed dominated by the intact water configurations.The solvation shifts due to a monolayer of explicit water molecules alone without the solvation models are smaller by B0.15 eV (0.30 eV) compared to those with CANDLE (SaLSA).In general, fluctuations of the water structure between multiple configurations under room temperature could render predictions from individual minimized structures unreliable, but this is apparently not an issue for TiO 2 because the water binds strongly forming a reasonably rigid structure.A direct application of the solvation models without any explicit water layers, which should in principle directly capture the statistical average of all water configurations, underestimates the experimental shift by 0.2 eV.As in the previous case of functionalized Si surfaces, the CANDLE predictions are systematically 0.1 eV smaller than the SaLSA predictions.

IrO 2 (110) surface
Next, we consider IrO 2 another rutile oxide, which is a highly active and the only known acid-stable catalyst for the oxygen evolution reaction (OER). 61,62IrO 2 is metallic 63 and it is critical to understand the effect of solvation on its work function, since that determines the interfacial energetics and charge transfer between IrO 2 catalysts and semiconductor photoelectrodes. 64e now apply our solvation protocol to estimate the work function shift for the most stable (110) surface 64 of IrO 2 .
Table 3 shows the solvation shifts of the IrO 2 (110) surface with and without a layer of the explicit water molecules.In this case, the water molecules on the surface dissociate to produce OH bound to the surface Ir and H to the surface O.The intact 'flat' configuration is not even meta-stable in this case because it dissociates without a barrier (comparing Fig. 3 to Fig. 2).Additionally, water binds more strongly to IrO 2 than to TiO 2 (110) surfaces, with a binding energy in vacuum of 2.1 eV per water molecule compared to 1.4 eV.(The same is true for binding energies in solution: compare Tables 2 and 3.) The dissociation of H 2 O to form OH groups at the surface of IrO 2 is likely related to its high catalytic activity for OER.
The work function of IrO 2 in liquid has not yet been measured but photoemission measurements in UHV have been    reported. 65We predict a solvation shift of 0.6 eV by directly applying both solvation models to the bare surface, and 0.75 eV for CANDLE (0.79 eV for SaLSA) after including an explicit water layer in the 'dissociated' configuration.These shifts are within 0.2 eV of the similar dissociated configuration for TiO 2 .The solvation shifts due to a monolayer of explicit water molecules alone without solvation models are smaller by B0.20 eV (B0.24 eV) compared to CANDLE (SaLSA).The stability of the intact configuration leads to a 0.4 eV larger shift in TiO 2 , and since that configuration is unstable in IrO 2 , we predict a lower solvation shift E0.75 eV for the IrO 2 (110) surface, as shown in Table 3.

c-WO 3 (001) surface
Stoichiometric WO 3 surface.Finally, we study the effect of solvation on the surface of g-WO 3 (the room temperature phase), another important photoanode material for PEC, which is earth abundant, acid stable, and has a suitable optical gap of 2.6 eV (which is 0.4 eV smaller than that of rutile TiO 2 ). 8,66,67At room temperature, WO 3 adopts a monoclinic phase and its (001) surface with the c(2 Â 2) (or p ffiffi ffi struction has the lowest surface energy. 64][70] We apply the solvation models with and without a layer of explicit water molecules to the stable WO 3 (001) surface.As before, we calculate the intact 'flat' and 'dissociated' water configurations shown in Fig. 4. Their binding energies (in vacuum and in solution) are comparable to those on TiO 2 but much smaller than those on IrO 2 (compare Tables 2-4), which indicates that WO 3 surfaces are not very reactive with H 2 O.In this case, the intact configuration is more stable by 0.23 eV per water molecule than the dissociated one, which is expected since the reconstruction leaves each surface W with two short WQO bonds and two normal W-O bonds, leaving a Lewis acid site for binding an intact H 2 O. Regardless, the solvation shifts for the two configurations, also shown in Table 4, differ only by 0.2 eV and the predictions without explicit water also agree within 0.2 eV.The solvation shifts due to a monolayer of explicit water molecules alone without solvation models are smaller by B0.2 eV (B0.3 eV) compared to CANDLE (SaLSA).Overall, we find theoretical solvation shifts for WO 3 to be B1.0 eV.
Experimentally, pristine WO 3 surfaces have very high ionization energies B9.6-9.8 eV, 71,72 which can decrease dramatically by B2 eV due to adsorption of molecules such as water when the sample is exposed to air. 73Electrochemical flat-band measurements find the VBM at BÀ3.0 to À3.1 eV relative to NHE (7.44-7.54eV below vacuum). 1 Therefore, the above experimental results combined with eqn ( 1) and ( 2), as well as pH PZC = 0-1 59 indicate that the experimental estimate of V dip ranges from 9.6-7.54+ 0.059 Â 0 = 2.06 eV to 9.8-7.44 + 0.059 Â 1 = 2.42 eV, which is more than twice that of the theoretical predictions!Based on our findings for Si and TiO 2 , we do not expect such large discrepancies.Additionally, our predictions conform with the 0.5-1.0eV range of solvation shifts observed for most oxides in literature, 10,57,74 but the experimental estimates of the solvation shift of WO 3 do not fall in this range.This indicates the possibility that additional phenomena are at play in the solvation of WO 3 surfaces.photoemission measurements. 71,72STEM and LEED studies also observed O-deficient surfaces, such as a p(2 Â 2) reconstruction 18 with O atoms on top of 1/4 surface W atoms (compared to 1/2 at the most stable c(2 Â 2) surface (Fig. 5)), or a p(1 Â 1) reconstruction with no on-top O atoms.‡ ‡ 70 Next to a liquid electrolyte, surface O vacancies could be further stabilized by bonding with water molecules.These bonds could contribute additional dipoles to the interface and influence the solvation shift.Direct experimental measurements of O vacancies on solvated surfaces would be extremely challenging and have not yet been reported.Consequently, in the remainder of this section, we discuss the formation energies of surface O vacancies in vacuum and solution, and their effect on the band positions, to explore this possibility.Table 5 compares the formation energy of O vacancies in bulk WO 3 as well as on the (001) surfaces in vacuum and solution at different O vacancy concentrations.The formation energies are presented for two limits: the O-rich limit where molecular oxygen provides O and the O-poor limit where the oxide (or water, when available) provides it.(See ESI, † for details on the vacancy formation energy calculations and additional discussions.)The formation energies of bulk O vacancies in WO 3 are positive and \1 eV over the entire range of O-rich to O-poor conditions.In comparison, the formation energies of surface vacancies in vacuum are much smaller than the bulk case (by B2 eV) and even become negative in the O-poor limit, indicating that O deficiency could be thermodynamically favorable in this limit. 75Interestingly, the formation energy of O vacancies on WO 3 (001) surfaces in solution are further stabilized by 0.1-0.4eV compared to those in vacuum, which indicates higher concentrations of O vacancies on the solvated surfaces.
Now that we have established the plausibility of high concentrations of O vacancies on solvated WO 3 surfaces, we explore the effect of these vacancies on the band edge positions in solution.Table 6 and Fig. 6 show the total band position shift due to interfacial dipoles, including the effect of vacancies and solvation (when present), for various concentrations of vacancies on the WO 3 (001) surface in vacuum and in solution.As expected, the potential shift varies linearly with vacancy concentration because the surface dipole per vacancy approaches a constant in the dilute limit, and the potential shift is proportional to the induced surface dipole per unit area.(Zero vacancy concentration corresponds to the ideal surface calculations from the previous section.) The surface dipole per vacancy, and hence the shift at identical vacancy configurations, is much larger in solution than vacuum because water molecules bind at the vacancy site and the extra H atoms (relative to the surface without vacancies) provide a large positive dipole and shift the bands towards the vacuum level.In this case, the shifts predicted using explicit water molecules alone underestimate those including the solvation model by up to 0.8 eV, unlike the previous underestimation B0.2 eV for the stoichiometric surfaces.This shows the importance of accounting for the effect of the liquid beyond the first layer using solvation models, especially for complicated surface structures including vacancies.
The total shift with the c(4 Â 4) surface vacancy configuration in solution of 2.1 eV (Table 6) matches the experimental  This journal is © the Owner Societies 2015 solvation shift of B2.0-2.4 eV (Table 4).This configuration corresponds to a surface vacancy concentration of 25%, which is high, but plausible as discussed above: the surface vacancy formation energies of WO 3 in solution are lower than in vacuum, and such high surface vacancy concentrations of WO 3 have been observed even in UHV STM and LEED experiments. 18,19xperimental characterization of the WO 3 surface structure in the electrochemical environment is necessary to confirm our predictions.
Optimum alignment of WO 3 band edge positions with the water redox potentials requires shifting both the VBM and CBM further towards the vacuum level.Based on the above results, this could potentially be achieved by increasing the O vacancy concentration at the surface.Several surface preparation techniques lead to highly O deficient WO 3 surfaces in vacuum. 18,70t may be possible to retain this O deficiency in the aqueous environment, especially under acidic conditions where WO 3 is stable.

Band position shifts in non-aqueous solvents
The previous section indicates the possibility of adjusting band positions by controlling surface preparation.Alternately, band positions would also vary with solvent, which could be particularly important for carbon dioxide reduction (CO2RR) in non-aqueous solvents.Here, we present a preliminary study of the variation of band positions with solvent, which has not been previously studied either experimentally or theoretically.For simplicity, we focus on the hydrophilic IrO 2 (110) surface since it has a non-negligible solvation shift in water that could be described by the solvation models accurately without including explicit solvent molecules.We include a range of non-polar to polar solvents for which the SaLSA or CANDLE solvation models have been previously constructed: carbon tetrachloride (SaLSA), chloroform (SaLSA) and acetonitrile (CANDLE).
We expect intuitively that the effect should increase with the dielectric constant, e, of the solvent and decrease with the distance, d, of the solute charge from the dielectric surface.In fact, for a point charge in a dielectric cavity of radius d, the solvent shift in the potential is proportional to (1 À e À1 )d À1 .The solvation shift V dip for a surface depends on its charge distribution in a complicated manner, but we can use the above spherical-cavity form to examine the correlations, with the solvent vdW radius 76 as an estimate of the distance d.Fig. 7 shows the solvation shifts in the band positions of the IrO 2 (110) surface.Indeed the shifts correlate well with the solvent dielectric function and solvent molecule size in the form (1 À e À1 )R vdW À1 , as discussed above.It also shows that the shift does not correlate well with dielectric constant alone, neither linearly nor inversely.Water sets the upper bound on solvation shifts, since it has a high dielectric constant and a small molecule size; increasing the dielectric constant further does not have an appreciable effect due to the (1 À e À1 ) dependence.However, it is possible to reduce the solvation shifts by selecting solvents with either larger molecules or smaller dielectric constants.

Discussion
The effect of solvation on band positions is an important consideration in determining the alignment of semiconductor band edges with redox potentials across the solid-electrolyte interface.The magnitude of the solvation shift varies tremendously from one system to another, but our results indicate that it is possible to estimate them based on the nature of the surface (specifically its hydrophobicity).Non-ionic hydrophobic surfaces, such as the hydrogen and methyl-terminated silicon surfaces in our study, interact weakly with water and exhibit small solvation shifts less than 0.5 eV.Vacuum calculations might suffice for such surfaces, especially for qualitative estimates seeking an accuracy B0.5 eV, as discussed in ref. 11.
Ionic hydrophilic surfaces interact strongly with water, and especially for oxides, the solvation shift of the band edge positions tends to fall in the 0.5-1.0eV range.In general, predictions using vacuum calculations alone can be qualitatively incorrect for these surfaces and therefore the effect of H 2 O should be included.
In some cases, defects on the surface can lead to larger solvation shifts, as we found for oxygen vacancies at the surface of WO 3 .In this case, water molecules bind to the oxygen vacancy sites with the hydrogen atoms pointing outwards, which contributes a large positive dipole that can shift the electronic states by as much as 2 eV.Such surfaces require particular care, and calculations should account for restructuring of the surface and should analyze concentrations of possible surface defects, in addition to accounting for direct solvation effects.
In terms of the performance of the solvation models for predicting the solvation effects of the band positions of various surfaces, we find that the CANDLE and SaLSA continuum solvation models can accurately describe the solvation shifts for hydrophobic surfaces, when applied directly to the surface with no explicit water molecules.They agree qualitatively with experimental results for hydrophilic surfaces, but miss the stronger bonds of the surface with water and consistently underestimate the solvation shift by B0.2 eV.The agreement is considerably improved by including a layer of explicit water molecules, where these strong interactions are now handled by electronic density functional theory instead, but this introduces the complexity of dealing with many possible water configurations.The solvation models are adequate to describe the effect of subsequent layers of water, which keeps the number of water configurations manageable and mitigates the need for explicit statistical sampling via molecular dynamics.Additionally, we could also reparametrize empirical solvation models such as CANDLE to remove this underestimate: this could be useful for accurate high-throughput theoretical screening of PEC materials.
In addition, we also find that the results of the empirical CANDLE solvation model agree very well with the non-empirical and more computationally-expensive SaLSA model.This is particularly important for studies of reaction mechanisms with charged species, since CANDLE corrects the systematic under-solvation of anions in SaLSA (and most other solvation models).We also note that SaLSA and CANDLE show a systematic offset of B0.1 eV in the predicted solvation shifts.In fact, correlating theoretical and experimental potentials of zero charge of single crystalline metal surfaces predicts E NHE 0 = 4.55 eV for the SaLSA model and 4.66 eV for the CANDLE model, 17 which also differ by B0.1 eV (and agree well with the experimental estimate B4.44 eV 37 ).This is, therefore, a systematic difference in the dipole layer of the liquid at the interface, independent of the solid surface.
Finally, our results indicate two avenues for engineering band positions of surfaces in solution.First, surface preparation techniques that produce oxygen-deficient surfaces may be used to shift band positions of oxides e.g.WO 3 further towards the vacuum level.An experimental investigation of whether such surfaces remain stable in water is necessary to determine the feasibility of this approach in order to optimize band alignment with water redox potentials.Second, solvents consisting of larger molecules with lower dielectric constants may be used to reduce the solvation shift.This is particularly relevant for carbon dioxide reduction, which may be carried out in non-aqueous solvents.Additionally, solvent effects can modify the interface energetics and charge transfer between photoelectrodes and catalysts by introducing a porous structure of the catalysts, where the water can penetrate the catalyst layer and be in contact with photoelectrodes, e.g. in WO 3 /IrO 2 interfaces. 62Such effects have been successfully predicted using atomistic models of the interface combined with solvation models, 64 and can be systematically exploited to improve the efficiency of the PEC.

Conclusions
We studied the solvation shifts of the electronic states of various surfaces using recently-developed continuum solvation models suitable for such systems.This computationally-efficient approach agrees well with much more expensive AIMD simulations of explicit solid-liquid interfaces, and is far more accurate than vacuum calculations alone while requiring minimal computational overhead.§ § Accounting for strong surface-liquid interactions such as formation of chemical bonds by using a single layer of explicit water molecules, we found good agreement with experimental band positions for highly ionic TiO 2 and WO 3 surfaces, and predicted the as-yet unmeasured Fermi level position of IrO 2 in liquid water.The solvation shifts vary considerably depending on the nature of the surface and the solvation model predictions are typically accurate to within B0.2 eV, as summarized in Table 7.Most importantly, we correlated the magnitude of the solvation shift with the surface type, the dielectric constant and molecule size of the solvent, and the strength of their mutual interaction, which provides valuable guidance for the discovery of new efficient photocatalyst materials.

Fig. 1
Fig. 1 Geometry of (a) hydrogen, (b) methyl and (c) carboxylic acid terminated Si(111) surfaces.Calculations used an inversion-symmetric (1 Â 1) surface unit cell with four layers and atom positions optimized self-consistently in both vacuum and solution.

Fig. 2
Fig. 2 Geometry of (a) bare TiO 2 (110) surface and surfaces with (b) flat adsorbed water molecules and (c) dissociated adsorbed water molecules.Calculations used an inversion-symmetric (1 Â 1) surface unit cell with six layers and atom positions optimized self-consistently both in vacuum and solution.

Fig. 3
Fig. 3 Geometry of (a) bare IrO 2 (110) surface and surface with (b) dissociated adsorbed water molecules.The intact flat configuration is unstable and dissociates spontaneously to configuration (b).Calculations used an inversion-symmetric (1 Â 1) surface unit cell with six layers and atom positions optimized self-consistently in both vacuum and solution.

WO 3
surface with oxygen vacancies WO 3 is intrinsically n-type due to oxygen deficiencies, which sets its Fermi level close to (0.2-0.3 eV below) its CBM in

Fig. 4
Fig. 4 Geometry of (a) bare WO 3 (001) surface and surfaces with (b) flat adsorbed water molecules and (c) dissociated adsorbed water molecules.Calculations used an inversion-symmetric (1 Â 1) four layer surface unit cell with self-consistent atom positions in solution.

Fig. 5
Fig. 5 Schematic top view of configurations of O vacancies on the WO 3 (001) surface ((a) 0%, (b) 50%, (c) 25%, (d) 12.5%).The solid yellow circles are the surface O atoms, empty yellow circles mark the O vacancies and black arrows indicate the lattice vectors of the supercell.For clarity, the monoclinic distortions (see Fig. 4) have been omitted in this illustration.

Fig. 6
Fig. 6 Variation of band position shift of WO 3 (001) surfaces due to vacancies and/or solvation as a function of vacancy concentration (per surface O atom of the stoichiometric surface) in vacuum and solution (see Table6).Vacancy concentrations of 0.125, 0.25 and 0.5 correspond to the p(4 Â 4), c(4 Â 4) and p(2 Â 2) configurations respectively.

Fig. 7
Fig. 7 Correlation of the predicted solvation shift of IrO 2 (110) band positions in various solvents with various combinations of solvent dielectric constant e and vdW radius R vdW .

Table 2
1,2,58,59 solvation shifts (V dip ) of the TiO 2 (110) surface with and without explicit H 2 O molecules, compared to experimental shift of the VBM due to liquid H 2 O,1,2,58,59and the binding energies of explicit water molecules to the surface (per molecule), E

Table 3
Predicted solvation shifts (V dip ) of the IrO 2 (110) surface with and without one layer of explicit water molecules, and binding energies per explicit water molecule, E

Table 4
1,71,72ed solvation shifts (V dip ) of the stoichiometric WO 3 (001) surface with and without explicit H 2 O molecules, compared to experimental shift of the VBM due to liquid H 2 O,1,71,72and binding energies per explicit water molecule, E H 2 O bind .The large discrepancy between our theoretical predictions E1 eV and experimental measurements 42 eV suggest that the experimental surfaces may not be the stoichiometric defect-free surface we studied

Table 5
Formation energies of various concentrations of oxygen vacancies in bulk WO 3 and WO 3 (001) surfaces in vacuum and solution (explicit water layer + CANDLE solvation model).The O-rich limit corresponds to deriving O from molecular oxygen whereas the O-poor limit derives O from WO 3 or H 2 O (when present)

Table 6
Total shift of the band positions of WO 3 (001) due to interfacial dipoles, V dip+vacancy including contributions due to O vacancies and/or solvation, for various configurations of surface O vacancies in vacuum and solution.The H 2 O-ML calculations include a monolayer of explicit water molecules without any solvation model, whereas the CANDLE calculations include the CANDLE solvation model in addition to the explicit monolayer

Table 7
Summary of solvation shifts[eV]in band positions of various systems: solvation model predictions compared to experiment and AIMD results.Here, SM refers to solvation models alone (range includes SaLSA and CANDLE predictions), ML is monolayer of explicit water alone and ML + SM is solvation models in addition to explicit monolayer.See Tables 1-4 and 6 for more details