Insights from density functional theory calculations on heteroatom P-doped ZnIn2S4 bilayer nanosheets with atomic-level charge steering for photocatalytic water splitting

ZnIn2S4 (ZIS) is an efficient photocatalyst for solar hydrogen (H2) generation from water splitting owing to its suitable band gap, excellent photocatalytic behaviour and high stability. Nevertheless, modifications are still necessary to further enhance the photocatalytic performance of ZIS for practical applications. This has led to our interest in exploring phosphorus doping on ZIS for photocatalytic water splitting, which has not been studied till date. Herein, phosphorus-doped ZnIn2S4 (P-ZIS) was modelled via Density Functional Theory to investigate the effects of doping phosphorus on the structural and electronics properties of ZIS as well as its performance toward photocatalytic water splitting. This work revealed that the replacement of S3 atom by substitutional phosphorus gave rise to the most stable P-ZIS structure. In addition, P-ZIS was observed to experience a reduction in band gap energy, an upshift of valence band maximum (VBM), an increase in electron density near VBM and a reduction of H* adsorption–desorption barrier, all of which are essential for the enhancement of the hydrogen evolution reaction. In overall, detailed theoretical analysis carried out in this work could provide critical insights towards the development of P-ZIS-based photocatalysts for efficient H2 generation via solar water splitting.

www.nature.com/scientificreports/ photocatalytic performances. Among all types of photocatalysts, two-dimensional (2D) bilayer zinc indium sulfide (ZnIn 2 S 4 or ZIS) has attracted much research interests due to its high chemical stability, suitable band gap energy, facile synthesis, ease of modification and excellent photocatalytic activity. By adopting appropriate modification strategies, ZIS has the potential to provide superior photocatalytic efficiencies for practical HER applications in the energy sector 21 . ZIS, a ternary metal chalcogenide, is a layered structure semiconductor with a direct band gap that could exist in three different crystal polymorphs (cubic, trigonal, and hexagonal) 22 . Hexagonal ZIS has been reported to be the most thermodynamically stable structure with high photoactivity and is therefore given more focus in photocatalytic applications [23][24][25][26][27] . It exhibits a band gap energy (E g ) of ~ 2.40 eV, with a conduction band minimum (CBM) of − 0.85 V (more negative than HER potential of 0 V) and a valence band maximum (VBM) of + 1.55 V (more positive than OER potential of + 1.23 V) against NHE potential at pH = 0 25 . Hence, ZIS is capable in utilizing a broad range of solar incident ray for the photoexcitation of electrons to drive photocatalytic water splitting. Despite the aforementioned merits, ZIS still suffers from several drawbacks such as high electron-hole recombination rate, sluggish diffusion of charge carriers to active sites and low exposure of active (110) facets for HER 21,[27][28][29][30] . Such limitations inhibit the photo-redox reaction of water and deteriorate photocatalytic performances. Therefore, appropriate modifications must be done to improve intrinsic electronic properties, increase charge separation efficiency and promote the migration of charge carriers to active sites for initializing the redox reactions. Elemental doping, either at substitutional or interstitial sites, involves the introduction of external heteroatoms into the crystal lattice of semiconductor. This addition could bring about new intermediate levels that reduce the band gap energy of the photocatalyst, thereby improving its photo-responsiveness towards solar spectrum 22,31 . Additionally, the presence of heteroatom dopants could alter the material's electronic structure, which could in turn promoting the separation of charge carriers and their migration rates [32][33][34] . Phosphorus doping (P-doping) has been reported to show significant improvements on numerous semiconductors by accelerating the separation and transmission of photoinduced electron-hole pairs [35][36][37][38][39][40][41][42] . However, to the best of our knowledge, no prior research has been conducted to investigate the intrinsic nature and effects of P-doping on ZIS. Hence, in this work, Density Functional Theory (DFT) calculations were performed to investigate and elucidate the effects of P-doping on ZIS. The doping nature as well as the structural and electronic properties of P-ZIS, including its charge density distribution, band structure and density of state (DOS) were studied in detail. Systematic investigations were also carried out to examine the performance of P-ZIS towards water interaction, HER and OER.

Computation detail and methodology
The theoretical structure was optimized using Vienna Ab initio Simulation Package (VASP) with projector augmented wave (PAW) method and created reciprocal space with plane-wave basis. DFT calculation was implemented in VASP, via generalization-gradient approximation (GGA) with exchange-correlation function of Perdew-Burke-Ernzerhof (PBE) 41 . Basic plane-wave settings were based on previously reported ZIS modelled structure in order to obtain consistent, comparable and convergence results. For instance, the energy cut-off was set at 500 eV. The position of all atoms and cell shapes were allowed to relax until an energy convergence of 1 × 10 −5 eV and force convergence of 0.01 eV/Å were obtained. During initialization calculation, 1-by-1 atomic bilayer with 2 Zn atoms, 4 In atoms and 8 S atoms was adapted. The Monkhorst-Pack k-point mesh was set at 3 × 3 × 1. Additional vacuum layer of approximately 15 Å was applied to eliminate possible interaction between periodic images 43 . Hybrid functional HSE06 was used to obtain the most stable structure during the relaxation calculation. For interstitial P-doping calculation, single P atom was added into possible interstitial sites of 1-by-1 atomic bilayer structure. Each intrinsic S atom was replaced by P atom during the substitutional P-doping calculation. In the calculation of the bandgap and DOS, high symmetry k-points of interest for ZIS structure were used including Ŵ, M, K and ŴA . All model lattice structures were visualized using Visualization for Electronic and Structural Analysis (VESTA).
For surface reaction study involving surface adsorption and active site interaction, single layer structure was utilized for investigation as recommended 28,41,44,45 . Herein, 1-by-1 atomic single layer ZIS (1 Zn atom, 2 In atoms and 4 S atoms) as well as P-ZIS (1 Zn atom, 2 In atoms, 3 S atoms and 1 P atom) structures were used for ease of convergence during surface reaction studies. 1.5 Å separation between active sites and adsorbed species were adopted for initialization calculation 42,46 . A number of active sites used in this study was chosen based on other reports 29 . During surface adsorption and interaction studies, additional Van der Waals (vdW) correction term utilizing D2 method of Grimme (DFT-D2) was accounted in computation of lattice relaxation to obtain accurate calculations in interatomic forces, potential energy and stress tensor 47 . Following that, free energy values calculated in DFT were utilized for analysis on both ZIS and P-ZIS structures towards HER, OER and water interactions, whereby adsorption energies and the corresponding Gibb's free energies were calculated to conclude potential performance improvement (see Supplementary Information online for computation and calculation details).

Results and discussion
Pristine structure optimization and phosphorus doping simulation. To ensure that the correct structure was formed, several key parameters were cross-checked with literature values including basic lattice parameters a = b = 3.85Å, c = 24.68Å, α = β = 90 • , γ = 120 • , bilayer thickness (~ 2.468 nm) and intermediate layer gap (~ 0.4 nm) 41,48 . Figures 1 and S1 show the simulated structures with different facets, which are consistent with those reported in the literature 41,48,49 . It should be noted that even though the simulated bilayer ZIS structure exhibited an identical arrangement of atoms at the top and bottom layers, the bond lengths of Zn1'-S1' (top layer) and Zn1-S1 (bottom layer) differed due to potential interactions across the interlayer spacing. A similar argument can be made for the other bonds present. To confirm that the simulated ZIS structure was fully DFT calculations were then conducted to analyse two different doping behaviours: (1) interstitial P-doping (IPD); and (2) substitutional P-doping (SPD), where the S atom was replaced due to its similar electronegativity and atomic radius of ~ 100 nm 41,42 . Figure 1c illustrates eight of the possible P-doping sites on ZIS. Table 1 summarizes the formation energy ( E f ) of the eight possible P-doping sites (see Supplementary Information online for calculation details). It should be noted that the lowest E f calculated corresponded to the most stable and energetically favourable structure 41,50 .
Based on Table 1, the formation energy calculated for each of the eight cases had a negative value, which inferred that all of the IPD and SPD processes were energetically favourable. Possible IPD scenarios were considered at 4 different sites with interstitial spaces larger than the P atom. After lattice relaxation computation, significant distortions were observed across all IPD structures (see Fig. S3). As shown in Table 1, the calculated formation energies for IPD structures (ranging from − 0.991 to − 1.172 eV) were much higher than those of SPD structures, which suggested that IPD generally created less favourable structures. On the other hand, all relaxed SPD structures experienced little to no lattice distortion due to the similar atomic radii (~ 100 nm) between S and P atoms (see Fig. S3). Consequently, the calculated formation energies for SPDs (ranging from − 2.011 to − 2.082 eV) were comparatively lower than IPDs, with SPD-3 exhibiting the lowest value. In short, it could be concluded that the substitutional doping of P-atom (replacing intrinsic S3 atom) on ZIS gave rise to the most stable structure, as evidenced by its lowest formation energy of − 2.082 eV. It was believed that the breaking of intrinsic In-S bonds and the formation of new In-P bonds around the In atoms in ZIS were favoured. Thus, all subsequent analyses and studies were performed on the SPD-3 structure. Figures 2, 3 and 4 show the comparison of lattice structures, DOS and charged density distributions between pristine ZIS and P-ZIS (SPD-3) to unravel the potential effects of P-doping towards photocatalytic water splitting.
Comparing the relaxed structures of ZIS and P-ZIS (SPD-3) in Fig. 2, no significant lattice distortions were found. However, a larger interatomic spacing was observed when P atoms substituted the intrinsic S atoms at the bottom layer. This was attributed to the replacement of smaller S atoms with larger P atoms in the P-ZIS structure. The observation also coincided with the negative shift in (110) and (120) peaks in the powder diffraction pattern of P-ZIS, which indicated an increment in lattice spacing as compared to that of pristine ZIS (see   Table 1.
Summary of E f for each interstitial and substitutional doping site.
Doping Sites www.nature.com/scientificreports/ Fig. S2). Contrastingly, the interatomic spacing around S1 and S1' atoms in P-ZIS were found to be reduced, with a slight positive shift of the (006) facet. Furthermore, the lattice bilayer thickness and interlayer spacing were also reduced to 2.477 nm (− 0.36%) and 0.396 nm (− 1.25%), respectively. Apart from those, all other characteristic peaks of P-ZIS were well-aligned to the pristine ZIS structure. The difference in the relative intensity between both patterns could be ascribed to the additional stress and strain imposed by the doping of P atoms 28 . The calculated DOS and band structure for both pristine ZIS and P-ZIS (SPD-3) are shown in Fig. 3. The band gap energy of pristine ZIS was found to be approximately 0.404 eV (see Fig. 3a, c). Although hybrid function HSE06 was used for band gap and DOS calculations, the results were clearly underestimated as the experimental band gap energy of ZIS should be roughly 2.40 eV 25 . It is known that using DFT to calculate the band gap of complex structures tends to have lower precision. Similar observations were also reported in other papers, where relatively smaller band gaps for ZIS (ranging from ~ 0 to 0.916 eV) were found using DFT simulation 29 48 , such spatial separation of electronic states near the VBM and CBM could potentially give rise to errors and an underestimation of band gap energies. Figure 3b shows the calculated DOS of P-ZIS, where a variation of band edge was detected following P-doping. It can be observed that the VBM upshifted by 0.242 eV, following by a significant increase in electron density near the VBM associated to the additional P-3p orbital (see orange shading in Fig. 3b). As the CBM remained unchanged, P-ZIS exhibited an overall reduction in band gap energy to 0.162 eV (− 59.9%). This finding was consistent with other DFT studies conducted over anion-doped (O and N) ZIS, where only an upshift of VBM was observed with no changes to the CBM and without the presence of any midgap state 31,51 . Additionally, P-ZIS exhibited different band structures from pristine ZIS (see Fig. 3c, d) owing to the contribution from the atomic orbitals 28 , which further demonstrates that P-doping could alter the electronic properties of ZIS. In short, the reduction in band gap energy could improve the photo-responsiveness of P-ZIS, by harvesting a broader range of the solar spectrum for photocatalytic applications. Moreover, the increase in electron density near the VBM indicates that more ground-state electrons would be available for photoexcitation and subsequent participation in HER upon light irradiation. Simultaneously, the additional states near the VBM close to the Fermi level could enhance the metallic conductivity of P-ZIS, thereby facilitating the mobility of holes and further inhibiting the recombination of photoexcited charge carriers for augmented performance in HER 51 . Figure 4 reveals the charge density distribution around the atoms in ZIS and P-ZIS at the doping layer. It is clear that the S3 atom possessed high charge distribution, which was replaced by P atom via substitutional doping. Because the electronic states near the CBM were primarily contributed by the S-3p orbital, such substitution step resulted in the decrease of DOS near the CBM of P-ZIS. Similarly, the introduction of P atom into ZIS contributed to higher charge densities near the VBM in the calculated DOS. It is evident that after P-doping, the charge density distribution around S atoms (potential HER active sites) increased significantly due to the delocalization of charges. This would potentially improve the interaction between hydrogen atoms with potential active sites, which would consequently enhance hydrogen adsorption and the photocatalytic performance towards HER.
Water interaction study. During photocatalytic water splitting, the adsorption of water molecules onto the photocatalyst surface serves as an important step for HER and OER. This is critical for providing H + atoms for the reductive generation of H 2 molecules, while simultaneously driving the formation of surface hydroxyl radicals to generate O 2 molecules. Therefore, the interaction and adsorption behaviour of water molecules on pristine ZIS and P-ZIS structures were analysed via DFT. The bonding strength of water adsorption energy ( E H 2 O * ) was calculated (see Supplementary Information online for calculation details). By definition, a more negative value of E H 2 O * would imply a more exothermic adsorption process, and thus, a stronger interaction between the H 2 O molecule and the photocatalyst surface 43 . In addition, the effect of P-doping on H 2 O interaction was also analysed on the basis of the changes in bond angle and bond length of the adsorbed H 2 O molecule.
All possible sites (Zn, In, S and P) for H 2 O adsorption in both pristine ZIS and P-ZIS structures were investigated. Based on a study by Li et al. 45 , H 2 O molecules have the tendency to be adsorbed onto positively charged atoms or atoms with lower localized electron density. As all S atoms are non-positively charged, there was no adsorption of H 2 O molecule on the S atoms observed, irrespective of the initial position of the H 2 O molecule. Similarly, no H 2 O molecule was found to be adsorbed onto the Zn atoms on both pristine ZIS and P-ZIS. This could be ascribed to the high localization electron cloud surrounding the Zn atom (see Fig. 5), which induced repulsion with the free H 2 O molecule and prevented its adsorption. These findings indicated that neither S nor Zn atom could serve as active sites for the adsorption of H 2 O. On the other hand, the adsorption and activation of H 2 O molecule were observed on In2 atoms in both pristine ZIS and P-ZIS, as well as P-atoms in P-ZIS (see Fig. 6). This suggested that both In2 and P could serve as activation regions for H 2 O and potentially act as reactive sites for OER. Interestingly, no H 2 O molecule was adsorbed onto In1, which could be due to insufficient interatomic     . Close-to-zero value of ∆G H * implies that reaction barriers in both adsorption and desorption steps are compromised which favours HER, serving as an indicator for good photocatalyst for HER 43 .
To evaluate the hydrogen adsorption behaviours and potential HER performances, the surface properties of monolayer pristine ZIS and P-ZIS were analysed in detail. The analysis was carried out to elucidate the effect of P-doping on the hydrogen adsorption behaviour over the least favourable site (S1), the most favourable site (S2) as well as the doping site (S3 for ZIS and P for P-ZIS). To facilitate the visualization of hydrogen adsorption, 6 different relaxed structures and their corresponding hydrogen interactive bond lengths obtained from DFT are depicted in Fig. 6. The free energy diagram for HER for each site on ZIS and P-ZIS is displayed in Fig. 7.
Based on the calculation conducted for pristine ZIS (Fig. 7), it can be confirmed that S1 was the least favourable HER active site (high G H * of 1.799 eV) while S2 was the most favourable HER active site (low G H * of 1.342 eV). After performing P-doping, the charge density distributions around the least (S1) and most favourable (S2) active sites were observed to have drastically improved (Fig. 4), which could collectively enhance the interaction and adsorption of H* for HER to take place 55 . From Figs. 6 and 7, the G H * on S1 was found to have reduced from 1.799 (for pristine ZIS) to − 0.324 eV (for P-ZIS). The latter, with a value much closer to zero, indicated that the adsorption-desorption barrier was significantly reduced, which rendered P-ZIS more favourable for H* adsorption. This affirmed that P-doping could result in the improvement of HER performance on the least favourable S1 site. Similarly, the increment of charge density distribution after P-doping on S2 further improved the favourability of H* adsorption onto the most favourable active site (as indicated by the close-to-zero G H * ). The significant reduction of the adsorption-desorption barrier for P-ZIS was attributed to the elongated S-H interactive bond, which in turn eased the desorption of H* to form H 2 . Likewise, the substitution of P atom favoured the adsorption of H* as indicated by its closer-to-zero G H * of − 0.596 eV as compared to the S3 atom (1.399 eV). The relatively longer P-H bond over S3-H bond could further promote the desorption of H* for H 2 formation and subsequently improve the HER activity. In short, it is clear that P-doping in ZIS enhances HER performance by (1) reducing H* adsorption-desorption barriers for the generation of H 2 , (2) increasing electron density near HER active sites and (3) inhibiting electron-hole pairs recombination as previously discussed.
Oxygen evolution reaction and adsorption study. In the overall photocatalytic water splitting process, OER is considerably more complicated than HER as it involves a four-electron transfer pathway. Firstly, H 2 O molecule is adsorbed onto the surface active sites, which is then followed by the formation of 3 different oxygenated reaction intermediates (HO*, O* and HOO*). Subsequently, molecular O 2 is formed from the HOO* species, in which the reaction occurs at coordinatively unsaturated surface active sites. In short, OER consists of 4 elementary reaction steps (OER1, OER2, OER3 and OER4) as listed in Eqs. (2)(3)(4)(5) 43 . For comparison, the overpotential η OER required to achieve a downhill trend for all OER free-energy steps at standard equilibrium potential (where external bias U = 1.23 V 56 ) was evaluated using Eq. (6). To be classified as a good OER catalyst, the calculated η OER should be as low as possible 43 .
Herein OER analysis was carried out for both pristine ZIS and P-ZIS to study the effect of P-doping on OER performance in In2 active site. The different reaction intermediate steps are shown in Fig. 8. The free energy diagram for OER pathway is also reflected in Fig. 9, with each coordinate Gibb's free energy of reaction and ∆G OER of corresponding pathway tabulated in Tables S1 and S2 respectively.
Focusing on pristine ZIS, it was previously mentioned that the In2-O bond for adsorbed H 2 O was 2.280 Å (Fig. 5a). When the O-H bond was firstly cleaved via deprotonation of H atom to form HO*, the In2-O bond length was found to decrease to 2.264 Å (Fig. 8a). The further cleavage of the 0.960 Å O-H bond via subsequent H deprotonation to form O* led to an additional reduction in the In2-O bond length to 2.227 Å (Fig. 8b). When O* reacted with free H 2 O molecule to form OOH*, the In2-O bond length decreased to 2.214 Å (Fig. 8c). The OER on pristine ZIS was completed following the cleavage of the 0.987 Å O-H bond via the last H deprotonation step. This was followed by the desorption of O 2 from the ZIS surface, which resulted in the production of a free O 2 molecule.
The initial In2-O bond for adsorbed H 2 O on the P-ZIS structure was previously found to be 2.312 Å (see Fig. 5b). After the first cleavage of the O-H bond to form HO*, the In2-O bond length decreased to 2.281 Å as shown in Fig. 8d. For P-ZIS, the subsequent deprotonation of H atom to form O* was deduced to be relatively easier compared to ZIS due to the longer O-H bond of 0.965 Å as shown in Fig. 8e. Next, O* on P-ZIS reacted with the free H 2 O molecule to form OOH* as shown in Fig. 8f. In order to produce free O 2 molecule, the O-H bond must be cleaved via deprotonation, followed by the desorption of O 2 from P-ZIS. Similarly, as both of the O-H (0.992 Å) and In2-O (2.221 Å) bond lengths for P-ZIS were comparatively longer than those in ZIS, easier H deprotonation and O 2 desorption were therefore expected over P-ZIS. These could collectively lead to improved OER performance of P-ZIS as compared to its pristine counterpart.
To further unravel the mechanism of OER, the free energy pathway diagrams for the process over pristine ZIS and P-ZIS were constructed and analysed. For pristine ZIS as presented in Fig. 9a, it can be seen that OER2 and OER4 (as in Eq. (3) and (5) respectively) presented an uphill process with no applied bias (U = 0 V). At the standard equilibrium potential for OER (U = 1.23 V), OER4 remained uphill with a η OER value of 5.147 V. This implied that the desorption of O 2 from the photocatalyst surface acted as the rate determining step (RDS). As shown in Fig. 8, the In2-O bond for HOO* was the shortest among all intermediate steps, which led to its strong adsorption onto In2. To achieve a downhill trend for each step, a bias U of 6.38 V was required. The relatively high η OER suggested that OER over pristine ZIS was energetically unfavourable. Figure 9b depicts the free energy diagram for OER over P-ZIS. It can be seen that OER3 and OER4 (as in Eqs. (4) and (5) respectively) presented an uphill process in the absence of an applied bias U. Similar to pristine ZIS, OER4 was observed to be the RDS www.nature.com/scientificreports/ at the standard equilibrium potential for OER (U = 1.23 V) with a η OER value of 5.123 V. It is noteworthy that even though P-doping reduced the overpotential by 0.024 V and out-performed pristine ZIS in the OER intermediate steps, the high η OER (> 1 V) indicated that P-ZIS was not energetically efficient in promoting OER for O 2 production via overall water splitting 56 . Owing to the infeasibility of OER4's uphill process, the HOO* radical is anticipated to undergo oxidative recombination with H + to produce H 2 O 2 instead 59,60 . In short, P-ZIS has the potential to efficiently drive solar water splitting to generate H 2 ; however, it does not result in the generation of the desired O 2 product.

Conclusion
In this work, detailed DFT calculations and computation analysis were successfully performed to fully unravel the nature of P-doping on ZIS. The most stable P-doped ZIS structure was obtained via the replacement of S3 atoms with substitutional doping of P.
The calculated E f of − 2.082 eV suggested that the formation of P-ZIS was energetically favourable. A few key points were observed following the doping of P atoms: (1) The band gap of P-ZIS was reduced, indicating enhanced photo-responsiveness; (2) The VBM was upshifted close to the Fermi level due to P-3p orbital contributions, which led to enhanced hole mobility and charge carrier separation; (3) The electron density near the VBM of P-ZIS was increased, which provided more ground-state electrons for photoexcitation and subsequent participation in HER; (4) The H* adsorption was drastically improved due to the increase in charge density distribution around the most favourable HER active site (S2); this resulted in the significant reduction in the H* adsorption-desorption barrier, which could in turn enhance HER performance. The P-ZIS photocatalyst not only displayed improved interaction on the intrinsic OER active site (In2), but also introduced P atoms as new reactive sites for H 2 O adsorption. It should be noted that although P-ZIS exhibited a lower overpotential for OER compared to ZIS, its high η OER value of 5.123 V indicated that OER was still energetically unfavourable over P-ZIS. All in all, we strongly believe that this work would provide critical insights into development of high performing P-ZIS-based photocatalyst for enhanced H 2 generation, which could ultimately bring about the successful commercialization of solar-driven water splitting in the near future. www.nature.com/scientificreports/