Polyatomic ab Initio Complex Potential Energy Surfaces: Illustration of

Article pubs.acs.org/JCTC

Polyatomic ab Initio Complex Potential Energy Surfaces: Illustration of Ultracold Collisions Debarati Bhattacharya,*,† Anael Ben-Asher,† Idan Haritan,† Mariusz Pawlak,‡ Arie Landau,† and Nimrod Moiseyev*,†,§,∥ †

Schulich Faculty of Chemistry, §Department of Physics, and ∥Russell-Berrie Nanotechnology Institute, Technion-Israel Institute of Technology, Haifa 32000, Israel ‡ Faculty of Chemistry, Nicolaus Copernicus University in Toruń, Gagarina 7, 87-100 Toruń, Poland S Supporting Information *

ABSTRACT: Resonances are metastable states that decay after a finite period of time. These states play a role in many physical processes. For example, in recent cold collision experiments, autoionization from a resonance state was observed. Complementing such observations with theory provides insight into the reaction dynamics under study. Theoretical investigation of autoionization processes is enabled via complex potential energy surfaces (CPESs), where the real and imaginary parts, respectively, provide the energy and decay rate of the system. Unfortunately, calculation of ab initio polyatomic CPESs are cumbersome; hence, they are not in abundance. Here, we present an ab initio polyatomic CPES utilizing a recently developed approach, which makes such calculations feasible. This CPES helps interpret the autoionization process observed in the He(23S) + H2 collision. From the behavior of the calculated CPES we can conclusively determine the nature of the autoionization process. Moreover, this CPES was used to generate reaction rates for the collision of He with ortho- and para-H2. These reaction rates are obtained from first principles. The results show a remarkable agreement with the cold collision experimental measurements, which demonstrates the robustness of our method. Hereby, we provide a computational tool for designing and interpreting new types of experiments that involve resonance states, e.g., in nucleobase damages (DNA or RNA) or in interatomic (intermolecular) Coulombic decay.

1. INTRODUCTION Autoionization is a process in which a molecular system in a resonance (metastable) state decays into a cation and a free electron. The system, however, does not ionize immediately, although it has sufficient energy to do so. Such a resonance state has a finite lifetime. Autoionization has been experimentally observed in molecular scattering,1 cold molecular collisions,2−4 and interatomic (intermolecular) Coulombic decay (ICD)5,6 processes. In addition, it has been indicated that such a process cause damages to DNA or RNA nucleobases due to the generation of slow moving electrons.5,7−9 Complementing such experiments with theoretical and computational studies is essential for analyzing and interpreting these observed processes. To aptly describe an autoionization process, in principle, the time-independent Schrödinger equation is solved by imposing outgoing boundary conditions.10 Under these conditions the Born−Oppenheimer Hamiltonian has a discrete electronic spectrum with real and complex eigenvalues. The real eigenvalues are associated with bound states, whereas the complex eigenvalues represent resonance states. The real and imaginary energy parts of the complex eigenvalues correspond to the position (relative to the ionization threshold) and lifetime (inverse decay rate) of the resonance states, respectively. Complex potential energy surfaces (CPESs) are © 2017 American Chemical Society

obtained from these complex eigenvalues as a function of the geometrical configuration of the molecular system. Specifically, these calculated CPESs are the ones required for the interpretation of the observed resonance phenomena. Practically, calculations of CPESs are not trivial. Imposing outgoing boundary conditions within an ab initio approach is very difficult. Conceptually, a simple approach to calculate atomic electronic resonances is complex scaling (CS),11−20 yet this procedure suffers from numerical and technical problems.21 Moreover, the description of molecular resonances within CS is not straightforward. A more convenient and popular approach is augmentation of the molecular Hamiltonian with a complex absorbing potential (CAP).22−37 Of note, because the CAP operator is artificial, its effect needs to be removed.38 Implementing CS and CAPs requires modifications of standard quantum chemistry packages (SQCPs). Nowadays, most computer implementations use CAPs but are limited to a few electronic structure methods.33,35,39 An alternative approach is based on the stabilization technique40 followed by an analytical continuation from the real to the complex plane.41−46 It makes use of the existing SQCPs without any need for modification. In the stabilization Received: January 26, 2017 Published: March 13, 2017 1682

DOI: 10.1021/acs.jctc.7b00083 J. Chem. Theory Comput. 2017, 13, 1682−1690

Article

Journal of Chemical Theory and Computation

in a proper manner for the general case of metastable decay rates. The reason for choosing this particular system is to demonstrate the accuracy of our method for which experimental cold-chemistry reaction rates are available. At very low temperatures, processes such as autoionization are well pronounced. Moreover, the sensitivity of these quantum resonances to the CPES structure poses a challenge to any state-of-the-art ab initio calculations. Therefore, to test our RVP method, we use the obtained ab initio CPES to compute the reaction rates for the collision of excited triplet helium with ortho- and para-hydrogen molecules. These reaction rates are calculated from first-principles, where only the Planck’s constant, charges, and masses of the electrons and nuclei were used as input parameters. Our approach yields remarkable agreement with the measured reaction rates. The ability to obtain the experimental result without a priori assumptions or the use of any fitting parameters paves the way for designing new kind of experiments, where the decay rates can also be large, such as in an interatomic (intermolecular) Coulombic decay (ICD) process.5,58

technique, the eigenvalue spectrum of the Hamiltonian is plotted with respect to the gradual variation of the basis set in real space.47 Thus, stabilization graphs can be obtained by any SQCP, with any existing Hermitian electronic structure method and at any level of theory. The stable region of this eigenvalue spectrum provides an estimation of the resonance position. Nevertheless, calculation of the resonance’ decay rates using this approach is more complicated. Recently, Landau et al.46 showed that the stabilization method, combined with an analytical continuation via the Padé approximant, can correctly account for atomic and molecular resonance positions as well as decay rates. This technique is abbreviated as resonances via Padé (RVP). In this report, we present a high-level electronic structure based RVP calculation of a highly excited CPES of a polyatomic supermolecule. Specifically, here we provide a theoretical investigation of the autoionization process that results from an atom−molecule collision in the milliKelvin temperature regime. Recent remarkable progress in cold-collision experiments allows for new observations of quantum phenomena such as, autoionization, anisotropy effect, and quantum tunnelling.2−4,48 These quantum effects were observed when an excited helium atom collides with molecular hydrogen at ultracold temperatures.2 Complementing such experiments with theoretical and computational approaches will help analyze and interpret these observations. Previous works3,48,49 have also analyzed this autoionization process by a theoretical approach. However, these theoretical analyses were obtained by some approximations, which hold in this particular case but are not suitable in the general case. The real and imaginary parts were not calculated concurrently.3,49 The real potential energy surface (PES) was obtained by a single-reference coupled cluster approach, designed to treat bound electronic states, while the imaginary part was fitted from the Tang−Toennies model (see ref 3). In addition, the reaction rates were calculated from the phase shift between the amplitudes of the outgoing and incoming waves in the S-matrix (see eqs 23 and 24 in ref 50 and section 3.5 in ref 10). This approach holds only when the phase shift is real. Thus, when the decay rate is very small, this approach serves as a qualitative approximation (again, see section 3.5 in ref 10). With these approximations, a quantitative agreement with the experiment was obtained by using a fitting parameter. Specifically, the correlation energy had to be shifted by 0.4% throughout the calculated real PES.48 Other CPESs were calculated in a fashion similar to the strategy employed in ref 3, in which the two components, real and imaginary, of the CPES are calculated by two separate techniques. For example in refs 51−55 the real part was approximated via a configuration-interaction calculation, while the imaginary part (the width) was extracted from complex Kohn scattering calculations. Another example is the FanoADC-Stieltjes approach, in which the real part is calculated by algebraic diagrammatic construction scheme up to second order [ADC(2)] and the decay rates by using the Fano theory of resonances.56,57 Alternatively, CPESs were also calculated by CAPs within ADC(2) (where the two complex parts were treated concurrently).26,27 With everything considered, ab initio CPES are not abundant. In this study we present a high-level electronic structure CPES for highly excited Feshbach resonance state of a polyatomic system, in which the real and imaginary energy components are treated on equal footing. Our work aims to demonstrate that interpretation of autoionization processes of polyatomic systems can be achieved

2. DETAILED THEORY WITH RESULTS Explicitly, we theoretically investigate the following molecular reaction, He(3S,1s2s) + H 2 → [He*−H 2] → He(1S,1s2) + H 2+ + e−

(1)

to gain more insight into the observed autoionization process.2 The reactants are an excited helium atom and ground-state molecular hydrogen. These two subsystems form a resonance (metastable) state, which decays upon collision. The products are helium in its ground state, an ionized hydrogen molecule, and a free electron. A schematic representation of the He*−H2 system (supermolecule) is shown in Figure 1, where He*

Figure 1. Schematic representation of the He*−H2 supermolecule.

corresponds to He(3S,1s2s). The molecular hydrogen bond length is r0 and the distance between He and the center of mass of H2 is represented by R. The angle between R⃗ and r0⃗ is ϕ. Here, we detail the different steps of the calculation. In the first step, the potential energy curve for the supermolecule was generated in the real plane. In the second step, the stabilization technique was used. This was followed by analytical dilation of an energy function from the real to the complex plane using the Padé approximant. Resonances were identified as stationary points in the complex plane. In the final step, the CPES was generated by calculating the complex eigenvalues as a function of the geometrical configurations. The calculated CPES is then used to evaluate the collision reaction rate. 2.1. Real Potential Energy Surface. Calculations of real or complex potential energy surfaces (PESs) require the consideration of different geometrical configurations. Figure 1 1683


Article


intensity of the red color corresponds to the intensity of the decay rates, where a darker shade corresponds to a faster decay, and a lighter shade to a slower decay. The area of interest, where molecular autoionization has been experimentally observed,2 is shown as an inset in Figure 2. The inset exposes a shallow potential well, which could have been overlooked on the larger scale (see the black rectangle on the blue curve in Figure 2). The depth of the well for T-shape and linear geometries appeared to be around 2.87 × 10−5 and 5.012 × 10−5 hartree (6.3 and 11 cm−1), respectively. This emphasizes the need for high accuracy in describing the supermolecule resonance state. Assurance regarding the accuracy of the calculated shallow well was provided from several convergence diagnostics of this result, which are described in details in the Computational Details and Supporting Information. Further confidence was provided by the comparison of the reaction rate, obtained from our CPES, which is in excellent agreement with the experiment (Figure 5 below). In the following subsection we shift from the real to the complex plane. 2.2. Resonances via Padé. In the previous subsection we showed that the state of interest is embedded in the continuum of the ionic species and, hence, is a resonance state. Therefore, its eigenvalues should be complex. To calculate the CPES of the state of interest, we used the RVP method.46 First, we generated the stabilization graphs, which were computed in the real space. Within this technique, the eigenvalues are evaluated when a given set of finite basis functions are scaled by a real factor. A plot of the eigenvalues (energies) as a function of the real scaling parameter is known as a stabilization graph. Continuum and resonance states comport differently when a finite basis is scaled by a small real factor. This behavior is attributed to the wave function associated with each state. The continuum states are associated with a delocalized wave function, whereas the resonance states are more localized in the interaction region. Therefore, the resonance eigenvalues would not be strongly affected by a small compression and/or expansion of the basis set. On the contrary, continuum eigenvalues depend strongly on the scaling parameter. The regions where the continuum eigenvalues attempt to cross the resonance eigenvalues are known as the avoided crossings. These features show up in the stabilization graph, which is presented in Figure 3a, for the system of interest in T-shape geometry at R = 1.5 Å. Such stabilization graphs were computed at each geometry. The black curves (Figure 3a) are the calculated eigenvalues as a function of the real scaling parameter α. In the given figure, three branches are shown, out of which one branch exhibit a relatively stable region of the eigenvalues. To calculate the resonance position as well as its decay rate, we then moved from the real plane (stabilization graphs) to the complex plane, by using analytical continuation. We used the Padé approximant, where an energy function, E (of a real scaling parameter, η), is fitted to a ratio between two polynomials (P and Q),

presents a schematic representation of the system under study. Numerous configurations of the supermolecule exist for different angles, ϕ and distances, R and r0. In this investigation, we assumed that the hydrogen molecule was a rigid rotor and, hence, kept r0 fixed at 0.74085 Å. The distance R varied over a wide range, while the angle was restricted to ϕ = 0 and ϕ = π/2. The computation began in the real plane by treating the He*−H2 state as if it was a bound state. To calculate this real PES, we used the implemented equation-of-motion coupled cluster (EOM-CC) method in the Q-Chem electronic structure package.59 The spin-flip (SF) variant of EOM-CC with singles, doubles, and perturbative triples correction [EOM-SF-CCSD(dT)]60 was employed by using the primitive 5ZP basis set,61 augmented with 3s, 3p and 6s, 6p diffuse basis functions on the hydrogen atoms and helium atom, respectively. Additional technical details and convergence tests are provided in section 3 and in the Supporting Information, respectively. Figure 2 displays the real potential curve for the He*−H2 supermolecule at ϕ = π/2 (T-shape). In addition, it presents

Figure 2. Real energy surfaces with schematic decay rates. Real potential energy curve of the neutral-excited (blue) and cationic species (green) at ϕ = π/2. The neutral-excited and cation total energies are in hartrees, with the intermolecular separation in angstroms. This is an approximation because resonance potential curves are complex. The intensity of the red color reflects the decay rate from the neutral-excited state to the cationic state, with stronger intensities reflecting higher decay rates. The inset zooms into the region around 6 Å, wherein a shallow well shows up. This is of interest as molecular autoionization has been experimentally observed in this region.

the PES of the He−H2+ product state. Clearly, the real PES of the supermolecule (in blue), lies above the ionization threshold, i.e., the PES of the ionized He−H2+ species (in green). This is an indication that He*−H2 is a resonance state embedded in the continuum of He−H2+. Thus, this real PES of the neutral excited system approximates the appropriate CPES. Similar trends and tendencies were observed for the ϕ = 0 (linear) geometry (Figure S1 in the Supporting Information). To properly describe the resonance state, we needed to move into the complex plane and identify the resonance complex eigenvalues. This process, which is described below, provides the decay rates that are directly proportional to the imaginary part of these eigenvalues. The decay rates as a function of R are schematically represented by the red arrows in Figure 2. The

E (η) =

P(η) Q (η)

(2)

The data for the fitting were only taken from the stable part of the stabilization graph. Before moving on to the complex plane, we verified the validity of the chosen set of data points and that it holds the relevant information. To do so, we reproduced the stabilization graph by using this chosen data set and the Padé 1684


Article


Figure 3. Stabilization curves, their reproduction, and cusps. (a) Exact and analytically dilated energy stabilization plot for the supermolecule in Tshape geometry at R = 1.5 Å. The blue area presents an example for the set of points used as an input for the Padé method. The magenta dashed line shows the analytically dilated energy values, as generated by the Padé method from the selected points. (b) The violet arrow points to the cusp, where both the θ and α trajectories meet. The figure represents the cusp at a distance of R = 1.5 Å for the system in the T-shape geometry. The α and θ trajectories are represented by the black and red curves, respectively. Cusps are associated with stationary points in the complex plane, i.e., resonances.

2.3. CPES of He*−H2. The procedure described in the previous subsections allows for calculation of a resonance state in a certain geometry. CPESs are obtained by repeating the same procedure at each molecular configuration. Specifically, this was executed at various distances, R, for both linear and Tshape geometries. The complex potential, real [Re E(R,ϕ)] and decay rate [Γ(R,ϕ) = −2 Im E(R,ϕ)], for the T-shape geometry (ϕ = π/2) of He*−H2 are given in Figure 4a,b, respectively. A similar behavior was also observed for the complex potential of the linear geometry (ϕ = 0, Figure S2a,b in Supporting Information). In Figure 4a it is clearly seen that the depth of the potential well remains unchanged after analytic continuation. This implies that the real PES provides a good approximation to the resonance position, which indicates that the approximation made in ref 48 is justified. The best fit for the calculated decay rate, Figure 4b, was found to be a single exponential curve (or linear in logarithmic scale, see inset). The relation between a single exponential function to the Penning ionization decay rate is discussed in ref 63. Therefore, we conclude that the investigated autoionization process is a Penning ionization. The CPES can be represented by a truncated interaction potential [E(R,ϕ) → V(R,ϕ)], expressed as a power series (in cos ϕ)50,64,65

approximant. A successful replication of the stabilization graph should be able to reproduce the stable branch of the curve and qualitatively predict the regions of avoided crossings as well as the other branches. The analytically dilated curve in Figure 3a (in magenta) successfully reproduced the original stable branch of the stabilization (in black) and correctly predicts the position of avoided crossings. This function was fitted with a small part of the original stabilization (marked in blue). This is a numerical proof that the fitted polynomial expansion is indeed an analytical energy function of the scaling parameter. Once an energy function was fitted (with eq 2), analytical dilation into the complex plane was performed by choosing η to be complex, i.e., η = αeiθ, where α and θ are real scaling parameters, with θ being the angle of rotation into the complex plane (as in CS10−20). Resonances are identified as stationary points in the complex plane which satisfy ∂E(η) =0 ∂α

and

∂E(η) =0 ∂θ

(3)

Graphically, these stationary points are associated with cusps in the α- and θ-trajectories.10,62 A θ-trajectory is generated by fixing α and varying θ over a range of values. Similarly, a αtrajectory is generated by fixing θ and varying α over a specified range. Convergence is achieved when these two trajectories form cusps that meet as shown by the purple arrow in Figure 3b. The real part of the cusp corresponds to the resonance position and the imaginary part to the decay rate (Γ). The occurrence of a cusp has been discussed in details in chapter 7 of ref 10. Finally, we checked for the stability of the complex resonance energies with respect to small modifications of the Padé input data set. It was observed that modification of the input data set (e.g., marked in blue in Figure 3a), resulted in minor variation of the resonance energies, which were orders of magnitude smaller than the calculated values. It is important to note that to obtain such a minor variation, only those data sets were chosen that successfully reproduced the original stabilization curve.

V (R ,ϕ) = Re(V (R ,ϕ)) −

i Γ(R ,ϕ) = 2

∑ Vj(R)Pj(cos ϕ) j

1 ≅ V0(R ) + V1(R ) cos ϕ + V2(R ) (3 cos2 ϕ − 1) 2

(4)

where, R⃗ and ϕ are the same as defined in Figure 1, Pj are the Legendre polynomials, Vj(R) are the amplitudes. Further simplification is enabled by the symmetric nature of the system, and thus the V1 containing term was neglected (ref 50). Therefore, only V0(R) and V2(R) have to be computed. By assigning, for example, ϕ to be 0 and π/2 in eq 4, we obtain the expressions for calculating V0(R) and V2(R): 1685


Article


Figure 4. Ab initio complex potential energy surface. (a) Real part of the complex energy for the system in T-shape geometry. The inset zooms in the region of the shallow potential well. The overlap of the orange curve (real component of the complex eigenvalues after analytic continuation) with the black curve (the real eigenvalues obtained from bound state calculations) implies that the resonance position is well approximated by the boundstate calculations for this state. (b) Rate of decay (inverse lifetime) from the imaginary part of the complex energy for the system in T-shape geometry. The inset highlights the decay rate in logarithmic scale. The decay rate fits into a single exponential curve, thus confirming the decay process as a case of Penning ionization.2,63 (c) The real part of the interaction potential in terms of V0 and V2 in Hartree. (d) Decay rate (inverse lifetime) of the interaction potentials in Hartree.

1 [V (R ,0) + 2V (R ,π /2)] 3 2 V2(R ) = [V (R ,0) − V (R ,π /2)] 3 V0(R ) =

ĤHe *−H2 = −

2μHe *−H

+ Re(V

(5)

The real parts of V0(R) and V2(R) are presented in Figure 4c where the depth of the wells are approximately 3.5 × 10

ℏ2 2

2 2 ĵ ∂2 l̂ + + ∂R2 2μHe *−H R2 2μH r0 2 2

2

i (R ,ϕ)) − Γ He *−H2(R ,ϕ) 2

He *−H 2

where R is the length between the He atom and the center of mass of H2, ϕ is the rotational angle between R⃗ and the H2 bond. j2̂ and l2̂ are the squared angular momentum operators of the diatomic and atom−diatom relative rotations, respectively. Then, the time-independent Schrödinger equation was solved,

−5

hartree (7.7 cm−1) and 2.28 × 10−5 hartree (5 cm−1). The decay rate components of the CPES are presented in Figure 4d. 2.4. First Principle Reaction Rates. The reaction rates

ĤHe *−H2 ΨiHe *−H2(R ,ϕ)

were calculated with the ab initio CPES (via V0(R) and V2(R)),

⎛ ⎞ i = ⎜ϵiHe *−H2 − Γ iHe *−H2⎟ΨiHe *−H2(R ,ϕ) ⎝ ⎠ 2

which allows direct comparison with the experimental results. To this end, the generated CPES, presented in eq 4, was

and the obtained eigenfunctions have complex eigenvalues, *−H2 *−H2 where ϵHe is the energy and ΓHe is the decay rate of the i i ith state, respectively. In addition, the time-independent

substituted into the nuclear Hamiltonian for the metastable He*−H2 state, 1686


Article


Figure 5. Reaction rates for the collision of He* with H2 (a) in the rotational ground state (para) and (b) in the rotational first excited state (ortho). The blue curves represent the experimental results reported in ref 48, and the red curves are the first-principles theoretical results, computed from the ab initio CPES, without using any fitting parameter.

Schrödinger equation of the nuclear Hamiltonian of the product He−H2+, ĤHe − H2+ = −

ℏ2 2μHe − H + 2

i. When He and H2 collide at a certain collision energy, Ecol, they form the resonance complex He*−H2 with *−H2 energy ϵHe , which is equal to Ecol. i

2 2 ĵ ∂2 l̂ + + ∂R2 2μHe − H + R2 2μH+ r0 2 2

*−H2 (=Ecol) ii. The supermolecule He*−H2 with energy ϵHe i +

2 , which is autoionizes to He−H2+ with energy EHe−H f

2

+

+ V He − H2 (R ,ϕ)

+

2 + Ekinetic . lower than Ecol, i.e., Ecol = EHe−H f e It is pertinent to mention here that the reaction rates were calculated on the basis of eq 6 and by using the adiabatic variational theory shown in refs 49 and 50. Comparison with the sub-Kelvin experimental reaction rates is presented in Figure 5. The blue and red curves in Figure 5, depict the experimental and our theoretical findings, respectively, which are in excellent agreement for both paraand ortho-H2 collisions with He*. Of note, our results are within the experimental uncertainty because the experimental absolute values in ref 48 were normalized according to the thermal rate measurement by Oskam and co-workers.67 The theoretical reaction rates were computed without any fitting parameters, and this manifests the accuracy of the calculated CPES. In addition, it illustrates the universality of the method in generating CPESs as well as reaction rates for any polyatomic system in any decay process.

was solved, +

+

+

− H2 − H2 ĤHe *−H2 ΨHe (R ,ϕ) = EfHe − H2 ΨHe (R ,ϕ) f f +

2 where EHe−H is the energy of the state f. The corresponding f eigenvalues and eigenfunctions of the metastable state and the product were integrated in the scattering theory to compute the reaction rates for the two collisions. In the first, the hydrogen molecule was in the rotational ground (para) state, and in the second, it was in the first excited (ortho) state. The equation for computing the reaction-rates (derived on the basis of Chapter 8 of ref 10 and ref 66) is given by

RR(Ecol) =

vcol kcol 2

P(Ecol)

(6)

3. COMPUTATIONAL DETAILS The EOM-SF-CCSD(dT)60 method was employed by using the primitive 5ZP basis set,61 augmented with 3s, 3p and 6s, 6p diffuse basis functions on the hydrogen atoms and helium atom, respectively. The 1S ground state of He−H2 was the reference configuration used to calculate the target 3S resonance state. The cationic ground-state PES was obtained by the same method and basis set, where the high-spin 4S state was the reference configuration. Most quantum chemistry codes and basis sets were designed to treat bound molecular states. To aptly describe a resonance state, a large number of diffuse functions are required. Details regarding the convergence with respect to the additional diffuse functions and cardinal number ζ are found in the Supporting Information, Tables S1 and S2. All the basis functions were atom centered Gaussian functions. The extra set of diffuse functions on the 5ZP basis set were added in an even tempered manner (with a spacing of 2.0). We also removed the tight binding g-functions from both the atoms, so as to reduce

where P(Ecol ) =

Ecol

∫−∞

+

+

dEfHe − H2 ρ(EfHe − H2 ) 2

×

− H2+ ΨHe (R ,ϕ) f

Γ He *−H2(R ,ϕ) *−H2 ΨHe (R ,ϕ) i 2π

is the dimensionless transition probability from the initial state i *−H2 = Ecol, to all the final states. The relative with energy ϵHe i

2Ecol μHe *−H

velocity of the reactants (He* and H2) is vcol =

, the

2

wave vector of He* and H2 is kcol =

2μHe *−H Ecol 2

ℏ

, and

He−H2+

ρ(Ef ) is the density of states of the product He−H2+. The prefactor vcol/kcol2 guarantees that the RR has cm3/s units. This expression takes into account the conservation of energy by the following assumptions: 1687


Article


theory. Analysis of the imaginary part of this CPES demonstrate that the decay rate for this process best fit a single exponential curve. This confirms that the investigated autoionization process is a Penning ionization.63 The ab initio CPES was then used to calculate the reaction rates for the collision of He with para- and ortho-H2. The excellent agreement with the experimental reaction rate validated that an appropriate level of theory was used, particularly because no fitting parameters and no a priori assumptions were invoked. Moreover, it shows that RVP allows for accurate CPES calculations. This innovative approach is crucial because CPESs are essential for interpreting and analyzing many quantum phenomena, particularly for cold collisions. The potential applicability of the RVP approach is thus manifold. First, it is not restricted to any particular standard quantum chemistry package and is straightforward to apply. Second, it is applicable to a large variety of molecules, ranging from small polyatomic molecules to medium-sized nucleobases. Moreover, it is not limited to a specific mechanism of the autoionization process, e.g., it can be applied to study ICD processes with lifetimes on the order of femtoseconds,58 as well as to nucleobase damages.7−9 Finally, the RVP method is extremely accurate and enables the calculation of experimental observables without the need of any fitting parameter.

computational cost without much loss of accuracy. The basis was also tested for basis set superposition error (BSSE), which can have large values if the basis/method is not appropriately chosen. The calculated BSSE values were negligible, on the order of 1 × 10−8 hartree for all distances. Other contributions, namely, relativity and diagonal Born−Oppenheimer corrections, were found to be negligible.48 Confidence regarding the accuracy of the calculated shallow well was provided by the minor contribution of the perturbative triples (on top of the singles and doubles model), which appeared to be about 4.55 × 10−6 hartree (1 cm−1) only. In addition, we find that the full configuration interaction (FCI) correction48 to our well depth is less than 1 cm−1. The FCI values were calculated with a smaller basis set, as currently FCI calculations of the CPES with a large basis set is practically unfeasible. Finally, the existence of the shallow well is essential for the excellent reproduction of the experimental reactionrates as demonstrated in the Supporting Information (Figure S3). Once the real PES was generated, we scaled our basis functions by a real parameter as the first step of the stabilization calculations. In our case, the scaling was restricted to the diffuse functions only; i.e., the functions whose exponents were smaller than a certain critical value (in this case

Polyatomic ab Initio Complex Potential Energy Surfaces: Illustration of

Recommend Documents