Theoretical Investigations of the Rotational Barrier in Anisole: An ab

sumed a 2-fold barrier, which implies a relatively sharp maximum for the perpendicular ... minimum energy conformation at 0 = 18'. Seip and Seip pro- ...
0 downloads 0 Views 1MB Size
Download PDF
J . Phys. Chem. 1990, 94, 4483-4491

4483

Theoretical Investigations of the Rotational Barrier in Anisole: An ab Initio and Molecular Dynamics Study David C. Spellmeyer: Peter D. J. Grootenhuis; Michael D. Miller,* Lee F. Kuyper,g and Peter A. Kollman**f Department of Pharmaceutical Chemistry, University of California, San Francisco, San Francisco, California 94143, Wellcome Research Labs, Burroughs Wellcome Company, Research Triangle Park, North Carolina 27709, and IBM Research Center, Almaden Research, San Jose, California 95120 (Received: December 1, 1989)

Theoretical investigations at the AMI, RHF/3-21G, RHF/6-31G*, MP2/6-3lG*//RHF/3-21G, and MP2/6-31G*// RHF/6-3IG* levels of theory have been used to study the gas-phase rotational barrier in anisole. Geometries at all levels are similar. Vibrational frequency analyses indicate that the coplanar structure is a local minimum at both the AMI and RHF/3-21G levels. The conformer with the methoxy substituent perpendicular to the phenyl ring is a transition structure at the AMI level but is a minimum at the RHF/3-21G level. A transition structure between these minima is located at both the RHF/3-21G and the RHF/6-31G* levels. The potential at either ab initio level is very flat between CP = 45O and 90' but is rather more steep with AMI. Calculations at the MP2/6-3lG*//RHF/3-21G and MP2/6-3lG*//RHF/6-31G* levels indicate that the perpendicular conformer is slightly higher in energy than the transition structure. The 4-fold term for rotation is small, but very significant, based on the ab initio energy profile. Molecular mechanics parameters derived for the AMBER force field have been used in an isothermal-isobaric simulation of liquid anisole to assess the effect of the liquid structure on the dihedral angle population. The population of the intermediate values of CP is slightly higher than in the liquid than is predicted from the gas-phase energy profile, suggesting a slightly decreased barrier to rotation in the liquid, in direct contrast to experiment.

Introduction The determination of the rotational barrier of the methoxy group in anisole, 1, has been the subject of a great deal of experimental and theoretical studies recently.'-'* These studies H14

I

F'

C8-H15 \H16

Hll 1

disagree not only on the magnitude of the gas-phase barrier height but also on the shape of the energy profile itself. In addition, the barrier is found to be higher in the solid than in the liquid and higher in the liquid than in the gas phase. Most of the recent studies have focused only on the gas-phase barrier and on the barrier in dilute solutions. Very little work has been done to rationalize the phase dependence of the rotational barrier. We present the first such theoretical study, focusing on the barrier in the gas phase. Most experimental studies of this rotational barrier have assumed a 2-fold barrier, which implies a relatively sharp maximum for the perpendicular conformation. Fewster measured the gas-phase barrier as 3.6 kcal/mol.' Friege has measured the barrier in the gas phase as 5.7 kcal/mol.* Owen and Hester have measured the barrier in liquid anisole as 6.0 kcal/moL3 Kerr constant measurements," however, suggest that there is a nonplanar minimum energy conformation at 0 = 18'. Seip and Seip proposed the existence of two stable conformers to rationalize this conflict dataas Konschin et al. find that two stable conformers exist at 130 K, separated by about 0.6 kcal/mol.6 The barrier between these two conformers has been estimated as about 12 kcal/mol in the solid.' Recent N M R studies arrive at the "inescapable" conclusion that the barrier is both 2-fold and 4-fold 'University of California, San Francisco. IBM,Almaden. 1 Burroughs Wellcome.

*

0022-3654/90/2094-4483$02.50/0

in nature, indicating the presence of a second stable conformer in the gas phase.E Assuming that this is the case, the barrier height is estimated as less than 3.1 kcal/mol, though the exact value has not been determined thus far. All of the theoretical studies to date tend to agree on the relatively low barrier to rotation of the methoxy group of about 1-2 kcal/mol. Unlike the experimental studies, most all agree on the existence of a second conformer. The estimates of the difference in energy between the conformers and on the barrier height do not concur, however. Nagy and JancBrova's PCILO calculations both show a barrier of 1.5 kcal/mol using fixed ge~metries,~ while Gerhards et al. find a barrier of 2.1 kcal/mol with the same method.1° Anderson, et al., employing STO-3G calculations, find a barrier of about 1 kcal/mol using standard geometries and a slightly higher barrier with partial optimization." Konschin et al. find the planar conformer to be about 2 kcal/mol more stable than the perpendicular conformer utilizing partial geometry optimization at the STO-3G leveL6 These authors state that the "V4 component for gas-phase barriers is small but not unimportant". Klessinger and Zywietz have also used STO-3G calculations to show that the presence of the second conformer is dependent upon optimization of the CO bond lengths and CCO and COC bond angles.I2 At the 4-31G//STO-3G level they find a barrier of 1.9 kcal/mol. The perpendicular conformer is 1.70 kcal/mol less stable than the planar at this level. They also suggest that at least some 4-fold component must be present in the rotational barrier. Finally, Schaefer and Penner have calculated ( I ) Allen, G.; Fewster, S. In Internal Rotations in Molecules; OrvilleThomas, W. J., Ed.; Wiley: New York, 1974; Chapter 8. (2) Friege, H.; Klessinger, M. Chem. Ber. 1979, 112, 1614. (3) Owen, N. L.; Hester, R. E. Spectrochim Acia, Pari A 1969, 25, 343. Lister, D. G.;Owen, N. L. J. Chem. SOC.,Faraday Trans. 2 1973,69, 1304,. (4) Aroney, J.; Lefoer, R. J. W.; Pierens, P. K.; The, M. G.N. J. Chem. SOC.B 1969, 666. ( 5 ) Seip, H. M.;Seip, R. Acra Chem. Scand. 1973, 27, 4024. (6) Konschin, H.. Tylli, H.; Grundfelt-Forsius, C. J. Mol. Struct. 1981, 77 . . , 51

(7) Tylli, H.; Konschin, H. J. Mol. Struct. 1977, 42, 7. (8) Schaefer, T.; Penner, G. H. Can. J . Chem. 1988, 66, 1635, 1641. Schaefer, T.; Penner G. H. J. Mol. Struct. (THEOCHEM)1987, 157, 179. Schaefer, T.; Sebastian, R. Can. J. Chem. 1989, 67, 1148. (9) Nagy, L. T.; JancBrovl, M.Chem. Papers 1985, 39, 289. (10) Gerhards, J.; Ha, TK.; Perlia, X . Helu. Chim. Acta 1982, 65, 105. (1 1) Anderson, G. M., 111; Kollman, P. A,; Domelsmith, L. N.; Houk, K. N. J. Am. Chem. SOC.1979, 101, 2344. (12) Klessinger, M.; Zywietz, A. J. Mol. Struct. 1982, 90, 341.

0 1990 American Chemical Societv

4484

The Journal of Physical Chemistry, Vol. 94, No. 11, 1990

Spellmeyer et al.

TABLE I: Bond Lengths (in A) of Anisole Calculated at Specific Values of 0 parameter CI-C2 C2-C3 C3-C4 C4-C5 C5-C6 C6-CI C1-07 07-C8 C2-H9 C3-HI0 C4-Hll C5-HI2 C6-H 13 C8-H I4 C8-H 15 C8-HI6

9 = Oo

IOo

1.3816 1.3877 1.3788 1.3892 1.3762 1.3892 1.3718 1.4356 1.0695 1.0722 1.0712 1.0720 I ,0696 I .083 1 1.0776 1.0831

1.3818 1.3876 1.3791 1.3889 1.3765 1.3890 1.3721 1.4360 1.0696 1.0722 1.0712 1.0720 I ,0696 I ,0832 1.0777 1.0829

20' 1.3820 1.3870 1.3797 1.3884 1.3772 1.3882 1.3736 1.4376 1.0698 1.0722 1.0713 1.0720 I ,0696 I ,0832 1.0778 1.0826

30' 1.3823 1.3862 1.3806 1.3876 1.3783 1.3870 1.3759 1.4399 1.0701 1.0721 1.0713 1.0720 1.0697 1.083 1 1.0780 1.0823

RHF/3-21G 40' 50' 60' 1.3826 1.3829 1.3830 1.3852 1.3842 1.3834 1.3817 1.3829 1.3839 1.3866 1.3855 1.3847 1.3796 1.3810 1.3821 1.3855 1.3839 1.3826 1.3789 1.3825 1.3855 1.4424 1.4448 1.4470 1.0704 1.0705 1.0707 1.0721 1.0721 1.0720 1.0713 1.0714 1.0714 1.0720 1.0719 1.0719 1.0697 1.0698 1.0699 1.0829 1.0828 1.0828 1.0780 1.0780 1.0779 1.0819 1.0817 1.0817

the 2-fold and 4-fold components of the energy barrier at the STO-3G level.* They find a barrier of about 1.5 kcal/mol with local minima at = 0' and 90'. To our knowledge, no theoretical estimates exist of the barrier height for the pure liquid, obtained from molecular dynamics or Monte Carlo calculations, or for the solid, obtained from crystal calculations. Because alkyl aryl ethers are of interest in both biological and in host-guest system^,'^-'^ we sought to derive molecular mechanical parameters for the rotational barrier in anisole. Despite the abundance of experimental and theoretical data on anisole, we have found it impossible to determine accurate molecular mechanical parameters from what has been presented to date. We have, therefore, undertaken an extensive theoretical study of the rotational barrier in anisole.

Computational Details: Ab Initio Preliminary ab initio geometry optimizations of the coplanar and the perpendicular conformations of anisole were carried out at the R H F level of theory with the 3-21G basis setT6and the GAUSSIAN~O-UCSF program.17 Geometry optimizations employing this program were performed in Cartesian space rather than with internal coordinate gradient optimization, requiring some 73 CPU hours on a VAX-8600 running the VMS operating system. Single point energy calculations at the RHF/6-31F*I8 level were performed on the RHF-3-21G geometries, also with the GAUSSIANSO-UCSF program. Further optimizations were required, but were considered to be prohibitively expensive on this machine. Therefore, optimizations were continued with the GAUSSIANS~ packageIg on the San Diego Supercomputer Center's Cray-X/ MP-48 running the CTSS operating system. Full gradient optimizations on internal coordinates20 were carried out for the coplanar and the perpendicular conformations at the RHF/3-21G level of theory. Geometries were found to be within machine precision of those obtained with Cartesian optimization. Analytical vibration frequency calculations were performed on these structures. Both were found to have no imaginary vibrational frequencies, indicating that minima exist for the coplanar and perpendicular conformers on this potential surface. (13) Makriyannis, A.; Fesik, S. J. Am. Chem. SOC.1982, 104, 6462. (14) Bacconari, D. P.; Doluge, S.; King, R. W. Biochemistry 1982, 21, 5068. ( I S ) Grwtenhuis, P. D. J.; Kollman, P. A. J. Am. Chem. SOC.,in press. Grootenhuis, P. D. J.; Gerden, J. v.; Dijkstra, P. J.; Harkema, R. S.; Reinhoudt, D. N. J. Am. Chem. SOC.1987, 109, 8044. Kollman, P. A,; Wipff, G.; Singh, U. C. J. Am. Chem. SOC.1985, 107, 2122. (16) Binkley. J. S.; Pople, J. A.; Hehre, W. J. J. Am. Chem. SOC.1980, 102, 939. (17) Singh, U. C.; Kollman, P. A. J. Compuf. Chem. 1984, 5, 129. (18) Hariharan. P. C.; Pople, J. A. Theor. Chim. Acto 1973, 28, 213. ( I 9) GAUSSIAN86; Frisch, M. J.; Binkley, J. S.; Schlegel, H. B.; Raghavachari, K.; Melius, C. F.; Martin, R. L.; Stewart, J. J. P.; Bobrowicz, F. W.; Rohlfing, C. M.; Kahn, L. R.; DeFrees, D. J.; Seeger, R.; Whiteside, R. A,; Fox, D. J.: Fleuder, E. M.; Pople, J. A. Carnegie-Mellon Publishing Unit, Pittsburgh, PA, 1984. (20) Schlegel. H. B. J. Comput. Chem. 1982, 3, 214.

70' 1.3828 1.3829 1.3845 1.3843 1.3829 1.3817 1.3873 1.4488 1.0709 1.0720 1.0715 1.0719 1.0702 1.0828 1.0779 1.0821

80' 1.3825 1.3829 1.3847 1.3843 1.3832 1.3815 1.3887 1.4502 1.0709 1.0719 1.0715 1.0719 1.0704 1.0828 1.0779 1.0825

90' TS 1.3818 1.3830 1.3831 1.3833 1.3845 1.3839 1.3845 1.3846 1.3831 1.3822 1.3818 1.3825 1.389 1.3857 1.4505 1.4471 1.0707 1.0707 1.0719 1.0720 1.0715 1.0715 1.0719 1.0719 1.0707 1.0700 1.0827 1.0828 1.0779 1.0779 1.0827 1.0817

RHF/6-3 1G* CP = Oo 90' TS 1.3858 1.3910 1.3794 1.391 1 1.3778 1.3936 1.3497 1.3982 1.0727 1.0757 1.0748 1.0756 1.0742 1.0851 1.0798 1.0851

1.3856 1.3854 1.3860 1.3860 1.3854 1.3856 1.3633 1.4050 1.0750 1.0755 1.0751 1.0755 1.0750 1.0854 1.0805 1.0854

1.3872 1.3855 1.3854 1.3859 1.3847 1.3860 1.3606 1.4041 1.0745 1.0756 1.0750 1.0755 1.0744 1.0852 1.0804 1.0842

. 2-

P

1.75' 1.50'

1.25'

.I

-f

1.00-

m

0.75' 0.50-

o 0 0.00

10

20 . 30

50

40

260

70

80

90 5

Q,

Figure 1. The relative energies of the rotational barrier in anisole as a function of the dihedral angle 9.

U

Figure 2. Geometry of the 3-21G-optimized transition structure for rotation of the methoxy group in anisole. The vectors represent the imaginary vibrational frequency calculated at the 3-21G level. Many of the vectors are not shown as they are very small relative to the movement of the methyl group and the hydrogens on the phenyl group.

Location of the transition structure for the interconversion of these two minima was performed with partial optimizations, varying the dihedral angle (C2-C,-07-C8;see 1 for numbering scheme) from 80' to 10' in 10' increments. All other variables were optimized at the RHF/3-21G level. At each point, the approximate Hessian from the previous point was used, as was the wave function, resulting in some computational savings. In most cases, three of four optimization cycles were sufficient for geometry convergence. Optimized geometries at all levels of theory are shown in Tables 1-111. The resulting energy curve (Figure 1) shows a weak maximum near the 60° conformer. A full transition state optimization employing MNDOZ'force constants as the initial guess of the Hessian provided the transition structure. Analytical vibrational frequency analysis indicates that this is indeed a transition structure at this level of theory. Figure 2 shows the vectors of the imaginary vibrational frequency at the 3-21G level. Finally, single point energy calculations at the RHF/6-3 lG* level were performed on seiected 3-21G geometries. Because the relative energies obtained at this level were similar to those obtained at the 3-21G level, only those of the coplanar, the 90°, and (21) Dewar, M. J. S.; Thiel, W. J. Am. Chem. SOC.1977, 99, 4899.

~

The Journal of Physical Chemistry, Vol. 94, No. I I, I990 4485

Rotational Barrier in Anisole

TABLE 11: Bond Angles (in deg) of Anisole Calculated at Specific Values of parameter

c I -c2-c3 c2-c3-c4 c 3-c4-c 5 C4-C5-C6 C5-C6-C1 C2-C 1-C6 c 2 - c 1-07 C6-C 1-07 C1-07-C8 C I-C2-H9 C3-C2-H9 C2-C3-H IO C4-C3-H I O C3-C4-H I 1 C5-C4-H 1 I C4-C5-H 12 C6-C5-H 12 C 5 C 6 - H 13 C I-C6-H 13 07-C8-H 14 07-C8-H 15 07-C8-H 16 H I4-C8-H I5 H 14-C8-H 16 H 15-C8-H 16

0 = 0' 119.8 120.7 119.1 120.5 120.2 119.6 124.4 116.0 120.8 120.9 119.3 119.2 120.1 120.5 120.4 120.0 119.5 121.7 118.1 111.4 105.6 1 11.4 109.6 109.3 109.6

IO" 119.8 120.7 119.2 120.5 120.2 119.7 124.3 116.1 120.7 120.8 119.4 119.2 120.1 120.5 120.4 120.0 119.5 121.7 118.1 111.3 105.6 111.4 109.7 109.3 109.5

20'

30"

119.8 120.7 119.2 120.4 120.2 119.7 124.0 116.3 120.5 120.7 119.5 119.2 120.1 120.5 120.4 120.0 119.6 121.7 118.1 1 1 1.2 105.7 11 1.4 109.8 109.3 109.3

119.8 120.6 119.2 120.4 120.1 119.7 123.7 116.5 120.1 120.6 119.6 119.3 120.1 120.4 120.3 120.0 119.6 121.7 118.2 111.2 105.7 111.4 109.9 109.3 109.3

40'

Q

RHF/3-21G 50" 60'

119.9 120.5 119.3 120.4 120.1 119.8 123.3 116.9 1 19.6 120.4 119.8 119.4 120.0 120.4 120.3 120.0 119.6 121.6 118.3 1 1 1.1 105.8 111.4 109.9 109.3 109.3

119.9 120.5 119.4 120.3 120.1 119.8 122.8 117.4 118.9 120.2 119.9 119.5 120.0 120.3 120.3 120.0 119.6 121.5 118.4 111.0 105.7 111.5 109.9 109.3 109.3

119.9 120.4 119.5 120.3 120.0 119.9 122.1 118.0 118.0 119.9 120.2 119.6 120.0 120.2 120.3 120.0 119.7 121.4 118.5 110.8 105.8 111.5 110.0 109.2 109.5

70' 119.9 120.3 119.5 120.2 120.0 120.1 121.2 118.7 117.1 119.6 120.5 119.7 120.0 120.2 120.2 120.0 119.1 121.3 118.7 110.7 105.9 111.3 110.0 109.2 109.7

80" 119.9 120.3 119.6 120.2 119.9 120.1 120.5 119.3 116.4 119.3 120.8 119.7 120.0 120.2 120.2 120.0 119.7 121.2 118.9 110.8 106.0 111.0 110.0 109.2 109.9

90" 119.9 120.2 119.6 120.2 119.9 120.2 119.9 119.9 116.2 119.1 121.0 119.8 120.0 120.2 120.2 120.0 119.8 121.0 119.1 110.9 106.0 110.9 110.0 109.1 110.0

70' 0.0 -0.2 0.4 -0.5 179.2 -178.9 -178.8 -179.7 -179.9 -179.7 -179.7 70.0 -110.8 61.2 -179.6 -60.4

80" 0.1 -0.3 0.4 -0.4 179.4 -179.3 -179.0 -179.6 180.0 -179.6 -179.5 80.0 -100.5 61.6 -179.1 -59.8

90' 0.2 -0.3 0.3 -0.2 179.3 -179.3 -179.2 -179.7 180.0 -179.6 -179.2 90.3 -90.3 60.7 180.0 -60.7

TS 119.9 120.4 119.5 120.3 120.0 120.0 122.0 118.0 118.0 119.9 120.2 119.6 120.0 120.2 120.3 120.0 119.7 121.4 118.5 110.8 105.8 111.5 1 10.0 109.2 109.5

RHF/6-3IGZ 0 =Oo 90' TS 119.5 119.7 119.7 121.1 120.4 120.6 119.0 119.6 119.4 120.6 120.4 120.4 120. I 119.7 119.9 119.7 120.3 120.1 119.8 121.9 124.5 119.8 118.0 115.7 119.8 115.5 117.2 121.2 119.2 120.1 119.3 121.1 120.3 118.9 119.6 119.4 120.0 120.0 120.0 120.6 120.2 120.3 120.5 120.2 120.3 120.0 120.0 120.0 119.4 119.6 119.5 121.4 121.1 121.3 118.5 119.2 118.8 111.2 111.0 1 1 1.5 106.9 106.6 106.2 111.5 1 1 1.2 111.9 109.1 109.4 109.4 109.3 108.9 109.0 109.1 109.4 108.9

TABLE 111: Dihedral Angles (in dep;) of Anisole Calculated at Specific Values of Q parameter 0 = 0' 10' c I-C2-C3-C4 0.0 -0.1 0.0 -0.1 c 2-c 3-c4-c5 C3-C4-C5-C6 0.0 0.0 0.0 0.0 C4-C5-C6-C1 -1 79.9 c3-c2-c 1-07 179.4 C 5 C 6 - C 1-07 179.9 -179.4 -179.9 -179.7 C6-C 1-C2-H9 180.0 180.0 CI-C2-C3-H I O 180.0 C2-C3-C4-H 1 1 180.0 180.0 180.0 C3-C4-C5-H 12 180.0 C4-C5-C6-H 13 180.0 10.0 C2-Cl-07-C8 0.0 C6-CI-07-C8 180.0 -170.8 61.2 55.7 C 1-07-C8-H 14 179.9 174.7 C 1-07-C8-H 1 5 C 1-07-C8-H 16 -61.1 -66.6

20' -0.2 -0.1 0.2 -0.1 178.8 -178.7 -179.4 179.9 180.0 -179.9 180.0 20.0 -161.5 50.8 169.9 -71.4

30' -0.3 -0.1 0.3 -0.2 178.4 -178.3 -179.2 179.9 -179.9 -179.9 -179.9 30.0 -152.0 47.2 166.4 -75.0

RHF/3-21G 40' 50' 60' -0.4 -0.3 -0.1 -0.1 -0.1 -0.1 0.4 0.4 0.4 -0.3 -0.4 -0.5 178.2 178.3 178.8 -178.1 -178.1 -178.5 -179.1 -178.9 -178.8 180.0 -179.9 -179.8 -179.9 180.0 -179.9 -179.9 -179.8 -179.8 -179.9 -179.9 -179.8 40.0 50.0 60.0 -142.2 -131.9 -121.3 46.5 50.7 57.6 165.8 169.9 176.8 -75.6 -71.4 -64.2

the transition structures were deemed necessary. Geometries of the coplanar, perpendicular, and transition structures were optimized at the RHF-6-31G* level to assess the effect of basis set upon the geometry. The coplanar and the perpendicular conformations were performed on the IBM-3090, while the transition structure was located on the SDSC Cray-XMP/48. The computational expense was diminished with the use of the full second derivative matrix from the 3-2 1G harmonic frequency analysis as the initial Hessian for the transition-state optimization. Single point energy calculations were carried out at the MP2/6-31G* level of theoryz2on the fully optimized RHF/3-21G and on the RHF/6-31G* geometries in order to assess the effect of electron correlation on the potential surface. These calculations were performed with the G A U S S I A N Bpackage ~ on an IBM-3090300-E running the VM operating system, and required about 3 h of CPU time each. In order to make the comparison between theory and experiment more complete, semiempirical calculations of the rotational barrier were performed with the reaction pathway feature of M O P A C . ~ ~ These calculations were carried out with the Hamiltonian (22) Msller, C.; Plesset, M. S. Phys. Reu. 1934, 46, 618. (23) Stewart, J. J . Quantum Chemistry Program Exchange, Program N.

455, Bloomington, IN. (24) Dewar, M. J . S.; Zoebisch, G.G.;Healy, E. F. J . Am. Chem. SOC. 1985, 107, 3902.

TS -0.1 -0.2

0.5 -0.5 178.8 -178.6 -178.8 -179.8 -179.9 -179.8 -179.8 60.9 -120.4 58.2 177.3 -63.7

RHF/6-3 IG* 90' TS 0.0 0.0 -0.4 0.0 0.2 -0.1 0.0 -0.2 0.5 0.0 0.0 -0.4 180.0 -178.6 178.3 180.0 178.6 -178.0 180.0 179.0 -178.6 180.0 179.7 -179.9 180.0 179.9 -179.8 180.0 179.5 -179.8 180.0 179.2 180.0 -0.1 -90.8 59.7 180.0 90.8 -122.4 -61.2 -60.7 54.9 180.0 180.0 173.9 61.3 60.7 -67.2 = 0'

contained in the MOPAC package and run on a VAX 8650 running BSD 4.3 Unix. Vibrational frequency calculations verify that the coplanar conformer is a minimum and that the perpendicular conformer is a maximum on this potential energy surface. Geometries of the coplanar and the perpendicular conformations are presented in Table V. The energies are presented in Table IV. Computational Details: Molecular Mechanics Molecular mechanics calculations were carried out with the Sun-3/50 version of the AMBER 3.1 package.2s Molecular dynamics calculations were performed with the AMBER 3.1 package on a Stellar GSlOOO workstation. The initial conformation for the molecular dynamics calculations was obtained in the following manner: A lattice of 125 monomers was constructed with a volume chosen to be 30% larger than the experimental density of 0.996.26 The resulting box has a size of 31.84 A per side, allowing a minimum number of high-energy nonbonded contacts. The monomer geometry that was replicated on the lattice was chosen (25) AMBER 3.0; Singh, U C.; Weiner, P. K.; Caldwell, J.; Kollman, P. A,; University of California, San Francisco. AMBER 3.1, internal release; Singh, U. C.;Weiner, P. K.; Caldwell, J.; Seibel, G. A.; Pearlman, D. A.; Kollman, P. A. University of California, San Francisco. (26) Dreisbach, R. R.; Martin, R. A. Ind. Eng. Chem. 1949, 41, 2875, as reported i n Timmermans, J. Physicochemical Constants of Pure Organic Compounds; Elsevier: New York, 1950.

4486

The Journal of Physical Chemistry, Vol. 94, No. I I , 1990

Spellmeyer et al.

TABLE IV: Energies of Optimization and Single-Point Energies" absolute energies optimized energies

0.0 10.0 20.0 30.0 40.0 50.0 60.0 TS 70.0 80.0 90.0

RHF/3-21G// RHF/3-21G

RHF/6-31GS// RHF/6-3 lG*

-342.67 21 4 -342.67 198 -342.67 153 -342.67 086 -342.67 01 2 -342.66953 -342.66 932 -342.66932 -342.66 939 -342.66 95 1 -342.66956

-344.58 326

AMI// AMI -15.85 -1 5.80 -15.66 -15.42 -15.1 I -14.83 -14.65

RHF/6-31G*// RHF/3-21G -344.58 165

-344.58073 -344.58 105

-14.49 -14.33 -14.23

single-point energies RMP2//6-31G*// RHF/3-21G -345.64251

RMP2/6-31G*// RHF/6-31G* -345.64 320

-344.57 883

-345.63 902

-345.63 986

-344.57 896

-345.63 889

-345.63 970

single-point energies RMP2//6-31G*// RHF/3-21G 0.0

RMP2/6-31Gt// RHF/6-31G* 0.0

Zero-Point Vibrational Energies

0.0 TS 90.0

RHF/3-2IG +O. 143 902 +0.143 162 +O. 143 I86 relative energies optimized energies

70.0 80.0 90.0

RHF/3-21G// RHF/3-21 G 0.0 0.10 0.38 0.80 I .26 1.63 I .77 1.77 1.72 1.65 1.61

0.0 TS 90.0

RHFJ3-21G 0.0 -0.46 -0.45

0.0 10.0 20.0 30.0 40.0 50.0 60.0

TS

RHF/6-31G*// RHF/6-3 1G* 0.0

AMI// AMI 0.0 0.05 0.19 0.43 0.74 1.02 1.20

RHF/6-31G*// RHF/3-2 1G 0.0

I .59 1.39

1.36 1.52 1.62

1.77

2.19

2.10

1.71

2.27

2.20

Zero-Point Vibrational Energies

Ab initio values are in hartrees; semiempirical and relative energies are in kcal/mol

TABLE V: Geometries of Coplanar and Perpendicular Conformers of Anisole Optimized with AM1 bond lengths, 8, bond angles, deg parameter CI-C2 C2-C3 c3-c4 C4-C5 C5-C6 C6-C I CI-07 07-C8 C2-H9 C3-H I O C4-H 1 1 C5-H 12 C6-H 13 C8-H 14 C8-H I5 C8-H 16

O0 1.3995 1.3955 1.3925 1.3972 1.3897 1.4082 1.3823 1.4223 1.098 1 1.1002 1.0991 1.1003 1.0988 1.1171 1 . 1 196 1.1 171

90' 1.4008 1.3936 1.3954 1.3953 1.3937 1.4010 1.3917 1.4262 1.0990 1.1001 1.0997 1.1000 1.0990 1.1 169 1.1183 1.1169

parameter C 1-C2-c3 C2-C3-C4 C3-C4-C5 C4-C5-C6 C5-C6-C 1 C2-CI-C6 C2-CI-07 C6-CI-07 C1-07-C8 C 1-C2-H9 C3-C2-H9 C2-C3-H10 C4-C3-H I O C3-C4-H1 I C5-C4-H 1 1 C4-C5-H 12 C6-C5-H 12 C5-C6-H 13 CI-C6-H 13 07-C8-H 14 07-C8-H 15 07-C8-H 16 H 14-C8-H I5 H 14-C8-H 16 H 15-C8-H 16

O0 119.1 120.6 119.9 120.5 119.1 120.8 124.5 114.7 116.2 121.1 119.8 119.4 120.0 120.1 119.9 119.9 119.6 121.4 119.5 110.7 103.5 110.7 I 10.6 110.6 110.5

90' 119.0 120.3 120.2 120.3 119.0 121.1 119.5 119.2 113.1 119.9 121.1 119.7 120.0 119.9 119.9 I 19.9 119.7 121.1 119.9 110.4 103.9 1 10.4 110.9 110.8 110.3

dihedral angles, deg parameter 00

C1-C2-C3-C4 C 2-c 3-C4-C 5 C3-C4-C5-C6 C4-C5-C6-C 1 C3-C2-C 1-07 C5-C6-Cl-07 C6-C I-C2-H9 CI-C2-C3-H10 C2-C3-C4-H11 C3-C4-C5-H12 C4-C5-C6-H I3 C2-C1-07-C8 C6-C 1-07-C8 C -07-C8-H 1 1 4 C -07-C8-H 1 1 5 C 1-07-C8-H 1 6

0.0 0.0 0.0 0.0 -1 80.0 180.0 0.0 180.0 180.0 -1 80.0 -1 80.0 0.0 180.0 -61.4 180.0 61.5

90° 0.2 -0.1 0.2 -0.3 174.6 -1 74.6 -0.2 -1 79.7 180.0 -179.8 -1 79.6 90.0 -94.9 -60.5 -179.4 61.7

Rotational Barrier in Anisole

The Journal of Physical Chemistry, Vol. 94, No. 11, 1990 4481

to be the minimized structure with set to IO'. This ensures that the initial bond lengths, which are removed from the calculation in the molecular dynamics simulation, have been minimized. It also ensures that there are nonzero forces on the molecule rotation for rotation about CP in hopes of aiding the conformational equilibration process. Despite the large volume, this lattice places many atoms very near to each other, resulting in very high energies and forces. Initial attempts to minimize this structure either with steepest descent minimization or with lowtemperature molecular dynamics were futile. Therefore, a 30 000-configuration Monte Carlo run was performed on rigid bodies with the volume, temperature, pressure, and the number of particles all held constant. The nonbonded parameters chosen were the ones described in the following section. However, no charge interactions were calculated. The resulting configuration had a much lower energy and forces for the initial configuration of the molecular dynamics study. All remaining molecular dynamics simulations were performed from this initial configuration. The molecular system is modeled with the AMBER energy equation27as foIIows: N,

Vtota~=

Na

E K r n ( r n -req,12 + n=EI Kon(en- eeqJ2 + n= I

+

= [ ~ * ~ c *and ~ ] lR*ij / ~ = R*i R*j. where The isothermal-isobaric simulations were run for 30 ps with a 1.5-fs time step. The leap-frog algorithmz8 is used for the integration of the equations of motion. The SHAKE algorithm29 was used to constrain all bonds to their initial lengths. All other interactions were calculated according to the AMBER energy equation. The initial temperature, from which initial velocities are assigned, was chosen to be 10 K. Constant temperature dynamics were maintained by coupling the system to a heat bath of 300 K with a 4 K target range and a coupling parameter of 0.2 ps-'. Constant-pressure dynamics were maintained by coupling to a pressure bath with a reference pressure of 1 atm. The inverse compressibility of the system was set to 44.6 X IOd bar-' and the pressure relaxation time was set to 0.5 ps-l. Center-of-mass motion was removed from the initial velocities only. The simulations were carried out through a series of twenty lo00 time-step runs. SHAKE was applied to all bonds in the system, with a relative tolerance of 0.0005 A. A residue based nonbonded pair-list was generated and was updated every five time steps. In the current version of AMBER, the minimum image is determined on an atom-atom basis, rather than on a molecule-molecule basis. Because of this, cutoff distances must be less than 1 /2 the box size minus 1/2 the van der Waals radius (8.0 A for anisole) of the solvent molecules. Therefore, the cutoff distance for nonbonded and charge-charge interactions was chosen to be 1 1 .O A for the first 7.5 ps and 9.5 A for the remainder of the simulations. A constant dielectric function was used with c = 1 .O. All I ,4 nonbonded interactions (27) Weiner, S. J.; Kollman P. A,; Nguyen, D. T.; Case, D. A. J . Compuf. Chem. 1986, 7, 230. (28) Verlet, L. Phys. Rev. 1967, 159, 98. (29) Van Gunsteren, W . F.; Berendsen, H. J. C. Mol. Pbys. 1977, 34, 1311.

were scaled by 1/2, as were all 1,4 electrostatic interactions. Trajectories were saved every 100 time steps (150 fs) and the energy, pressure and volume were saved every 50 time steps (75 fs) for later analysis. the first 7.5 ps were used for equilibration purposes only. Results and Discussion Ab Initio Studies-Gas Phase Barrier. The geometries of the conformers of anisole at all levels of optimization are presented in Tables 1-111 (see 1 for the numbering scheme). Our optimized structures agree with the structures optimized by Shaefer and Sebastian,s and with Seip and Seip's electron diffraction study,5 but are somewhat different from the microwave structure.M Onda et al. fixed the phenyl ring in a planar conformation and fix most of the remaining parameters at values determined from STO-3G optimizations. Several features are of interest in our structures. In all cases, the internal C-C-C angles are nearly 120' and the phenyl ring remains nearly planar. The bond lengths of the benzene ring alternate slightly in the coplanar conformer with the Cl-C2,C3-C4, and c5-c6 bonds slightly shorter than the remainder (see Table I). As 4 increases, the bond alternation decreases, until all bonds are nearly equal in the 90' conformer. In the 90' conformer, the cl-c2and c& bonds are 1.3818 A. As the overlap of the p-type lone pair on the oxygen and the T orbital on the phenyl group is decreased by rotation of the methoxy from coplanar toward perpendicular, the CI-C7 bond length decreases. In the 0' conformer the C1-C7 bond length is 1.3718 A a t the 3-21G level and 1.3497 A at the 6-31G* level. These are similar to the electron diffraction value of 1.361 A,5 but significantly shorter than the microwave value of 1.399 A.30 In the 90' conformer, the C1-07 bond length is 1.3890 A in the 3-21G structure and 1.3633 A in the 6-31G* structure. At the same time, rotation of the methoxy group increases the overlap of the T orbital of the phenyl ring and the low-lying u * ~orbital. , This serves to lengthen the 07-Csbond upon rotation from 1.4356 A in the coplanar conformer to 1.4505 8, in the 90' conformer in the 3-21G structures. This effect is diminished somewhat in the 6-31G* structures where the 0,-C bond lengthens from 1.3982 A in the 0' conformer to 1.4050 in the 90' conformer. Both values are somewhat shorter than the electron diffraction value of 1.423 A.S Notice that both 0-C bond lengths are considerably shorter in the 6-31G* structures than in the 3-21G structures. The lengthening of the CI-07 bond is reproduced in the AM1 geometries, but the lengthening of the 07-C8 bond is not. The C,-O7 bond length is 1.3823 8, in the 0' conformer and 1.3917 A in the 90" conformer (see Table V). The 07-C8 bond length is 1.4223 A in the 0' conformer and just slightly longer (1.4262 A) in the 90" conformer. Steric interactions between the methyl group and H9 are thought to be the cause of the rather large c1-0& angle (124.4' in the 3-21G structure; see Table 11) and the small c6-cI-07 angle (1 16.0') in the coplanar conformer relative to the perpendicular conformer. The 3-21G value of 124.4' and the 6-31G* value of 124.5' for the cI-o7-c8 angle are both in excellent agreement with experiment ( 124°).5 These angles approach the idealized value of 120' as the methoxy group rotates, but are still some 4' different at the transition structure. As the steric interaction between the methyl group and H9 decreases, the angles approach 120'. This same steric interaction causes a rotation about the 07-Cs bond as the methoxy group rotates out of plane. Rotation must occur to allow H9 and the hydrogens of the methyl group to avoid each other. The CI-O7-C8-HI5 dihedral angle is 180' in the coplanar conformer, quickly decreases to -165.8' in the 40' conformer as the point of closest contact, and increases again to 179.6' in the 70' conformer. Relief of this steric interaction results in the decrease of the CI-C,-H9 angle from 120.9' in the coplanar conformer to 119.1' in the 90' conformer, and in the slight increase in the C I - C ~ - H I angle ~ in the 90' conformer

8:

(30) Onda, M.; Toda, A,; Mori, S.; Yamaguchi, 1. J . Mol. Sfrucf.1986, 144, 47.

The Journal of Physical Chemistry, Vol. 94, No. 11, 1990

4488

I

0

1

SLLMO

10

20

30

40

50

60

70

80

90

@ in degrees

Figure 3. Orbital energies (in electronvolts) of selected orbitals as a function of rotation about a. All energies are from fully optimized structures at the RHF/6-3IG* level.

( 1 19.1') relative to the coplanar conformer ( 1 18.1').

This steric interaction is also partly responsible for the rather open CI-07-CB angle of 120.8' at the 3-21G level ( I 19.8' with 6-31G*) in the coplanar conformer. Both of these are in agreement with the electron diffraction value of 120.0°,5but differ significantly from the microwave value of 113.8'. Onda et al. admit that this value is dependent upon the choice of structural parameters used in their study.jO Although this angle decreases as 4 increases to 116.2" in the 3-21G 90' conformer ( 1 15.5' at the 6-31G* level), it is large for an sp3-hybridized oxygen atom, indicating that considerable sp2 character exists. The decrease in this angle is more likely a result of the increase in the overlap of the T orbital and the u*O,-Cs orbital rather than the relief of steric interactions. The methoxy group tends to twist toward the phenyl ring as the overlap increases. These geometric changes are paralleled in the orbital energies. Figure 3 presents the orbital energies for selected orbitals as a function of the rotation about 4. Energies are from the 6-31G* optimized geometries of the coplanar, transition structure, and from the perpendicular conformation. In all cases, the orbitals on the oxygen have been omitted for clarity. The SHOMO and the LUMO cannot mix with the oxygen lone-pair orbital or with the uws or u*orc8 orbitals. Thus, the energy of these two orbitals remains relatively unchanged in all three structures. The energy of the oxygen nonbonded pair rises rapidly from the coplanar to the transition structure, and less rapidly between the transition structure and the perpendicular conformation. This is expected. The stabilizing two-electron interaction with the S L U M 0 is decreased upon rotation of the methoxy group out of plane, while the destabilizing four-electron interaction with the HOMO is decreased. Rotation beyond 60' does little to further decrease overlap, and thus, the energy rises less rapidly. As is seen in the figure, the energy of the HOMO decreases between the coplanar and the perpendicular structures. This is a reflection of a decreased four-electron destabilizing interaction with the nonbonded oxygen lone pair and an increase of the stabilizing two-electron interaction with the orbital. Though of no major significance to this study, the donation of electrons from the uCH orbital into the U*CO bond strengthens the C8-HISbond and decreases the O7-C8-HISangle somewhat. The CS-Hl5 bond (1.0776 A at the 3-21G level, and 1.0798 A with 6-31G*) is shorter than are the C8-HI4and C8-Hl6 bonds (1.0831 A at the 3-21G level and 1.0851 8, with 6-31G*), and the 0,Cs-HIs angle of 105.6' (106.2', 6-31G*) is almost 6' smaller than the 07-C8-Hf4 and O7-CS-Hl6 angles of 11 1.4' ( 1 1 IS', 6-31G*). Good overlap of the ucH orbital with the uc,+ orbital is attained only when the hydrogen is anti to the phenyl ring. This trend is maintained throughout the rotation of C#I from 0' to 90' (see Table II), and is also found in the AMI geometries (see Table VI.

Spellmeyer et al. A very interesting phenomenon is that, as the methoxy group rotates, the oxygen moves out of the plane of the phenyl ring, yet the phenyl ring itself remains nearly planar. As measured by the sum of the dihedral angles C3-c2-cI-07 and cs-C6-c,-07, the degree of nonplanarity reaches a maximum of 3.7', in the 40' conformer. The maximum is 2.2' when the nonplanarity is measured by the sum of the C2-CI-07-Cs and c6-c~-o7-cs dihedral angles. Klessinger found similar results, but his data indicate that the oxygen is 3.1 'out of plane in the 90' conformer.12 This is likely a result of the increasing steric interaction of the methyl group at Cs and H9. The resulting force would push the oxygen away and down from the methyl group, and consequently, out of the plane of the phenyl group. The geometries of all conformers at the RHF/3-21G and at the RHF/6-3 lG* levels are very similar. As is generally true, the bond lengths in the 6-31G* structure are slightly longer than in the 3-21G structure.)' The notable exception is the 0,-Cs bond length that is 0.043 8, shorter in the 6-31G* transition structure than in the 3-21G structure. Presumably, this is due to the better description of the T-u*o,+ orbital interaction with the larger basis set. Likewise, the CI-07-C8angle is 0.8" smaller in the 6-31G* structure than in the 3-21G structure. Also, not surprisingly, the 6-3 lG* transition structure occurs 1.2' closer to the coplanar structure than in the 3-21G calculations. Consequently, the steric interaction between the methyl group and H9 are more significant, and the methyl group is 3.3' less staggered than in the 3-21G transition structure. Otherwise, the 6-31G* structures compare favorably with those calculated at the lower level. It should be noted that the 6-31G* optimized coplanar conformation had the additional constraint that the phenyl ring be held planar. This is not a severe approximation, as the degree of nonplanarity at the 3-21G level is at most 0.2'. Unlike the a b initio profiles, the energetic profile at the AMI level is definitely 2-fold in nature (see Table I V and Figure 1). The energy rises slowly from the coplanar conformer to the perpendicular transition structure. Semiempirical methods are known to reproduce rotational barriers quite poorly,32and it seems as if this case is no exception. In contrast to the AM1 profile, the 3-21G profile has a sharp rise between 10' and the transition state, while remaining relatively flat between the transition and the 90' conformer. As has been suggested earlier, there must be a substantial 4-fold barrier to rotation about the C1-07 bond.s At the RHF/3-21G level, the transition state occurs at 4 of 60.9' and is 1.8 kcal/mol higher in energy than the coplanar structure. This is only 0.2 kcal/mol higher in energy than the 90' conformer. Inclusion of the zero-point vibrational energies for all three decreases the barrier to 1.2 kcal/mol and the 90' conformer remains 0.1 kcal/mol lower than the transition state. Thus, the overall shape of the curve does not change upon inclusion of zero-point vibrational energies. The relative energies calculated at the RHF/631G*//RHF/3-21G level are nearly coincident with the RHF/3-21G values, differing by no more than 0.05 kcal/mol. Not unexpectedly, the energy of the fully optimized transition structure (4 = 59.7') at the RHF/6-31G* level is 1.6 kcal/mol higher that that of the coplanar structure and is 0.2 kcal/mol above the 90' conformer. An interesting difference appears when the energies are calculated at the RMP2/6-3lG*//RHF/3-2lG level. With this level of theory, the "transition state" is 2.2 kcal/mol higher than the coplanar conformer, while the 90' conformer is 2.3 kcallmol higher. Thus, it would appear that the 90" conformer is a maximum on the energy profile, rather than a minimum at this level. This too is not surprising. Nobes et al. have found similar results with meth~xyethene.~~ The difference in energy at the RMP2/6-31G*//RHF/6-31G* level is also 0.1 kcal/mol. Further optimizations at the MP2/6-31G* level are needed to (31) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Theory; Wiley: New York, 1986. (32) For example, the eclipsed versus staggered energy difference in ethane is 1.2 kcal/mol with the AMI method. (33) Nobes. R. H.; Radom, L.; Allinger, N. L. J . Mol. Srrucr. (THEOCHEM) 1981,85, 185.

The Journal of Physical Chemistry, Vol. 94, No. 11, 1990 4489

Rotational Barrier in Anisole H

0.0812

I CT-

'

0.4183 0.1745 H

H 0.1745

0.1408 H

-0.1905

H 0.1453

2.5

HC

os I I

I

HC

'HC

H 0.1408

II

"."0

HC

Figure 4. The partial atomic charges and the AMBER atom types used in the molecular mechanics and dynamics studies.

clarify the nature of the perpendicular conformation. These are currently beyond our resources. Despite the fact that there is still some ambiguity about the presence of a second conformer based on our highest level calculations, the following conclusions are justified: ( I ) There is a substantial 4-fold component to the rotational barrier about the CI-07 bond. (2) This 4-fold component apparently is invariant to the level of the ab initio calculation and is thus likely to be steric in nature, rather than electronic. (3) The 2-fold component is dependent upon basis set and is sensitive to post-SCF methods. Thus, it is likely electronic in nature, rather than steric. (4) The barrier to rotation is nearly 2.1 kcal/mol in the gas phase. (5) The energy profiles at the AMI level are significantly different than those with the a b initio methods. (6) The energy profiles with the ab initio methods are all qualitatively similar, with an almost flat potential between 4 = 50' and 90' and a rapid rise in energy between C$ = Oo and SOo. The rapid rise in the energy can be attributed to a combination of events. As the methoxy group rotates out of plane, the following must take place: (1) Steric interactions between H9 and the methyl group initially increase. To avoid this, the methyl group rotates about 14' from a staggered conformation. (2) The oxygen moves out of the plane of the phenyl ring. (3) The stabilizing overlap of the A* orbital and the oxygen in-plane lone-pair orbital decreases. These all reach a maximum between 40 and 60°. Meanwhile, the following occur between 60° and 90° and serve to decrease the energy, at the R H F level, and to maintain the flat barrier at the MP2 level: ( 1 ) The stabilizing A - O * ~ , + overlap increases. (2) The methyl group reestablishes a staggered conformation. (3) Further dimunition of the steric interaction between H9and the methyl group. (4) Relief of the angle strain the the CI-o7-c8 and c 6 4 7 - c 8 bonds. The existence or absence of the second minimum hinges upon the subtle interplay of these events. the few tenths of a kcal/mol differences in the R H F and the MP2 energies is probably due to the overestimation of the contribution of the u-u*o,+ overlap to stabilization of the energy at the R H F level, relative to the destabilization brought about by decreasing the overlap of the oxygen lone pair with the A* orbital. The exact role of basis set size limitations and the lack of electron correlation during optimization add to the uncertainty in the relative energy of the perpendicular form relative to the transition state. Parameter Development. Partial atomic charges were determined at the RHF/6-31G*//RHF/3-21G level with the ESP method of Singh and K ~ l l m a using n ~ ~ the perpendicular conformer as the basis. The first choice of geometry for the charge fit is the coplanar conformation, as the molecule will assume that conformation most frequently. However, the charge distribution in this conformation is rather unsymmetrical, which could lead to unsymmetrical conformation populations in the simulations. Thus, the charges were symmetrized by averaging the values for, e.g., the C2 and (26 atoms, etc. The resulting charges are shown in Figure 4. The 6-31G* basis set was chosen because it has been shown to provide reliable charges and is consistent with the charges for a new force-field under development at UCSF.34

60

120

180

240

300

360

@ Figure 5. A plot of the relative energies of the rotational barrier in anisole as obtained with AMBER in comparison with MP2/6-31G*// RHF/6-31G* single point energies.

The intramolecular parameters for anisole were taken from Grootenhuis and Kollman.'s Several nonstandard parameters were also developed for this simulation. We can assume that the 4-fold component in the rotational barrier is invariant with respect to the 2-fold component. Therefore, the torsional parameters for the 07-C8 bond were initially developed to reproduce the rotational barrier obtained at the RHF/3-21G level. The resulting energy profile deviates from the ab initio profile by not more than 0.1 kcal/mol. This is must lower than the generally accepted value of the error in potential functions of this type. On the basis of the MP2/6-31G* calculations, we know that the 2-fold component of this barrier is too low. The 2-fold component was increased to reproduce the energy profile at the MP2/6-31G*//RHF/631G* level. The V2/2 and V4/2 terms are 2.60 and 0.325 kcal/mol, respectively. It should be noted that these are not directly related to the experimentally determined V2 and V4 terms, as many other energy components are included in the rotational barrier. The final barrier is shown in Figure 5. There is a very slight asymmetry in the barrier with respect to the perpendicular conformation. Attempts to locate the source of this asymmetry proved fruitless. We decided that the error introduced in this manner was insignificant because the magnitude of the largest difference in energies is 0.1 kcal/mol and that this occurs prior to reaching the plateau in the barrier. Attempts were made to reproduce the ab initio potential using only the VI.V,, and V, terms in the Fourier series for the dihedral angle terms. Though not exhaustive, these terms proved to be insufficient to reproduce the rotational barrier. This phenomenon is unusual but points up a deficiency in truncation of the torsional potential at only three terms. It has become clear that the standard AMBER nonbonded parameters for carbon and hydrogen do not accurately reproduce the liquid properties for hydrocarbon^.^^ Therefore, nonbonded parameters for carbons and hydrogen were taken from liquid simulations of hydrocarbons. These have been shown to reproduce the structure and density of benzene and for the alkane series of methane to butane.34 The parameters for the ether oxygen were taken from Jorgensen's OPLS parameter^.^^ The resulting parameter set is shown in Table VI. Molecular Dynamics-Liquid Phase. The molecular mechanical parameters described above were used to perform a 30-ps molecular dynamics simulation of liquid anisole. A box of 125 anisole monomers was used with periodic boundary conditions at 300 K and 1 atm of pressure in the NPT ensemble. The system was equilibrated for 7.5 ps. That the system was in equilibrium is shown by the convergence of the total energy and of the volume. These are shown in Figures 6 and 7, respectively. A further 22.5 ps of simulation was used for the analysis of the run. Although the system is relatively small, it is quite adequate for the determination of the conformational population about the Cl-O7 bond. The choice of nonbonded parameters seems to be (35) Jorgensen, W . L.; Tirado-Rives, J. J . Am. Chem. SOC.1988, 110,

(34) Spellmeyer, D. C.; Kollman, P. A., unpublished results.

1657

4490 The Journal of Physical Chemistry, Vol. 94, No. 11, 1990

Spellmeyer et al. TABLE VI: AMBER All-Atom Parameters for Anisole and Other Aryl Alkyl Ethers atom name CA CT HC

atomic mass 12.01 12.0 I 1.oo

os

16.00 Bond Parameters

650 . . . . . . . .' . 2 5 1d.o l i s 1 i . o i j . 5 26.0 2i.5 2i.o 2i.5 3d.o Time (in picoseconds)

Figure 6. The total energy (in kcal/mol) as a function of time (in ps) for the molecular dynamics study.

bond type CA-HC CA-CA CT-HC CA-OS CT-OS

force const 340.0 469.0 331.0 300.0 320.0

equilibm length 1.08 1.40 1.09 1.40 1.41

Angle Parameters angle type C A-CA-C A CA-OS-CT CA-CA-OS CA-CA-HC OS-CT-HC HC-CT-HC

7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0 27.5 30.0 Time (in picoseconds)

Figure 7. The volume (in A') as a function of time (in ps) for the

dihedral type X-CA-C A-X X-C A-OS-X X-CA-OS-X X-CT-OS-X

equilibm angle 120.0 1 13.0 120.0 120.0 109.5 109.5

Dihedral Angle Parameters division peak phase factor height angle 4 5.3 180.0 2 2.600 180.0 2 0.325 180.0 3 1.15 0.0

periodicity 2 2 4 3

Improper Dihedral Angle Parameters

molecular dynamics study.

quite good; the calculated average volume per particle is 191 A3, while the experimental value is 179.3 The energy fluctuations as a function of time are presented in Figure 6, and the volume fluctuations are shown in Figure 7. It is clear that these properties have reached an equilibrium value. Because the torsional barrier was derived from very high level ab initio calculations, and the volume of the liquid is in good agreement with experiment, we believe we have an adequate representation of the pure liquid with which to assess the conformational profile of the methoxy rotation in anisole. Averaging of the two dihedrals C2-C,-O7-C8and C6-C1-07-C8 for the 22.5 ps of data collection is shown in Figure 7. It is clear that the populations of the conformations of the methoxy group are being sampled inadequately and that they are not yet in equilibrium at the end of 30 ps. Further simulations were deemed too expensive to continue the simulation beyond this point. Therefore, the following results are qualitative in nature and are not to be taken as quantitative. The expected population density based on a Boltzmann distribution of the gas-phase energy at the 3-21G level and at the MP2/6-3lG*//RHF/6-3lG* level is shown in Figure 8. Notice that the conformation population of dihedreal angle, which was initially near 180' the C6-C,-07-C8 in the starting conformation, is very similar to the predicted gas-phase population in the vicinity of 180'. Likewise, the conformation population of the C2-C,-07-C8 dihedral angle is similar to that of the gas-phase population near 0". The intermediate angles from about 40' to 140' and from 220' to 320' are populated more in the liquid phase than in the gas phase. Despite the fact that the methoxy group is not quite equilibrated, it is clear that the barrier in the liquid is very similar to the gas-phase barrier. It is not very likely that the liquid barrier is twice that of the gas phase. The estimate of the barrier to rotation based upon the dihedral angle population densities from the molecular dynamics simulation of 1.5 kcal/mol is substantially lower than the gas-phase estimate of 3.1 kcal/mol.8 This value is also lower than that obtained from the gas-phase estimate of 3.1 kcal/mol.8 This value is also lower than that obtained from the gas phase barrier of 2.3 kcal/mol from which the parameters were developed.

force const 85.0 46.5 70.0 35.0 35.0 35.0

dihedral type CA-CA-CA-CA

peak height 20.0

phase angle 180.0

periodicity 2

Lennard-Jones Parameters atom HC

os

CT CA

R* 1.4592 1.6837 1.9082 1.9082

L

0.015 0.170 0.1 IO

0.086

0.07 I

0 Figure 8. The population density of the dihedral angles

+

obtained from molecular dynamics simulation in comparison with the Boltzmann distribution as predicted from the gas-phase rotational barrier as predicted from the AMBER gas-phase barrier.

The most plausible explanation of the apparent contradiction between previous experiments and our liquid results is that the experimental results were obtained under the assumption that the barrier is 2-fold in nature. Recent studies by Schaefer would indicate that the barrier would be considerably diminished if the 4-fold nature of the barrier was taken into account in the evaluation of the experimental data.8 These authors state that the nmr data for liquid anisole would indicate a barrier of nearly 9.6

J . Phys. Chem. 1990, 94, 4491-4493 kcal/mol if a 2-fold barrier is assumed. However, if the barrier is assumed to be 4-fold in nature, the rotational barrier estimate from experiment would be in the range of 3.1 kcal/mol.8 Our a b initio calculations would suggest that there is a substantial 4-fold component to the barrier in the gas phase. Our molecular dynamics simulations would indicate that the barrier in the liquid phase is not different from the gas phase, in qualitative agreement with the results of Shaefer.8

Conclusion In this study, we have presented definitive high-level ab initio calculations showing the energy profile of the rotational barrier in anisole. The potential is composed of a large 2-fold term and a small, but very important, 4-fold term. The overall rotational barrier is about 2.1 kcal/mol, and is somewhat sensitive to the level of optimization. At all levels of theory, the coplanar structure is the minimum energy arrangement. However, it appears that there is not a second stable structure at the MP2/6-31G* level of theory, in contrast to the results at the RHF/6-31G* level. In either case, the barrier between 60" and 90" is very flat (