Vibrational relaxation rate coefficients in CS–He collisions up to the dissociation limit: mixed quantum-classical calculations and neural network predictions

Abstract

Rate coefficients for the process CS(v) + He → CS(v–Δv) + He with Δv = 1, 2, 3 and for v up to 45 were calculated in the 80–5000 K temperature using a mixed quantum-classical (MQC) method. The dataset was then completed using a Neural Network (NN) model, trained on the MQC rate coefficients. The reliability of the MQC method was first verified by comparing the new MQC results with initially state-selected ro-vibrational rate coefficients up to v = 2 computed with the vibrational close-coupling infinite-order sudden (VCC-IOS) method [F. Lique and A. Spielfiedel, Astron. Astrophys., 2007, 462, 1179] employing the same potential energy surface (PES). To enable calculations for higher vibrationally excited states, a new analytical PES was developed that is suitable even under large bond distortions.

---

1 Introduction

The study of inelastic molecular collisions is essential for modelling non-equilibrium physical environments, from plasmas and (exo)planetary atmosphere composition and dynamics to the characterisation of the interstellar medium (ISM).^1–8 Non-local thermal equilibrium (non-LTE) conditions require models that rely on extensive sets of state-to-state collisional data, which are usually impossible to obtain from experiments alone.⁹

Therefore, the theoretical approach is the necessary way to proceed: collisional information is obtained from molecular dynamics simulations using accurate potential energy surfaces (PESs), which describe the interactions between colliding particles. When available, these theoretical results can be compared with experimental data to assess the reliability of the adopted computational methods. A variety of techniques can be used to describe nuclear motion and compute state-to-state collisional cross sections and rate coefficients. Full quantum methods represent the gold standard as they incorporate all possible quantum effects. However, these methods are computationally expensive and their cost increases dramatically with the number of molecular quantum states considered in the calculation, restricting their application to systems where few quantum states are accessible or when the excitation of small and medium sized rigid molecules is induced by light colliders such as He or H₂, unless approximations are introduced. Methods based on classical mechanics, mostly quasi-classical trajectory approaches (QCT), are often the alternative. They can be readily scaled for large systems and are easily parallelisable. Nevertheless, these methods inherently neglect quantum effects and schemes for the assignment of the final quantum states need to be applied, whose accuracy, particularly for vibrational states, is still debated.^10,11

The mixed quantum-classical (MQC) approach^12,13 represents an effective compromise between fully quantum and classical methodologies. In this framework, the degrees of freedom exhibiting the most significant quantum behaviour are treated quantum mechanically, while the remaining ones are treated classically. This enables the inclusion of the quantum effects of interest with a consistent reduction in computational cost. In recent years, the original MQC formalism has been further developed and refined, in applications where rotations (e.g. ref. 14 and references therein) or vibrations (e.g. ref. 15 and references therein) constitute the quantum subset.

In the modelling of ISM or exoplanet atmospheres, there is a growing need to extend the database of available rotationally resolved rate coefficients data for collision-induced transitions between ro-vibrational levels across a large number of species. Indeed, results from full quantum simulations have been obtained for a variety of systems,^16–21 but are often limited to cold temperature conditions, mainly addressing the quenching of rotational states in the ground vibrational state. However, recent years have seen the detection of a growing number of vibrationally excited molecules, making the availability of comprehensive vibrational quenching and excitation rate coefficients an urgent requirement. Indeed, in the past decade, we have undergone a revolution in our ability to observationally study the molecular Universe. The high resolution of new telescopes such as Atacama Large Millimeter/submillimeter array (ALMA) and James Webb Space Telescope (JWST) enables very detailed observations of small astrophysical bodies such as planetary atmospheres, cometary comae, and exoplanets as well as interstellar environments such as molecular clouds and protoplanetary disks. In particular, infrared observations of vibrational molecular signatures, boosted by the MIRI and NIRSPEC instruments on JWST, are flourishing.

Using MQC methods, where the vibrational motions are treated quantum mechanically while rotations and translations are described classically, can be highly effective in such cases. Indeed, this approach has been used to successfully reproduce and extend experimental results on the vibrational quenching of Earth's atmospheric species (N₂, CO, and O₂) through collisions with atoms and diatomic molecules.^22–26

Among diatomic molecules frequently detected in astrophysical environments, carbon monosulfide (CS) is of particular interest, as it has been detected in interstellar molecular clouds, circumstellar envelopes and protoplanetary disks.^27–33 More recently, CS has also been identified in exoplanetary atmospheres as an intermediate species in the photochemical formation of SO₂.³⁴ Molecular collisions involving CS have been studied by Lique and coworkers (CS–He collisions),^35,36 and by different groups (CS–H₂ collisions).^37–40 The results from Lique and coworkers for CS–He collisions are especially relevant, as they provide detailed insights into ro-vibrational excitation processes: first they developed an intermolecular PES to study rotational quenching of CS,³⁵ and then extended their work to include the first two excited vibrational states (this work will be referred as 07LiSp hereafter).³⁶ Rate coefficients for inelastic collisions were calculated using the vibrational close-coupling infinite-order sudden approximation (VCC-IOS) quantum method^41,42 and made available through the BaseCol database.⁴³

In this work we focus on CS–He inelastic collisions. We take on the PES developed in ref. 36, briefly described in Section 3, and apply the MQC method, described in Section 2, to calculate the rate coefficients for the same vibrational transitions for which VCC-IOS results are available. This enables us to directly compare the performance of the MQC and VCC-IOS scattering approaches (Section 5.1). We then develop a new intermolecular analytical PES, based on the Improved Lennard-Jones model⁴⁴ (Section 3.2), which can describe larger deformations of the CS molecule. This new PES can therefore be used to obtain vibration-to-translation/rotation (V–T/R) rate coefficients for highly vibrationally excited states of CS, approaching its dissociation limit. Using the ILJ PES we construct a comprehensive dataset of V–T/R rate coefficients covering all CS vibrational levels up to 45, explicitly calculated by the MQC method and completed using a Neural Network (NN) model (Section 4).

It is worth noting that, while the MQC method can yield state-to-state rate coefficients for specific ro-vibrational transitions (see Section 5.1), the primary objective of this study is to calculate a comprehensive dataset of rotationally averaged rate coefficients for the collisional quenching of vibrationally excited CS by He (see Section 5.2). This dataset includes all rate coefficients for the following processes:

CS(v) + He → CS(v–Δv) + He, with Δv = 1, 2, 3

for v = 1,…, 45 in the temperature interval 80–5000 K.

2 MQC collisional cross sections and rate coefficients

Our MQC implementation for vibrational energy transfer in atom–diatom collisions simultaneously solves the time-dependent Schrödinger equation for the vibrational degree of freedom and ro-vibrational coupling, and the classical equations of motion for the other degrees of freedom (rotation and translation). The total Hamiltonian (H) of the system in Jacobi coordinates (see Fig. 1) can be written as:where r is the diatom bond length, p_r is the molecular linear momentum, m is the molecular mass, v(r) is the vibrational oscillator function, j(t) is the time dependent molecular rotational angular momentum, μ is the reduced mass of the total atom–diatom system, R is the distance between the atom and the center of mass of the diatom, P_R is the system linear momentum, l is the orbital angular momentum, and V(R, r, Θ) is the interaction potential.

Fig. 1 Jacobi coordinates used in this work to describe the CS–He dimer interaction potential. Θ = 0° corresponds to a collinear configuration with the He atom approaching the C atom moiety and Θ = 180° sees the He atom approaching the S atom.

The quantum Hamiltonian is composed by the first three terms of eqn (1). The rotational contribution is present in the quantum and in the classical Hamiltonian, and, together with the effective potential V_eff introduced in the following, defines their interaction. and the remaining terms compose the classical Hamiltonian. The Ehrenfest averaged potential, V_eff(R, Θ), is used instead of V(r, R, Θ) in the classic Hamiltonian and also couples the two Hamiltonians:⁴⁵

V_eff(r, Θ) = 〈Ψ|V(R, r, Θ)|Ψ〉 (2)

where Ψ is the vibrational time-dependent wavefunction expressed in terms of the rotationally distorted vibrational basis ϕ_v:with a_v representing the expansion coefficients for a generic vibrational level v, E_v its energy, and ħ the reduced Planck constant. ϕ_v(r, t) is expressed as:where ϕ⁰_v(r) denotes eigenfunctions of the radial one-dimensional Schrödinger equation based on the rotationless intramolecular potential, and H_v′v is the first order centrifugal stretching term:with r_eq being the diatom equilibrium distance.

The coefficients a_v are obtained by plugging eqn (3) and (4) into the time-dependent Schrödinger equation, and propagating in time the following set of coupled equations together with the classical equations:

The determination of both H_v'v and 〈ϕ_v′|V(R, r, Θ)|ϕ_v〉 (for which V(R, r, Θ) is expanded in a Taylor series around r_eq) requires the calculation of the matrix coupling elements between vibrational states M_v'v:where Δr = r − r_eq. Here we considered terms up to k = 3.

The probability for the vibrational transition v_i → v_f can be obtained as P_vi→vf = |a_vf(∞)|² and depends on the initial conditions (i.e. the initial random values of angle and action) of the particular trajectory considered.

The total energy conservation and micro-reversibility principle imposing that P_vi→vf(E) = P_vf→vi(E) are guaranteed by the utilisation of a symmetrisation procedure¹³ which assigns a total energy E to a specific transition, so that:

with ΔE = E_vf − E_vi and E_cl is the classical energy, i.e. the energy pertaining to the classical degrees of freedom, which in the present case is defined by the sum of the initial translational and rotational energies.

If we consider a specific initial ro-vibrational state (v_i, j_i) we can obtain the probability of the collisional transition to a final vibrational state v_f, by averaging over a number of trajectories N allowing a statistical sampling over the initial conditions:

The corresponding cross sections σ_vi,ji→vf can then be obtained:where J is the total angular momentum, J_max is its maximum value, E − E_vi − E_ji is the kinetic energy E_k. Note that the above cross sections correspond to initially state selected rotational states only, i.e. no resolution for the final rotational quantum number has been considered. However, cross sections resolved for the final rotational states can also be obtained by applying one of the binning schemes commonly used for this purpose in QCT calculations. The rate coefficients can be calculated by averaging over the Boltzmann distribution of kinetic energy as:with β = (k_BT)⁻¹ where k_B is the Boltzmann constant and T the temperature.

In many applications, we are interested in rate coefficients for vibrational transitions averaged over the initial rotational quantum numbers j_i. This approach is especially suitable for heavy molecules such as CS, for which rotational spacings are small and many rotational states contribute to the relaxation process over the temperature range considered here. Under these conditions, statistical effects are expected to reduce the sensitivity of vibrational relaxation rates to the specific initial rotational state. It is therefore important to evaluate the validity of this approximation under the conditions relevant to astrochemical observations. While laboratory experiments have accessed rotational states as high as j = 113,^46,47 astrochemical detections remain limited to much lower excitation. In high-excitation regions such as Orion KL, the highest detected transition of CS is j = 26–25,⁴⁸ and in colder dense cores rotational levels with j < 10 are typically observed. In these environments, the gas density often exceeds the critical density for rotational transitions, ensuring that the rotational population follows a Boltzmann distribution. Nevertheless, in extreme non-LTE or very low-density environments (e.g., cometary comae) state-specific rotational effects may become more important. However, as recently discussed by Bowesman et al.,⁴⁹ approaches assuming rotational LTE and vibrational (or vibronic) non-LTE are expected to capture the dominant non-LTE effects observable with the James Webb Space Telescope, particularly in view of the computational infeasibility of treating the full set of rovibrational transitions explicitly. Furthermore, radiative-transfer models generally require the inclusion of higher-lying rotational levels that may be weakly populated, even if not directly observed, as neglecting them can affect the excitation balance and lead to biased physical conditions.

To determine rotationally averaged rate coefficients, it is convenient to calculate rotationally averaged “cross sections” σ_vi→vf(E_cl, T_ref):

I is the moment of inertia of the diatom. In this formulation a reference arbitrary temperature T_ref = 300 K is used, which cancels out when rate coefficients are computed:¹²with ε = E − E_vi. The above expression is valid for de-excitation collisions.

In the following, we calculated the rate coefficients using eqn (11) to compare the coefficients obtained by the MQC method with those obtained by the VCC-IOS approach,³⁶ and then using eqn (13) to produce the dataset of V–T/R rate coefficients for vibrational quenching in a wide temperature range. In both cases, adequate coverage of the energy spectrum is important for good sampling, which can be achieved by using denser sampling of low energies to describe low temperatures and including a representative number of points at higher energies.

We solved the combined classical and quantum coupled equations using a variable-order, variable-step Adams predictor–corrector integrator, imposing an accuracy of 1 × 10⁻⁸.

3 Potential energy surfaces

3.1 CS intramolecular potential and vibrational wavefunctions

The evaluation of the matrix coupling elements in eqn (7), the vibrational wavefunction (eqn (3)) and the vibrational energy levels requires the use of the rotationless vibrational wavefunctions ϕ^j=0_v(r, t) for the isolated CS molecule.

Initially, to enable comparison with the results of ref. 36, we used the same matrix elements and vibrational energy levels. These were obtained using the vibrational wavefunctions evaluated by the Fourier Grid Hamiltonian (FGH) method,⁵⁰ which is based on a CS potential energy function derived from a Rydberg–Klein–Rees (RKR) curve⁵¹ and using the experimental spectroscopic constants from Bogey et al.⁵² The work reports matrix elements for the ground state and the first two excited vibrational states.

To determine collisional vibrational transition rate coefficients for highly excited states, we employed CS vibrational wavefunctions obtained from the recent spectroscopic model, developed by Paulose et al. (hereafter referred to as JnK).⁵³ This model starts from Coxon & Hajigeorgiou potential energy curve,⁵⁴ with refined parameters that yield improved agreement with experimental data. Experimental vibrational energy levels were used when available through the measured active rotation–vibration energy levels procedure.⁵⁵ The model consists of 50 vibrational states up to the dissociation limit, with rotational states ranging up to j = 257 for v = 0. The energy levels for the first three vibrational states obtained by the two models agree within 0.1 cm⁻¹.

The matrix elements were generated using the vibrational wavefunction from Level16 program,⁵⁶ on a grid spanning 0.75 < r (Å) < 8.0 with a step size of 0.05 Å. The resulting values differ from those reported by Lique and Spielfiedel³⁶ by less than 4%, which is well within the expected 10% accuracy of the method.

3.2 CS–He intermolecular potential: a new analytical PES

The quality of scattering cross sections and rate coefficients critically depends on the quality of the PES describing the interaction potential both at short and long range. In particular, long-range interactions play a fundamental role at low temperatures, where orientation effects—which originate at large intermolecular separations—can significantly influence the dynamics and outcomes of the collisions.

We again start by using the intermolecular PES developed by Lique and Spielfiedel,³⁶ obtained through an analytical fit to ab initio points calculated at the coupled cluster level with single and double excitations and perturbative triples [CCSD(T)], using Dunning's augmented correlation-consistent quadruple-zeta basis set⁵⁷ (aug-cc-pVQZ) further augmented with diffuse [3s3p2d2f1g] bond functions. Basis set superposition error (BSSE) corrections were included following the procedure of Werner et al.⁵⁸ The dependence on the carbon–sulfur bond length r was explicitly incorporated by evaluating the potential at three geometries, 2.4, 2.9, and 3.4 bohr. This makes the potential energy surface reliable for vibrational levels up to v = 2.

Throughout this work, the CS–He potential is described using a Jacobi coordinate system, where r represents the carbon–sulfur bond length, R is the distance between the He atom and the CS centre of mass, and Θ is the angle between these two vectors (Fig. 1).

To address higher vibrationally excited states, we have also developed a new potential energy surface designed to be reliable for larger distortions of the CS bond.

As in previous works,^15,26,59 where we determined large datasets of rate coefficients for vibrational energy transfer, we adopt an analytical model that provides a reliable representation of the interaction, both at equilibrium and for large distortions of the CS bond. This model is expressed as a combination of effective contributions, whose basic parameters have clear physical interpretations in terms of the fundamental properties of the interacting species—such as the static dipole polarizability, the number of outer-shell electrons determining its magnitude, and the permanent electric dipole moment. The numerical values of these parameters were adjusted against accurate ab initio interaction energies for the CS–He ground singlet state computed for a wide range of intermolecular distances using the CCSD(T) approach in the complete basis set (CBS) limit.

In such computations, counterpoise-corrected⁶⁰ interaction energies obtained using the aug-cc-pVQZ and aug-cc-pV5Z⁵⁷ basis sets have been properly combined^61,62 to obtain reliable extrapolated CBS values. To also probe the dependence of the intermolecular potential on the stretching of the CS bond length beyond its equilibrium value r^CS_eq = 2.9 bohr, two further values (3.19 and 3.335 bohr) were considered. Furthermore, ab initio estimations of the CS permanent dipole moment and static dipole polarizability, which are quantities needed for the analytical representation of the global PES (see below), have been obtained as previously done^63–66 for different monomers at the multireference ACPF (Averaged Coupled Pair functional) level in combination with the aug-cc-pV5Z basis set and as a function of the internuclear distance. In particular, the natural molecular orbitals in the ACPF calculations were those from the complete active space self-consistent field reference wave functions. The active space (CAS) considered here is defined by distributing 18 electrons in the orbitals indicated as (2, 3, 4)σ_g (2, 3, 4)σ_u (1, 2)π_u (1, 2)π_g. The 1σ_g 1σ_u core molecular orbitals were constrained to be doubly occupied and excluded from the used CAS while also being fully optimised. Additional estimations of the CS dipole polarizability have been carried out at the CCSD(T)/aug-cc-pV5Z level of theory. All electronic structure calculations were carried out by using the Molpro 2012.1 package.⁶⁷

The complete formulation of the analytical PES is provided in Appendix A, and the corresponding computational subroutine is available in the SI.

Fig. 2 and 3 present several comparisons between the new analytical PES and the ab initio PES derived in ref. 36 for CS at its equilibrium bond length. The contour plot shown in Fig. 2 indicates that the two PESs exhibit a very similar overall topology, with only minor differences in the shape of the potential minimum for certain orientations. The potential cuts shown in Fig. 3 demonstrate excellent agreement between the current ab initio estimations (calculated at the CCSD(T)/CBS level) and the two PESs for the three dimer configurations considered (top three panels). Both PESs accurately reproduce the predicted position and depth of the potential well, with minor discrepancies, below 0.5 meV, occurring for V_anal. The two bottom panels of the same figure illustrate the corresponding angular dependence of the interaction energy at fixed intermolecular distance R values (R = 4 Å and R = 5 Å). In all cases, the quantitative differences in the attractive region are very small (ca. 0.1 meV) and limited to a few specific orientations. This also constitutes a validation of the new analytical PES for small distortions of CS bond length. Indeed, such good agreement between the analytical PES and the ab initio points is remarkable. While this level of agreement is expected for the V07LiSp PES, which was constructed by fitting a large number of ab initio points at the CS equilibrium bond length, the analytical PES relies solely on a limited number of physically motivated parameters. The last condition is basic to formulate a PES reliable in the full range of accessible molecular orientations and in an extended range of molecular deformations.

Fig. 2 Contour plots of the new analytical potential energy surface, V_anal (right panel), and of the V_07LiSp PES³⁶ (left panel), both evaluated at the equilibrium bond length r = r^CS_eq. The minimum of V_anal is located at Θ = 78° and R = 3.90 Å, with a well depth of −2.61 meV, while for V_07LiSp the minimum occurs at Θ = 75°, R = 3.84 Å, and −2.70 meV.

Fig. 3 Potential energy cuts at the equilibrium bond length r = r^CS_eq comparing the results from the V_07LiSp PES,³⁶ the present analytical potential V_anal, and the CCSD(T)/CBS ab initio points.

Direct comparison with the ab initio data computed at the CCSD(T)/CBS level is feasible only for small deviations from the equilibrium bond length (see Fig. S1 in the SI, corresponding to a 15% CS bond elongation). A more extensive comparison would require accounting for the increasing multireference character of the electronic wavefunction at larger distortions, which lies beyond the scope of the present work. The figure shows that the position of the potential minimum is satisfactorily reproduced, with a maximum discrepancy in the well depth of approximately 0.5 meV.

Overall, the accuracy of the analytical PES, based on comparison with ab initio points, is expected to be within 10% for small deformations of the CS bond and within 30% for the largest deformations.

4 The neural network model

In recent years, machine learning (ML) techniques have been increasingly employed to predict cross sections and rate coefficients from a limited set of computed values. The requirement for extensive state-to-state datasets, often spanning wide temperature ranges, makes data-driven prediction from a restricted number of computationally expensive calculations an appealing approach. Consequently, rate coefficients derived from QCT and MQC simulations for both reactive and inelastic scattering have been modelled using Neural Networks (NN) and Gaussian Process Regression (GPR) methods.^25,68–71 Both techniques have demonstrated comparable predictive performance, provided they are applied within the interpolation domain of the training data. However, differences in training efficiency can arise depending on the size of the dataset and the dimensionality of the feature space.⁶⁹

In this work, we construct the dataset directly starting from rate coefficients rather than on cross sections, since the rotationally averaged cross sections (eqn (12)), which depend on the reference temperature, do not represent intrinsic collision properties. In contrast, rate coefficients can be straightforwardly related to physically relevant observables, such as collisional excitation and relaxation rates. Besides, we verified that building the dataset starting from cross sections or from rate coefficients leads to negligible differences in the final results.

Here, an NN model was selected over a GPR one, as a comparative evaluation demonstrated its superior performance on our specific datasets. To ensure a fair comparison, both models first underwent hyperparameter optimization via an extensive grid search. For the GPR model, we optimized parameters including the value of ν (nu) when using the (Matérn kernel). For the NN model, we optimised the architecture by varying the number of dense layers and the number of neurons per layer. The performance of the optimised models was then evaluated using a Leave-One-Out Cross-Validation (LOOCV) strategy. This was applied to the three datasets of rotationally averaged rate coefficients for the CS(v) + He → CS(v–Δv) + He transitions, corresponding to Δv = 1, 2, and 3. In this procedure, the model is trained on all data points in a set except one, which is kept out for testing. This process is repeated for every data point, enabling us to judge the overall effectiveness of each model. The ANN model consistently produced better predictive accuracy in these tests, resulting in its selection as the optimal model for this study.

Following this initial step, the optimal neural network architecture was identified through a two-stage grid search process designed to balance the accuracy of the predictions with the physical plausibility of the resulting output surface. In the first stage, a conventional grid search was performed on a wide range of model hyperparameters, using a standard LOOCV strategy to evaluate performance and select a smaller subset of promising candidates. The second stage involved a more refined evaluation focused on the smoothness of the final interpolated surface. To achieve this, the training data points – consisting here of all the computed rate coefficients values for each Δv – were combined with the model predictions to create a complete surface representation. The smoothness of this surface was then assessed quantitatively by calculating the mean absolute Laplacian. This metric was derived by first interpolating the scattered data points onto a regular grid and then applying a Laplacian kernel via 2D convolution to approximate the surface curvature at each point. Ultimately, models that produced a lower mean absolute Laplacian value were favoured, as this indicates a smoother and less erratic surface. This ensures that the final selection is accurate and generates physically realistic outputs. All codes can be retrieved at ref. 72.

This two-stage optimisation process led to the selection of our final model configuration. The optimal architecture identified was a three-layer network with 256 neurons per layer ([256, 256, 256]). The model was trained using the Adam optimiser with the mean squared error (MSE) loss function. Training was run for 1000 epochs for all datasets, and the model weights from the epoch that achieved the minimum training MSE were selected for the final model. The optimisation process also indicated different optimal batch sizes: a batch size of 50 was used for the Δv = 1 dataset and a batch size of 256 was employed for both the Δv = 2 and Δv = 3 datasets.

5 Results and discussion

5.1 Comparison between MQC and VCC-IOS rate coefficients

The availability of the results for the CS–He ro-vibrational de-excitation cross sections and rate coefficients in ref. 36 gives us the opportunity to compare these quantities using two completely different methods: the time-independent quantum VCC-IOS of ref. 36, and the time-dependent MQC approach of the present work.

The results by Lique and Spielfiedel³⁶ cover ro-vibrational transitions for v = 0, 1, and 2, ranging from j = 0 to j = 37. The reported accuracy is within a factor of 2–3 for the ro-vibrational excitation and within a factor of 1.5 for pure rotational excitation within one vibrational manifold.

Those rate coefficients were compared with the corresponding MQC values calculated using eqn (11) obtained by running 5000 trajectories for each of the 34 initial classical energy values E_cl, in the interval 50–10 000 cm⁻¹, selected on an unevenly distributed energy grid, with a higher density of points at lower energies. These settings should ensure an accuracy of the MQC method within 20%. An initial separation of R = 20 Å was adopted which reflects the coordinate space extension used in Lique and Spielfiedel.³⁶ As mentioned in Section 2, the present implementation of the MQC method does not explicitly select the final rotational state j_f of the CS molecule. However, since the VCC-IOS data include rotational states only up to j = 37, a histogram binning procedure was applied to enable a meaningful comparison. Specifically, the cross sections (eqn (10)) were computed by including only those trajectories leading to a final rotational energy equal to or lower than that corresponding to j_f = 37. Accordingly, the rate coefficients can be compared for the following process:

Calculations were performed for initial rotational states j_i = 0, 1, 2, 5, 10, and 30 with v_i = 1, and for j_i = 0, 5, 10, and 30 with v_i = 2, coupling only the first three vibrational levels. The results for these selected transitions are presented in Fig. 4, while a detailed quantitative comparison of the corresponding rate coefficients is provided in Table S1 of the SI. Fig. S2 in the SI further illustrates the comparison in terms of the ratio between the MQC and VCC-IOS rate coefficients.

Fig. 4 MQC (dashed lines) and VCC-IOS³⁶ (solid lines) rate coefficients for the transition as a function of temperature.

For processes with Δv = 1 (upper panels in Fig. 4), the agreement is excellent for the lowest initial rotational quantum numbers, regardless of the initial vibrational state and across the entire temperature range. For j_i ≤ 5, the maximum deviation remains below 50%, which is well within the combined numerical uncertainties of both methods. The main difference between VCC-IOS and MQC data is that VCC-IOS rate coefficients are basically independent on the initial rotational excitation whereas MQC values show an increase as j_i grows, an effect which decreases with temperature. As a result, MQC and VCC-IOS rate coefficients match very well (i.e. differences are below the methods accuracy) for at least j_i ≤ 10 in the low temperature range and nearly for all the considered j_i's for the highest temperature. This behaviour may be partly attributed to the imposed constraint of j_f ≤ 37. Considering the rotational state distribution within the investigated temperature range, and the known propensity of vibrational relaxation processes to favour transitions with Δj = j_f − j_i = 0, with the probability decreasing rapidly as |Δj| increases, particularly at low temperatures, this limitation is expected to have only a minor effect for the lowest j_i values. In contrast, for the highest j_i levels, where transitions populating high rotational states is more likely, an accurate description of the final rotational distribution becomes increasingly important.

It is worth noting, however, that ro-vibrational excitation rate coefficients calculated using the full close-coupling (CC) method for CO–He collisions⁷³ show a similar trend, with a slight increase as the initial rotational excitation of the molecule rises in the 1000–1500 K temperature range, even though those calculations are limited to j_f ≤ 14. This behaviour suggests that a comparable dependence may also apply to the CS–He system as captured by the MQC-calculated rate coefficients.

The process involving the loss of two vibrational quanta shows a comparable behaviour; however, in this case, the best agreement is observed for j_i > 10 and discrepancies tend to be larger.

Although the level of agreement with the results of Lique and Spielfiedel³⁶ varies with the initial rotational quantum number, the fact that the present work focuses on producing a dataset of rotationally averaged vibrational de-excitation rate coefficients ensures that these differences remain within the expected accuracy of the VCC-IOS results.

5.2 Dataset of vibrational excitation rate coefficients

The rotationally averaged de-excitation rate coefficients for the process CS(v) + He → CS(v–Δv) + He, for vibrational levels up to v = 45, were calculated using eqn (13) and the new analytical PES. The accuracy of the results, considering the combined uncertainties of the PES and of the MQC method, is expected to lie within 25–30% for the transitions involving lower vibrational states and within 50% for the highest ones. Fig. 5 compares the MQC-calculated rate coefficients for the CS(v = 1) + He → CS(v = 0) + He transition obtained with the new V_anal and the 07LiSp PESs. Both surfaces yield similar temperature dependencies, while the new PES systematically predicts values larger by approximately a factor of two. This systematic increase may be related to the slightly earlier onset of the repulsive wall in the analytical PES compared to the V07LiSp PES for collisions occurring within the acceptance angular cone Θ < 90° (Fig, 3), corresponding to He approaching CS from the carbon side (Θ = 0°). Given the larger distance of the carbon atom from the molecular center of mass, such collisions are expected to be more effective in inducing vibrational energy transfer.

Fig. 5 MQC rate coefficients for CS(v = 1) + He → CS(v = 0) + He calculated using the new V_anal (black curve) and the 07LiSp PESs (red curve) as a function of temperature.

MQC de-excitation rate coefficients were calculated for the first three vibrational states and, starting from v_i = 5, every five levels up to v_i = 45. Above this level, the predissociative character of the vibrations caused numerical instabilities in the simulations. For each calculation, all vibrational levels within a 7000 cm⁻¹ interval from v_i were coupled. The temperature range of 80–5000 K was chosen to ensure good coverage of environments ranging from the interstellar medium (ISM) to hot exoplanetary atmospheres. A total of 5000 trajectories were computed for each of the 37 values of E_cl, spanning the 50–20 000 cm⁻¹ interval on an uneven energy grid with a higher point density at low energies. All trajectories were initiated at a separation of R = 50 Å. Vibrational energy levels were taken from the JnK model,⁵³ and matrix coupling elements were obtained from the vibrational wavefunctions following the procedure described by Hong et al.¹⁵

Fig. 6 presents the collisional relaxation rate coefficients corresponding to the loss of one (top panel), two (middle panel), and three (bottom panel) vibrational quanta over the studied temperature range. As expected, in all cases the rate coefficients increase with both the initial vibrational level—owing to the decreasing energy gap—and temperature. The latter dependence is more pronounced at low temperatures and for low vibrational states: the rate coefficients reach a plateau at progressively lower temperatures as v_i increases. Consequently, vibrational relaxation proceeds more rapidly for highly excited vibrational states at low temperatures (as expected from the stronger direct coupling between such states), while becoming nearly independent of v_i at high temperature. This temperature trend becomes increasingly pronounced as Δv grows, and the corresponding rate coefficients decrease in magnitude. The combined effect of these behaviours produces a larger variation in the rate coefficients across the temperature range for the lowest vibrational states, as illustrated in Fig. S3 in SI.

Fig. 6 MQC rate coefficients for CS(v) + He → CS(v–Δv) + He calculated for selected v as a function of temperature.

The dataset was completed by using the NN model as described in Section 4. Panel A of Fig. 7 shows both the MQC calculated and the predicted NN rate coefficients for the CS(v_i) + He → CS(v_i − 1) + He transitions, while panel B shows the corresponding contour plot, illustrating that the NN predicted values smoothly fill in the calculated ones across the full vibrational range. Similar plots for the v_i → v_i − 2 and v_i → v_i − 3 transitions are provided in the SI (Fig. S4 and S5) and exhibit the same overall behaviour. All the rate coefficients are available, together with the NN models, at ref. 74.

Fig. 7 MQC calculated and NN predicted rate coefficients for CS(v_i) + He → CS(v_i − 1) + He. Panel A shows the calculated values (in blue) and the predicted values (in red) as functions of the initial vibrational level and temperature. Panel B presents the corresponding contour plot, displaying both the calculated and predicted rate coefficients.

6 Conclusions and further remarks

In this work, we employed a mixed quantum–classical (MQC) approach combined with a neural network (NN) scheme to construct a complete dataset for the vibrational relaxation of the CS molecule in collisions with He atoms over the 80–5000 K temperature range. Selected vibrational transitions were explicitly calculated using the MQC method, whose performance was preliminarily assessed by comparison with ro-vibrational rate coefficients obtained from the VCC-IOS approach. The agreement between the two methods is generally good, particularly for low initial rotational states. The main difference lies in their dependence on the initial rotational excitation: while VCC-IOS rate coefficients are largely independent of the initial rotational quantum number, the MQC results exhibit a moderate increase with rotational excitation.

A new analytical potential energy surface (PES) was developed using physically meaningful parameters and is expected to remain reliable even for highly excited vibrational states. The MQC-calculated rate coefficients display the expected physical behaviour, increasing with both temperature and initial vibrational level. The NN approach was then applied to extend and complete the MQC dataset, yielding a smooth and physically consistent representation across the full vibrational and temperature ranges.

The resulting MQC + NN dataset provides a coherent description of CS–He collisional vibrational relaxation under a wide range of conditions relevant to astrophysical and (exo)planetary environments. This work highlights the potential of hybrid MQC (or semiclassical) and machine-learning strategies for producing accurate and comprehensive collisional datasets for complex molecular systems.

Future studies will explore the application of this approach to generate large datasets of rotationally state-resolved rate coefficients for vibrational transitions, particularly in regimes where the rotational dynamics exhibit predominantly classical behaviour.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data supporting this article, including the routine for the analytical PES, are available as part of the supplementary information (SI) or can be found in the following repositories: Neural Network codes: CurveFittingML, 2022, https://github.com/lstorchi/CurveFittingML. Trained Neural Network models and rate coefficient datasets: M. Pezzella, F. Pirani, M. Bartolomei, Q. Hong, F. Lique, L. Storchi and C. Coletti, DATA related to: Vibrational relaxation rate coefficients in CS–He collisions up to the dissociation limit: mixed quantum-classical calculations and neural network predictions, 2025, https://doi.org/10.5281/zenodo.17554162. Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp04424c.

Appendix

Appendix A: Formulation of the analytical PES

The analytical semi-empirical potential, V_anal, is made of three components:

V_anal = V_ind + V_ILJ + V_rep. (14)

The first term describes the induction contribution (V_ind) that arises from the interaction between the permanent dipole of CS (μ_CS) and the induced dipole on helium that depends on its polarizability (α_He). This contribution is represented by the following semiempirical formula:⁷⁵

with P₂(Θ) being the second Legendre polynomial. The dipole is analytically described using a polynomial expression of the form:

μ_CS = 0.7758 − 1.45Δr + 0.3Δr² + 0.45Δr³ − 0.17Δr⁴ (16)

where Δr = r − r_eq.

V_ILJ represents the contribution to the interaction arising from the balance between the short-range atom-molecule size repulsion and the long-range dispersion attraction. It has been modelled according to the improved Lennard-Jones (ILJ) formulation.⁴⁴

The V_ILJ contribution is expressed as the sum of pairwise interactions between He with C and S atoms. For each component, it takes the following form:

where R_i describes the distance between the He atom and the C or S atom, which depends on the intermolecular R and the intramolecular r separation distances and on the angle Θ (see Fig. 1). Moreover, ε and R_m represent the depth and the minimum location, respectively, of the pair interaction potential well. These potential parameters have been evaluated from the effective polarizability of each atom within the CS molecule in accordance with the general guidelines reported in ref. 76. The effective atomic polarizability of interest (which is different from that of the isolated atom) is obtained by partitioning the CS molecular polarizability, given in Fig. 8, at each intramolecular distance r, while the polarizability of He is taken from ref. 77. The analytical profile of the CS molecular polarizability, reported in the lower panel of Fig. 8 as a solid line, is expressed in a similar way as previously obtained for the CO molecule (see Appendix A of ref. 66).

Fig. 8 CS dipole moment (upper panel) and molecular polarizability (lower panel) as a function of the intramolecular distance. Solid lines correspond to the analytical expressions (see text).

The term n(R_i) is:

where β_i(r), depending on the hardness of the interacting partners, a property that relates to their polarizability,⁷⁸ is modulated by the CS stretching Δr. Accordingly, the following empirical formula has been adopted to provide reliable β_i(r) values even for large CS deformations, associated with high vibrational levels close to the molecular dissociation:

CS is characterized by a very large equilibrium distance. Its outer electronic charge distribution is highly anisotropic, defining the strength of the repulsion at each orientation angle Θ, and it cannot be represented as a combination of two spheres contributions related to the size of the two atoms. Accordingly, the sum of the two V_ILJ(R_i) components, characterizing previous formulations, has been corrected here with the addition of the term V_rep(r, Θ), which varies significantly with both the intermolecular distance and the molecular orientation:

V_rep(r, Θ) = c_repe^−3R sin⁶(Θ). (20)

For the two interacting pairs, the adopted ε and R_m, evaluated at the equilibrium distance r_eq, and the other potential parameters are given in Table 1. Note that the zeroth order values of ε and R_m, predicted by the atomic polarizabilities of the interacting pairs, have been fine-tuned (R_m by 2–3% and ε by 10–15%) to achieve the best comparison with the ab initio points.

Table 1 Parameters for the analytical potential energy surface at r = r_eq

Parameter	Value
r^CSeq (Å)	1.5346
μr^CSeq (ea0)	0.7758
R^S–Hem (Å)	3.74
R^C–Hem (Å)	3.85
ε^S–He (meV)	2.00
ε^C–He (meV)	1.28
β^eqS–He	8.0
β^eqC–He	8.2
β^1S–He	2.0
β^1C–He	6.0
crep (meV)	70 000

---

The new PES behaves correctly upon deformation of the CS bond thanks to the physically meaningful parameters of its formulation. This is clear from the trend of the analytical expressions for the CS polarizability and dipole moment as a function of r (Fig. 8). These values are in good agreement with the ab initio predictions.

Acknowledgements

L. S., M. P. and C. C. acknowledge support by Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations. National Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) Project ID CN_00000013 “Italian Research Center on High-Performance Computing, Big Data and Quantum Computing” funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di “campioni nazionali di R&S (M4C2-19)” – Next Generation EU (NGEU) and the CINECA award under the ISCRA C initiative, grant IsCc7, for granting the high performance computing resources. M. B. acknowledges the Grant No. PID2024-155666NB-I00 of the Ministerio de Ciencia, Innovación y Universidades (Spain) for funding. Q. H. acknowledges funding from the Chinese Academy of Sciences Project for Young Scientists in Basic Research (Grant No. YSBR-107). F. L. thanks financial support the PNP program (INSU CNRS).

References

1. I. V. Adamovich and W. R. Lempert, Plasma Phys. Controlled Fusion, 2014, 57, 014001.
2. R. Celiberto, I. Armenise, M. Cacciatore, M. Capitelli, F. Esposito, P. Gamallo, R. K. Janev, A. Laganà, V. Laporta, A. Laricchiuta, A. Lombardi, M. Rutigliano, R. Sayós, J. Tennyson and J. M. Wadehra, Plasma Sources Sci. Technol., 2016, 25, 033004.
3. C. D. Pintassilgo and V. Guerra, J. Phys. Chem. C, 2016, 120, 21184–21201.
4. A. V. Pavlov, Geomagn. Aeron., 2011, 51, 143–169.
5. V. Yankovsky and E. Vorobeva, Atmosphere, 2020, 11, 116.
6. M. Keidar, Plasma Sources Sci. Technol., 2015, 24, 033001.
7. E. Roueff and F. Lique, Chem. Rev., 2013, 113, 8906–8938.
8. J. Tennyson, S. N. Yurchenko, J. Zhang, C. A. Bowesman, R. P. Brady, J. Buldyreva, K. L. Chubb, R. R. Gamache, M. N. Gorman, E. R. Guest, C. Hill, K. Kefala, A. Lynas-Gray, T. M. Mellor, L. K. McKemmish, G. B. Mitev, I. I. Mizus, A. Owens, Z. Peng, A. N. Perri, M. Pezzella, O. L. Polyansky, Q. Qu, M. Semenov, O. Smola, A. Solokov, W. Somogyi, A. Upadhyay, S. O. Wright and N. F. Zobov, J. Quant. Spectrosc. Radiat. Transfer, 2024, 326, 109083.
9. M. Capitelli, I. Armenise, E. Bisceglie, D. Bruno, R. Celiberto, G. Colonna, G. D’Ammando, O. De Pascale, F. Esposito and C. Gorse, et al., Plasma Chem. Plasma Process., 2012, 32, 427–450.
10. M. González-Martínez, W. Arbelo-González, J. Rubayo-Soneira, L. Bonnet and J.-C. Rayez, Chem. Phys. Lett., 2008, 463, 65–71.
11. D. S. Lambrecht, G. N. I. Clark, T. Head-Gordon and M. Head-Gordon, J. Phys. Chem. A, 2011, 115, 5928–5935.
12. G. D. Billing, Comput. Phys. Commun., 1984, 32, 45–62.
13. G. D. Billing, Comput. Phys. Rep., 1984, 1, 237–296.
14. B. Mandal, D. Bostan, C. Joy and D. Babikov, Comput. Phys. Commun., 2024, 294, 108938.
15. Q. Hong, L. Storchi, Q. Sun, M. Bartolomei, F. Pirani and C. Coletti, J. Chem. Theory Comput., 2023, 19, 8557–8571.
16. M. Kirste, L. Scharfenberg, J. Kłos, F. Lique, M. H. Alexander, G. Meijer and S. Y. T. van de Meerakker, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 82, 042717.
17. Y. Ajili, K. Hammami, N. E. Jaidane, M. Lanza, Y. N. Kalugina, F. Lique and M. Hochlaf, Phys. Chem. Chem. Phys., 2013, 15, 10062–10070.
18. F. Daniel, M. Gérin, E. Roueff, J. Cernicharo, N. Marcelino, F. Lique, D. C. Lis, D. Teyssier, N. Biver and D. Bockelée-Morvan, Astron. Astrophys., 2013, 560, A3.
19. T. Stoecklin, O. Denis-Alpizar and P. Halvick, Mon. Not. R. Astron. Soc., 2015, 449, 3420–3425.
20. O. Denis-Alpizar and T. Stoecklin, Mon. Not. R. Astron. Soc., 2015, 451, 2986–2990.
21. O. Denis-Alpizar, T. Stoecklin, A. Dutrey and S. Guilloteau, Mon. Not. R. Astron. Soc., 2020, 497, 4276–4281.
22. Q. Hong, M. Bartolomei, C. Coletti, A. Lombardi, Q. Sun and F. Pirani, Molecules, 2021, 26, 7152.
23. Q. Hong, Q. Sun, M. Bartolomei, F. Pirani and C. Coletti, Phys. Chem. Chem. Phys., 2020, 22, 9375–9387.
24. Q. Hong, Q. Sun, F. Pirani, M. A. Valentín-Rodríguez, R. Hernández-Lamoneda, C. Coletti, M. I. Hernández and M. Bartolomei, J. Chem. Phys., 2021, 154, 064304.
25. Q. Hong, L. Storchi, C. Coletti, J. Li, Q. Sun and J. Li, J. Chem. Phys., 2024, 160, 084305.
26. Q. Hong, M. Bartolomei, F. Pirani, Q. Sun and C. Coletti, J. Chem. Phys., 2025, 162, 114308.
27. A. Penzias, P. Solomon, R. Wilson and K. Jefferts, Astrophys. J., 1971, 168, L53.
28. R. N. Martin and A. H. Barrett, Astrophys. J. Lett., 1975, 202, L83–L86.
29. C. Henkel and J. Bally, Astron. Astrophys., 1985, 150, L25–L27.
30. K. Drdla, G. R. Knapp and E. F. van Dishoeck, Astrophys. J., 1989, 345, 815.
31. C. Sánchez Contreras, V. Bujarrabal and J. Alcolea, Astron. Astrophys., 1997, 327, 689–698.
32. M. Agúndez, J. P. Fonfría, J. Cernicharo, C. Kahane, F. Daniel and M. Guélin, Astron. Astrophys., 2012, 543, A48.
33. R. Le Gal, K. I. Öberg, R. A. Loomis, J. Pegues and J. B. Bergner, Astrophys. J., 2019, 876, 72.
34. S.-M. Tsai, E. K. Lee, D. Powell, P. Gao, X. Zhang, J. Moses, E. Hébrard, O. Venot, V. Parmentier and S. Jordan, et al., Nature, 2023, 617, 483–487.
35. F. Lique, A. Spielfiedel and J. Cernicharo, Astron. Astrophys., 2006, 451, 1125–1132.
36. F. Lique and A. Spielfiedel, Astron. Astrophys., 2007, 462, 1179–1185.
37. O. Denis-Alpizar, T. Stoecklin, P. Halvick, M.-L. Dubernet and S. Marinakis, J. Chem. Phys., 2012, 137, 234301.
38. O. Denis-Alpizar, T. Stoecklin, P. Halvick and M.-L. Dubernet, J. Chem. Phys., 2013, 139, 204304.
39. O. Denis-Alpizar, T. Stoecklin, S. Guilloteau and A. Dutrey, Mon. Not. R. Astron. Soc., 2018, 478, 1811–1817.
40. B. Yang, P. Zhang, C. Qu, P. C. Stancil, J. M. Bowman, N. Balakrishnan and R. C. Forrey, Phys. Chem. Chem. Phys., 2018, 20, 28425–28434.
41. R. Goldflam, D. J. Kouri and S. Green, J. Chem. Phys., 1977, 67, 5661–5675.
42. R. Goldflam, S. Green and D. Kouri, J. Chem. Phys., 1977, 67, 4149–4161.
43. M. L. Dubernet, C. Boursier, O. Denis-Alpizar, Y. A. Ba, N. Moreau, C. M. Zwölf, M. A. Amor, D. Babikov, N. Balakrishnan, C. Balança, M. Ben Khalifa, A. Bergeat, C. T. Bop, L. Cabrera-González, C. Cárdenas, A. Chefai, P. J. Dagdigian, F. Dayou, S. Demes, B. Desrousseaux, F. Dumouchel, A. Faure, R. C. Forrey, J. Franz, R. M. García-Vázquez, F. Gianturco, A. Godard Palluet, L. González-Sánchez, G. C. Groenenboom, P. Halvick, K. Hammami, F. Khadri, Y. Kalugina, I. Kleiner, J. Klos, F. Lique, J. Loreau, B. Mandal, B. Mant, S. Marinakis, D. Ndaw, P. Pirlot Jankowiak, T. Price, E. Quintas-Sánchez, R. Ramachandran, E. Sahnoun, C. Santander, P. C. Stancil, T. Stoecklin, J. Tennyson, F. Tonolo, R. Urzúa-Leiva, B. Yang, E. Yurtsever and M. Zóltowski, Astron. Astrophys., 2024, 683, A40.
44. F. Pirani, S. Brizi, L. Roncaratti, P. Casavecchia, D. Cappelletti and F. Vecchiocattivi, Phys. Chem. Chem. Phys., 2008, 10, 5489–5503.
45. P. Ehrenfest, Z. Phys., 1927, 45, 455–457.
46. R. Ram, P. Bernath and S. Davis, J. Mol. Spectrosc., 1995, 173, 146–157.
47. H. Uehara, K. Horiai and Y. Sakamoto, J. Mol. Spectrosc., 2015, 313, 19–39.
48. N. R. Crockett, E. A. Bergin, J. L. Neill, C. Favre, P. Schilke, D. C. Lis, T. A. Bell, G. Blake, J. Cernicharo, M. Emprechtinger, G. B. Esplugues, H. Gupta, M. Kleshcheva, S. Lord, N. Marcelino, B. A. McGuire, J. Pearson, T. G. Phillips, R. Plume, F. van der Tak, B. Tercero and S. Yu, Astrophys. J., 2014, 787, 112.
49. C. A. Bowesman, S. N. Yurchenko, A. Al-Refaie and J. Tennyson, TIRAMISU: Non-LTE radiative transfer for molecules in exoplanet atmospheres, 2025, https://arxiv.org/abs/2508.12873.
50. C. C. Marston and G. G. Balint-Kurti, J. Chem. Phys., 1989, 91, 3571–3576.
51. R. J. Le Roy, J. Quant. Spectrosc. Radiat. Transfer, 2017, 186, 158–166.
52. M. Bogey, C. Demuynck and J. Destombes, J. Mol. Spectrosc., 1982, 95, 35–42.
53. G. Paulose, E. J. Barton, S. N. Yurchenko and J. Tennyson, Mon. Not. R. Astron. Soc., 2015, 454, 1931–1939.
54. J. A. Coxon and P. G. Hajigeorgiou, Chem. Phys., 1992, 167, 327–340.
55. T. Furtenbacher, A. G. Császár and J. Tennyson, J. Mol. Spectrosc., 2007, 245, 115–125.
56. R. J. Le Roy, J. Quant. Spectrosc. Radiat. Transfer, 2017, 186, 167–178.
57. R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
58. H.-J. Werner, B. Follmeg and M. H. Alexander, J. Chem. Phys., 1988, 89, 3139–3151.
59. J. Yang, Q. Hong, M. Bartolomei, C. Coletti, F. Pirani, Q. Sun and J. Li, Phys. Rev. A, 2025, 111, 032804.
60. S. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
61. A. Halkier, T. Helgaker, P. Jorgensen, W. Klopper, H. Koch, J. Olsen and A. K. Wilson, Chem. Phys. Lett., 1998, 286, 243–252.
62. A. Halkier, T. Helgaker, P. Jorgensen, W. Klopper, J. Olsen and A. K. Wilson, Chem. Phys. Lett., 1999, 302, 437–446.
63. M. Bartolomei, E. Carmona-Novillo, M. I. Hernández, J. Campos-Martínez and R. Hernández-Lamoneda, J. Comput. Chem., 2011, 32, 279–290.
64. M. Bartolomei, F. Pirani, A. Laganà and A. Lombardi, J. Comput. Chem., 2012, 33, 1806–1819.
65. A. Lombardi, F. Pirani, A. Laganà and M. Bartolomei, J. Comput. Chem., 2016, 37, 1463–1475.
66. A. Lombardi, F. Pirani, M. Bartolomei, C. Coletti and A. Laganà, Front. Chem., 2019, 7, 309.
67. H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby, M. Schütz, P. Celani, T. Korona, G. Rauhut, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, C. Hampel, G. Hetzer, A. W. Lloyd, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, P. Palmieri, R. Pitzer, U. Schumann, H. Stoll, A. J. Stone, R. Tarroni and T. Thorsteinsson, MOLPRO, Version 2012.1, a Package of Ab Initio Programs, 2012, see https://www.molpro.net.
68. J. Arnold, J. C. San Vicente Veliz, D. Koner, N. Singh, R. J. Bemish and M. Meuwly, J. Chem. Phys., 2022, 156, 034301.
69. Q. Hong, L. Storchi, M. Bartolomei, F. Pirani, Q. Sun and C. Coletti, Eur. Phys. J. D, 2023, 77, 128.
70. J. Yang, J. Li, J. Li and J. Li, J. Phys. Chem. A, 2024.
71. D. Bossion, G. Nyman and Y. Scribano, Artif. Intell. Chem., 2024, 2, 100052.
72. L. Storchi, CurveFittingML, 2022, https://github.com/lstorchi/CurveFittingML.
73. D. R. Flower, Mon. Not. R. Astron. Soc., 2012, 425, 1350–1356.
74. M. Pezzella, F. Pirani, M. Bartolomei, Q. Hong, F. Lique, L. Storchi and C. Coletti, DATA related to: Vibrational relaxation rate coefficients in CS–He collisions up to the dissociation limit: mixed quantum-classical calculations and neural network predictions, 2025,  DOI:10.5281/zenodo.17554162.
75. G. C. Maitland, M. Rigby, E. B. Smith and W. A. Wakeham, Intermolecular Forces: Their Origin and Determination, Clarendon Press, Oxford, 1987.
76. R. Cambi, D. Cappelletti, G. Liuti and F. Pirani, J. Chem. Phys., 1991, 95, 1852–1861.
77. T. N. Olney, N. Cann, G. Cooper and C. Brion, Chem. Phys., 1997, 223, 59–98.
78. M. Capitelli, D. Cappelletti, G. Colonna, C. Gorse, A. Laricchiuta, G. Liuti, S. Longo and F. Pirani, Chem. Phys., 2007, 338, 62–68.

---