The Construction of Surface-Frustrated Lewis Pair Sites to Improve the Nitrogen Reduction Catalytic Activity of In2O3

The construction of a surface-frustrated Lewis pairs (SFLPs) structure is expected to break the single electronic state restriction of catalytic centers of P-region element materials, due to the existence of acid-base and basic active canters without mutual quenching in the SFLPs system. Herein, we have constructed eight possible SFLPS structures on the In2O3 (110) surface by doping non-metallic elements and investigated their performance as electrocatalytic nitrogen reduction catalysts using density functional theory (DFT) calculations. The results show that P atom doping (P@In2O3) can effectively construct the structure of SFLPs, and the doped P atom and In atom near the vacancy act as Lewis base and acid, respectively. The P@In2O3 catalyst can effectively activate N2 molecules through the enzymatic mechanism with a limiting potential of −0.28 eV and can effectively suppress the hydrogen evolution reaction (HER). Electronic structure analysis also confirmed that the SFLPs site can efficiently capture N2 molecules and activate N≡N bonds through a unique “donation-acceptance” mechanism.


Introduction
Ammonia (NH 3 ) is not just an irreplaceable nitrogen-containing chemical raw material in the synthesis of traditional chemical products such as fertilizers, nitric acid, and explosives.In recent years, ammonia fuel has emerged as a clean energy source of worldwide interest due to its ease of transport and storage and its lack of carbon emissions in operation [1][2][3].At present, large-scale industrial ammonia synthesis is still based on the Haber-Bosch (H-B) reaction, which requires high temperatures (400-500 • C) and high pressures (100-200 bar), which not only results in extremely high energy consumption but also in extremely stringent equipment requirements.Many researchers are searching for new and clean ammonia synthesis technology to replace the traditional industrial ammonia synthesis technology in order to alleviate the problems of high resource consumption and severe environmental pollution caused by industrial ammonia synthesis [4].Among them, electrochemical nitrogen reduction reaction (NRR) technology is a method with great application potential and research value, which has outstanding advantages such as high efficiency, low energy consumption, and simple reaction devices [5][6][7].Electrochemical nitrogen reduction reaction technology is based on the synthesis of ammonia (or NH 4 + ) from N 2 and water (H 2 O or H + ) at ambient temperature and pressure, driven by electricity with the applied voltage.The conversion and storage of intermittent energy is facilitated when the electricity consumed comes from clean energy sources such as solar, wind, and tidal power.
The development of efficient catalysts is a central task in the commercialization of electrocatalytic nitrogen reduction ammonia technology.It is well known that transition metals (TMs) with both d-orbital electrons and unoccupied d-orbital electronic structures Molecules 2023, 28, 7130 2 of 12 can form N 2 -M σ-donation and M-N 2 π-bonding configurations with nitrogen molecules, which can lead to the activation of nitrogen via the π-bonding pathway [6].Therefore, most of the electrocatalysts reported so far contain transition metals such as Fe, Mo, W, etc.However, HER is more likely to occur in most transition metal-based electrocatalysts, leading to problems such as the poor selectivity of the ammonia synthesis reaction on the catalysts and low Faradaic efficiency.
Compared to transition metal catalysts, elemental materials in the p-block are inherently weak in hydrogen precipitation activity, which is feasible for selective electrocatalytic ammonia synthesis.Légaré et al. [8] found that the Lewis acid-containing elemental boron in the boranylidene group can also inject electrons into the N 2 molecule via back-donation due to the simultaneous presence of empty sp 2 orbitals and occupied p-orbitals and πantibonding orbitals to achieve N 2 activation.By doping porous carbon with elemental fluorine with higher electronegativity, Liu et al. designed and prepared highly active (ammonia yield: 197.7 µg h −1 mg cat −1 ) p-block elemental-based electrocatalysts (F-doped carbon) enriched with Lewis acid sites and programmed elevated temperature desorption in the nitrogen atmosphere [9].The N 2 desorption peak of the F-doped carbon was found to be located at 436 • C, which is 115 • C higher than that of the undoped F sample, indicating that the Lewis acid sites can give the catalyst a stronger nitrogen binding capacity.Hu et al. improved the NRR performance of the catalyst in neutral electrolyte by simultaneously modifying oxygen vacancies and hydroxyl groups on the surface of Bi 4 O 5 I 2 to induce the creation of unoccupied orbitals around the Bi atoms, which drastically lowered the protonation barriers of N 2 molecules [10].The above studies show that the presence of either acidic or basic Lewis sites can, to some extent, improve the ability of the catalyst to interact with N 2 .
The SFLPs system has acidic and basic Lewis active centers that are not quenched with each other, and the construction of this surface structure is expected to break the limitation of the single electronic state of the catalytic center of the p-block elemental materials, which can not only target the inert molecules through the synergistic effect of Lewis acids and bases, but also activate the adsorbed molecules efficiently through the unique "donationacceptance" process [11].Lin et al. took advantage of the existence of empty p-orbitals after the sp 2 hybridization of B atoms and paired electrons after the sp 3 hybridization of N atoms [12].The construction of local environments by spatially hindered B-N atom pairs can play a "pulling" role for N 2 molecules, and inertness can be achieved by weakening the N≡N bond molecular activation by weakening the N≡N bond.This work provides strong support for the construction of SFLPs sites on the surface of p-block element-based electrocatalysts and the elucidation of the mechanism of nitrogen reduction.
Indium-based oxides are a class of p-block elemental materials with simple compositions that are widely used in the catalytic reduction of CO 2 at room temperature to synthesize CO, CH 3 OH, and other C1 chemicals [13][14][15][16].The surface atomic arrangement of indium-based oxides is easily tunable, and SFLPs sites can be constructed using simple modifications such as the introduction of oxygen vacancies.Moreover, compared with BCN 2D nanosheets, the SFLPs sites of indium-based oxides can be constructed around the oxygen vacancies on the surface of the material, which also helps to construct high-density and high-activity surface SFLPs sites [17].
In this work, defective In 2 O 3 with one oxygen vacancy (V-In 2 O 3 ) was selected as the research object, and the construction of oxygen defects and heteroatom doping means were used to fine-tune the design of the spatial configuration of the Lewis acid/base center of the SFLPs sites on the surface of the material, the electronic structure, to explore the form of acid/base site combinations in the SFLPs (In/doped elements), and to deeply and systematically investigate the influence of the electronic structure of the doped elements, the spatial site-resistance configuration, and other factors on the FLPs sites on the activation of nitrogen molecules.At the same time, the influence of the constructed SFLPs sites on the nitrogen reduction reaction pathway was analyzed, which provided experimental and theoretical guidance for the preparation of p-region element-based electrocatalysts with excellent performance in nitrogen reduction for ammonia synthesis.

Results and Discussion
The surface of In 2 O 3 (110) is one of the major exposed crystal planes of c-In 2 O 3 observed in experiments [18][19][20].It has also been shown to be thermodynamically stable under electrocatalytic conditions in theoretical studies [21].Therefore, we chose the surface of In 2 O 3 (110) for our study.The top layer structure of the optimized perfect In 2 O 3 (110) surface is shown in Figure 1a, and consists of 8 In atoms and 12 oxygen atoms forming a repeating unit.In addition, previous reports have shown that oxygen vacancies (V O ) on the surface of In 2 O 3 can generate SFLPs sites, which serve as active sites for reactions [22][23][24].At the same time, Qin et al. calculated the vacancy formation energies of different oxygen vacancies using density functional theory, and the calculations showed that the V O4 site had the highest formation energy [25].
the spatial site-resistance configuration, and other factors on the FLPs sites on the activation of nitrogen molecules.At the same time, the influence of the constructed SFLPs sites on the nitrogen reduction reaction pathway was analyzed, which provided experimental and theoretical guidance for the preparation of p-region element-based electrocatalysts with excellent performance in nitrogen reduction for ammonia synthesis.

Geometric Structure of Doping V-In2O3
The surface of In2O3 (110) is one of the major exposed crystal planes of c-In2O3 observed in experiments [18][19][20].It has also been shown to be thermodynamically stable under electrocatalytic conditions in theoretical studies [21].Therefore, we chose the surface of In2O3 (110) for our study.The top layer structure of the optimized perfect In2O3 (110) surface is shown in Figure 1a, and consists of 8 In atoms and 12 oxygen atoms forming a repeating unit.In addition, previous reports have shown that oxygen vacancies (VO) on the surface of In2O3 can generate SFLPs sites, which serve as active sites for reactions [22][23][24].At the same time, Qin et al. calculated the vacancy formation energies of different oxygen vacancies using density functional theory, and the calculations showed that the VO4 site had the highest formation energy [25].Indium oxide is a kind of P region elemental material with simple composition which is often used in the catalytic reduction of CO 2 at room temperature to synthesize C1 chemicals such as CO and CH 3 OH.The surface atomic structure of indium-based oxides is easily modified, such as by heteroatom doping to build SFLPs sites.Meanwhile, as far as we know, there have been no investigations on the construction of SFLPs by doping nonmetallic elements on the surface of In 2 O 3 .To further enhance the catalytic ability of the SFLPs sites and avoid the shielding effect of oxygen vacancies in different electronic states, we modified the surface by replacing O 5 with various nonmetal elements (named n-M@In 2 O 3 , nonmetal = C, N, F, Si, P, S, Cl,) and screened for efficient and stable acid-base site combinations in SFLPs.The optimized structure of several nonmetal elements doped in In 2 O 3 (110) is shown in Figure S1.The results of the geometry optimization showed that the different non-metallic elements selected could be used to replace the O 5 atom for the modification of the SFLPs site.Subsequently, the stability of the constructed SFLPs site at the actual applied temperature was further tested using 10 ps AIMD simulation for all doped structures.As shown in Figure 2, the curves of the total energy change of the different structures were calculated as the simulation time increased.For C@In 2 O 3 and F@In 2 O 3 , the total energy of the system became significantly smaller within the first 1ps of the simulation.Correspondingly, the structure of the surface SFLPs sites changed.The C atom was oxidized to CO and then left the surface, while the F atom moved to the O vacancy and combined with In 3 and In 4 atoms, also changing the surface SFLPs site structure.Furthermore, although the total energy of B@In 2 O 3 and Si@In 2 O 3 did not change significantly during the simulation, the SFLPs site structure changed.The doped B atoms bound to the oxygen atoms on the surface to break the SFLPs site structure, while the Si atoms filled the original O 4 vacancy defect, as well as the doped F atoms.The diagram of the radial distribution function also supported the above results, as can be seen in Figure S2.In general, AIMD simulation results showed that only the doping and substitution of P, S, and Cl atoms could construct relatively stable SFLPs site structures.Therefore, in the following discussion, P@In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 were selected as research objects to further screen potential and efficient NRR catalysts.

Adsorption and Protonation of Nitrogen
N 2 adsorption is the first and most essential step of the whole NRR process.As shown in Figure S3, the most stable structure of N 2 adsorbed on V-In 2 O 3 , P@In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 was optimized.The results show that among the four materials, N 2 could be stably adsorbed only on the P@In 2 O 3 surface.Simultaneously, the bond length data of N 2 on different materials were summarized in Table S1.The length of the adsorbed N-N bond on P@In 2 O 3 was 1.17 Å, while the length of the free N 2 bond was 1.108 Å, indicating that the N-N bond was activated.In addition, the adsorption energy of N 2 on different materials was calculated, as shown in Figure 3a.The adsorption energies of N 2 on V-In 2 O 3 , P@In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 were −0.15, −0.10, −0.46, and −0.11 eV, respectively.The more negative the adsorption energy, the stronger the adsorption of N 2 on the catalysts.The results further showed that P@In 2 O 3 can effectively adsorb N 2 molecules.It was also interesting to find that there was a good linear relationship between the bond length of In-N and the adsorption energy of N 2 , where the shorter the bond length, the stronger the adsorption.Variations in temperature and energy against the time for the AIMD simulations of (a) V-In2O3, (b) B@In2O3, (c) C@In2O3, (d) Cl@In2O3, (e) F@In2O3, (f) P@In2O3, (g) S@In2O3, and (h) Si@In2O3.The simulation was run under 400 K for 10ps with a time step of 1 fs.

Adsorption and Protonation of Nitrogen
N2 adsorption is the first and most essential step of the whole NRR process.As shown in Figure S3, the most stable structure of N2 adsorbed on V-In2O3, P@In2O3, S@In2O3, and Cl@In2O3 was optimized.The results show that among the four materials, N2 could be stably adsorbed only on the P@In2O3 surface.Simultaneously, the bond length data of N2 on different materials were summarized in Table S1.The length of the adsorbed N-N bond on calculated, as shown in Figure 3a.The adsorption energies of N2 on V-In2O3, P@In2O3, S@In2O3, and Cl@In2O3 were −0.15, −0.10, −0.46, and −0.11 eV, respectively.The more negative the adsorption energy, the stronger the adsorption of N2 on the catalysts.The results further showed that P@In2O3 can effectively adsorb N2 molecules.It was also interesting to find that there was a good linear relationship between the bond length of In-N and the adsorption energy of N2, where the shorter the bond length, the stronger the adsorption.Following the capture and adsorption of N2 molecules, the first hydrogenation step of N2* is critical and is often the potential determining step (PDS) in the protonation step.In addition, the Gibbs free energy (ΔG) of the first hydrogenation of N2* is often used as a descriptor for selecting catalysts with excellent catalytic performance.Therefore, we calculated the ΔG values of N2*→N2H* on V-In2O3, P@In2O3, and S@In2O3.Following the capture and adsorption of N2 molecules, the first hydrogenation step of N2* is critical and is often the potential determining step (PDS) in the protonation step.In addition, the Gibbs free energy (ΔG) of the first hydrogenation of N2* is often used as a descriptor for selecting catalysts with excellent catalytic performance.Therefore, we calculated the ΔG values of N2*→N2H* on V-In2O3, P@In2O3, S@In2O3, and Cl@In2O3 to screen potential NRR catalysts.As can be seen in Figure 3c, the results of the ΔG calculations show that the protonation of N2* on the surfaces of V-In2O3, S@In2O3, and Cl@In2O3 is extremely difficult.However, the free energy of the protonation of N2* on the surface of P@In2O3 is only 0.19 eV, which indicates that P@In2O3 can effectively activate the N-N triple bonds and lower the energy barrier of the protonation of N2*.Therefore, P@In2O3 is selected as a candidate NRR electrocatalyst for further investigation.Under the real conditions of electrocatalysis (in Following the capture and adsorption of N 2 molecules, the first hydrogenation step of N 2 * is critical and is often the potential determining step (PDS) in the protonation step.In addition, the Gibbs free energy (∆G) of the first hydrogenation of N 2 * is often used as a descriptor for selecting catalysts with excellent catalytic performance.Therefore, we calculated the ∆G values of N 2 *→N 2 H* on V-In 2 O 3 , P@In 2 O 3 , and S@In 2 O 3 .Following the capture and adsorption of N 2 molecules, the first hydrogenation step of N 2 * is critical and is often the potential determining step (PDS) in the protonation step.In addition, the Gibbs free energy (∆G) of the first hydrogenation of N 2 * is often used as a descriptor for selecting catalysts with excellent catalytic performance.Therefore, we calculated the ∆G values of N 2 *→N 2 H* on V-In 2 O 3 , P@In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 to screen potential NRR catalysts.As can be seen in Figure 3c, the results of the ∆G calculations show that the protonation of N 2 * on the surfaces of V-In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 is extremely difficult.However, the free energy of the protonation of N 2 * on the surface of P@In 2 O 3 is only 0.19 eV, which indicates that P@In 2 O 3 can effectively activate the N-N triple bonds and lower the energy barrier of the protonation of N 2 *.Therefore, P@In 2 O 3 is selected as a candidate NRR electrocatalyst for further investigation.Under the real conditions of electrocatalysis (in aqueous solution), HER is inevitable as the main side reaction of the NRR system [5][6][7].As show in Figure S4, the ∆G values of HER for V-In 2 O 3 , P@In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 are −0.91,−0.85, 0.15, and −0.08 eV, respectively.In addition, the potential difference between the N 2 * protonation and HER processes (U NNH* − U H* ) was calculated to evaluate the selectivity of electrocatalysis reactions with different catalysts, as shown in Figure 3d.The results show that P@In 2 O 3 has excellent catalytic selectivity for NRR and can effectively inhibit the HER, which is considered to be a promising electrocatalyst for NRR.

Catalytic Mechanism in Solution
To explore the reaction mechanism of NRR on P@In 2 O 3 , Gibbs free energies of different mechanisms on P@In 2 O 3 and V-In 2 O 3 were calculated and the results are shown in Figure 4.The implicit solvent model was used to consider the influence of the solvation effect on the reaction when calculating the Gibbs free energy diagram.The adsorption and desorption processes are not sensitive to the solvation effect [26], so the calculation started from N 2 * protonation.The possible NRR pathway on P@In 2 O 3 is shown in Figure 4a, and since N 2 * is adsorbed on the surface of P@In 2 O 3 in a side-on configuration, the enzymes and distal mechanisms were mainly considered for the NRR process.The results show that the PDS of the enzyme mechanism was NHNH* + H + /e − →NH 2 NH* and ∆G was 0.278 eV, where the PDS of the distal mechanism was NHN*+ H + /e − → NH 2 N* with ∆G of 0.285 eV.It was found that the PDS of the above two mechanism species was not the typical first hydrogenation process, and the protonation steps after PDS were spontaneous processes.At the same time, the limiting potential difference between the two mechanisms was only 0.07 eV, suggesting that NRR can occur on P@In 2 O 3 through both the enzyme and distal mechanisms.
selectivity of electrocatalysis reactions with different catalysts, as shown in Figure 3d.Th results show that P@In2O3 has excellent catalytic selectivity for NRR and can effectively inhibit the HER, which is considered to be a promising electrocatalyst for NRR.

Catalytic Mechanism in Solution
To explore the reaction mechanism of NRR on P@In2O3, Gibbs free energies of differ ent mechanisms on P@In2O3 and V-In2O3 were calculated and the results are shown in Figure 4.The implicit solvent model was used to consider the influence of the solvation effect on the reaction when calculating the Gibbs free energy diagram.The adsorption and desorption processes are not sensitive to the solvation effect [26], so the calculation started from N2* protonation.The possible NRR pathway on P@In2O3 is shown in Figure 4a, and since N2* is adsorbed on the surface of P@In2O3 in a side-on configuration, the enzyme and distal mechanisms were mainly considered for the NRR process.The results show that the PDS of the enzyme mechanism was NHNH* + H + /e − →NH2NH* and ΔG was 0.278 eV, where the PDS of the distal mechanism was NHN*+ H + /e − → NH2N* with ΔG of 0.285 eV.It was found that the PDS of the above two mechanism species was not the typical firs hydrogenation process, and the protonation steps after PDS were spontaneous processes At the same time, the limiting potential difference between the two mechanisms was only 0.07 eV, suggesting that NRR can occur on P@In2O3 through both the enzyme and dista mechanisms.As shown in Figure 4b, the N 2 molecule was adsorbed in an end-on configuration on the surface of V-In 2 O 3 , so the alternating, distal, and mixed mechanisms are considered.The results show that the PDS of these three reaction mechanisms was the first hydrogenation process (N 2 *+ H + /e − →NNH*) with the ∆G is 1.18 eV.In the case of the mixing mechanism, the protonation steps after the PDS were all processes of energy decrease and could be formed spontaneously.Meanwhile, for the other two mechanisms, the protonation steps after the formation of the first NH 3 were all spontaneous processes.It is noteworthy that the NNH* intermediate was not adsorbed at the constructed SFLPs site.Instead, it formed In-N bonds with the two In atoms (In 3 and In 4 ) at the oxygen vacancy and was adsorbed on the surface of V-In 2 O 3 by end-on configuration.This also shows that it is not sufficient to build SFLPs sites by oxygen vacancies alone.

Origin of Catalytic Activity
In order to identify the origin of catalytic activity at the P-doped SFPLs sites on P@In 2 O 3 , the ELF levels of V-In 2 O 3 , P@In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 were calculated and the results are shown in Figures 5 and S6.It can be clearly seen from the figure that there are obvious localization electrons at the oxygen vacancy (between In 3 and In 4 ) on the surfaces of V-In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 , while the electron density of the oxygen vacancy on P@In 2 O 3 equals the uniform density of the electron gas.These oxygen vacancies with localization electrons shield the catalytic activity of the constructed SFLPs sites by making it difficult to induce empty orbitals (Lewis acids) from the surrounding metal atoms.Interestingly, the results of partial charge density calculations for the valence band (VB) and conduction band (CB) support this opinion.The partial charge density values of the VB and CB of V-In 2 O 3 , P@In 2 O 3 , S@In 2 O 3 , and Cl@In 2 O 3 are shown in Figure S7.It was found that only the VB of P@In 2 O 3 was mainly contributed by the In atoms (In 3 and In 4 ) and the P atoms, whereas the VB of the other results was mainly contributed by the In 3 and In 4 atoms.As the major constituent in VB, the electron-rich P atom acts as a Lewis base site.It can interact with the empty π orbital of N 2 , providing electrons to activate N 2 molecules.In the meantime, the CB of P@In 2 O 3 is mainly contributed by the In 3 atoms.The empty orbital of the In 3 atom acts as a Lewis acid site and is able to accept the electron of the N 2 molecules.The apparent electron "donation-acceptance" process between the N 2 molecules and the SFLPs sites on P@In 2 O 3 promotes the effective activation of the N≡N bonds.
In order to verify the "donation-acceptance" process of the electrons, the difference charge density difference of the N 2 adsorbed by the different structures was calculated, as shown in Figure S8.The results show that electron transfer to N 2 occurred at the P@In 2 O 3 SFLPs site, which confirms the "donation-acceptance" mechanism of electrons.The Bader charge analysis also confirms this phenomenon and the calculated values in Table S3 indicate that the In 3 and P atoms are injected into the N 2 * molecule.In addition, the integrals of projected crystal orbital Hamilton population (ICOHP) data and bond length data suggest that the SFPLs sites can effectively activate the N≡N bonds, as shown in Tables S1 and S4.At the same time, the partial densities of states (PDOS) and the projected crystal orbital Hamilton population (pCOHP) of N 2 adsorbed on different structures were calculated to further investigate the "donation-acceptance" mechanism.As shown in Figure 5c, the 2p orbital of N 2 and the 3p orbital of P have a large hybrid region at the high-energy region (−2 eV), and a 2π* orbital of N 2 shift to the high-energy region near the Fermi level, indicating strong electron donation interaction and orbital interaction between the P atom and N 2 .As shown in Figure S9, the bonding orbitals are mainly contributed by the mutual hybridization of the s and p orbitals of the In 3 and N atoms between the −6 and −8 eV levels.This shows that the empty 5p orbital of the In atom can effectively accept the σ electron of the N 2 molecule and activate the N≡N bonds.

Materials and Methods
The DS-PAW program integrated into the Device Studio program [27] was used to perform the calculations under the framework of density functional theory (DFT).The interaction between valence electrons and the ionic core was described using the projection-augmented wave (PAW) base set [28,29], within the exchange-correlation function described by Perdew, Burke, and Ernzerhof (PBE) [30].The D3 correction method proposed by Grimme [31,32] was used to make up for the deficiency of the GGA functional description of van der Waals interaction.One-electron Kohn−Sham orbitals were expanded with a kinetic energy cutoff of 630 eV.The geometry optimization adopted the 3 × 2 × 1 Monkhorst-Pack grid of the K-points mesh to sample the Brillouin zone.The convergence criterion for self-consistent iteration was set to 1 × 10 −5 eV and the structures were fully relaxed until the final force on each atom was less than 0.03 eV Å −1 .To test the stability of the doped structures, Ab initio molecular dynamics (AIMD) simulations of all doping systems were performed using a canonical ensemble (NVT) with Nosé thermostat at 400K temperature with a total time of 10 ps and a time step of 1 fs.In addition, the density of states (DOS) was calculated using the Vienna Ab initio Simulation Package (VASP) [33,34].DOS calculations using 4 × 3 × 1 Monkhorst-Pack grids of k-points mesh.At the same time, the LORBSTER [35] program was used to analyze the properties of COHP [36].
In this work, the most thermodynamically stable body-centered cubic bixbyite (c-In2O3) crystal structure was selected, which had undergone extensive theoretical and experimental investigation [19,37,38].The surface of In2O3 (110) was modeled in a (1 × √2) supercell with five periodic atomic layers containing 40 In atoms and 60 O atoms.The bottom three layers of the model were fixed to simulate bulk properties, while the rest of

Materials and Methods
The DS-PAW program integrated into the Device Studio program [27] was used to perform the calculations under the framework of density functional theory (DFT).The interaction between valence electrons and the ionic core was described using the projectionaugmented wave (PAW) base set [28,29], within the exchange-correlation function described by Perdew, Burke, and Ernzerhof (PBE) [30].The D3 correction method proposed by Grimme [31,32] was used to make up for the deficiency of the GGA functional description of van der Waals interaction.One-electron Kohn−Sham orbitals were expanded with a kinetic energy cutoff of 630 eV.The geometry optimization adopted the 3 × 2 × 1 Monkhorst-Pack grid of the K-points mesh to sample the Brillouin zone.The convergence criterion for self-consistent iteration was set to 1 × 10 −5 eV and the structures were fully relaxed until the final force on each atom was less than 0.03 eV Å −1 .To test the stability of the doped structures, Ab initio molecular dynamics (AIMD) simulations of all doping systems were performed using a canonical ensemble (NVT) with Nosé thermostat at 400 K temperature with a total time of 10 ps and a time step of 1 fs.In addition, the density of states (DOS) was calculated using the Vienna Ab initio Simulation Package (VASP) [33,34].DOS calculations using 4 × 3 × 1 Monkhorst-Pack grids of k-points mesh.At the same time, the LORBSTER [35] program was used to analyze the properties of COHP [36].
In this work, the most thermodynamically stable body-centered cubic bixbyite (c-In 2 O 3 ) crystal structure was selected, which had undergone extensive theoretical and experimental investigation [19,37,38]

Figure 1 .
Figure 1.Optimized structure of (a) the perfect In2O3 (110) surface and (c) defective In2O3 (110) with one oxygen vacancy.Electron localization function of (b) the perfect In2O3 (110) surface and (d) defective In2O3 (110) with one oxygen vacancy.The black numbers in the diagram are the Bader charges of the atoms.

Figure 1 .
Figure 1.Optimized structure of (a) the perfect In 2 O 3 (110) surface and (c) defective In 2 O 3 (110) with one oxygen vacancy.Electron localization function of (b) the perfect In 2 O 3 (110) surface and (d) defective In 2 O 3 (110) with one oxygen vacancy.The black numbers in the diagram are the Bader charges of the atoms.As shown in Figure 1, in the absence of oxygen vacancies on the surface, In 3 and In 4 atoms (Lewis acids) can form localized structures similar to SFLPs with O atoms (Lewis bases) from the other chain, but the O 4 attached to In 3 and In 4 atoms shields them from their functionality as Lewis acids.The unsaturated ligand In atoms (In 4 and In 5 ) and O atoms (O 5 ) on the other chain form SFLPs sites when oxygen vacancies are present on the surface of In 2 O 3 (110) [22,24].Meanwhile, Bader charge results show that the atomic local charges of In 3 , In 4 , and O 5 are +1.193 and +1.267, −1.153, respectively.The results show that the creation of oxygen vacancies makes two In atoms (In 3 and In 4 ) on the In 2 O 3 (110) surface into potential electron donors, further suggesting that oxygen vacancy defects can construct SFLPs sites compared to a perfect In 2 O 3 (110) surface.However, the O 4 vacancy defects on the surface of In 2 O 3 (110) still hold a small charge, which has a shielding effect on

Figure 2 .
Figure 2. Variations in temperature and energy against the time for the AIMD simulations of (a) V-In 2 O 3 , (b) B@In 2 O 3 , (c) C@In 2 O 3 , (d) Cl@In 2 O 3 , (e) F@In 2 O 3 , (f) P@In 2 O 3 , (g) S@In 2 O 3 , and (h) Si@In 2 O 3 .The simulation was run under 400 K for 10ps with a time step of 1 fs.

Figure 3 .
Figure 3. (a) Adsorption energy of N2 on different materials.(b) The adsorption energy of N2 on different materials is a function of the In-N bond length.(c) The free energy of the first hydrogenation of N2 on different materials.(d) Diagram of the potential difference between the NRR and HER on different materials.

Figure 3 .
Figure 3. (a) Adsorption energy of N 2 on different materials.(b) The adsorption energy of N 2 on different materials is a function of the In-N bond length.(c) The free energy of the first hydrogenation of N 2 on different materials.(d) Diagram of the potential difference between the NRR and HER on different materials.

Figure 4 .
Figure 4. Gibbs free energy diagrams and intermediate structure of N 2 reduction process on (a) P@In 2 O 3 and (b) V-In 2 O 3 .

Figure 5 .
Figure 5. (a) Electron localization function and (b) partial charge density of the VB and CB of P@In2O3.The value of the isosurface is 0.002 e/Å −3 .(c) The pDOS and -pCOHP between N and P atoms on P@In2O3-adsorbed N2.

Figure 5 .
Figure 5. (a) Electron localization function and (b) partial charge density of the VB and CB of P@In 2 O 3 .The value of the isosurface is 0.002 e/Å −3 .(c) The pDOS and -pCOHP between N and P atoms on P@In 2 O 3 -adsorbed N 2 .
. The surface of In 2 O 3 (110) was modeled in a (1 × √ 2) supercell with five periodic atomic layers containing 40 In atoms and 60 O atoms.The bottom three layers of the model were fixed to simulate bulk properties, while the rest of the model was allowed to relax during the geometry optimization process.A vacuum gap