Cu4 Cluster Doped Monolayer MoS2 for CO Oxidation

The catalytic oxidation of CO molecule on a thermodynamically stable Cu4 cluster doped MoS2 monolayer is investigated by density functional theory (DFT) where the reaction proceeds in a new formation order of COOOCO* (O2* + 2CO* → COOOCO*), OCO* (COOOCO* → CO2 + OCO*), and CO2 (OCO* → CO2) desorption with the corresponding reaction barrier values of 0.220 eV, 0.370 eV and 0.119 eV, respectively. Therein, the rate-determining step is the second one. This low barrier indicates high activity of this system where CO oxidation could be realized at room temperature (even lower). As a result, the Cu4 doped MoS2 could be a candidate for CO oxidation with lower cost and higher activity without poisoning and corrosion problems.

Scientific RepoRts | 5:11230 | DOi: 10.1038/srep11230 for CO oxidation 8 . This structure is stable due to the strong chemical bonds among Cu 4 and monolayer MoS 2 25 while the triangular Cu 3 active site is acted for CO oxidation by adsorbing O 2 and more CO molecules and the Cu 4 cluster completely inserts into the sandwich of MoS 2 . We have found a new OCO* intermediate state (MS) with a small E bar of 0.370 eV on Cu 4 cluster during the CO oxidation process. Our calculations suggest that the Cu 4 cluster embedded in a monolayer MoS 2 is a good candidate for CO oxidation.

Results and discussion
The experimental fabrication of the Cu 4 doped MoS 2 could be intricate since there are two different doping sites on monolayer MoS 2 , Mo vacancy and S vacancy 40,41 . It is known that Re atoms and Co atoms occupying Mo sites in monolayer MoS 2 have been synthesized by the chemical vapor transport (CVT) and chemical vapor deposition (CVD) method, respectively 42,43 . Since the value of Pauling electronegativity of Cu (1.90) is almost the same of Re (1.90) and Co (1.88), the Cu atom could substitute Mo atom in the MoS 2 by means of CVT or CVD, which has been proved to be feasible through the density functional calculations 25 . Then, the S vacancies are prepared by low-energy argon sputtering or electron irradiation. Last, the Cu 3 atoms are embedded into these vacant sites through the physical vapor deposition. The corresponding structure of Cu 4 doped MoS 2 is shown in Fig. 1(a) where one Cu atom substituting the Mo atom is denoted as Cu 1 , and the other three Cu atoms replacing three S atoms on the surface of monolayer MoS 2 are named as Cu 3 . The flat triangular Cu 3 active site on the surface plays a vital role for CO oxidation 8 . The Cu 3 atoms and surface S atoms are at the same plane, making it less active to be oxidized than the metal atoms above the graphene surface. The bond lengths of Cu 3 -Cu 3 and Cu 3 -Cu 1 are 2.53 Å and 2.57 Å, being almost the same of Cu-Cu (2.55 Å) in Cu bulk. The bond lengths of Cu 3 -Mo and Cu 1 -S are 2.66 Å and 2.28 Å. The latter is shorter than the bond length of Cu-S of 2.41 Å in Cu-doped MoS 2 25 . Thus, Cu 1 -S bond is stronger than Cu-S bond. From Hirshfeld charge analysis, the electron transfer of Cu 1 and Cu 3 are 0.146 e and 0.076 e, respectively. The direction of electron transfer is in agreement with the values of Pauling electronegativity of Cu (1.90), S (2.58) and Mo (2.16) 44 . There are about 0.374 e transfer from the Cu 4 cluster to monolayer MoS 2 . The electron transfer can also be verified by the charge density difference (CDD) for Cu 4 doped MoS 2 . As shown in Fig. 1(b), the blue and red regions represent the areas of electron accumulation and depletion, respectively. Obviously, different electron affinities of Cu, S and Mo determine the electron distribution. The pronounced charge density The partial density of states (PDOS) projected on the 3d orbitals of Cu 1 , the 2p orbitals of its neighboring S, 3d orbitals and 4s orbitals of Cu 3 and 4d orbitals of its neighboring Mo are plotted, as shown in Fig. 1(c),(d). The strong interaction between Cu 1 and S can be further confirmed by the PDOS in Fig. 1(c) where significant hybridization between Cu 1 -3d and the adjacent S-2p is present, denoting the stability of the system. Furthermore, the hybridization between 4s and 3d orbitals of Cu 3 and the adjacent Mo-4d is also found in Fig. 1(d).
To gain more insight into the stability of Cu 4 doped MoS 2 , we also calculated the E b of Cu 4 doped MoS 2 , being 3.262 eV, showing strong interaction between Cu 4 cluster and its neighboring S and Mo atoms. For the possible diffusion problem of Cu atom, we computed the energy barrier and reaction energy of Cu diffusion. Since the Cu 1 strongly bonds with three S atoms and it is in the middle layer of the Cu 4 doped MoS 2 , the diffusion of the Cu 1 atom is difficult. Therefore, only the surface Cu 3 atoms are considered here. Because all three Cu 3 atoms are the same in symmetry, we study only one Cu 3 atom diffusing to the neighboring hollow site or Mo top site [ Fig. 2(a)]. The diffusion energy barriers for the two cases are the same with a value of 1.41 eV and the reaction energies are 1.02 eV and 1.22 eV, respectively. Considering the high diffusion barriers and the endothermic reactions, the diffusion of Cu 3 is absent and the Cu 4 doped MoS 2 system thus is an energetically stable structure.
To further prove the stability of the system, first principle molecular dynamics (MD) simulation at a constant temperature of T = 500 K in the NVT ensemble (i.e., constant particle number, volume and temperature condition) has been carried out for 5 ps with the time step of 1 fs. Three structures from MD calculation are present in Fig. 2(b). It is found that Cu 4 cluster is fixed in the vacancies of MoS 2 and the Cu atoms are located at the original sites after 5000 dynamics steps at 500 K. Thus, the stability of the studied Cu 4 doped MoS 2 system at room temperature is expected.
The adsorption process of O 2 molecule on Cu 4 doped MoS 2 is considered for possible side-on and end-on configurations. With the former configuration, two found adsorption structures of O 2 molecule are shown in Fig. 3    For CO oxidation reaction, the adsorption of CO and CO 2 also should be considered. In this system, the E ad-CO values are − 1.105 eV on top site,− 0.990 eV on bridge site and − 0.957 eV on hollow site, respectively. Thus, CO prefers adsorbing on the top site, as illustrated in Fig. 3(a). Since E ad-O2 (− 1.749 eV) is stronger than E ad-CO (− 1.105 eV), O 2 is preferentially adsorbed on the hollow site, which indicates there is no CO poisoning problem. And the E ad-CO2 value is − 0.140 eV, which proves that the CO 2 is easy to leave the surface. Note that as long as the first O 2 is adsorbed on the hollow site, sites for other O 2 adsorptions are absent. Because the adsorption energy value of two O 2 on the catalyst is − 1.709 eV, which is weaker than that of only one O 2 . As a result, the Cu 4 doped MoS 2 without poisoning and oxidation problems could be a good catalyst with a long cycle life. To further more comprehensively understand the adsorption of O 2 and CO on the catalyst, the other sites near the Cu 4 cluster are considered. Because MoS 2 surface is inert at the basal planes terminated by S atoms, both E ad-O2 value and E ad-CO value are about − 0.1 eV, which imply that O 2 and CO only are adsorbed on Cu 4 cluster of the catalyst.
All configurations and their E ad-nCO values are given in Table 1. For the co-adsorption of O 2 + nCO, the E ad-nCO are − 0.401 eV,− 0.350 eV and − 0.186 eV respectively for n = 1, 2, 3. The corresponding PDOS of Cu n and Cu n -CO are shown in Fig. 5. When n = 1 and 2, the orbitals below Femi level are away from the Femi level, denoting more stable states. If n = 3, however, the orbitals change less, denoting weaker interaction of the third CO with Cu atom. To simplify the latter discussion, we neglect the adsorption of the third CO on Cu 4 cluster in the following.
Now we begin to consider CO oxidation on Cu 4 doped MoS 2 . It is well known that there are two mechanisms for CO oxidation: Langmuir-Hinshelwood (LH) mechanism and Eley-Rideal (ER) mechanism. The former involves all the reacting intermediates on the surface, whereas the latter does species from the direct reaction with a surface intermediate 46,47 . Both will be discussed in details in the following.
Firstly, the LH mechanism of the one CO co-adsorption with O 2 molecule is considered, which is denoted as mLH shown in Fig. 6 (top views in Fig. S1). The O atom in O 2 molecule approaches the CO molecule and bonds with the C atom of the CO molecule to form the OCOO* intermediate state (MS1 in Fig. 6) with the E bar = 0.302 eV. Then, the O-O bond breaks and CO 2 molecule is released from the catalyst with E bar = 0.292 eV. Therein, the rate-determining step for the mLH mechanism is the formation of OCOO* intermediate (TS1 in Fig. 6). On the other hand, for mER mechanism, the first un-adsorbed CO molecule directly reacts with the activated O 2 molecule. Due to the presence of two adsorption configurations of O 2 molecule, there are two kinds of mER mechanisms in the reaction path (mER in Fig. 6). The migration of CO molecule toward the pre-adsorbed O 2 molecule is determined as the rate-determining step with E bar = 0.385 eV (TS3 and TS4 in Fig. 6).
When the O 2 molecule and first CO molecule have been adsorbed on the triangular Cu 3 active site, the second CO molecule can be further adsorbed, which reacts with the adsorbed O 2 . We assume that the reaction could follow bLH mechanism or bER mechanism, we define b denoting the case where two CO molecules are involved in the reaction path. In the bLH mechanism, the adsorption structures of the two CO molecules and O 2 molecule are shown in Fig. 3(e),(f). The initial structure of Fig. 3(f) Fig. 6). Their E bar values are 0.364 eV and 0.370 eV, respectively. Although both values are pretty much the same, the product of FS5, OCO*, releases 0.22 eV more energy for the formation of first CO 2 than that of FS4 while OCO* is beneficial for the next oxidation reaction. As a result, FS5 tends to happen relative to FS4. In the case of bER mechanism (Fig. 6), E bar = 0.381 eV for the release of the first CO 2 .
After releasing the first CO 2 , the following structures are shown in Fig. 6(FS1~FS6). In addition to OCO* (FS5), others have one O atom on the triangular Cu 3 site. The O atom needs larger E bar to form CO 2 with CO due to the strong interaction between the Cu 3 atoms and O atom, as shown in Fig. S2. However, the formation of CO 2 from the OCO* intermediate state has the smallest E bar (0.119 eV) to escape from the surface of the catalyst, implying desorption of the second CO 2 .
The above results show that the LH mechanism is better than the corresponding ER mechanism as the O 2 molecule itself is not activated enough without the cooperative adsorption of CO. Then, the co-adsorbed CO molecules affect the E bar while the rate-determining step changes from the first step to the second one. The E bar of ER mechanism is affected by the adsorbed number of CO molecule too. Last, it should be noted that in the bLH mechanism, the OCO* (FS5 in Fig. 6) is found as the last product which is particularly favorable for the release of the second CO 2 with E bar = 0.119 eV. As a result, the most optimal reaction path is given in Fig. 7. Among three E bar values in the reaction path, the largest E bar value is 0.370 eV where the rate-determining step is the formation of the OCO* intermediate state.
Thus, no matter from the point of view of dynamics or thermodynamics, the reaction path (Fig. 7) is the most optimal process. Because of the complexity and diversity of the experimental conditions, there are a few other reaction paths in CO oxidation reaction. All kinds of reaction paths and their E bar values are shown in Fig. S3 and Fig. S4. Now we compare the corresponding E bar values at the rate-determining reaction step for our system and related systems which are listed in Table 2. As shown in Table 2, Zn-embedded graphene and Au-embedded graphene have smaller E bar than this system. However, it is noteworthy that our system is more stable than Zn-embedded graphene while Au-embedded graphene has CO poisoning problem, thus their applications are limited. Thus, the overall performance of Cu 4 doped MoS 2 for the CO oxidation should be the best among the considered systems.
Under the biggest doping concentration, Cu 4 doped 3 × 3 monolayer MoS 2 supercell is also considered for CO oxidation. For CO and O 2 adsorption, their adsorption energy values are shown in Table 1 where the both supercells have almost the same values. The stability of 3 × 3 supercell is also determined, as shown in Fig. S5. The results show that even Cu 4 doped monolayer MoS 2 has the biggest doping concentration, it still possesses good catalytic activity and stability for CO oxidation as our results from the 4 × 4 supercell studied above.
In summary, our comprehensive DFT studies of CO oxidation on Cu 4 doped monolayer MoS 2 suggest that the protruded triangular Cu 3 site is the main active site for CO catalytic oxidation while the number of CO adsorbed molecule produces a significant effect on the energy barriers of the CO oxidation reaction. During the reaction, an OCO* intermediate state is found, which leads to the energy barrier of CO oxidation of 0.370 eV. As a result, the Cu 4 doped monolayer MoS 2 with outstanding catalytic activity without poisoning and oxidation problems could be a good candidate for CO oxidation with low cost and high activity.

Methods
In this work, all calculations are performed using the spin-unrestricted density functional theory (DFT) as implemented in the DMol 3 code 48 . Exchange-correlation functions are taken as a generalized gradient approximation (GGA) with Perdew-Wang correlation (PWC) 49 . DFT semi-core pseudo potentials (DSPPs) core treatment 50 is implemented for relativistic effects, which replaces core electrons by a single effective potential. In addition, double numerical plus polarization (DNP) is chosen as the basis set and the quality of orbital cutoff is fine. The convergence criteria of the geometrical optimization are set to be 1.0 × 10 −5 hartree for the energy change, 2.0 × 10 −3 hartree/Å for the gradient, and 5.0 × 10 −3 Å for the displacement, respectively. The smearing parameter is set to be 0.005 hartree in the geometric optimization. For transition states (TS) searching, the calculation firstly performs a linear synchronous transit (LST) 51 maximum, which is followed by an energy minimization in directions conjugating to the reaction pathway. TS approximation obtained via LST/optimization is then used to perform a quadratic synchronous transit (QST) 51 maximization to find more accurate transitional states. The convergence tolerance of the root mean square (RMS) force is 2.0 × 10 −3 hartree/Å and the maximum number for QST step is set as 10. In the simulation, three-dimensional periodic boundary conditions are taken. The simulation cell consists of a 4 × 4 monolayer MoS 2 supercell with a vacuum width of 18 Å, which leads to negligible interactions between the system and their mirror images. In order to prove the effect of doping concentration, 3 × 3 monolayer MoS 2 supercell is also considered. For geometric optimization and the search for the transition state (TS), the Brillouin zone integration is performed with 3 × 3 × 1 k-point sampling. After structure relaxations, the density of states (DOS) are calculated with a finer k-point grid  where E(MoS 2 ), E(Cu) and E(Cu 4 /MoS 2 ) are the total energies of the monolayer MoS 2 with three S vacancies and a Mo vacancy, the free Cu atom, and the Cu 4 doped MoS 2 , respectively. For one molecule (CO, O 2 , CO 2 ) adsorbed on catalyst, the adsorption energy values of E ad-M (the subscript M denotes the corresponding molecule) are determined by,