Computational screening study of double transition metal carbonitrides M′2M″CNO2-MXene as catalysts for hydrogen evolution reaction

Two-dimensional (2D) transition metal carbonitrides (MXene) have attracted growing interest in electrocatalytic hydrogen production due to its structural and electronic properties. In this work, the hydrogen evolution reaction (HER) activity of all 64 O-terminated ordered double transition metal carbonitrides in the form of M′2M″CNO2 (M′ = Ti, V, Cr, Zr, Nb, Mo, Hf, Ta; M″ = Ti, V, Cr, Zr, Nb, Mo, Hf, Ta) has been investigated by well-defined density functional theory (DFT) calculations. The results indicate that there are 11 M′2M″CNO2-MXene candidates whose HER performance is superior to that of Pt. Moreover, according to the stability screening, it is proved that Ti2NbCNO2, Mo2TiCNO2, and Ti2VCNO2 are more stable than other candidates. Especially, Ti2NbCNO2 have the potential to be perfect HER catalyst with the small Gibbs free energies of hydrogen adsorption (ΔGH) value of 0.02 eV, abundant catalytic sites on the C-side, and better stability. This work paves the way on designing excellent HER catalyst candidates based on M′2M″CNO2-MXenes.


INTRODUCTION
Electrochemical water splitting through hydrogen evolution reaction (HER) associated with catalysis is an attractive method for hydrogen production 1,2 , but it is still challenging due to the lack of effective and low-cost HER catalyst. Noble metals such as platinum (Pt) have been verified as excellent HER catalysts, but the corresponding scarcity and high price limit its large-scale applications. Therefore, finding a cost-effective alternative to platinum is crucial. Computational screening is an effective method for exploring catalysts. First-principles calculations can be used to evaluate the performance of a mass of potential candidates before selecting catalysts suitable for experimental synthesis. In comparison to traditional materials, two-dimensional (2D) materials have more superiorities as HER catalysts. The large specific surface area of 2D materials can provide abundant active sites, and its diverse chemical properties also supply convenience for adjusting its catalytic performance. In recent years, some efficient and 2D HER catalysts, including MoS 2 3-5 , metal phosphides 6 , g-C 3 N 4 7-9 , and chalcogenides 10 , have been extensively studied. However, the poor charge transfer performance 3 and restricted active sites 4,11 are major defects of these catalysts. Wherefore, it is still necessary to search efficient and economical 2D HER catalysts.
In recent years, 2D materials involving transition metal carbides, nitrides, and carbonitrides known as MXene have been widely studied in both experiment [12][13][14][15][16] and theory [17][18][19][20] . MXenes usually have a general formula M n+1 X n T x , where M denotes the early transition metals, X represents C and/or N, T x represents surface functional groups (O, OH, or F), and n takes values from 1 to 3. Owing to its electronic properties, excellent stability, hydrophilicity, extensive chemical and structural adjustments, MXenes have possessed broad application prospects in nanodevices 21,22 , hydrogen storage 23,24 , supercapacitors 25 , photolytic or electric catalysis 26,27 and other fields. At the same time, MXenes are also effective HER catalysts, and have recently received extensive attention in electrocatalytic water splitting to produce hydrogen. Huang et al. 28 29 . Therefore, we tried to construct "sandwich" structure of transition metal carbonitrides M 3 CN, that is, double transition metal carbonitrides M′ 2 M″CNO 2 -MXene to improve its HER activity.
In this work, theory studies based on well-defined DFT calculations were used to systematically investigate the HER activity and geometric structure of the O-terminated double transition metal carbonitrides M′ 2 M″CNO 2 -MXene (M′ = Ti, V, Cr, Zr, Nb, Mo, Hf, Ta; M″ = Ti, V, Cr, Zr, Nb, Mo, Hf, Ta), as well as the stability of some selected MXene candidates. By high-throughput 1 computational screening, 11 M′ 2 M″CNO 2 -MXene candidates whose HER performance is even better than that of Pt were found. Further, the catalytic active sites and stability of the selected candidates were taken into consideration. Ultimately, Ti 2 NbCNO 2 , Mo 2 TiCNO 2 , and Ti 2 VCNO 2 were identified as highly efficient and more stable HER catalysts, among which Ti 2 NbCNO 2 with the ΔG H value of 0.02 eV and a great deal of catalytic sites on the C-side is the most ideal HER catalyst candidate.

RESULTS AND DISCUSSION
The bare structure of ordered double transition metal carbonitrides Figure 1a shows the geometric structure of the studied bare ordered double transition metal carbonitrides. The M′ layers are located in the outermost layer of the MXenes monolayer, and M′′ layer fills the middle layer of MXenes, and the carbon layer and nitrogen layer are distributed on both sides of M′′ layer. In other expression, the atomic layer of another type of transition metal element was used to replace the original intermediate layer of the M 3 CN MXenes 28 . Due to the high chemical activity of the bare MXenes monolayer surface, the generated MXenes are easily functionalized by surface terminal groups (O, OH, and F) 30 . The bare MXenes have been proved to be thermodynamically prone to be functionalized with oxygen terminals 31 . Therefore, we only focus on the properties of O-terminated MXenes in the subsequent researches.

The O-terminated structure of ordered M′ 2 M″CNO 2 -MXenes
For O-terminated MXenes, in order to maximize the coordination number of the outermost transition metals, the oxygen groups are normally located in the hollow position of the outermost layers. The basal planes of bare MXenes have two different hollow sites, denoted as fcc and hcp (Fig. 1b). Therefore, there are four different structures for M′ 2 M″CNO 2 -MXene (Fig. 1c), marked as fcc-fcc, fcc-hcp, hcp-fcc, and hcp-hcp. Based on the calculation results, we selected the configurations of M′ 2 M″CNO 2 with the lowest energy as the most stable structures for further research, as shown in Fig Table 1; the data are in good agreement with the previous results 28 .
HER activity screening of M′ 2 M″CNO 2 Figure 3 represents the flow chart for selecting the most ideal HER catalyst from all 64 M′ 2 M″CNO 2 -MXenes. As we know, the Gibbs free energy of hydrogen adsorption is usually used as the indicator of the HER activity 32 and defined as where ΔE H represents the adsorption energy of hydrogen, ΔE ZPE and ΔS H are the difference in zero point energy (ZPE) and the entropy between the hydrogen adsorbed state and gas phase, respectively. The ΔE H is defined as where E(MXene + nH) denotes the energy of MXene with n adsorbed hydrogen atoms, and E(MXene + (n−1)H) represents the energy of MXene with n−1 adsorbed hydrogen atoms. E(H 2 ) is energy of H 2 in the gas phase. The ΔE ZPE is computed by  (Fig. 4).
Notably, when the second transition metal element is introduced to form M′ 2 M″CNO 2 , the HER activity of MXenes will be greatly improved. Obviously, the M′ compared with M′′ have a greater impact on the HER activity. The outermost transition metal M′ is Ti or Nb, most of M′ 2 M″CNO 2 -MXenes exhibit desired HER activity on the C-side, while the M′ is Ta, the HER activity is better on the N-side. Overall, the number of shadow region of the C-side is more than that of N-side, indicating that the C-side of M′ 2 M″ CNO 2 -MXenes is more active for HER. It is well known that ΔG H is connected with the bonding energy between H atoms and the surface of MXenes. The bonding energy is relevant to the electronegativity of all atoms of MXenes (such as M′, M″, C, and N atoms), which results in different electronic structure for the Cand N-side of M′ 2 M″CNO 2 . Therefore, it is hard to summarize the regularity of HER activity for whole M′ 2 M″CNO 2 -MXenes.
The volcano curve can be used to intuitively characterize the HER activity of electrocatalysts. The abscissa of the volcano curve is the ΔG H while the ordinate is the overpotential (η r = ±ΔG H /e − ). The closer the coordinate is to the top of the volcano, the more perfect HER performance of the catalyst. We divide the volcano curve into two branches, the points on the left branch represent the strong adsorption between H atom and the catalyst, while the  right branch is the opposite. Moreover, the volcano diagram is employed to characterize the HER activity of the C-and N-side of all 64 M′ 2 M″CNO 2 -MXenes, as shown in Fig. 5a, b. Here some representative HER catalysts are marked in Fig. 5c, d, such as Pt 33 , MoS 2 (ref. 34 ), WS 2 (ref. 35 ), C 3 N 4 (ref. 36 ). The profiles indicate that there are 11 M′ 2 M″CNO 2 -MXenes on the C-or N-side have the better electrocatalytic HER activity than Pt, for the C-side such as Ti 2 NbCNO 2 , Ti 2 VCNO 2 , Ti 2 HfCNO 2 , Mo 2 TiCNO 2 , Ti 2 ZrCNO 2 , Ti 2 TaCNO 2 , and Nb 2 TaCNO 2 , and for the N-side such as Ta 2 NbCNO 2 , Ta 2 TiCNO 2 , Ta 2 ZrCNO 2 , and Zr 2 CrCNO 2 . The calculated |ΔG H | ranges from 0.071 to 0.003 eV at 1/4 coverage. Furthermore, the coordinate for the N-side of Ta 2 NbCNO 2 with the ΔG H value of 0.003 eV is very close to the peak of the volcano, meaning that very little overvoltage is required by this material during the HER electrochemical reaction, which indicates its excellent HER activity. Of course, the rest of M′ 2 M″CNO 2 appeared in Fig. 5c, d also have quite good prospects as HER catalysts, but we will not do much research on them in the follow-up work.
Furthermore, for these 11 M′ 2 M″CNO 2 -MXenes we preliminary screened, the effect of H coverage on the HER activity is investigated. The C-or N-side of the 2 × 2 × 1 M′ 2 M″CNO 2 monolayer has four O atoms that can be used as hydrogen  adsorption sites. As shown in Fig. 6a, taking the C-side of Ti 2 ZrCNO 2 as an example, the H coverage (θ) of the C-side of Ti 2 ZrCNO 2 is 1/4, 1/2, 3/4, and 1, respectively. The relationship between the ΔG H and H coverage for the 11 M′ 2 M″CNO 2 -MXenes has been shown in Fig. 6b. ΔG H presents an overall upward tendency with the increase of H coverage. The points in the shadow (|ΔG H | < 0.2 eV) stand for the MXenes with fit coverage considering to have outstanding HER activity. It can also be understood as the HER active sites of MXenes. Most of the M′ 2 M″ CNO 2 -MXenes we screened show HER activity only when the H coverage is 1/4. It is worth noting that the C-side of Ti 2 NbCNO 2 and Ti 2 ZrCNO 2 still have HER activity when the H coverage is 2/4, which indicates that there are more HER active sites for them.

Thermodynamic and dynamic stability of M′ 2 M″CNO 2 -MXenes
The stability of the materials is the prerequisite for its synthesis and application. Therefore, the stability of 11 M′ 2 M″CNO 2 -MXenes we have screened is investigated in order to prove that our research objects can be synthesized experimentally. Firstly, the formation energy (E formation ) of 11 M′ 2 M″CNO 2 -MXenes is calculated, which could be obtained by subtracting the energy of the 2D MXene from the energy of the most stable competitors with the same number of atoms and stoichiometry 37 .
As shown in Table 1, there are 3 M′ 2 M″CNO 2 -MXenes with negative formation energy, namely Ti 2 NbCNO 2 , Ti 2 VCNO 2 , and Mo 2 TiCNO 2 , which indicates that these candidates are thermodynamically stable. Furthermore, to check the dynamic stability of the aforementioned 3 M′ 2 M″CNO 2 -MXenes, the canonical ensemble (NVT) ab initio molecular dynamics (AIMD) simulation is performed at 300k for 10 ps, as shown in Fig. 7. The energy of the MXenes fluctuate slightly above and below the average value, and the crystal structures have no obvious deformation (Fig. 7d), indicating the prominent dynamic stability of them. Therefore, we believe that these 3 M′ 2 M″CNO 2 -MXenes are stable and can be synthesized experimentally.
In summary, the HER activity of all 64 M′ 2 M″CNO 2 -MXenes has been performed by well-defined density functional theory calculations. The crystal structure, Gibbs free energy of hydrogen adsorption, and the stability of several MXene candidates are systematically investigated. According to the calculated ΔG H , the outermost transition metal M′ is Ti or Nb, the C-side for most of M′ 2 M″CNO 2 -MXenes exhibit the desired HER activity while the M′ is Ta showing better HER activity on the N-side. In addition, the  In addition, to verify the stability of Ti 2 NbCNO 2 in an acidic environment, the NVT ensemble AIMD simulation is performed at 300 K for 4 ps. The energy oscillate slightly near the equilibrium configuration, and no notable geometric deformation is found ( Supplementary Fig. 1), suggesting the excellent stability of Ti 2 NbCNO 2 in an acidic environment.

Computational details
All calculations were carried out by employing the Vienna Ab initio Simulation Package (VASP) 38 . The projector augmented wave (PAW) 39 method was employed to explain the interaction between ions and valence electrons, and the exchange-correlation interaction of electrons was handled by the generalized gradient approximation 40 with Perdew-Burke-Ernzerhof 41 potential. All calculation parameters were selected after the convergence test of the total energy of the system. The cut-off energy of plane wave was set to be 520 eV, the energy and force convergence were selected as 10 −5 eV and −0.01 eV/Å, respectively. The Grimme's DFT-D3 (ref. 42 ) method was exerted to effectively rectify the van der Waals interactions. The Brillouin zone of 5 × 5 × 1 k-points mesh created by Monkhorst Pack Scheme were used for the structural optimization, and the Gamma point (1 × 1 × 1) was employed for the calculation of the zero point energy and AIMD simulations. The 2 × 2 × 1 supercell was applied in computation screening research of HER activity for all MXenes monolayers, and the 4 × 4 × 1 supercell was taken for the AIMD simulations. In order to avoid interference between MXenes monolayers, at least 15 Å vacuum layer was required in the z-axis direction of all MXenes structural models.

DATA AVAILABILITY
Data will be available from the corresponding author upon reasonable request.