1 Introduction

Cycloaddition reaction is a pericyclic chemical reaction, in which two or more unsaturated molecules or parts of the same molecule combine with the formation of a cyclic adduct in which there is a net reduction of the bond multiplicity. Cycloadditions are usually described by the backbone size of the participants, which include the [4+2] cycloaddition and the [3+2] cycloaddition. The [3+2] cycloaddition is one of widely used organic chemical reaction, and has been extensively studied for the synthesis of 5-membered heterocycles^[1-2].

5-Amino-1H-tetrazole (5-AT), which was discovered by Thiele in 1893, is one of the most important intermediates for preparation of the high-energy density materials (HEDMs). In the past few years, a number of experiments have been carried out to explore the synthetical methods of 5-AT and derivatives thereof^[3-12]. The reaction of hydrazoic acid (HN₃) combing with the cyanamide (NH₂CN) is a useful method to synthize 5-AT, which is a typical [3+2] cycloaddition (Scheme 1).

To the best of our knowledge, less theoretical study on the title reaction has been reported ^[13-17]. The 1, 3-polar azide anion cyloaddition to nitriles was studied by Jursic, and the calculations agreed with the fact that electron-withdrawing substituents on the nitrile group decreased the activation barrier and facilitated the reaction ^[13].The transition states of both HCN+HN₃→H₂CN₄ and HCN+N₃^-→HCN₄^- reaction were investigated by Chen, and it was found that the latter reaction was more favored than the former one in view of the chemical kinetics and thermodynamics, tetrazole (H₂CN₄) and tetrazolate anion (HCN₄^-) were formed more easily in an alkali environment than in other systems ^[14]. The mechanisms of tetrazole formation by addition of azide to nitriles and why the reaction catalyzed by zinc(Ⅱ) salts were studied by Fahmi. The calculations indicated that coordination of the nitrile to the zinc ion was the dominant factor affecting the catalysis, this coordination substantially lowered the barrier for nucleophilic attack by azide ^[15-16]. The mechanism of the azide-nitrile cycloaddition mediated by the dialkylltin oxide-trimethylsilyl azide catalyst system and a new vilsmeier-haack type organocatalyst was investigated by David, and it was shown that as compared to the dialkylltin oxide-trimethylsilyl azide method, the organocatalytic system presented herein had the advantage of higher reactivity, in situ generation from inexpensive materials, and low toxicity ^[17]. In this present work, we attempt to perform an intensively theoretical calculation for the solvent effect of the title reaction using B3LYP and high-electron-correlation QCISD and MP2 methods with 6-311+G^* basis set.

2 Computational Methods

The geometries of all reactants (HN₃ and NH₂CN), transition state, and product (5-AT) have been fully optimized by B3LYP, QCISD, and MP2 methods at 6-311+G^* level. B3LYP is a DFT method using Becke′s three-parameter nonlocal exchange functional with the nonlocal correlation of Lee, Yang, and Parr ^[18-19]. QCISD is an extension of configuration interaction that corrects for size-consistency errors in the all singles and double excitation configuration interaction methods^[20-22]. MP2 stands for the second-order Moller Plesset perturbation theory ^[23]. 6-311+G^* is a split-valence triple-ζ polarization basis set augmented with diffuse functions^[23]. The structures and imaginary frequencies of transition states were confirmed by the vibration analysis and the intrinsic reaction coordinate (IRC) method at the same level. The changes of the bond lengths for the reaction as a function of the intrinsic reaction coordinate s were obtained at B3LYP/6-311+G^* level from -3.00 to 3.00 (amu)^1/2 bohr using the IRC method with a step size of 0.1 (amu)^1/2 bohr ^[24-26]. The single-point energies of all stagnation points were calculated at same levels. The statistical thermodynamic method [Eq. (1)] and Eyring transition state theory with Wigner correction [Eq. (2)] were used to determine the thermodynamic functions and rate constants of all reactions from 200 to 450 K.

$ \begin{array}{l} \Delta H = \sum {{H_{{\rm{product}}}}-\sum {{H_{{\rm{reactant}}}}} } \\ \Delta {H^ \ne } = {H_{{\rm{TS}}}}-\sum {{H_{{\rm{reactant}}}}} \\ \Delta S = \sum {{S_{{\rm{product}}}}-\sum {{S_{{\rm{reactant}}}}} } \\ \Delta {S^ \ne } = {S_{{\rm{TS}}}} - \sum {{S_{{\rm{reactant}}}}} \\ \Delta G = \Delta H - T\Delta S\\ \Delta {G^ \ne } = \Delta {H^ \ne } - T\Delta {S^ \ne }\\ {K^{\rm{ \mathsf{ θ} }}} = {e^{\left( { - \Delta G/RT} \right)}} \end{array} $

(1)

$ \begin{array}{l} k\left( T \right) = g\left( {{k_{\rm{b}}}T/h} \right)\exp \left( {{\Delta _{\rm{r}}}^ \ne {S_{\rm{m}}}/R-\Delta _{\rm{r}}^ \ne {H_{\rm{m}}}/RT} \right)\\ g = 1 + {\left( {h{v^ \ne }/{k_b}T} \right)^2}/24\\ A = g\left( {{k_{\rm{b}}}T/h} \right)\exp \left( {\Delta _{\rm{r}}^ \ne {S_{\rm{m}}}/R} \right) \end{array} $

(2)

where k(T), reaction rate constant. g, Wigner emendation factor. A, reaction frequency factor. k_b, Boltzmann factor, and h Planck′s constant. Δ^≠_rH_m and Δ^≠_rS_m are the mole activation enthalpy and mole activation entropy of a reaction system, and v^≠ is the imaginary frequency of a transition state.

To further explore the solvent effects, the above calculations were carried out using self-consistent reaction field (SCRF) approach with the polarizable continuum model (PCM) of the Tomasi′s group ^[27-32] in four solvents (carbon tetrachloride, dimethylsulfoxide, aceton, and water). All calculations are carried out using the Gaussian 09 software package ^[33].

3 Results and discussion 3.1 Geometries

The geometric parameters of the reactants (HN₃ and NH₂CN), products (5-AT), and TS optimized at B3LYP/6-311+G^* level of theory in the gas phase and four different solvents (carbon tetrachloride, acetone, dimethylsulfoxide, and water) are shown in Fig. 1 and Table 1. It is shown that the bond lengths participated in the reaction in TS increases in gas phase and in four solvents compared with those of the corresponding bonds in the reactants. The lengths of C1—N2 bond in TS are longer than that in HN₂CN. The other bond lengths, such as N6—N7 and N7—N8, have similar trend. All above bond lengths are further enlongated in 5-AT compared with those in TS. The distance of C1—N8 and N2—N6 have reversed trend in order to form the 5-membered heterocycle. It is also seen that the bond angles participated in the reaction in TS and product are decreased in comparison with those of the corresponding bonds in the reactants.

3.2 Reaction Path Properties

The variations of the bond lengths with the IRC are shown in Fig. 2. The trends of all bond lengths changing with the intrinsic reaction coordinate s are similar in the gas phase and in four solvents. The bond lengths of C1—N8 and N2—N6 change dramatically with s. In contrast, the bond lengths of all the other bonds change little. As the reaction proceeds from reactants to products, the length of the forming bond C1—N8 and N2—N6 remains minor changed until the intrinsic reaction coordinate reaches -0.1 (amu)^1/2 Bohr, where they start to fast decrease almost linearly with an increasing s, and the decreasing trend, however, slow after reaching 0.1 (amu)^1/2 Bohr. Synchronously, the elongating bond length, such as C1—N3, C1—N2, N6—N7, and N7—N8, increase slowly and almost linearly with s.

It is clear that C1 and N2 are in one molecule (NH₂CN), and N6 and N8 are in another molecule (HN₃) at the starting of cycloaddition. With the reaction proceeding, C1 (or N2) and N8 (or N6) are slowly closed with each other until the distance of C1 (or N2) and N8 (or N6) match the condition of forming a single bond. The elongating bond lengths, however, remain almost unchanged comparing with the forming bond during the reaction due to they are not broke but only changing the kind of bonds.

3.3 Energies

The zero-point energies (E_ZPE) and total energies (E) of all species are calculated by using B3LYP, QCISD, and MP2 method with 6-311+G^*level of theory (Table 2). The potential energies profile of the reaction (corrected by the zero-point energies, emendation factors were 0.96 ^[34]) are shown in Fig. 3. The activation energy barriers are increased with the solvent polarity. It indicates that the reaction become more difficult in polar solvent than in non-polar solvent.

3.4 Thermodynamics and Kinetics

The gas-phase thermodynamic and kinetic properties of the reaction are calculated from 200 K to 450 K with a step size of 50 K (Table 3). It is shown that ΔH and ΔS of the reaction decreases with increasing of temperature, and the equilibrium constant (K^Θ) is similar to them. ΔG of the reaction increase with temperature, and it will become positive when the temperature is higher than 450 K. It indicates that the reaction easily occurs at low temperature, and it will not carry out spontaneously when the temperature exceeds 450 K. According to Van′t Hoff equal pressure equation, the equilibrium constants of the reaction is increased with the decreasing of temperature because it is an exothermic reaction. So the reaction has advantage to be taken place at low temperature thermodynamically.

The reaction rate is remarkably increased analyzing from the kinetics, and the reaction is favored with increase of temperature. There is very little change in the reaction frequency factor (A), and A is considered to be a constant, so the reaction belongs to the Arrhenius-type reaction. Synthetically considering the thermodynamics and kinetics, the best reaction temperature is from 300 K to 350 K.

The thermodynamic and kinetic properties of the reaction in gas phase and solvents are shown in Table 4. All reactions are exothermic, entropy-decreasing, and spontaneous, and the equilibrium constants (K^Θ) are large. Analyzing from the thermodynamics, ΔH and ΔG of the reaction in dimethylsulfoxide solvent are the most negative, and K^Θ is the largest, so the reaction easier occurred spontaneously in dimethylsulfoxide solvent. Eeach rate constant in the solvents is smaller than that in gas phase, so the reaction is favored in gas phase kinetically.

The reaction rate in this calculation is smaller than the experimental data of the [3+2] reaction, because catalysts such as Rh, zinc bromide, aluminium chloride, lewis acid are used to accelerate the reaction process, and increase the reaction rate in experiment. The reaction energy barrier can be decreased about 6~16kcal/mol by catalysts, and the reaction rate can be increased about 10³~10^8[15-17].

4 Summary

In this work, we have analyzded the reaction paths, thermodynamic and kinetic properties and solvent effects for the cycloaddition reaction of hydrazoic acid with cyanamide. The B3LYP, QCISD and MP2 methods with 6-311+G^* level of theory are employed to optimize the geometries of all stationary points in gas phase and four solvents.

By analyzing MEP of the reaction, it is found that the regions of occurrence of the cycloaddition reaction in gas phase mainly occurs in the range of s=-0.1 to 0.1 (amu)^1/2 bohr on the MEP of the reaction. The trends of all bond lengths changing with s in the gas phase and four solvents are similar.

The equilibrium constant in dimethylsulfoxide solvent is the largest, so the reaction easier occurred spontaneously in dimethylsulfoxide solvent thermodynamically. Rate constant in the solvents is smaller than that in gas phase kinetically. In gas phase, the reaction is favored with increase of temperature, but has advantage to take place at low temperature thermodynamically, and 300 K to 350 K is the most feasible temperature to the reaction.