Nonrigidity effects — a missing puzzle piece in the description of low-energy anisotropic molecular collisions

Mariusz Pawlak, ∗ Piotr S. Żuchowski, † Nimrod Moiseyev, and Piotr Jankowski Faculty of Chemistry, Nicolaus Copernicus University in Toruń, Gagarina 7, 87-100 Toruń, Poland Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University in Toruń, Grudziądzka 5, 87-100 Toruń, Poland Schulich Faculty of Chemistry and Department of Physics, Technion–Israel Institute of Technology, Haifa 32000, Israel

The breakthrough in controlling the movement and internal degrees of molecules with external fields [1], which started about 20 years ago, currently allows to study the collisions, reactivity and properties of molecules in very cold regime. Such molecules provide new and propitious prospects in precision spectroscopy, fundamental physics, astrochemistry, and quantum engineering. Understanding of quantum effects, resonance phenomena, and reaction dynamics in a low-energy range opens the gate to design and create materials with unusual functionality and elements of quantum computers [2,3]. Cold collision experiments realized by merging two supersonic beams have become an important technique for studying chemical reactions in temperatures near 1 K [4][5][6][7][8][9] in which unseen earlier quantum features, such as resonances or interference, are revealed. The group of Narevicius performed the first experiment in that field focused mainly on the Penning ionization (PI) process of colliding hydrogen isotopologues with excited metastable helium atoms [4,5,8,9]. In such a reaction, an electron is moved from the molecule to the only partially-occupied orbital of excited atom, then the initially excited electron of the atom is kicked out of the system. Three products occur: the atom in the ground state, the molecular ion in the ground state, and a free electron. These colliding systems are of great interest to astrophysicists and astrochemists studying conditions and reactions in outer space. Hydrogen and helium are the most abundant elements, whereas molecular hydrogen is the most common molecular species in the Universe [8,10]. The first observation of metastable helium atom in the atmosphere * e-mail: teomar@chem.umk.pl † e-mail: pzuch@fizyka.umk.pl of one of exoplanets [11] has boosted its importance for astrochemistry, and one can expect that its interaction with the omnipresent hydrogen molecule will be carefully investigated.
Recently, Narevicius and co-workers directly probed the anisotropy in atom-molecule interactions through orbiting resonances by changing the rotational state of the molecule [9]. That work reveals a crucial role of the anisotropy of the interaction energy, due to various orientations of the H 2 molecule in the complex of He(1s2s, 3 S 1 ) (≡ He * ) with H 2 , in the dynamics in sub-kelvin regime. To elucidate physical phenomena presented in their novel experiment, the authors of Ref. [9] used the state-ofthe-art first-principles calculations. The interaction energy of the complex comprised of the metastable helium atom and the hydrogen molecule was first calculated using the supermolecular approach employing the coupled cluster method with singles, doubles and perturbative triples (CCSD(T)), and is further denoted as E The CCSD(T) method is known as the gold standard in quantum chemistry and in many cases provides values of properties accurate enough to predict experimental results [12]. However, in the case of He * +H 2 , the values of reaction rate coefficients calculated with the potential based on E CCSD(T) int did not agree with experiment even qualitatively. Thus, the correction to the interaction energy, denoted here by δE FCI int , was calculated from the full configuration interaction (FCI) method and added to the E CCSD(T) int energy. The agreement of the rate coefficients calculated from the E CCSD(T) int +δE FCI int energies with the experimental values was improved but still qualitatively incorrect, since one additional resonance, not present in the experiment, was predicted for very low collision energies. To obtain quantitative agreement of the calcu-lated rate coefficients with the experiment, the authors scaled the correlation part of the interaction energy by a factor of 1.004, that can be viewed as adding the correction 0.4%E corr int to the E CCSD(T) int +δE FCI int energy, suggesting that the basis set incompleteness was responsible for the inaccuracy of the rigorous ab initio interaction energy surface.
The motivation of the present studies was to find a reason why the ab initio interaction energy surface, even obtained at the FCI level of theory, was not accurate enough to precisely reproduce the experimental results. Here, we show that the missing piece of the puzzle is the nonrigidity of the H 2 molecule, neglected in the theoretical study of Ref. [9]. It has been recently shown that taking into account the monomer nonrigidity effects is necessary to obtain very precise agreement with the spectroscopic or scattering experiments [13][14][15][16][17]. In the present paper, we demonstrate for the first time the importance of the monomer flexibility in low-energy molecular anisotropic collisions. Only by inclusion of vibrations of the molecule in description of the complex, we are able to correctly predict results of very subtle scattering experiments with no fine-tuning to the experimental data whatsoever. Moreover, we present how to incorporate, in a simple and effective way, the nonrigidity effects into theoretical studies for any colliding molecules.

RESULTS AND DISCUSSION
Theoretical consideration of our problem can be divided into two steps: first, the preparation of the most reliable interaction energy surface and second, scattering calculations. Since the complex comprises the metastable atom He * and the diatomic molecule H 2 , the positions of the nuclei can be described by three coordinates: the distance R between He * and center of mass (COM) of H 2 , the angle θ between the H 2 bond and the COM-He * direction, and the distance r between the hydrogen nuclei. Thus, in principle one has to use the three variables (R, θ, r) to parametrize the interaction energy surface and next to perform the scattering calculations. However, even for relatively small atom-diatom systems the full-dimensional treatment is rather rare. In most applications, the quantum scattering calculations are performed within the rigid-rotor approximation, i.e., assuming that the molecules in the complex are rigid. This widely used approximation is well physically motivated, since the internal vibrations of interacting molecules are much more energetic than the intermolecular modes [13]. Usually, in the rigid-rotor calculations the rigid monomer interaction energy surfaces are used, obtained from ab initio calculations for the monomers with fixed geometries. Such calculations have been employed also in Ref. [9] to study the PI reaction of He * +H 2 . However, it has been shown very recently [16] that if in the rigidrotor scattering calculations one uses the interaction energy averaged over the vibrations of the monomers, the results are closer to the full-dimensional calculations and experimental data. It has been demonstrated that the vibrationally averaged surfaces perform better than the rigid-monomer ones also in predictions of other physical properties, like rovibrational spectra [13-15, 18, 19] or virial coefficients [20,21]. A main drawback of such approximation is that, in principle, to obtain the vibrationally averaged surface, one has to know the corresponding full-dimensional one. To avoid construction of a full-dimensional surface and minimize the number of geometries for which ab initio calculations have to be performed, we use the method developed in Ref. [22]. The details of this method are presented in the Methods section. Thus, in the current study we try to capture the nonrigidity effect by using the interaction energy surface averaged over vibrations of H 2 and rigid-rotor scattering calculations. We have to emphasize that although the averaged surface depends only on the intermolecular coordinates (R, θ), it cannot be regarded as the rigid one, since to obtain it one has to perform ab initio calculations for various values of the intramolecular distance r. Such a surface includes in an effective way information about the nonrigidity effects of the complex.
The most common technique to obtain PI rate coefficients is to use the complex potential in which the imaginary part describes the losses due to the ionization process [23,24]. With such potential one can solve the Schrödinger equation for the nuclear coordinates, for instance, using the close-coupling scattering method. More recently, two of us developed a new approach based on adiabatic theory and scattering theory for cold collision experiments [25,26] dubbed as adiabatic variational theory (AVT). AVT together with the dedicated new closedform expression for PI rate coefficients [27] provides a relatively simple method to implement and was used in our investigation. This technique, where the diatom is treated as a rigid rotor, allows to uncouple the rotations of the diatom and the complex from the atom-molecule separation.
In Fig. 1 we present the reaction rate coefficients of He * with para-H 2 in the ground rotational state (j = 0) and ortho-H 2 in the first excited rotational state (j = 1). The strong effect of anisotropy on the resonant structure was discussed in details in Refs. [9,26], however it is worth mentioning that in the interaction between He * and para-H 2 (j = 0) only the isotropic part of the potential is probed, because the wavefunction of the molecular hydrogen in the lowest rotational state is spherically symmetric. When the interacting H 2 molecule is in the j = 1 rotational state, the leading term of anisotropy of the potential contributes directly into the effective interaction and can firmly affect the positions of resonances. As demonstrated in Ref. [26], by excluding from calculations the orientation-dependent part, the low-temperature resonance at collision energy of k B × 0.27 K is missed in the theoretical results. It shows that this peak, as opposed to the peak at 2.4 K, arises totally from the anisotropic interaction. It takes place when the excited helium atom The results have been convoluted with the experimental energy spread. Neither scaling nor fitting parameters have been used in the calculations. The experimental values are taken from Ref. [9] (black points with error bars). Theoretical results are slightly below the experimental data, however, the latter are burdened with systematic normalization error larger than the vertical discrepancy. collides with ortho-H 2 (j = 1). Hence, the rate structures are entirely different for different rotational states of the hydrogen molecule.
The evidence of monomer nonrigidity effects in cold anisotropic collisions is demonstrated in Fig. 1. It is clearly seen that the results with the interaction potential from the CCSD(T) method supplemented by the FCI correction are still not satisfactory. We found out that the discrepancy between the reported theoretical and experimental data is not due to the incompleteness of the used basis set as the authors of Ref. [9] stated, but due to the assumed stiffness of the molecule. By adding the correction corresponding to the flexibility of the diatom, Reaction rate coefficients of He(1s2s, 3 S1) with H2 with respect to relative energy (in K) between the colliding subsystems. The theoretical rate coefficients have been calculated based on the interaction potential obtained using CCSD(T) with additional corrections resulting from FCI and inclusion of the effect of the hydrogen molecule nonrigidity. The black solid curve represents results for 100% of para-H2 (j = 0), whereas the magenta dashed-dotted curve represents results for the mixture of 98% of para-H2 (j = 0) and 2% of ortho-H2 (j = 1). The inset exhibits the region with a small bump in millikelvin regime. The experimental data are shown as black points with error bars [9].
termed as δE flex int , where the vibrations are averaged, we obtained an excellent agreement with the measurements over the whole range of temperatures. The nonrigidity correction shifts the energy of the resonance at 0.27 K to the position matching the experimental data. Also the magnitude of the rate coefficient curve is about 25% larger than the one calculated without the nonrigidity correction. These results are very similar to those obtained in Ref. [9] with the artificial 0.4% increase of the correlation part of the interaction energy. Note that in our entire calculations we did not apply any scaling or fitting parameters as well as we did not shift the final results to adjust to the experiment. Our results are slightly below the experimental ones, but the latter have been normalized to the absolute scale according to thermal rate observations at 300 K, see the Methods in Refs. [8,9], and this procedure introduced a systematic error to the experimental data much larger than the vertical discrepancy. In the Supplementary Information we provide a figure corresponding to Fig. 1, with the rate coefficient curves shifted by a constant value to match the normalized experimental data at the collision energy around k B ×2.4 K, as it was done in Ref. [9]. After such additional "normalization", the theoretical resonance structure around 0.27 K perfectly agrees with the experiment.
In the experimental data for para-H 2 , given in the upper panel of Fig. 1, one can see a very flat bump around the collision energy of k B ×0.3 K, not discussed in Ref. [9]. On the basis of our theoretical considerations we can explain an origin of that feature. One of the beams used in the Narevicius' experiment is formed of the para-H 2 molecules. However, a purity of para-H 2 was limited to 98% [9]. Therefore, we have added 2% of the reaction rate coefficients of ortho-H 2 (j = 1) given in the lower panel of Fig. 1 to 98% of the reaction rate coefficients of para-H 2 (j = 0) given in the upper panel of Fig. 1; the resulting rates are presented in Fig. 2. Now one can see, in the inset of this figure, that the bump appears on the theoretical curve, properly predicting the position and the shape of its experimental counterpart.

CONCLUSIONS
Low-energy collision experiments allow to directly and clearly reveal the true nature of interaction. Only at a high level of theory we are able to predict or confirm experimental data as well as understand quantum phenomena in chemical reactions. In this paper, we have demonstrated a significant role of monomer nonrigidity effects on the position and intensity of scattering resonances in anisotropic cold molecular collisions. We have investigated excited metastable helium atoms colliding with hydrogen molecules in the temperature range from a few dozen kelvins to 1 millikelvin. We have provided the most accurate interaction energy surface that can be treated as the reference one for novel quantum chemistry methods for molecular systems in the resonance state. The calculated rate coefficients are in remarkable agreement with measurements [9]. We have demonstrated that the discrepancy between the experimental and theoretical results discussed in Ref. [9] is due to the assumption that the H 2 molecule is rigid. We have complemented the theoretical description of the interaction energy of the complex by inclusion of its dependence on the flexibility of the molecule. It shows that the approach beyond the commonly used rigid-rotor approximation is indispensable even when rigorous state-of-the-art computations are performed at the FCI level of theory. The correction to the interaction energy due to the nonrigidity of the monomer is a few times larger than the uncertainties caused by basis-set incompleteness or generated by the Born-Oppenheimer approximation. Our finding is not limited to a specific complex considered here nor restricted to the PI reaction process. We believe that it can be vitally important in precisely controlled cold chemistry, where quantum effects in chemical reactions dominate.

Vibrationally averaged surface
To construct an interaction energy surface averaged over the vibrations of the monomers, V , one can take advantage of the fact that the molecules in the complex preserve their identity and frequencies of internal vibrations are much higher than those of the intermolecular modes, and perform averaging over the intramonomer coordinates similar to the adiabatic approximation in the electronic-structure theory [13]. Since we are interested in the He * +H 2 complex, let us limit our consideration to the atom-diatom case. If the Jacobi coordinates are used, the interaction energy surface V (R, θ, r) can be represented as the truncated Taylor expansion around some reference geometry r c [22] r=rc . The higher order terms can be neglected if only modest deformations of the monomer are allowed, as those corresponding to a few lowest vibrational states of the H 2 molecule. The f 0 function is in fact the rigid monomer two-dimensional surface calculated for the H-H separation equal r c , while the remaining terms account for the nonrigidity effects, i.e., the dependence of the surface on the intramolecular coordinate. The potential V TE from Eq. (1) can be easily averaged over the v vibrational state of the monomer and the resulting formula reads We can use the V TE surface as an approximation to the V one, i.e., assume V ≈ V TE . The values of r v and r 2 v can be calculated from the theoretical properties of the monomer or even estimated from the empirical spectroscopic constants. The surfaces averaged according to the approximation of Eq. (2) and the corresponding formula for the diatom-diatom case, were used in the rigid-rotor dynamical calculations, both the bound-states and the scattering, and provided the results in excellent agreement with the experimental ones [14-16, 18, 19]. In practical applications, there is no need to know the surface V TE for any values of the (R, θ) coordinates, but it would be enough to calculate it on the grid points, for instance the ones used in the scattering calculations. For each grid point in (R, θ), we can compute the interaction energy f 0 and the numerical values of the f 1 and f 2 derivatives. Of course, one can also calculate V TE for a given grid of geometries and then fit an analytical function to obtain the surface. From the expression (2), one can see a very useful feature, namely that the values of f 1 and f 2 may be calculated on a different level of theory than f 0 , e.g., the leading term f 0 on the highest possible level, and the values of f i defining the higher order terms, specifying the dependence of V TE on r, may be calculated at a somehow lower level of theory. Since the first and second order terms are much smaller than the leading term f 0 , such an additional approximation only slightly increases the uncertainty of V TE or V TE , whereas it may significantly reduce the computational effort required. Such a strategy has been applied, for instance, to the H 2 -CO complex in Refs. [14,15] and led to the very accurate rovibrational spectra.
To prepare the vibrationally averaged surface for the He * +H 2 complex, we have used the formula of Eq. (2) in the following way. The leading term f 0 was set to be equal to the rigid-monomer interaction energy of Ref. [9]. That interaction energy can be written as E CCSD(T) int +δE FCI int , using our notation, and was obtained as a sum of the interaction energy calculated at the CCSD(T) level and the FCI correction for the intramolecular distance 1.4487 bohr. Thus, we have to set r c to be equal to the same value to make our choice of f 0 consistent with Eq. (2). Since we consider the ground vibrational state of H 2 , v = 0, the values of r 0 and r 2 0 for j = 0 were set to the round-off values from Ref. [28] equal to 1.4487 bohr and 2.1270 bohr, respectively. The equality of r c and r 0 is of course not accidental, because the interaction energy of Ref. [9] was calculated for the rigid H 2 molecule with the vibrationally averaged geometry. With such a choice of r c , the first-order term in Eq. (2) vanishes and the nonrigidity effect is fully described by the second-order term. The values of r 0 and r 2 0 for the molecule in the first excited rotational state (j = 1) are slightly different: 1.4509 bohr and 2.1334 bohr, respectively. Therefore, for this case the first-order term in Eq. (2) contributes to the nonrigidity correction.
In the standard, supermolecular approach the interaction potential is very difficult to obtain due to the fact that our potential is not the ground-state one and is coupled to the scattering state of the He+H + 2 +e − system. In Ref. [29] it was shown that using a carefully tailored start guess it is possible to converge the CCSD(T) interaction potential and that also the symmetry-adapted perturbation theory (SAPT) [30][31][32] provides a good representation of the short-range potential. Nonetheless, it is very difficult to stabilize first and second derivatives of the interaction potential with respect to the nuclear coordinate motion in the supermolecular method as it inherently relies on subtraction of big numbers and loss of accuracy is unavoidable. We previously stated that the derivative of the interaction energy can be obtained at a lower level of theory and it still catches the essential physics. As a matter of fact, the application of SAPT greatly facilitates the calculation of nonrigidity effects. Since in SAPT the interaction energy is obtained directly from the wavefunction of monomers, the interaction energy is stable and very inexpensive. This is due to the fact that in SAPT we obtain the interaction energy directly as a sum of the perturbation theory terms in which the expansion parameter is the interaction potential between monomers. Here we use the interaction energy which is the sum of first two terms of perturbation series in the interaction operator between H 2 and helium analogously to Ref. [29]. The values of f 1 and f 2 were calculated from the four-and five-point central-difference formula, respectively, with h = 0.05 bohr. For each grid point in (R, θ), the interaction energy was calculated for five separations r c + kh, k = −2, −1, 0, 1, 2. The SAPT calculations were carried out with the d-aug-cc-pVTZ basis set. We estimated that uncertainty of calculation of flexi-bility correction in the interaction well is underestimated by about 2% for linear geometry and about 10% for Tshape geometry. In absolute numbers these uncertainties are below 0.01 cm −1 . For discussion of uncertainty of δE flex int see the Supplementary Information. In Ref. [9], the effect on the interaction energy surface beyond the Born-Oppenheimer approximation was neglected on the basis of the analysis of the value of the diagonal Born-Oppenheimer (DBO) correction [33] calculated at the minimum of the surface. Here we performed an extended analysis for two angular orientations, θ = 0 • and 90 • , several values of R, and the H 2 geometry fixed at r c . The DBO correction to the interaction energy, δE DBO int , was obtained by subtracting from the DBO correction calculated for the complex the value calculated for the monomer at large separations of interacting species. Finally, let us discuss the effect of basis set incompleteness. It is very difficult to estimate since the basis set extrapolation techniques are questionable to use in the present case, where the system is not in its ground state and the system concerned is essentially a resonance. The basis set used in the present calculations at the CCSD(T) level, aug-cc-pV6Z with bond functions, is already well saturated for the dispersion energy which dominates the interaction energy for the considered system. A computationally expensive increase of the basis set to aug-cc-pV7Z shifts the interaction energy in the global (linear) minimum by about −0.022 cm −1 , while for the local minimum at a T-shape geometry (for 10.85 bohr) the shift is about −0.027 cm −1 . Close to the inner turning points at 9.5 bohr, these values are −0.034 and −0.045 cm −1 , respectively. Unfortunately, one cannot perform an extrapolation to the complete basis set, since the resonant nature of the interaction affects the stability of ab initio calculations at the level of accuracy of the order of 0.01 cm −1 . Nonetheless, one should bear in mind that such an interaction energy should not be extrapolated in the usual sense but rather than that, stabilized. For this particular system, the coupling between the continuum and the bound state is weak, hence, the potential obtained by the standard quantum chemistry methods gives right answer for right reasons. An improvement of theory to go below the 0.01 cm −1 accuracy is a formidable task and a new approach would be needed to address such demands. Since the real potential is dominated by the dis- persion energy which is a variational quantity, the more complete basis set implies the deeper potential. However, if one proceeds from the aug-cc-pV6Z basis set to the augcc-pV7Z one, the shift of the interaction energy surface by 0.02-0.03 cm −1 at the minimum separation is an order of magnitude smaller than in the case of the δE flex int correction. Thus, the effect of the basis set incompleteness on the positions and intensities of the resonances is negligible in comparison to the effects caused by the δE flex int correction and is also much smaller than the effect of the δE DBO int correction. According to the works of Refs. [34,35], we may expand the vibrationally averaged surface V in Legendre polynomials For the collision of an atom with a homonuclear diatomic molecule the index η is even due to symmetry reasons, V (R, −θ) = V (R, θ). In other words, terms for odd η vanish. Thus, the two leading terms are given by V 0 (R) and V 2 (R)(3 cos 2 θ − 1)/2, where V 0 (R) and V 2 (R) are radial isotropic and anisotropic interaction potential terms, respectively. It should be emphasized that the considered complex is not in the bound state but in the resonance one. The total potential energy surface is above the ionization threshold, thus the electronic state of the He(1s2s, 3 S 1 )+H 2 system is embedded in the continuum of the He(1s 2 , 1 S 1 )+H + 2 +e − system. Consequently, the total potential energy surface has to be complex where the imaginary part represents the ionization width (inverse lifetime). Two new approaches have been lately developed for PI widths: one is based on Fano-algebraic diagramatic construction method [36], the next one uses the stabilization method with an analytical continuation via the Padé approximant [37,38]. In our studies, we took the imaginary part of V 0 and V 2 from Ref. [38], obtained by the latter technique.
The calculated isotropic, V 0 , and anisotropic, V 2 , radial interaction potential terms obtained on three levels of theory are presented in the upper panel of Fig. 3. One can see that adding the δE FCI int correction to E

CCSD(T) int
apparently changes both V 0 and V 2 , but in opposite directions: V 0 becomes deeper and V 2 slightly shallower. The subsequent addition of δE flex int to the E

CCSD(T) int
+δE FCI int interaction energy makes the resulting V 0 even deeper, whereas for the V 2 term the effect of δE flex int almost cancels the effect of δE FCI int . The lower panel of Fig. 3 presents how the δE FCI int and δE flex int corrections to the interaction energy enter the V 0 and V 2 terms of the interaction potential.

Adiabatic variational theory
To solve the Schrödinger equation for the complex with the previously prepared interaction potential and consequently to calculate rate coefficients, we used the AVT approach that has been developed for cold atommolecule collision experiments [25,26]. This technique has been recently successfully applied for the He(1s2p, 3 P 2 )+H 2 system [39]. It enables one to reduce the complexity of the problem enhancing the computational performance without losing physical essence. Within AVT, we represent the potential (3) in a basis set consisting of many angular functions. The single angular function in our case is a product of two spherical harmonics -one is responsible for the description of rotations of the molecule, whereas the other one for the description of rotations of the whole complex. Such a matrix constructed for a given intermolecular distance R needs to be diagonalized providing a set of eigenvalues. The process has to be repeated for different values of R. The obtained eigenvalues, after ordering, create so-called adiabats (effective potentials) depending on R. Then, from the practical reasons all adiabats are shifted to get asymptotes at zero. Therefore, for each of them the dissociation threshold is at zero. Next, we solve the one-dimensional Schrödinger equation many times, each time with a different adiabat treating R as a variable. Here, any technique can be used, e.g., by spanning the wave function space in a basis of trial functions [40,41]. Finally, we apply the simple and easy to implement formula for reaction rate coefficients that has been derived based on AVT and non-Hermitian scattering theory [27]. Only the information about complex eigenenergies, the reduced mass of the atom-diatom system, and the rotational state of the molecule are required. In calculations, we used 21 partial waves corresponding to end-over-end angular momenta of the complex (l = 0, 1, ..., 20). The Schrödinger equation was solved by the DVR with the box size of 500 bohr and with 2000 basis sine-functions. The results are fully converged with respect to the number of partial waves and of basis functions. The calculated reaction rate coefficients are convoluted over the experimental collision energy spread (8 mK). Neither scaling nor fitting parameters have been used in our calcula-tions.
In Table S1  In Table S1

S4
Let us now discuss the uncertainty of the flexibility correction δE flex int . In Table S2 we gathered second derivatives of the components of the interaction energy for 10.5 bohr for Tshape and linear geometry. As one can see, the dispersion energy by far dominates the total flexibility correction. Given how good overall performance of SAPT was for the classically allowed region [3], one can safely assume that the effect of higher order SAPT corrections will be marginal. There are two main uncertainties related to the dispersion derivative used in this paper: basis set incompleteness and the time-dependent Hartree-Fock (TDHF) approximation [2,3]. To address the first uncertainty, we performed test calculations using the d-aug-cc-pVQZ basis set for a few geometries around the minimum in T-shape and linear configurations. Using standard basis set extrapolation technique [5] for the dispersion and the exchange-dispersion (note that only these components of the interaction energy should be extrapolated, since the others do not explicitly depend on excitations between subsys-TABLE S2. Second derivatives of the interaction energy components for the He * +H 2 system for the atom-molecule distance of 10.5 bohr. To estimate the basis set limit, we used the d-aug-cc-pVQZ results for E (1) elst , E exch , E (2) ind , E exch−ind , and extrapolated values of E (2) disp and E (2) exch−disp . For the induction and dispersion energies we used the TDHF model [2][3][4], exchange and electrostatic interactions were obtained from the Hartree-Fock densities. The derivatives with respect to r were calculated at r = r c , in units cm −1 /(bohr) 2 .
d-aug-cc-pVTZ estimated basis set limit S5 tems) we found that for the linear geometry the second derivative is underestimated by about 2%, while for T-shape by about 10%. These numbers translate to about 6% underestimation and small 3% overestimation for the V 0 and V 2 parts of the δE flex int correction, respectively. In absolute numbers these values are well below the uncertainty of basis set incompleteness for the E CCSD(T) int part of the total interaction energy. To address the uncertainty of the TDHF method let us note that this model only very slightly overestimates the dispersion energy for the metastable helium dimer by about 3% [2]. Similarly, the comparison of longrange isotropic C 60 coefficient obtained with TDHF [4] (112.8 E h a 6 0 ) and the accurate value of Bishop and Pipin [6] (108.24 E h a 6 0 ) suggests that, indeed, the dispersion energy can be slightly overestimated. Assuming that is proportional to overall performance of the dispersion, one can conclude that the inaccuracy due to the TDHF method is marginal and contributes to less than 0.01 cm −1 in the minimum range.