key: cord-0058667-p7jidfy6
authors: Valença Ferreira de Aragão, Emília; Faginas-Lago, Noelia; Rosi, Marzio; Mancini, Luca; Balucani, Nadia; Skouteris, Dimitrios
title: A Computational Study of the Reaction Cyanoacetylene and Cyano Radical Leading to 2-Butynedinitrile and Hydrogen Radical
date: 2020-08-20
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58808-3_51
sha: 97d7fea5deb432549c9a7acadbf631c980e5ceb3
doc_id: 58667
cord_uid: p7jidfy6

The present work focuses on the characterization of the reaction between cyanoacetylene and cyano radical by electronic structure calculations of the stationary points along the minimum energy path. One channel, leading to C[Formula: see text] N[Formula: see text] (2-Butynedinitrile) + H, was selected due to the importance of its products. Using different ab initio methods, a number of stationary points of the potential energy surface were characterized. The energy values of these minima were compared in order to weight the computational costs in relation to chemical accuracy. The results of this works suggests that B2PLYP (and B2PLYPD3) gave a better description of the saddle point geometry, while B3LYP works better for minima.

Cyanopolyynes are a family of carbon-chain molecules that have been detected in numerous objects of the Interstellar medium (ISM), such as hot cores, star forming regions and cold clouds [1] [2] [3] [4] . They are all linear molecules, with alternating carbon-carbon triple and single bonds. The simplest cyanopolyyne, HC 3 N, has been among the first organic molecules to be detected in the ISM [5] and up to date also HC 5 N, HC 7 N, HC 9 N and HC 11 N have been detected at least once in the ISM [6, 7] (the detection of HC 11 N, however, has been recently disputed by Loomis et al. [8] and Cordiner et al. [9] ). HC 3 N and HC 5 N are also abundant in solar-type protostars (see for instance a recent work on IRAS 16293-2422 by Jaber Al-Edhari et al. [10] ).

The shortest and most abundant member of the cyanopolyyne family, HC 3 N (cyanoacetylene), is a precursor of chain elongation reactions: the successive addition of C 2 H molecules generates the other members of its family, as summarised by Cheikh et al. [11] . According to the same authors, however, addition of CN radical instead of C 2 H would result in a chain termination reaction by the formation of dicyanopolyynes, NC−(CC) n -C N. These species have not been observed in the ISM so far because they lack a permanent electric dipole moment and cannot be detected through their rotational spectrum. However, it has been suggested that they are abundant in interstellar and circumstellar clouds [12] . The recent detection of NCCN [13] via its protonated form NCCNH + seems to corroborate the suggestion by Petrie et al. [12] .

In 2004, Petrie and Osamura had explored through computational means the formation of C 4 N 2 (2-Butynedinitrile) in Titan's atmosphere through many pathways [14] . They had found in particular that the cyano radical addition to cyanoacetylene leads to C 4 N 2 . In order to characterize all the stationary points of the potential energy surface, the authors had used at the time the hybrid density functional method B3LYP in conjunction with triple-split valence Gaussian basis set 6-311G** for geometry optimizations and vibrational frequency calculations. Moreover, they had also performed single-point calculations with CCSD(T) in conjunction with aug-cc-pVDZ basis-set, a choice made at the time due to computational costs. As computational power has risen in the last 15 years and new methods have been implemented into quantum chemistry software, more accurate results for the geometries and energies can be obtained.

In our laboratory we have already investigated several reactions of astrochemical interest providing insightful results for the understanding of processes observed in the ISM [15, 16] including those involving or leading N-bearing organic molecules [17] [18] [19] . Recently, we have focused on studying the reaction between HC 3 N and CN in collaboration with the experimental part of the Perugia group. The preliminary investigation of the potential energy surface of the system HC 3 N + CN showed the presence of a number of reactive channels. The search for intermediate and product species always involves the computation of four carbon atoms, two nitrogen atoms and one hydrogen atom, i.e. 39 electrons in total. Ab initio calculations can become prohibitively expensive with a rising number of electrons, therefore a balance between chemical accuracy and computational cost must be reached. The focus of this paper is the exit channel of the reaction between cyanoacetylene and cyano radical that leads to the formation of 2-Butynedinitrile and hydrogen. Four different computational methods were benchmarked in order to check if relatively cheap methods can get accurate results for this system. A comparison to Petrie and Osamura's results was also done.

The paper is organized as follows: in Sect. 2, the methods and the construction of the potential energy surface are outlined. Preliminary results are reported in Sect. 3 and in Sect. 4 concluding remarks are given.

The Potential Energy Surface (PES) of the system was investigated through the optimization of the most stable stationary points. The overall multiplicity of the system was 2 (doublet) and calculations were performed adopting an unrestricted formalism. Following a well established computational scheme [20] [21] [22] [23] [24] [25] [26] , we optimized the geometry of the stationary points, both minima and saddle points, using a less expensive method with respect to the one employed in order to get accurate energies. Calculations for geometries were performed in order to benchmark three methods: Unrestricted-Hartree-Fock (UHF) [27, 28] , density functional theory (DFT), with the Becke-3-parameter exchange and Lee-Yang-Parr correlation (B3LYP) [29, 30] , and double-hybrid DFT method B2PLYP [31] combined or not with Grimme's D3BJ dispersion [32, 33] . All methods were used in conjunction with the correlation consistent valence polarized basis set aug-cc-pVTZ [34] . In each level of theory, a vibrational frequency analysis was done to confirm the nature of the stationary points: a minimum in the absence of imaginary frequencies and a saddle point if one and only one frequency is imaginary. The assignment of the saddle points was performed using intrinsic reaction coordinate (IRC) calculations [35, 36] . Then for each stationary point for all methods, the energy was computed with coupled cluster including single and double excitations and a perturbative estimate of connected triples (CCSD(T)) [37] [38] [39] . CCSD(T) is a more accurate method than the ones used for the optimizations, but prohibitive computational costs restricts the use of the method in this particular system to fixed geometry calculations. Finally, zero-point energy correction obtained through the frequency calculations were added to energies obtained from all methods to correct them to 0 K. All calculations were performed using the Gaussian 09 code [40] .

The calculation of the different electronic structures shows that the attack of cyano radical on the cyanoacetylene is an energetically favorable process that leads to the formation of an adduct intermediate. Figure 1 gather the geometries at the minimum energy path, starting from the HC 3 N and CN reactants and leading to 2-Butynedinitrile and hydrogen. These geometries were optimized at UHF/aug-cc-pVTZ, B3LYP/aug-cc-pVTZ, B2PLYP/aug-cc-pVTZ and B2PLYPD3/aug-cc-pVTZ levels. The reported interatomic distances are in red for UHF, blue for B3LYP and green for both B2PLYP and B2PYLP with Grimme's D3 dispersion, since no difference in their geometries was recorded. The geometries obtained with the methods above are very similar. The differences between bond lengths and angles are small, with the exception of the transition state geometries. In the saddle point structure, B2PLYP seems to provide a more reasonable distance for weak interaction as the one between carbon and hydrogen: it is 0.162Å shorter than the distance obtained with B3LYP and 0.241Å shorter than the one proposed by UHF. B2PLYP is a double-hybrid density functional that combines Becke exchange and Lee, Young and Parr correlation with Hartree-Fock exchange and a perturbative second-order correlation obtained from Kohn-Sham orbitals [31] . According to the author of the method, B2PLYP reports very good results for transition state geometries.

Most of the geometries might be similar, but this does not necessarily means that all energies also are. Table 1 gathers the relative energies of all geometries computed with every method. The reactants were taken as the reference for the energy. Within the B2PLYP method, the account of Grimme's D3 dispersion changes slightly the energy value, even for an identical geometry. As every method returns different absolute energy values, a direct comparison is not pertinent. However, the trend in the evolution of the relative energies can be compared. First, it can be observed that all stationary points are below the energy of the reactants in all the methods that have been employed. Table 1 . Energies (kJ.mol −1 , 0 K) of the different geometries relative to the reactants. The energies are computed at UHF/aug-cc-pVTZ, B3LYP/aug-cc-pVTZ, B2PLYP/aug-cc-pVTZ, and B2PLYPD3/aug-cc-pVTZ levels of theory. In parentheses are the values for the same geometry computed at CCSD(T)/aug-cc-pVTZ level of theory.

Second, in all methods the adduct intermediate MIN 1 is the most energetically stable structure. In addition, a barrier between the adduct and the products is well characterized for B3LYP and both B2PLYP methods, but it is not the case for the UHF method. Though the saddle point was identified with UHF, the energy of the products returned a higher value. It will be shown, however, that when a single-point energy calculation with these same geometries optimized with UHF is done with the CCSD(T)/aug-cc-pVTZ method, the energy of the products is below the energy of the saddle-point.

In order to compare the geometries optimized with all those methods, singlepoint CCSD(T) calculations were also performed. The results are reported again in Table 1 , in parenthesis. Firstly, it can be observed that the energies of all stationary points, i.e. minima and transition state, are below the energy of the reactants in all the methods used. Secondly, as B2PLYP and B2PLYPD3 optimized geometries were mostly identical, the energies computed with CCSD(T)/augcc-pVTZ were also the same. Thirdly, it can be noticed that the energies of the B3LYP optimized geometries are very close to the ones of B2PLYP methods, with the exception of the transition state geometry: for this saddle point, energy difference is around 19.6 kJ.mol −1 . At last, the energies of the UHF optimized geometries are the most different from the others methods. In particular, the energies of the intermediate and transition state gave the lowest energy at CCSD(T) level. In contrast, the energy of the products was the highest one using the UHF geometry. Nevertheless, the energy barrier between the intermediate and the products is now well characterized for all methods.

It is also interesting to look at the barrier height values obtained for each method. Table 2 reports the energy changes and barrier heights, computed at the same levels of theory, for the process leading to 2-Butynedinitrile and hydrogen. The interaction of HC 3 N and CN gives rise to the adduct MIN 1 (or H(CN)C 3 N), more stable than the reactants in all levels of theory. This adduct evolves, through a barrier leading to the transition state, to the products of the reactive channel C 4 N 2 and radical H. In all levels of theory, the products are less stable than the adduct. Making a comparison between methods, UHF estimates the largest enthalpy changes and barrier height. UHF is followed by B3LYP, B2PLYPD3 and B2PLYP in this order. In relation to CCSD(T), a higher-level method, UHF and B3LYP energies are systematically overestimated. B2PLYP and B2PLYPD3 overestimates the CCSD(T) in the enthalpy variation attributed to the formation of the adduct, but underestimate both the barrier and enthalpy change for the reaction that leads the intermediate to the product. Table 2 . Enthalpy changes and barrier height (kJ.mol −1 , 0 K) computed at UHF/augcc-pVTZ, B3LYP/aug-cc-pVTZ, B2PLYP/aug-cc-pVTZ, and B2PLYPD3/aug-cc-pVTZ levels of theory. In parentheses are the values for the same geometry computed at CCSD(T)/aug-cc-pVTZ level of theory.

While in this work the aug-cc-pVTZ basis set was employed for every method, Petrie and Osamura [14] carried geometry optimizations at B3LYP level in conjunction with the 6-311G** basis-set and computed the single-point energies at CCSD(T) level in conjunction with aug-cc-pVDZ basis-set. In respect to the geometries obtained, only the transition state and the intermediate geometries were published in their paper. The bond distance and angle values are similar to the ones in the B3LYP/aug-cc-pVTZ optimized geometries. The difference is smaller than when the B3LYP geometries were compared to the ones obtained with other methods. The energies computed at CCSD(T)/aug-cc-pVDZ level are listed in Table 2 Table 2 .

UHF, DFT (with B3LYP functional) and double-hybrid DFT (with B2PLYP and B2PLYPD3 functionals) were benchmarked on four points of the PES of the astrochemically relevant reaction between cyanoacetylene (HC 3 N) and cyano radical. As far as optimized geometries for stationary points are concerned, the UHF method seems to be inadequate, while DFT methods seem to be more reliable. In particular, B3LYP functional seems to work well for minima, while better functionals like the B2PLYP seem to be necessary for transition state geometries when van der Waals interactions are present. As far as energies are concerned, however, only more correlated methods like CCSD(T) seems to provide accurate results.

A more general conclusion concerns the reaction mechanism suggested by our calculations, which is also in line with the previous determination by Petrie and Osamura. Similarly to the case of other reactions involving CN and other species holding a triple C−C bond (such as ethyne, propyne and 2-butyne [41] [42] [43] ), the CN radical interacts with the electron density of triple bond to form an addition intermediate without an activation energy. This is in line with the large value of the rate coefficient derived for the reactions also at very low temperature [11] and makes this process a feasible route of dicyano-acetylene (2-Butynedinitrile) even under the harsh conditions of the interstellar medium.

After completing the derivation of the potential energy surface for the title reaction, we will run kinetic calculations to derive the rate coefficient and product branching ratio. The calculated rate coefficient will be compared with the experimental values derived by Cheikh et al. [11] while the reaction mechanism and product branching ratio will be compared with those inferred by the crossed molecular beam experiments which are now in progress in our laboratory. A thorough characterization of the title reaction will allow us to establish its role in the nitrogen chemistry of the interstellar medium. To be noted that both cyano-and dicyano-acetylene have a strong prebiotic potential [44] .

Vibrationally excited HC3N toward hot cores

Survey observations to study chemical evolution from high-mass starless cores to high-mass protostellar objects I: HC3N and HC5N

A search for cyanopolyynes in L1157-B1

Observations of 13 C isotopomers of HC3N and HC5N in TMC-1: evidence for isotopic fractionation

Detection of interstellar cyanoacetylene

The detection of HC9N in interstellar space

Detection of HC11N in the cold dust cloud TMC-1

Non-detection of HC11N towards TMC-1: constraining the chemistry of large carbon-chain molecules

Deep Kband observations of TMC-1 with the green bank telescope: detection of HC7O, nondetection of HC11N, and a search for new organic molecules

History of the solar-type protostar IRAS 16293-2422 as told by the cyanopolyynes

Low temperature rate coefficients for the reaction CN + HC3N

NCCN in TMC-1 and IRC+ 10216

Probing non-polar interstellar molecules through their protonated form: detection of protonated cyanogen (NCCNH + )

NCCN and NCCCCN formation in Titan's atmosphere: 2. HNC as a viable precursor

Silicon-bearing molecules in the shock L1157-B1: first detection of SiS around a Sun-like protostar

Interstellar dimethyl ether gas-phase formation: a quantum chemistry and kinetics study

Dimerization of methanimine and its charged species in the atmosphere of Titan and interstellar/cometary ice analogs

A theoretical investigation of the reaction between the amidogen, NH, and the ethyl, C2H5, radicals: a possible gas-phase formation route of interstellar and planetary ethanimine

Low temperature kinetics and theoretical studies of the reaction CN + CH3NH2: a potential source of cyanamide and methyl cyanamide in the interstellar medium

Stereoselectivity in autoionization reactions of hydrogenated molecules by metastable noble gas atoms: the role of electronic couplings

Crossed-beam and theoretical studies of the S( 1 D) + C2H2 reaction

The intermolecular potential in NO-N2 and (NO-N2) + systems: implications for the neutralization of ionic molecular aggregates

The proton affinity and gas-phase basicity of sulfur dioxide

Observation of organosulfur products (thiovinoxy, thioketene and thioformyl) in crossed-beam experiments and low temperature rate coefficients for the reaction S( 1 D) + C2H4

SSOH and HSSO radicals: an experimental and theoretical study of [S2OH] 0/+/− species

Theoretical study of reactions relevant for atmospheric models of Titan: interaction of excited nitrogen atoms with small hydrocarbons

New developments in molecular orbital theory

Self-consistent orbitals for radicals

Density functional thermochemistry. III. The role of exact exchange

Ab Initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields

Semiempirical hybrid density functional with perturbative secondorder correlation

Effect of the damping function in dispersion corrected density functional theory

Efficient and accurate double-hybrid-meta-GGA density functionals-evaluation with the extended GMTKN30 database for general main group thermochemistry, kinetics, and noncovalent interactions

Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen

An improved algorithm for reaction path following

Reaction path following in mass-weighted internal coordinates

Many-body perturbation theory and coupled cluster theory for electron correlation in molecules

A fifth-order perturbation comparison of electron correlation theories

Full configurationinteraction and state of the art correlation calculations on water in a valence double-zeta basis with polarization functions

Gaussian 09, Revision A. 02

Crossed beam reaction of the cyano radical, CN (X 2 Σ + ), with methylacetylene, CH3CCH (X 1 A1): observation of cyanopropyne, CH3CCCN (X 1 A1), and cyanoallene, H2CCCHCN (X 1 A )

Crossed beam reaction of cyano radicals with hydrocarbon molecules. II. Chemical dynamics of 1-cyano-1-methylallene (CNCH3CCCH2

A ) formation from reaction of CN (X 2 Σ + ) with dimethylacetylene CH3CCCH3 (X 1 A 1 )

Crossed beam reaction of cyano radicals with hydrocarbon molecules. IV. Chemical dynamics of cyanoacetylene (HCCCN; X 1 Σ + ) formation from reaction of CN (X 2 Σ + ) with acetylene, C2H2 (X 1 Σ + g )

Elementary reactions of N atoms with hydrocarbons: first steps towards the formation of prebiotic N-containing molecules in planetary atmospheres