J. Phys. Chem. 1993,97, 12752-12759
12752
A Force Field for Simulations of 1,2-Dimethoxyethane and Poly(oxyethy1ene) Based upon ab Initio Electronic Structure Calculations on Model Molecules Grant D. Smith Eloret Institute, I 1 78 Maraschino Drive, Sunnyvale, California 94087
Richard L. Jaffe NASA Ames Research Center, Moffett Field, California 94035
Do Y. Yoon' IBM Almaden Research Center, 650 Harry Road, San Jose, California 95120-6099 Received: March 12, 1993'
The force field parameters for simulations of 1,2-dimethoxyethane (DME) and poly(oxyethy1ene)were determined mainly utilizing the geometries and energies of the conformational minima and low-lying barriers in DME and diethyl ether (DEE) as determined from a b initio electronic structure calculations. A set of partial atomic charges, used to describe electrostatic interactions, was parameterized based upon the dipole moment and the partial atomic charges for DEE as determined from a b initio electronic structure calculations. The force field parameters thus determined reproduce well the experimental second virial coefficients of gas-phase DEE vs temperature. Furthermore, a stochastic dynamics simulation of DME based on our force field leads to an electron diffraction radial distribution function in good agreement with experiment. Detailed comparison of the electron diffraction radial distribution function as a function of conformer distribution with experiment corroborates the previous finding of Astrup that the order of conformer populations in DME is tg*gr > tgt > ttt.
Introduction Simulation of polymers via molecular mechanics, Monte Carlo, molecular dynamics and related techniques requires an accurate description of the conformational energy of the polymers of interest. In the case of poly(alky1 ethers) such as poly(oxyethylene) (POE), where strong polar interactions plus the oxygen gauche effect' must be carefully described, the use of general conformational force fields,parameterizedso as to describe reasonably a large number of molecules, is problematic. An example of the uncertainties incurred when using general conformational force fields is seen by comparing the conformational energies of the tgt and ttt conformers of l ,2-dimethoxyethane (DME), a model molecule for POE, as determined from various frequently used force fields. These conformers are illustrated in Figure 1. The popular MM2 force field yields an energy difference (the energy of the tgt conformer relative to the ttt conformer) of 0.70 kcal/mol* for a dielectric constant of 1.O. When the recommended dielectric constant of 1.5 is utilized, the MM2 force field yields an energy difference of 0.31 kcal/moP, while the MM3 force field yields an energy difference of 0.05 kcal/mol3. The Quanta force field4 yields an energy difference of 1.55 kcal/mol for a dielectric constant of unity. From the wide range of energies predicted by the various force fields for the relative conformer energies of DME, it is apparent that in order to be assured of achieving a quantitative description of properties for polymers containingthe types of interactionspresent in DME, it is best to parameterize force fields based upon model molecules representing the specific polymers to be studied. Prediction of equilibrium polymer properties requires accurate geometries and energies for the low-energy conformers of the model molecules, while prediction of dynamic properties addi Abstract published in Advance ACS Absrracts, October 15, 1993.
0022-3654/93/2097- 12752$04.00/0
tionally requires accurate descriptions of barriers to transitions between these conformers (rotational energy barriers). In a separate paper*we present results of an extensive ab initio electronic structure analysis of the conformational energies of DME. We began our analysis by obtaining the optimized geometries of the ten unique low-energy conformers of DME at the self-consistent-field (SCF) level using a D95** basis set, a Dunning6 double-{ basis set with polarization functions of the form (9s5pld/4s 1p)/ [4s2pld/2s 1p] . Theseoptimizedgeometries were then used in MP2 (second-order Maller-Plesset perturbation treatment of electron correlation) calculations of conformer energies using the larger D95+(2df,p) basis set, the Dunning double-fbasis set augmented with diffuse functionsand additional polarization functions for the carbon and oxygen atoms denoted (lOs6p2dlf/4slp)/[5s3p2dlf/2slp]. With this larger basis set, our calculations yielded an energy for the gauche conformation of the DME C-C bond of +0.1 kcal/mol relative to the trans conformation, much greater than the value of -0.4 kcal/mol predicted previously from a second-order rotational isomeric state model (RIS) parametrized to reproduce NMR vicinal coupling constants in gas-phase DME.' A salient result of our investigation was the determination that the energy of the tg*gi conformer (see Figure 1) lies only ca. 0.2 kcal/mol above that of the ttt conformation, due to strong 0.-H attractions. The conformer populations of DME obtained from the ab initio conformational energies were found to be in good agreement with those derived from electron diffractionexperiments and were used to successfully reproduce experimental values of dipole moments and NMR vicinal coupling constants in DME.5 Additionally, we find that a third-order RIS model for POE based upon the ab initio DME conformer energies*reproduces well the experimental values of the characteristic ratio and mean-square dipole moment ratio of POE. In this paper, using as a starting point an empirical force field for model alkyl ethers which has been applied successfully in modeling of crystalline poly(oxymethy1ene) (POM),Qwe detail 0 1993 American Chemical Society
Force Field for 192-Dimethoxyethaneand Poly(oxyethy1ene)
0
The Journal of Physical Chemistry, Vol. 97, No. 49, 1993 12753
8
TABLE I: Force Field Parameters Nonbonded Dispersion and Repulsion Parameters"
6
interacting pair (ij)
Aii (kcal/mol)
c.-c
14976.0 75845.0 2650.0 4320.0 33702.0 14176.0
o--o
m
H-*H C-H
0-.c 0.-H
Cii (A6 kcal/mol) 640.8 398.9 27.4 138.2 505.6 104.5
Bii (A-l) 3.090 4.063 3.740 3.415 3.577 3.902
Partial Atomic Chargesb atom (i) methyl carbon methylene carbon methyl hydrogen methylene hydrogen oxygen
charge (Qi) -0.163 -0,066 0.097 0.097 -0.256
Bonded Parameters: Stretches" stretch (ij) C-C C-H O-C
Figure 1. The lowest energy conformers of 1,2-dimethoxyethane. The sequence t$l&t$3 refers to the conformations of the O-C,C-C, and C-O bonds, with t = trans and g = gauche. A five-center P interaction is
Results lad Discussion Force Field Parameterization. In our force field we express the conformational energy E(r) of a molecule or ensemble of molecules as a function of nonbonded interatomic separations plus internal coordinates (bond lengths, valence angles, and dihedral angles) as
where r is the vector of atomic coordinates, PB(q) is the nonbonded energy associated with the atom pair i and j, b(q) is the potential energy due to stretching of the covalent bond between atoms i and j, @(OijL) is the potential energy due to bending of the valence angle defined by atoms i, j, and k, and E(@$u) is the potential energy due to twisting of the dihedral angle defined by atoms i, j, k, and 1. The summations run over all interactions of each type present in the molecule or ensemble of molecules. Coupling between variousvalence stretches, bends, and torsions is introduced only through the nonbonded terms. As a result, this force field may not reproducethe detailedvibrational
Pii (kcal/mol/A*)
1.51 1.09 1.39
618.0 655.0 739.0
Bonded Parameters: Bends
indicated for the tg*f conformer (see Appendix).
the parameterization of a force field specific to DME and POE based upon the results of our ab initio electronic structure calculations on model molecules. In addition to the ab initio values of the DME conformer energiesand geometriesdescribed above, we utilize ab initio energies and geometries of rotational barriers between the lowest energy conformations. These were calculated in a manner analogous to that for the low-energy conformers:s geometries of the saddle points were determined at the SCF level with a D95** basis set and the energiesat the MP2 level with the larger D95+(2df,p) basis set. In parameterizing a set of partial atomic charges for use in the force field, dipole moments, and partial atomic charges for diethyl ether (DEE) were used. These were determined from ab initio electronic structure calculationson DEE analogous to those described above for DME. Finally, we test the force field by performing a stochastic dynamics (SD) simulation of gas-phase DME for the purpose of comparing conformer populations with estimates from electron diffraction experiments. In addition, the electron diffraction radial distribution function, determined from the simulation, is compared with experiment.
roii(A)
bend (iik)
O'iik
C-C-C" C-C-Ha H-C-H' C-C-Ob O-C-H" C-O-CO
(radians)
kbiik
1.9373 1.9109 1.8902 1.9031 1.9211 1.9471
~
(kcal/mol) 107.0 86.0 77.0 172.0 112.0 149.0
Bonded Parameters: Torsions torsion (ijkl)
fold
O-C-C-H' C-C-C-H" H-C-C-H" C-O-C-Ha
3
O-C-C-Ob O-C-C-Ob O-C-C-ob C-O-C-Cb C-O-C-Cb C-oc-Cb a
$"iju (radians)
0 0 0 0 r *I2 0 0 712 0
k * ~(kcal/mol) u 0.28 0.28 0.28 0.81 0.05 2.55 0.00 1 .OO 0.70 0.32
Taken from ref 9. From this work.
properties of polyether chains. The nonbonded terms include intramolecular and intermolecular contributions while the stretching, bending and torsional terms are strictly intramolecular. The energy functions and their parameterization,based primarily upon the conformational properties of DME and DEE as determined by ab initio electronic structure calculations, are described below for each type of interaction (nonbonded, stretches, bends, and torsions). Nonbonded Parameters. In our force field we describe interactions between nonbonded atoms i and j, in terms of dispersion and repulsion interactions, and coulombic interactions, given as
~
~
+ ~= ( exp(-rijBi) r ~ ~ - cijrij4 )
332.08q,q
(2)
r.. U
where rij is the separation between atoms i and j (in A), and qi and qj are the partial atomic charges of atoms i and j, respectively. Energies are in units of kcal/mol. A pair of atoms i and j are considered nonbonded if they are in different molecules or are atoms in the same molecule separated by two or more atoms. The parameters for the nonbonded dispersion and repulsion interactions, Aij, Bij, and Cij, for our force field were adopted from an empirical force field successfullyused to describe crystal structures and energetics of a closely related polymer to POE, poly-
Smith et al.
12754 The Journal of Physical Chemistry, Vol. 97, No. 49, 1993 0.122 0.122
Q?
0.122
CkB
0.122
W
0.W
0.097
tr Figure 2. The tt conformer of diethyl ether. Partial atomic charges are given for each atom. The charges were obtained by adjustment of the Mulliken charges as explained in the text. Partial charges for 1,2dimethoxycthane,given in Table I, were determined from these charges as explained in the text.
(oxymethylene) (POM).g For covenience these parameters are reproduced in Table I. It is a well-established practice to use the pairwise coulomb energy based on atom-centered partial electric charges to model the electrostatic contribution to the nonbonded energy in polar molecules. We successfully used such an approach recently in a force field for poly(viny1 chloride).lo In such an approach, the set of partial atomic charges used in the force field description of a polymer must result in both charge neutral chains and repeat units. For straightforward application of the force field in molecular mechanics and dynamics simulations,equivalent atoms (e.g., methylene hydrogen atoms) should have identical charges that are invariant with respect to conformational changes. We chose to use the results of a Mulliken population analysis of DEE as the starting point for obtaining these partial charges for our force field for POE. DEE can be viewed as a POE monomer, -[CH20CH2]-, with methyl end groups. In determining the partial charges, the optimized geometry and dipole moment of the tt conformation of DEE was first established via ab initio molecular orbital calculations at the SCF level using a D95** basis set. The tt conformation of DEE is illustrated in Figure 2. The resulting Mulliken atomic charges were adjusted as follows. First, the methyl hydrogen charges were replaced by the average value for all methyl hydrogen atoms. Second, the methyl carbon charge was reduced by 0.015 and the oxygen charge increased by 0.03 in order to obtain a charge neutral repeat unit, -[CH2OCH2]-, while maintaining overall charge neutrality for the molecule. Finally, all charges were scaled proportionally to yield agreement between the calculated dipole moment of tt DEE, PDEE~given by
where the sum of position vector multiplied by the partial charge is taken over all atoms of the molecule, and the ab initio SCF value of 1.312 Debye. (The experimentally determined dipole moment for DEE is 1-06ll-l. 15 Debye.'*) The resulting charges are illustrated in Figure 2. For force field computations involving DME and POE polymers of the form H[CH20CHJnH we treat the terminal methyl groups as methylene groups with an additional methylene hydrogen and set the methyl carbon charge equal to a methylene carbon charge minus the charge of the additional hydrogen in order to maintain charge neutrality. These charges are given in Table I. As a test of the nonbonded parameters, the gas-phase second viral coefficient as a function of temperature was calculated for DEE using the nonbonded parameters given in Table I and the partial charges as shown in Figure 2. A fixed tt molecular geometry as given by the ab initio calculations was used. For a fixed molecular geometry, the second virial coefficient is given by the expression13 B2(T)
J ."Jrf'sin
ec[exP(-v(rc,ec,~,,e,4,$)/kT)1] d$d&ied&dOcdrc
(4)
-2000
250
300
350
400
450
T (K) Figure 3. The gas-phase second virial coefficient of diethyl ether. The calculatedvalues (solid curve) are for a fixed tt conformationalgeometry. Experimental data (filled circles) are from ref 13.
where the polar coordinates r,, &, and 4c describe a vector connecting the centers of mass of two DEE molecules and the Eulerian angles 0,4, and describe the relative orientations of the molecules. Here, U(rc,ec,4c,e,4,+) is the intermolecular potential energy of the molecules, k is Boltzmann's constant, and T is temperature. The integration of eq 4 was performed numerically. The number of function evaluations for each of the six integrations was doubled until each integral converged to a relative tolerance of 0.05., i.e., II,, - Idd/&wl I0.05, where I,,, and [old are the new and old values of the integral, respectively. A comparison of calculated and experimental values13 for the second virial coefficient as a function of temperature is shown in Figure 3. Agreement is good, especially considering that for a given force field using a fixed geometry results in a reduction in magnitude of the second virial coefficient relative to that which would result if flexibility were taken into account. A further test of the nonbonded force field will be to compare calculated thermodynamic properties of POE from molecular dynamics simulations with experimental values. Bonded Parameters: Stretches and Bends. Bond stretching functions of the form
+
1 - rij)' P ( r i j ) = -Pij(roij 2
and bending functions of the form
where rij is the bond length and Bijk the valence angle, are used to describe two center and three center bonded interactions, respectively. The diagonal harmonic stretching and bending force constants ksij and kbijk were adopted from the alkyl ether force field? with a single exception. These force constants are similar to those obtained in our ab initio electronic structure calculations for DME and DEE with the exception of the C-C-O bend. Here the force constant from the alkyl ether force field is significantly smaller than the ab initio value. We believe that this coordinate is described better using the larger ab initio force constant (scaled by 0.8),especially given the relatively large force constants for the off-diagonal terms involving this coordinate in the full alkyl ether force field. The scaled ab initio force constant is given Table I. The valence geometry parameters of the force field were calibrated by comparing the predicted geometries of the two lowest energy conformers of DME, ttt and tgt, illustrated in Figure 1, with results of the ab initio calculations. The predicted geometries were determined by energy minimization (minimization of eq 1 with respect to the atomic positions), using the nonbonded, stretching and bending functions given by q s 2, 5 , and 6,respectively, and fixing the torsional angles at the ab initio values for the respective conformers. The geometric parameters
Force Field for 1,2-Dimethoxyethane and Poly(oxyethy1ene)
The Journal of Physical Chemistry, Vol. 97, No. 49, 1993 12755
TABLE II: DME and DEE Valence Geometries C-CO
case
C-O
C-H
C-O-Cb
C-C-O
O-C-H
C-C-He
H-C-He
107.9 108.8 109.7 108.6
110.9 110.7 110.5 110.6
109.6 109.4 109.1 109.5
107.9 108.8 108.1 108.1
108.9 107.8 108.6
1 10.0 110.6 110.2
110.3 110.0 110.2
107.3 108.1 107.5
DME rrr ab initio trr force field rgr ab initio rgr force field
1.52 1.52 1.51 1.52
1.40 1.40 1.40 1.40
1.09 1.09 1.09 1.09
rr ab initio rr force field rr microwaved
1.52 1.52 1.52
1.40 1.39 1.41
1.09 1.09 1.09
114.0 113.8 113.9 113.9 DEE 114.3 114.1 1 12.4
a Bond lengths are in A. b Valence angles are in deg and are averages if nonequivalent angles with a given label are present. c Angles involving methylene hydrogen atoms only. From microwave spectroscopy experiments of ref 11.
Pu and Boijk were adjusted so as to give the best agreement between force field and ab initio geometries of the ttt and tgt conformers of DME. The geometric parameters for stretches and bends so determined are shown in Table I. A comparison of force field and ab initio geometries for the ttt and tgt conformers of DME ispresentedinTable 11. Agreement isquitegood, with thegreatest discrepancy being in the C-C-O bond angle, where the force field fails to reproduce the 2 O difference in the angle between the ttt and tgt states predicted by the ab initio calculations. The ab initio and force field geometries for DEE, using the nonbonded and valence parameters given in Table I and the partial charges given in Figure 2, are also shown in Table 11. Agreement between the ab initio and force field geometries is good. These geometries also agree well with the experimental geometry of tt DEE determined from microwave spectroscopyll as indicated in the table. We therefore conclude that the use of diagonal force constants alone is sufficient for representing valence geometries of DME and DEE and will therefore be adequate for simulations of POE. However, we note that the neglect of off-diagonal force constants limits the ability of the force field to quantitatively reproduce vibrational spectra.9 Bonded Parameters: Torsions. In order to successfully model conformation dependent equilibrium and dynamic properties of any polymer, it is critical that the energies of the low-energy conformations and the rotational energy barriers between the low-energy conformations be accurately represented. The three lowest energy conformations of the model molecule DME are illustrated in Figure 1, and the low-energy transitions between these conformers are illustrated in Figure 4. Energies for these low-energy DME conformationsand associated rotational energy barriers, determined from ab initio molecular orbital calculations at the MP2 level using a D95+(2df,p) basis set and geometries optimized at the SCF level using a D95** basis set, are given in Table 111. At reasonable temperatures, these conformers constitute a large majority of the DME conformer population. Also shown in Table I11 is the ab initio energy for the ttg conformer, which is an important intermediate in one of the low-energypaths between the ttt and tg*gi conformers, as illustrated in Figure 4. In our force field, the intrinsic torsional energies are expressed by functions of the form n
where &ju is the dihedral angle formed by the atoms, i, j, k, and I of a four atom consecutive sequence. The torsional energy parameters k' and function minimum angles, 4oiju, were adopted f r o m x e alkyl ether force field? except for the 0-CC-O and C-O-C-C torsions. Parameters for these torsions, each represented by a sum of 1-, 2-, and 3-fold terms (n = 1, 2, and 3), were determined so as to optimize agreement between force field and ab initio energies for the DME conformations and rotational energy barriers given in Table 111. The torsional parameters so determined are given in Table I, and the resulting conformational energies and rotational energy barriers are listed
2.5 2.0 1.5 1.o 0.5
2.5 2.0 1.5 1.o
0.5
rotational isomerization path F w e 4. Schematic representations of the low-energy paths between the low-energy conformers of 1,2-dimethoxyethane. Illustrated are conformational energies and barriers from ab initio calculations, from force field calculations, and from calculations with the modified force field (see Appendix).
in Table I11and shown in Figure4 for the transition paths between low-energy conformations. Examination of Table I11 indicates that the force field does a good job in representing the energies of the low-energy DME conformers. For the rotational energy barriers, however, we found that better agreement between ab initio and force field values than about 0.4 kcal/mol for both the ttg i-tg*gi and t t t t g t barriers could not be obtained simply by adjusting the 1-, 2-, and 3-fold barriers for the 0-C-C-O and C-O-C-C torsions while maintaining a reasonable energy for the tg*gi conformer. An attempt to improve agreement for the rotational barriers within the framework of the simple functional form of the force field presented here, by introduction of an additional simple energy term, is discussed in the Appendix. ConformationalEnergies and Geometries. Ab initio electronic structure calculations indicate that the tgt and tg*gi conformer energies lie only slightly higher than the ttt energy (Table 111). Our force field calculations indicate that the nonbonded energies in these conformers are considerably higher than in ttt, as shown in Table 111, where the force field energies for low-energy conformers and barriers are broken down into contributions from nonbonded terms and bonded terms. These differences in nonbonded energies are due primarily to strong 0-0 coulombic repulsion, which is offset somewhat in the tg*gi conformer by strong 0-H coulombic attraction. In our force field these coulombic repulsion is offset by favorable bonded interactions,
Smith et al.
12756 The Journal of Physical Chemistry, Vol. 97,No. 49, 1993
TABLE IIk Energies of Low-Energy DME Conformers and Barriers energies" force field coulombic
stretch
bonded bend
torsion
0.00 2.27 0.99 0.22
0.00 0.02 0.07 0.04
0.00 0.01 0.52 0.34
0.00 -2.06 -1.60 0.46
0.62 0.02 0.21 1.85
0.04 0.01 0.07 0.02
0.02 0.1 1 0.23 0.13
0.43 1.67 1.16 -0.49
nonbonded
conformer
ab initio
total
dispersion/repulsion Conformers
ttt tgl t P S ttg
0.00 0.15 0.23 1.43
0.00 0.14 0.19 1.38
0.00 -0.12 0.22 0.31
ttt-tgt trt-ttg trgi-tBg' tPtJBS
2.31 2.04 2.03 1.36
1.88 1.97 2.48 1.41
0.76 0.16 0.82 -0.10
Barriers
0
Energies are in kcal/mol. All values are relative to the ttt conformer.
TABLE IV Energies of Higher Energy DME and DEE Conformers and Barriers conformation ab initio energy force field energy DME 1.74 3.19 1.70 2.07 2.64 2.66 7.84
1.51 1.64 1.86 2.41 3.08 3.13 8.90
DEE 1.45 (1.37 2.94 2.55 7.24
0.10)*
1.17 2.23 1.90 6.71
I, Energies are in kcal/mol, relative to the ttt conformer (DME) or tt conformer (DEE). Experimentally determined value from ref 14.
TABLE V Torsional Geometries of DME and DEE Conformers and Barriers ab initio force field 4 c o k 4 O c b k 4oc conformation (deg) (deg) (deg) (deg) (deg) (deg) Low-Energy DME Conformers and Barriers ttt tgt tgfg' ttg ttr-tgr ttr-ttg
tfSfPf tPtfPf
0.0 106.4 104.9 4.9 60.0 1.3 67.5 108.9
0.0 -4.7 -95.8 101.3 -1.4 62.6 -105.5 -53.5
0.1 -2.4 -1.5 -0.9 -1.2 -0.5 0.2 -2.4
-0.1 102.8 94.9 4.8 56.0 -0.1 46.0 103.7
0.1 -2.4 -100.1 95.0 -1.2 59.0 -95.0
106.6 115.6 102.8 -105.1 0.9 12.1
94.0 97.7 -98.6 73.5 -93.8 96.6
-60.0
Other DME Conformers t P P
PBP Pkgi gfSgf
primarily a strong 2-fold function for the 0-C-C-O torsion (see Table I) which favors the gauche state. The need to include such a strong function is indicativeof the strength of the oxygen gauche effect in DME. Ab initio and force field energies for the six remaining higher energy conformers of DME, which were not considered in parameterizing the force field, are presented in Table IV. As the populationofeachofthesehigherenergyconformers islow (
