Quasiclassical trajectory study of the C(¹D) + HD reaction

Abstract

The isotopic product CD/CH branching ratios, thermal rate coefficients, as well as other dynamical quantities of the C(¹D) + HD → CD(H) + H(D) reaction are investigated by detailed quasiclassical trajectory calculations on the highly accurate singlet ground-state (ã¹A′) and the first excited-state (b̃¹A′′) global ab initio potential energy surfaces (PESs) recently constructed by us. The calculated CD/CH branching ratios are in reasonable agreement with experiment. The thermal rate coefficients in the temperature range of 200–1500 K are calculated, and the obtained values at room temperature are in very good agreement with available experimental data. The distinct topographical features between the present and previous PESs, which influence the CD/CH branching ratio, are also discussed. In addition, the effect of the b̃¹A′′ PES is investigated, and the results show that the contribution from the b̃¹A′′ PES to the total reactivity is rather noticeable.

---

1 Introduction

The isotopic product DX/HX branching ratio for the X + HD (X is an atom or radical) reaction serves as a very sensitive probe of the reaction dynamics because of the large mass disparity of the H and D atoms, and it has been used to acquire understanding of the chemical reaction dynamics for many years. For example, a combined theoretical and experimental investigation of the DCl/HCl branching ratio of the Cl + HD reaction has revealed that the van der Waals (vdW) forces in the entrance valley play a decisive role to the outcome of this reaction at low collision energies, and it has become widely recognized that the weak vdW forces could have substantial effects on the reaction dynamics.^1,2

Unlike the Cl + HD reaction that proceeds with an activation barrier, the C(¹D) + H₂ reaction is characterized by a deep potential well and regarded as a prototype in the study of the dynamics of the complex-forming reactions.^3–5 There are a number of kinetic and dynamical studies for the C(¹D) + H₂ and C(¹D) + D₂ reactions,^6–16 but the C(¹D) + HD reaction has not been sufficiently studied. In 1997, Sato et al.¹⁷ reported the product CD/CH branching ratio of 1.61 ± 0.10 for the title reaction at the collision energy of 3.7 kJ mol⁻¹ and measured the thermal rate coefficients using laser-induced fluorescence (LIF) technique. There was also an earlier experiment by Fisher et al.¹⁸ that reported a value of 1.69 ± 0.07 for the product CD/CH branching ratio by using LIF technique at collision energies less than 4.2 kJ mol⁻¹. To our knowledge, the detailed dynamical experimental investigations such as crossed molecular-beam (CMB) experiments are still lacking for the C(¹D) + HD reaction. On the theoretical side, several research groups^19–22 studied the title reaction using a global ab initio potential energy surface (PES) for the singlet ground state (ã¹A′) constructed by Bussery-Honvault et al. (BHL PES)²³ and its modified version (RKHS PES).²⁴ In particular, Defazio et al.¹⁹ obtained the cross sections, isotopic branching ratios and rate coefficients using the real wave packet (WP) method, Kang et al.^20,21 reported the product CD/CH branching ratio as well as the vector correlation between products and reagents using the quasiclassical trajectory (QCT) method, and Lin et al.²² investigated the state-to-state integral cross sections (ICSs) and thermal rate coefficients by the statistical model method. In 2011, Varandas and coworkers carried out QCT calculations on a double many-body expansion (DMBE) PES²⁵ constructed by themselves and obtained the rate coefficients as well as isotopic product branching ratios for the title reaction.²⁶ In 2013, Sun et al.²⁷ performed quantum wave packet calculations on a global ab initio PES for the C(¹D) + H₂(D₂) reaction, which yielded rate coefficients in good agreement with experiment. Recently, a highly accurate global ab initio PES for the ã¹A′ state has been reported by our group,²⁸ referred to as the ZMB-a surface. This surface is unique in the accurate description of the regions of vdW interactions and around conical intersections (CIs). The subsequent accurate quantum dynamics²⁹ (QD) and QCT³⁰ calculations on this surface for the C(¹D) + H₂ reaction yielded some interesting results. Besides the singlet ground state, the first excited (b̃¹A′′) state is also expected to play an important role.³¹ Up to now, only Defazio et al.¹⁹ have investigated the contribution of the b̃¹A′′ state to the overall ICSs and rate coefficients for the C(¹D) + HD reaction using the BJHL PES constructed by Bussery-Honvault et al.³² Most recently, a highly accurate global ab initio PES for the excited b̃¹A′′ state has been constructed by our group.³³ Accurate QD calculations on our ã¹A′ and b̃¹A′′ surfaces yielded various dynamical quantities for the C(¹D) + D₂ reaction in very good agreement with results from crossed-beam experiments, which was not achieved before on the PESs constructed by other groups. More importantly, it was found³³ that the weak long-range forces cause vdW saddles in this kind of complex-forming reaction which have very different dynamical effects from vdW wells at low collision energies, and the importance of the excited b̃¹A′′ PES was also revealed. In addition, a general definition of vdW saddle in a chemical reaction was given in that work,³³ which is expected to become a useful concept in the study of complex-forming reactions. However, the dynamics for the title reaction has not been studied on the newly built and more accurate PESs. To gain more insights into the dynamics of the C(¹D) + HD reaction, detailed QCT calculations are performed on our ã¹A′ and b̃¹A′′ surfaces in this work.

Meanwhile, despite the fact that it is affordable to perform full-dimensional accurate quantum dynamics calculations at the state-to-state level for tri-atomic reactive systems nowadays, the QCT method remains a valid and efficient tool for the study of reaction dynamics,^34–37 which can give good approximation to various scattering dynamical properties and is useful for dissecting the various reaction mechanisms. In many cases, however, the QCT calculations with the conventional histogram binning (HB) approach suffer from lacking the constraint of the zero-point energy (ZPE) in each vibrational mode of the products. To deal with this problem, the Gaussian binning (GB) type approaches^38–41 have been proposed, where each reactive trajectory is assigned to a weight coefficient according to a Gaussian function centered on the integer vibrational action of the products. Because the C(¹D) + HD reaction has two isotopic product channels, the present work also provides us a good chance to verify the effectiveness of the GB approach by comparing the HB and GB results.

In the present work, detailed QCT calculations for the C(¹D) + HD reaction are performed on the most recent ã¹A′ (ZMB-a) and b̃¹A′′ ab initio PESs. Various dynamical and kinetic quantities are obtained with both the HB and GB approaches, and compared with the available experimental measurements as well as the previous theoretical results. Moreover, the important effect of the PES topographical features and the contribution of the b̃¹A′′ PES are analyzed. This paper is structured as follows. In the next section, the details for the QCT calculations are elucidated; then the results are presented and discussed; the conclusions are given in the last section.

2 Methods and computational details

The QCT method has been described in detail elsewhere,^42–45 here we just give some details. The present QCT calculations were performed by using a modified version of the VENUS96 code^46,47 customized to incorporate our ã¹A′ and b̃¹A′′ PESs.^28,33 All trajectories were initiated at the C(¹D) + HD asymptote with a separation of 8.0 Å, and a time step of 0.02 fs was chosen, which guarantees an energy conservation of better than 1 in 10⁶. The maximum impact parameter b_max was estimated at different initial conditions by a scheme similar to that described in our previous work.³⁰ In this paper, we firstly calculated batches of 10 000–15 000 trajectories at collision energies of 0.01, 0.02, 0.03, 0.038, 0.05, 0.08, 0.1, 0.2, 0.3, 0.4 and 0.5 eV with the initial HD molecule set on its ground rovibrational state (v = 0, j = 0), where the impact parameter b is randomly sampled between 0 and b_max with the b_max values of 5.9, 5.0, 4.5, 4.3, 3.9, 3.4, 3.2, 2.7, 2.5, 2.4 and 2.3 Å, respectively; we also calculated batches of 10 000 trajectories at the temperatures of 200, 298, 500, 1000, 1500 K with b randomly sampled between 0 and 7.0 Å, in which the relative translational energy was randomly selected to mimic the Boltzmann distribution and the initial HD was set on different rovibrational states related to the rovibrational Boltzmann distribution at the given temperature. In the GB approach, the full-width-half-maximum (fwhm) of the Gaussian function was set to 0.1.

The ICS is derived by the traditional expression:

where N_r is the number of the reactive trajectory, N_t is the number of the total trajectory, and w_i is the Gaussian weight of the i-th reactive trajectory in the GB approach and is set to 1 in the HB approach.

The thermal rate coefficient can be calculated from

where g_e = 1/5 is the electronic degeneracy factor, k_B is the Boltzmann constant, μ is the reduced mass of the reagents, E_vj is the rovibrational energy of the HD (v, j), Q_vj is the corresponding partition function, σ(E_rel,v,j) is the reactive cross section, and E_rel is the relative translational energy.

The DCS is obtained by the method of moment expansion in Legendre polynomials with the expression of^42,43

where P_n(cos θ) is the n-degree Legendre polynomial, Ω is the solid angle (dΩ = 2π sin θdθ), and θ is the scattering angle defined as the angle between the velocity vector of the product H/D and that of the reagent HD in the center-of-mass system (θ = 0 degrees is taken as the forward direction). The truncation of the polynomial series is addressed by performing the Smirnov–Kolmogorov test.⁴⁸ Significance levels higher than 94% can be achieved using 4 to 12 Legendre moments.

3 Results and discussion

3.1 Excitation functions

Fig. 1 shows the excitation functions for both the CH + D and CD + H channels of the state-specific C(¹D) + HD (v = 0, j = 0) reaction on the ZMB-a surface, together with the previous theoretical results on the BHL and DMBE surfaces.^19,26 As can be seen, for both channels the present HB and GB ICSs decrease sharply from very large values with the increase of the collision energy and then keep nearly constant at high collision energies, indicating that the title reaction proceeds through the insertion mechanism on the ground surface. The similar trend can also be found in the previous theoretical excitation functions,^19,26 however, some differences should be noted. At the collision energy of 0.01 eV, both the ICSs on the ZMB-a and DMBE PESs are much larger than those obtained from the BHL PES. Considering that the deep well regions of the three PESs are broadly similar and the CI regions are not accessible at low energies, the much smaller ICSs on the BHL PES at low energies could be attributed to the deficiencies of the BHL PES in the long-range region. In particular, as mentioned in our previous work,³³ the vdW saddle on the ZMB-a PES enhances reactivity at low collision energies, which is absent on the BHL PES. In addition, the present CH + D ICSs are much larger than those on the DMBE PES at the collision energies of 0.02–0.3 eV, and the present CD + H ICSs are somewhat larger than those on the DMBE PES at the collision energies of 0.02–0.1 eV, which could be attributed to the different topographical features of the two PESs in the CI and vdW regions. The different PES topographical features between the ZMB-a and DMBE PESs are further discussed in the following section.

Fig. 1 Excitation functions for both the CH + D (upper panel) and the CD + H (lower panel) channels of the C(¹D) + HD (v = 0, j = 0) reaction obtained from various singlet ground surfaces. The black squares display the QCT-HB results on the ZMB-a surface, the blue triangles indicate the QCT-GB results on the ZMB-a surface, the red circles give the QCT-HB results on the DMBE surface,²⁶ and the pink dashed line presents the quantum wave packet (WP) results on the BHL surface.¹⁹

By comparing the present HB and GB ICSs in Fig. 1, we can find that for collision energies larger than 0.02 eV the difference is not evident. However, at the collision energy of 0.01 eV, the GB ICS is larger than the HB ICS in the CH + D channel, while in the CD + H channel the GB ICS is smaller than its HB counterpart. To investigate the cause of this phenomenon, we check the product vibrational state resolved ICSs at this collision energy. For the CH + D channel, the present HB ICS (v′ = 0), 30.3 ± 0.5 Ų, is smaller than the GB ICS (v′ = 0), 33.6 ± 1.5 Ų, and the difference between the HB ICS (v′ = 1), 0.6 ± 0.1 Ų and the GB ICS (v′ = 1), 0.0 Ų, is very small; for the CD + H channel, the present HB ICS (v′ = 0), 49.3 ± 0.5 Ų, is nearly the same with the GB ICS (v′ = 0), 48.6 ± 1.6 Ų, while the HB ICS (v′ = 1) is 7.6 ± 0.3 Ų, which is much larger than the GB ICS (v′ = 1), 2.1 ± 0.4 Ų. Obviously, the CD + H ICS (v′ = 1) is overestimated in the HB approach, which results from that many reactive trajectories with v′(CD) being in the range from 0.5 to 1.5 are treated as effective reactive trajectories in the HB approach but assigned very small weights in the GB approach. Thus the ICSs, especially those for the vibrational excitation product channel, are very sensitive to the ZPE leakage at low collision energies, where the GB approach should be used to correct the ZPE leakage.

In addition, we investigate the contribution of the excited b̃¹A′′ PES. As shown in Fig. 2, for both channels, the excitation functions on the b̃¹A′′ PES show similar trend with the collision energy to those on the ground PES; the contribution of the b̃¹A′′ PES to the total reactivity is about 25–45%. In particular, at the collision energy of 0.038 eV, for the CH + D channel the present GB ICS values are 16.0 ± 0.7 Ų on the ã¹A′ surface and 5.6 ± 0.5 Ų on the b̃¹A′′ surface, which can be compared with the WP values¹⁹ of 13.6 Ų on the ã¹A′ BHL surface and 7.0 Ų on the b̃¹A′′ BJHL surface by Defazio et al.; for the CD + H channel, the present GB ICS values are 25.8 ± 0.9 and 19.0 ± 0.8 Ų on the ã¹A′ and b̃¹A′′ surfaces respectively, and the corresponding WP values¹⁹ are 17.6 and 14.5 Ų.

Fig. 2 Excitation functions for both the CH + D (upper panel) and the CD + H (lower panel) channels of the C(¹D) + HD (v = 0, j = 0) reaction obtained from our ã¹A′ and b̃¹A′′ surfaces. The corresponding contribution of the b̃¹A′′ surface is also displayed.

3.2 Isotopic branching ratios

The calculated product CD/CH branching ratios obtained by ICS_CD+H/ICS_CH+D for the C(¹D) + HD (v = 0, j = 0) reaction as a function of the collision energy are presented in Fig. 3. As can be seen, at the collision energy of 3.7 kJ mol⁻¹, the present HB value of the product CD/CH branching ratio is 1.58 ± 0.03 and the GB value is 1.61 ± 0.10 on the ZMB-a surface, whereas the HB result from the DMBE surface is 1.91 (ref. 26) and that from the BHL surface is 1.78.²⁰ Although the CD/CH branching ratios on the ZMB-a surface are in very good agreement with the experimental value of 1.6 ± 0.1 reported by Sato et al.¹⁷ and an earlier experimental value¹⁸ of 1.69 ± 0.07 which is measured at collision energies less than 4.2 kJ mol⁻¹, when the b̃¹A′′ contribution is included, the overall HB product CD/CH branching ratio at the collision energy of 3.7 kJ mol⁻¹ is 1.93 ± 0.07 and the corresponding GB value is 2.07 ± 0.23, which become somewhat larger than experiment. However, Hickson et al.^49,50 have recently reported new experimental rate coefficients for the C(¹D) + H₂ and C(¹D) + D₂ reactions, and they concluded that several key secondary reaction rates required in the fitting procedure used by Sato et al.¹⁷ may be incorrectly estimated, indicating the potential uncertainty of their reported experimental branching ratio,¹⁷ thus more precise experiments on the C(¹D) + HD reaction are required.

Fig. 3 CD/CH product branching ratios for the C(¹D) + HD (v = 0, j = 0) reaction as a function of the collision energy on our ã¹A′ ZMB-a surface (upper panel) and both the ã¹A′ and b̃¹A′′ surfaces (lower panel). On the upper panel, the black squares display the HB results, the blue triangles give the GB results, the violet line gives the QCT-HB results on the DMBE surface,²⁶ and the blue dot line presents the QCT-HB results on the BHL surface.²⁰ On the lower panel, the pink squares display the HB results, and the navy triangles give the GB results. The experimental value of Sato et al.¹⁷ (in solid red circle) is also presented.

More interestingly, as seen from the upper panel of Fig. 3, at collision energies smaller than 0.1 eV, the present GB branching ratios on the ZMB-a surface do not change drastically, while it seems the corresponding HB results decrease slightly with the increase of the collision energy and the trend is similar to the previous QCT-HB results on the BHL and DMBE PESs,^20,26 which may result from that in the HB approach the reactivity of the CD (v′ = 1) channel is overestimated at low collision energies; for collision energies larger than 0.3 eV, both our HB and GB results tend to increase slightly with the rise of the collision energy, while the previous QCT results are nearly collision energy independent. Nevertheless, it shows both the present and previous QCT results on the ground surfaces for the CD + H channel are more populated than the CH + D channel in the whole energy range. This can be explained from a statistical point of view as follows. The title reaction proceeds via a deep well and the yielding of product CD(H) can be regarded as the decomposing of the D–C–H complex, the long lifetime of which leads to a statistical population of the internal states. Then more states of CD related are populated than those of CH related as a result of the difference in the reduced masses. Because of that each state of the D–C–H complex has equal probability of decomposing, the CD product will possess more open channels than the CH product at a given collision energy and is statistically favored. Similar phenomenon has also been found in other insertion reactions such as O(¹D) + HD.⁵¹ When the b̃¹A′′ contribution is included, both the HB and GB branching ratios display almost negative collision energy dependence at collision energies smaller than 0.3 eV, indicating that the dynamics picture on our b̃¹A′′ PES is more complex.

It should be noted that the present branching ratios on the ZMB-a surface are appreciably smaller than those obtained from the DMBE surface over the whole energy range, which can be explained by the distinct PES topographical features. In the following we give a qualitative analysis based on the comparison of the contour plots for the ZMB-a (Fig. 4) and DMBE (Fig. 10 in ref. 25) surfaces, where the dissociated process from intermediate to product can be regarded as the H (or D) atom leaving from the deep well to the product region in various directions (corresponding to various potential energy curves). We use red dot lines to distinguish two kinds of regions: for the region above the red dot lines, the energy is rising monotonously from the deep well to the product region without any barrier in each potential energy curve; whereas for the region under the red dot lines, there will be a barrier caused by CIs along each potential energy curve, and this barrier may inhibit the D + CH channel because it is more difficult to surmount the barrier for the heavier D-atom than for the lighter H-atom. It can be easily found that the regions under the red dot lines on the DMBE PES are much more widespread than those on ZMB-a, which may lead to an underestimation of the reaction probability of the D + CH channel and result in larger CD/CH branching ratios.

Fig. 4 Contour plot of the ZMB-a PES for the H atom moving around the CH product with the CH bond length fixed at its equilibrium bond length of R_CH = 2.116 bohr, where the CH product lies along the X-axis with the center of the bond located at the coordinate origin. Contours are in kcal mol⁻¹ and the energy of the C(¹D) + HD asymptote is taken as zero. The blue lines are contours equally spaced by 5.0 kcal mol⁻¹, while the black lines are contours equally spaced by 0.02 kcal mol⁻¹.

3.3 Thermal rate coefficients

The thermal rate coefficients of the CH + D and CD + H channels as well as the thermal CD/CH branching ratios (i.e., k_CD+H/k_CH+D) calculated with the QCT-HB approach are presented in Table 1, where the previous QCT-HB results on the DMBE surface and available experimental data are also listed for comparison. As seen, the rate coefficients for both channels are weakly dependent on temperature from 200 to 1500 K. Besides, the present ã¹A′ rate coefficients are larger than those computed on the DMBE surface over the whole temperature range, resulting from different topologies of the two PESs in the CI and vdW regions. The b̃¹A′′ surface contributes noticeably to the overall rate coefficient, which accounts for about 30–43%. At room temperature, the present thermal rate coefficient for the title reaction is (2.26 ± 0.14) × 10⁻¹⁰ cm³ s⁻¹, which is in good agreement with the experimental value¹⁷ of (1.7 ± 0.4) × 10⁻¹⁰ cm³ s⁻¹. In particular, our value of the CD + H channel is (1.50 ± 0.06) × 10⁻¹⁰ cm³ s⁻¹, which agrees well with the experimental value¹⁷ of (1.2 ± 0.5) × 10⁻¹⁰ cm³ s⁻¹; as for the CH + D channel, our value of (0.76 ± 0.04) × 10⁻¹⁰ cm³ s⁻¹ is also in very good agreement with the experimental value¹⁷ of (0.7 ± 0.3) × 10⁻¹⁰ cm³ s⁻¹.

Table 1 Thermal rate coefficients k and CD/CH branching ratios Γ_CD/CH for the C(¹D) + HD reaction

T (K)	State	10^10k (cm^3 s^−1)									ΓCD/CH
		CD + H			CH + D			Total
		This work	Others^a	Exp^b	This work	Others^a	Exp^b	This work	Others^a	Exp^b	This work	Others^a
a Ref. 26 QCT-HB/DMBE.b Ref. 17.
200	ã^1A′	0.87 ± 0.02	0.79 ± 0.02		0.54 ± 0.01	0.38 ± 0.02		1.41 ± 0.04	1.17		1.61 ± 0.05	2.08
	b̃^1A′′	0.65 ± 0.01			0.23 ± 0.01			0.88 ± 0.04			2.83 ± 0.13
	ã^1A′ + b̃^1A′′	1.52 ± 0.04			0.77 ± 0.04			2.29 ± 0.13			1.97 ± 0.11
298	ã^1A′	0.88 ± 0.02	0.73 ± 0.02	1.2 ± 0.5	0.53 ± 0.01	0.35 ± 0.01	0.7 ± 0.3	1.41 ± 0.03	1.08	1.7 ± 0.4	1.67 ± 0.04	2.08
	b̃^1A′′	0.62 ± 0.02			0.23 ± 0.01			0.85 ± 0.05			2.70 ± 0.15
	ã^1A′ + b̃^1A′′	1.50 ± 0.06			0.76 ± 0.04			2.26 ± 0.14			1.97 ± 0.12
500	ã^1A′	0.83 ± 0.02	0.66 ± 0.02		0.52 ± 0.02	0.36 ± 0.04		1.35 ± 0.06	1.02		1.61 ± 0.07	1.83
	b̃^1A′′	0.58 ± 0.02			0.26 ± 0.01			0.84 ± 0.04			2.23 ± 0.12
	ã^1A′ + b̃^1A′′	1.41 ± 0.06			0.78 ± 0.04			2.19 ± 0.15			1.81 ± 0.12
1000	ã^1A′	0.80 ± 0.03	0.67 ± 0.03		0.50 ± 0.02	0.35 ± 0.02		1.30 ± 0.07	1.02		1.59 ± 0.08	1.91
	b̃^1A′′	0.52 ± 0.02			0.30 ± 0.02			0.82 ± 0.06			1.73 ± 0.13
	ã^1A′ + b̃^1A′′	1.32 ± 0.07			0.80 ± 0.06			2.12 ± 0.20			1.65 ± 0.16
1500	ã^1A′	0.83 ± 0.03	0.76 ± 0.03		0.50 ± 0.02	0.41 ± 0.02		1.33 ± 0.08	1.17		1.68 ± 0.10	1.85
	b̃^1A′′	0.52 ± 0.02			0.34 ± 0.02			0.86 ± 0.06			1.53 ± 0.11
	ã^1A′ + b̃^1A′′	1.35 ± 0.07			0.84 ± 0.06			2.19 ± 0.19			1.61 ± 0.14

---

We can also see from Table 1 that the present thermal CD/CH branching ratio is nearly temperature-independent on the ZMB-a surface, which is consistent with the previous QCT results on the DMBE surface. It is interesting that the difference between the CD/CH branching ratios under thermal conditions and those at the initial state-specific (v = 0, j = 0) conditions is not significant on the ZMB-a surface, indicating that the CD/CH branching ratios are not sensitive to the initial reagent HD rotational excitation. This is understandable because the long-lived intermediate makes almost all memory of the initial conditions of the reagents be lost. However, the total thermal CD/CH branching ratio shows negative temperature dependence, resulting from that the thermal branching ratio on our b̃¹A′′ surface decreases with the increasing temperature, indicating nonstatistical dynamics behaviours on the b̃¹A′′ surface. Obviously, the addition of the b̃¹A′′ contribution can not only increase the absolute rate coefficient of the two product channels, but also influence the temperature dependence of the thermal CD/CH branching ratio, thus it is rather necessary to investigate reactions which happen on the b̃¹A′′ surface.

3.4 Product state and angular distributions

Fig. 5 shows the product vibrational distributions for the C(¹D) + HD (v = 0, j = 0) → CH(CD) (v′, j′) + D(H) reaction at different collision energies. As can be seen, for both channels, the ICS (v′ = 1)/ICS (v′ = 0) ratios are smaller than 0.7 over the whole collision energies, indicating that most products (CH or CD) are in the ground vibrational state; the inclusion of the b̃¹A′′ contribution further reduces the total ICS (v′ = 1)/ICS (v′ = 0) ratios evidently. For each channel, both the HB and GB ICS (v′ = 1)/ICS (v′ = 0) ratios tend to increase with the rise of the collision energy, and the GB ratios are much smaller than the HB ratios at low collision energies. At given collision energy, the ICS (v′ = 1)/ICS (v′ = 0) ratio of the CH + D channel is smaller than that of the CD + H channel, which consists with the consequence of decomposing the D–C–H complex. In particular, at low collision energies, the GB results display no threshold for the CD (v′ = 1) product in the CD + H channel while predict a reaction threshold of 0.1 eV for the CH (v′ = 1) product in the CH + D channel. Therefore at collision energies less than 0.1 eV, only CD (v′ = 1) product can be found, this interesting phenomenon would be observed experimentally.

Fig. 5 HB and GB results of product vibrational branching ratio defined as ICS (v′ = 1)/ICS (v′ = 0) as a function of collision energy for both the CH + D (upper panel) and CD + H (lower panel) channels of the C(¹D) + HD (v = 0, j = 0) reaction.

The calculated rotationally state-resolved ICSs for the C(¹D) + HD (v = 0, j = 0) → CH(CD) (v′ = 0, j′) + D(H) reaction at the collision energy of 0.038 eV are displayed in Fig. 6. For each channel, both the HB and GB rotational state distributions on the ZMB-a surface are highly inverted and peak at high rotational states, which supports the complex-forming mechanism. In particular, the GB rotational state distribution of the CH + D channel shows a sharp peak at j′ = 9 and extends to j′ = 13, while that of the CD + H channel peaks at j′ = 13 and reaches j′ = 18. For both channels, the inclusion of the b̃¹A′′ contribution makes the total rotational distributions a little cooler. The rotational state distribution of the CD + H channel is hotter than that of the CH + D channel, which consists with the experiment.¹⁷

Fig. 6 HB and GB results of the rotationally state-resolved integral cross sections calculated at the collision energy of 0.038 eV for both the CH (v′ = 0, j′) + D (upper panel) and CD (v′ = 0, j′) + H (lower panel) channels of the C(¹D) + HD (v = 0, j = 0) reaction.

The DCSs for both the CH + D and CD + H channels of the C(¹D) + HD (v = 0, j = 0) reaction at the collision energy of 0.038 eV are shown in Fig. 7. As expected, for each channel, both the HB and GB DCSs on the ZMB-a surface are nearly forward–backward symmetric with considerable polarization between forward/backward and sideways scattering, which is a clear indication of that the reaction proceeds via a long-lived complex. The present HB DCS of the CH + D channel on the ZMB-a surface is consistent with the previous QCT result on the BHL PES,²¹ which is also forward–backward symmetric with a similar polarization between forward/backward and sideways scattering. When the b̃¹A′′ contribution is included, the overall DCSs for both the CH + D and CD + H channels are also nearly forward–backward symmetric, and in contrast to this, a forward bias appears in the overall DCS when the contribution of the b̃¹A′′ BJHL PES is added.⁵² Unfortunately there has been no experimental data for the DCS of the title reaction up to now, and more experimental studies are desired to verify the theoretical predictions.

Fig. 7 HB and GB results of the differential cross sections calculated at the collision energy of 0.038 eV for both the CH + D (upper panel) and CD + H (lower panel) channels of the C(¹D) + HD (v = 0, j = 0) reaction.

4 Conclusions

In this work, we have reported the detailed QCT calculations for the C(¹D) + HD reaction on our recently constructed ã¹A′ and b̃¹A′′ ab initio PESs. The isotopic product CD/CH branching ratios, thermal rate coefficients, and other dynamical quantities have been investigated. The present ICSs on the ZMB-a surface are much larger than those on the previous singlet ground surfaces at low collision energies, which can be attributed to the improvement of the ZMB-a surface in the vdW and CI regions. The calculated overall CD/CH branching ratios show almost negative collision energy dependence at collision energies smaller than 0.3 eV; at the collision energy of 3.7 kJ mol⁻¹, our results are in reasonable agreement with the experimental value. The calculated thermal rate coefficients at room temperature are in very good agreement with the available experimental data. For each product channel, the DCS is dominated by scattering in both forward and backward directions with almost forward–backward symmetry; most products are in the vibrational ground state; the rotational state distribution for the CD + H channel is hotter than that for the CH + D channel, which is consistent with the available experimental measurements. The contribution of the b̃¹A′′ PES accounts for around 25–45% to the total ICSs and rate coefficients; the addition of the b̃¹A′′ contribution reduces the ICS (v′ = 1)/ICS (v′ = 0) ratio and leads to a somewhat colder rotational distribution, while the overall DCSs remain nearly forward–backward symmetric. We have also compared the various dynamical quantities obtained using the HB approach with their GB counterparts; although the GB DCSs are similar to the HB ones, the GB approach yields much smaller ICS (v′ = 1), colder product vibrational and rotational distributions for both channels than those obtained by the HB approach, indicating that the ZPE problem could be effectively mitigated by the GB approach. It is our hope that the present study would stimulate further experimental interests in the title reaction.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This work is supported by the Chinese Ministry of Science and Technology (No. 2013CB834601), the National Natural Science Foundation of China (No. 21473218 and 21303217), and the Beijing National Laboratory for Molecular Sciences and the CAS Key Research Project of Frontier Science. Many computations are carried out at National Supercomputing Centers in Tianjin (on TianHe-1(A)) and Shenzhen.

References

1. D. Skouteris, D. Manplopoulos, W. Bian, H.-J. Werner, L. Lai and K. Liu, Science, 1999, 286, 1713–1716.
2. W. Bian and H.-J. Werner, J. Chem. Phys., 2000, 112, 220–229.
3. T. González-Lezana, Int. Rev. Phys. Chem., 2007, 26, 29–91.
4. K. Liu, Int. Rev. Phys. Chem., 2001, 20, 189–217.
5. P. Casavecchia, Rep. Prog. Phys., 2000, 63, 355–414.
6. F. J. Aoiz, L. Bañares and V. J. Herrero, J. Phys. Chem. A, 2006, 110, 12546–12565.
7. D. C. Scott, J. de Juan, D. C. Robie, D. Schwartz-Lavi and H. Reisler, J. Phys. Chem., 1992, 96, 2509–2518.
8. K. Mikulecky and K.-H. Gericke, J. Chem. Phys., 1993, 98, 1244–1251.
9. M. R. Scholefield, S. Goyal, J.-H. Choi and H. Reisler, J. Phys. Chem., 1995, 99, 14605–14613.
10. Z. Zhang, H. Ma and W. Bian, J. Chem. Phys., 2011, 135, 154303.
11. A. Bergeat, L. Cartechini, N. Balucani, G. Capozza, L. F. Philips, P. G. Casavecchia, G. Volpi, L. Bonnet and J.-C. Rayez, Chem. Phys. Lett., 2000, 327, 197–202.
12. N. Balucani, G. Capozza, L. Cartechini, A. Bergeat, R. Bobbenkamp, P. Casavecchia, F. J. Aoiz, L. Bañares, P. Honvault, B. Bussery-Honvault and J.-M. Launay, Phys. Chem. Chem. Phys., 2004, 6, 4957–4967.
13. N. Balucani, P. Casavecchia, F. J. Aoiz, L. Bañares, J.-M. Launay, B. Bussery-Honvault and P. Honvault, Mol. Phys., 2010, 108, 373–380.
14. N. Balucani, G. Capozza, E. Segoloni, R. Bobbenkamp, P. Casaveechia, T. Gonzalez-Lezana, E. J. Rackham, L. Bañares and F. J. Aoiz, J. Chem. Phys., 2005, 122, 234309.
15. P. Defazio, C. Petrongolo, B. Bussery-Honvault and P. Honvault, J. Chem. Phys., 2009, 131, 114303.
16. P. Defazio, B. Bussery-Honvault, P. Honvault and C. Petrongolo, J. Chem. Phys., 2011, 135, 114308.
17. K. Sato, N. Ishida, T. Kurakata, A. Iwasaki and S. Tsuneyuki, Chem. Phys., 1998, 237, 195–204.
18. W. H. Fisher, T. Carrington, C. M. Sadowski and C. H. Dugan, Chem. Phys., 1985, 97, 433–448.
19. P. Defazio, P. Gamallo, M. González, S. Akpinar, B. Bussery-Honvault, P. Honvault and C. Petrongolo, J. Chem. Phys., 2010, 132, 104306.
20. L. Kang and M. Zhu, J. Mol. Struct., 2010, 945, 116–119.
21. L. Kang, Int. J. Quantum Chem., 2011, 111, 117–122.
22. S. Lin and H. Guo, J. Chem. Phys., 2004, 121, 1285–1292.
23. B. Bussery-Honvault, P. Honvault and J.-M. Launay, J. Chem. Phys., 2001, 115, 10701–10708.
24. L. Bañares, F. J. Aoiz, S. A. Vázquez, T.-S. Ho and H. Rabitz, Chem. Phys. Lett., 2003, 374, 243–251.
25. S. Joseph and A. J. C. Varandas, J. Phys. Chem. A, 2009, 113, 4175–4183.
26. S. Joseph, P. J. S. B. Caridade and A. J. C. Varandas, J. Phys. Chem. A, 2011, 115, 7882–7890.
27. Z. Sun, C. Zhang, S. Lin, Y. Zheng, Q. Meng and W. Bian, J. Chem. Phys., 2013, 139, 014306.
28. C. Zhang, M. Fu, Z. Shen, H. Ma and W. Bian, J. Chem. Phys., 2014, 140, 234301.
29. Z. Shen, J. Cao and W. Bian, J. Chem. Phys., 2015, 142, 164309.
30. Y. Wu, C. Zhang, J. Cao and W. Bian, J. Phys. Chem. A, 2014, 118, 4235–4242.
31. X. Liu, W. Bian, X. Zhao and X. Tao, J. Chem. Phys., 2006, 125, 074306.
32. B. Bussery-Honvault, J. Julien, P. Honvault and J.-M. Launay, Phys. Chem. Chem. Phys., 2005, 7, 1476–1481.
33. Z. Shen, H. Ma, C. Zhang, M. Fu, Y. Wu, W. Bian and J. Cao, Nat. Commun., 2017, 8, 14094.
34. S. J. Greaves, E. Wrede, N. J. Goldberg, J. Zhang, D. J. Miller and R. N. Zare, Nature, 2008, 454, 88–91.
35. J. Cao, Z. Zhang, C. Zhang, K. Liu, M. Wang and W. Bian, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 13180–13185.
36. M. Wang, X. Sun and W. Bian, J. Chem. Phys., 2008, 129, 084309.
37. A. J. C. Varandas, ChemPhysChem, 2005, 6, 453–465.
38. L. Bonnet and J.-C. Rayez, Chem. Phys. Lett., 1997, 277, 183–190.
39. L. Bañares, F. J. Aoiz, P. Honvault, B. Bussery-Honvault and J.-M. Launay, J. Chem. Phys., 2003, 118, 565–568.
40. L. Bonnet and J.-C. Rayez, Chem. Phys. Lett., 2004, 397, 106–109.
41. A. J. C. Varandas, Chem. Phys. Lett., 2007, 439, 386–392.
42. F. J. Aoiz, V. J. Herrero and V. S. Rábanos, J. Chem. Phys., 1992, 97, 7423–7436.
43. F. J. Aoiz, L. Bañares and V. J. Herrero, J. Chem. Soc., Faraday Trans., 1998, 94, 2483–2500.
44. T. D. Sewell and D. L. Thompson, Int. J. Mod. Phys. B, 1997, 11, 1067–1112.
45. G. H. Peslherbe, H. Wang and W. L. Hase, Adv. Chem. Phys., 1999, 105, 171–201.
46. W. L. Hase, R. J. Duchovic, X. Hu, A. Komornicki, K. F. Lim, D.-H. Lu, G. H. Peslherbe, K. N. Swamy, S. R. V. Linde and A. J. C. Varandas, VENUS96: A General Chemical Dynamics Computer Program; Bulletin Number 16 of the Quantum Chemistry Program Exchange, Indiana University, Bloomington, IN, 1996, p. 671.
47. X. Hu, W. L. Hase and T. Pirraglia, J. Comput. Chem., 1991, 12, 1014–1024.
48. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in Fortran 77, Cambridge University Press, Cambridge, U.K., 1992.
49. K. M. Hickson, J.-C. Loison, H. Guo and Y. V. Suleimanov, J. Phys. Chem. Lett., 2015, 6, 4194–4199.
50. K. M. Hickson and Y. V. Suleimanov, Phys. Chem. Chem. Phys., 2017, 19, 480–486.
51. M. Hankel, G. G. Balint-Kurti and S. K. Gray, J. Phys. Chem. A, 2001, 105, 2330–2339.
52. P. Honvault, B. Bussery-Honvault, J.-M. Launay, F. J. Aoiz and L. Bañares, J. Chem. Phys., 2006, 124, 154314.

---