Dynamical importance of van der Waals saddle and excited potential surface in C(1D)+D2 complex-forming reaction

Encouraged by recent advances in revealing significant effects of van der Waals wells on reaction dynamics, many people assume that van der Waals wells are inevitable in chemical reactions. Here we find that the weak long-range forces cause van der Waals saddles in the prototypical C(1D)+D2 complex-forming reaction that have very different dynamical effects from van der Waals wells at low collision energies. Accurate quantum dynamics calculations on our highly accurate ab initio potential energy surfaces with van der Waals saddles yield cross-sections in close agreement with crossed-beam experiments, whereas the same calculations on an earlier surface with van der Waals wells produce much smaller cross-sections at low energies. Further trajectory calculations reveal that the van der Waals saddle leads to a torsion then sideways insertion reaction mechanism, whereas the well suppresses reactivity. Quantum diffraction oscillations and sharp resonances are also predicted based on our ground- and excited-state potential energy surfaces.

T he van der Waals (vdW) interactions are important in many chemical and biological processes 1-3 such as the formation of tertiary structures of biopolymers, electron tunnelling in protein crystals and urea's preferential solvation of peptides, to name only a few. The possible role of weak vdW interactions on reaction dynamics is an issue of great interest 2,[4][5][6][7][8][9][10][11] . Although the importance of orientation or stereodynamical effects of the long-range anisotropic interactions was realized earlier [12][13][14][15][16] , and has recently become apparent in the emerging field of ultracold bimolecular reactions 17 , it is quite remarkable that the shallow vdW well 2 in the entrance channel of a classic activated chemical reaction (Cl þ H 2 ) can have an evident effect on the outcome of the reaction at low collision energies, which has been well recognized since the work of Skouteris et al. 2 Besides the stereodynamical effects on reactivity, the resonances caused by the vdW well are also detected; in particular, it is just found that the shallow reactant and product-valley vdW wells can support quasibound states that give rise to resonances in the F þ H 2 reaction 18 .
However, despite the exciting advances in revealing the effects of the vdW well in direct reactions dominated by the activation barrier 2,4-11 , the role of the vdW structure in the dynamics of a complex-forming reaction dominated by a deep well is still unclear. Currently, people tend to believe that the vdW wells also exist in other types of chemical reactions. Here we show that the weak vdW forces in the entrance channel of the C( 1 D) þ D 2 reaction (Fig. 1a), which is a prototypical complex-forming reaction [19][20][21][22] , give rise to entrance channel vdW saddle (see below for a definition and Supplementary Table 1), which has very different dynamical effects from vdW well at low collision energies. This is achieved with extensive dynamical calculations on the highly accurate ab initio potential energy surfaces (PESs) with vdW saddles constructed by us, and our findings are important for acquiring a deeper understanding of complex-forming reactions. For many basic bimolecular reactions such as H þ H 2 , Cl þ CH 4 and C þ D 2 , the vdW structure results from the dispersion 1 and quadrupole-quadrupole interactions between the reactants, which are rather weak, and an accurate description on such detailed structures remains a challenge nowadays in the construction of PESs, as is the construction of an excited-state PES.
The C( 1 D) þ D 2 reaction is much more difficult to characterize quantum mechanically than the activated reaction type, and has been the subject of many experimental and theoretical works. Experimentally, several kinetic and dynamical studies on the C( 1 D) þ D 2 reaction have been reported [23][24][25] , including a few crossed-beam experiments 24,25 . In particular, Liu 24 has reported the excitation function, and Balucani et al. 25 have obtained the product angular distributions and vibrational branching ratio. Theoretically, Bussery-Honvault et al. 26 have developed a global ab initio PES (called Bussery-Honvault-Honvault-Launay (BHL)) for theã 1 A 0 state; later, the ab initio data are refitted with the reproducing kernel Hilbert space (RKHS) method to remove some spurious features of the BHL PES, and the modified PES is denoted as the RKHS PES 27 . Although there are other PESs and subsequent dynamical calculations 28,29 , the RKHS PES is the most widely used in dynamical studies 27,[30][31][32][33][34] . Most recently, a highly accurate global ab initio PES (or Zhang-Ma-Bian-a (ZMB-a)) has been constructed by our group 35 , which is unique in the accurate description of the regions of vdW interactions and around conical intersections (CIs), and further dynamical calculations 36 performed on this PES for the C( 1 D) þ H 2 reaction confirm its accuracy. On the other hand, theb 1 A 00 excited-state PES may play an important role, and the only reported globalb 1 A 00 PES is the Bussery-Honvault-Julien-Honvault-Launay (BJHL) surface constructed by Bussery- Honvault et al. in 2005 (ref. 37), but the inclusion of the contribution of the BJHL surface worsens the agreement of the theoretical results with the crossed-beam experiments 38 , indicating the deficiency of the BJHL surface.
The main purpose of this work is to study the effects of weak vdW interactions and excited-state PES on a complex-forming reaction. The C( 1 D) þ D 2 reaction is chosen, since the crossedbeam scattering data by Liu are available for comparison. To include the contribution of theb 1 A 00 excited state, we have constructed an ab initiob 1 A 00 PES. To make it fully convincing, we have performed extensive quantum dynamics (QD) and quasiclassical trajectory (QCT) calculations on the singletground-(both ZMB-a and RKHS) and first-excited-state PESs for the C( 1 D) þ D 2 reaction.

Results
Excited PES and vdW saddle. The construction of accurate global ab initio ground-state PESs has played a key role in deepening our understanding of many polyatomic reactions; however, extension to the electronically excited state is rather difficult, and so far much less experience is available. Here a highly accurate global analytical PES for theb 1 A 00 excited state of the title system is constructed by fitting 46,600 symmetry unique ab initio energy points (see Methods and Supplementary Methods). The main advantage of ourb 1 A 00 PES is the accurate  description of the regions of the vdW interactions and those along the perpendicular insertion and around CIs, and further details are given in Supplementary Discussion. Stationary-point properties on ourb 1 A 00 surface are presented in Supplementary  Table 2. As can be seen, ourb 1 A 00 surface has a shallow well at a highly bent C 2v geometry with a depth of 14.4 kcal mol À 1 , which is absent on the BJHL surface 37 , and may play an important role in reaction dynamics. Meanwhile, saddle point instead of minimum is located in the vdW interaction region on ourb 1 A 00 surface, supported by our geometry optimization and frequency analysis with very accurate ab initio procedures (Supplementary Table 1). In other words, the presentb 1 A 00 PES constructed by us does not have any vdW well; instead, it has vdW saddles in the entrance and exit channels. This is remarkable, since vdW saddle could be regarded as the counterpart of vdW well, but is a concept unknown before for a chemical reaction. In addition, our ab initio calculations (Supplementary Table 1) and further analysis indicate that our ZMB-a surface also has similar vdW saddles. As shown in Fig. 1b, the local potential surface around the vdW stationary point as a function of two variables (C-D-D bending angle and C---DD distance) resembles a saddle very much, and that is why we call it vdW saddle.
Here we further induce a general definition of vdW saddle in a chemical reaction as: in case that the vdW interactions cause a saddle point instead of minimum (vdW well is absent), the PES region surrounding the local saddle point is referred to as vdW saddle. The key property of vdW saddle is its depth, which is here defined as the energy difference between the vdW saddle point and the corresponding asymptote. Then, the depth of the entrance-channel vdW saddle on ourb 1 A 00 PES is B0.22 kcal mol À 1 , somewhat shallower than that on ourã 1 A 0 PES (ZMB-a). It is interesting to reveal the dynamical effects of the vdW saddles, and the comparison of dynamical results on the two available PESs (ZMB-a with vdW saddle and RKHS with vdW well) makes it possible, as we will see below.
Excitation function and reaction mechanism. Figure 2 shows various integral cross-sections (ICSs) calculated by us together with the experimental measurements using the crossed-beam technique by Liu 24 . We see that the ICSs yielded from our QD calculations on ourã 1 A 0 andb 1 A 00 PESs are in very good agreement with experiment. As far as we know, this is the first effort to compare the theoretical results with the crossed-beam relative excitation function by Liu 24 and the agreement achieved on our surfaces is satisfactory in view of the degree of uncertainty in the measurements. The crossed-beam technique can determine the relative excitation function, but in principle it is unable to determine absolute ICSs 2 . We can also notice from Fig. 2a that theb 1 A 00 excited PES has an evident contribution to the absolute ICSs, which accounts for 430%. However, as shown in Fig. 2b, the shape of the decreasing ICS as the collision energy rises is quite similar to that obtained on the single ZMB-a surface. On the other hand, it can be easily noticed from Fig. 2b that, although the ICSs calculated on the ZMB-a and RKHS PESs display a similar trend and have close values at higher collision energies, the ICSs from ZMB-a have much larger values at low collision energies (E c o0.05 eV) and are much more consistent with the crossedbeam excitation function. This phenomenon is attributed to the improvement of our PES (ZMB-a) in the vdW regions in the following analysis. Furthermore, the ICSs calculated with the QCT method on ZMB-a are shown to be in very good agreement with the QD results on the same PES, which justifies the following analysis with the QCT method based on the different topologies of the ZMB-a and RKHS PESs.
The deep well regions of the ZMB-a and RKHS PESs are broadly similar, and both PESs have barriers at similar linear C-D-D geometries and with similar heights. In addition, both PESs support entrance-channel vdW complexes with similar linear equilibrium geometries and energies. However, the two PESs have very different topologies around the entrance-channel stationary points; as mentioned above, the RKHS PES has a vdW well, whereas the ZMB-a PES exhibits a vdW saddle, which is rather flat around the stationary point or local saddle point. It is interesting to mention that the depth of the vdW saddle on ZMB-a is B0.3 kcal mol À 1 (Fig. 1b), which is very close to the depth of the vdW well on RKHS. It should be noted that the ab initio calculations for ZMB-a are much larger than those for RKHS, and we think it is important to treat the five states simultaneously (see Methods) due to the mixing of five electronic states in the entrance channel, but we notice that the single-state multireference configuration interaction calculations are performed for the RKHS PES 26,27 .
To account for the much larger ICSs on ZMB-a than on RKHS at low collision energies, the QCT method is used. The QCT calculations could not only yield dynamical quantities in good accord with accurate QD results for the complex-forming reaction 39 , but also provide an intuitive picture of the reaction mechanism. Reactive classical trajectories leading to the D and CD products on ZMB-a are compared with trajectories from the same initial conditions on RKHS in Fig. 3. We find that most of the trajectories that react to produce CD þ D on ZMB-a are   Fig. 3 for two particular classical trajectories is thus considered to be quite general. On the other hand, detailed analyses reveal that, when the C atom approaches a non-rotating D 2 reagent, the reorientation effects induced by the vdW interactions of ZMB-a lead to a D-D torsion towards the linear C-DD vdW saddle point, and later, the C atom inserts into the D-D bond from a sideways direction. This is illustrated in Supplementary Fig. 2 and called as the D-D torsion then C sideways insertion reaction mechanism here. This mechanism is confirmed by direct examination of trajectories, and the analysis of the relative population of reactive trajectories indicates that it contributes most to the reactivity (Supplementary Fig. 3; Supplementary Discussion) at low collision energies. In addition, Fig. 3 shows that the left part of the vdW well on RKHS rises more steeply as R decreases than that of the vdW saddle on ZMB-a, which exhibits a platform structure around R ¼ 5.0 bohr; we consider this delicate structure on ZMB-a to be realistic since the fitting reproduces our highly accurate ab initio calculations very well in this region. Furthermore, we can also see from Fig. 3 that the long-range vdW interactions of ZMB-a extend to longer distances than those of RKHS, thus allowing reactive trajectories with larger impact parameter, which also contributes to the larger ICSs. At higher collision energies where the effects of vdW interactions are evidently weakened, the other small differences between the two surfaces make the ICSs on RKHS even become slightly larger than those on ZMB-a, which further highlights the suppression role of the vdW well at low collision energies.
There are five electronic states that correlate to the reagents C( 1 D) þ D 2 , but the lowest two (ã 1 A 0 andb 1 A 00 ) are the most relevant. As for the comparison with the excitation function measured in experiment, the inclusion of Renner-Teller effects is not expected to improve the agreement, since similar agreement is obtained (as shown in Fig. 2a,b) when the contribution of thẽ b 1 A 00 state is excluded. The three higher excited states can only participate in the reaction via the non-adiabatic pathways through the CIs, and their contributions are considered to be small in the collision energy range of this work. In particular, the minimum energy crossing point on the CI seam between theã 1 A 0 andc 1 A 0 states is located at 12.4 kcal mol À 1 as shown in Fig. 1a, below which thec 1 À0:08 cm A 0 andẽ 1 A 0 states are not accessible. Furthermore, although the rotational temperature of the experimental D 2 beam 24 is 100 K and only the D 2 (j ¼ 0) is considered here, our QCT results show that the difference between the ICS(j ¼ 0) and ICS(j ¼ 1, 2) is very small, thus the ortho-para nuclear spin effects could be neglected. Consequently, we think that, without further high-precision experimental work, it is unlikely to clarify the remaining small differences between the QD calculations on our PESs and experiment at low energies in Fig. 2.
Clearly the larger ICSs on the ZMB-a PES at low collision energies (Fig. 2b) will lead to larger rate coefficients. Supplementary Fig. 4 shows the initial state-specified rate coefficients as a function of temperature computed by us with the QD method on ourã 1 A 0 (ZMB-a) andb 1 A 00 PESs. We see that, theb 1 A 00 state makes a contribution of B38% to the overall reactivity, and our calculations suggest a larger rate coefficient than the experimental value of Sato et al. 23 at room temperature. Our result is indirectly supported by a most recent experiment of Hickson et al 34 . In this experiment, they 34 obtained a room temperature rate coefficient for the C( 1 D) þ H 2 reaction, which is larger than the upper limit reported by Sato et al. 23 even when the experimental uncertainty is considered. Further experimental work is required to clarify the divergence 23,34 in the measurements of rate coefficient.
Product angular and state distributions. We also perform largescale QD calculations on ourã 1 A 0 andb 1 A 00 PESs with vdW saddles to predict more detailed scattering properties encoded in differential cross-sections (DCSs). In particular, dense resonance structures and quantum diffraction oscillations are revealed, which manifest themselves as two kinds of different oscillation patterns in calculated DCSs (Fig. 4). We present the energy dependence of DCS in Fig. 4a, and we see that the DCSs near the two extreme angles (0°and 180°) display similar highly oscillatory structures, indicating densely distributed sharp resonances, which result from the numerous long-lived quasibound states supported by the deep well. We can also notice that the amplitude of the oscillations is reduced quickly when the scattering angle deviates from 0°or 180°. Meanwhile, we show the angular dependence of the product rotational state-resolved DCS at a representative collision energy in Fig. 4b, which exhibits an interesting oscillation pattern with peaks. This phenomenon may be caused by quantum interferences between partial waves, which appears to be similar to the diffraction oscillations just observed experimentally in the H þ D 2 reaction by Zare and coworkers using the photoloc technique 40,41 . However, they are of different nature, and as demonstrated by Zare and coworkers 41 , the latter is caused by well separated sets of total angular momentum, leading to scattering at the same angles.
In addition, we present a three-dimensional CD product flux surface plot at the collision energy of 0.038 eV in Fig. 5, in which the concentric circles correspond to different ro-vibrational states (n 0 , j 0 ) of the CD product. By integrating over the angular distribution shown in Fig. 5, we can obtain the rotational state distribution of the CD product, which is found to be in very good agreement with the available experimental measurements 23 ( Supplementary Fig. 5). We can also see that the product angular distribution is dominated by scattering in both the forward and backward directions with roughly forward-backward symmetry; this trend persists over the whole energy range, according with the complex-forming mechanism. In particular, the calculated results on ourã 1 A 0 andb 1 A 00 PESs at higher collision energies are consistent with the observations of the crossed-beam experiments 24,25 , in which the product angular distributions are found to be forward-backward symmetric. It should be noted that the calculated DCSs on theã 1 A 0 PES, either RKHS or ZMB-a, are nearly forward-backward symmetric 25,27,36 ; however previous calculations show that the addition of the contribution of theb 1 A 00 BJHL PES 37 leads to a forward bias in DCS 38,42 , in contradiction to experiment, indicating that the presentb 1 A 00 PES constructed by us is more accurate than the BJHL PES. PES is presented in Fig. 6, which is about the product CD vibrational branching ratio defined as ICS(n 0 ¼ 1)/ICS(n 0 ¼ 0). We see that the present QD result calculated on ourã 1 A 0 andb 1 A 00 PESs is in very good agreement with the available experimental value obtained by Balucani et al. 25 using the crossed-beam technique; meanwhile, the addition of the contribution of our excitedb 1 A 00 PES evidently reduces the total vibrational branching ratio and improves the agreement with experiment. In contrast, the present QD result calculated on the RKHS PES shows a much larger value, and the addition of the contribution of the BJHL (b 1 A 00 ) PES would evidently further increase the vibrational branching ratio according to the previous calculations 38,42 . A schematic lower limit of the QD branching ratio including contributions from both the RKHS (ã 1 A 0 ) and BJHL (b 1 A 00 ) PESs is also indicated in Fig. 6, and the contribution from the BJHL PES is indirectly estimated from the previous QD calculations for the C( 1 D) þ H 2 reaction 38 . Actually, the previous calculations on the BJHL excited PES yielded significantly larger vibrational branching ratio 38,42 than that on the singlet-groundstate PES (BHL or RKHS), and it was concluded by Casavecchia and coworkers 38 that the participation of the BJHL surface leads to an excessive contribution from the product CH(n 0 ¼ 1) channel, well above the possible experimental uncertainty. As shown in Supplementary Fig. 6, the topological feature of our b 1 A 00 PES is very different from that of the BJHL PES (Supplementary Discussion), and this should be the main cause of the large divergence. Our classical trajectory analysis further indicates that the nonstatistical and complex reaction mechanisms related to the shallow well structure on ourb 1 A 00 PES are responsible for the highly cold CD vibrational distribution. It should be mentioned that the D 2 beam used in the experiment 25 has a relative rotational population, and slightly larger theoretical values may be obtained than those shown in Fig. 6 when the initial D 2 rotational excitation is considered. We have calculated the vibrational branching ratio using the QCT method by averaging the contributions from different rotational states according to the experimental condition, and the obtained value is somewhat larger, but still in very good agreement with experiment. We think that the present theoretical values on our a 1 A 0 andb 1 A 00 PESs are within the experimental error limit according to the relevant statements of Casavecchia and coworkers 25,38 .

Discussion
Accurate QD and QCT calculations studying vdW interactions and excited-state PES contributions in the C( 1 D) þ D 2 reaction have been reported, using theã 1 A 0 andb 1 A 00 ab initio PESs constructed by us, the latter of which is reported in this work. Our PESs have entrance and exit-channel vdW saddles, and are unique in the accurate description of the vdW regions. The QD calculations using ourã 1 A 0 andb 1 A 00 PESs yield dynamical quantities in very good agreement with experiment, including the excitation function, product vibrational branching ratio and rotational state distribution; the obtained product angular distributions are consistent with those found in the crossedbeam experiments, which are forward-backward symmetric. These were not achieved on theã 1 A 0 andb 1 A 00 PESs reported by others. Furthermore, interesting diffraction oscillations and numerous sharp resonances are predicted. The importance of both vdW saddle and excitedb 1 A 00 PES in reaction dynamics is underscored. Our calculations indicate that theb 1 A 00 excited-state PES plays a significant role in the C( 1 D) þ D 2 reaction dynamics. The present QD calculations on our ZMB-a PES yield reasonable cross-sections supported by the recent crossed-beam experiment, whereas the same calculations on the earlier RKHS PES with the vdW well produce much smaller cross-sections at low collision energies, at which detailed analyses with the QCT method reveal that the anisotropic interactions around the vdW saddle point on ZMB-a lead to a D-D torsion and then C sideways insertion reaction mechanism, whereas the vdW well on RKHS inhibits the reactivity.
The concept of 'vdW saddle' proposed in this work may have general importance, and our conjecture is, the vdW saddle originates from the nature of the polyatomic complex-forming reactions that involve deep wells. The shallow entrance-channel vdW well has played an important role in acquiring a deeper understanding of the direct reactions that are dominated by activation barriers; it is our hope that the vdW saddle will prove to be a useful concept in the study of complex-forming reactions that are dominated by deep wells.

Methods
Ab initio calculations and PES. All ab initio calculations were performed with the MOLPRO package 43 . The internally contracted multireference configuration interaction 44 , and five-state-averaged complete active space self-consistent field methods 45 were used. Five reference states were used for generating the internally contracted pairs. The active space consists of six electrons distributed among seven orbitals. The basis set used is the Dunning's correlation consistent quadruple-zeta basis augmented by diffuse functions (aug-cc-pVQZ) 46 . A highly accurate global PES for theb 1 A 00 state of the title system was constructed based upon 6,639 symmetry unique energy points obtained from our ab initio calculations. Accurate analytical fits were generated using many-body expansions with the permutationally invariant polynomials, and the scheme of PES construction is similar to that described in our recent work 35 , where a highly accurate global ground-state (ã 1 A 0 ) PES called ZMB-a was constructed. The root-mean-square error of the global fit is 0.42 kcal mol À 1 for energy points below 15 kcal mol À 1 relative to the C( 1 D) þ D 2 asymptote, though the fitting errors in dynamically important regions are much lower. To illustrate the global feature of ourb 1 A 00 surface, several typical contour plots as functions of Jacobi coordinates R and r at different chosen angles are shown in Supplementary Fig. 6. Further details about ourb 1 A 00 surface and discussion are given in Supplementary Methods and Supplementary Discussion.
QD and QCT calculations. The QD calculations including all Coriolis couplings were performed using the real-wave packet approach as implemented in the DIFFREALWAVE code 47 . Numerical details are given in Supplementary Table 3, and numerous test calculations were performed to check the convergence. All partial wave contributions up to J ¼ 50 were included, yielding converged ICSs and DCSs for collision energies up to B0.3 eV. About 120,000 iteration steps were implemented for each partial wave calculation, which were reduced at high J values. The QCT calculations were performed using a modified version of the VENUS96 code 48,49 customized to incorporate the PESs. Batches of 15,000-30,000 trajectories were calculated for each set of initial conditions with an integration time step of 0.02 fs, which guarantees an energy conservation of better than 1 in 10 6 . The trajectories were initiated at the C( 1 D) þ D 2 asymptote with an atom-diatom separation of 15.0 bohr. To fix the zero-point energy problem, the Gaussian binning 50 approach was used with a full width at half maximum of 0.1 in the Gaussian weighting function. Further details about the methods and calculations are given in Supplementary Methods and Supplementary References. Data availability. The data that support the findings of this study are available within the article (and its Supplementary Information files) and from the corresponding author on request.