High-Performance of Electrocatalytic CO 2 Reduction on Defective Graphene-Supported Cu 4 S 2 Cluster

: Electrochemical CO 2 reduction reaction (CO 2 RR) to high-value chemicals is one of the most splendid approaches to mitigating environmental threats and energy shortage. In this study, the catalytic performance of CO 2 RR on defective graphene-supported Cu 4 S 2 clusters as well as isolated Cu 4 X n (X = O, S, Se; n = 2, 4) was systematically investigated based on density functional theory (DFT) computations. Calculation results revealed that the most thermodynamically feasible product is CH 3 OH among the C1 products on Cu 4 X 2 clusters, in which the Cu 4 S 2 cluster has the best activity concerning CH 3 OH synthesis with a limiting potential of − 0.48 V. When the Cu 4 S 2 cluster was further supported on defective graphene, the strong interaction between cluster and substrate could greatly improve the performance via tuning the electronic structure and improving the stability of the Cu 4 S 2 cluster. The calculated free energy diagram indicated that it is also more energetically preferable for CH 3 OH production with a low limiting potential of − 0.35 V. Besides, the defective graphene support has a signiﬁcant ability to suppress the competing reactions, such as the hydrogen evolution reaction (HER) and CO and HCOOH production. Geometric structures, limiting potentials, and reduction pathways were also discussed to gain insight into the reaction mechanism and to ﬁnd the minimum-energy pathway for C1 products. We hope this work will provide theoretical reference for designing and developing advanced supported Cu-based electrocatalysts for CO 2 reduction.


Introduction
Global warming and energy shortage have become two major issues in the 21st century. The excessive CO 2 emission caused by rapid consumption of fossil resources calls for the instant requirement of clean and renewable energy [1]. CO 2 itself is an inert molecule ∆G 0 f = −394.4 KJ/mol [2], so it is urgent to explore a robust catalyst for its reduction and conversion. In recent years, electrochemical CO 2 reduction reaction (CO 2 RR) has attracted widespread attention as one of the most excellent methods to alleviate energy shortage and reduce the global carbon footprint [3,4]. In CO 2 RR, various metal electrodes play a dual role, not only as an electron conducting medium but also as catalytically active sites. Among them, Cu-based material has excellent activity and selectivity to reduce CO 2 into hydrocarbons such as methanol, methane, formic acid, ethylene, etc. [5,6]. In these products, methanol is a desired hydrocarbon because it can be used as feedstock for direct methanol fuel cells [7]. More importantly, the metallic Cu has been verified as an active catalyst for methanol synthesis with a turnover frequency even comparable with that of industrial catalysts [8].
With the development from single-crystal metal electrodes to supported metal nanoparticle electrodes [9][10][11], the supported sub-nanometric Cu clusters [12,13] have received considerable attention in the past few years. The superior catalytic performance comes from three factors: The unique electronic structure of the nano-size cluster; the high specific surface area; and the significant effect of the substrate support. In the catalytic experiments, the support has an important role in immobilizing clusters and preventing them from aggregation. For instance, defect engineering of graphene [14,15] has been used to immobilize the cluster against sinter by providing strong anchoring sites which could effectively regulate the electronic properties and catalytic performance in CO 2 reduction. In particular, theoretical results have confirmed that the defective graphene supported Cu (Cu 55 [16] and Cu 4 [17]) nanocluster can improve the CO 2 RR selectivity for hydrocarbon fuels at low overpotential. Additionally, Ma et al. reported that the experimental realization of the Cu@Cu x O nanoparticles decorated on defective graphene achieved a high CO 2 reduction Faradaic efficiency of 60.6% [18]. By analyzing experiments and DFT calculations, Hu et al. demonstrated that the defective graphene-supported Cu 13 cluster exhibits excellent performance for CO 2 RR with a maximum Faradaic efficiency towards methane of 81.7% and an outstanding stability of 40 h [19].
It is noteworthy that the activity and selectivity of CO 2 RR are highly sensitive to the structure of the supported Cu-based cluster [20]. To date, the size-selected Cu 4 cluster is found to present good catalytic activity for producing hydrocarbons at near atmospheric pressure by comparing with other larger-size catalysts [21,22]. In the experiment, the subnanometer-sized Cu nanoclusters have been successfully prepared. For example, the size-controlled naked Cu cluster [23][24][25] can be synthesized by the electrochemical oxidation-reduction method and the one-pot procedure based on wet chemical reduction, or combination of gas-phase cluster ion sources, mass spectrometry, and soft-landing techniques. However, the impurity-doping effects could further enhance the durability, selectivity, and chemical activity of clusters. The modified Cu nanoclusters [26,27] [28]. Similarly, Zhao et al. reported that the CuS nanosheet arrays supported on nickel foam are a robust catalyst for CO 2 RR toward methane, which has high Faradaic efficiency of 73 ± 5% and a low Tafel slope of 57 mV dec −1 [29]. Lu et al. have synthesized Cu 3 -X clusters (Cu 3 -CI, Cu 3 -Br, Cu 3 -NO 3 ) which have high selectivity for the electrocatalytic reduction of CO 2 to C 2 H 4 with a maximum Faradaic efficiency of 55.01% [30].
It can be seen that control type, location, and content of the dopant atoms are important in terms of affecting the catalytic performance for the modified cluster [31,32]. The electronic structure of the nanocluster can be easily regulated by per-oxidization or heteroatom modulation, providing remarkable possibilities to tune the catalytic performance in CO 2 RR. What is the effect of doped clusters on the catalytic activity and selectivity of CO 2 reduction with different valence states Cu as electrocatalyst? How does the defective graphene support affect their performance? To answer these questions, high-precision computational methods are first used to give an in-depth understanding of the CO 2 reduction mechanism of the isolated Cu 4 X n clusters, to identify the promising candidate with optimal catalytic performance. On this basis, the CO 2 RR mechanism of the Cu 4 S 2 cluster adsorbed on defective graphene is considered as a representative model to boost its stability and application. The catalytic activity and product selectivity are evaluated via analyzing the limiting potential of CO 2 RR towards different C1 products. The geometric structure, charge distribution, and all possible reaction pathways were also discussed to help further understand the reaction mechanism and catalytic properties. We also hope that this work can give theoretical guidance for experimental research and provide reference for the design of high-performance supported Cu-based catalysts in the future. Potential energy surfaces (PESs) of the Cu 4 X 2 clusters were probed via global optimization algorithms to locate the global minima (GM) structure and low-lying energy isomers (details in Figure S1). The GM of the Cu 4 O 2 , Cu 4 S 2 , Cu 4 Se 2 cluster exhibit a similar configuration (in Table 1), in which four Cu atoms and one doped atom form a double trigonal pyramidal structure, and the other doped atom binds with two Cu atoms on the side edge of Cu 4 . With the decrease of electronegativity of the doped atom, the Cu 4 X 2 configuration exhibits small stretch in the xy plane; as for the vertical direction, the lower half (triangular pyramid of Cu 4 structure) presents a slight compression, while the upper half presents a small stretch. It is noteworthy that the structure of the Cu 4 S 2 cluster predicted in this paper is partly consistent with the lowest energy structure of the Cu n S (n = 4) cluster proposed by Li et al. [33], verifying the solidity of our results. The natural bond orbital (NBO) analysis presented in Table S2 shows that the higher electronegativity dopant atoms can form the positively charged Cu sites, which tends to attract the electronegative oxygen in CO 2 , thus facilitating the initial step of CO 2 RR.

Results and Discussion
It is easy to understand that the ability of capturing CO 2 molecules is essential as the introductory step in CO 2 RR. The electrochemical adsorption of CO 2 molecules can be typically represented as The electrochemical adsorption of CO 2 involves one (H + + e -) pair transfer to the CO 2 species, where CO 2 is first hydrogenated to generate hydroxyl or carboxyl groups; it is usually controlled via an applied potential. After testing the possible electrochemical adsorption sites of CO 2 , the calculated reaction free energy of COOH * formation is 0.41 eV, 1.02 eV, and 0.73 eV lower than that of OCHO * formation on the Cu 4 O 2 , Cu 4 S 2 , Cu 4 Se 2 clusters, respectively. This means that the intermediate COOH * is more preferable in this case, and the lowest-energy configuration is located and presented in Table 1. As seen in Table 1, all reaction free energies are negative, indicating that the electrochemical absorption of CO 2 can occur spontaneously. In these three clusters, E react of COOH * on Cu 4 S 2 and Cu 4 Se 2 are close, but still larger than that of Cu 4 O 2 , suggesting the former two clusters are more active for CO 2 activation. In the adsorbate COOH, both C and O are binding with the Cu atoms. The binding distances d (Cu−C) and d (Cu−O) enlarge slightly with the decrease of electronegativity of the chalcogen; that is in line with the trend of charge transfer between cluster and adsorbed species. The charge redistributions before and after electrochemical CO 2 adsorption are listed in Table S3. The NBO analysis revealed that the adsorption process is accompanied with partial charge transfer of about −0.5 from the Cu 4 X 2 cluster (electron donor) to the COOH adsorbate (electron acceptor). Meanwhile, the transferred charge decreases from Cu 4 O 2 to Cu 4 S 2 (or Cu 4 Se 2 ), resulting in the weakness of Cu-O and Cu-C bonding between COOH * and Cu 4 S 2 (or Cu 4 Se 2 ). In addition, the structures of COOH * on three different clusters are similar, in which ∠ O−C−O bends significantly compared with the linear CO 2 molecule. Besides, both C-O bonds are elongated to ensure that COOH * retains sufficient activity to continue the subsequent hydrogenation reactions. Notice that the adsorption of COOH * does not change the main structure of doped clusters, which is important to maintain the stability of the Cu 4 X 2 cluster in the CO 2 RR process. Table 1. The GM of the Cu 4 X 2 cluster and structural details of electrochemical CO 2 adsorption on the Cu 4 X 2 cluster.

Cluster Electrochemical Adsorption of CO
The red, yellow, blue, orange, gray, and white spheres denote O, S, Se, Cu, C, and H atoms, respectively. The unit for E react is eV, d isÅ, and ∠ O−C−O is • .

Reaction Pathway of CO 2 RR on the Cu 4 X 2 Cluster
In order to obtain a better understanding of the reaction mechanism of electrocatalytic CO 2 reduction, the reaction networks for each elementary step were analyzed to identify the most thermodynamically favorable pathway and major side reactions of each doped cluster. The reaction schemes of all possible pathways for the reduction of CO 2 to methanol or methane, as well as their competing reactions, are given in the Supplementary Materials. The lowest energy pathway of CO 2 RR to CH 3 OH and CH 4 is shown in Figure 1, which helps us to explore the catalytic performance explicitly.
As shown in Figure 1, the free energy diagram of CO 2 RR pathways on the Cu 4 X 2 clusters starts from the intermediate COOH * . The lowest energy pathway for CH 3 OH formation is identified (pathway 1, in blue), and the other pathway (pathway 2, in red) represents CO 2 RR towards CH 4 formation. The electrochemical adsorption of CO 2 on the Cu 4 O 2 , Cu 4 S 2 , and Cu 4 Se 2 clusters are accompanied with a spontaneous exothermic process with ∆G of −0.24 eV, −0.65 eV, and −0.57 eV, respectively, which means the initial step can effectively adsorb COOH * (state 1). Subsequently, the hydroxyl group is hydrogenated on C to from HCOOH * accompanied by a Cu-C bond cleavage or a Cu-O bond cleavage caused by hydrogenation on O; these are both endothermic reactions, while the ∆G of the former is 0.58 eV, 0.50 eV, and 0.12 eV lower than the latter on the Cu 4 O 2 , Cu 4 S 2 , and Cu 4 Se 2 clusters, respectively. Then, the intermediate HCOOH * (state 2) undergoes three degrees of hydrogenation to generate CHO * (state 3), CH 2 O * (state 4), and CH 3 O * (state 5). Finally, the adsorbate CH 3 O is hydrogenated at the oxygen atom to produce CH 3 OH. The rate-determining steps of pathway 1 on the Cu 4 O 2 and Cu 4 S 2 cluster are both the hydrogenation of COOH * to form HCOOH * (from 1→ 2), which is endergonic, and require potentials of −0.56 V and −0.48 V, respectively. However, the rate-determining step for the Cu 4 Se 2 cluster is created concerning CHO * (from 2→ 3) via a (H + + e -) pair transfer from a solution with a larger potential of −0.82 V. It is obvious that the Cu 4 S 2 cluster has the minimum rate limiting reaction energy in methanol formation among these clusters, indicating that it can be considered as a promising electrocatalyst for CH 3 OH formation during the CO 2 reduction process. In pathway 2 for methane formation, the COOH * undergoes a series of hydrogenations and takes the form of At this point, desorption of CH 4 is the most difficult step in pathway 2 with limiting potentials of −0.56 V, −0.79 V, and −1.27 V for the Cu 4 O 2 , Cu 4 S 2 , and Cu 4 Se 2 clusters, respectively. A free energy profile shows that there are partially overlapping intermediates between pathways 1 and 2, while the sixth step is the product-determining step for producing methanol or methane. The kinetic barrier diagram of the sixth step on the Cu 4 O 2 cluster, the Cu 4 S 2 cluster, and the Cu 4 Se 2 cluster is shown in Figure S5. The corresponding results show that the energy barrier for CH 3 OH formation is 1.73 eV, 0.99 eV, and 1.49 eV lower than that of CH 4 formation, respectively. Combined with the free energy diagram, the methanol production step is a spontaneous exothermic process, while the release of methane is endothermic, suggesting the methanol product is both thermodynamically and kinetically favored on the Cu 4 X 2 clusters.

The Catalytic Performance on Cu 4 X 2 Cluster
The activity of the CO 2 -to-CH 3 OH reaction is governed by the limiting potentials, which are −0.56 V, −0.48 V, and −0.82 V for the Cu 4 O 2 , Cu 4 S 2 , and Cu 4 Se 2 clusters, respectively. The results indicated that both the Cu 4 O 2 and Cu 4 S 2 clusters exhibit excellent electrocatalytic activity for methanol synthesis. The presence of two heteroatoms can tune the electronic structure of the adjacent Cu atom, leading to partial charge transfer from Cu to chalcogen. Therefore, the oxidized Cu + site can effectively promote CO 2 activation, which is an important intermediate for triggering an electroreduction reaction. For the selectivity of the Cu 4 X 2 cluster, the CH 3 OH formation step is a spontaneous exothermic process, while the branch of CH 4 production requires an endothermic process to proceed, in which the intermediate CH 3 O * is considered as a bifurcation for these two pathways. Therefore, the Cu 4 X 2 cluster exhibits high CO 2 selectivity for CH 3 OH synthesis, instead of producing CH 4 via CHO * → CH 2 O * → CH 3 O * → CH 4 on the Cu 4 X cluster [34]. In order to analyze the variation of product selectivity with dopant content, we then examine the bond length, bond energy, and adsorption energy of the adsorbate CH 3 O on Cu 4 X 2 and Cu 4 X clusters. As shown in Table 2, the CH 3 O * binds to the cluster via an oxygen atom; product selectivity at this point depends on how easily the Cu-O or C-O bond can be cleaved. With the dopant content increasing, the bond energy of Cu-O becomes smaller and that of C-O becomes larger, which is consistent with the changes in bond length. The adsorption energy of CH 3 O * on the Cu 4 X 2 cluster is smaller than that on the Cu 4 X cluster, which implies the subsequent reaction prefers the desorption of methanol on the Cu 4 X 2 cluster. These calculations further confirmed that the production of CH 3 OH is more favorable than that of CH 4 on the Cu 4 X 2 cluster. In addition, the doping of four O, S, and Se atoms on the Cu 4 cluster has also been tested as an electrocatalyst in CO 2 RR for comparison. The activity of the Cu 4 X 4 cluster will not be explained in detail as it is not improved, and the discussion can be found in the Supplementary Material.  In view of above discussions, the Cu 4 S 2 cluster is found to be an excellent catalyst due to its being more favorable for CH 3 OH production and having low limiting potential in CO 2 RR. Based on the clear understanding of the reduction mechanism of doped clusters, we are committed to investigating the defective graphene supported Cu 4 S 2 cluster as a promising electrode material for CO 2 RR, so as to broaden its practical applications. The geometric structure and stability of the Cu 4 S 2 cluster which is adsorbed on defective graphene in different orientations are first tested in the case of the fully relaxed system. Figure 2 shows the optimized structure and the charge distributions with the lowest-energy binding configuration (denoted as Cu 4 S 2 /SV). The removal of one carbon atom from the graphene sheet leads to the undercoordination of three C atoms around the vacancy, which exerts strong attraction for the Cu atom and allows the Cu 4 S 2 cluster to be stably anchored to the defective graphene. From Figure 2a, the Cu atom at the bottom forms three covalent bonds with the adjacent C atoms on the defective graphene, and the average height of the cluster from the defective surface is 1.61Å. The electrons of the Cu 4 S 2 cluster prefer to transfer to the substrate, resulting in the positive charges on the cluster with the value of 3.89 e. The density of states (in Figure S7) reveals the metallic behavior that makes the structure unique and have good conductivity, and the C and Cu form a new bond corresponding to the overlap of the C p and Cu d orbitals. Moreover, the binding energy of anchoring the Cu 4 S 2 cluster at the single vacancy site of graphene (V C ) is to be defined as: where the E total , E V C , and E cluster represent the total energy of the Cu 4 S 2 /SV system, the graphene with a single vacancy, and the Cu 4 S 2 cluster, respectively. From the calculation result, the binding energy of the Cu 4 S 2 /SV system is +5.07 eV, indicating that the supported Cu 4 S 2 cluster can be stably located on the defective graphene and continue the followup reactions.
(a) (b) Figure 2. The optimized structure (a) and charge density difference plot (b) for the Cu 4 S 2 /SV catalyst. The yellow and blue color regions mean the charge accumulation and depletion, respectively.

Reaction Pathway of CO 2 RR on Cu 4 S 2 /SV
As the initial step of Cu 4 S 2 /SV in CO 2 RR, the first hydrogenation reaction involves a (H + + e -) pair transfer to CO 2 . There are two possible adsorbed species, COOH * or OCHO * , in competition with the H * for HER, as illustrated by Equations (1), (2), and (4), respectively.
In both cases of CO 2 activation, the linear structure of the CO 2 molecule transfers into a bent structure for further CO 2 RR. The calculated reaction free energy of the CO 2 hydrogenation occurs on the O site to form intermediate COOH * (0.24 eV) and is found lower than that on C to form OCHO * (0.40 eV). To examine the selectivity of HER on CO 2 RR, as shown in Figure 3a, the reaction free energy for H * formation is 0.58 eV larger than that of COOH * formation. This means that the Cu 4 S 2 /SV catalyst displays higher CO 2 RR activity as well as lower HER activity. Thus, the CO 2 reduction occurs predominantly before the HER process.
In this part, the reaction mechanism for electroreduction of CO 2 to CH 3 OH, CH 4 , HCOOH, and CO was investigated, beginning with the intermediate COOH * . Among these, the formations of CO and HCOOH are accompanied by the two transfer steps of (H + + e -) pairs, and the reaction pathways are presented in Figure 3b and 3c, respectively. For the Cu 4 S 2 /SV catalyst, the optimum pathway of CO production goes through the pathway CO 2 → COOH * → CO * → CO. The rate-determining step is the hydrogenation of the hydroxyl group to generate the CO * , whose limiting potential is calculated to be −1.43 V. In the process of HCOOH generation, the last step of HCOOH desorption is the rate-determining step with a high free energy change of 1.15 eV. This is due to the stable adsorption of HCOOH * that hinders its desorption from the catalyst surface. The results clearly showed that these high limiting potentials exhibit poor selectivity for producing CO and HCOOH. In addition, the basic steps of all possible intermediates of CO 2 hydrogenation to C1 products on the Cu 4 S 2 /SV catalyst and their corresponding reaction free energies are shown in Figure 3d. There are several thermodynamically feasible pathways for the reduction of CO 2 to CH 3 OH and CH 4 , in which the lowest free energy diagrams are listed in Figure 4. The blue line in Figure 4 shows the lowest energy reaction pathway for producing CH 3 OH accompanied by six (H + + e -) pair transfer processes. The first hydrogenation favors the formation of COOH * (state 1); the next step is to generate an adsorbed HCOOH * (2) through transfer of another (H + + e -) pair. Subsequently, HCOOH * is successively hydrogenated to produce CHO * (3) and CH 2 O * (4). Notably, the fifth step of the proton and electron transfer leads to generation of CH 2 OH * (5), where the free energy of the evolution of CH 2 OH * is 0.66 eV lower than the formation of CH 3 O * . Finally, the CH 2 OH * can be hydrogenated to form CH 3 OH (6). In the pathway of CH 3 OH formation, the ratedetermining step is CHO * → CH 2 O * with a low limiting potential of −0.35 V. The great activity is caused by the strong interaction between Cu 4 S 2 and defective graphene, which could further tune the electronic structure and stability of the supported cluster. The charge transfer from the cluster to the substrate indicates the Cu d orbitals become more vacant and available to be adsorbates. Another red line represents the pathway which takes eight (H + + e -) pairs of transfer steps to produce methane. In Figure 4, the pathway to generate methanol and methane diverges in the sixth step: one is a pair proton-electron transfer to the C end of CH 2 OH * and eventually produces methanol; the other is to form the H 2 O molecule first, and then the remaining CH 2 * undergoes continuous hydrogenations to produce methane at last. Obviously, the step of CH 3 OH formation is thermodynamically downhill, rather than the branch of CH 4 formation requiring the endothermic sort for the next hydrogenation reaction. Hence, the Cu 4 S 2 /SV electrode material is more energetically preferable for CH 3 OH production compared with other C1 products in CO 2 reduction. For the pathway of CH 4 formation, the calculation shows that removing the OH * from the catalyst is the most endothermic step during the whole reaction process. In this case, the step of releasing H 2 O becomes the rate-determining step, of which the limiting potential is −0.51 V. Due to the limiting potential of CH 3 OH formation being lower than that of CH 4 formation, the former may be preferred from a thermodynamic point of view. Based on the above discussions, the calculation results indicate that the Cu 4 S 2 /SV catalyst is not only more favorable for generating CH 3 OH but also can effectively suppress the competing reactions, such as HER, CO, and HCOOH production. In addition, recent experiments [35][36][37] have shown that the supported Cu-based electrode as a electrocatalyst for CO 2 RR has higher CH 3 OH selectivity, suggesting that the free energy changes for generating CH 3 OH are smaller than those of other products. This is consistent with our calculation.

Computational Details
Density functional theory (DFT) [38] calculations were used to investigate the reaction mechanism of isolated Cu 4 X n clusters in CO 2 reduction. The first section in this study is to search for local or even global minima of molecular clusters using ABCluster software [39] which adopts an artificial bee colony algorithm to conduct the global optimization and conformation search. Full details about the global minima nanocluster searching process are available in the Supplementary Material. After that, the subsequent geometry optimizations and vibration frequencies were implemented at the level of B3LYP/6-31g(d, p) [40,41] by the Gaussian16 program [42]. Next, the various active sites for electrochemical CO 2 adsorption and hydrogenation reactions on the doped clusters were tested. Following that, the electronic energies were refined with the single-point energy calculation employing the high-level DLPNO-CCSD(T) method [43], along with the cc-pVTZ basis set [44] using the ORCA software [45]. The Gibbs free energies were then obtained through electronic energies at the DLPNO-CCSD(T)/cc-pVTZ level of theory corrected with the zero-point energies (ZPEs) computed by the B3LYP/6-31g(d, p) method. Besides, the conductor-like polarizable continuum model (CPCM) [46] with toluene as the solvent was used in this study. Finally, the reasonable intermediates were sieved using the calculation mentioned above, after that we could obtain the most feasible pathway for CO 2 reduction for different products on the doped clusters.
The calculations about the defective graphene supported doped cluster were conducted using the Vienna ab initio simulation package (VASP, version 5.3.2) [47], which is based on the DFT with the projector-augmented wave method (PAW) [48]. The PBE exchange-correlation functional [49] and van der Waals (vdW) interactions described via a pair-wise force field by the DFT-D2 method of Grimme [50] were adopted for all the defective graphene-supported system calculations. In this work, the model structure comprises 5 × 5 single-layer graphene unit cells with a single vacancy as the substrate, and the doped cluster was placed above the vacancy site. The vacuum region was set as 15Å for preventing the mirror images effect. The MonKhorst-Pack mesh k-points 4 × 4 × 1 was used for calculation, and a plane-wave energy cutoff is used of 520 eV after convergence test (see Figure S8). All the atoms were allowed to relax during the model optimization process. The total energy converged to 10 -7 eV and the maximum ionic force was less than 0.01 eV/Å. In addition, VASPsol [51] was used to describe the effect of electrostatics, cavitation, and dispersion on interactions between solutes and solvents. The VASPsol model considering the solvent molecule as the continuum solvent, and the dielectric constant were used to represent the solvent effect. Here, the dielectric constant of toluene = 2.37 was used to simulate the implicit solvent environment.
The computational hydrogen electrode (CHE) model was described by Nørskov et al. in 2004 [52], and is combined with energetic from DFT simulations to calculate the free energy changes (∆G) for various basic steps of CO 2 reduction. In the CHE model, the chemical potential of a proton-electron pair is defined as half of the gaseous H 2 at equilibrium, under any pH values and temperatures, the reaction is equilibrated at 0 V, 1325 Pa reversible hydrogen electrode (RHE). The free energy for a basic electrochemical reduction step is defined as Here, " * " denotes an active site on the catalyst. Then the potential-dependent free energy change is thus defined as in which µ is the chemical potential, e is the elementary positive charge and U is the electrode potential versus RHE. The chemical potential is shifted by −eU when an external potential U is applied. While at zero applied potential, ∆G becomes −U L /e, in which U L means the limiting potential in the CO 2 RR process.

Conclusions
To sum up, the reaction mechanism of CO 2 electroreduction to C1 products on a defective graphene-supported Cu 4 S 2 cluster and isolated Cu 4 X n clusters was investigated. DFT calculations indicated the CH 3 OH is the most feasible product among the C1 products on Cu 4 X 2 clusters, and the limiting potentials for producing CH 3 OH are in the order of Cu 4 S 2 < Cu 4 O 2 < Cu 4 Se 2 . The Cu 4 S 2 cluster has excellent selectivity and the best catalytic activity (−0.48 V) during the conversion of CO 2 to CH 3 OH among these isolated clusters. For the Cu 4 S 2 /SV catalyst, defect-engineered graphene could induce a strong interaction between cluster and substrate, which facilitates the improved catalytic performance through adjusting the electronic structure and stability of the Cu 4 S 2 cluster. The Cu 4 S 2 /SV catalyst not only shows superior activity to produce CH 3 OH with a low potential of −0.35 V, but can also significantly suppress other competing reactions, such as HER, CO and HCOOH production. We hope the combination of modification effect and defective support engineering can provide a valuable strategy for the Cu-based catalyst design with highly electrocatalytic performance in CO 2 reduction.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/catal12050454/s1, Figure S1: The ground minimal structures and low-lying energy isomers of the gas-phase Cu 4 X 2 clusters. Figure S2: Global minimum geometry composition and some selected bond lengths of Cu 4 X 2 clusters and Cu 4 X 4 clusters; Figure S3: Mechanistic free energies diagram of electrochemical CO 2 reduction on Cu 4 O 2 cluster, Cu 4 S 2 cluster, and Cu 4 Se 2 cluster with no applied potential; Figure S4: Free energies diagram for producing H 2 , CO, and HCOOH with no applied potential; Figure S5: Kinetic barrier diagram of product-determining step on Cu 4 O 2 cluster, Cu 4 S 2 cluster, and Cu 4 Se 2 cluster; Figure S6: Mechanistic free energies diagram of electrochemical CO 2 reduction on the Cu 4 O 4 cluster, Cu 4 S 4 cluster, and Cu 4 Se 4 cluster with no applied potential; Figure S7: Total density of states (DOS) of Cu 4 S 2/ SV catalyst, Projected density of states (PDOS) of p orbitals of C, and PDOS of d orbitals of Cu; Figure S8: Variation the reaction free energy of COOH * on Cu 4 S 2/ SV catalyst as a function of the cutoff energy value in the plane-wave calculations; Table S1: Geometric parameters , the maximum, and minimum harmonic vibrational frequencies of the GM structure of gas-phase Cu 4 X 2 clusters; Table S2: The charge distribution of elements in the neutral Cu 4 X, Cu 4 X 2 , and Cu 4 X 4 clusters by the natural bond orbital charge analysis; Table S3: The charge distribution of elements in the neutral Cu 4 X, Cu 4 X 2 , and Cu 4 X 4 clusters with electrochemical CO 2 adsorption by the natural bond orbital charge analysis. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study and not reported in the Supplementary Materials are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.