Wei Yan*a,
Rui Shan Tan
a and
Shi Ying Lin
b
aSchool of Science, Shandong Jianzhu University, Jinan 250101, China. E-mail: yanwei19@sdjzu.edu.cn
bSchool of Physics, Shandong University, Jinan 250100, China
First published on 22nd May 2023
A global potential energy surface (PES) for the electronic ground state of the Na + HF reactive system is constructed by three-dimensional cubic spline interpolation of 37000 ab initio points obtained using the multireference configuration interaction method including the Davidson's correction (MRCI + Q) with auc-cc-pV5Z basis set. The endoergicity, well depth and properties of the separated diatomic molecules are in good agreement with experimental estimations. Quantum dynamics calculations have been performed and compared with those of the previous MRCI PES as well as experimental values. The better agreement between theory and experiment indicates the accuracy of the new PES.
Many experiments had been done to gain insights into the characteristics of Na + HF reaction.2–9 They found that the ground state reaction, namely Na(3S) + HF, was endoergic by about 0.87–0.93 eV5,7 and no reaction took place with low vibrational state of HF (v < 2).2,4–6 However, the reactivity was highly enhanced with the vibrational excitations of HF (v ≥ 2) and the rate constants increased linearly with v, which is a typical characteristic of the late-barrier reaction.2 Meanwhile, the influences of HF(j) initial rotational excitations had also been measured and the reaction rates showed initially decreasing then increasing dependence on j.3 In addition, from the rainbow structures of the nonreactive differential cross sections, the existence of van der Waals wells were predicted to locate in the entrance channel of both ground and first excited potential energy surfaces (PESs) and the estimations for the well depths were also given.5,6,8
As for theoretical studies, the PES is crucial in reaction dynamic calculations. For the Na + HF reaction system, the PESs of the ground as well as excited electronic states were constructed using semiempirical10 or ab initio8,9,11–17 methods. The ground state PESs and relevant dynamics studies had successfully reproduced some of the experimental results, such as the existence of van der Waals well and late barrier, the large endoergicity and the effect of initial vibrational excitation.13–15,17,23–29 The research on excited electronic states had been carried out by Sevin12 and Truhlar et al.15–17 using different methods. They confirmed the involvement of harpooning process during the reactions whether the reaction started in the ground or the excited electronic state. They also confirmed the experimental prediction of the existence of a deep well in the first excited state. Therefore, a metastable complex could be formed and dissociated to ground state via harpooning process.
We had also done some dynamics studies both in the adiabatic and nonadiabatic representations31,32 based on the 2 × 2 diabatic PES17 (denoted as TTCPP-B). In our previous work, for the Na(3P) + HF reaction, the agreement between the nonadiabatic results and the experiment implied that the interstate coupling term of the employed PES was quite accurate. For the Na(3S) + HF reaction, the probabilities obtained in the two representations were almost indistinguishable, which indicated the negligibly small nonadiabatic effect. So, if we focused on the dynamic characteristics of the ground state, the adiabatic PES was sufficient. The theoretical results were in good agreement with experimental ones as a whole. However, theory–experiment disparity still existed. The reaction rate for the v = 2 state was much smaller than the experimental data and the effects of the initial rotational excitations were only in qualitative agreement with the experiment, which implied that the PES probably overestimated the true barrier height. Further refinement of the ground state PES was needed to close the theoretical and experimental gap.
Following our previous work, in this paper, we construct a new PES for the ground electronic state of NaFH system. The MRCI + Q level of theory is used with aug-cc-pV5Z basis set. The global PES was obtained by fitting a larger number of ab initio points (37000 points) which cover not only the interaction region but also the asymptotic regions of Na + HF and NaF + H channels. Quantum wave packet calculations are performed based on the new PES and compared with the experimental data to examine its accuracy. The new PES has a lower barrier than TTCPP-B PES and the dynamic results show better agreement with the experiment.
This letter is organized as follows. In Section 2, the relevant theoretical and numerical methods are briefly described. In Section 3, the new developed PES and corresponding dynamic results are displayed and discussed. Final conclusions are made in the last section.
The reference configuration state functions for the MRCI calculations were taken from state-averaged complete active space self-consistent field (CASSCF)38,39 calculations for equally weighted 12A′ and 22A′ states. In the CASSCF calculations, the active space consists of seven electrons distributed in eight valence orbitals, including the 3s and 3p orbitals of Na, the 2p orbitals of F and the 1s orbital of H. The inner orbitals were constrained to be doubly occupied but fully optimized. The total number of contracted configurations in the MRCI calculations was about 2.0 × 106 corresponding to about 6.6 × 107 uncontracted configurations.
A nonuniform and truncated direct product grid was used in the internal coordinates (RHF, RNaF, θNaFH), where RAB stands for the distance between A and B atoms and θNaFH stands for the bond angle enclosed by Na–F and H–F bonds. We chose 50 grids in the RHF coordinate ranging from 1.0 to 20.0 Bohr, 37 grids in the RNaF coordinate ranging from 2.2 to 25.0 Bohr and 20 grids in the θNaFH coordinate within 0°–180°, which gives rise to a total of 37000 grids. The ab initio energies computed on these grids were then used to generate the PES by fast three-dimensional cubic spline interpolation method.40
|ψ1〉 = DĤscaled|ψ0〉, |
|ψk+1〉 = D(2Ĥscaled|ψk〉 − D|ψk−1〉), k ≥ 1 |
The Hamiltonian and wave packet are expressed in the reactant Jacobi coordinates and discretized in terms of a mixed DVR/FBR (discrete variable representation/finite basis representation) scheme.54 The equidistant Fourier grid representation is utilized for the two radial degrees of freedom and the actions of the radial kinetic energy operators onto the wave packet are evaluated using the fast sine Fourier transform. For the angular degrees of freedom, the associated Legendre function and the Wigner rotation matrix are used. In this FBR, the rotational kinetic energy operators can be expressed as a diagonal or a tridiagonal matrix, so actions of them onto the wave packet can be evaluated in a straightforward manner. In order to calculate the action of the PES operator, the wave packet is transformed to a grid representation through a pseudospectral transformation55,56 so that the operator is diagonal in this case.
Feature | Property | This work | TTCPP-B17 | Exp57 |
---|---|---|---|---|
a The separated atom and diatomic molecule were set to 60 a.u. apart when performing the geometry optimization. The calculated results did not change with the further increase of this distance, indicating that there was no interaction between them, as was the case at infinity. | ||||
Reactant (Na + HF): | ||||
RHF | 1.73 | 1.73 | 1.73 | |
RNaF | 60.00a | ∞ | ∞ | |
D0 (HF) | 6.00 | 5.66 | 5.90 | |
![]() |
||||
van der Waals well (Na⋯FH): | ||||
RHF | 1.75 | 1.74 | — | |
RNaF | 4.60 | 4.68 | — | |
θNaFH | 116.59° | 118.0° | — | |
![]() |
||||
Transition state: | ||||
RHF | 2.54 | 3.49 | — | |
RNaF | 3.85 | 3.66 | — | |
θNaFH | 72.35° | 90.8° | — | |
![]() |
||||
Product (NaF + H): | ||||
RHF | 60.00a | ∞ | ∞ | |
RNaF | 3.69 | 3.64 | 3.64 | |
D0 (NaF) | 4.98 | 4.58 | 5.0 ± 0.1 |
As can be seen, the equilibrium geometries as well as the dissociation energies of HF and NaF molecules are in excellent agreement with the experimental data. The endoergicity, defined as ΔE = D0(HF) − D0(NaF), which is a key factor of the reaction, is 1.02 eV in this PES. Some uncertainty exists in the experimental value of D0 for NaF molecule. So, the experimental suggestions of endoergicity are about 0.87–0.93 eV.5,7 The location of the van der Waals well in the new PES is similar with TTCPP-B PES,17 but the well depth of the former is a little deeper than of the latter (−0.082 eV vs. −0.074 eV, with respect to the energy minimum of the Na(3S) + HF asymptote). The geometries of the transition state in the two PESs are quite different. The saddle point in the new PES has shorter F–H distance and less bent structure, although the Na–F distance is similar. Whereas, our value of θNaFH is in good agreement with that obtained by Laganà et al. (75°),14 and similar with Shapiro and Zeiri's result (60°).10 As we mentioned in our previous work,31 the noticeable difference between theory and experiment shown on the rates for Na + HF (v = 2) may due to the overestimated barrier height of the PES. The barrier height resulting from our ab initio calculation is 1.002 eV relative to the Na(3S) + HF asymptote which is lower than that on the TTCPP-B (1.270 eV)17 or semiempirical (1.266 eV) PES10 but close to Laganà's result (0.811 eV) in which some artificial corrections to the original ab initio energies are made.14 The lower barrier may have a significant effect on the reactivity and we will discuss this in detail in the next section.
To further characterize the reactive properties of the reaction, the schematic one-dimensional minimum energy path is displayed in Fig. 1. The eigen-energies of several vibrational states of the reactant and the ZPEs of the stationary points are also given to help illustrate the results of dynamic calculations. The detailed rovibrational eigen-energies of the HF/NaF (v = 0–5, j = 0–7) molecules which are obtained by solving the one-dimensional Schrödinger equation are listed in the ESI.† As the figure shows, the v ≤ 1 states of HF energetically lie below while the higher vibrational states above the saddle point. It is worth noticing that the v = 2 state energetically lies above the saddle point, but below the ground state of product by about 0.041 eV. Even though the barrier height is lower than TTCPP-B PES, the reaction probability of v = 2 state still has threshold energy which will be illustrated in detail below.
The fitted PES is exhibited by three-dimensional surfaces as a function of the RNaF and RHF bond coordinates in Fig. 2 and the corresponding two-dimensional contour plots are displayed in the ESI.† The bond angles (θNaFH) are equal to 60°, 72°, 117° and 180°, respectively. Among them, the surfaces for 72° and 117° outline the topography near the transition state and van der Waals well, respectively. As the figure shows, the PES is smooth all through the configuration space from the asymptotic to interaction regions. The well is visible in the surface for 117° in the entrance channel. The PES flattens so the saddle point is very inconspicuous which parallels the observation made by Sevin et al. at MRCI level12 and by Truhlar et al. at MBPT(2) level.15
Parameters | Na(3S) + HF (v, j = 0) |
---|---|
Grids/basis | NR = 445, R ∈ [0.5, 22] |
Nr = 115, r ∈ [0.5, 12] | |
Nθ = 45, θ ∈ [0°, 180°] | |
j = 0, 1, 2, … , 44 | |
![]() |
|
Damping; exp(−dx(x − xd)2), x = R or r | Rd = 18.0, dR = 0.005; rd = 8.0, dr = 0.005 |
![]() |
|
Initial wave packet; exp(−(R − R0)2/2δ2) × cos(k0R) | E0 = ℏ2k02/2μR = 0.2 eV; δ = 0.15, R0 = 16.0 |
![]() |
|
Dividing surface | rf = 7.9 |
![]() |
|
Propagation steps | 40![]() |
![]() |
|
Jmax | 130 (v = 2); 170 (v = 3); 180 (v = 4); 220 (v = 5) |
The Coriolis coupling effect is fully considered in the J > 0 partial waves. In our calculation, the body-fixed z axis is chosen to coincide with the vector which joins the Na atom and the center of mass of HF. Because the H atom is relatively light, we can assume that the z axis nearly coincides with the principal axis of inertia of the system. The Coriolis term should be rather small in this case. Therefore, we truncate the maximum helicity quantum number (Ω) to save computational costs. We set Ωmax = min(3, J) for J < 100 and Ωmax = 5 for J > 100 and it can provide indistinguishable results from the ones with Ωmax = J. In order to let the ICSs converge to more than 0.7 eV and the rate constants up to more than 1000 K, a large number of J partial waves contributions should be calculated. We calculated explicitly at intervals of ΔJ = 5 and 20 for J ∈ [0, 20] and [20, Jmax] respectively. The probabilities of uncalculated partial waves were obtained by means of interpolation over the J space. The accuracy of this method has been proven in previous work.50,58,59
In the previous studies, the disparity between experimental and theoretical results mainly exists in the v = 2 state.31 Therefore, we focus on the reactivity of this state first. The reaction probabilities for three different partial waves, J = 0, 40, 100, for the v = 2 state are displayed in Fig. 3 and compared with the results obtained on the TTCPP-B PES. Although the new PES has a lower barrier which energetically lies below the eigen-energy of v = 2 state, a threshold still appears at about 0.04 eV due to the endothermic nature of the reaction and the reaction probabilities show monotonically increasing energy dependence. On the other hand, in the TTCPP-B PES, the effective barrier height for the v = 2 state is 0.08 eV which is in good correlation with its threshold energy. The sources of the threshold energies are different with the two PESs. The better agreement between theory and experiment of the relative rate constants (which shown below in Fig. 6) indicates that the new developed PES is more accurate. From this point of view, the threshold energy of this work is more accurate. Because of the lower barrier and endothermicity, the probability obtained using the new PES is larger than the old one obviously for every partial wave. The threshold energy of each newly calculated partial wave probability is smaller than the corresponding value based on TTCPP-B PES and the difference between the corresponding two thresholds increases with J. Otherwise, all the reaction probability curves show oscillatory energy dependence and the oscillation is more pronounced in the probability curves obtained with the new PES, which is probably due to the deeper well in the reactant valley.
![]() | ||
Fig. 3 Initial state specified HF (v = 2, j = 0) total reaction probabilities for several selected Js with two different PESs (solid lines for the new PES and dash lines for TTCPP-B PES31). |
For the v ≥ 3 states, whose eigen-energies are greater than the endoergicity, the reaction probabilities are threshold-less and remain large throughout the energy range as shown in Fig. 4. Due to the deeper well, the oscillatory structures are much more pronounced in the probability curves calculated with the new PES similar to the v = 2 case and the probabilities of J = 0 obtained with the new PES are smaller than those obtained with TTCPP-B PES. The probabilities shift towards higher energies as J increases due to the additional centrifugal barrier. As the figure shows, the threshold energies increase more slowly with J in the results obtained with the new PES, leading to lower threshold energies for J > 0 probabilities compared to the results obtained with TTCPP-B PES. It may due to the relatively long R (3.81 vs. 3.67 a.u.) of the transition state in the new PES which is called “corner cutting” for reactions of heavy–light-heavy systems.60,61 So, we need more J partial waves to make the ICSs converge using the new PES.
![]() | ||
Fig. 4 Initial state specified HF (v = 3, 4, 5, j = 0) total reaction probabilities for several selected Js with two different PESs (solid lines for the new PES and dash lines for TTCPP-B PES31). |
In Fig. 5, the initial state specified ICSs are displayed and compared with the previous values. At a first glance, the initial vibrational excitation is effective in promoting the reactivity for this late-barrier reaction. Although the new PES has a lower barrier height, the calculated ICSs show similar characteristic to the previous results. The ICS for v = 2 state increases monotonically from the threshold energy while those for v ≥ 3 states are threshold-less and decrease sharply from very large values. However, differences still exist between the two sets of data. The new obtained ICSs for a given vibrational state are greater than the corresponding one obtained with TTCPP-B PES, especially in the low energy region. This enhancement of reactivity may be due to the increase of the cone of acceptance. In addition, the resonance structures are not fully canceled out by partial wave summation at low collision energies and the locations of the peaks are consistent with each other for different initial vibrational states which presumably suggests the same origin of them. On the other hand, the differences between the two PESs seem to have slight effects on the reactivity of the threshold-less reactions at high collision energies and their ICSs nearly the same at the collision energies greater than 0.4 eV.
![]() | ||
Fig. 5 Initial state specified HF (v = 2, 3, 4, 5, j = 0) integral cross sections for the Na(3S) + HF reaction (Solid lines for the new PES and dash lines for TTCPP-B PES31). |
To test the accuracy of the newly constructed PES, in Fig. 6, the initial state specified rate constants are computed and compared with the experimental data2 as well as our previous results using TTCPP-B PES.31 The data shown in the figure is relative rate constants with respect to kv=3 (i.e. kv/kv=3) at T = 940 K to mimic the experimental conditions. As displayed in the figure, the newly calculated rate constants for v > 3 states are in excellent agreement with the experimental results. The kv=4/kv=3 ratio is exactly equal to the average experimental value and the kv=5/kv=3 ratio is exactly equal to one of the experimental values which is closest to the average value. Unfortunately, the new PES still underestimates the kv=2/kv=3 ratio. In our previous work, we had artificially modified the TTCPP-B PES to reduce the barrier height and found it could narrow down the theory–experiment disparity effectively.31 In this work, the barrier height of the new PES is lower than that of the TTCPP-B PES by 0.268 eV. However, the kv=2/kv=3 ratio increases only by a factor of 6.36 (0.089 vs. 0.014) and is still much smaller than the experimental value (kv=2/kv=3 = 0.65–0.96), partly because of the non-zero threshold energy of v = 2 state. Another possible reason for the disparity is the role of initial rotational excitation. The experimental measurement is carried out at T = 940 K at which the maximally populated rotational state is j = 4. In order to clarify the effect of rotational excitation, we perform additional dynamics calculations for other four initial rotational states (j = 2, 4, 6, 7) of the v = 2 manifold. The eigen-energies of these state are 1.250, 1.284, 1.336 and 1.369 eV, respectively. The eigen-energy of j = 2 state is still smaller than that of the product. The calculated ratios for these states also displayed in Fig. 6. As the figure shows, the initial rotational excitations promote the reactivity as the ratios increase with j and up to 0.24 at j = 7. However, the ratios are still smaller than the experimental result. In the experimental side, the possibility of error in the measurement cannot be completely excluded. If we accept the experimental result, the theory–experiment disparity shown for v = 2 state calls for further study.
![]() | ||
Fig. 6 Initial state specified relative rate constants (kv/kv=3) of the Na(3S) + HF (v = 2–5) → NaF + H reaction with respect to the v = 3 state (kv=3) at T = 940 K. + and red open symbols denote the experimental values measured at different conditions, and the solid and dashed lines denote the average experimental values.2 The theoretical results calculated on the new and TTCPP-B31 PESs are denoted as solid circles and solid squares, respectively. Black hollow triangles denote the values of different initial rotational states (j) of v = 2 state on the new PES. |
Quantum dynamics studies are carried out using Chebyshev wave packet method to verify the accuracy of the new PES. The initial vibrational excitations enhance the reactivity significantly in line with the experimental observations. The oscillation is more pronounced in the newly calculated probabilities of every initial state due to the deeper well. The centrifugal barriers of J partial waves increase more slowly using the new PES due to the “corner cutting” for reactions of heavy–light-heavy systems. For v ≥ 3 states, the reactions are threshold-less and the newly calculated ICSs are obviously greater than the old ones in low energy area then tend to be the same in the high collision energy region. The difference in the details of the PESs seems to have no effect on the ICSs of the threshold-less reaction at high collision energy. The better agreement between the relative reaction rates and the experimental values for v > 3 states indicates that the new PES is more accurate. We pay more attention to the v = 2 state because of the theory-experiment disparity shown in the previous work. Although the new PES has a lower barrier, the probability of v = 2 state still has threshold because v = 2 state energetically lies below the product. For v = 2 state, the new PES gives enhanced reactivity compared to the TTCPP-B PES. However, the relative reaction rate is still much smaller than the experimental value. The consideration of rotational state excitation, which improves the agreement though, is not enough to remove the theory-experiment disparity. The reason for this theory–experiment disparity is still unclear and calls for further studies.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ra01885g |
This journal is © The Royal Society of Chemistry 2023 |