Quantum-Chemical Investigations on the Structure and Stability of Mixed Trimers Containing HC₃N in Combination with H₂C₂ and/or HCN Analyzed by QTAIM, NBO and SAPT Methods

Abstract

The present work deals with the computational study of HC₃N$··$HCN$··$H₂C₂-, (HC₃N)₂$··$H₂C₂-, and HC₃N$··$(H₂C₂)₂-mixed trimers. The different equilibrium structures of the different low-lying minima on the corresponding potential energy surface (PES) were accurately determined, and the relative stabilities were computed by extrapolation procedures to the complete basis set limit. For each mixed trimer, the non-covalent interactions ruling the structure of the most stable isomer were analyzed using the QTAIM (Quantum Theory of Atoms in Molecules) approach. Additional insights into these interactions were provided by the Natural Bond Orbital (NBO) and Symmetry-Adapted Perturbation Theory (SAPT) methods. These results can be used to assist further theoretical investigations and experimental studies on the formation of larger molecules potentially relevant in astrochemistry.

1. Introduction

Hydrogen cyanide (HCN), cyanoacetylene (HC₃N), and acetylene (H₂C₂) are very simple neutral organic molecules that are very relevant in astrochemistry. The first one has been discovered in many different regions of the interstellar medium (ISM) [1,2,3,4]; a similar pervasion in the ISM was also found for cyanoacetylene [5,6,7] and acetylene [8,9]. Due to their abundance in the interstellar medium, they have been extensively investigated (see Refs. [10,11,12,13,14] and references therein). In addition, recent studies [15,16] highlighted that neutral weakly bound complexes made by these molecules may be involved in the first steps of the formation path of larger compounds. In the literature, clusters made up of two or more units of the same molecule were investigated for hydrogen cyanide (see Refs. [17,18] and references therein), cyanoacetylene (see Refs. [15,19] and references therein), and acetylene (see Refs. [20,21,22] and references therein). Moving to the corresponding mixed clusters (i.e., made up of different molecules), the studies available in the literature are mainly focused on dimers. In the case of cyanoacetylene, its mixed dimers with acetylene [23], di-acetylene [24], ammonia [25], carbon monoxide [26] and carbon dioxide [27] were analyzed. As for the mixed dimers involving at least one unit of acetylene, their structures with hydrogen cyanide were investigated [28,29]. Indeed, in addition to dimers, larger mixed clusters are also becoming increasingly important in astrochemistry due to their potential connection to the formation path of molecules [30], and therefore, we decided to focus our attention on weakly bound mixed aggregates containing at least one unit of cyanoacetylene (HC₃N) and one unit of acetylene (H₂C₂).

In the present work, we report on the results obtained by a computational study carried out on HC₃N$··$HCN$··$H₂C₂-, (HC₃N)₂$··$H₂C₂-, and HC₃N$··$(H₂C₂)₂-mixed molecular complexes. For each trimer, we accurately characterized the equilibrium geometries corresponding to the different minima on the corresponding potential energy surface (PES), and we computed their relative stabilities; in addition, by employing different approaches, we analyzed the non-covalent bonding interactions ruling the most stable structure of each trimer. These results provide the basic data necessary to guide future spectroscopic experimental investigations on these weakly bound clusters; in addition, these data could assist searches in remote environments and could also demonstrate usefulness in modeling the several steps leading to the formation of more complex molecules.

2. Materials and Methods

To obtain the geometries corresponding to the different minima for each of the mixed trimers investigated, we followed the procedure established in a previous work [31], which is only briefly summarized here. The efficient sampling of the PES is obtained by carrying out three distinct steps, each one involving density functional theory (DFT) calculations performed at an increasingly higher level of theory. At the end of each step, redundant structures were excluded using a mix of criteria based on the mean squared deviations of their geometries and their energy difference. The first two steps employed first the GFN2-xtB method [32] and then the B3LYP functional [33]. Finally, in the last step, the unique isomers were computed and optimized using the recently proposed double hybrid functional revDSD-PBEP86 [34] in conjunction with the “calendar” jul-cc-pVTZ basis set [35]; the combination of this functional method with the triple-zeta calendar basis set, hereafter labeled simply as revDSD, has been proven to yield accurate predictions for geometries, dipole moments, and spectroscopic parameters (see Refs. [36,37,38] and references therein). All the DFT calculations were performed with the inclusion of D3 corrections [39] and Becke–Johnson damping [40] since, despite the widespread use of DFT methods in different chemical applications (see, for example, Refs. [41,42,43] and references therein), the accurate modeling of van der Waals adducts requires the inclusion of dispersion contributions [44,45,46]. Subsequent frequency calculations (carried out at the same level of theory) confirmed that all these isomers were the real minima on the PES (i.e., no imaginary frequencies). At the same level of theory (revDSD), the structures of HC₃N-, HCN-, and H₂C₂-isolated monomers were optimized; the corresponding data (equilibrium geometries, harmonic frequencies, and intensities of the fundamentals) are reported in the Supplementary Material (see Table S1).

Accurate predictions of the binding energies (BEs) were computed by extrapolation to the complete basis set (CBS) limit using two composite schemes, each one relying on a series of single-point calculations carried out using both MP2 [47] and CCSD(T) [48] levels of theory, thus efficiently accounting for both the basis set superposition effects and the basis set truncation errors. In the first one (hereafter labeled as CBS-1), the BEs were derived by correcting the MP2/CBS values obtained using the aug-cc-pVnZ (n = T, Q and 5) basis sets [49,50] with CCSD(T)/aug-cc-pVTZ energies; core-valence (CV) corrections were obtained at the MP2 level using the cc-pCVTZ basis set [51]. As a second extrapolation approach, we employed the so-called jun-ChS scheme [52], which was recently designed to accurately evaluate the non-covalent interactions ruling molecular complexes [53,54]. For comparison, in the present work, we also reported the interaction energies (counterpoise corrected) obtained at the revDSD level of theory, hereafter labeled as IE/CP-revDSD). A more realistic determination of the binding energies of these clusters should consider at least the correction for zero-point energy; however, following previous work [53], in the present paper we decided not to include these kinds of contributions.

When dealing with weakly bound complexes, the analysis of the non-covalent interactions ruling their structures is generally performed by investigating the corresponding electron density $\rho$ (see Ref. [55] and references therein). So, the most stable isomer for each of the three mixed clusters was studied by two different methods, namely the QTAIM (Quantum Theory of Atoms in Molecules) [56] and NCIs (non-covalent interactions) [57]; some of their recent applications to molecular complexes can be found in Refs. [58,59] (and the references therein). Within the framework of the QTAIM, key information about the nature of bonding interactions can be obtained by the characterization of so-called critical points (CPs) locations in space where the gradient norm of the electron density is zero; these critical points are usually classified into four different kinds according to the number of the negative eigenvalues ${\lambda }_{\mathit{i}}$ of the corresponding Hessian matrix. We focused our attention on CPs corresponding to second- and first-order saddle points that are associated with the presence of a bond between two neighboring atoms and several bonds involved in a ring, respectively. NCI instead employs the reduced density gradient (RDG) function, together with the sign of the second largest eigenvalue of the Hessian matrix of electron density (${\lambda }_2$) to plot the regions associated with non-covalent interactions; the strength and type of each interaction can be inferred from the value of $\rho$ and the sign of ${\lambda }_2$, respectively. In the present work, we used both the color-filled isosurface 3D plot and the color-mapped 2D scatter graph (both using the proper color map to graphically represent the different kinds of interactions). Besides the results yielded by the QTAIM and NCI methods, further insights into these interactions were obtained using the Natural Bond Orbital (NBO) approach, thus computing, for each donor–acceptor couple, the corresponding stabilization energy E⁽²⁾ as yielded by second-order perturbation analysis [60]; the greater the interaction between the donor and acceptor NBO orbitals, the larger the corresponding value of E⁽²⁾. Eventually, the most stable structures of (HC₃N)₂$··$H₂C₂- and HC₃N$··$(H₂C₂)₂-mixed molecular complexes were analyzed using the Symmetry-Adapted Perturbation Theory (SAPT) energy decomposition method [61]; within this approach, the non-covalent interactions ruling these structures are broken down into physically meaningful components, namely the electrostatics, exchange-repulsion energy, induction, and dispersion terms. These kinds of calculations were carried out using the so-called gold standard level (i.e., SAPT2+(3)δMP2/aug-cc-pVTZ [62,63]. Some recent applications of NBO and SAPT to the investigations of molecular complexes are reported in Ref. [64] (and the references therein).

All the electronic calculations were carried out using Gaussian16 software (version C.01) [65], while QTAIM and NCI analysis were performed using the Multiwfn (rev. 3.8) program [66]; SAPT energy decompositions were made by employing the PSI4 (version 1.9.1) software package [67].

3. Results and Discussion

3.1. HC₃N$··$HCN$··$H₂C₂ Clusters

By employing the procedure previously outlined, we identified nine different low-lying isomers for the clusters containing one unit of cyanoacetylene, one unit of hydrogen cyanide, and one unit of acetylene (the corresponding equilibrium geometries, dipole moments, and BE values are reported as Supplementary Data, see Table S2); in Figure 1, the structures of these nine lowest energy geometries are shown.

The structure labeled A (which resulted in being the most stable) presents the three different units aligned into a linear structure; the distance between the nitrogen atom of cyanoacetilene and the hydrogen atom of HCN is 2.163 Å, while the distance between the nitrogen atom of HCN and the hydrogen atom of H₂C₂ is 2.312 Å. In B, we present planar cyclic geometry, with the two units of acetylene and cyanoacetylene almost parallel. In C again, a cyclic geometry is given, but in this case, the units of cyanoacetylene and hydrogen cyanide are parallel (and the hydrogen and nitrogen atoms are distributed in a sort of head–tail arrangement), with the acetylene unit almost perpendicular to both. In D, we present acetylene and cyanoacetylene almost parallel to each other, with the hydrogen cyanide unit almost perpendicular to both. E, F, and G show planar structures, H corresponds again to a linear arrangement of the monomer units, and I is characterized by the HC₃N and HCN units both being perpendicular to the acetylene monomer.

Table 1. Equilibrium rotational constants (in MHz), overall dipole moments (in Debye), and components along the axes, binding energies (BEs, in kcal mol⁻¹), and interaction energies (IEs, in kcal mol⁻¹) for the isomers of HC₃N$··$H₂C₂$··$HCN clusters.

	Isomer A	Isomer B	Isomer C	Isomer D
A	$-$	2851	2461	2508
B	248	858	953	877
C	$-$	659	687	650
Overall dipole moment (|μ|)	8.82	1.06	0.980	6.36
Dipole moment component (μx, μy, μz)	0.0, 0.0, 8.82	−1.06, 0.77, 0.0	−6.34, −0.75, 0.0	6.36, −0.27, 0.23
BE/CBS-1 a	−7.85	−7.55	−7.52	−6.98
BE/jun-ChS b	−7.65	−7.54	−7.47	−6.94
IE/CP-revDSD c	−7.82	−7.57	−7.60	−6.97
	Isomer E	Isomer F	Isomer G	Isomer H
A	2310	2164	1851	$-$
B	658	739	654	229
C	512	551	484	$-$
Overall dipole moment (|μ|)	7.92	5.61	7.62	0.889
Dipole moment component (μx, μy, μz)	7.75, −1.65, 0.0	−5.61, 0.09, 0.0	−7.3, −2.27, 0.0	0.0, 0.0, −0.889
BE/CBS-1 a	−6.43	−6.16	−5.62	−4.44
BE/jun-ChS a	−6.34	−6.13	−5.54	−4.25
IE/CP-revDSD b	−6.38	−6.14	−5.59	−4.41
	Isomer I
A	35184
B	264
C	262
Overall dipole moment (|μ|)	0.755
Dipole moment component (μx, μy, μz)	0.76, 0.0, 0.0
BE/CBS-1 a	−3.68
BE/jun-ChS a	−3.71
IE/CP-revDSD b	−3.76

^a Obtained using the CBS-1 extrapolation scheme (see the text for details). ^b Obtained using the jun-Chs extrapolation scheme (see the text for details). ^c Obtained using the CP-revDSD method (see the text for details).

reports the equilibrium rotational constants for these different isomers, together with the corresponding predicted values of their dipole moments; in the same Table, the binding energies (BEs) computed by employing both CBS-1 and jun-ChS approaches are also listed (together with the predictions given by revDSD calculations).

Both these two extrapolation schemes identified the linear isomer A to be the most stable one; it is worthwhile to note that the other linear isomer we found (H, see also the geometry in Table S2), with only the position of acetylene and hydrogen cyanide units swapped, is predicted to be less stable by about 3.4 kcal mol⁻¹. The predictions of the BEs yielded by these two approaches are in very close agreement (their mean absolute difference, MAD, is very small, 0.08 kcal mol⁻¹). Other spectroscopic parameters (harmonic frequencies and intensities) useful to guide future spectroscopic investigations on isomer A are listed in

Table 2. Harmonic frequencies and intensities for isomer A of HC₃N$··$H₂C₂$··$HCN.

Freq (cm−1)	Intensity (km mol−1)	Freq (cm−1)	Intensity (km mol−1)	Freq (cm−1)	Intensity (km mol−1)
9.7	1.58	511.8	5.90	2111.5	27.97
9.7	1.58	511.8	5.90	2115.8	15.04
30.4	0.01	647.2	10.95	2309.7	49.09
30.4	0.01	647.2	10.95	3336.1	655.37
69.2	2.06	694.8	36.31	3373.9	215.92
100.0	5.18	694.8	36.31	3456.8	81.35
100.0	5.18	813.2	62.91	3492.3	1.35
136.1	0.35	813.2	62.91
144.4	19.34	830.7	41.15
144.4	19.34	830.7	41.15
230.5	0.52	894.2	0.24
230.5	0.52	2001.7	13.76

.

At this point, it is worthwhile to list the most relevant changes in the monomers after the formation of isomer A (see Tables S1 and S2): in the H₂C₂ unit, we observed an increase (by 0.005 Å) in the C-H bond pointing toward the N atom of HC₃N; in the HCN unit, we saw an elongation of the C-H bond by 0.008 Å; and the HC₃N unit was characterized by the slight contraction (by 0.002 Å) of its proton–acceptor C-N triple bond. Comparing the harmonic data of the isomer A (Table 2) with the data of the isolated monomers (Table S1), we can see that, associated with these structural changes, there are some relevant shifts in the frequencies as well.

The structure of the most stable isomer (A) has been further investigated using the QTAIM method. The topological analysis pointed out that the geometry of this isomer is due to the presence of two attractive bond (3,−1) critical points labeled as CP1 and CP2 and highlighted by the arrows between the nitrogen and hydrogen atoms of the different monomers, as shown in Figure 2.

At these two CPs (highlighted by the two arrows), the values of the two most relevant QTAIM parameters, namely the electron density, $\rho \left(r\right)$, and its Laplacian, ${\mathsf{\nabla }}^2\rho \left(r\right)$, can be used to assess the type of interactions observed. For CP1, we obtained 0.017 and 0.060 au for $\rho \left(r\right)$ and ${\mathsf{\nabla }}^2\rho \left(r\right)$, respectively, while CP2 was characterized by values of 0.012 au for $\rho \left(r\right)$ and of 0.046 au for ${\mathsf{\nabla }}^2\rho \left(r\right)$. The commonly accepted ranges for hydrogen bonding are 0.002–0.034 au for $\rho \left(r\right)$, and 0.024–0.139 au for ${\mathsf{\nabla }}^2\rho \left(r\right)$ [68], so both CP1 and CP2 can be associated with hydrogen bonds between the different units, with the former showing a larger value of ${\mathsf{\nabla }}^2\rho \left(r\right)$. The NCI analysis confirmed the existence of these non-covalent interactions among the different units; see Figure S1 (where these hydrogen bonds are pointed out by the green regions in the 3D isosurface plot and the green color spikes in the 2D scatter graph). In the NBO framework, the CP1 highlighted by QTAIM is mainly described as charge transfer (with an E⁽²⁾ of 6.38 kcal mol⁻¹) from the lone pair of the N atom of HC₃N to a nonbonding orbital between the C-H pair of the HCN unit. In the same way, CP2 is mainly related to charge transfer (with a stabilization energy of 3.02 kcal mol⁻¹) from the lone pair of the N atom of HC₃N to a nonbonding orbital of the C-H bond of the H₂C₂ unit.

3.2. (HC₃N)₂$··$H₂C₂ Clusters

In the case of the mixed trimers containing two units of cyanoacetylene and one unit of acetylene, we found seven different low-lying structures (the values of the corresponding equilibrium geometries, dipole moments, and BEs can be found in the Supplementary Data, see Table S3); in Figure 3, we report their geometries.

The structure labeled A (the most stable one) presents cyclic geometry among the three units; the two units of cyanoacetylene are almost parallel between them (in a sort of head–tail arrangement), and one of them has one hydrogen pointing towards the pi-cloud of the acetylene unit. The distance between the two centers of mass of the two HC₃N units is 3.578 Å, while the distance between the center of mass of H₂C₂ and that of the whole complex is 4.014 Å. In B, a linear structure can be seen. In C, we find a T-shaped geometry, but in this case, the two units of cyanoacetylene are colinear, and the hydrogen atom of cyanoacetylene points toward the p-cloud of the acetylene unit. In D, another cyclic geometry is observed; this time, the unit of acetylene is almost parallel to one unit of cyanoacetylene, and one of its hydrogen atoms is pointing toward the other cyanoacetylene unit. The other remaining structures (E, F, and G) are all planar. The corresponding equilibrium rotational constants and overall dipole moments (and components along the axes corresponding to the geometries listed in Table S3) are listed in

Table 3. Equilibrium rotational constants (in MHz), overall dipole moments (in Debye), and components along the axes, binding energies (BEs, in kcal mol⁻¹), and interaction energies (IEs, in kcal mol⁻¹) for the lowest-energy isomers of (HC₃N)₂$··$H₂C₂ clusters.

	Isomer A	Isomer B	Isomer C	Isomer D
A	1505	$-$	35219	2484
B	670	153	159	379
C	464	$-$	158	329
Overall dipole moment (|μ|)	0.15	9.84	9.69	7.75
Dipole moment component (μx, μy, μz)	0.14, −0.05, 0.0	−9.84, 0.0, 0.0	−9.69, 0.0, 0.0	−7.75, −0.06, 0.0
BE/CBS-1 a	−8.29	−7.12	−6.76	−6.11
BE/jun-ChS b	−8.33	−6.90	−6.63	−6.07
IE/CP-revDSD c	−8.39	−7.13	−6.86	−6.11
	Isomer E	Isomer F	Isomer G
A	2310	2164	1851
B	658	739	654
C	512	551	484
Overall dipole moment (|μ|)	8.87	8.75	3.97
Dipole moment component (μx, μy, μz)	8.77, −1.41, 0.0	8.74, −0.21, 0.0	1.23, −3.77, 0.0
BE/CBS-1 a	−5.94	−5.70	−5.19
BE/jun-ChS b	−5.87	−5.65	−5.16
IE/CP-revDSD c	−5.94	−5.72	−5.15

^a Obtained using the CBS-1 extrapolation scheme (see the text for details). ^b Obtained using the jun-Chs extrapolation scheme (see the text for details). ^c Obtained using the CP-revDSD method (see the text for details).

.

Both CBS-1 and jun-ChS schemes were used to rank the different isomers according to their stability, and isomer A was identified as the most stable one (see the binding energies reported in Table 3). Also, for this kind of cluster, these two extrapolation methods gave results that are in very close agreement, including their MAD at only 0.08 kcal mol⁻¹.

Table 4. Harmonic frequencies and intensities for isomer A of (HC₃N)₂$··$H₂C₂.

Freq (cm−1)	Intensity (km mol−1)	Freq (cm−1)	Intensity (km mol−1)	Freq (cm−1)	Intensity (km mol−1)
24.7	3.55	241.2	1.06	805.6	137.46
41.7	5.50	512.2	6.36	889.4	0.61
47.8	1.05	512.8	0.99	891.3	0.12
61.4	6.61	514.2	8.35	2002.0	3.26
75.4	1.21	514.6	8.43	2098.5	14.94
83.6	3.75	648.1	6.56	2107.4	3.02
89.1	0.12	656.9	1.86	2294.4	32.76
102.5	4.32	694.0	44.40	2302.0	49.50
133.7	1.72	704.6	31.89	3378.2	158.15
230.9	0.12	742.6	14.34	3387.8	314.70
232.3	0.10	756.9	34.41	3452.2	76.60
236.2	0.34	799.6	82.67	3488.0	1.90

lists other spectroscopic parameters (harmonic frequencies and intensities) that can be used to assist future experimental works on isomer A.

Looking at the changes in the structural parameters of the monomers (see Table S1) upon the formation of this isomer (see the data in Table S3), we noted that the most relevant ones concern the unit of HC₃N with the C-H bond pointing towards H₂C₂ with the elongation (by 0.005 Å) of the C-H bond, and the H₂C₂ unit, where its C-H bond towards the N of the other HC₃N unit experiences a similar change. The other cyanoacetylene unit shows smaller structural changes. The shifts in the harmonic frequencies of the complex (Table 4) from the values of the individual units (Table S1) are in line with the structural changes.

The analysis of isomer A carried out within the QTAIM framework relates its geometry to the interplay between three attractive bond (3,−1) CPs (labeled as CP1, CP2, and CP3) and one (3,+1) ring critical point (labeled RCP and involving steric effects), among the three molecular units, as shown in Figure 4.

On the basis of the corresponding electron density and Laplacian values, both CP1 and CP3 can be classified as hydrogen bonding (CP1: $\rho \left(r\right)=$ 0.0099 au, ${\mathsf{\nabla }}^2\rho \left(r\right)=$ 0.030 au; CP3: $\rho \left(r\right)=$ 0.011 au, ${\mathsf{\nabla }}^2\rho \left(r\right)=$ 0.043 au). CP2 is characterized by the smaller value of its Laplacian parameter (${\mathsf{\nabla }}^2\rho \left(r\right)=$ 0.020 au), as expected in the case of van der Waals interactions. In Figure S2, we report the NCI analysis carried out on the same isomer, showing the non-covalent interactions among the different units (see the green regions in the 3D isosurface plot and the green color spikes in the 2D scatter graph). The NBO analysis correlates CP1 mainly to a very strong interaction (with an E⁽²⁾ of 2.42 kcal mol⁻¹) between a bonding orbital of H₂C₂ (donor) and a nonbonding orbital of HC₃N (acceptor); similarly, the CP2 between the two cyanoacetylene units is mainly associated with an interaction (with stabilization energy of 0.52 kcal mol⁻¹) between a bonding orbital of one HC₃N and a nonbonding orbital of the other. Finally, CP3 is mainly correlated to charge transfer (with stabilization energy of 1.86 kcal mol⁻¹) from the lone pair of the N atom of HC₃N (donor) to the nonbonding orbital of the C-H bond of the H₂C₂ unit (acceptor). Following the procedure of previous studies on mixed trimers [69,70,71], SAPT calculations were carried out by considering two different fragments, which, in this case, included the acetylene and the cyanoacetylene subunits, thus focusing on the interactions between these combined cyanoacetylene clusters and the acetylene. The corresponding results (see Table S4) show that isomer A is characterized by the largest contributions of electrostatic induction and dispersion factors to the overall stability; the most relevant contributor to the stability of isomers B and C is the electrostatic term, while for isomer D, it is the dispersion contributor.

3.3. HC₃N$··$(H₂C₂)₂ Clusters

Concerning the mixed trimers containing one unit of cyanoacetylene and two units of acetylene, our research procedure led to the characterization of five different low-lying geometries corresponding to minima in the PES (the data of the corresponding equilibrium geometries, dipole moments, and BEs are listed in the Supplementary Data, see Table S5); in Figure 5, we report these structures.

The structure labeled A (the most stable one) presents cyclic geometry among the three units where one unit of cyanoacetylene and one of acetylene are almost parallel (their centers of mass are separated by about 3.7 Å), and the other unit of acetylene is almost perpendicular to them (and its distance from the center of mass of the whole cluster is 3.548 Å). In B, we observe a T-shaped structure, where one unit of acetylene and cyanoacetylene are colinear, and the hydrogen atom of cyanoacetylene is pointing toward the p-cloud of the other acetylene unit. In C, a cyclic geometry is again presented, but in this case, the units of cyanoacetylene and acetylene are almost parallel, while the hydrogen atom of acetylene is pointing toward the p-cloud of the second acetylene unit. In D, another cyclic geometry is evident, but this time, the hydrogen of cyanoacetylene is pointing toward the p-cloud of one unit of acetylene, while one of the hydrogen atoms of the second unit of acetylene is pointing toward the cyanoacetylene unit. Finally, in E, we have another T-shaped structure. From the point of view of relative stability, both extrapolation procedures (CBS-1 and jun-ChS) clearly identified isomer A as being the favorite one; even in this case, the agreement between the predictions yielded by the two methods is remarkable (MAD of only 0.03 kcal mol⁻¹). For these isomers, we report the corresponding bonding energies (together with the predictions obtained at the revDSD levels of the theory), their equilibrium rotational constants, and their dipole moments in

Table 5. Equilibrium rotational constants (in MHz), overall dipole moments (in Debye), and components along the axes, binding energies (BEs, in kcal mol⁻¹), and interaction energies (IEs, in kcal mol⁻¹) for the different isomers of HC₃N$··$(H₂C₂)₂ clusters.

	Isomer A	Isomer B	Isomer C	Isomer D
A	2573	35221	1729	1908
B	842	279	1057	716
C	634	277	656	521
Overall dipole moment (|μ|)	3.88	5.24	3.78	4.12
Dipole moment component (μx, μy, μz)	3.64, 1.39, 0.0	5.24, 0.0, 0.0	−3.24, −1.94, 0.0	3.95, 1,11, 0.0
BE/CBS-1 a	−5.33	−4.85	−4.79	−4.67
BE/jun-ChS b	−5.32	−4.79	−4.80	−4.71
IE/CP-revDSD c	−5.31	−4.90	−4.82	−4.76
	Isomer E
A	35242
B	247
C	245
Overall dipole moment (|μ|)	4.25
Dipole moment component (μx, μy, μz)	4.25, 0.0, 0.0
BE/CBS-1 a	−3.84
BE/jun-ChS b	−3.74
IE/CP-revDSD c	−3.84

^a Obtained using the CBS-1 extrapolation scheme (see the text for details). ^b Obtained using the jun-Chs extrapolation scheme (see the text for details). ^c Obtained using the CP-revDSD method (see the text for details).

.

Table 6. Harmonic frequencies and intensities for isomer A of HC₃N$··$(H₂C₂)₂.

Freq (cm−1)	Intensity (km mol−1)	Freq (cm−1)	Intensity (km mol−1)	Freq (cm−1)	Intensity (km mol−1)
14.5	3.51	509.8	7.07	800.8	107.45
41.5	3.49	512.4	4.38	889.9	0.08
53.4	0.18	638.2	5.22	2002.8	1.85
69.4	0.25	639.7	10.44	2005.5	8.06
70.3	0.76	646.6	6.80	2110.3	1.19
78.7	0.24	652.7	0.62	2305.3	25.25
88.5	0.00	684.5	35.82	3384.8	162.00
102.8	3.21	684.6	39.51	3391.6	233.70
135.3	0.43	781.4	24.04	3461.0	77.84
229.8	1.68	788.9	106.58	3492.3	2.27
230.1	0.03	799.1	117.80	3497.1	1.67

lists additional spectroscopic parameters (harmonic frequencies and intensities) that can be useful for the future experimental characterization of isomer A.

The comparison of the isomers’ structural data (Table S5) with the corresponding data of the isolated monomers (Table S1) shows that the most relevant changes are to be found in the two H₂C₂ units, with the elongation (by about 0.004 Å) of the C-H bonds pointing toward the acetylene and the N atom of HC₃N, respectively. Looking at the harmonic data reported in Table 6 and comparing them with the values listed in Table S1, we can see that the shifts in the vibrational frequencies upon the formation of the complex are in line with these structural changes.

The structure of the most stable isomer (A) was then investigated within the QTAIM framework; the analysis (see Figure 6) revealed that, as in the case of (HC₃N)₂$··$H₂C₂, the geometry of the most stable isomer of HC₃N$··$(H₂C₂)₂ is associated with the interplay between three attractive bonds (3,−1) critical points (labeled as CP1, CP2, and CP3), and one (3,+1) ring critical point localized among the different units (labeled as RCP), associated with steric effects.

Looking at CP3 ($\rho \left(r\right)=$ 0.011 au, ${\mathsf{\nabla }}^2\rho \left(r\right)=$ 0.043 au), the value of its Laplacian parameter clearly identifies it as capable of hydrogen bonding, while both CP1 ($\rho \left(r\right)=$ 0.0076 au, ${\mathsf{\nabla }}^2\rho \left(r\right)=$ 0.024 au) and CP2 ($\rho \left(r\right)=$ 0.0042 au, ${\mathsf{\nabla }}^2\rho \left(r\right)=$ 0.014 au) should be associated with weaker interactions. The main results of the NCI analysis are reported in Figure S3, whereas before, the non-covalent interactions were visualized as green regions in the 3D isosurface plot and green color spikes in the 2D scatter graph. The NBO analysis provides further insights, showing that there is a very strong interaction (having an E⁽²⁾ of 1.39 kcal mol⁻¹) between a bonding orbital of one unit of H₂C₂ (donor) and a nonbonding orbital of the second unit of H₂C₂ (acceptor), as highlighted by the QTAIM in terms of CP1. Between the second unit of acetylene and the cyanoacetylene molecule, the NBO analysis revealed weak interactions, mainly described as the charge transfer between the bonding orbital of H₂C₂ (donor) and a nonbonding orbital of HC₃N (acceptor) and between the lone pair of the nitrogen atom of HC₃N (donor) and the H₂C₂ unit (acceptor), correlating with the CP2 localized by QTAIM. Between the cyanoacetylene molecule and the first acetylene unit, NBO identified a significant charge transfer (with a stabilization energy of 1.91 kcal mol⁻¹) from the lone pair of the N atom of HC₃N (donor) to a nonbonding orbital of the C-H bond of H₂C₂ (acceptor), in line with the CP3 pointed out by the QTAIM. SAPT analysis was performed by considering the two acetylene molecules as one whole unit; the results (see Table S6) show that the overall stability of isomers A, C, and D is mainly due to the combined effect of both electrostatic and dispersion interactions, while isomer B is ruled by the dominant role of the former ones.

4. Conclusions

The present work carried out on mixed trimers made by at least one molecule of cyanoacetylene (HC₃N) and one molecule of acetylene (H₂C₂) demonstrated that these isomers possess a very rich PES (especially in the case of HC₃N$··$HCN$··$H₂C₂ clusters), with several different minima, the corresponding equilibrium geometries of which were accurately determined. The analysis carried out within the framework of the QTAIM theory highlighted that the structure of the most stable isomer (as identified by the accurate extrapolation schemes to the CBS limit here employed) is ruled by the interplay among different attractive critical points and, in the case of (HC₃N)₂$··$H₂C₂ and HC₃N$··$(H₂C₂)₂ trimers, is also ruled by the presence of one steric RCP. The stability of these isomers is due to several interactions, which have been thoroughly assessed by NBO analysis, characterized by significant values of the corresponding stabilization energies, and by SAPT calculations. This work clearly shows how the interplay of these different methods of analysis and the efficient extrapolation techniques for obtaining their binding energies can be combined to investigate mixed clusters involving different monomers. These results will help us to better understand the formation path of larger molecules and, providing reliable estimates of equilibrium geometries, can also be useful to guide and assist future spectroscopic works focused on these weakly bound complexes.