Electronic structure of aqueous two-dimensional photocatalyst

The electronic structure, in particular the band edge position, of photocatalyst in presence of water is critical for photocatalytic water splitting. We propose a direct and systematic density functional theory (DFT) scheme to quantitatively predict band edge shifts and their microscopic origins for aqueous 2D photocatalyst, where thousands of atoms or more are able to be involved. This scheme is indispensable to correctly calculate the electronic structure of 2D photocatalyst in the presence of water, which is demonstrated in aqueous MoS2, GaS, InSe, GaSe and InS. It is found that the band edge of 2D photocatalysts are not rigidly shifted due to water as reported in previous studies of aqueous systems. Specifically, the CBM shift is quantitatively explained by geometric deformation, water dipole and charge redistribution effect while the fourth effect, i.e., interfacial chemical contact, is revealed in the VBM shift. Moreover, the revealed upshift of CBM in aqueous MoS2 should thermodynamically help carriers to participate in hydrogen evolution reaction (HER), which underpin the reported experimental findings that MoS2 is an efficient HER photocatalyst. Our work paves the way to design 2D materials in general as low-cost and high-efficiency photocatalysts.


INTRODUCTION
Generating hydrogen by photocatalytic water splitting through photoelectrochemical cells (PECs), as first demonstrated in the pioneer work by Fujishima and Honda 1 , is a promising approach for producing environmentally friendly and scalable solar fuel. To assist photocatalytic reactions, highly efficient sunlight absorbers and suitable catalysts are required. The key is therefore to find an earth-abundant semiconductor for sunlight absorption, charge carrier separation, transport, as well as surface redox reaction. The photocatalytic process is initiated by photo-excitation of electrons from the valence band of the semiconductor to the conduction band. The semiconductor should have an optimal band edge that straddles the water redox potential, namely, its conduction band minimum (CBM) should be higher than the hydrogen reduction potential (H 2 /H 2 O) to accomplish hydrogen evolution reaction (HER), and its valence band maximum (VBM) should be lower than the water oxidation potential (O 2 /H 2 O) to accomplish oxygen evolution reaction (OER). Then, the electrons (holes) can migrate to the surface active sites to reduce (oxidize) water into hydrogen (oxygen). There are two popular cells for photocatalytic water splitting, the photochemical cell (PC) 2,3 and the photoelectrochemical cell (PEC) [2][3][4] . In the former, both HER and OER occur on the same photocatalyst and the product is a mixture of H 2 and O 2 ; in the latter, the two processes occur on different semiconductor electrodes and H 2 and O 2 are separated automatically. Thus, for PEC we concerned in this work, each semiconductor electrode provides part of the water splitting potential, i.e., either HER one, or OER one 5 .
Recently, the advent of two-dimensional (2D) materials presents fascinating properties and new possibilities for photocatalytic water splitting (i.e., MoS 2 ,GaS, InSe, GaSe, InS, C 3 N 4 , BiOCl, and ZnIn 2 S 4 ). In 2D materials, the surface-to-volume ratio is maximized to provide more reaction sites for photocatalytic reaction. The 2D nature shortens the distance for photo-generated electrons and holes to migrate to the reactive surface, greatly lowering the electron-hole recombination 6 . To find and/or design 2D materials for water splitting, a critical problem is to establish the band edge positions in the presence of water, since interaction with water can dramatically affect electronic properties of the 2D materials 7-11 . However, due to the complexity of calculating the electronic structure of the solid surface in contact with water, density functional theory (DFT) calculations of the band edges of the aqueous photocatalysts are usually performed without any water, accordingly, neglecting the structural and chemical effect brought by the solid/water interface. So far, very few studies have simulated the entire solid/water interface at DFT level to obtain the solid's band edge in the presence of water 7,10-13 . These few existing DFT studies proceed by first calculating VBM and CBM of the solid without water using their primitive cell, then shifting the VBM and CBM by the same value, i.e., the average electrostatic potential (AEP) offset calculated from the solid/water interface 14 . The final result is taken as the band edge of the combined solid/ water system 10,11,15 . This indirect approach significantly reduces computation demand since the electrostatic potential can be easily obtained in the self-consistent loops of the DFT calculation, where small supercells with limited number of water molecules are used. And the primitive cell employed in the first step allows high level exchange-correlation to be employed 10 . However, noticeably, this indirect approach assumes that the VBM and CBM of the solid are rigidly shifted by the same amount of energy due to the presence of water. For 2D photocatlysts/water systems, an interesting issue is if the assumption in such an indirect procedure still holds true. Additionally, the microscopic origins of these band edge shifts need to be studied quantitatively.
To determine band edge positions of 2D materials in the presence of water and answer the above issues, as well as to design 2D photocatalysts in general, direct DFT calculations to determine the band edge shifts of the combined 2D-photcatlysts/ water system are indispensable. To achieve this goal, it is necessary to solve large number of atoms in the combined system. In this work, we meet this challenge by using our recently developed state-of-the-art real-space Kohn-Sham (KS) DFT solver that can self-consistently simulate thousands or more atoms 16 . To be more specific, we choose MoS 2 as an example of 2D photocatlysts to demonstrate our theoretical approach and illustrate the aforementioned key issues. MoS 2 is a particularly interesting catalyst because it is less costly than noble-metals and very stable under acidic conditions 17 . Experimentally, MoS 2 showed excellent intrinsic HER catalytic activity 18,19 . Moreover, the band edge of monolayer MoS 2 without water is typically considered to straddle the water redox potential 6 .
In this work, we propose a direct and systematic DFT approach to quantitatively predict band edge positions, band edge shifts and microscopic origins of these shifts for realistic 2D photocatalyst/water systems, where thousands of atoms or more are involved.In previous work 7, 11 , the origins of the band edge shift induced by water are classified into three parts: geometrical distortion, water dipole and electron redistribution.However, the first two origins are just qualitatively discussed while only electron redistribution effect is quantitatively calculated. Here, we establish a set of elegant procedure to quantitatively and directly calculate the contributions to the band edge shift from the aforementioned three microscopic origins, which is applicable to a wide range of 2D material systems. By directly calculating the electronic structure of the combined MoS 2 /water system, we discover that the CBM and VBM of MoS 2 are shifted by different values due to water. The value of CBM shift perfectly matches with the sum of the contributions from geometrical distortion, water dipole and electron redistribution. However, these three elements could not explain the VBM shift. Thus, the fourth effect, the chemical contact between 2D photocatalyst and water, is disclosed in our work. Four more promising 2D photocatalysts, GaS, InSe, GaSe and InS 6,20-22 are investigated, and we found that it is rather general that under aqueous condition, the VBMs of 2D materials whose surface atoms contain out-of-plane orbitals will shift differently from their CBMs, and the fourth origin, chemical hybridization dominates this difference.

RESULTS
Direct DFT scheme for band edge shift analyses Before presenting the calculated data, we would like to discuss our direct DFT approach (shown in Fig. 1) to analyze the band alignment, the band edge shift and microscopic origins for these shifts, which underpins our following findings. First, we calculate the pristine 2D material without water whose band edge relative to the vacuum level is shown in Fig. 1a; then we calculate water without 2D material, and the AEP in the bulk region of water, V H , relative to the vacuum level by ΔV 0 is shown in Fig. 1b. Second, we calculate the combined 2D material/water system in which the bulk water part has its own AEP, V 0 H , as shown in Fig. 1c, and we do a global shift by setting V 0 H ¼ V H . This way, because V 0 H is aligned with V H (Fig. 1c), and V H relative to the vacuum level is known from Fig. 1b, the band edge of the 2D material in the 2D material/ water system is correctly aligned to the vacuum level. Finally, the influence of water to the band edge shift of 2D material is obtained by comparing Fig. 1a, c.
To understand the influence of water on the band edge shift of 2D material in the combined 2D material/water system, we shall quantitatively decompose the effect into three sources: the geometric deformation of 2D material due to water, the water dipole potential and the charge redistribution. Figure 1d is used to calculate the band edge of the distorted 2D material, where the distortion is induced by water and the distorted configuration is extracted from the fully relaxed 2D material/water interface system. By comparing the band edge calculated from Fig. 1a, d, the effect of geometric deformation is obtained. Similarly, owing to the interaction with 2D material, the arrangement of water molecules surrounding 2D material is greatly changed, usually presenting a well-organized double layer structure, i.e., the wellknown water double layer 23 . Combined with the fact that intrinsically water molecule has a strong electric dipole, the dipole orientation of the reorganized water molecules should dramatically shift the band edge of 2D material. To quantitatively estimate the effect of water dipole, Fig. 1e shall be employed. By removing 2D material from the combined 2D material/water system, the electrostatic potential profile and ΔV 1 of the remaining water part can be calculated. ΔV 1 is the difference between the average electrostatic potential (labeled by dashed line in Fig. 1e) in the marked vacuum region (corresponding to the 2D material region in the combined 2D material/water system) and bulk water region. The band edge shift due to water dipoles is thus determined by ΔV 1 − ΔV 0 (Fig. 1b). With respect to the third effect, the charge redistribution brought by the interfacial interaction between 2D Water Vacuum level  The averaged electrostatic potential of bulk water V H relative to the vacuum level, resulting in ΔV 0 . c The band edge of 2D material in the combined 2D material/water system, where the AEP in the bulk water region is V 0 H , which is aligned with V H . d The band edge of deformed 2D material relative to the vacuum level. e The AEP of the water part extracted from the combined 2D material/water system. ΔV 1 is the corresponding averaged electrostatic potential difference between 2D material region and bulk water region. f The AEP generated from the differential charge density of the combined 2D material/water system. V 000 H is the AEP of bulk water region and ΔV DCD is the difference between the AEP of 2D material region and V material and water. we firstly quantify this charge redistribution by differential charge density (DCD), which is defined as ρ i corresponds to the charge density of combined 2D material/ water system, ρ water is the charge density of the water part extracted from the combine system, ρ 2D represents the charge density of the extracted 2D material part. As is illustrated in Fig. 1f, by solving possion equation, electrostatic potential change induced by such charge redistribution is obtained. ΔV DCD denotes the difference between the averaged electrostatic potential of 2D material region and that of bulk water region where DCD equals zero. Following the direct DFT scheme for band edge shift analyses we proposed here, specifically, a systematic and detail investigation on electronic structures of aqueous MoS 2 , a typical example of 2D photocatlyst, are performed. A finite slope of electrostatic potential distribution could be induced by the asymmetry of water layers contacting the two surfaces of MoS 2 in a single sampled MoS 2 /water configuration. Note that in Fig. 1, the averaged results of electrostatic potential distribution over sufficient samples are plotted as schematics for analyzing the band alignment and band edge shift of semiconductors in the presence of water.  Fig. 1b we construct a water/vacuum system and calculate V H , which is relative to the vacuum level by ΔV 0 . Some further details of this calculation is given in the Supplementary Note 2, in particular we found ΔV 0 = 2.73 eV ( Supplementary Fig. 4). While the value of ΔV 0 is calculated to bẽ 3.7 eV as reported in refs 24,25 due to the use of different pseudopotentials, the choice of pseudopotential does not change the MoS 2 /water level alignment.
Finally, we calculate the combined MoS 2 /water system as shown in Fig. 1c. To this end, we construct supercells using atomic configurations from the LAMMPS MD simulation as discussed in the method section, in particular we sample 20 snapshots equally spaced in 10 ns time span after equilibration of the MD simulation. We adopt periodic boundary condition in calculating the electronic structure of the combined system. The V H (Fig. 1b) of the bulk water sets the energy reference, and we align the V combined MoS 2 /water system to the vacuum level. The following equations describe this analysis: where Eq. (2) is used to compute V H to determine the AEP of bulk water with respect to vacuum level, Eq. (3) is employed to align AEP of bulk water in MoS 2 /water system with respect to AEP in the water/vacuum system, Eqs. (4) and (5) align the CBM and VBM with respect to AEP of bulk water in the MoS 2 /water system. In this way, the CBM and VBM of the MoS 2 /water system is aligned with respect to vacuum level. Figure 2a plots the CBM, VBM of the 20 snapshots of atomic configurations, obtained using Eqs. (4) and (5). The CBM (red square) and VBM (blue circle) fluctuate in unison as a function of the MD simulation time, reflecting the dynamics of the water molecules. The timescale of photocatalytic chemical reactions is typically longer than μs [26][27][28] . The fluctuation of the band edges is on the timescale of ns (Fig. 2). It is thus appropriate to work with the averaged edge potentials for band alignment studies on aqueous 2D photocatalysts. Figure 2b gives the band edge shift of MoS 2 due to water in the combined MoS 2 /water system for the 20 snapshots. It is obtained by subtracting the CBM/VBM values of pristine MoS 2 from those of the combined MoS 2 /water system. Both the CBM and VBM of the combined MoS 2 /water system display a significant upshift (the shifted values are positive in Fig.  2b). This is very important for MoS 2 as a photocatalyst: the value of the upshift suggests this material to be effective for HER but not OER.
To understand the microscopic origin of the band edge shift of MoS 2 due to water, we quantitatively calculated the band edge shift induced by geometric deformation, water dipole and charge redistribution effect following the direct DFT scheme discussed above.

Geometric deformation effect
We start by analyzing the CBM of MoS 2 in the combined MoS 2 / water system. The upshift value of CBM of the 20 snapshots is plotted in Fig. 2b (red squares) whose average is 0.42 eV. Here we investigate the contribution toward this value due to geometric deformation of MoS 2 in the presence of water. Turns out this contribution adds a negative contribution.
Starting from a given LAMMPS MD configuration of the combined MoS 2 /water system, we remove the water in the supercell, and compute the CBM of the deformed MoS 2 using Fig.  1d. The CBM position of the distorted MoS 2 is plotted in Fig. 2c (denoted by red squares, CBM-d). We found that the CBM of the distorted MoS 2 slightly shifts downwards with an average value of −0.04 eV, with respect to the CBM of pristine MoS 2 . The downshift value as a function of MD simulation time is plotted in Fig. 2d (as red square). We conclude that the geometric deformation of MoS 2 in the combined MoS 2 /water system contributes a negative value to the total band edge shift of CBM.
The formation of ripples in freestanding MoS 2 can have dramatic effect on its band edges 29 . In our MD simulation, the ripple effect is included since all atoms are dynamically evolved. Because MoS 2 is inside water, different from the freestanding situation, the rippling is largely suppressed in the MoS 2 /water system.
Water dipole effect Next, the dipole orientation of the water molecules are analyzed for the combined MoS 2 /water system, which turns out to have a dramatic effect on the band edge of MoS 2 . First, the number density of hydrogen and oxygen atoms averaged over the 10 ns MD simulation times is plotted in Fig. 3a where the sulfur surface of MoS 2 towards water is located at~5 Å. Due to the Coulomb interaction, water molecules near the MoS 2 surface are organized into the well-known water double layer. About 10 Å away from the MoS 2 surface, the density of bulk water is recovered. The density in the first peak of the double layer, which is nearest to the surface is about three times of that in bulk water region, consistent with previous literature 30 . The structure of water near the MoS 2 surface is analyzed by the statistical information of dipole orientation of the water molecules, i.e., the probability of dipole orientation shown in Fig. 3b defined as the H-O-H angle bisector. The angle θ in Fig. 3b is defined as the angle between the dipole of water molecule and the Z-direction, which is perpendicular to the MoS 2 surface, thus the dipole of water molecule is pointing to the MoS 2 surface if θ is greater than 90°and pointing out of the MoS 2 surface if θ is less than 90°. In the water double layer region, it is clear that the dipole orientation is no longer distributed randomly like that of bulk water region. For water molecules nearest to the MoS 2 surface (at about 7.5 Å in Fig. 3b), the dipoles tend to point to the MoS 2 surface with θ = 150 ∘ . This is consistent with the DFT calculation for a single water molecule adsorbing on the surface of MoS 2 31 . Note that the plot in Fig. 3a is a statistically averaged result of 10 6 configurations. In a specific configuration, the water molecule nearest to the MoS 2 surface tends to point their hydrogen atom to MoS 2 . Actually, in Fig. 3a, the blue line (hydrogen) begins to increase at a z-direction value of 6.7 Å, while the red line (oxygen) begins to increase at a larger value of 7.7 Å. This means that atoms with the shortest distance to MoS 2 are hydrogen atoms, i.e., the nearest water molecules point their hydrogen atom towards MoS 2 surface. The region from 7.9 to 9.1 Å corresponds to the first density peak in Fig. 3a, where the gradual transition of water dipole orientation from pointing to the MoS 2 surface to pointing away of the surface is observed. A similar transition is found in the region from 10.9 to 12.6 Å, which corresponds to the second density peak in Fig. 3a. Such ordered arrangement of water dipoles in the water double layer induces an important electrostatic potential shift in MoS 2 .
The band edge shift of MoS 2 caused by water dipole is then evaluated. After the combined system is converged by DFT, we remove MoS 2 from the supercell and calculate the electrostatic potential profile and ΔV 1 of the remaining water part, as schematically shown in Fig. 1e, which includes the contribution from the water double layer. The potential shift due to water dipoles is thus determined by ΔV 1 − ΔV 0 where ΔV 0 = 2.73 eV as discussed above. For the 20 snapshots, the dipole induced shifts are presented in Fig. 2d (cyan triangles), showing a prominent upshifts for all the atomic configurations with an average of 0.21 eV. In Fig. 4a, we plot this dipole induced band edge shift and the total CBM shift together, which allows us to conclude that the fluctuated behavior of CBM is largely due to the water dipole effect.

Charge redistribution effect
A third factor causing the CBM shift is charge redistribution induced by the interaction between MoS 2 and water, which results in a remarkable influence on the observed band edge shift of MoS 2 . Following Eq. (1), the DCD of this system is calculated, where the 2D material is replaced by MoS 2 specifically. Fig. 3c plots the DCD for a sampled configuration at the MD simulation time of 1 ns. Clear charge accumulation (blue) is seen near sulfur atoms at the MoS 2 /water interface, and charge reduction (purple) around hydrogen atoms. This indicates non-negligible electron D. Kang et al.
transfer from water to MoS 2 in the combined MoS 2 /water system. The average DCD along the direction perpendicular to MoS 2 /water interface is plotted in Fig. 3d, which exhibits positive electron changes (0.52 e) in the MoS 2 and negative changes in the water. As illustrated in Fig. 1f, by solving the Poisson equation, the potential shift derived from DCD, denoted as ΔV DCD , is calculated. Figure 2d (green triangle) gives this potential shift due to charge redistribution in the combined MoS 2 /water system for the 20 snapshots, the average is found to be 0.26 eV.
When summing the potential shifts due to the calculated three microscopic sources, i.e., the water dipole effect (0.21 eV), the charge redistribution effect (0.26 eV), and the geometric deformation effect (−0.04 eV), the total potential shift is 0.42 eV, which is exactly the upshift of CBM of MoS 2 in the combined MoS 2 /water system, as shown in Fig. 4b. The quantitative match not only indicates the accuracy of our direct DFT approach to evaluate each microscopic contribution to the band edge shift, but also states that for CBM of MoS 2 in the presence of water, its band edge shift is largely due to the interfacial electrostatic interaction with a small negative contribution from the geometric structure deformation of the atoms.
The fourth effect: chemical contact So far we focused on the CBM shift which is critical for HER. Using the same procedure, the VBM shift for MoS 2 in the combined MoS 2 /water system can be analyzed. First, in the presence of water, the VBM of MoS 2 is upshifted with an averaged value of 0.53 eV (Fig. 2b), which is different from the value of CBM shift (0.42 eV). Then, with respected to the three origins of the band edge shift owing to water: geometric deformation, water dipole and charge redistribution, their contributions to the upshift of VBM are 0.04, 0.21, and 0.26 eV, respectively. Obviously as plotted in Fig. 4c, the VBM shift does not match perfectly with the sum of these three effects. An additional upshift of about 0.03 eV is found for VBM of MoS 2 in the combined MoS 2 /water system. In short, due to water, compared with the CBM shift, the VBM of MoS 2 is upshifted differently with a value of 0.11 eV. This difference comes from two parts: the geometric deformation with a difference of 0.08 eV and an additional unknown effect causing 0.03 eV. The converged shift difference is about 0.10 eV, which is ten times greater than the statistical error, i.e., 0.01 eV To determine the additional unknown effect, the wavefunction distributions of VBM, CBM in real space for the combined MoS 2 / water system, are plotted in Supplementary Fig. 3. Different from CBM, the wavefunction of VBM not only distributes on Mo atoms and S atoms but also on water molecules. There are orbital contributions from water molecules for VBM but not for CBM. It can be well understood from the orbital projection of VBM and CBM of MoS 2 . Since S atoms locate at the interfacial sites closer to the water molecules, we focus on the orbital projection from S atoms. The CBM of MoS 2 is composed of in-plane orbitals (p x , p y ) of S atoms, and VBM composed of the out-of-plane orbital (p z ). The interactions between the p z orbital of S atoms in MoS 2 and p z orbital of O atoms in water give rise to an emerging hybridized VBM state in the combined MoS 2 /water system, which mimics a chemical-like contact. On the other hand, for CBM, there is no outof-plane orbitals of S atoms to interact with p z orbtials of water. We conclude that in the presence of water, the VBM of the combined system is provided by this hybridized VBM state sitting Fig. 3 Water dipole and differential charge density for the combined MoS 2 /water system. a Laterally averaged number density of oxygen and hydrogen atoms for 10 6 configurations over a 10 ns time span of the MD simulation. b Dipole orientation of water molecules in the combined MoS 2 /water system. c Differential charge density (DCD) for the sampled configuration at the MD simulation time of 1 ns. The isosurface value is set at 6 × 10 −9 eÅ −3 . Charge accumulation (blue) is seen near the sulfur atoms and reduction (purple) near the hydrogen atoms at the MoS 2 /water interface. d The laterally averaged DCD along the Z-direction. at 0.03 eV above the original VBM. As a result, the band edge shift of MoS 2 due to water is different for VBM and CBM, not only on the quantitative number but also microscopic origin, i.e., the fourth effect: chemical contact.
Four other promising 2D photocatalysts GaS, InSe, GaSe, and InS 6,20-22 are calculated, in order to infer the generality of the qualitative findings of the combined MoS 2 /water system. As plotted in the left panel of Supplementary Fig. 6, the contributions of the p z orbital of S/Se atoms are proportional to the size of the red circles since S/Se atoms locate at the interfacial sites closer to the water molecules when these 2D photocatalysts contact with water. Obviously, for their VBMs, significant projection from the p z orbital of S/Se atoms are observed while negligible contributions from these out-of-plane orbitals are displayed for their CBMs. These indicate that in the presence of water, similar to MoS 2 , there are emerging hybridized VBM states sitting above the original VBMs for GaS, InSe, GaSe, and InS, thus the shifts of their VBMs differ from those of their CBMs. Accordingly, quick verification calculations are performed by adsorbing one layer of water molecules on both surfaces of these four monolayer 2D photocatalysts and wavefunction distributions of VBM, CBM are plotted in the middle and right panel of Supplementary Fig. 6. As expected, the wavefunctions of their VBMs not only distribute on these 2D photocatalysts but also on water molecules, which exhibits chemical hybridization of the p z orbitals between the 2D materials and water. The bandgaps of these four promising 2D photocatalysts with and without water are summarized in Supplementary Table 1 (See Supplementary Note 4). We found that in the presence of water, all the bandgaps become smaller, which means that the shifts of their VBMs differ from their CBMs due to the existing chemical contact between 2D photocatalysts and water. Notably, the band gap of GaS is reduced by 0.11 eV even with only one layer of water molecules. We conclude that it is rather general that under aqueous condition, the VBMs of 2D materials whose surface atoms contain out-of-plane orbitals will shift differently from their CBMs, and the different shifts are caused by the fourth microscopic origin: the emerging hybridized states appear between the out-of-plane orbitals of 2D materials and water.

DISCUSSION
Because our DFT approach directly calculates the combined 2D material/water system self-consistently, it is able to capture the correct band edge shifts and the complete microscopic origins of these shifts for cases when CBM and VBM are shifted by different values due to the presence of water, e.g., MoS 2 /water, GaS/water, GaSe/water, InSe/water, and InS/water systems. The reported direct DFT scheme is systematic to quantitatively calculate band edge positions, band edge shifts and microscopic origins of the shifts, for a wide range of aqueous material systems where thousands of atoms or more are able to be involved. The findings of the band edge shifts is of great importance for artificial photosynthesis, e.g., MoS 2 is good for HER but not OER. The CBM of MoS 2 receives an averaged upshift of 0.42 eV due to water. The predicted upshift of CBM band edge in MoS 2 due to water should thermodynamically help photon activated electrons to participate in HER. Furthermore, in a recent experiment 19 of MoS 2 grown on GaN nanowires which sit on Si wafer, MoS 2 is found to provide excellent electronic pathways for photo-electrons to traverse toward water molecules absorbed on MoS 2 . These microscopic discoveries underpin the experimental findings that MoS 2 is an efficient photocatalyst for hydrogen evolution reaction. This upshift is due to three microscopic contributions: the geometric deformation effect contributes −0.04 eV, the water dipole effect contributes 0.21 eV, and the charge redistribution effect contributes 0.26 eV. Concerning VBM of MoS 2 , apart from these three contributions, an additional upshift is due to hybridiation of the pz orbitals in MoS 2 and in water, which generates a VBM state mimicking a chemical-like contact. Therefore, for any other semiconductors, especially those promising photocatalysts, whose surface atoms contain out-of-plane orbitals, such as GaS, InSe, GaSe, and InS shown in Supplementary Fig. 6, the interfacial chemical interaction plays an indispensable role in the band edge shifts for semiconductor/liquid systems. The chemical contact effect could be sensitive to the level of functionals in DFT, which is so far at the GGA level, and a higher level analysis such as using GW may be in order. And when microscopic physics of the water molecules and water/solid interface are not important, implicit solvent models via reaction fields can be a very powerful approach. The reaction field method is much more computationally efficient, it provides a long-time average of the solvent by a dielectric function and is not limited by the number of atoms one can handle. We believe combining the reaction fields with the large scale DFT/MD approach of this work into a multiscale approach, should be an interesting future direction.
In this work we focus on investigating the electronic structure of aqueous 2D photocatalysts in neutral environment, which provides benchmark results for future studies on pH effect. For example, in the acidic environment, the following formula could be employed to estimate the potential change induced by pH change 32 : Where Δϕ is the potential change,k B is Boltzmann's constant, T is the temperature and e is the elementary charge. α = 1 gives the Nernstian pH dependence. Therefore, for a decrease in pH of one unit amounts to a potential increase of 59 meV. In summary, we proposed a direct DFT approach to quantitatively evaluate band edge positions, band edge shifts and microscopic origins of these shifts for 2D materials in water, which is applicable for a wide range of 2D material systems. Taking monolayer MoS 2 as an example, Kohn-Sham DFT is applied to directly calculate the combined MoS 2 /water system containing over 1800 atoms. The direct calculation avoids using the rigid shift approach thus the CBM and VBM do not need to shift by a same amount. The CBM shift of the combined MoS 2 /water system can be well described by quantitative origins including the geometric deformation effect, the water dipole effect and the charge redistribution effect. For the different shift found in VBM, the chemical contact effect contributes a value of 0.03 eV while a value of 0.08 eV derived from the geometric deformation effect. Our findings on the property of aqueous photocatalyst interfaces with 2D materials provide useful information about low-cost and high-efficiency photocatalysts.

Molecular dynamic simulation
Since the photocatalyst is surrounded by water, the interaction between the semiconductor and water must be taken into account to predict the band edge shift due to water. The water molecules continue to adsorb and desorb on the solid surface, we account for this dynamical nature by calculating many time-evolution configurations. The supplementary Fig. 1 plots a snap-shot of an atomic configuration of MoS 2 /water system containing 1872 atoms, including 560 water molecules, 64 Mo atoms and 128 S atoms. These atoms are inside an orthorhombic simulation supercell of 21.91 × 25.29 × 36.88 Å 3 .
We determine the atomic configurations of the MoS 2 /water system by molecular dynamics (MD) using the LAMMPS package 33 . The calculations are carried out in two steps. First, in an isothermal-isobaric (NPT) ensemble with pressure 1 atm and temperature 298.15 K, the MoS 2 /water system is equilibrated until the density of water in the bulk region reaches its equilibrium density of 1 g/cm 3 . MD simulation of 1 ns is adequate for this process. Second, the resulting configuration is placed in a canonical (NVT) ensemble for a further 10 ns MD simulation, during which the configurations used in the electronic structure calculation are sampled. In these two steps, Nose-Hoover thermostat is used to maintain the system temperature at 298.15 K while the SHAKE algorithm used to keep the bonds rigid in the liquid water molecules. The force-field parameters in the LAMMPS MD are adopted from those correctly describing the cohesive interactions in MoS 2 as well as between MoS 2 and water 34 . Equation (7) is the energy expression in the force field, where U bond , U angle , U lj , and U c represents bond potential, angular potential, Lennard-Jones potential and Coulomb potential, respectively. k r and k θ are the harmonic bond constant and angle constant respectively, r 0 and θ 0 are the equilibrium bond length and angle. r ij is the distance between the ith and jth atoms. θ ijk is the angle between the ith, jth, and kth atoms with the jth atom at the top corner. ϵ is the interaction strength between two atoms and σ the distance at which the potential energy between two atoms vanishes. ϵ 0 is the dielectric constant of free space. q i and q j are the partial charges on the ith and jth atoms. The harmonic bond, angle force constants and Lennard-Jones parameters are obtained by fitting to the experimental lattice parameters and mechanical properties 34 . The Lennard-Jones parameters describing the interaction between water and MoS 2 is derived from the geometricmean combining rules, i.e., ϵ ij ¼ ffiffiffiffiffiffi ffi where the interaction cutoff is set to 1.2 nm. The long-range electrostatic interactions are calculated with the Particle-Particle-Particle-Mesh (PPPM) method 35 where the partial charges for Mo and S atoms are obtained from the DDAP formalism 36 . For water molecules, the SPC/E force field 37 is adopted. The equilibrium bond length and angle are optimized based on the relaxed geometric structure of MoS 2 obtained by DFT. These parameters are summarized in Table 1. DFT-MD calculations are also performed to confirm the accuracy of these adopted parameters to describe the MoS 2 /water interface structure (see Supplementary Note 5). The DFT-MD results show that the overall structure derived from the force field is satisfactory and can serve as the starting point of further electronic structure calculations.

DFT computations
Having determined the atomic configurations of MoS 2 /water by MD, their electronic structure is calculated by our recently developed real-space KS-DFT method RESCU 16 . With advanced computational mathematics and high parallelization efficiency, RESCU is capable to perform accurate KS-DFT self-consistent calculations for over ten thousand atoms, and we refer interested readers to the original work for technical details of this method 16 . In this work, due to the large system size, a single Γ-point ksampling is adequate for both total energy and band energy calculations, which is verified by convergence tests. The atomic cores are defined by the optimized norm-conserving Vanderbilt (ONCV) pseudopotential 38 . A double-zeta polarized (DZP) atomic orbital basis is used for H and O elements. For Mo element, the basis set contains single-zeta 4s/4p, doublezeta 4d/5p, triple-zeta 5s, and one orbital with f-symmetry. For S element, the basis set contains double-zeta 3s/3p, one orbital with d-symmetry, and one orbital with f-symmetry. The real-space grid resolution is set at 0.3 Bohr. The exchange-correlation is treated at the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) level. The calculated GGA band gap of prinstine monolayer MoS 2 is 1.75 eV, which agrees with known GGA values 39 . GW calculations, which takes into acount of electron many body interaction, indicate the quasi-particle energy gap of monolayer MoS 2 to be about 2.8 eV 40 . However, the experimental photoluminescence (PL) spectra gives an excitonic gap of 1.8-1.9 eV 41 . The experimental results tend to prove that the exciton binding energy of about 1.0 eV cancels out the electron many body interaction in monolayer MoS 2 . This makes our GGA gap of 1.75 eV reasonable for monolayer MoS 2 .
In the present work, we perform direct DFT calculations of electronic structure of the combined MoS 2 /water system where thousands or more atoms are necessary. To this end, our calculations are performed at the GGA level. Our KS-DFT calculations are carried out on a computer cluster of 4-nodes with a total of 64 computing cores each having 4-GB memory. RESCU took about 2 wall-clock hours to converge one electronic structure calculation for the MoS 2 /water system of 1872 atoms. To calculate the averaged band edge shift of MoS 2 in the presence of water accurately, 60 configurations in the MD trajectory of 0.1-6.0 ns with a time step of 0.1 ns are sampled to calculated the CBM and VBM of the combined MoS 2 /water systems. The averaged CBM and VBM values versus the number of sampled configurations are plotted in Supplementary Fig. 8, the results converge by sampling beyond 30 configurations, and the variation is about 0.01 eV in the sampling number interval [40,60], which can be taken as the error bar of the band edges.

DATA AVAILABILITY
The data supporting the findings of this work are available from the corresponding author on reasonable request.