J. Phys. Chem. 1995,99, 5873-5882
5873
Ab Initio Calculations on Small Molecule Analogues of Polycarbonates H. Sun, S. J. Mumby, J. R. Maple, and A. T. Hagler" Biosym Technologies, Inc., 9685 Scranton Road, San Diego, Califomia 92121 Received: Janualy 11, 1994; In Final Form: December 29, I994@
In order to develop an ab initio force field for polycarbonates, extensive quantum mechanical calculations were carried out on several model compounds: carbonic acid, methyl and dimethyl carbonates, phenyl carbonate, and 2,2-diphenylpropane. The calculations were performed with Hartree-Fock (HF), MgllerPlesset second order perturbation (MP2), and density functional theory (DFT). Full geometry optimizations were performed to characterize the global and local minimum energy structures and transition states of internal rotations. These geometry optimizations reveal large differences in valence coordinates between different conformational states, indicating the need to take molecular flexibility into account when calculating properties of these systems. The equilibrium structures and energies are discussed in terms of substituent group effects and steric crowding. Finally, atomic partial charges were calculated via three methods: Mulliken population analysis, calculation of atomic polar tensors, and by fitting the electrostatic potential surface. The latter calculations begin to uncover weaknesses in the isotropic partial atomic charge model and quantitatively demonstrate the possible importance of atomic dipoles and even quadrupoles in these systems.
I. Introduction Computer simulation and molecular modeling of organic and biological molecules have been performed for many years, and have provided considerable insight into the structural order and interactions within these systems at an atomistic level. More recently, an increasing number of such investigations of polymeric systems have been reported in the literature. This is due, in part, to the extensive utilization of advanced experimental techniques such as NMR and IR spectroscopy to probe the properties and structures of polymers that have created a significant challenge for computational chemists. In addition, the development of faster computers affords the possibility of performing large scale, quantitative molecular dynamics (MD) and Monte Carlo (MC) simulations to investigate polymeric systems. However, a fundamental underpinning of all such simulations is a realistic molecular potential energy function that accurately accounts for both inter- and intramolecular interactions. As a part of a serial investigation, polycarbonates are the focus of the study described here. One of the most important polycarbonates is bisphenol A polycarbonate (BPAPC; see Figure 1). The outstanding optical and mechanical properties of this polymer have drawn much attention during the past one or two decades, and many experimental investigations have been conducted to probe its microscopic structure, dynamics, and properties. One of the most active areas of these studies has been the attempt to understand the phenyl ring motions and how these motions affect the properties of BPAPC. In addition to several dynamic mechanical spectroscopy (DMS) investigations,' solid state NMR studies by Schaefer et aL2 and Spiess3 suggest that the dominant motion of BPAPC at room temperature is a 180" flipping of the aromatic rings about their C;, symmetry axis. These flips are superimposed on small amplitude ring oscillations and main chain wiggles. This suggestion has been supported by recent I3C NMR studies conducted over a range of temperat~res.49~ The results of the low temperature NMR ~ t u d i e sare ~ . ~consistent with X-ray structural analysis6s7 for the diphenyl carbonate of 2,2-bis(4-hydroxyphenyl)propane, @
Abstract published in Advance ACS Abstracts, March 15, 1995.
Figure 1. Bisphenol A polycarbonate (BPAPC).
a crystalline analogue of BPAPC. Variable temperature FTIR studiess of BPAPC have also provided valuable information regarding the dynamic properties of this polycarbonate. The cooperative motion in polycarbonates has been studied recently? and it has been found that the motion involves several repeat units. In addition to the large amount of experimental work, many theoretical investigations have been reported. As early as the 1960s, Williams and FloryIo studied the conformational energies of polycarbonates by using an empirical force field. Tonellil performed conformational energy calculations and found that the phenylene ring motion occurs quite easily. This result was obtained later by Bendler,I2 using MNDO calculations. Johansen, Rettrup, and JensenI3 calculated electrostatic potentials and deformation densities of methyl carbonate. In recent years, Bicerano and Clarki4conducted semiempirical MNDO, AM1, and approximate ab initio PRDDO (partial retention of diatomic differential overlap) calculations for several small-molecule analogues of BPAPC and the poly(ester carbonate) of bisphenol A (BPA). Laskowski et al.I5 performed ab initio Hartree-Fock (HF) calculations on small model molecules of BPAPC. More calculated results of this group's work can be found in a recent review article written by Jaffe, Yoon, and McLean.I6 Sung et al.I7 have also conducted CNDO calculations on similar molecular analogues. Labrenz and Schroer'* carried out a thorough analysis of geometries, energies, and dipole moments of the conformers of carbonic acid and its methyl, ethyl esters by using various semiempirical and ab initio calculation methods, and experimental dielectric measurements. Very recently, Hutnik, Argon, and SuterI9 developed a classical force field for BPAPC based on the latest published theoretical and experimental results. However, owing to the limitations imposed by computational resources, most of these calculations are empirical or semiempirical in nature, from which controversial structural and energetic results have arisen. The ab initio
0022-3654/95/2099-5873$09.00/00 1995 American Chemical Society
5874 J. Phys. Chem., Vol. 99, No. 16, 1995
calculation^^^^^^^^^ have been primarily focused on minimumenergy structures, and the information obtained from such calculations is often insufficient to develop a general force field. In fact, most force fields described above invoke a rigid valence geometry, an approximation we shall address in some detail. Several “relaxed” geometry empirical force fields have been developed for molecular modeling of organic and biological systems, such as CVFF?O CHARMM,21AMBER?2 MM3,23and CFF93.24-26 The commonly used approach to develop such a force field has been to start with the experimental geometries and vibrational frequencies for a set of related model molecules and to find force field parameters that provide the best fit to these data. However, this approach suffers from the sparsity of available experimental data. The quality of the potential function is generally dependent on the amount of available experimental data: the more experimental data that are available, the better the force field may be determined. Alternatively, one can enhance the available experimental data by combining it with an ab initio force field.24-28 To develop such an ab initio force field, we start with quantum mechanical calculations on several small model molecules, which sample the molecular interactions in which we are interested. The calculations should provide detailed and accurate information about molecular structures, conformational energies, electrostatic properties, and vibrational frequencies. In the current investigation, more extensive ab initio calculations have been performed on small molecule analogues of polycarbonates. These studies extend the previous work in terms of the size of the model molecules considered, the number of conformations examined, the basis sets employed, and the level of theory. The results provide direct insight into the properties of BPAPC, and the data are the basis for the subsequent development of a polycarbonate force field. Also, they provide valuable information on the conformational energies and spectroscopic properties of small organic molecules such as dimethyl carbonate, which have been the subject of extensive experimental work for several d e c a d e ~ ? ~ - ~ but l on which theoretical calculations are sparse at best. In the second stage, a force field is developed based on these calculation^.^^ In addition, we take this opportunity to test the applicability of density functional theory to these small molecule analogues. Density functional theory (DFT) was developed in the early 1960’s (Kohn and Sham, 1965) for solid state physics and since then it has slowly gained acceptance among chemists owing, in part, to the recent advances in the so-called nonlocal gradient corrected functionals (for example, see Perdew, 1991). Since the exchange and correlation energies are included within the density functional formalism, it has the potential to provide accurate results at a fraction of the cost over traditional postHartree-Fock methods such as MP2 and CI. 11. Calculation Details All ab initio calculations reported in this paper were conducted on Silicon Graphics SGI-4D240-GTX and IBM RISC6000/550 workstations using the ab initio software package TURBOMOLE.32 In keeping with our research most calculations for later force field development were performed using the Hartree-Fock method with the 6-31G* basis set.33 The direct SCF method was used in all calculations. MBller-Plesset second order perturbation method (MP2)34with the 6-31G* and 6-311G** basis sets35 was used on smaller molecules in order to investigate the impact of using larger basis sets and electron correlation. For comparison, calculations using density functional theory36with the recently developed ACM functional37in conjunction with DZVP basis functions3*(DFT-
Sun et al. (ACM)/DZVP) were carried out. These calculations are particularly interesting for large molecules because recent studies39show that the DFT results are generally comparable with those of MP2 which becomes computationally expensive as the size of molecule increases to the level of monophenyl carbonate or diphenylpropane. Full geometry optimization using analytical gradients was performed to characterize global and local minimum-energy structures as well as transition-state structures. Total energies were calculated, as were the analytical gradients and Hessians (first and second derivatives of the total energies) in order to test the optimized structures and to characterize the energy surface of these molecules for the derivation of analytical force fields.27 The electrostatic properties such as atomic monopole moments (atomic charges), molecular dipole moments, and their derivatives were calculated. Atomic partial charges were calculated from the €IF wave functions by using three methods: Mulliken population analysi~,4~,~’ by fitting to electrostatic p o t e n t i a l ~ , 4 ~and - ~ ~analysis of atomic polar tensors48(derivatives of molecular dipoles). Five molecules were selected to explore the energy surface of polycarbonates: carbonic acid, methyl and dimethyl carbonate, phenyl carbonate, and 2,2-diphenylpropane. These model compounds include the most important molecular structures or fragments comprising polycarbonates. The illustrations of the optimized structures of these molecules and the numbering of atoms are given in Figures 2 to 4. Among these molecules, 2,2-diphenylpropane is especially important for understanding the phenyl ring motions in BPAPC, because it is representative of the smallest segment that contains the two phenyl rings. Carbonic acid and its methyl and phenyl derivatives were chosen to probe substitution effects, intemal barriers to rotation about backbone bonds, C1-03 and 03-C5, and the coupling between two functional groups separated by the OCOO group. Dimethyl carbonate is one of the few stable small carbonate molecules for which experimental data are available in the literature. Thus, this molecule provides an opportunity to compare the theoretical results with experimental data. Diphenyl carbonate would be another model molecule of interest, but practically it is extremely expensive to perform accurate ab initio calculations for this molecule. In addition, based on the comparison of results calculated for methyl and dimethyl carbonates, the coupling between the two phenyl groups in diphenyl carbonate is expected to be weak. Evidence to this effect was obtained by Laskowski, et a l . I 5 who calculated the trans-trans and trans-cis energy difference for diphenyl carbonate with the STO-3G basis set, and reported that the energy difference is only 0.04 kcal/mol lower than that calculated for monophenyl carbonate. In this study, we performed similar calculations at the HF/6-31G* level and obtained 0.06 kcal/mol difference between these two conformers. Thus, there is good evidence to support the assumption that the interactions between the phenyl and carbonate groups are represented well by those in the monophenyl carbonate, and that the coupling between the two phenyl groups separated by the OCOO group can be represented solely by nonbonded interactions.
111. Results and Discussions We first studied the consequence of using different basis sets and computational methods. Since Hartree-Fock calculations with smaller basis sets such as 3-21G and 4-21G have been reported in the l i t e r a t ~ r e ,our ~ ~ calculations ,~~ started with HF/ 6-31G* which in most cases yields reasonable geometric and energetic results for organic compounds.33 MP2 calculations with the 6-31G* basis set were performed, and the comparison
J. Phys. Chem., Vol. 99, No. 16, 1995 5875
Ab Initio Calculations on Polycarbonates Carbonic acid
I1 H6 I
Methyl carbonate
111
I
Dimethyl carnonate
I
111
Figure 2. The optimized structures of carbonic acid and its methyl derivatives. Trans-trans, trans-cis, and cis-cis conformers are shown.
0'
w
1I
I
I1
cm cm
Figure 3. Four optimized structures of monophenyl carbonate. Atoms are numbered in structure I only, but apply to all conformers. Structure I is the minimum energy structure; the phenyl ring is twisted about the 03-C5 bond out of plane by about 64". Structures I1 and I11 were calculated by constraining the phenyl ring to be twisted by 90" out of plane and in plane, respectively. Structure IV is one of the trans-cis conformers.
with HF/6-3 lG* results provides direct information about the influence of electron correlation. MP2 calculations with a larger basis set, 6-31 1G**, are more reliable in general. Unfortunately, these calculations were limited to small molecules (carbonic acid, dimethyl carbonate) in this study owing to limitations of computational resources. Since DFT results have been found to be comparable with those of MP2 calculations,40DFT(ACM)/ DZVP calculations were conducted both for the information they could provide, and to further test their applicability on systems where all methods could be applied. In Table 1, we list the optimized geometrical parameters for carbonic acid, dimethyl carbonate, monophenyl carbonate, and 2,2-diphenylpropane.The molecular structure of dimethyl carbonate has been resolved
cH3 cH3
111
IV
Figure 4. Four optimized structures of 2,2-diphenylpropane. Atoms are labeled by numbers in structure I only, but apply to all of the conformers. The minimum energy structure I has C2 symmetry, with the phenyl rings twisted in an opposite sense from the C5-Cl-C4 plane by 50'. Structure I1 has C, symmetry; one ring is in the plane while the other is twisted 90' out of the plane. Structures I11 and IV exhibit CzVsymmetry; both rings are either twisted 90' or 0'. based on gas electron diffraction data.31 For comparison, the experimental data are also given in the table. It is apparent that there are pronounced differences among the results of the four calculation methods. In particular, HF/ 6-3 lG* structural parameters of carbonates are significantly different from those of MP2 and DFT calculations. The calculated bond lengths of carbon-oxygen are too short at the HF/6-31G* level. Inclusion of the MP2 correction with the same basis functions (6-3 1G*) elongates these bonds by 0.02 to 0.03 A. However, with a larger basis set, the MP2/6-3 11G** bond lengths are reduced slightly, which yields the best agreement with experimental data, as shown in this table.
5876 J. Phys. Chem., Vol. 99, No. 16, 1995
Sun et al.
TABLE 1: Comparison of Optimized Geometric Parameters of Model Molecules of Polycarbonate, in Their Most Stable Conformation@ HF/6-31G* MP2/-31G* MP2/6-311G** bond length (A) c1-02 C1-03 03-H5 bond angle (") 02-C1-03 C 1-03-H5 bond length (A) c1-02 C1-03 03-C5 bond angle (") 02-C1-03 C1-03-C5 bond length (A) c1-02 C1-03 C1-04 03-C5 04-H6 bond angle (") 02-C1-03 02-C1-04 C1-03-C5 C1-04-H6 C3-05-C7 C3-05-C8 torsion (") C7-C5-03-C( 1) bond length (A) Cl-C2 Cl-C4 bond angle (") C3-CZ-Cl C4-Cl-C5 C1-04-C6 C1-04-C8 torsion (") C4-Cl-C5-C9 out-of-plane ( o ) b C1-04-C6-C8
125.1 105.4
1.217 1.345 0.977 125.9 105.4
1.205 1.339 0.963 125.9 104.9
1.212 1.340 0.968 125.6 106.2
125.4 116.4
1.219 1.344 1.440 126.2 113.3
1.207 1.337 1.431 126.2 113.2
1.213 1.339 1.434 126.0 114.4
Monophenyl Carbonate 1.185 1.319 1.318 1.385 0.951
1.214 1.349 1.352 1.405 0.977
1.208 1.344 1.345 1.399 0.968
126.9 124.7 119.8 107.5 120.9 117.3
127.6 125.7 116.7 105.2 121.3 116.4
127.5 125.4 118.6 106.0 121.6 116.5
63.5
59.1
59.2
2,2-Diphenylpropane 1.543 1.541
trans-trans trans-cis cis -cis trans-trans trans-cis cis-cis
MP2/6-31G*
MP2/6-311G**
DFT
0.0 2.17 13.04
0.0
Carbonic Acid 0.0 2.19 13.59
0.0 2.44 14.19 Dimethyl Carbonate 0.0 0.0 3.80 3.36 24.53 26.73
1.98 12.29
0.0 3.16 21.79
Monophenyl Carbonate
Dimethyl Carbonate 1.190 1.313 1.417
HF/6-31G*
DFT
Carbonic Acid 1.188 1.315 0.951
TABLE 2: Comparison of Conformational Energies (in kcaVmo1F
1.541 1.536
107.3 109.6 119.3 123.1
107.7 109.4 119.3 122.9
50.5
50.8
-1.4
-1.6
Illustration of the structures and the numbering of atoms are given in Figures 2-4. The out-of-plane angle of a-b-c-d is defined as the angle between the a-b bond and the b-c-d plane.
Comparison of the data in this table indicates that, as in previous systems, the DFT calculations are in excellent agreement with the MP2 methods. For the molecules we studied, the results of DFT(ACM)/DZVP fall in the region between those of MP2/ 6-31G* and MP2/6-311G** calculations. Conformational energies calculated using the different methods are summarized in Table 2. Comparison of the data given in this table shows that the MP2 correction with both 6-31G* and 6-311G** basis sets has a small effect on the conformational energies for the molecules we calculated, although the MP2 correction generally yields slightly lower conformational energies than the HF calculations. It is of interest to point out that the DFT energies are again generally in better agreement with the correlated HF results although in this case they exhibit a significant trend toward lower values. In order to develop an accurate force field representing the potential energy surface, it is important to model the changes of internal coordinates (bond lengths, angles, etc.) for different conformations. In Tables 3-7, we give the optimized valence coordinates for the calculated conformations. Only HF/6-3 lG*
optimized translskewed wandplanar cislskewed
I (optimized) I1 (0,90) 111 (90,90) IV (0,O)
0.0
0.0
0.04 1.01 2.81
0.20 1.06 2.49
0.00 0.05 0.37 2.08
2,2-Diphenylisopropane 0.0 2.17 3.93 17.21
0.00 1.80 3.50 17.53
Illustration of these structures and numbering of atoms is given in Figures 2-4.
TABLE 3: Comparison of Structural and Energetic Parameters Calculated for Conformers of Carbonic Acid, at the HF/6-31G* Level" structure bond length (A) c 1-02 C1-03 C1-04 03-H5 04-H6 bond angle (") 02-C1-03 02-Cl-04 H5 - 0 3 -C 1 H6-04-C 1 nonbonded distances (A) H5. * -02 H6. * - 0 2 H6. -03 H5. * -04 H5- * -H6 conformational energies (kcal/mol) HF MP2
trans-trans 1.188 1.315 1.315 0.95 1 0.95 1 125.1 125.1 107.9 107.9
trans-cis 1.180 1.316 1.333 0.950 0.952 124.4 124.6 108.3 110.8
2.291 2.291 2.966 2.966 3.670
2.299 2.975 2.993 2.159 3.096
0.0
2.44 2.14
0.0
cis-cis 1.172 1.333 1.333 0.947 0.947 122.1 122.1 115.0 115.0 2.996 2.996 3.340 3.340 2.027 14.19 13.51
The MP2 energies are calculated for the HF/6-31G* optimized structures.
data are given in these tables. However, our calculations indicate that the changes in these values for different conformations are not sensitive to the methods employed, although the absolute values are relevant to the methods used. Carbonic Acid and Its Methyl Derivatives. Three optimized structures of each of the molecules were examined. Following the literature, we refer to these structures as trunstrans ( I ) , trans-cis (11),and cis-cis (HI) conformers; they are schematically illustrated in Figure 2. In addition to geometric and energetic parameters given in Tables 1 and 2, the optimized geometric parameters for different conformers are summarized in Tables 3, 4,and 5 for carbonic acid, methyl carbonate, and dimethyl carbonate, respectively. Focusing on the geometrical parameters, several interesting points become apparent. That the sp2 hybridized carbon atom C1 forms three CJ bonds with two of the adjacent oxygens, and an additional JC bond with the carbonyl oxygen, is clearly shown by the short length of the C1=02 bond. However, the p
Ab Initio Calculations on Polycarbonates
J. Phys. Chem., Vol. 99, No. 16, 1995 5877
TABLE 4: Comparison of Structural and Energetic Parameters Calculated for Conformers of Methyl Carbonate at the HF/6-31G* Levelu structures trans--trans trans-cis cis-cis
TABLE 7: Comparison of Structural and Energetic Parameters Calculated for Conformers of 2,ZDiphenylpropane at the HF/6-31G* Levelu
bond length (A) c1-02 C1-04 C1-03 03-C5 04-H6 bond angle (") 02-C1-03 02-C1-04 C5-03-C 1 H6-04-C1
torsion angle (") c 4 - c 1-c5-c9 C5-Cl-C4-C8 bond length (A) Cl-C2 Cl-C4 Cl-C5 bond angle (") c3-c2-c1 C4-Cl-C5 C4-Cl-C2 c 5 -c 1- c 2 C 1-C4-C6 Cl-C4-C8 c 1 -c5-c7 c 1 -c5-c9 out-of-planeangle (") C1 -C4/C4-C6-C8 c 1-c5/c5-c7-c9 conformational energies (kcaVmol) HF MP2
1.189 1.319 1.308 1.419 0.95 1 126.0 124.3 116.2 107.4
1.183 1.329 1.310 1.418 0.95 1 123.5 123.6 121.1 107.4
1.175 1.330 1.330 1.407 0.946 121.4 121.7 123.4 114.8
conformational energies (kcaymol) 0.0 0.0
HF MP2
4.01 3.55
14.21 13.07
The MP2 energies are calculated for the HF/6-31G* optimized structures. TABLE 5: Comparison of Structural and Energetic Parameters Calculated for Conformers of Dimethyl Carbonate at the HF/6-31G* Levela structures trans-trans trans-cis cis-cis bond length (A) c1-02 C1-03 C1-04 03-C5 04-C6 bond angle (") 02-C1-03 02-C1-04 C5-03-C1 C6-04-C1
1.190 1.313 1.313 1.417 1.417 125.4 125.4 116.4 116.4
1.184 1.322 1.315 1.419 1.416
1.179 1.328 1.328 1.409 1.409
122.7 124.7 116.3 121.2
19.0 19.0 31.4 31.4
conformational energies (kcal/mol) 0.0 0.0
HF MP2
3.80 3.25
24.53 23.58
The MP2 energies are calculated for the HF16-3 lG* optimized structures. TABLE 6: Comparison of Structural and Energetic Parameters Calculated for Conformers of Monophenyl Carbonate at the HF/6-31G* Level" structures bond length
c 1-02
I
I1
I11
IV
1.184 1.319 1.318 1.386 0.951
1.187 1.313 1.319 1.385 0.951
1.182 1.321 1.323 1.385 0.952
(A)
C1-03 C1-04 03-C5 04-H6 bond angle (") 02-C1-03 02-C1-04 C5 -03-C 1 H6- 04-C 1 C7-05-C3 C8-05-C3 torsion (") C7-C5-03-C(1) out-of-plane angle (") 03-C5/C5-C7-C8 nonbond distances
(A)
H(Ph). * -02 conformational energies (kcaUmo1) HF MP2
1.185 1.319 1.318 1.385 0.951 126.9 124.7 119.8 107.5 120.9 117.3
126.4 124.8 118.8 107.7 119.1 119.1
128.1 124.0 126.0 107.6 112.9 125.8
122.5 124.3 122.7 107.4 119.1 119.1
63.5
91.7
180.0
91.9
3.3
3.1
0.0
3.5
2.838
3.503
2.227
4.163
0.00 0.00
0.04 0.23
1.01 1.03
2.81 2.47
a The MP2 energies are calculated for the HF/6-31G* optimized structures.
electron of the carbon atom seems also to be shared by the other two oxygens to some extent. This is evident because the lengths of the C1-03 and C1-04 bonds are about 0.1 8, shorter than
structures
I
I1
111
IV
50.5 50.5
90.0 0.0
90.0 90.0
0.0 0.0
1.543 1.541 1.541
1.547 1.546 1.542
1.545 1.548 1.548
1.554 1.563 1.563
107.3 109.6 107.5 112.5 119.3 123.1 119.3 123.1
107.1 111.1 109.2 110.1 119.7 123.0 121.4 121.4
104.5 109.0 110.8 110.8 121.5 121.5 121.5 121.5
110.4 121.8 106.1 106.1 115.9 129.0 115.9 129.0
1.4 1.4
0.0 1.7
1.1
1.1
0.0 0.0
0.0 0.0
2.17 1.80
3.93 4.01
17.21 16.81
The MP2 energies are calculated for the HF/6-31G* optimized structures. the 03-C5 or 04-C6 bond lengths, irrespective of the method used. Our calculations incidate that the barrier heights for rotation about the C1-03 or C1-04 bonds (about 12 kcal/ mol) are significantly larger than those typically found for single bonds (about 3 kcaymol). When the C1-03 or C1-04 bond is elongated, the C1=02 bond is shortened, and vice versa. In terms of electronic structure, when the C1-03'or C1-04 bonds are elongated, either by twisting the bonds or by substituting the oxygen by different groups, electron density migrates from the C1-03 or C1-04 region to the C1=02 region, strengthening the n character of C1=02. By comparing the lengths of the C-0 bonds in Table 1 and across Tables 3, 4, and 5, one sees that the substitution of a hydrogen by a methyl group slightly decreases the bond lengths of C1-03 and C1-04, and increases the bond length of C1=02. But the major change of the internals is the bond angle C1-03-C(Me). For example, it changes from roughly 108" to 116' at the HF/6-31G* level, when the hydrogen is replaced by a methyl group. The trans-cis conformation is important because it is only a few kcaymol less stable than its corresponding trans-trans counterpart. At the HF/6-31G* level of theory, this energy difference is only 2.44 kcdmol for carbonic acid and 3.80 kcal/ mol for dimethyl carbonate. These energy differences decrease by some 0.2-0.5 kcal/mol to 2.19 and 3.36 kcal/mol, respectively, at the MP2/6-31G* level of theory. With the DFT(ACM)/DZVP method, these values are slightly lower still at 1.98 and 3.16 kcal/mol, respectively, close to the MP2 values. Katon and C ~ h e nreported ~~ the enthalpy difference of the trans-cis and trans-trans conformations in liquid dimethyl carbonate to be 2.6 f 0.5 kcal/mol, as determined from the temperature dependence of the infrared spectrum. The theoretical results determined in this study fall at the higher end of these experimental observations. The cis-cis conformers of the substituted carbonates are much less stable relative to those of the trans-trans. For carbonic acid, the cis-cis conformer is more than 14 kcal/mol higher in energy than that of the trans-trans conformer at the HF/6-3lG* level, indicating strong repulsion between the two hydrogens. When one of the hydrogens is replaced by a methyl group, evidently the steric crowding and electrostatics compensate in part because the carboxylic hydrogen is placed in such
5878 J. Phys. Chem., Vol. 99, No. 16, 1995
a way that it bisects the H-C-H angle of the opposite methyl group. However, if both hydrogens are replaced by methyl groups, the repulsion between these two methyl groups causes a much more significant increase in the total energy to 24.53 kcal/mol at the HF/6-31G* level. In carbonic acid the DFT results agree most closely with the highest level MP2/6-31G* results. For dimethyl carbonate, although all methods yield extremely high energies (> 20 kcal/mol), the D I T calculations lie slightly lower than the HF while the MP2 calculations indicate a slight rise. Monophenyl Carbonate. Laskowski et al.15316performed calculations on phenyl formate using several different basis sets from STO-3G to 6-31G* and concluded that, in order to generate the correct geometrical parameters for such a molecule, the split Gaussian basis set 6-3 1G with polarization functions is required. These authors calculated the phenyl ring rotation barriers in phenyl formate using the 6-31G* basis set. Also calculations were performed on the trans-cis and trans-trans conformations of phenyl carbonate. We extended these calculations by performing full geometry optimization on phenyl carbonate. The results are summarized in Table 6, which lists geometric parameters for the four optimized structures. A schematic illustration of the conformations calculated can be found in Figure 3. The minimum-energy structure I of monophenyl carbonate has Cl symmetry. The phenyl ring was found to be twisted out of the OCOO plane by 63.5" at the HF/6-31G* level, 59.1" at the MP2/6-31G* level, and 59.2" with the DFT(ACM)/DZVP method. Other geometrical parameters, such as the bond angle C1-03-C5, are also in good agreement with the results of these a ~ t h o r s . ' ~ ,The ' ~ energy barrier to rotation about the 03-C5 bond was found to be a little higher in monophenyl carbonate than in phenyl formate. Specifically, the HF/6-3 lG* barrier heights were found to be 0.04 kcal/mol for structure II and 1.01 kcaVmol for structure 111, compared with 0.01 kcal/mol and 0.72 kcal/mol, respectively. Focusing on the MP2 energies, we see that both barrier heights are slightly increased for I1 and I11 although the cisskewed structure is somewhat lower. Again, DFT also yields small energy differences between these states, somewhat lower than MP2. Note that these results do not agree with the semiempirical result^.^^.^^ The trans-cis conformer is 2.81 kcaV mol less stable than the trans-trans at the HF level and 2.49 kcal/mol at the MP2 level. These values agree well with previous ab initio result^'^ but are significantly higher than the value of 1.3 kcaVmo1 estimated by Williams and Flory'O from their empirical force field. Once again, we note significant conformational dependence of valence coordinates in Table 6. For example, both the C1-03-C5 and C8-C5-03 angles open significantly to relieve the steric strain in the planar configuration. In fact, the lower energy of this configuration in light of these deformations is somewhat unexpected and provides valuable information in determining both valence and nonbonded terms of analytical force fields. 2,ZDiphenylpropane. The four conformers calculated are shown in Figure 4. The minimum-energy structure (I) has C2 point group symmetry. At the HF/6-31G* level, both phenyl rings are twisted 50.5" from the C4-Cl-C5 plane in an opposite sense, such that the relative position of the two rings is roughly perpendicular, as illustrated in Figure 4. MP2 calculations were not performed for this molecule because of the computational limitations. DFT calculations yield very similar torsion angles (50.8"). Laskowski et al,'5J6*reported an average value of 5 1" for these dihedral angles based on their STO-3G and 4-31G SCF calculations. The value of these dihedral angles seems insensitive to the basis sets. Other
Sun et al. literature values for this dihedral angle include 60" from CNDO calculations, 50" computed using the MNDO and AM1 methods, and 48.1" from the PRDDO technique.I4 The crystal structure of the diphenyl carbonate of 2,2-bis(4-hydroxyphenyl)propane has been resolved,6s7 and it has been found that the intrinsic molecular C2 symmetry of the diphenylisopropane segment is broken. Consequently, two values have been reported for the corresponding torsion angles, 61" and 36". This is clearly due to the packing forces which make any rigorous quantitative comparison with calculations of isolated molecules meaningless. This distortion emphasizes the difficulty of comparing conformations of crystal structures with isolated molecules. However, we note that these distorted values are not far from those found in this and other theoretical studies. The calculated bond lengths and angles are also in the range of those obtained from X-ray crystal The C1-C2 bond length is longer than the Cl-C4 and Cl-C5 bonds, reflecting the fact that the latter two bonds are between sp3 and sp2 carbons, while the first bond is between two sp3 atoms. The C1 atom is the center of a distorted tetrahedral arrangement of substituent groups. The different bond angles C4-Cl-C2 (107.5") and C5-Cl-C2 (112.5") demonstrate that the plane of the C3-Cl-C2 atoms is not perpendicular to the plane of the C4-Cl-C5 atoms. This results from the steric crowding between methyl and phenyl groups. Note that the phenyl rings are also bent slightly for the same reason. The large difference in equivalent valence angles also begins to indicate the inherent problems in fixed, rigid geometry approximations and at the same time provides a test of analytical representations of the energy surfaces. Structure I1 was optimized with C, symmetry constraints. Consequently, the plane of the C2-Cl-C3 atoms is perpendicular to the plane of the C4-Cl-C5 atoms. One phenyl ring is twisted 90" from the C4-Cl-C5 plane while the other is in the plane. We have shown21that this structure, known as the Morino structure, is the most favorable transition state during the ring-flip process. The energy difference between this conformation and structure I provides a measure of the energy barrier to the ring flipping process. This barrier height was found to be 2.17 kcdmol at the HF/6-31G* level, which is slightly higher than the value of 1.9 kcaYmol that was reported by Laskowski et al.'5316 However, calculations at higher levels indicated that the energy barrier may be even lower. The MP2 energy calculated at the optimized 6-31G* geometry was found to be 1.80 kcal/mol. Full geometry optimization with the DFT(ACM)/DZVP method yields 1.79 kcal/mol. Note that the geometrical parameters for this structure are similar to those for structure I, which is consistent with the small energy differences between these two conformers. The results of several calculations of this energy barrier height have been published in the literature; the values are 5 kcaYmol from CND0,ls 10.1 kcaYmol from CNDO/2,I4 and 9 kcaVmol from PRDD0.I4 We note that these semiempirical and approximate ab initio results seem closer to the measured activation energy of 12 kcal/mo14 for glassy polycarbonate of bisphenol A (PCBPA) than those of ab initio SCF calculations. However, the problem of a direct comparison of calculations performed on an isolated molecule with experimental data derived from measurements on the condensed phase arises. It has been ~ h o w n ' that ~ . ~the ~ large difference found between theory and experiment for the higher level calculations is most likely due to strong intermolecular interactions, especially the steric crowding between phenyl rings in the condensed phase. Structure 111, with both rings twisted by 90" from the C4Cl-C5 plane, is another possible transition state for the phenyl ring flipping process, but this will result in an energy barrier
J. Phys. Chem., Vol. 99, No. 16,1995 5879
Ab Initio Calculations on Polycarbonates twice as high as for structure II as the transition state. A related phenomenon is that several intemals in this structure are distorted from the minimum-energy conformation. In particular, the bond angle C4-C1-C5 has decreased from 107.3' to 104.5", revealing that the two methyl groups are squeezed together by some 3'. Note that the bond lengths are slightly elongated too. Structure IV has the highest incremental energy over that of structure I. With both phenyl rings in the plane of the C4C1-C5 atoms, it has a total energy 16.81 kcal/mol higher than that of the minimum-energy structure (I) at the MP2/6-31G*/ /HF/6-31G* level, and 17.53 kcal/mol for the fully optimized DFT structure. This structure is associated with a highly distorted bond angle C4-Cl-C5 and bond lengths C1-C4 and Cl-C5. Obviously, the repulsion between the two phenyl groups is the dominant cause of this energy increase. These structures are important for obtaining an accurate description of the intramolecular interactions and consequently for developing a force field. Atomic Charges and Dipoles. Partial atomic charges are used as a standard model of the electrostatic interactions in computer simulations of macromolecular systems. It is clear that electrostatic interactions are important for polar systems, and the carbonate group OCOO certainly belongs to this category. The questions are, how does one derive reasonable charge representations, and how accurate are these? In this study, we compared three approaches: (1) Mulliken population analysis (MPA),4°,41which is perhaps the simplest and most widely used approach in the literature; (2) atomic partial charges determined by fitting the electrostatic potentials (ESP) calculated from point charges to those calculated from ab initio wave function^;^^-^^ and (3) the atomic polar tensor (AFT)48approach for planar molecules. In the ESP approach, the electrostatic potentials are represented by a grid around the molecule. The charges are then treated as parameters to fit the electrostatic potential calculated from ab initio wave functions. The resulting charges one obtains in this method often depend on the grid used.42 However, one can get around this problem to some extent if the grid points are placed far enough from the nuclei (an indication that higher multipoles may be significant resulting in a deviation from a l/r dependence over the entire range). After testing several options, we performed the calculations by using a spherical grid model in which the grid points were placed symmetrically on several sets of concentric spherical shells around each of the nuclei. The first shell was placed at a distance 1 8, greater than the van der Waals radii (1.7 8, for C, 1.4 A for 0, and 1.2 8, for H) from each of the nuclei in the molecule. Nine layers of shells separated by 0.25 8, were used in the calculation so that the thickness of the sampling region was 2.0 A. In the region where two spheres overlap, distances between all pairs of grid points were calculated. If any two grid points were closer than 0.25 A, one of the points was removed from the region to avoid oversampling. The ESP fits were also constrained by symmetry and dipole moments. The quality of the partial charge representation can be evaluated by calculating the root-meansquare (RMS) deviation, relative root-mean-square (RRMS) deviation and maximum absolute deviation (MaxD) between calculated from the atomic partial the electrostatic potential, charges, and V:, calculated directly from ab initio wave functions. RMS and RRMS deviations are given by the expressions
vc,
/ M
and M
M
where M is the number of grid sites. The APT method is applicable to planar molecules. Dinur and Haglefl8 showed that, for such a molecule, the atomic partial charges used to calculate electrostatic energies can be defined by the Cartesian out-of-plane derivatives of the molecular dipole moment. For a molecule in the X-Y plane, the partial charge of atom A is given as: (3) It has been shown48 that these charges, q A , reproduce the molecular dipole moment exactly for the planar system. These charges also reproduce most features of the electrostaticpotential energy surface for the molecule, because the molecular dipole moment dominates the electrostatic interactions for an electrically neutral molecule. These authors48went on to show that higher order atomic multipoles can be defined from analogous derivatives of the next higher molecular multipole moments. Despite the limitation of planarity, the advantages of this approach are clear: (1) there are no adjustable parameters, the charges do not depend on the way of setting grids, and other similar problems usually faced by numerical method; (2) the decomposition of electrostatic energies into contributions from atomic multipoles provides insight into the mechanism of electrostatic interactions. In particular, one can analyze the errors introduced by truncating higher order terms of the multipole expansion. Table 8 summarizes the results for carbonic acid and planar monophenyl carbonate (conformer 111) calculated with MPA, APT, and ESP methods. Table 9 contains data calculated for dimethyl carbonate that cannot be treated by this implementation of the ATP method because of the out-of-plane hydrogens of the methyl groups, and, consequently, only MPA and ESP calculations were performed. The first thing to note from Table 8 is the overall degree of fit. As one would expect, the electrostatic potential gives the best fit. This occurs because the charges are treated as adjustable parameters to fit exactly this observable-the electrostatic potential. The atomic polar tensor method produces a slightly less good fit. There are several noteworthy points to be made in comparing the ESP results and the APT results. First of all, we note that, in both cases, the dipole moment is fit exactly. Thus, this points out the fact that there are an infinite set of charges which will reproduce the dipole moment and result in very different electrostatic potentials in the surrounding space. Secondly, the deviations noted in the atomic polar tensor method indicate the failing of the representation of the continuous electronic distribution of the molecule by point charge^.^^,^^ This, of course, is well k n 0 w n , 4 ~ -and ~ HirschfeldS0and Coppens5' have also discussed this at length. One needs higher order moments such as atomic dipoles and perhaps quadrupoles to adequately represent such distributions. Systems such as those containing oxygen atoms with their lone-pair orbitals which are highly anisotropic would be expected to require higher order electric moments. We will retum to this shortly. The results also give
Sun et al.
5880 J. Phys. Chem., Vol. 99,No. 16, 1995
TABLE 8: Comparison of Calculated Atomic Partial Charges for Carbonic Acid and Monophenyl Carbonate STO-3G MPAa
ESPb
Carbonic Acid c1 0.470 0.977 02 -0.332 -0.545 03 -0.313 -0.561 H5 0.244 0.345 iu (D) 0.055 0.368 RMS(kcal/mol) 1.9 0.13 RRMS (%) 43.7 3.1 MaxD (kcal/mol) 5.8 0.5
6-31G* APTC
ESP
MPA
APT
(trans-trans) 0.651 1.025 0.937 0.707 -0.407 -0.594 -0.647 -0.528 -0.450 -0.706 -0.593 -0.487 0.328 0.491 0.448 0.397 0.368 0.876 0.176 0.176 0.9 2.5 0.3 0.8 22.6 48.2 5.6 15.6 3.2 7.6 1.1 2.8
Carbonic Acid (cis-cis) c1 0.483 1.040 0.660 1.054 0.984 0.721 02 -0.281 -0.499 -0.378 -0.531 -0.549 -0.505 03 -0.324 -0.629 -0.477 -0.727 -0.672 -0.508 H5 0.224 0.358 -0.336 0.466 0.455 0.400 P (D) 2.917 4.199 4.195 5.411 5.711 5.711 RMS (kcaymol) 3.1 0.1 1.1 0.7 0.3 1.0 RRMS (%) 33.4 1.6 12.0 5.7 2.0 7.7 MaxD (kcal/mol) 10.1 0.5 3.6 2.3 1.1 3.2
Phenyl Carbonate (Conformer 111) 0.482 0.998 0.645 1.000 -0.323 -0.533 -0.392 -0.580 -0.312 -0.572 -0.447 -0.678 -0.270 -0.494 -0.333 -0.644 0.143 0.444 0.245 0.366 -0.075 -0.226 -0.121 -0.239 C8 -0.094 -0.266 -0.132 -0.238 c9 -0.055 -0.022 -0.049 -0.195 c10 -0.053 0.001 -0.046 -0.197 c11 -0.071 -0.117 -0.087 -0.205 HC7 0.080 0.109 0.090 0.227 HC8 0.098 0.127 0.104 0.227 HC9 0.070 0.071 0.065 0.210 HClO 0.069 0.062 0.063 0.210 HCll 0.069 0.074 0.070 0.206 H6 0.243 0.344 0.324 0.479 P (D) 0.055 0.248 0.248 2.330 RMS (kcaymol) 1.9 0.16 1.0 4.6 RRMS (%) 58.2 4.9 30.9 94.0 MaxD(kcal/mol) 6.0 0.6 3.9 13.6
c1 02 04 03 c5 c7
1.013 -0.669 -0.644 -0.488 0.686 -0.475 -0.579 -0.004 0.034 -0.280 0.219 0.297 0.140 0.127 0.162 0.460 0.554 0.3 6.3 1.0
0.686 -0.492 -0.485 -0.286 0.190 -0.172 -0.162 -0.116 -0.121 -0.148 0.153 0.165 0.130 0.129 0.139 0.390 0.554 0.97 19.9 4.3
with our knowledge of the significantly anisotropic character of the oxygen atom.53 One of the functions of Tables 8 and 9 is to illustrate the basis set dependency, by comparing calculations with two basis sets: 6-31G* and STO-3G. It is clear that the MPA charges are strongly dependent on the basis set employed. Indeed, several values of those charges increase over 100% on changing from the STO-3G to the 6-31G* basis set. However, both ESP and APT charges appear less sensitive to the choice of the basis set. The dependence of the partial charges on the molecular configuration and conformation can also be examined from comparison of different molecules and different conformations. The data in Tables 8 and 9 show that both ESP and APT charges are slightly more sensitive to the substituent group and conformation changes. In contrast, the MPA charges appear less dependent on the configuration and conformation than the ESP and APT charges. Since the ESP charges provide a reasonable representation of the electrostatic potential, this dependence simply restates the fact that the electron distribution changes as the local chemical environment changes. Comparison of the charges listed in Table 8 indicates that, accompanied with the larger deviations, the size of APT charges are systematically smaller than those derived from the ESP approach. This is presumably due to the higher order moment effects which are rigorously excluded from the APT charges. As noted above, and shown by Dinur and Hagler,48higher order multipole terms in many systems are significant, and it follows that they will contribute significantly to the electrostatic potentials and intermolecular interactions. For planar molecules in the XY plane, the atomic point dipole is defined by:48
Mulliken population analysis. Fit to electrostatic potential. Calculated by using atomic polar tensor method.
TABLE 9: Atomic Partial Charges of Dimethyl Carbonate HF/STO-3G c1 02 03
c5
H5 1 H52 P (D)
RMS (kcaymol) RRMS (%) MaxD (kcal/mol)
HF/6-3 lG*
MPA
ESP
MPA
ESP
0.471 -0.330 -0.266 -0.056 0.088 0.082 0.195 2.0 39.8 5.9
1.217 -0.569 -0.477 -0.150 0.137 0.083 0.048 0.6 11.1 1.9
1.042 -0.599 -0.596 0.187 0.186 0.190 1.376 5.9 90.8 15.9
1.274 -0.694 -0.536 -0.03 1 0.073 0.132 0.389 0.6 8.7 1.7
us a feeling for the level of error in the various approximations. Thus, the electrostatic potential has RMS deviations of 0.3, but maximum deviations of as much as 1 kcaVmol in the case of the 6-31G* results. Since we are at distances greater than the van der Waals distance we would expect that such interactions would be sampled, and thus these errors represent errors that we might encounter in calculations. This suggests that, as we suspected from the above results, the higher order moments are indeed significant. In this case, unlike in other systems that have been treated, for example, hydrogen fluoride,52not only are the dipoles important, but higher order moments including quadrupoles seem to be important. Again, this is consistent
where the subscript A denotes atom A, MXZand Myz are the second moments relative to the center of mass, and g A is the ratio between the mass of atom A and the total mass of the molecule. The atomic point dipoles defined above exactly reproduce the molecular quadrupole moments. Consequently, by including the atomic point dipoles one can reproduce the electrostatic potential of the molecular charge distribution in the range where the contribution of molecular multipoles higher than the quadrupole moment is unimportant. We calculated the electrostatic potentials by including the atomic point dipoles for trans-trans and cis-cis conformers of carbonic acid; the results are summarized in Table 10. The size of the atomic dipoles are listed in the table. It is interesting to note that they are even more sensitive to the change of conformation. Clearly, the quality of reproducing the electrostatic potential is significantly improved if the atomic dipoles are included. The RRh4S deviations drop from 15.6% to 9.6% for trans-trans, and from 7.7% to 3.9% for cis-cis, showing improvements of approximately one-half in both cases. However, the influence of higher order moments are still detectable, particularly for the trans-trans conformer. But in general, by including atomic monopole and dipole moments into the calculations, one can recover the electrostatic potentials with an accuracy of 95% or higher.
Ab Initio Calculations on Polycarbonates
J. Phys. Chem., Vol. 99, No. 16, 1995 5881
(e)
TABLE 10: Atomic Monopole and Dipole (M) Moments of Carbonic Acid, Calculated Using the APT Method" trans-trans
c1 02 03 04
H5 H6
cis-cis
Q
IMl
Q
0.7070 -0.5277 -0.4867 -0.4867 0.3970 0.3970
0.2996 0.2697 0.5796 0.5796 0.4308 0.4308
0.7205 -0.5045 -0.5078 -0.5078 0.3998 0.3998
Q 0.8 15.6 2.8
Q+M 0.5 9.6 2.4
so (A) ESP APTQ APTQM
trans -trans
RMS (kcaVmol) RRMS (%) MaxD (kcaVmo1)
(MI 0.2332 0.1334 0.6252 0.6252 0.3407 0.3407
TABLE 11: Deviations of the Electrostatic Potentials Calculated from ESP, APTQ, APTQM, and MPA Multipoles with Respect to Those Calculated from HF/6-31G* Wave FunctioW
cis-cis
Q 1.o 7.7 3.2
Q+M 0.5 3.9 2.2
The deviations are calculated by comparing electrostatic potential values calculated from the atomic multipole moments and those calculated from the ab initio HF/6-3 1G* wavefunctions.
MPA
ESP
APTQ
The sampling grids for the electrostatic potentials reported in Tables 8, 9, and 10 are placed in such a way that the first layer of grid points is at 2.2 A, 2.4 A, and 2.7 8, from the centers on the hydrogen, oxygen, and carbon atoms, respectively. It is known that the point charges cannot provide a good representation of the electrostatic potential in regions too close to the n u ~ l e u s . ~ *We - ~ ~investigated the quality of reproducing the electrostatic potential at different distances from the nuclei, by using the atomic multipole parameters obtained. Three different regions around the molecules, with the first layer of the grids placed at SO= 2, 3, and 4 A from the nuclei, were investigated. The same grid setting as reported in Tables 8, 9, and 10 was used. The results, in terms of deviations of the calculated electrostatic potentials between those calculated from the atomic multipole expansions and those calculated from ab initio wave functions, are summarized in Table 11. Four sets of atomic multipole representations were used: point charges determined by the ESP method, point charges derived from the APT (APTQ), point charges plus point dipoles derived from the APT approach (APTQM), and Mulliken charges (MPA). Comparing the data of SO= 2 A with the data in Table 8 and 10, we see that all deviations are generally larger. In the best case, the ESP model, R R M S deviations are 9.03% for trans-trans and 2.97% for cis-cis; and the maximum deviations are 3.19 and 2.38 kcdmol, respectively. The quality of the fit increases quickly as the distance SOincreases, and all three deviations drop off significantly from SO= 2.0 A to 3.0 A. Finally, the energy surface becomes smooth and the deviations gradually converge. The large deviations at regions close to the nuclei (SO = 2 A) have a strong impact on short range intermolecular interactions, e.g., hydrogen bonding. However, this is not of concern for the carbonates studied here. Our results demonstrate in yet another way the importance and difficulty of representing highly anisotropic electron distributions in order to produce reasonably accurate electrostatic potentials and energies. Our study, along with other^$^-^^ points out the deficiencies of using atomic point charges and possible ways to improve the quality of representing the electrostatic field. Colonna et a1.46also recently reported their studies on modeling the electric field of formamide by using an overlap multipole expansion (OME) technique. However, in practice, the point-charge approximation is still the current standard because this model provides a simple, fast approach to calculate the electrostatic interactions, and the efficiency of energy computation is still a
APTQM
2.0 trans-trans RMS (kcallmol) 0.63 RRMS (%) 9.03 MaxD ( k c d m o l ) 3.19 RMS (kcaVmo1) 1.38 RRMS (%) 19.72 MaxD (kcdmol) 8.47 RMS (kcaVmol) 0.85 12.11 RRMS MaxD ( k c d m o l ) 4.49 RMS (kcal/mol) 3.08 RRMS (%) 44.03 MaxD (kcaVmol) 11.23
(a)
cis-cis RMS (kcdmol) RRMS (%) MaxD (kcaVmo1) RMS (kcal/mol) RRMS (a) MaxD (kcdmol) RMS (kcdmol)
RRMS (%) MPA
MaxD (kcaVmo1) RMS (kcal/mol) RRMS (%) MaxD (kcaUmol)
0.43 2.97 2.38 1.38 9.49 7.24 0.82 5.68 4.36 0.93 6.39 3.76
3.0
4.0
0.15 4.45 0.65 0.5 1 15.31 1.90 0.28 8.20 1.22 1.88 55.90 6.10
0.05 2.84 0.16 0.25 15.50 0.64 0.09 5.78 0.38 1.10 68.44 2.73
0.14 1.57 0.58 0.62 7.14 2.18 0.27 3.13 1.05
0.05
0.50 5.71 1.46
0.92 0.16 0.32 5.87 0.85 0.09 1.68 0.35 0.30 5.53 0.85
a The comparisons are for trans-trans and cis-cis carbonic acids, with the first grid-shell at different distances from the nucleus.
major concem in molecualr dynamics simulations. The importance of more accurate representations of electrostatic interactions and its effect on physical properties of interest is a subject of ongoing research.
IV. Summary Full geometry optimizations were performed on polycarbonate model molecules using ab initio Hartree-Fock, MP2, and DFT methods. The minimum energy, local minimum energy, and transition states of intemal rotations about covalent bonds in the backbone of polycarbonate were calculated, and these calculations reveal, on a higher level of theory than those reported previously, important information about the conformational structures and energies of such polymers. Our results show that the minimum-energy barrier height for synchronized phenyl ring rotations in 2,2-diphenylpropane is less than 2 kcal/mol. Complete rotation of a phenylene ring relative to the carbonate group costs only about 1.0 kcal/mol. Large fluctuations of the dihedral angle between the phenylene ring and the carbonate group, from -60" to 60", requires only a fraction of a kcdmol. The trans-cis conformational energies of carbonate molecules depend on the substituent group; it is about 3 kcal/mol for dimethyl carbonate, and 2 kcdmol for phenyl carbonate. It was found that large variations occur in valence coordinates, especially bond angles, for the same molecule in different conformations. These provide a challenge for analytical potential functions which must reproduce these trends if they are to accurately represent the energy surfaces of these molecules. As with other studies, it was found that density functional calculations result in energies and geometries that are compatible with the higher order of post Hartree-Fock results. The results in the carbonates are consistent with the MP2 results, and, because of less demanding computational resources, larger
5882 J. Phys. Chem., Vol. 99, No. 16, 1995
molecules such as the 2,2-diphenylpropane molecule can be treated with full geometry optimization. Based on the ab initio calculations, atomic partial charges are calculated by using three methods: Mulliken analysis (MPA), fitting electrostatic potential (ESP), and the atomic polar tensor approach (APT). The comparison begins to indicate weaknesses in the isotropic partial atomic charge model. The comparison of APTQ potentials with the ab initio potentials results in reduction in the fit, and the AFT charges are, in general, about 20% smaller than the ESP charges. These differences can be accounted for in part by including the atomic dipole moments which are found to be significant (APTQM). Regarding the transferability and sensitivity of the charges calculated, the APTQ,APTQM, and ESP parameters were found to be less sensitive to the basis set employed, but more sensitive to the conformations and configurations than the MPA analysis. Acknowledgment. This work was supported in part by Biosym Technologies’ Polymer Project, which is funded by a consortium of companies and government laboratories, and a grant from the Advanced Technology Program of the National Institute of Standards and Technology. H.S. thanks Professor B. E. Eichinger and Professor R. 0. Watts for beneficial discussions and encouragement. References and Notes (1) A review: Bubeck, R. A.; Smith, P. B.; Bales, S. E. Order in the Amorphous “state” of Polymers; Keinath, S . E., Miller, R. C., Rieke, J. K., Eds.; Plenum: New York, 1987; pp 347-358. (2) Schaefer, J.; Stejskal, E. 0.; McKay, R. A.; Dixon, W. T. Macromolecules 1984, 17, 1479. (3) Spiess, H. W. In Advances in Polymer Science; Kausch, H. H., Zachmann, H. G., Eds.; Springer-Verlag: Berlin, 1984; Vol. 66. Schmidt, C.; Kuhn, K. J.; Spiess, H. W. Macromolecules 1985, 18, 71. (4) Roy, A. K.; Jones, A. A.; Inglefield, P. T. Macromolecules 1986, 19. 1356. ( 5 ) Henrichs, P. M.; Nicely, V. A. Macromolecules 1990, 23, 3193. (6) Perez, S.; Scaringe, R. P. Macromolecules 1987, 20, 68. (7) Henrichs, P. M.; Luss, H. R.; Scaringe. R. P. Macromolecules 1989, 22, 2731. (8) King, J. A,, Jr.; Codella, P. J. Macromolecules 1990, 23. 345. (9) Ngai, K. L.; Rendell, R. W.; Yee, F. F.; Macromolecules 1988. 21, 3396. Xiao, C.; Yee, A. Macromolecules 1992, 25, 6800. (10) Williams, A. D.; Flory, P. J. J . Polym. Sci., Polym. Phys. Ed. 1968, 6 , 1945. (11) Tonelli, A. E. Macromolecules 1972, 5, 558. (12) Bendler, J. T. Ann. N . Y. Acad. Sci. 1981, 371, 299. (13) Johansen, H.; Rettrup, S.; Jensen, B. Theor. Chim. Acta 1980, 55, 267. Additional references relating to calculations on small carbonate molecules can be found in: Quantum Chemistry Literature Data Base (QCLDB); Hosoya, H., Morokuma, K., Yamabe, S., Ohno, K., Eds.; Release 92. (14) Bicerano, J.; Clark, H. A. Macromolecules 1988,21,585. Bicerano, J.; Clark, H. A. Macromolecules 1988, 21, 597. (15) Laskowski, B. C.; Yoon, D. Y.; McLean, D.; Jaffe, R. L. Macromolecules 1988, 21, 1629. (16) Jaffe, R. L.; Yoon, D. Y.; McLean, A. D. In Computer Simulation of Polymers; Roe, R. J., Ed.; Prentice-Hall: Englewood Cliffs, NJ, 1991.
Sun et al. (17) Sung, Y. J.; Chen, C. L.; Su, A. C. Macromolecules 1990.23, 1941. (18) Labrenz, D.; Schrikr, W. J . Mol. Struct. 1991, 249, 327. (19) Hutnik, M.; Argon, A. S.; Suter, U. W. Macromolecules 1991, 24, 5956. (20) Dauber-Osguthorpe, P.; Roberts, V. A,; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Proteins: Struct., Funct., Genet. 1988, 4 , 31. (21) Nilsson, L.; Karplus, M. J . Comput. Chem. 1986, 7, 591. (22) Weiner, S. J.; Kollman, P. A.; Ngnyen, D. T.; Case, D. A. J . Comput. Chem. 1986, 7, 230. (23) Allinger, N. L.; Rahman, M.; Lii, J. J . Am. Chem. SOC.1990, 112, 8293. (24) Maple, J. R.; Dinur, U.; Hagler, A. T. Proc. Natl. Acad. Sci. U.S.A. 1988, 85, 5350. (25) Maple, J. A.; Hwang, M. J.; Stockfisch, T. P.; Dinur, U.;Waldman, M.; Ewig, C. S.; Hagler, A. T. J . Comput. Chem. 1994, 15, 162. (26) Hwang, M. J.; Stockfisch, T. P.; Hagler, A. T. J . Chem. SOC. 1994, 116, 2515. (27) Sun, H.; Mumby, S. J.; Maple, J. R.; Hagler, A. T. J. Am. Chem. Soc. 1994, 116, 2978. (28) Sun, H. J . Comput. Chem. 1994, 15, 752. (29) Katon, J. E.; Cohen, M. D. Can. J . Chem. 1975, 53, 1378. (30) Collingwood, B.; Lee, H.; Wilmshurst, J. K. Aust. J . Chem. 1966, 19, 1637. (31) Mijlhoff, F. C. J . Mol. Struct. 1977, 36, 334. (32) Haeser, M.; Ahlrichs, R. J . Computer. Chem. 1989, 10, 104. Ahlrichs, R.; Baer, M.; Haeser, M.; Horn, H.; Koelmel, C. Chem. Phys. Lett. 1989, 162, 165. (33) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J . Chem. Phys. 1972, 56, 2257. Hariharam, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213. The 6-31G* basis set we used, which has been implemented in TURBOMOLE program, is slightly different from the conventional one. Instead of six, five d polarization functions (d3++, d;-+, d-, dyL.d a ) are employed. (34) Pople, J. A.; Seeger, R.; Krishnan, R. Int. Quantum Chem. 1977, IIS, 149. Krishnan, R.; Pople, J . A. Int. Quantum Chem. 1978, 14, 91. (35) Krishnan, P. C.; Pople, J. A. J . Chem. Phys. 1980, 72, 4244. (36) Parr, R. G.; Yang, W. Oxford University Press, 1988. Labanowski, J. A.; Andzelm, J. W. Springer-Verlag, 1991; Ziegler, T. Chem. Rev. 1991, 91, 651. (37) Becke, A. D. J . Chem. Phys. 1993, 98, 5648. (38) Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Can. J . Chem. 1992, 70, 560. (39) Ewig, S.; Liang, C.; Polavarapu, P. L.; Stouch, T. R.; Hagler. A. T. J . Am. Chem. SOC., in press. (40) Mulliken, R. S. J . Chem. Phys. 1955, 23, 1833. (41) Davidson, E. R. J . Chem. Phys. 1967, 46, 3320. (42) Williams, D. E.; Yan, J.-M. Advances in Atomic and Molecular Physics; Academic Press: New York, 1988; Vol. 23, p 87. (43) Williams, D. E. J . Comput. Chem. 1988, 9, 745. (44) Williams, D. ,E. Biopolymers 1990, 29, 1367. (45) Colonna, F.; Angyh, J. G.; Tapia, 0. Chem. Phys. Lett. 1990,172, 55. (46) Colonna, F.; Evleth. E.; Angyin, J. G. J . Comput. Chem. 1992, 13, 1234. (47) Urban, H. J.; Famini, G. R. J . Comput. Chem. 1993, 14, 353. (48) Dinur, U.; Hagler, A. T. J . Chem. Phys. 1989, 91, 2949. (49) Hagler, A. T.; Lapicirella, A. Peptides: Chemistry, Structure and Biology; Walter, R., Meienhofer, J. Eds.; Ann Arbor Science Publishers: Ann Arbor, MI, 1975. (50) Hirschfeld, F. L. Theor. Chim. Acta 1977, 44, 129. (51) Coppens, P.; Sevens, E. D. Adv. Quantum Chem. 1977, 10, 1. (52) Maple, J. R. and Hagler, A. T. Paper in preparation. (53) Hagler, A. T.; Lapicirella, A. Biopolymers 1976, 15, 1167. JW40082B
