Investigation of dual atom doped single-layer MoS 2 for electrochemical reduction of carbon dioxide by first-principle calculations and machine-learning

first-principle


INTRODUCTION
The increasing consumption of fossil fuels has induced massive release of carbon dioxide (CO 2 ) in the atmosphere and caused severe energy crisis and environmental pollution on a global scale [1,2] .One sustainable approach is to decrease CO 2 emissions while transforming CO 2 into value-added products.Nevertheless, the CO 2 molecule is very stable, which requires a high activation energy to activate and break the C=O bond [3][4][5][6][7] .Among those developed methods [8,9] , the electrochemical CO 2 reduction reaction (CO 2 RR) is one promising solution and has received lots of experimental and theoretical attention owing to its simple operation, controllable selectivity, and practical potential for industrial applications [10,11] .
In particular, the single-atom catalysts (SACs) have been a rapidly developing field in recent years owing to their powerful potential for heterogeneous catalysis [12] .The well-defined active sites provide a great platform for investigating the reaction mechanism and establishing the correlation between structures and activity [3,11,[13][14][15][16][17][18][19][20] .Significant progress has been made in applying SACs for single-intermediate electrochemical reactions, i.e., the hydrogen evolution reaction (HER) [21][22][23][24] .The SACs also exhibit promising electrocatalytic applications in other types of multi-intermediate reactions, including oxygen reduction reactions (ORR) [25][26][27][28] , CO 2 RR [29][30][31] , and N 2 reduction reactions (NRR) [32,33] .The catalytic activity of SACs, however, is usually limited to the simple electronic structure and low density of metal active sites [34] .Meanwhile, the SACs tend to form metal clusters during experimental synthesis, which causes great challenges in the usage of SACs efficiently [18,35] .Moreover, due to the presence of only one type of active site, it is difficult to break the inherent linear relationship of adsorption strength of intermediates by SACs [36][37][38] .The catalytic activity can be regulated by balancing the adsorption of reaction intermediates on the catalyst surface [39,40] .
In this case, a promising strategy to regulate the adsorption of intermediates is via introducing a secondary metal site, as indicated by prior studies [41][42][43] .We have termed it as dual atom catalysts (DACs) [44] .On account of the synergetic effects of dual active sites, DACs can better maximize the catalytic potential of SACs for various multi-step reactions, leading to boosted catalytic performance [45][46][47][48][49] .For example, Yan et al. experimentally synthesized the Pt 2 dimer dispersed on graphene, which catalyzes the hydrolytic dehydrogenation of ammonia and boron at a reaction rate nearly 17-fold faster than that of a single Pt atom [50] .Ren et al. synthesized Fe-Ni DACs embedded in N-doped porous carbon as a highly efficient catalyst for CO 2 reduction [13] .In theory, Zhao et al. predicted that Cu 2 dimer loaded on porous C 2 N nanosheets has high selectivity for CH 4 production [51] .The Co-, Ni-, and Cu-based DACs are predicted to exhibit higher activity for O 2 reduction than the corresponding single-atom counterparts [45,52] .In order to obtain excellent transition metal (TM)-based DACs, an appropriate stabilizing substrate is essential, which not only offers the coordination sites for stably capturing metal atoms but also acts synergistically with the active center during the electrocatalytic process.Currently, a lot of experimental and theoretical studies focus on the N-doped carbon or graphene as the stabilizing substrate.Notably, during the synthesis of 2D nanosheets of molybdenum disulfide (MoS 2 ), inherent vacancy defects are very common and easy to form, mostly S vacancies [53][54][55] .Not only the single S vacancy but also the double S vacancies and clustered S vacancy line can be realized experimentally [56] .These S vacancies can be used as the anchor sites for catalytic atoms due to their high binding affinity for atoms and molecules.Experimentally, many isolated metal atoms, such as Co and Pt, have been successfully anchored at the single S vacancy of MoS 2 [57,58] .Thus, we hypothesized that the MoS 2 with available double S vacancies could also be a potential substrate to anchor DACs [59] .
In this work, we theoretically explored the CO 2 RR performance of a series of dual-metal (M: Cu, Fe, Ni, Mn, Cr, Co)-doped single-layered MoS 2 (denoted as MoS 2 -M 2 or MoS 2 -M 1 M 2 , Figure 1) via density functional theory (DFT) calculations.In the optimized homonuclear and heteronuclear DACs, some of the adjacent metal atoms will form the metal-metal bond [Figure 1B and D], while others will not [Figure 1A and C].The results showed that water molecules would first occupy the active site, which is difficult to desorb, and would stabilize MoS 2 -M 2 /M 1 M 2 for further CO 2 reduction.Among the examined 21 DAC compositions, MoS 2 -NiCr is identified as a highly promising candidate for catalyzing CO 2 reduction to CH 4 .More importantly, we incorporated random forest regression prediction in a machine learning approach by training the DFT calculated data to identify important feature factors that influence the activity of CO 2 RR and the adsorption of the key *CO intermediate, where the distance between the two metal centers and the number of electrons in the outer layers of the metal atoms play a significant role.

COMPUTATIONAL DETAILS
All the spin-unrestricted DFT calculations are performed in the DMol 3 code [60,61] .The exchange-correlation effect is described via the Perdew-Burke-Ernzerhof (PBE) [62] functional of the generalized gradient approximation (GGA) [63] .The double numerical plus polarization (DNP) basis is employed using the DFT semi-core pseudopotential (DSPP) for the core treatment.The van der Waals interaction between CO 2 RR intermediates and DACs is described by empirical dispersion-corrected DFT (DFT-D3).To simulate the aqueous solvent environment, a conductor-like screening model (COSMO) is adopted [64][65][66] .In geometric optimization, the convergence threshold of energy is 2 × 10 -5 Ha; the maximum displacement is set as 0.005 Å, and the force applied to each atom is 0.004 Ha/Å.A 4 × 4 × 1 rectangular 2H-MoS 2 supercell with adjacent double S vacancies was constructed, and a 3 × 3 × 1 Monkhorst-Pack k-mesh was used to sample the Brillouin zone.Moreover, an 18 Å vacuum space was set to avoid interactions of adjacent images.The canonical ensemble (NVT) ab initio molecular dynamics (AIMD) simulations are performed in a Nosé-Hoover thermostat at 300K for five picoseconds (ps) in a time step of one femtosecond (fs).
The formation energies of homonuclear and heteronuclear DACs, E f , are calculated by the following equation [67] : where N represents the number of doped atoms, E total is the total energy of DACs, E MoS2 denotes the energy of MoS 2 substrate with double S vacancies, and E coh is the cohesive energy of the dopant.
According to the computational hydrogen electrode (CHE) model [68] , the Gibbs free energy change (ΔG) of each elementary reaction step of CO 2 RR is calculated by ΔG = ΔE + ΔZPE -TΔS, where ΔE is the reaction energy change calculated by DFT calculations, while ΔZPE and TΔS represent the difference in zero-point and entropy change at 298 K.For gas phase molecules, the entropy is derived from the NIST database and details are provided in the Supporting Information [Supplementary Table 1].
The limiting potential (U L ) of the reaction is calculated as U L = -ΔG max /e, where ΔG max corresponds to the maximum free energy change among all the CO 2 RR elementary steps.

Stability
According to the above equations, we calculated the formation energy (E f ) [Figure 2A] to assess the thermodynamic stability of the six kinds of homonuclear DACs and 15 kinds of heteronuclear DACs.The E f of all DACs were negative, ranging from -4.92 to -6.16 eV [Supplementary Table 2], indicating high thermodynamic stability.In addition, the AIMD simulations are carried out to verify the dynamic stability of the DACs.From Figure 2B to E (MoS 2 -MnCr, MoS 2 -FeMn, MoS 2 -CrCo, and MoS 2 -NiCr), the temperature slightly fluctuates around 300 K, and the energy changes within ± 0.01 eV.No obvious deformation occurs in these frameworks during the AIMD simulation, further confirming the high dynamic stability of the catalysts.

Activation of CO 2
The activation of CO 2 over the active center is the first step during electrocatalytic CO 2 RR.However, from Supplementary Table 3, water adsorption is energetically more preferable than CO 2 except for MoS 2 -MnCr.From the optimized structures [Supplementary Figure 1], the O atom of the adsorbed H 2 O is coordinated to one single metal center or the metal-metal bridge site.Note that the adsorbed water molecule cannot split spontaneously due to its highly endothermic dissociation process (0.28~1.07 eV, Supplementary Table 4).Therefore, we subsequently considered CO 2 adsorption and reduction after water molecules first occupy the active site.The binding interaction between CO 2 and MoS 2 -embedded DACs ranges from -0.45 to -0.98 eV [Supplementary Figure 2].Especially, the binding strength of CO 2 on MoS 2 -NiCr is strong with high adsorption free energy of -0.98 eV and curved O−C−O bond angle of 139.106°, which is accompanied by considerable charge transfer of around 0.6 |e| from catalyst to CO 2 [Figure 3A].Note that in the case of MoS 2 -NiCr, the two O atoms of CO 2 are coordinated to the Ni and Cr centers, respectively.In other DAC systems, only one of the O atoms of CO 2 is coordinated to one metal center, which is accompanied by weaker adsorption strength (-0.45~-0.88eV) and less charge transfer (0.01~0.02 |e|) between the anchored metal dimer and CO 2 (MoS 2 -CrCo is shown as an example in Figure 3B).However, the favorable adsorption proves that CO 2 molecules can be effectively captured and activated by the metal dimers embedded in sulfur vacancies of MoS 2 .

Scaling relations
In most cases, the potential limiting step for CO 2 electroreduction is the reduction of *COOH to *CO (twoelectron reduction) or the reduction of *CO to *CHO (deep reduction).Thus, the overall catalytic efficiency depends strongly on the adsorption energies of *COOH [E ads (COOH)], *CO [E ads (CO)], and *CHO [E ads (CHO)] [69,70] .The reduction of CO 2 to CO involves a two-step electroreduction, i.e.,  [36,71] .Therefore, we first examined the adsorption of *CHO, *COOH and *CO.The detailed adsorption energy and the adsorption geometries are given in Supplementary Table 5 and Supplementary Figure 3. From Figure 4, the scaling relations are completely broken compared to those observed in pure metal surfaces.In the linear diagram of E ads (COOH) vs. E ads (CO) [Figure 4A], the scattered points are distributed in the whole region, indicating that the DAC electrocatalysts can effectively break the linear relationship.Note that most of these investigated DACs have strong *CO adsorption, which means that the generated CO would undergo further hydrogenation to form deep reduction products.For the relationship between E ads (CO) and E ads (CHO) [Figure 4B], the scaling relationships are slightly weakened with scattered points compared with those of the pure metal surfaces.In addition, NiCr and CrCo are two special cases that deviate greatly from the linear relationship of pure metal surfaces and show small differences between E ads (CO) and E ads (CHO); thus, they can be used as candidates to achieve the desired low overpotential for deep reduction products.Consequently, we selected these two systems for our subsequent calculations.We also analyzed the projected density of states (PDOS) of the NiCr and CrCo candidates after CO adsorption.From Supplementary Figure 4, within the energy region from -8 to -4 eV below the Fermi level, it can be clearly seen that the p orbital of C in the adsorbed *CO is strongly hybridized with the d orbital of the metal, proving that *CO has strong adsorption with the metal active site, which is beneficial to the deep reduction reaction of CO [72] .

The pathway of CO 2 RR
In the following, we systematically investigated the reduction pathway of CO 2 RR on MoS 2 -NiCr and MoS 2 -CrCo after the formation of strongly bound *CO (*CO is firstly produced through the two-electron DACs and the transition metal surfaces.The linear proportional relationships between the adsorbents were obtained on Ni, Cu, Ag, Pd, Au, Pt, and Rh surfaces [70] .DACs: Dual atom catalysts.pathway: CO 2 → *COOH → CO).The *CO can be further reduced to other C1 products, such as HCHO, CH 3 OH, and CH 4 .The free energy diagrams of all the possible C1 products are shown in Figure 5, the structural schematics are shown in Supplementary Figure 5, and the detailed data of free energy are provided in Supplementary Table 6

Selectivity of CO 2 RR vs. HER
In CO 2 RR, the HER always competes with CO 2 reduction in an aqueous solution [73] .Firstly, the occupation of sites was initially considered, and as shown in Supplementary Table 7, H adsorption in most dual-atom systems is not as strong as H 2 O and CO 2 adsorption.Therefore, the diatomic sites are more likely to take the CO 2 RR path.Secondly, it is necessary to assess the selectivity of CO 2 RR by comparing its limiting potential (U L ).In the CO 2 reduction process, we consider the comparison between the limiting potential of the electrochemical steps and that of HER.Accordingly, a more positive value of ΔU L [U L (CO 2 RR) -U L (HER)] implies higher reaction selectivity for CO 2 reduction.From Figure 7, the ΔU L of the NiCr dimer (0.26 V) is located in the upper right corner, indicating its high CO 2 RR selectivity, while the ΔU L of the CrCo dimer is close to 0, indicating its poor selectivity.Furthermore, the ideal electrocatalysts should be well accompanied by effective CO 2 activation.In other words, the strong adsorption of CO 2 over the catalyst can inhibit H on the catalyst surface, thus hindering the competitive HER as the CO 2 will occupy the available active sites [6,74,75] .The calculated adsorption free energies of *H on MoS 2 -NiCr and MoS 2 -CrCo metal sites are -0.84 and -0.45 eV, respectively (inset in Figure 7), while the adsorption free energies of CO 2 are -0.98 and -0.68 eV, respectively.This indicates that CO 2 adsorption is more favorable than H* adsorption.Hence, the adsorption of CO 2 is preferred, while the adsorption of *H is hindered.By comparing the reaction activity and selectivity, the MoS 2 -NiCr is screened to be a promising dual-site electrocatalyst to promote the CH 4 formation with moderate rate-determining step (RDS) barriers and high CO 2 RR selectivity over HER.Our theoretical prediction will provide useful insights for future experimental verification of the high electrocatalytic performance of Ni-Cr DACs supported on MoS 2 substrates.

Machine learning analysis
From the data calculated above, the CO 2 RR activity of DACs and the binding strength of the intermediate *CO are closely related, with weak CO binding favoring CO gas desorption and strong CO binding facilitating adsorption of added H to *CHO.At present, the underlying factors affecting CO 2 RR activity remain to be discovered.Furthermore, DAC systems are much more complex than TM surfaces.Therefore, it is difficult to accurately describe the CO 2 RR activity of DACs with only one descriptor.Without performing intensive DFT calculations, there is a strong need to identify more readily available variables to describe the CO 2 RR activity of DACs.Thus, by using a machine learning approach, we investigated the correlation between ΔG *CO or U L (CO) and the intrinsic factors of six homonuclear and 15 heteronuclear catalysts.Proper feature selection is essential for machine learning models to identify the hidden rules behind the input data.In our work, we considered seven very basic parameters to describe the geometric and electronic properties of DACs, including the distance between two metal atoms (d M-M ), the average distance between two metal atoms and six Mo atoms (d M-Mo ), the radii of two metal atoms (R 1 and R 2 ), the number of outer electrons of two metal atoms (Ne 1 and Ne 2 ), the Pauling electronegativity (P 1 and P 2 ), the first ionization energy (I 1 and I 2 ), and the electron affinity (A 1 and A 2 ).Importantly, we examined the correlations between the factors, and as can be seen in Supplementary Figure 6, most combinations of factors are poorly related to each other.Some of the factors vary with the regularity of the periodic table, e.g., Ne, R, etc.Thus, these factors and coefficients are variables that can be approximated as independent of each other.It is important to note that we augmented the data for all the DACs studied because MoS 2 -M 1 M 2 and MoS 2 -M 2 M 1 correspond to two different sets of variable combinations [Supplementary Table 8].In this way, each DAC possesses two sets of input features.
We used a random forest regression algorithm from the scikit-learn toolkit [76] .The DFT computed ΔG *CO values were then compared with the values predicted based on the random forest study.The DFTcomputed ΔG *CO input data were randomly perturbed and divided into a training set and a test set in a ratio of 6:1.As shown in Figure 8A, the predicted values based on the random forest have a similar trend to the values calculated by DFT, with a lower mean square error of 0.058.There is a high R 2 value, 0.93 for the training score and 0.91 for the test score, indicating that the random forest prediction algorithm adequately trained the model by learning the factors inherent in the model to reach an accurate prediction.The importance of the seven features on ΔG *CO was also evaluated.In Figure 8B, the distance between the metals (d M-M ) is the most influential on ΔG *CO , with a feature importance value of 0.34, while the sum of feature importance values of the radius of the metal atoms (R 1 and R 2 ) and the distance between the metal and the Mo atoms (d M-Mo ) is only 0.01.That means that the synergistic effect between the DACs has a strong influence on the catalytic efficiency.In addition, the outer electron number (Ne) of the metal atom also plays an important role, with a sum of feature importance (Ne 1 + Ne 2 ) of 0.20, which can be interpreted as forming a metal-metal bond between DACs that cannot efficiently bind the CO 2 RR intermediates.However, the importance of the remaining three features (P, I, and A) was relatively insignificant.We also predicted the limiting potential in the CO 2 → CO process based on the Random Forest algorithm, and the predictions were highly similar to the DFT [Figure 8C and D].The feature importance pie charts show similarities to those described above.Machine learning links the correlations between the intrinsic structure and the properties, providing powerful insights into the understanding of the CO 2 RR activity of DACs.Particularly, since the activity of the dual-atom catalyst in the CO 2 RR process is closely correlated with these important factors, we can apply these descriptors to predict the activity of other dual-atom compositions.

Potential limitations
There is one thing that needs to be added: our work is based on first-principles calculations to investigate the activity of electrochemical reduction of CO 2 by dual atom doped single-layer MoS 2 .From a theoretical point of view, the DACs predicted by us have relatively negative formation energies (E f ) and stable structures through AIMD, which indicates that it is feasible to synthesize such structures.Recently, an ingenious approach has successfully assembled DACs of Ni and Fe into the interlayer of MoS 2 [77] .These DACs exhibit higher catalytic activity toward acidic water splitting.Our predicted MoS 2 -FeNi structure was confirmed through this experiment.Therefore, these structures that we predict, namely the doping of different dual atoms (Cu, Co, Cr, Mn, etc.) in the single-layer MoS 2 , are expected to be realized in the future.
Furthermore, our computations rely on a traditional CHE model that neglects the display potential and display solvation factors, which do affect the precision of the performance evaluation to some extent.
Although the method has some limitations in the evaluation of activity due to the significant computational cost savings and relatively accurate simulation accuracy of the CHE model, this method is very popular for large-scale prediction and performance screening of new materials [78][79][80][81][82][83] .In other words, while taking into account the calculation speed and accuracy, the performance evaluation at the same atomic level is also of great reference significance.

CONCLUSION
In summary, the reaction activity of various dual atoms embedded in defective MoS 2 monolayers, named MoS 2 -M 2 /M1M2, for CO 2 reduction was systematically studied using computational DFT approaches.We theoretically studied 21 dimer electrocatalysts.Our results demonstrate that the defective MoS 2 monolayer with double S vacancies can anchor the two TM atoms stably.Through the analysis of the adsorption relationship of key intermediates, it was found that MoS 2 -CrCo and MoS 2 -NiCr candidates significantly deviated from the linear relationship; thus, they were selected for further analysis of deep reduction.We found that MoS 2 -CrCo shows a lower barrier energy for CH 4 production (0.44 eV), but its selectivity (ΔU L = 0.02 eV) over competitive HER is low.In contrast, the MoS 2 -NiCr system possesses superior selectivity (ΔU L = 0.26 eV) and catalytic activity for CH 4 production with a low rate-determining electrochemical barrier of 0.58 V.In the whole reaction process, water exists in the form of coordination in the formation process of C1 products.Moreover, our machine learning study demonstrated that adsorption of the key *CO intermediate and CO 2 RR activity can be well described by some basic parameters, such as the distance between the center of metal atoms and the number of outer electrons of the metal atoms.This work presents an atomic-level investigation of the screening and design of novel DACs supported on defective MoS 2 , providing useful insights for further investigations, including theoretical and experimental attempts.

Figure 1 .
Figure 1.The geometric structure of MoS 2 -M 2 and MoS 2 -M 1 M 2 .Some of the adjacent metal atoms will form metal-metal bonds (B and D), while others will not (A and C).The dark cyan and yellow balls represent Mo and S atoms, respectively, and the dark blue and purple balls represent the two TM atoms.TM: Transition metal.

Figure 3 .
Figure 3.The adsorption geometry and charge density difference of CO 2 adsorption over (A) MoS 2 -NiCr and (B) MoS 2 -CrCo DACs, with an isosurface level of 0.002 e•Å -3 .The green and red regions represent electron depletion and accumulation, respectively.DACs: Dual atom catalysts.

Figure 4 .
Figure 4. Relationship between the binding energies (A) E ads (COOH) vs. E ads (CO) and (B) E ads (CHO) vs. E ads (CO) of MoS 2 embeddedDACs and the transition metal surfaces.The linear proportional relationships between the adsorbents were obtained on Ni, Cu, Ag, Pd, Au, Pt, and Rh surfaces[70] .DACs: Dual atom catalysts.

Figure 7 .
Figure 7.The limiting potential difference between CO 2 reduction [U L (CO 2 RR)] and HER [U L (HER)] over MoS 2 -NiCr and MoS 2 -CrCo catalysts.The inset in the figure is the free energy diagram of HER.CO 2 RR: Carbon dioxide reduction reaction; HER: hydrogen evolution reaction; U L : limiting potential.

Figure 8 .
Figure 8.Comparison of (A) ΔG *CO and (C) [U L (CO)] obtained by DFT with values predicted by machine learning; (B and D) feature importance based on a random forest regression.DFT: Density functional theory.
. One can see that the hydrogenation of *CO to *CHO is energetically more favorable than the formation of *COH.Moreover, for both MoS 2 -NiCr and MoS 2 -CrCo, the generation of CH 4 needs lower energy input than the generation of HCHO and CH 3 OH, indicating that CH 4 would be the main reduction product of CO 2 RR.From Figure5A, the potential limiting step of CH 4 formation on MoS 2 -NiCr corresponds to *CO reduction to *CHO and *OH reduction to H 2 O, which need comparable endothermic free energy of 0.56 and 0.58 eV, respectively.While on MoS 2 -CrCo [Figure5B], the potential limiting step of CH 4 formation corresponds to *CO 2 reduction to *COOH or *CO reduction to *CHO, which needs comparable endothermic free energy of 0.44 and 0.43 eV, respectively.Figure6shows the detailed geometry of reaction intermediates during CH 4 formation on MoS 2 -NiCr and MoS 2 -CrCo catalysts.On MoS 2 -NiCr [Figure 6A], the various intermediates (*CHO, *CH 2 O, *CH 3 O, *O, and *OH) from a deep reduction of *CO are all coordinated to both the Ni and Cr atoms.In the case of MoS 2 -CrCo [Figure 6B], the reaction intermediates are mainly singly coordinated to the Cr atom.This indicates that the Ni and Cr centers in MoS 2 -NiCr work synergistically as dual active sites to affect the adsorption and bonding of CO 2 RR intermediates, while in MoS 2 -CrCo, only the Cr center plays the key role and functions as the single active site for CO 2 RR.