Ab initio and density functional theory study of the diazene isomerization

11 October 1996 .

.

.

:

.

CHEMICAL PHYSICS LEI"rERS

<,:

ELSEVIER

Chemical Physics Letters 261 (1996) 13-17

Ab initio and density functional theory study of the diazene isomerization Branko S. Jursic Department of Chemistry, University of New Orleans, New Orleans, LA 70148, USA Received 27 February 1996; in f'mal form 12 July 1996

Abstract The HF, MP2, QCISD(T), G1, G2, and G2MP2 ab initio and twenty-two density functional theory (DFT) methods were used to perform the computational studies of the syn- and cis-diazene transformation into trans-diazene. In all calculations the Gaussian type 6-311 + G(2d,2p) basis set was used. The DFT computed geometries and energies are compared with ab initio results and the validity of several DFI" methods for studying this and similar molecular systems are evaluated.

1. Introduction There is considerable interest in the chemical community to establish reliable computational methods for the study of reaction mechanisms, especially for systems where experimental results cannot be easily obtained. In these cases, one usually combines the optimization of transition state structures and the estimation of postulated reaction barriers for determining properties of interest. Many computational methods correctly predict the geometries of the reactants and products of studied reactions but these methods become problematic when studying systems where bonds are broken or formed [1,2]. This is even more profound for molecules with heteroatoms [1,2]. Transition state structures are good examples of problematic chemical systems that do not have fully formed bonds and, thus, they require computational methods that include electron correlations [3]. The ab initio methods which use high level correlation interactions are computationally expensive and can be practically employed for only small chemical systems. On the other hand, Db-T methods have electron

correlations and require only a portion of the time that is necessary for ab initio optimization. It has become well known that the physical and chemical properties of molecules can be predicted through electron density. This is now known as the Hohenberg-Kohn theorem [4]. It becomes computationally applicable with a Kohn-Sham [5] type theory, which makes use of orbitals for the representation of electron density. Many density functional theory (DFF) studies demonstrate the usefulness of this approach [6,7]. We have successfully employed DFT methods for computing the geometries [8-10] and activation barriers [11-13] for many different chemical systems. Here we are presenting our computational study of the diazene isomerization by ab initio and DFT methods.

2. Computational methods All calculations were performed with the GAUSSIAN 94 computational package [14]. Tradi-

0009-2614/96/$12.00 Copyright © 1996 Elsevier Science B.V. All fights reserved PII S0009-261 4(96)00887- l

14

B.S. Jursic / Chemical Physics Letters 261 (1996) 13-17

tional ab initio methods were used, including Hartree-Fock (HF) [15], the second, third, and fourth order of the M011er-Plesset energy correction (MPn) [16,17], and the quadratic configuration interaction with singly and doubly substituted configurations and with the incorporation of contributions of triply substituted configurations (QCISD(T)) [ 18]. G 1, G2, and G2MP2 high level ab initio computational studies were also performed. These methods require GAUSSIAN-1 (G1) and GAUSSIAN-2 (G2) [14] for computing very accurate energies. The G2MP2 is a modified version of G2 which uses MP2 instead of MP4 for the basis set extension correlation. It is nearly as accurate as the full G2. The full configuration interaction (CI) is the most complete nonrelativistic treatment of the molecular system possible, within the limitations imposed by the chosen basis set [20]. It represents a possible quantum state of the system while modeling the electron density in accordance with the definition of the basis set in use. Practical configuration interaction methods augment the Hartree-Fock method by adding only a limited set of substitutions, truncating the CI expansion at the same level of substitution. For example, the CISD method adds single and double excitations to the HF determinant. The disadvantage of this limited CI variation is that it is not size-consistent. The QCISD method adds terms to CISD to restore size consistency. On the other hand the G2 procedure is known to generate highly accurate energies [19] and in many of the computational studies it is used as a reference together with experimental data Three different groups of DFF calculations were employed: exchange functional, hybrid, and exchange functional in combination with correlation functional methods. These are the Slater exchange functional p4/3 with a theoretical coefficient of 2 / 3 commonly referred to as the local spin density exchange (HFS) [21], the Xot exchange functional p4/3 with an empirical coefficient of 0.7 (XALPHA) [21], and Becke's 1988 functional [22] which includes the Slater exchange along with correlations involving the density gradient. Three hybrid DFF methods were applied including Becke's three parameter functional [23] which has the f o r m AExSlater -I- ( l - a)Enx F q- B A E x Becke qE~VWN+ CA E~°nl°ca~, where the nonlocal correlation is provided by the LYP [24] expression (B3LYP).

The constants A, B, and C are those determined by Becke by fitting to the G1 molecular set. Becke's three parameter functional [23], with the nonlocal correlations provided by the Perdew 86 expression [25] (B3P86) and by the Perdew/Wang 91 expression [26] (B3PW91), was also used. A combination of exchange and correlation functional DFT methods was also used. Three exchange functionals, Slater (S) [21], Xet (XA) [21], and Becke's 88 (B) [22], were combined with Vosko, Wilk, and Nusair's 1980 (VWN) [27], LYP [24], P86 [25], PW91 [26], PL [28], and VWN5 [27] correlation functionals. For all calculations an extended Gaussian type basis set, 6-311 + G(2d,2p), was employed, and in one instance a variation of the Gaussian basis sets was used. Explanation and abbreviation of the basis sets are found in the GAUSSIAN computational package [20].

3. Results and discussion

Although only the trans- and cis-diazenes have been observed experimentally, the other three structures corresponding to the two minima and two transition structures on the potential energy surface have been postulated in several theoretical investigations [29-31]. There are photochemical studies of organic derivatives of l,l-diazene [32]. Several theoretical groups have also calculated the order of the electronic states of 1,l-diazene [33,34]. Goddard III and Davis [34] used HF, generalized valence bond (GVB), and configuration interaction (GVB-CI) calculated geometries and energies of both singlet and triplet 1,1-diazenes. Here we are presenting DFF and ab initio computational studies of the diazene isomerization. The computed transition structures, TS1 and TS2, representing the isomerization of cis- and syn-diazene into trans-diazene are presented in Table 1. By definition, HF does not include an electron correlation and does not perform well for systems that require an electron correlation such as the one studied here. The computed bond distances are considerably shorter, which is a common problem with HF

B.S. Jursic / Chemical Physics Letters 261 (1996) 13-17

calculations [1,2]. The MP2/6-311 + G(2d,2p) theory model predicts longer bond distances while CISD geometries are somewhere between the HF and MP2 computed parameters (Table 1). On the other hand, QCISD(T) ab initio methods predict the longest bond distances or the most loose transition states of all of the methods used in this study. The multi-configurational self-consistent field, MCSCF, computational study performed by Helgaker and coworkers [30]

~ ( ~ TS1

15

computes a considerably longer N-N bond distance than does MP2. Contrary to this, the limited CI (or CISD) computational study produces geometries of these two transition state structures that are somewhere between the HF and MP2 computed transition state structures (Table 1). Considering that QCISD(T)/6-311 + CA2d,2p) is the highest level ab initio computational method, the computed structures should be the most accurate.

~TS2

Table 1 The geometric parameters of transition state structures for the transformation of cis- and syn-diazene into trans-diazene as computed by using 6-311 + g(2d,2p) basis set Transition state structure TSI a

Transition state structure TS2 b

Method

r21 ( ~ )

r31 (,~)

r42 ( ~ )

a312 (o)

a421 (o)

r21 (m)

r32 (,~)

r43 (,~k)

a321 (o)

a234 (o)

HF MP2 MCSSCF CI QCISD(T) B3LYP B3P86 B3PW91 HFS HFB Xct BPL BP86 BPW91 BLYP BVWN BVWN5 XAPL XAP86 XAPW91 XALYP XAVWN XAVWN5 SPL SP86 SPW91 SLYP SVWN

1.201 1.237 1.250 1.221 1.243 1.223 1.219 1.220 1.232 1.250 1.219 1.236 1.233 1.231 1.236 1.234 1.236 1.207 1.204 1.203 1.206 1.205 1.205 1.220 1.217 1.216 1.219 1.218

1.024 1.046 1.057 11037 1.062 1.055 1.055 1.055 1.096 1.089 1.08 1.069 1.074 1.072 1.073 1.066 1.068 1.060 1.065 1.063 1.063 1.057 1.059 1.076 1,082 1.079 1.080 1.073

0.979 0.992 0.990 0,985 0.995 0.993 0.992 0.992 1.017 1.009 1.005 0.996 1.001 0.999 1.000 0.994 0.996 0.993 0.997 0.995 0.996 0.991 0,992 1.004 1.009 1.007 1.008 1.002

110.5 109.4 103.4 109.8 109.3 110.2 110.1 110.1 110.5 110.0 110.6 110.0 110.0 110.0 110.1 110.0 110.0 10.7 10.5 10.5 10.7 10.6 10.7 10.5 10.4 10.4 10.5 10.5

179.4 178.6 178.9 179.2 178.1 178.5 178.5 178.5 187.8 187.6 188.1 187.9 187.9 187.9 188.0 188.0 188.0 188.4 188.3 188.4 188.4 188.4 188.4 188.1 188.1 188.1 188.1 188.2

1.016 1.043 1.066 1.030 1.055 1.052 1.051 1.052 1,097 1.088 1.079 1.068 1,074 1,071 1,072 1,065 1.067 1,060 1,066 1,063 1,064 1.057 1.059 1.076 1.082 1.080 1.081 1.073

1.246 1.284 1.299 1.268 1.288 1.273 1.268 1.269 1,284 1.207 1.268 1.291 1.285 1.284 1.290 1.289 1.290 1.254 1.250 1.248 1.253 1.252 1.253 1.269 1.264 1.263 1.268 1.267

1.284 1.362 1.367 1.319 1.374 1.334 1.333 1.334 1.378 1.376 1.357 1.253 1.259 1.356 1.358 1.349 1.352 1.337 1.345 1.342 1.342 1.334 1.336 1.356 1.364 1,361 1.361 1.353

120.5 122.1 121.0 121.0 122.1 120.9 120.9 121.0 121.2 120.8 121.3 120.9 121.0 121.1 120.9 120.9 120.9 121.3 121.4 121.4 121.3 121.3 121.3 121.3 121.4 121.4 121.3 121.2

52.3 49.1 50.9 50.9 48.6 50.9 50.7 50.7 50.7 50.4 50.7 50.4 50.3 50.2 50.5 50.4 50.4 50.7 50.5 50.5 50.8 50.7 50.7 50.6 50.4 50.4 50.7 50.6

16

B.S. Jursic / Chemical Physics Letters 261 (1996) 13-17

Therefore, we will use these as a reference. All calculations agree that the NNH bond angle for the cis-trans isomerization is nearly 180 ° and the system is planar with the other HNN bond angle around 110 ° (Table 1). Major disagreement is observed when the N - N bond distances are compared. The closest values to the QCISD(T) ab initio results were obtained with Becke's 88 exchange functional based DFT methods and with hybrid DFT methods. The computed values are in the range of 0.01 ~,. A major disagreement was obtained with Siater's and Xet exchange functional based DFT methods. A similar behavior was observed for the syn-trans isomerization transition state (Table 1). Again, the best agreement is between geometries computed with the QCISD(T)/6-311 + G(2d,2p) ab initio theory model and with Becke's exchange functional based DF-I" methods. In computational studies, it is appropriate to compare the energies of identical chemical systems computed with different methods as a measure of their accuracy. The total energies for the two isomerization reactions are presented in Table 2. For the isomerization reactions, neither the reaction barriers nor the heats of reaction are known experimentally because both syn- and cis-diazene have not been observed. The transition state TS 1 represents a bond distortion and not a bond breakage. One can predict that a low level ab initio calculation should not estimate an acceptable activation barrier. The HF computed barrier is 49.7 kcal/mol. It is well known that the HF ab initio method overestimates reaction barriers. Usually MP2 underestimates the same bartier, but in the case of the cis-trans isomerization there is only a slight difference. To our disappointment, MCSCF generates an activation barrier through TS1 that is higher than the one computed by HF (Table 2). This is also the case for QCISD(T). This is not a usual pattern and we have no explanation for this behavior. Interestingly, all of the ab initio studies, G1, G2, and G2MP2, compute almost identical reaction barriers. We believe that the most accurate value is 44.6 kcal/mol which was computed with G2MP2 [19]. The hybrid methods compute identical values while Becke's 88 exchange functional based DFT methods generate a reaction barriers that differ by ~ 1 kcal/mol. Two other groups of DFI" methods, Slater's and Xot exchange functional based

Table 2 The barriers of the diazene isomerization (kcal/mol) and relative energies (kcal/mol) of cis and syn-diazene in regard to trans-diazene Method

AEI

AEH

AEm

A Ely

HF MP2 CI MCSSCF QCISD(T) GI(0 K) G1 G2(0 K) G2 G2MP2(0 K) G2MP2 B3LYP B3P86 B3PW91 HFS HFB Xa BPL BP86 BPW91 BLYP BVWN BVWN5 XAPL XAP86 XAPW91 XALYP XAVWN XAVWN5 SPL SP86 SPW91 SLYP SVWN

49.7 48.1 48.6 60.0 58.7 44.5 44.6 44.7 44.7 44.6 44.6 44.4 44.6 44.6 40.9 44.0 40.9 44.2 43.7 43.7 43.4 44.2 44.2 41.1 40.8 40.8 40.5 41.2 41.1 41.1 40.8 40.8 40.4 41.1

69.0 50.0 58.3 51.3 50.9 45.3 45.2 45.6 45.5 45.2 45.1 53.8 52.7 52.8 45.2 49.8 46.6 51.4 48.5 48.8 49.9 51.7 51.5 47.9 45.1 45.4 46.6 48.1 48.0 46.6 44.1 44.1 45.2 46.8

6.6 6.1 6.1 7.0 6.2 5.5 5.5 5. I 5.1 5.0 5.0 5.6 5.7 5.7 5.3 4.9 5.5 5.2 5.2 5.2 5.2 5.2 5.2 5.8 5.9 5.9 5.8 5.8 5.8 5.6 5.7 5.7 5.6 5.6

19.7 28.5 24.4 35.0 26.4 24.0 24.0 24.1 24.1 24.3 24.4 22.0 22.3 22.3 21.0 21.9 21.0 22.1 22.2 22.2 21.9 22.1 22.1 21.4 21.4 21.5 21.2 21.4 21.4 21.3 21.3 21.4 21.1 21.4

A Et = barrier for transformation of cis-diazene to trans-diazene through TS 1; A Eli = barrier for transformation of syn-diazene to trans-diazene through TS2; AE m = relative energy of cis-diazene; A E l v = relative energy of syn-diazene.

DFT methods underestimate the reaction barrier by 3 - 4 kcal/mol. In the transition state structure TS2, breaking the N - H bond and forming the H - N bond is occurring concertedly. Thus, an electron correlation should be more important than in the case of TS1. This is clearly demonstrated in going from 69 kcal/mol, computed with H F / 6 - 3 1 1 + G(2d,2p), to 50 kcal/mol, computed by the MP2/6-311 + G(2d,2p) theory model (Table 2). Because better results with

B.S. Jursic / Chemical Physics Letters 261 (1996) 13-17

QCISD(T) were not obtained, one would conclude that saturation was achieved. However, G1, G2, and G2MP2 ab initio methods compute the reaction barrier to be around 45 k c a l / m o l . The hybrid DFT methods estimate a reaction barrier ( ~ 53 k c a l / m o i ) that is much closer to QCISD(T) ( ~ 51 k c a l / m o l ) than to the G2MP2 (45.1 k c a l / m o l ) computed value. Here Becke's 88 exchange functional overestimates the barrier by around 3 k c a l / m o l . To our surprise the best values were obtained by Slater's exchange functional (Table 2). The enthalpies of the isomerization follow a similar trend. The cis-trans isomerization enthalpy computed with HF, MP2, and QCISD(T) are overestimated by 1 - 2 k c a l / m o l when compared with the G2 and G2MP2 values (Table 2). The best values were computed with Becke's 88 exchange functional based DFT methods, which differ by 0.2 k c a l / m o l from the G2 values. In the case of the s y n - t r a n s isomerization, the computed values with both hybrid and Becke's 88 exchange functional DFT methods differ by ~ 2 k c a l / m o l (Table 2).

4. C o n c l u s i o n Thus, the activation barrier for the cis-trans isomerization computed with the G2MP2 ab initio method is 44.6 k c a l / m o l with the enthalpy of the reaction at 5.0 k c a l / m o l . The reaction barrier for the sy n- t r a n s isomerization is 45.2 k c a l / m o l and the enthalpy of the reaction is 24.2 k c a l / m o l . The closest computed values to the G2 and G2MP2 energies are obtained with BPW91 and BLYP, which are nonlocal methods, followed by the B3PW91 and B3LYP hybrid DFT methods. In the best case the energies differ by less than 3 k c a l / m o l . Considering our previous investigations of reaction mechanisms with DFI" methods, we suggest these four D F T methods are reliable for the computational study of isomerization reactions of these and similar molecular systems.

References [I] W.J. Hehre, L. Radom, P. yon R. Schleyer and J.A. Pople, Ab initio molecularorbital theory (Wiley, New York, 1986). [2] W.J. Hehre, Experimentsin computationalorganic chemistry (Wavefunction,Inc., lrvine, CA, 1995).

17

[3] T.H. Dunning,ed., Advancein molecularelectronicstructure theory (JAI Press, Greenwich, CT, 1990); C. Gonzalez and H.B. Schlegel, J. Phys. Chem. 94 (1990) 5523. [4] P. Kohn and L.J. Sham, Phys. Rev. B 136 (1964) 864. [5] W. Kohn and L.J. Sham, Phys. Rev. A 140 (1965) 1133. [6] R.G. Parr and W. Yang, Density functional theory of atoms and molecules(Oxford Univ. Press, New York, 1989). [7] J.M. Seminario and P. Politzer, eds., Modern density functional theory a tool for chemistry (Elsevier, Amsterdam, 1995). [8] B.S. Jursic, Int. J. QuantumChem. 57 (1996). [9] B.S. Jursic, Chem. Phys. Lett. 236 (1995) 206, [10] B.S. Jursic, J. Mol. Struct. THEOCHEM 358 (1995) 145; Int. J. QuantumChem. 57 (1996) 213. [11] B.S. Jursic and Z. Zdravkovski, J. Chem. Soc. Perkin 2 (1995) 1223. [12] B.S. Jursic, Chem. Phys. Lett. 244 (1995) 163. [13] B.S. Jursic, J. Chem. Phys. 104(1996) 4151. [14] M.J. Frisch, G.W. Trucks, H.B. Schlegel, P.M.W. Gill, B.G. Johnson, M.A. Robb, J.R. Cheeseman,T. Keith, G.A. Petersson, J.A. Montgomery, K. Raghavachari, M.A. A1-Laham, V.G. Zakrzewski,J.V. Ortiz, J.B. Foresman,C.Y. Peng, P.Y. Ayala, W. Chen, M.W. Wong, J.L. Andres,E.S. Replogle,R. Gomperts, R.L. Martin, D.J. Fox, J.S. Binkley, D.J. Defrees, J. Baker, J.P. Stewart, M. Head-Gordon, C. Gonzalez, and J.A. Pople, GAUSSIAN 94, Revision B.3 (Gaussian, Inc., Pittsburgh, PA, 1995). [15] R. McWeeny and G. Dierksen, J. Chem. Phys. 49 (1968) 4852. [16] C. M~ller and M.S. Plesset, Phys. Rev. 46 (1934) 618. [17] S. Saebo and J. Almlof, Chem. Phys. Lett. 154 (1989) 83. [18] J.A. Pople, M. Head-Gordonand K. Raghavachari, J. Chem. Phys. 87 (1987) 5968. [19] L.A. Curtiss, K. Raghavachari and J.A. Pople, J. Chem. Phys. 98 (1993) 1293. [20] M.J. Frisch, A. Frisch and J.B. Foresman, GAUSSIAN 94 User's Reference (Gaussian,Inc. Pittsburgh,PA, 1995). [21] J.C. Slater, Quantumtheory of molecularand solids, Vol. 4. The self-consistentfield for molecular and solids (McGrawHill, New York, 1974). [22] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [23] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [24] C. Lee, W. Yang and R.G. Parr, Phys. Rev. B 37 (1988) 785. [25] J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [26] J.P. Perdew and Y. Wang, Phys. Rev. B 45 (1992) 13244. [27] S.H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58 (1980) 1200. [28] J.P. Perdew and A. Zunger, Phys. Rev. B 23 (1981) 5048. [29] L. Fan and T. Ziegler, J. Chem. Phys. 92 (1990) 3645. [30] H.J.A. Jensen, P. JCrgensenand T. Helgaker, J. Am. Chem. Soc. 109 (1987) 2895. [31] C.A. Parson and C.E. Dykstra, J. Chem. Phys. 71 (1979) 3025. [32] P.G. Schultzand P.B. Dervan,J. Am. Chem. Soc. 104 (1982) 6660. [33] G. Wagniere,Theor. Chim. Acta 31 (1973) 269. [34] J.H. Davis and W.A. Goddard III, J. Am. Chem. Soc. 99 (1977) 7111.