Ab initio calculations of the rotational barriers in formamide and acetamide: The effects of polarization functions and correlation

Ab initio calculations of the rotational barriers in formamide and acetamide: The effects of polarization functions and correlation

Journal of Molecular Structure (Theochem), 139 (1986) 197-206 Elsevier Science Publishers B.V., Amsterdam -Printed in The Netherlands AB INITIO CALiX...

730KB Sizes 0 Downloads 54 Views

Journal of Molecular Structure (Theochem), 139 (1986) 197-206 Elsevier Science Publishers B.V., Amsterdam -Printed in The Netherlands

AB INITIO CALiXJLATIONS OF THE ROTATIONAL BARRIERS IN FORMAMIDE AND ACETAMIDE: THE EFFECTS OF POLARIZATION FUNCTIONS AND CORRELATION

PAUL G. JASIEN*,

WALTER J. STEVENS and MORRIS KRAUSS

Molecular Spectroscopy 20899 (U.S.A.)

Division, National Bureau of Standards, Gaithersburg, MD

(Received 10 February 1986)

ABSTRACT Ab initio calculations have been used to determine the gas-phase rotational barrier about the CN bond in formamide and acetamide. The results indicate that the inclusion of polarization functions in the basis set leads to a substantial decrease (ca. 5 kcal mol-‘) in the calculated barrier height at the SCF level. Electron correlation effects decrease the barrier by less than 1 kcal mol-‘, while the addition of zero point energy corrections changes the barrier height only slightly. Based upon the current calculations, the 0 K rotational barriers for isolated formamide and acetamide are predicted to be 14.2 and 12.5 kcal mol.‘, respectively. INTRODUCTION

The energy barrier to rotation about an amide bond is of fundamental importance in the modeling of protein conformation. As such, there has been a great deal of study, both experimental [l-4] and theoretical [ 5- 131, to determine the magnitude of this barrier in simple model compounds. The compounds most often viewed as appropriate models for the amide bond in proteins are formamide and the various methyl substituted formamides. The experimental data on the formamide system have revealed a significant dependence of the rotational barrier on the solvent in which the experiment is done. Experimental determinations of this barrier seem to indicate that the height increases as the strength of the solvent hydrogen bonding increases. These results give values of 21.3, 16.8, and 16.9 kcal mol-’ [l] for the rotational barrier when the solvents are water, dioxane and acetone, respectively. Further work on neat formamide has yielded a value for this barrier of 18.9 kcal mol-’ [l] ; however, this result may not be directly comparable to the others due to the tendency of formamide to form linear chains, and therefore should not be viewed as a normal solute- solvent system. The experimental barrier for acetamide has been determined to be of the order of 1 kcal mol-’ lower than that of formamide [4] . *NBS-NRC

Postdoctoral

fellow.

198

It is therefore of interest to determine the inherent in vacua rotational barrier for molecules of this type. Towards this end, there have been numerous theoretical calculations on this problem. Unfortunately, space does not permit all of these studies to be mentioned, so only selected studies will be referred to. An early study by Radom et al. [ 51, using a 4-31G basis at the SCF level, yielded a value for the rotational barrier in formamide of 24.7 kcal mol-‘. This was shown by Nalewajski [8] to be an overestimate due to certain geometric constraints which were imposed. His calculations showed that geometry optimization led to a decrease of about 6 kcal mol-’ in the barrier height, demonstrating the importance of geometry relaxation in this problem. Radom and Riggs [ 91, performing full gradient-optimization calculations at the minimal basis STO-3G SCF level, found an extremely low barrier of 8.2 kcaI mol-‘. Calculations with the larger 4-31G basis set, at the STO-3G optimized geometry, yielded a value of 20.0 kcal mol-‘. In addition, they found a second saddle point corresponding to a transition state at higher energy (15.2 kcal mol-l) in the minimal basis set calculations, which disappeared when the larger basis set was used. Williams et al. [12], using a 4-21G basis set at the SCF level with full geometry optimization, obtained a barrier height of 18.1 kcal mol-‘. A calculation of the formamide barrier which is worthy of a short discussion is a study made in 1970 by Christensen et al. [6]. In this study, a large basis set of approximately triple zeta (TZ) quality was used and a barrier of 20 kcal mol-’ was obtained. Although full optimizations were not performed, a large number of structures were looked at. Additional calculations were done in which polarization functions were added to the heavy centers. (The previous calculations had already included p polarization functions on the hydrogens.) From their results they concluded that the calculated barrier is insensitive to the addition of polarization functions. Harding and Goddard [7], in a paper dealing with the low-lying excited states of formamide, reported a value of 14.2 kcal mol-’ for the barrier to rotation, based upon generalized valence bond (GVB) calculations with a double zeta (DZ) basis set. Using a three-pair GVB wavefunction as the reference, they went on to perform a very limited configuration interaction (CI) calculation within the GVB space and obtained a barrier height of 18.2 kcaI mol-‘. To the best of our knowledge, these calculations are the only ones in the literature which have investigated the effect of correlation on the barrier height. Unfortunately, the basis set which was used was probably not sufficient to give an accurate indication of the correlation effect. Furthermore, no geometry optimizations were performed since their major goal was to look at different electronic states and not the barriers to rotation. Theoretical calculations on the various methyl substituted formamides have been reported in the literature [lO-121. In this paper, discussion is confined to the case of acetamide. For this molecule, Radom and Riggs [ll] found barriers of 14.7 and 20.3 kcal mol-’ at the SCF level with STO-3G and 4-31G basis sets, respectively. Williams et al. [l2] reported a value of 18 kcal mol-’ based upon their 4-21G SCF calculations.

199

It can be seen that all previous ab initio calculations, done at approximately the DZ SCF level, have obtained results which are within the range of the experimentally determined values. However, care must be taken when making these direct comparisons, since the theoretical and experimental numbers are not directly comparable. First, the theoretical results have not been corrected for any difference in zero-point energies between the transition state and the equilibrium structure; it may be argued that such effects are small when compared to the overall energy of activation, and indeed this may be the case. A large inconsistency exists, however, when comparing the experimental and theoretical results. The ab initio calculations yield in vacua values for the rotational barriers, and because these barriers are seen experimentally to decrease as the ability of the solvent to hydrogen bond decreases, it would be expected that the solvent free value would be a few kcal mol-’ lower than the lowest experimental value. In the case of the GVB results of Harding and Goddard [7] this is true, but even these results gave a barrier in line with the small basis SCF calculations when approximations to the correlation energy were added. Only in the case of the semi-empirical PCILO calculations by Hofmann and Peinel [ 131 has the calculated in vacua result been predicted to be lower than those for the solvated species. In this paper, it will be demonstrated that the apparent inconsistency in the previous ab initio calculated rotational barriers may be accounted for by the lack of polarization functions in the basis sets. In addition, correlated results for these rotational barriers from extensively correlated wavefunctions will be reported, as will an estimate of the zero-point energy difference for the planar and twisted structures. With this data available, we will provide an estimate for the inherent barrier to rotation for formamide and acetamide at 0 K. CALCULATIONS

All SCF and GVB calculations in this study were performed with the HONDO [14] series of programs. The DZ basis sets used for heavy atoms were 3,l contractions of those given by Stevens et al. [15]. The compact effective potentials used to replace the 1s core in C, N and 0 were also taken from ref. 15. Basis sets of TZ quality were 2, 1, 1 recontractions of the DZ sets. The DZ set which was used for hydrogen was taken from Dunning [16] , while the TZ set was from Huzinaga [ 171. Calculations which included polarization functions used exponent values given as follows: Q(C) = 0.75, Q(N) = 0.77, ad(O) = 0.80, and CQ,(H)= 1.00. Calculations which included polarization functions on the heavy atoms, but not on the hydrogens are labelled as “+d”, while those utilizing a full set of polarization functions on all atoms are labelled with a “+P”. Test calculations were performed supplementing the basis set with only a single polarization function on the nitrogen and are labelled as DZ + d(N). Other calculations were done

200

by using a set of polarization functions for all atoms and an additional diffuse d function (Q = 0.10) on the nitrogen. This calculation is labelled as TZ + P + d’(N), the prime indicating that the function was diffuse. Full gradient optimizations were performed for both the equilibrium and transition-state structures under a C, symmetry constraint. All optimizations were done at the SCF level and the gradients were converged to better than 0.0005 a.u. for all coordinates. The C, symmetry constraint which was imposed on the planar structure of formamide and acetamide may result in a final structure which is not an absolute minimum in all coordinates. However, as has been discussed previously in the literature [ 181, the non-planarity of formamide due to pyramidalization of the nitrogen leads to no significant energetic change, because of the extreme flatness of the potential for this motion. In the case of the transition states, imposing C, symmetry on the molecule does not place any constraint on the pyramidalization of the NH2 group. GVB calculations, analogous to those of Harding and Goddard [7] were also performed in this study. The GVB pairs which were considered corresponded to a CO u bond, a CO 7~ bond, and a CN u bond. Calculations with the DZ basis set included all three GVB pairs, while those done with the DZ + P basis set included only the important CO 71bond GVB pair. The effect of correlation on the calculated energies for the rotational barrier in formamide was also studied. These effects were investigated via single-point energy evaluations with both the DZ + P and TZ + P + d’(N) basis sets. Calculations were done at the DZ + P SCF optimized geometries and were based upon two types of correlated wavefunctions. The first type included all single and double excitations from the single determinant SCF reference (CIDS), while the second included, in addition to single and double excitations, all higher order even excitations (quadruples, hextuples, etc.). These higher order effects were calculated via the approximate coupled cluster doubles (ACCD) method developed by Bachrach et al. [ 191. For both types of correlated wavefunctions, the calculations were carried out using the SCEP80 program [20]. The zero-point energy (ZPE) corrections for the rotation of formamide were estimated from calculated vibrational frequencies. These frequencies were obtained via finite differencing of the analytic first derivatives of the SCF wavefunction. The basis sets were of DZ and DZ + d quality and the calculations were carried out at the optimized structure for the particular basis set All zero-point energies were obtained by summing the 3N-6 (3N-7 for the transition state) vibrational frequencies. RESULTS

AND DISCUSSION

The calculated rotational barriers for formamide and acetamide at a number of different calculational levels are presented in Tables 1 and 2. All results are reported for fully optimized structures, except where indicated.

201 TABLE 1 Calculated rotational barriers for formamide (kcal mol-‘) Calculation level

Energy

DZ/SCF DZ + d( N)/SCF DZ + d/SCF DZ + P/SCF DZ + P/CIDS= TZ + P/SCFa TZ + P + d’(N-)/SCF’ TZ + P + d’(N)/CIDS’ TZ + P + d’( N-)/ACCDa

19.1 16.2 14.5 14.9 14.5 15.4 15.2 15.0 14.3

‘These calculations used the DZ + P SCF optimized geometries.

From Table 1, it can be seen that the DZ results obtained in this work are in good agreement with those done at a comparable level in other studies. With the addition of polarization functions on the heavy atoms (DZ + d), there is a dramatic decrease in the calculated barrier height for the formamide and acetamide molecules (4.6 and 5.2 kcal mol-‘, respectively). Further addition of polarization functions to the hydrogen atoms in formamide modifies the barrier only slightly (0.4 kcal mol-l) from that found at the DZ + d SCF level. Calculations done with the DZ + P basis at the optimized structure from the DZ calculations give a value of 17.2 kcal mol-’ for the barrier, halfway between the DZ and DZ + P results. This result indicates the importance of further geometry optimization after the addition of polarization functions. Extension of the basis set to TZ + P or TZ + P + d’(N), with no further geometry optimization beyond DZ + P, has only relatively minor effects on the energy barrier, yielding values of 15.4 and 15.2 kcal mol-’ compared to the DZP result of 14.9 kcal mol-‘. The effect of correlation on the rotational barrier in the case of formamide is also presented in Table 1. For the CIDS calculation utilizing the DZ + P basis, the effect of correlation is small, lowering the barrier by only 0.4 kcal mol-l, compared with the SCF result. In the larger basis CIDS calcuTABLE 2 Calculated rotational barriers for acetamide (kcal moldI) Calculation level

Energy

DZ/SCF DZ + d(N)/SCF DZ + d/SCF

17.9 11.6 12.7

202

lation the lowering is comparable, being only 0.2 kcal mall’. The near coincidence of the changes in energy due to correlation for two calculations with different sizes of basis set indicates that the small correlation effect is unlikely to be an artifact due to basis-set restrictions. Estimates of the differential higher order correlation, via the Davidson correction indicated that these effects may be as large as 1 kcal mol-l. Since an energetic effect of this size is larger than that found from the inclusion of single and double substitutions into the wavefunction, further investigation was considered to be worthwhile. In order to test these effects, the explicit determination of the higher order substitutions via an ACCD wavefunction was performed. Inclusion of these terms in the CI wavefunction leads to a further lowering of the barrier by 0.7 kcal mol-‘, which is indeed larger than that from single and double substitutions alone. Although no CI calculations were performed for the acetamide system, the results would be expected to be analogous to those in formamide. This expectation is based upon the similarity in the basic change in electronic structure in the two systems upon rotation (i.e., breaking the NC0 a bond conjugation). The results in Table 3 present the optimized O-C-N-H torsional angle at the transition state for both formamide and acetamide. The value of this angle is indicative of the extent of pyramidalization about the nitrogen atom and it may be observed that these values are quite different from those found for the equilibrium structures for which this angle is 0”. A complete tabulation of the optimized geometric parameters for these systems for the DZ, DZ + d(N), and DZ + d basis sets is given in Tables 4 and 5. The general trends which we find, for the changes in bond lengths and angles on going from the unrotated to the rotated structure, are in line with those found in other studies. No further discussion of these trends will be given here; however, the interested reader may wish to refer to the paper by Williams et al. [12]. From the data in Table 3, it is clear that there is a large change in the apparent hybridization of the nitrogen on the addition of polarization TABLE 3 O-C-N-H

Torsional angle for transition state structures (degrees) Formamide

DZ/SCF DZ+ d(N) DZ+d DZ+P 4-21G/SCFa aRef. 12.

67.7 57.7 56.7 56.9 64.0

Acetamide 65.1 56.0 56.0 62.3

203 TABLE 4 Calculated structural parameters for equilibrium and transition states of formamide@ DZ + d(N)

DZ

WC01 R(CN) R(CH) R(NH,)C R(NHt)C
DZ+d

Equi.

TS

Equi.

TS

Equi.

TS

1.234 1.370 1.087 0.998 0.996 124.3 121.6 119.3 121.4 0.0

1.223 1.435 1.085 1.003 1.003 124.3 120.5 115.5 115.5 67.7

1.235 1.365 1.088 1.002 0.999 124.3 121.5 118.9 120.9 0.0

1.218 1.437 1.092 1.009 1.009 122.6 119.7 108.6 108.6 57.7

1.206 1.364 1.096 1.001 0.998 124.7 122.4 119.0 121.5 0.0

1.197 1.440 1.093 1.010 1.010 124.9 121.4 107.8 101.8 56.7

aBond lengths in .k ; angles in degrees. bCalculated parameters using the DZ + P basis set give results in agreement with the DZ + d basis set to 0.002 I%in bond lengths and 0.2” in bond angles. ‘H, labels the hydrogen cis to the oxygen, while Ht labels the hydrogen trans to it.

functions, as can angle when going seems to be well barrier for the two

be seen by the approximately 9” change in the reported from a DZ to a DZ + P basis set. Furthermore, this angle correlated with the calculated height of the rotational molecules.

TABLE 5 Calculated structural parameters for equilibrium and transition states of acetamidea DZ

WW WCN) NCC)

RW-L~ RW-4 1

R( CH’) R( CH”)
DZ + d(N)

DZ+d

Equi.

TS

Equi.

TS

Equi.

TS

1.240 1.377 1.524 0.998 0.996 1.082 1.088 121.4 122.8 118.6 122.2 0.0 120.1 0.0

1.228 1.447 1.518 1.004 1.004 1.083 1.087 122.0 122.9 114.1 114.1 0.0 121.6 65.1

1.241 1.371 1.526 1.002 0.999 1.082 1.087 121.4 122.6 118.3 121.8 0.0 120.1 0.0

1.229 1.449 1.518 1.010 1.010 1.083 1.087 122.2 123.2 107.2 107.2 0.0 121.6 56.0

1.211 1.370 1.529 1.000 0.998 1.084 1.089 122.0 123.0 118.6 122.3 0.0 120.3 0.0

1.201 1.449 1.519 1.010 1.010 1.085 1.089 122.7 123.6 107.1 107.1 0.0 121.7 56.0

aBond lengths in A; angles in degrees. b& labels the hydrogen cis to the oxygen, while Ht labels the hydrogen trans to it.

204

This need for larger basis sets including polarization functions, in order to describe accurately the hybridization of the nitrogen atom as well as inversion barriers in amines, has been noted previously in the literature. Carlsen et al. [ 181, in their study of the inversion barriers in ammonia and formamide, found that the calculated values may be extremely basis-set dependent. More recently, Boggs and Niu [21] studied the effect of basis set on amine out-of-plane angles. They concluded that “the amino group angle is the result of a delicate balance between lone pair character of the nitrogen electrons and delocalization of these electrons into other portions of the molecule”, and further concluded that polarization functions are necessary in order to describe the lone-pair character. One further basis-set test was performed in this study in which calculations were performed with a DZ + d(N) basis. These results seem to indicate that with only d polarization functions on the nitrogen, the hybridization of this atom is described adequately compared with the fully polarized basis. However, although the barrier height is dramatically improved over the non-polarized basis, it still differs by l-2 kcal mol-’ from the fully d polarized basis-set results. This is probably due to an inadequate description of the delocalized nature of the electron density onto the carbonyl group, which necessitates the inclusion of polarization functions on the carbon and oxygen. At this point a short discussion of the relationship of the GVB calculations of Harding and Goddard [7] with the present correlated results is in order. We have performed calculations analogous to their three-pair GVB calculations utilizing our DZ basis set and find reasonable agreement with their results (12.8 versus 14.2 kcal mol-‘). Analysis of the resulting wavefunction indicates that the CO x bond GVB pair has the largest correlation effect of the three; its CI coefficient being 0.22 in the twisted form and 0.19 in the planar structure, as compared with the coefficients for the other two GVB pairs in both structures, which are less than 0.09. With this result in hand, we performed single-point calculations with the DZ + P basis set at the DZ + P SCF optimized geometry, with only the CO 7~bond in the GVB space. This calculation yielded a value of 11.0 kcal mol-’ for the rotational barrier and gave a CI coefficient for the GVB pair of 0.19 and 0.16 for the twisted and planar structures, respectively. These GVB results for the torsional barrier are certainly too low when compared to our correlated results. That this type of calculation overestimates the differential correlation effect is evidenced by the results of Harding and Goddard [7] on the inclusion of a slightly larger degree of correlation into their wavefunction. Their results indicated that the additional correlation led to an increase in the barrier of 4 kcal mol-I. The fairly significant coefficient for the CO 71 bond GVB pair, which has been calculated in the present study, may explain the reason for the relatively large effect on the energy barrier of the inclusion of higher order substitution effects in the ACCD calculations. The apparent higher order energy effects may be an artifact of the lack of a multi-configuration reference

205

for the CIDS wavefunction. In the present ACCD calculations, however, the constraint of a single configuration reference should not cause an appreciable error in the final calculated result, since the reference configuration is still the dominant term in the wavefunction and since the coupled cluster treatment of the correlation may correct for some of the deficiencies of the single reference configuration. In the present calculations, such an error is probably on the order of a few tenths of a kcal mol-’ at the most. Lastly, the calculated vibrational frequencies for the planar and rotated formamide molecules indicate that the ZPE correction is extremely small, of the order of 0.1-0.2 kcal mol-‘. Such a small correction is found despite the fact that the planar structure has an additional frequency, corresponding to the NH2 twist, added to its ZPE. This is due to the fact that in the transition state the vibrational mode corresponding to the NH2 inversion mode has increased dramatically in frequency. In addition, many of the intermediate energy modes (1000-1600 cm-‘) have slight shifts to higher energy. Frequencies for planar formamide have been reported previously, and the current values are in fair agreement with those of Sugawara et al. [22]. CONCLUSION

The results of this study demonstrate that previously calculated values for the rotational barriers in formamide and acetamide were overestimates of the true gas-phase values. The addition of polarization functions to the heavy centers, followed by full geometry optimization, results in a large decrease of about 5 kcal mol-’ in the calculated barriers. Correlation effects due to single and double substitutions from the reference are shown to have modest effects (<0.5 kcal mol-‘) in the case of the formamide barrier, although higher order effects lower the barrier even further. The large correlation effect in the CO 7~ bond may necessitate a multi-configurational reference in the case of a CIDS calculation. Given the current available data, our estimated value for the barrier to internal rotation in formamide is 14.2 kcal moT1, while that of acetamide is 12.5 kcal mol-‘. The differences in energies due to zeropoint effects have been included in these values; however, such effects are estimated to be quite insignificant. Because no other thermodynamic corrections have been included, these reported barriers correspond to a temperature of 0 K. REFERENCES 1 2 3 4 5 6

H. Kamei, Bull. Chem. Sot. Jpn., 41 (1968) 2269. T. Drakenberg and S. Forsen, J. Phys. Chem., 74 (1970) 1. T. Drakenberg, K. I. Dahlgvist and S. Forsen, J. Phys. Chem., 76 (1972) 2178. T. Drakenberg, Tetrahedron Lett., 18 (1972) 1743. L. Radom, W. A. Lathan, W. J. Hehre and J. A. Pople, Aust. J. Chem., 25 (1972) 1601 D. H. Christensen, R. N. Kortzeborn, B. Bak and J. J. Led, J. Chem. Phys., 53 (1970) 3912. 7 L. B. Harding and W. A. Goddard III, J. Am. Chem. Sot., 97 (1975) 6300.

206 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

R. F. Nalewajski, J. Am. Chem. Sot., 100 (1978) 41. L. Radom and N. V. Riggs, Aust. J. Chem., 33 (1980) 249. L. Radom and N. V. Riggs, Aust. J. Chem., 33 (1980) 2337. L. Radom and N. V. Riggs, Aust. J. Chem., 35 (1982) 1071. J. 0. Williams, C. Van Alsenoy and L. Schafer, J. Mol. Struct., (Theochem), 76 (1981) 171. H. J. Hofmann and G. Peinel, J. Mol. Struct., (Theochem), 91 (1983) 289. M. Dupuis and H. F. King, Int. J. Quantum Chem., 11 (1977) 613; J. Chem. Phys., 68 (1978) 3998. W. J. Stevens, H. Basch and M. Krauss, J. Chem. Phys., 81 (1984) 6026. T. H. Dunning, J. Chem. Phys., 53 (1970) 2823, S. Huzinaga, J. Chem. Phys., 42 (1965) 1293. N. R. Carlsen, L. Radom, N. V. Riggs and W. R. Rodwell, J. Am. Chem. Sot., 101 (1979) 2233. S. M. Bachrach, R. A. Chiles and C. E. Dykstra, J. Chem. Phys., 75 (1981) 2270. C. E. Dykstra, R. A. Chiles and M. D. Garrett, J. Comput. Chem., 2 (1981) 266. J. E. Boggs and Z. Niu, J. Comput. Chem., 6 (1985) 46. Y. Sugawara, Y. Hamada, A. Y. Hirakawa, M. Tsuboi, S. Kato and K. Morokuma, Chem. Phys., 50 (1980) 105.