A Spectroscopic Validation of the Improved Lennard–Jones Model

The Lennard–Jones (LJ) and Improved Lennard–Jones (ILJ) potential models have been deeply tested on the most accurate CCSD(T)/CBS electronic energies calculated for some weakly bound prototype systems. These results are important to plan the correct application of such models to systems at increasing complexity. CCSD(T)/CBS ground state electronic energies were determined for 21 diatomic systems composed by the combination of the noble gas atoms. These potentials were employed to calculate the rovibrational spectroscopic constants, and the results show that for 20 of the 21 pairs the ILJ predictions agree more effectively with the experimental data than those of the LJ model. The CCSD(T)/CBS energies were also used to determine the β parameter of the ILJ form, related to the softness/hardness of the interacting partners and controlling the shape of the potential well. This information supports the experimental finding that suggests the adoption of β≈9 for most of the systems involving noble gas atoms. The He-Ne and He-Ar molecules have a lifetime of less than 1ps in the 200–500 K temperature range, indicating that they are not considered stable under thermal conditions of gaseous bulks. Furthermore, the controversy concerning the presence of a “virtual” or a “real” vibrational state in the He2 molecule is discussed.


Introduction
The detailed characterization of several equilibrium and non-equilibrium properties of matter (in condensed and gaseous phases) is often obtained through the proper formulation of force fields associated with non-covalent intermolecular interactions [1]. The adoption of simple and accurate models of this type of interactions, to be easily used in molecular dynamics simulations of both ionic and neutral aggregates, still represents a basic question. In particular, such models must be given in the analytical form, from which the first and second derivative of the interaction, defining force, and force constant must be easily obtained and must present continuity of behavior. Moreover, they must involve few parameters having a defined physical meaning that can be used as proper scaling factors when the extension to systems at increasing complexity is attempted. This target can be achieved by investigating in detail prototype systems for which accurate experimental and theoretical information on the intermolecular interaction is easily obtainable.
Some years ago the adoption of an Improved Lennard-Jones (ILJ) function [20] permitted to obtain, for many noble gas pairs, the most accurate experimental values of well depth D e and equilibrium distance R e from the combined analysis of scattering experiments, with the resolution of fundamental quantum interference effects, spectroscopy, and transport properties. In the same paper, it has been also demonstrated that while ILJ provides asymptotically a dipole-dipole dispersion attraction coefficient equal to C 6 = D e · R 6 e , which is in good agreement with the most accurate theoretical and experimental values, the LJ model predicts C 6 = 2D e · R 6 e a factor 2 larger, with also a poor reproduction of the experimental observables. Moreover, values in the range of 7 to 9 (depending on the softness of the interacting partners) of the additional parameter β in ILJ formulation (See next section) work well for several neutral-neutral and ion-neutral cases [20,21]. Such values allowed a proper assessment of the role of the van der Waals interaction component in the formation of the weak hydrogen and halogen intermolecular bonds [21,22]. The excessive long-range attraction of LJ can be a strong limitation when the model is applied, as often made, to describe the behavior of big molecules, where many interaction centers are involved, several of them separated by large distances. For the application of ILJ function to systems at increasing complexity, like those involving biomolecules, the selective passage of chemical species in cellular channels and pores, the physical adsorption on single and multiple layers, further tests with a possible generalization of its formulation are desirable and probably necessary. The achievement of this target can be pursued through a sequence of steps.
In this paper, we test in detail the shape of the potential well predicted by ILJ and LJ on accurate ab initio values of the interaction and the combined analysis of spectroscopic features of both symmetric and asymmetric noble gas dimers. This study confirms that ILJ with β ≈ 9 provides the best representation for most of the investigated systems. The following steps should involve an accurate analysis of short-range repulsion of atom-atom systems with great difference in the polarizability (or in softness) to test the modulation of the repulsion by varying the β value. Moreover, other important information can be provided by a further accurate study of systems, formed by neutral partners interacting with both negative and single or multiple charged positive ions. For such systems, the depth, location, and shape of the potential well and the steepness of the first part of the repulsion arise from a more critical balance of the attraction, stronger respect to the neutral-neutral cases, and of the repulsion.
This work is organized as follows. The methodologies adopted in this study are summarized in Section 2, while obtained results and discussion are presented in Section 3, and, finally, some conclusions are provided in Section 4.
To correctly describe the electronic structure of heavier elements, it is necessary to include the relativistic effects [33][34][35]. As examples, (i) relativistic effects account for 1.7-1.8 V in a standard 2 V lead-acid battery cell [36], and (ii) the non-relativistic gold is white (like silver) so the yellow color of gold comes from relativity [34]. A typical way of including relativity in electronic structure calculations is through the use of pseudopotentials [37,38] (PP). The relativistic pseudopotentials used in this investigation for Xe and Rn were the small core energy consistent developed by Peterson et al. [39], which were adjusted to multiconfiguration Dirac-Hartree-Fock data based on Dirac-Coulomb-Breit Hamiltonian, with the accompanying aug-cc-pVQZ-PP and aug-cc-pV5Z-PP basis sets. Note that errors due to the used pseudopotentials are found negligible [39] (they are expected to provide a maximum contribution to D e of about 1.3 KJ/mol). In addition, the selected PP with the matched basis sets exhibits the systematic convergence and accuracy characteristic of their all-electron counterparts also used in this investigation for lighter He, Ne, Ar, and Kr [28,29]. All these PECs were calculated through the Gaussina09 computational code [40].
Through the R e and D e calculated values, combined with an accurate investigation of the radial dependence of the interaction, effective isotropic PECs have been constructed exploiting the LJ and ILJ analytical forms. The general formulation of the classical LJ model is given by the following equation: that for neutral-neutral systems, with n = 12 and m = 6, this equation turns into the following well-known form: For the ILJ function it has been proposed that where n(R) = β + 4 R R e 2 and β parameter describes the softness/hardness of the elements involved in the complex and β is experimentally set to 9 for systems involving noble gases [20]. For neutral-neutral systems, m assumes the value of 6 and the ILJ form becomes Rovibrational energies of each diatomic molecule were determined by solving the nuclear Schrödinger equation. To solve this equation, the Discrete Variable Representation method [41] was employed. Rovibrational spectroscopic constants, such as ω e , ω e x e , ω e y e , α e , and γ e , were calculated using the following expressions [42]: In Equation (5), E v,j represents the rovibrational energy, where the indices v and j indicate the vibrational and rotational quantum numbers, respectively. To verify the accuracy of spectroscopic constants, the Dunham method [43] was also used. This approach depends on the derivatives of PECs in the equilibrium configuration.
For each Ng-Ng molecules, the lifetime as a function of temperature was calculated using Slater's method which is described by the equation [44,45]: In the Equation (6), T is the temperature, R the universal gas constant, and E 0,0 the zero-point energy. This equation provides the lifetime for decomposition of the systems and it is a description purely dynamical with a vibrational analysis of the complexes, referring to the low or high rate of unimolecular decay and it is supposed to occur when the interaction coordinate reaches the dissociation threshold (D e ). In general, this approach is suitable for regions of intermediate pressure in the bulk. Table 1 shows the CCSD(T)/aug-cc-pVQZ, CCSD(T)/aug-cc-pV5Z, CCSD(T)/CBS, and experimental equilibrium distances for all Ng-Ng diatomic molecules (with Ng = He, Ne, Ar, Kr, Xe, and Rn). From this table it is possible to note that the equilibrium distances calculated with CCSD(T)/CBS level agree more effectively with experimental data [20,46]. CCSD(T)/aug-cc-pVQZ, CCSD(T)/aug-cc-pV5Z, CCSD(T)/CBS results and experimental dissociation energies, for the Ng-Ng molecules, are compared in Table 2. These results also show that the best agreement between theoretical and experimental data happens with the CCSD(T)/CBS level, mainly when compared with the data available in [20,47]. Table 1. CCSD(T)/aug-cc-pVQZ, CCSD(T)/aug-cc-pV5Z, CCSD(T)/CBS, and experimental equilibrium distances (Å) for the Ng-Ng molecules (Ng = He, Ne, Ar, Kr, Xe, and Rn).  The twenty-one complete PECs for all systems with CCSD(T)/aug-cc-pVQZ, CCSD(T) /aug-cc-pV5Z, and CCSD(T)/CBS levels are shown in the supplementary material (from Tables S1-S21). These PECs were built calculating the ground state electronic energies for different values of the internuclear distances (R) that ranged from the region of the strong interaction (R less than R e ) to the asymptotic region (R much larger than R e ). For R less and greater than equilibrium distance (R e ), it was used a step of 0.1 Å, while for R near to R e was considered a step of 0.01 Å. With these steps, it was determined approximately a hundred electronic energies for all Ng-Ng molecules (except for the Kr-Rn, Xe-Rn, and Rn 2 systems).

Molecules aug-cc-pVQZ aug-cc-pV5Z CBS
The β parameter of ILJ model (Equation (4)) was determined for each molecule by fitting, via Powell method [48], the set of CCSD(T)/CBS electronic energies as shown in Table 3. From this table, it is possible to note that β parameter, which describes the softness/hardness of the elements involved in the molecule, for each molecule is very close to the experimental value. This fact supports the experimental prediction that this parameter is close to 9 for most molecules formed with noble gases. The root means square deviation of the performed fitting varied from 3.97 × 10 −5 Hartree (for Ar-Kr system) to 1.05 × 10 −7 Hartree (for He 2 system) for all considered molecules (see Table S22 of Supplementary Information). Table 3. β parameter values adjusted using the CCSD(T)/aug-cc-pVQZ, CCSD(T)/aug-cc-pV5Z, and CCSD(T)/CBS electronic energies for the Ng-Ng molecules (Ng = He, Ne, Ar, Kr, Xe, and Rn).

Rovibrational Energies, Spectroscopic Constants, and Lifetime
Once the CCSD(T)/CBS ILJ PEC of the 21 studied molecules were obtained, their rovibrational energies were calculated using the reduced mass showed in Table S23 of Supplementary Material and they can be found in Tables S24 and S25 of Supplementary Material. The experimental vibrational energy spacings for the Ne 2 (1 transition), Ar 2 (5 transitions), Kr 2 (9 transitions), and Xe 2 (10 transitions) systems [20] were compared with the present results. From this comparison, it was found a difference of 0.38 cm −1 for Ne 2 (1-0 transition) and a maximum and minimum difference of 0.47 cm −1 (1-0 transition) and 0.09 cm −1 (2-1 transition) for Ar 2 , 0.60 cm −1 (1-0 transition) and 0.00 cm −1 (8-7 transition) for Kr 2 , and 0.27 cm −1 (3-2 transition) and 0.01 cm −1 (8-7 transition) for Xe 2 , respectively. Furthermore, from the point of view of the CCSD(T)/CBS calculation (D e = 0.914 meV, R e = 2.97Å, and β = 8.68), the He 2 system does not present a vibrational level within the PEC, i.e., the He 2 first rovibrational state is considered virtual (unbound). This fact agrees with that found by Wang et al. [49], but disagrees with the calculations conducted by Aziz et al. [50] (based on the LM2M2 semiempirical potential) and Tang et al. [51] that predict one weakly bound state for 4 He 2 dimer. This controversy is expected because 4 He 2 dimer interactions are composed by the combination of small mass and small atomic polarizability and it makes that the rovibrational energy of the lowest state places very close to that of the separated atoms. In particular, the small potential well arises from the critical combination of a limited repulsion with the weakest attraction existing in nature and it controls many peculiarities of He in gaseous [52] and condensed phases, as its anomalous phase diagram. This issue was resolved experimentally by Luo et al. [53,54] with the mass spectrometric observation of bound 4 He 2 in an extreme pulsed supersonic beam of He at temperatures less than 1mK. This fact was confirmed by Schöllkopf et al. [55] by using a novel diffraction experiment employing a nanoscale transmission grating. Finally, the current results suggest that, in addition to the He 2 system, also the He-Ne, He-Ar, He-Kr, and He-Rn molecules are less stable, as they only showed two vibrational levels within their PECs.
Moreover, to emphasize the ILJ sensitivity to the potential parameters and to cast further light on the controversy concerning the presence of a "virtual" or a "real" rovibrational state in He 2 molecule, confined at the dissociation limit of a very small potential well, we modulated slightly shape and depth of the potential well of the ILJ potential formulation. In particular, the shape has been adjusted by lowering β (maintaining always its value within the limit 7-9 typical of van der Waals forces for neutral-neutral systems [20,21] and accompanying it by a maximum D e increase of 0.1 meV (0.01 KJ/mol) respect to the CCSD(T)/CBS result, to include in this increase any possible uncertainty of ab initio calculations. The new ILJ formulation, adopting β = 7.6, D e = 0.988 meV (7.9687 cm −1 ) and R e = 2.974 Å, provides results still consistent with the experimental determination [46] and with those of other more sophisticated potential models [51,52] in an extended range of 2.0 to 6.0 Å of internuclear distances, and its well contains here a "real" vibrational state with a value of 7.9685 cm −1 .
To verify the influence of the R e , D e , and β parameters on the quality of the ILJ PEC, the rovibrational spectroscopic constants (RSC) of the 21 molecules were calculated considering these parameters determined at CCSD(T)/CBS level and with an experimental value of β. Tables 4-6 show the RSC calculated by using both Equation (5) (whose rovibrational energies were calculated using the DVR method) and Dunham method. It is important to mention that Equation (5) can only be used for systems that have at least four vibrational levels within the PEC well. Thus, the RSC for the He-Ne, He-Ar, He-Kr, He-Xe, He-Rn, and Ne-Ne molecules were only determined by the Dunham method. From these Tables, note that the RSC determined with ILJ model agrees more with the experimental data than LJ representation for twenty of the twenty-one studied molecules (except the He-Ne molecule). For almost all of the 21 studied molecules (except for He-Ar, He-Kr, He-Ne, Ne-Xe, and Ar-Kr systems), the RSC agrees more with experimental data when an ILJ PEC with β = 9 (experimental value) is used. This fact suggests that β = 9 is an accurate choice to describe molecular systems involving noble gases.
To specify each type of calculation, the following nomenclatures were used: D-ILJ- . For the He-Xe system, the RSC were calculated using the β parameter obtained at aug-cc-pV5Z basis set because the β adjustment for the CBS base did not converge and it is indicated in the Table 5 by the D-ILJ-5z-βFIT symbol. Figure 1 shows the lifetime as a function of temperature for all studied molecules, except for the He 2 system that has no vibrational level within its PEC well. From this figure, one can see that He-Ne and He-Ar molecules have a lifetime of over 1.0 picosecond for all considered temperature ranges (200-500 K) and that the He-Kr lifetime is slightly larger than 1 picosecond within the same temperature range. Following the recommendations of wolfgang [56], which states that a lifetime over 1.0 picosecond means that the PEC well is not deep enough to exclude the intermediate complex, these compounds can be considered unstable. It is not possible to determine the He 2 lifetime, because no vibrational level or only one bound state at the dissociation limit was found within the He 2 PEC.

Molecules
Methods ω e ω e x e ω e y e α e γ e Exp.

Conclusions
In this work, an accurate test of the β parameter value, defining the strength of both attraction and repulsion in the ILJ model (See Equations (3) and (4)), was obtained exploiting CCSD(T)/CBS electronic energies calculated for the complete family of the diatomic molecules formed by the He, Ne, Ar, Kr, Xe, and Rn noble gas atoms. For all considered molecules, β was found to be close to 9 and this feature is supported by both theoretical and experimental findings. To verify the influence of the shape of the potential well on the observables, ILJ PEC has been adopted to predict rovibrational energies, spectroscopic constants, and lifetime as a function of temperatures. The results suggest that ILJ analytical form with β ≈ 9 provides rovibrational spectroscopic constants (RSC) that agree more effectively with experimental than RSC determined with LJ PEC. This fact confirms that most of the LJ inadequacies at large and short intermolecular distances are overcome by the ILJ model. Predicted lifetimes indicate that the He-Ne and He-Ar molecules are not stable under temperature confined in the 200 to 500 K range.
We found that an increase in He 2 well depth of less than 0.1 meV (0.01 KJ/mol) accompanied by a slight change in its shape and position of the well (some fraction of a hundredth of Angstrom), all the characteristics that arise from a very critical balance of weak attraction with the repulsion, leads to the existence of a real vibrational level. Note that these changes are within the errors of any ab initio calculation, even of the CCSD(T)/CBS type that extends to long-range asymptotic regions.
Finally, as an important conclusion, further investigations, carried out combining accurate theoretical and experimental information and focused on the critical balance of attraction and repulsion controlled by the n(R) term, are expected to be crucial for the correct modulation of β parameter and of the numerical coefficient 4 (See Equations (3) and (4)) when systems with completely different nature and size are involved in noncovalent interactions.