The intermolecular potential energy, V_inter, of the two interacting molecules, CO and N₂, is formulated as a combination of two effective interaction components:

where V_vdW and V_elect represent the van der Waals (size repulsion plus dispersion-attraction) and the electrostatic interaction components, respectively. V_elect originates from the anisotropic molecular charge distributions of the two bodies, which asymptotically tend to the sum of the (permanent) quadrupole-(permanent) quadrupole and dipole-quadrupole interactions. Both V_vdW and V_elect depend on the distance R between the centers of mass of the two molecules (say a for CO and b for N₂), and on the Jacobi angular coordinates Θ_a, Θ_b and Φ describing the a − b mutual orientation as well (see Figure 1). In the present work we will consider a series of limiting configurations of the interacting molecules, specifically (Θ_a,Θ_b,Φ)=(90°, 90°, 0°), (90°, 90°, 90°), (90°, 0°, 0°), (0°, 90°, 0°), (180°, 90°, 0°), (0°, 0°, 0°) and (180°, 0°, 0°), which will be referred to as H, X, T_a, T_b1, T_b2, I₁ and I₂.

The van der Waals term, V_vdW of Equation (1), is expressed as a sum of the noncovalent contributions VvdWi as follows,

where r_i is the distance between different interaction centers of the involved molecules, which here are chosen to coincide with the constituting atoms (see Figure 1). Accordingly, the sum of Equation (2) runs over all four atom pairs of the CO-N₂ complex (see Figure 1). It must be emphasized here that in the present work the bond-bond formulation is replaced by the pseudoatom-pseudoatom one in order to properly account for the variation of the anisotropy with the atom-atom distance and the polarizability consistent with the diatom (Pirani et al., 2019).

The explicit form of VvdWi is obtained using the Improved Lennard-Jones (ILJ) potential function (Pirani et al., 2008):

where x_i is the reduced distance defined as

and ε and R_m are, respectively, the well depth and position of the corresponding pair interaction. Note that the ILJ function (Pirani et al., 2008) is more realistic than the original Lennard-Jones (12,6) one, with a much more accurate size repulsion (first term within the square brackets) and the long range dispersion attraction tail (second term within the square brackets) (Pirani et al., 2004; Lombardi and Palazzetti, 2008).

The exponent n in Equation (3) is expressed as a function of x_i using the following empirical equation (Pirani et al., 2004):

in which β is a parameter depending on the nature and the hardness of the interacting particles. For the present system, β has been set equal to 8 (a value typical of van der Waals interactions in neutral-neutral systems) for all atom-atom pairs.

The relevant ε and R_m parameters are directly related to the molecular polarizability of the involved partners, and a zeroth-order estimate of them can be obtained from correlation formulas (Pirani et al., 2001, 2019; Cappelletti et al., 2008) which exploit “effective” atomic polarizability.

In this way a tentative full dimensional PES is generated and related ε and R_m parameters can be fine-tuned by fitting experimental data and by comparing model predictions with accurate ab initio electronic structure calculations (see below).

The V_elect term of Equation (1) has been formulated as a sum of Coulomb potentials as follows:

with q_ja and q_kb being point charges (located on the monomers a and b, respectively, and having values consistent with the corresponding calculated molecular dipole and quadrupoles) and r_jk being the distance between them.

Such a formulation of V_elect must be used for cases in which the molecular dimensions are not negligible with respect to the intermolecular distance R (Maitland et al., 1987). For both monomers a linear distribution of charges as that used in a previous work, see Lombardi et al. (2016b), has been adopted, which consists of three charges placed on the atoms and on the molecule center of mass.

For charge values q_i (as well as for the corresponding position r_i), we took for N₂ the values reported in Lombardi et al. (2016b) while for CO we obtained them by exploiting corresponding dipole and quadrupole moment calculations (see below) together with simple geometrical considerations (see Appendix A). A comparison of the van der Waals and electrostatic interaction contributions to the intermolecular potential energy, for the various configurations of the two interacting molecules is reported in Figure S1 of the Supplementary Material.

The above-mentioned adopted values of ε and R_m (hereinafter called “predicted”) were fine-tuned by carrying out a comparison with ab-initio estimates of the intermolecular part of the interaction and an analysis of the second virial coefficient data (see next section). Predicted and optimized values for the case of rigid monomers are given in Table 1 together with the other parameters used to formulate the intermolecular potential.

In Figure 2 the main features of the intermolecular potential V_inter calculated for seven selected geometries of the interacting system at the CCSD(T) level of theory (see next section) are reported. Although such results will be discussed in more detail later, it is worth pointing out here that, as shown by the figure, the stability ranking of the investigated geometries on the empirical PES model globally agrees well with that of the ab initio values.

In order to carry out the fine tuning of the parameters of the semiempirical functional representation of the PES, we followed the procedure already employed in Bartolomei et al. (2012); Lombardi et al. (2016b), and we first investigated the dimer formation for two rigid monomers (with bond lengths and angles fixed at their equilibrium values). In this perspective, we performed CCSD(T) intermolecular potential calculations for the H, X, T_a, T_b1, T_b2, I₁, and I₂ configurations of the dimer (see section 2.1), as a function of the distance R of the two monomers. Results are reported in Figure 2 where they are compared with the corresponding empirical PES energy profiles. The supermolecular CCSD(T) energies were calculated using the MOLPRO package (Werner et al., 2006) and corrected using the counterpoise method (Boys and Bernardi, 1970; van Lenthe et al., 1987) in order to remove the basis set superposition error. For all supermolecular calculations, the Dunning's aug-cc-pVTZ basis set (Kendall et al., 1992) was used together with the bond function set [3s3p2d1f] developed by Tao Tao and Pan (1992) and placed on the midpoint of the intermolecular distance R that was varied in the range 2.5–8 Å. In order to assess the convergence of the supermolecular energies with the basis set, additional calculations using the aug-cc-pVQZ plus bond functions (see above) basis set were also performed for the usual dimer configurations around their minima. The corresponding results, compared in Table 2 with those obtained using the less extended basis set, show deviations of about 0.1–0.2 meV. The C-O and N₂ bond lengths of the linear monomers were set equal to 1.1283 Å (Linstrom and Mallard, 2018) and 1.1007 Å (Bartolomei et al., 2011), respectively.

In order to carry out a comparison with the flexible monomer empirical potential (to be described below), further CCSD(T) calculations were carried out by considering a maximum variation of 10% in the bond length of one of the two monomers. Such small deformations did not alter the single determinant character (Lee and Taylor, 1989) of the used wavefunctions and allowed the use of the CCSD(T) method also for the deformed arrangements. The ab initio interaction profiles related to flexible monomers, reported in the following, have been compared with the results of the model potential in Figures 6, 7, to prove the validity of the model predictions.

In order to introduce into the empirical model potential the dependence of the electrostatic component V_elect (see Equation 1) on monomer deformations, the variation of the CO and N₂ molecular multipole moments with the internal coordinates (and consequently, also that of the related charge distributions) was introduced. As for N₂, the molecular multipoles and the related point charge dependence as those adopted in Lombardi et al. (2016b) were used. In the case of CO, as previously done for CO₂ (Bartolomei et al., 2012) and N₂ (Lombardi et al., 2016b), multireference ACPF (Averaged Coupled Pair Functional) calculations were performed as a function of the stretching for the molecular dipole and quadrupole moments by following the guidelines reported in Bartolomei et al. (2011). In particular, the molecular orbitals in the ACPF calculations were all taken to be these natural ones from the complete active space self-consistent field (CASSCF) reference wave functions. Therefore, the considered active space (CAS) was assumed to distribute 10 electrons in the orbitals indicated as (2, 3, 4)σ_g(2, 3, 4)σ_u(1, 2)π_u(1, 2)π_g. The 1σ_g1σ_u core molecular orbitals were fully optimized, while being constrained to be doubly occupied and excluded from the used CAS. The Dunning's aug-cc-pV5Z basis set (Kendall et al., 1992) was employed and the calculations were performed using the MOLPRO package (Werner et al., 2006).

The parameters of the rigid monomers defining the empirical PES were optimized in order to best fit at the same time the measured second virial coefficients and the ab initio interaction energies. It is worth pointing out here that, given the physical ground of the procedure, both the number of parameters allowed to vary in the fit and their interval of variation are rather small. As an example, the long-range dispersion attraction coefficient values, defined as ε·Rm6 were allowed to vary within a 10% interval of their initial value (see also Appendix A of Cappelletti et al., 2008). Moreover, it is also worth pointing out here that, due to some inter-dependencies, the final best fit values of the adjustable parameters usually differ only by a few percent from the initial ones (see Table 1 for ε and R_m). The main features (equilibrium distance R_e and binding energy D_e for selected dimer geometries) of the optimized and predicted empirical PES's are reported in Table 2, and it can be seen that the small variations of the optimized ϵ and R_m (see Table 1) lead in general to a less attractive potential. Table 2 shows also that while optimized and predicted parameter values are both able to satisfactorily reproduce the relative stability of the limiting configurations (see Table 2), the former are closer to ab initio results, especially for the most attractive dimer geometries. However, it can also be noticed that the model PES is not capable of predicting the T_a configuration as the most stable of those here considered; nevertheless, it has to be stressed that the ab initio energy difference between the T_a, T_b1, and T_b2 configurations is quite low and less than 1 meV (see Table 2) and that the very few parameters of the model PES do not allow the description of this very fine behavior.

As anticipated above, to validate and optimize the rigid rotor PES, second virial coefficient (B(T)) values, including first quantum correction B_q1(T) to the classical estimate B_cl(T) (Pack, 1978) (B(T) = B_cl(T)+B_q1(T)), were also computed. The calculations evidence that quantum corrections are smaller than 1% even for temperature values as low as 200 K. A comparison of calculated B(T) values with experimental measurements (Jaeschke et al., 1988; McElroy and Buchanan, 1995) over the temperature range 273 < T < 350 K is given in Figure 3.

The figure shows that the predicted PES results underestimate the experimental data in the considered range of temperature, while those referring to the optimized PES lie between the two sets of the measured values. This behavior suggests that, in agreement with the above reported analysis for the ab initio energy values, the predicted PES provides an interaction that is too attractive, which indeed can be properly corrected by tuning the involved ε and R_m parameters.

A further check of the reliability of the empirical PES is reported in Figure 4, where the related spherical averaged potential is compared with an accurate ab initio estimate at the MP4 level and obtained from Karimi-Jafari et al. (2011). Again, it has to be pointed out that although the predicted PES shows a 15% deeper well, it becomes practically negligible after optimization.

The parametrization of the empirical PES has been generalized by properly introducing the dependence of V_vdW and V_elec (see Equation 1) on the molecular elongation. Specifically, the variation of the electric multipole affecting V_elec has been obtained as indicated above, while the change of the molecular polarizability has been taken into account for the modulation of the V_vdW potential parameters.

For the N₂ monomer deformation, the empirical bond length dependence of the molecular polarizability α as reported in Appendix B of Cappelletti et al. (2008) is used, while for the corresponding point charges distribution dependence the analytic formulas given in Appendix A of Lombardi et al. (2016b) are employed.

In the case of the CO monomer we used the empirical bond length dependence of the polarizablity α as reported in Appendix A, which is shown here in the lower panel of Figure 5 together with previous ab initio estimates. In addition, for the CO permanent dipole quadrupole dependence on the bond length, an appropriate representation can be obtained from the analysis of present ab initio results reported in the upper panels of Figure 5. To this end, ab initio data were fitted using suitable analytic functions providing the radial dependence of point charges on the CO bond length r, as shown in detail in Appendix A.

The effect of the stretching of a single monomer on the complex interaction predicted by the adopted potential formulation is illustrated in Figures 6, 7. In all cases the plots of the CCSD(T) values, obtained for rigid and flexible molecules, are also reported for comparison: it can be appreciated that the elongation of one monomer provokes the same trends in both PESs for the stability of the considered interaction profiles, confirming therefore the reliability of the adopted empirical model.

According to the QC approach, the quantum mechanical time-dependent Schrödinger equation for the nuclear motion is solved for the degrees of freedom of the system playing the most relevant role in vibration-to-vibration (VV) quantum energy exchange, vibrations and roto-vibrational couplings, by a coupled equations method to obtain the quantum transition amplitudes avv′(t), where v and v′ are the initial vibrational quantum numbers of the diatoms. The vibrational wavefunction is initialized as the product of the Morse functions ϕv0(rCO)ϕv′0(rN2) for the two infinitely separated diatoms, whose parameters are given in Table 3, and is expanded as:

where E_v and Ev′ are the vibrational energies of the oscillators in their initial states. The remaining degrees of freedom are treated classically by integrating the corresponding set of Hamilton equations of motion in an effective potential defined as the Ehrenfest average of the interaction potential.

The quantum transition amplitudes avfvf′(t) can then be used to calculate either specific state-to-state vibrational/rotational transitions cross sections, σvijivi′ji′→vfjfvf′jf′, or Monte Carlo averaged cross sections over the Boltzmann distribution of the initial rotational angular momenta j and j′ for CO and N₂, respectively:

where T₀ is an arbitrary reference temperature, which cancels out in the formulation of rate constants (Equation 9), I_CO and I_N2 are the moments of inertia of the diatoms, U is the classical energy, obtained by subtracting from total energy the vibrational energy of the two diatoms, U=E-Ev-Ev′, j_max and jmax′ are the upper limit for the randomly chosen rotational quantum numbers for the diatoms, l_max the upper limit for the angular momentum and μ is the reduced mass of the system.

From such averaged cross sections, rate coefficients for vibrational relaxation kvv′(T) (see e.g., Coletti and Billing, 2002; Billing et al., 2003) can be derived as follows:

where ϵ_min = 0 for exothermic and ϵ_min = ΔE for endothermic processes and U¯ is the symmetrized classical energy (Billing, 1984b):

which has been introduced to restore, in an approximate fashion, the quantum mechanical detailed balance principle (Billing, 1984a,b, 1987).

In the present treatment the anharmonic vibrational energy is formulated as:

where ω_ei is the wavenumber for the i-th oscillator and x_ei and y_ei are the anharmonicity constants (see Table 3 for the values employed in the calculation for N₂ and CO).

QC rate coefficients have been computed for the exothermic exchange of a single vibrational quantum of energy between the first excited vibrational level vi′=1 of N₂ (N₂(1)) and the ground vibrational level v_i = 0 of CO (CO(0))

and for the endothermic exchange of a single vibrational quantum of energy between the first excited vibrational level v_i = 1 of CO (CO(1)) and the ground vibrational level vi′=0 of N₂ (N₂(0))

and in the temperature range 80–3,000 K, by running trajectories at 15 initial values of total classical energy comprised between 50 cm⁻¹ and 15,000 cm⁻¹, with a more frequent sampling driven toward lower energies. For each energy value, 2,000 trajectories were considered, which should ensure an accuracy for rate constants of ~ 20% at lower temperatures and ~ 15% at higher ones.

For each trajectory, the impact parameter was randomly chosen between 0 and 10 Å, and the initial separation between the colliding partners was set equal to 15 Å. In the expansion (7) a total of 36 vibrational states was included, i.e., those consisting of a band of energy ΔE up to 14,000 cm⁻¹.

Figures 8, 9 show the dependence of the computed QC rate coefficients on the temperature for the exothermic CO(0)+N₂(1) → CO(1)+N₂(0) and the endothermic CO(1)+N₂(0) → CO(0)+N₂(1) processes, respectively.

For comparison, the same Figures show the experimental data (Sato et al., 1969; von Rosenberg et al., 1972; Mastrocinque et al., 1976; Allen and Simpson, 1980), when available, together with previous theoretical results (Kurnosov et al., 2003).

Rate coefficients computed at some selected temperature values are also reported in Tables 4, 5.

The behavior of the QC rate coefficient plots computed on our PES shown in Figures 8, 9 for both processes agrees well with the experimental ones over the whole temperature range. More in detail, deviations by a maximum factor of 3 are found at the lowest values of T, where the experimental data are more bound to be affected by lower accuracy, although the corresponding slope is well reproduced. Indeed the calculated and experimental value of β in the Arrhenius expression k(T)=Aexp(-βRT) for the T range 80–300 K differs for little less than 10%.

In the case of the exothermic VV exchange, where experimental data are available in a wider range of temperatures up to 3,000 K, the behavior at low temperatures is the same as for the endothermic process, but the deviation between calculated and experimental rate coefficients is found to decrease to ~ 20% at temperatures larger than 1,000 K (Figure 8). Thus, the agreement between theoretical and experimental data is rather good in the whole temperature range. Agreement at higher temperatures (superior than in previous theoretical determination) is an indication that the basic contributions to the interaction are well described. The slightly larger discrepancy at low temperature can be attributed either to the lower accuracy of low-temperature experiments or to the neglecting of a proper quantum treatment for rotations which might play a role in the vibrational energy exchange process at low collision energies.

Networking and assembling networked software applications for the Molecular science community begun within COST (www.cost.eu/) Action D23 (METACHEM) and D37 (Grid Computing in Chemistry: GRIDCHEM). In METACHEM the activities of various Molecular Science research laboratories were networked on a shared computing platform made of a geographically distributed cluster of heterogeneous computers operating as a single virtual parallel machine (Foster and Kesselman, 1999). In the following Action GRIDCHEM Grid solutions and paradigms for molecular science research developed by D23 were consolidated on the grid. GRIDCHEM leveraged the creation and the use of distributed computing infrastructures (the “Grid”) to drive collaborative computer modeling and simulation in chemistry toward “new frontiers in complexity and a new regime of time-to-solution” (GRIDCHEM, 2006–2010). At more infrastructural level, the European projects EGEE (Enabling Grids for E-sciencE, https://cordis.europa.eu/project/rcn/87264_en.html), first, and the EGI (European Grid Infrastructure, https://en.wikipedia.org/wiki/European-Grid-Infrastructure), next, provided various disciplines, including Molecular Science, with a world class level platform for computational collaborations. In particular, during EGEE-III, the first pilot computational application called GEMS (Laganà et al., 2010)—designed to enable the distributed calculation of cross sections and rate coefficients starting from the ab initio computation of the electronic structure of the molecular system followed by the fitting of the computed ab initio points using a combination of analytic formulae into a proper functional representation of the interaction in short, intermediate and long-range regions—was implemented. This allows at present routine computation of reactive and nonreactive properties of elementary systems based on the desired number of molecular geometries (Storchi et al., 2006) and collisional paths (Gervasi and Laganà, 2004) by distributing them on the grid. The scheme of the most recent evolution of GEMS toward the atomistic simulation of the kinetics of more complex systems, together with the related “data analysis and validation” and “open archive and reuse” in a cloud service perspective, is sketched in Figure 10.

The molecular science community has been accordingly organized into the COMPCHEM (Laganà et al., 2010) VO first and in the CMMST (Chemistry, Molecular and Materials Science a Technologies) VRC (Virtual Research Community) (https://wiki.egi.eu/wiki/Towards_a_CMMST_VRC) later. On a national basis, then, the different disciplinary communities were gathered in Joint Research Units (JRU), for example in Italy with IGI (http://www.italiangrid.it), in order to give a higher local momentum to networked applications. In this spirit, the following activities were pursued in networked Molecular Sciences:

orchestrate the initiatives of the e-infrastructure experts and of the disciplinary researchers so as to enable an effective intra- and trans-community networked implementation and coordination of a collaborative/competitive (synergistic) research environment by allowing a selection of compute resources based on quality parameters;

The cloud-oriented computing infrastructure developed for our calculations is an embryonic platform of the Beowulf type running under the OpenStack platform sketched in the lower side of Figure 3 of Vitillaro and Laganà (2018) named HERLA (https://en.wikipedia.org/wiki/Beowulf_cluster). HERLA has been established at the Dipartimento di Chimica, Biologia e Biotecnologia (DCBB) of the University of Perugia. The platform consists of a couple of HPC clusters, (CG/training) and (FE/research) running Scientific Linux 6.x with two distinct access nodes. The clusters are connected using NIS in a single-image system and are used first for students' training (CG) and second for scientists' research (FE). The management of Herla is presently carried out by the CMS² Consortium (http://www.cms-2.org/index.php) of the University of Perugia, CNR-ISTM - Perugia unit and the two companies Master-UP s.r.l. and Molecular Horizon s.r.l. Recently, cloud images of HERLA (VHERLA) have been created and deployed in a storage system (CEPH located at DCBB). This effort was meant to support as well the activities of the School on Open Science Cloud (SOSC17) held in Perugia on June 2017 in collaboration with the Department of Physics and Geology (DFG) and INFN Perugia (running under the INFN OpenStack platform sketched in the lhs upper side of Figure 3 of Vitillaro and Laganà, 2018).

Later, as sketched in the rhs side of Figure 4 of Vitillaro and Laganà (2018), the images of VHERLA were allocated to the OpenStack GARR Cloud platform in Palermo within the project “cnr-istm” and were used to install the version VHERLA(GARR-CLOUD) hscw (http://hscw.herla.unipg.it/ganglia/?p=2&c=FrontEnd) bearing the following features: Access node hscw (2 core, 4Gb RAM, 200Gb storage) and Cluster (7 nodes, 96 cores, 360Gb RAM, 700Gb storage). The Access node hscw, can be reached at the IP address [90.147.189.20] via SSH and the following 7 nodes, 96core, 380Gb, 512Gb scratch are defined at Torque(PBS)/MAUI as [Intel Xeon E3-12xx v2(Ivy Bridge)/2.6Ghz]. This allowed the generation of a virtual cluster for Molecular Sciences that has also been used for the training of the Students of the XIII EM TCCM (European Master in Theoretical Chemistry and Computational Modeling) Intensive Course (http://www-old.chm.unipg.it/chimgen/mb/theo2/TCCM2018/EM- TCCM2018/EM-TCCM/Welcome.html) with details at the following reference web page: http://hscw.herla.unipg.it URL.

As to the operational model we have already developed a sustainable open collaborative user/producer (Prosumer), whose prototype has been first implemented for running the ECTN (European Chemistry Thematic Network (http://ectn.eu/)) EChemTest e-tests based on the use of the educational LibreEOL (https://echemtest.libreeol.org) and GLOREP (https://glorep.unipg.it/) services Laganà et al., 2018a. An important feature of the Prosumer model adopted for EChemTest was the containment of costs by leveraging the fact that the ECTN member HEIs running e-test Self Evaluation Sessions (SES)s for the assessment of the Chemistry competences of their own students already act at the same time as consumers of EChemTest services and as producers of Question and Answers (Q&A)s, assessors of students' competences, designers and developers of e-learning materials (a typical cluster of the horizontal type) for the harmonization at the European level of the assessment of Molecular Science competences. This behavior is quite usual in education and research activities in which knowledge is a common good to be at the same time produced and consumed.

The only additional action needed to the end of making the prosumer model sustainable was to assign to a company the role of market spinner. For EChemTest this role was taken by Master-Up s.r.l. thanks to its nature as a former spinoff of the University of Perugia (started in the year 2004 out of the aggregation of some members of the Chemistry Department and the Mathematics and Informatics Department who were experts in molecular dynamics simulations and computer science) devoted to design, production and marketing services for technological innovation as its main goal. As a matter of fact, the mission of Master-Up has been since the very beginning the design and development of cloud services for molecular sciences and technologies. In the particular case of EChemTest, the Prosumer model has led to the production in the year 2018 of 2622 SESs, with an increase of 5 percent over the previous year. Further efforts have also been spent in the period up to January 2019 for the implementation of the already mentioned MOlecular Simulator Enabled Cloud Services (MOSEX), a European Open Science Cloud Pilot designed as a follow-up of the EGI COMPCHEM VO activities for the following applications:

Molecular electronic structure and dynamical properties programs of GEMS (including the QC ones reported reported in this paper)

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.