TraPPE-UA Force Field for Acrylates and Monte Carlo Simulations for

J. Phys. Chem. B 2009, 113, 6415–6425

6415

TraPPE-UA Force Field for Acrylates and Monte Carlo Simulations for Their Mixtures with Alkanes and Alcohols Katie A. Maerzke,† Nathan E. Schultz,‡ Richard B. Ross,‡ and J. Ilja Siepmann†,* Departments of Chemistry and of Chemical Engineering and Materials Science, UniVersity of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455, and Corporate Research Materials Laboratory, 201-2E-23, 3M Company, St. Paul, Minnesota 55144 ReceiVed: December 1, 2008; ReVised Manuscript ReceiVed: February 15, 2009

An extension of the transferable potentials for phase equilibria-united atom (TraPPE-UA) force field to acrylate and methacrylate monomers is presented. New parameters were fit to the liquid density, normal boiling point, saturated vapor pressure, and (where experimentally available) critical constants of 1,3-butadiene, isoprene, methyl acrylate, and methyl methacrylate using Gibbs ensemble Monte Carlo simulations. Excellent agreement with experiment was obtained for the parametrization compounds and seven additional acrylate and methacrylate compounds, with average errors in liquid density and normal boiling point of approximately 1%. The TraPPE-UA force field also predicts accurate heats of vaporization at 298 K. In addition, Gibbs ensemble Monte Carlo simulations of binary vapor-liquid equilibria for the mixtures methyl acrylate/1butanol and methyl acrylate/n-decane show that the TraPPE-UA acrylate force field performs well in mixtures with both polar and nonpolar molecules. These simulations also indicate structural microheterogeneity in the liquid phase of these mixtures. 1. Introduction Acrylate monomers are the building blocks for many important polymers which are used in numerous industrial applications, such as films and adhesives. The properties of acrylatebased polymers can often be tailored through blending with different polymers and with addititives, such as plasticizers. Despite the widespread use of these polymer blends, little progress has been made in developing a transferable force field to predict relative miscibilities of monomers, oligomers, polymers, and additives that are candidates for a blend. Previous simulation studies of acrylates have focused on the solubility of poly(methyl acrylate) in supercritical CO21 using the MM3 force field,2 the arrangement of molecules within acrylate and methacrylate clusters3,4 using the MM2 force field,5,6 methods for studying amorphous polymers7 using the AMBER force field,8,9 and ordering in liquid crystals10 using the Dreiding force field.11 However, none of these force fields was parametrized specifically for the thermophysical properties of acrylates or methacrylates; hence, their accuracy in these applications is uncertain. Therefore, in order to accurately model these systems a transferable force field which has been validated for acrylates is necessary. To this end, we have extended the transferable potentials for phase equilibria-united atom (TraPPE-UA) force field to include parameters for acrylate and methacrylate monomers. Beginning with a generalized force field for linear alkanes,12 it has been an ongoing effort of this research group to develop computationally efficient transferable force fields which accurately predict a wide range of thermophysical properties. These force fields, known collectively as the transferable potentials for phase equilibria, or TraPPE, force field, have been developed for linear and branched alkanes,12-14 alkenes,15 * Corresponding author. E-mail: [email protected]. † University of Minnesota. ‡ 3M Company.

alcohols,16 ethers,17 aldehydes,17 ketones,17 glycols,17 carboxylic acids,18 amines,19 nitroalkanes,19 nitriles,19 amides,19 thiols,20 sulfides,20 disulfides,20 carboxylate esters,21 benzene,22,23 hetereocyclic aromatics,23 substituted aromatics,24 and fused-ring aromatics,24 as well as small molecules such as H2O,25 CO2,26 N2,26 O2,27 and NH3.28 For organic compounds, the TraPPE force field employs two levels of sophistication. The simplest and most efficient, known as TraPPE-united atom (UA), uses a united atom representation for CHx groups. TraPPE-explicit hydrogen (EH)14,19,23 is more accurate, but also more computationally expensive. Hence, for the acrylate and methacrylate force fields, which will be used in the future to simulate large systems involving polymers, we have chosen to use the UA representation. Here we should mention two other transferable force fields, namely the OPPE (optimized potentials for phase equilibria)29-31 and SPEADMD (step potential equilibria and discontinuous molecular dynamics)32-34 force fields, that are available for a wide range of organic molecules and are also parametrized to vapor-liquid coexistence curves. The OPPE force field uses interaction sites that are displaced from the nuclear sites for many units (e.g., CHx units), thereby introducing additional parameters and computational complexity, but achieving somewhat more accurate predictions of the saturated vapor pressures than the TraPPE-UA version. The SPEADMD force field uses discontinuous step potentials and places greater weight in the parametrization on saturated vapor pressures than liquid state properties. The general philosophy behind the development of the TraPPE force field is not only to predict a wide variety of thermophysical properties as accurately as possible, but also to minimize the number of different atom types in order to ensure the maximum amount of transferability. Transferability ensures that the simulation of larger molecules with multiple functional groups is possible with a straightforward application of existing force field parameters; that is, fitting new parameters is not necessary. Furthermore, parameters are in general fit to smaller

10.1021/jp810558v CCC: $40.75  2009 American Chemical Society Published on Web 04/09/2009

6416

J. Phys. Chem. B, Vol. 113, No. 18, 2009

Maerzke et al.

molecules and then can be transferred to larger molecules; i.e., the CH3 pseudoatom parameters are fit to ethane, then transferred to larger molecules such as n-decane. The TraPPE force field has been used recently to predict viscosities and vapor-liquid equilibria for polyhydric alcohols35 as well as Hildebrand solubility parameters for various organic molecules.36 In addition, the TraPPE force field has been used extensively in simulations of complex chemical systems, such as octanol-water partitioning,37 gas chromatography,38-40 reversed-phase liquid chromatography,41-43 interfacial properties of aqueous systems,44,45 adsorption in zeolites and metal-organic frameworks,46-48 and various polymer systems.49-52 In the present work, the TraPPE-UA force field is extended to acrylate and methacrylate monomers. The remainder of this paper is organized as follows: the details of the force field and parametrization are discussed in section 2. The simulation details are provided in section 3. In section 4.1 the results for singlecomponent systems are presented, followed by binary systems in section 4.2. The conclusions of this work can be found in section 5. 2. Model In order to increase computational efficiency, the TraPPEUA force field employs pseudoatoms for all CHx groups (where 0 e x e 4) which are located at the position of the carbon atom. The nonbonded interactions are described by the simple pairwise-additive Lennard-Jones (LJ) and Coulomb potentials:

unonbonded(rij) ) 4ij

[( ) ( ) ] σij rij

12

σij rij

6

qiqj + 4π0rij

(1)

where rij, ij, σij, qi, and qj are the separation, LJ potential well depth, LJ diameter, and partial charges, respectively. The unlike LJ interactions are determined using the Lorentz-Berthelot combining rules:53

1 σij ) (σii + σjj) 2

and

ij ) (iijj)1/2

(2)

Nonbonded interations are calculated for intermolecular interactions and intramolecular interactions of atoms separated by four or more bonds. Intramolecular Coulombic interactions of beads separarated by three bonds are also calculated, but reduced in magnitude by a scaling factor of 1/2. Additionally, for intramolecular interactions between a hydroxyl hydrogen (which is not protected by a LJ interaction site) and an oxygen separated by four bonds we include a short-range repulsive term, given by:

urepulsive )

ax rij12

(3)

where for an ether oxygen aether/kB is 4.0 × 107 K Å12.17 Rigid bond lengths are used for bonded interactions. Bond angles are treated using a harmonic potential:

ubend(θ) )

kθ (θ - θeq)2 2

(4)

where θ is a bond angle, θeq is the equilibrium value for that bond angle and kθ is the force constant. For beads separated by

Figure 1. Schematic drawing of the four molecules used in the parametrization process. The squares and diamonds indicate sites for which this work introduces new Lennard-Jones parameters and partial charges, respectively.

TABLE 1: TraPPE-UA Force Field Parameters for Bond Lengths bond type

r (Å)

bond type

r (Å)

CHx-CHy12 CdC15 CdO21

1.54 1.33 1.20

C(carbonyl)-O21 C(carbonyl)-C(sp2)17 CHx-O(ether)17

1.344 1.52 1.41

three bonds, torsional interactions are computed using a cosine series:

utor ) a0 + a1 cos(φ) + a2 cos(2φ) + a3 cos(3φ) + a4 cos(4φ) (5) For the acrylate and methacrylate models, the transferability of the TraPPE force field was exploited to produce most of the necessary parameters. Bond lengths and bending parameters were taken from existing TraPPE force fields (see Tables 1 and 2). Lennard-Jones parameters for sp3 carbon pseudoatoms are taken from the TraPPE-UA alkane force field12 and the oxygen and carbonyl carbon parameters from the ether and ketone force fields,17 and CH2(sp2) is taken from the alkene force field.15 As the CHx pseudoatoms represent not only the repulsive and dispersive interactions of the valence electrons for the x C-H bonds, but also a portion of the interactions of valence electrons belonging to the C-C bonds, the LJ parameters need to reflect both contributions (see refs 12 and 13 for a more detailed description). Hence, new parameters were fit for CHx(sp2) pseudoatoms with single bonds to another CHx(sp2) pseudoatom (i.e., a conjugated double bond), a bonding environment not previously parametrized for the TraPPE force field. Following the philosophy of TraPPE, missing LJ parameters for CHx(sp2) pseudoatoms were fit to the simpler molecules 1,3-butadiene (CH) and isoprene (C) as well as methyl acrylate and methyl methacrylate, respectively. Figure 1 graphically shows the sites included in the parametrization procedure. Parameters were chosen to simultaneously reproduce the liquid density, normal boiling point, and critical temperature and density (where available). Partial charges were computed using the CM4 charge model of Kelly et al.54 with 1-octanol as the universal continuum solvent for the TraPPE charge parametrization.23 1-Octanol was chosen because of its intermediate dielectric constant and combination of polar and nonpolar moieties. To determine the partial charges, the gas-phase methyl acrylate structure was first optimized using Kohn-Sham density functional theory at the B3PW91/6-31++G(d,p) level of theory in the Gaussian 03

TraPPE-UA Force Field for Acrylates

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6417

TABLE 2: TraPPE-UA Force Field Parameters for Bond Angles angle type

θ0 (deg)

CHx-CH2-CHy CsCHdCH215 C-CdCH2 CHx-C(sp2)-CHy15 12

114 119.7 119.7 119.7

kθ/kB (K) 62500 70420 70420 70420

TABLE 3: TraPPE-UA Force Field Parameters for Nonbonded Interactions site 3 12

CH3(sp ) CH2(sp3)12 CH3(ether)12 CH2(ether)12 CH2(sp2)15 O(ether)17 C(carbonyl)17 O(carbonyl)17 CH(sp2) C(sp2)

σ (Å)

/kB (K)

q (e)

3.75 3.95 3.75 3.95 3.675 2.80 3.82 3.05 3.71 3.85

98 46 98 46 85 55 40 79 52 22

+0.25 +0.25 -0.25 +0.40 -0.40 -

program.55 The optimized gas-phase geometry was then used to optimize the solute’s charge distribution in the continuum solvent calculations at the same level of theory using the Minnesota Gaussian solvation model (MN-GSM) implemented in Gaussian 03.56 Partial charges for CHx pseudoatoms were determined by combining the CM4 charges on the carbon and hydrogen atoms. Nonbonded parameters for the TraPPE-UA acrylate force field can be found in Table 3. The CHx-CH2-CH2-CHy torsional parameters were taken from the TraPPE-UA alkane force field.12 The O-CH2CH2-O, O-CH2-CH2-CHx, and H-O-CH2-CHx torsional parameters were taken from the glycol17 and alcohol16 force fields, but the parameters were refit to the functional form used for the acrylate torsions (see eq 5). The remaining torsional parameters (see Table 4) were fit to the torsional energy calculated using Gaussian 0355 at the B3PW91/6-31++G(d,p) level of theory. For torsions involving an sp2 carbon connected to two additional interaction sites, e.g., CH2dCHsC(sp2)dO and CH2dCHsC(sp2)sO (see Figure 1), changing one torsional angle automatically changes the other. For simplicity, these torsions use the same potential parameters, but with an offset of 180° reflecting the planarity of the sp2 site. 3. Simulation Details Gibbs ensemble Monte Carlo57,58 (GEMC) simulations in the canonical (NVT) ensemble were carried out in order to compute the vapor-liquid coexistence curves (VLCC), saturated vapor pressures, and heats of vaporization for pure compounds. The canonical version of the Gibbs ensemble57 employs two separate periodic simulation boxes in thermodynamic contact, but without an explicit interface. In order to equilibrate the chemical potential, molecules are swapped between the boxes with the aid of coupled-decoupled configurational-bias Monte Carlo (CBMC).13,59-61 Volume exchanges are performed in order to equilibrate the pressure, and translational and rotational moves along with CBMC regrowths allow the system to reach thermal equilibrium. GEMC simulations were also performed in the isobaric-isothermal (NpT) ensemble to compute binary temperature-composition phase diagrams for the mixtures methyl acrylate/1-butanol and methyl acrylate/n-decane. The NpT Gibbs ensemble differs from the NVT version in that volume moves for each phase are performed separately in order for the system to equilibrate with an external pressure bath. Additionally,

angle type CHx-O-C(carbonyl) OsCdO18 O-CHx-CHy18 OdCsCHx18

21

θ0 (deg)

kθ/kB (K)

115 123 111 126

62500 40300 35300 40300

isobaric-isothermal (NpT) Monte Carlo simulations were performed at room temperature in order to compute the liquid density, cohesive energy, and radial distribution functions for single-component systems. For the acrylates and methacrylates, the system size consisted of 300 molecules, while for the alkenes (1,3-butadiene and isoprene) the system size was 512 molecules. A spherical cutoff of rcut ) 14 Å with analytic tail corrections62 was employed for the LJ interactions. The Coulombic interactions were computed using the Ewald summation technique with tin foil boundary conditions.63-65 The real-space cutoff was equal to rcut and the Ewald sum convergence parameter was set to 3.5/rcut. Liquidphase simulation box sizes were larger than twice rcut. Vaporphase box sizes were adjusted so that on average at least 20 molecules were present in the vapor phase. In cases where the saturated vapor pressure is low the vapor-phase box must be quite large in order for 20 molecules to be present. For large simulation boxes the Ewald sum becomes quite expensive; hence, once the vapor box length exceeds 50 Å we set rcut to be equal to 80% of half the box length. As there are few molecules in the vapor phase, this increases the efficiency by reducing the number of vectors used in the Ewald sum. Volume moves were accepted on average once every 10 Monte Carlo (MC) cycles (1 MC cycle ) N moves, where N is the number of molecules in the system). The fraction of swap moves was adjusted so that on average one swap move was accepted every 10-60 MC cycles. The higher values are used at lower temperatures, where the acceptance rate of swap moves is generally much lower. Of the remaining moves, approximately 1/3 were CBMC regrowths, with the rest divided equally between translational and rotational moves. For the binary systems, some trial and error is necessary in order to find a suitable overall composition and volume as the composition changes dramatically from one state point to the next as a function of temperature. For pure compounds two independent simulations were performed at every state point, while for binary mixtures four independent simulations were performed. At every state point at least 200 000 MC cycles of equilibration were performed for the acrylates and at least 100 000 MC cycles for the alkenes before starting the production runs. For the acrylates, the production runs consisted of 240 000 MC cycles divided into four blocks in order to compute statistical uncertainties, while for the alkenes the production runs consisted of 200 000 MC cycles divided into four blocks. Production runs for the binary mixtures consisted of 250 000 MC cycles divided into five blocks. Calculated statistical uncertainties are given as the standard error of the mean. The critical temperature (Tc) and density (Fc) were estimated using the scaling law66

Fliq - Fvap ) B(T - Tc)β

(6)

where β ) 0.325 and the law of rectilinear diameters67

1 (F + Fvap) ) Fc + A(T - Tc) 2 liq

(7)

6418


Maerzke et al.

TABLE 4: TraPPE-UA Force Field Parameters for Dihedral Angles torsion type

a0/kB (K)

a1/kB (K)

a2/kB (K)

a3/kB (K)

a4/kB (K)

CHxsOsCdO CHx-O-C(sp2)-C(sp2) CH2dCHx(sp2)sCdO CH2dCHx(sp2)sCsO OdCsC(sp2)sCH3 O-C(sp2)-C(sp2)-CH3 CHx-CHy-O-C(sp2) CH2dCHsCH)CH2 CH2dCHsC(sp2)sCH3 O-CH2-CH2-O17 O-CH2-CH2-CHx16 H-O-CH2-CHx16

1820.74 1820.74 823.03 823.03 195.19 195.19 2029.99 2034.58 1861.27 1258.09 893.21 368.58

-417.41 +417.41 -47.91 +47.91 -149.3 +149.3 -751.83 531.57 -349.97 0.0 176.62 209.80

-1373.14 -1373.14 -773.13 -773.13 164.38 164.38 -538.95 -1239.35 -1048.70 251.62 53.34 29.17

-30.19 +30.19 -1.99 +1.99 +24.12 -24.12 -22.10 460.04 -580.54 1006.47 769.93 187.93

0.0 0.0 0.0 0.0 28.45 28.45 -51.27 196.38 117.92 0.0 0.0 0.0

where Fliq and Fvap are the liquid and vapor densities, respectively. In the case of neutral molecules, the system size used in the present simulations is sufficiently large to reduce the error from finite size effects in the estimated critical properties to less than 1%.12 The normal boiling points (Tboil) were calculated using the Clausius-Clapeyron equation68

d ln Pvap ∆Hvap ) dT RT2

(8)

for the three temperatures nearest to Tboil. 4. Results and Discussion 4.1. Pure Compounds. GEMC simulations were performed for 1,3-butadiene and methyl acrylate in order to fit the CH(sp2) interaction site parameters. The Clausius-Clapeyron plot (see Figure 2) shows that the TraPPE force field accurately predicts the saturated vapor pressures for methyl acrylate and slightly overpredicts the saturated vapor pressures for 1,3-butadiene. Simulations of isoprene and methyl methacrylate were performed in order to fit the C(sp2) interaction site parameters. In this case, the saturated vapor pressure is overpredicted for isoprene and underpredicted for methyl methacrylate. As it is

Figure 2. Clausius-Clapeyron plot for 1,3-butadiene, isoprene, methyl acrylate, and methyl methacrylate. TraPPE 1,3-butadiene (triangles), isoprene (squares), methyl acrylate (diamonds), methyl methacrylate (circles), and experiment (solid lines). Experimental data for 1,3butadiene from the CRC Handbook,69 and experimental data for isoprene, methyl acrylate, and methyl methacrylate from the Antoine equation data from the NIST Chemistry WebBook.70

TABLE 5: Normal Boiling Points (K)a molecule 1,3-butadiene isoprene methyl acrylate ethyl acrylate n-butyl acrylate n-octyl acrylate 2-ethylhexyl acrylate 2-hydroxyethyl acrylate methyl methacrylate ethyl methacrylate n-butyl methacrylate a

TraPPE-UA CRC69 NIST70 Hoy72 avg exptl 260.82 297.23 354.77 375.21 422.31 505.21 493.04 461.11 377.44 395.24 437.02

268.7 307.2 353.9 372.6 418.2 502.2 490.2 464.2 373.7 390.2 433.2

268.6 307.3 353.5 372.8 420.1 490.2 373.5 391.7 434.7

283.8 353.3 372.6 421.8 491.5 375.0 -

268.65 307.25 353.57 372.66 420.03 502.2 490.63 464.2 374.06 390.95 433.95

Subscripts indicate uncertainties in the final digit.

not possible to simultaneously lower the vapor pressure for isoprene and raise the vapor pressure for methyl methacrylate, this is an acceptable compromise. It is well-known that the united-atom representation with accurate C-C bond length leads to a slight overprediction of the saturated vapor pressure for alkanes and alkenes.12-15,29,30,71 In order to validate the force field, additional GEMC simulations were performed for seven other acrylate and methacrylate monomers. A comparison of the normal boiling points with experimental values can be found in Table 5. The normal boiling points for 1,3-butadiene and isoprene are underpredicted by 3% relative to the average experimental value, which is entirely consistent with the TraPPE alkene model.15 The experimental data from Hoy72 for 1,3-butadiene were not included in any of the average values as the data from Hoy for 1,3-butadiene differ significantly from the other experimental data available for every quantity measured. The mean unsigned percent error (MUPE) for the acrylate and methacrylate monomers was 0.7%, which indicates excellent agreement with experiment. Liquid densities, shown in Table 6, were computed in the NpT ensemble at 293 K and 1 atm. Although TraPPE underpredicts the liquid density slightly for 6 of the 12 compounds, acceptable agreement with experiment is obtained with a MUPE (relative to the average experimental density) of 0.9%. Table 6 also includes values for the Dreiding force field with Mulliken (D/MUL) and ESP (D/ESP) partial charges.73 The MUPE for both D/MUL and D/ESP is 1.9% These results are graphically displayed in Figure 3, where we have plotted the simulated liquid densities against the average experimental values. If the simulated and experimental values agree, then the points in Figure 3 form a straight line with a slope of 1 and an intercept of 0. Again, it is evident that the TraPPE-UA force field yields more accurate predictions than the Dreiding force field. As the TraPPE result for 2-hydroxyethyl acrylate is considerably higher than the value reported by the CRC Handbook, some additional checks were performed. There is some confusion as



TABLE 6: Liquid Densities at 293 K (g/mL)a

a

molecule

TraPPE-UA

D/MUL73

D/ESP73

CRC69

Hoy72

avg exptl

1,3-butadiene isoprene methyl acrylate ethyl acrylate n-butyl acrylate n-octyl acrylate 2-ethylhexyl acrylate 2-hydroxyethyl acrylateb methyl methacrylated ethyl methacrylate n-butyl methacrylate

0.62271 0.67702 0.93732 0.90892 0.89233 0.88344 0.88413 1.0731 0.91824 0.89842 0.88463

0.64 0.93 0.9 0.84 0.95 0.93 0.89

0.65 0.92 0.89 0.85 0.92 0.90 0.88

0.6211 0.679 0.9535 0.9234 0.8898 0.881 0.880 1.011 0.9377 0.9135 0.8936

0.6460 0.9495 0.9156 0.8949 0.8810 1.106c 0.9303 -

0.6211 0.679 0.952 0.920 0.892 0.881 0.881 1.059 0.934 0.9135 0.8936

Subscripts indicate uncertainties in the final digit. b 296 K. c BASF MSDS.74 d 298 K.

TABLE 7: Relative Energies (in K) of Three 2-Hydroxyethyl Acrylate Conformers configuration

TraPPE

M06-2X/6-311+G(d,p)75

H-bond w/carbonyl oxygen linear chain (no H-bond) H-bond w/ether oxygen

0.00 1214.02 4133.36

0.00 1154.27 3865.18

TABLE 8: Critical Temperatures (K) and Densities (g/mL)a 1,3-butadiene n-butyl acrylate a

Figure 3. Scatter plot of the near-ambient liquid densities comparing the simulated results for the TraPPE force field (red triangles), Dreiding force field with Mulliken charges (green squares),73 and Dreiding force field with ESP charges (blue circles)73 against the average experimental values.69,72 The dashed, dot-dashed, and dotted lines are the linear regression for TraPPE, D/MUL, and D/ESP force fields. Data for 2-hydroxyethyl acrylate is not included in the regression due to the large uncertainty in the experimental value.

to the correct value for the liquid density, with the CRC Handbook reporting 1.011 g/mL69 and other sources74 reporting 1.106 g/mL. Judging from the TraPPE-UA predictions for the other acrylates, it appears likely that the value in the CRC Handbook is in error. Additionally, 2-hydroxylethyl acrylate has the added complication of intramolecular hydrogen bonding. In order to ensure that the TraPPE-UA force field correctly predicts the minimum energy configuration of 2-hydroxyethyl acrylate, the energies of the three possible conformers were computed using the Gaussian 03 program and compared to TraPPE-UA energies.55 The structures were first optimized at the M06-2X/MIDI!75 level of theory with the relevant dihedrals fixed. Next the dihedrals were unfrozen and the structures were further optimized at the M06-2X/6-311+G(d,p)75 level of theory. To compare the Gaussian results with TraPPE-UA, the energies were scaled with the lowest-energy structure set to zero. The results are summarized in Table 7. Clearly, the TraPPE-UA force field predicts the correct relative intramolecular energies for the 2-hydroxyethyl acrylate conformers. GEMC simulations were performed in order to compute the vapor-liquid coexistence curve (VLCC) for pure compounds.

TraPPE expt76 TraPPE expt77

Tc

Fc

4272 4251 6244 6443

0.2483 0.2455 0.2888 0.292


However, very little experimental data is available for acrylate monomers. In fact, the only data available are the critical points for 1,3-butadiene and n-butyl acrylate, where the critical point for n-butyl acrylate has been determined from vapor pressure measurements. The TraPPE critical temperature and density for 1,3-butadiene agree to within 0.5% and 1% of the experimental values, respectively. However, it is not surprising that TraPPE predicts the correct critical properties for 1,3-butadiene as these properties were considered when fitting the force field parameters. A more interesting case to consider is n-butyl acrylate,

Figure 4. Vapor-liquid coexistence curves for 1,3-butadiene (triangles) and n-butyl acrylate (circles). Experimental critical points (stars).76,77 Note that the experimental critical point for n-butyl acrylate has been determined from vapor pressure data, not a direct measurement.

6420


Maerzke et al.

TABLE 9: Heats of Vaporization at 298 K (kJ/mol)a molecule

TraPPE-UA

D/MUL73

D/ESP73

NIST70

PH78

Hoy72

Askadskii79

avg exptl

1,3-butadiene isoprene methyl acrylate ethyl acrylate n-butyl acrylate n-octyl acrylate 2-ethylhexyl acrylate 2-hydroxyethyl acrylate methyl methacrylate ethyl methacrylate n-butyl methacrylate

19.26 24.301 34.342 37.383 45.503 62.919 59.688 68.72 37.61 40.781 48.344

21 59 64 77 62 62 64

21 52 56 69 45 53 57

21.47 26.8 29.2 39.2 47.31 40.1 -

21 25 32 36 49 56 37 38 47

23.6 35.9 38.0 47.1 56.7 40.8 -

19 25 32 34 50 36 39 47

20.5 25.6 33.3 36.8 48.4 56.4 38.5 38.5 47

a


which was not one of the molecules used during parametrization. In this case, TraPPE predicts the critical temperature to within 3% of experiment, and the critical density to within 0.7% of experiment (see Figure 4). This is reasonable given that the experimental data is not a direct measurement, and the critical properties of methyl acrylate, the acrylate monomer used for parametrization, were not even considered during parameter fitting as there are no experimental values available. Figures showing the VLCCs and Clausius-Clapeyron plots and a table listing the critical constants for all acrylates and methacrylates are provided in the Supporting Information. Heats of vaporization were determined for the alkenes, acrylates, and methacrylates. For the low molecular weight compoundss1,3-butadiene, isoprene, and methyl acrylatesthe heat of vaporization at 298 K was determined via a GEMC simulation, where the heat of vaporization is calculated on the fly according to

j〉 ∆Hvap ) 〈Uvap - Uliq + Psat∆V

(9)

j are the instantaneous values of where Uvap, Uliq, Psat, and ∆V the molar internal energy of the vapor and liquid phases, the saturated vapor pressure, and the difference in molar volume of the liquid and vapor phases, respectively. For the remainder of the acrylate and methacrylate monomers, the vapor pressure is extremely low, making a GEMC simulation at 298 K impractical. For these molecules the heat of vaporization was calculated from NpT simulations of the liquid phase and simulations for an isolated molecule according to

∆Hvap ) Uiso - Uliq + RT

(10)

where Uiso and Uliq is the energy of an isolated molecule and the liquid phase, respectively. Given the importance of the solubility parameter for the description of the cohesiveness of polymers, many compilations report the values of the solubility parameter and not of the heat of vaporization for monomer units. Experimental heats of vaporization from Hoy,72 the Polymer Handbook (PH),78 and Askadskii79 as well as simulated heats of vaporization for the Dreiding model73 were determined by converting available Hildebrand solubility parameter data, where

j liq + RT ∆Hvap ) δH2V

(11)

j liq is the where δH is the Hildebrand solubility parameter and V molar volume of the liquid phase.

Figure 5. Scatter plot of the heats of vaporization comparing the simulated results for the TraPPE force field (triangles), Dreiding force field with Mulliken charges (squares),73 and Dreiding force field with ESP charges (circles)73 against the average experimental values.70,72,78,79 Values for the Dreiding model convertetd from solubility parameters; experimental values from the Polymer Handbook (PH),78 Hoy,72 and Askadskii79 also converted from solubility parameter data. The dashed, dot-dashed, and dotted lines are the linear regression for TraPPE, D/MUL, and D/ESP force fields.

The heats of vaporization from all experimental and simulation sources are listed in Table 9. Data for methyl acrylate from NIST70 is not included in the average experimental value as the vapor pressures do not agree with the Antoine equation. Again, data for 1,3-butadiene from Hoy72 is not included in the average experimental value. Clearly, the TraPPE force field performs much better at predicting heats of vaporization than the Dreiding force field with either charge model. The TraPPE force field has a MUPE of around 4% (relative to the average experimental value) whereas the Dreiding force field with Mulliken and ESP partial charges yields MUPEs of 42% and 23%, respectively. These results are graphically displayed in exp Figure 5 where ∆Hsim vap is plotted against ∆Hvap. Again, it is clear that the TraPPE force field performs much better than the Dreiding force field with either charge model. Radial distribution functions (RDFs) were computed from simulations in the NpT ensemble at 298 K and 1 atm. Although experimental RDFs are not available for these complex molecules, knowledge of their liquid structure can help to rationalize the differences in the properties of different acrylates. Examination of Figure 6 reveals that in many cases the RDF does not depend on the length of the alkyl side chain. For instance, the


Figure 6. Radial distribution functions for methyl acrylate (black), ethyl acrylate (red), n-butyl acrylate (green), n-octyl acrylate (blue), and ethylhexyl acrylate (magenta).

CH3(alkyl)-CH3(alkyl), CH3(alkyl)-CH2(sp2), and CH3(alkyl)O(ether) RDFs do not show a significant difference between methyl, ethyl, n-butyl, n-octyl, and even the branched 2-ethylhexyl acrylate. However, the CH3(alkyl)-O(carbonyl) RDF has a larger peak for smaller side chains. This could be the result of steric effects, the parametrization of the model, or a combination of both. The carbonyl oxygen has a partial negative charge, while for methyl acrylate the CH3 group has a partial positive charge. As the side chain becomes longer, the methyl group moves further away from any charge sites and is no longer attracted to the carbonyl oxygen. To see the effects of chain length on the charge site, we computed the CHx(R-ether)O(carbonyl) RDF, where CHx(R-ether) is the CHx group with a partial positive charge adjacent to the ether oxgyen. In this case, steric effects are clear. The peak height is greatly diminished for all of the acrylates with longer side chains, and there is a noticeable decrease in peak height for 2-ethylhexyl acrylate versus n-octyl acrylate, its linear isomer. Also of interest is the difference between the CH3(alkyl)-O(carbonyl) and CH3(alkyl)-O(ether) RDFs. As both of these oxygens carry a partial negative charge, one might expect that both CH3(alkyl)-O RDFs would be similar. However, a significantly larger peak is found for the carbonyl oxygen. This is likely caused by the larger accesssible surface area of the carbonyl oxygen, as also concluded for polymeric solvation.50 Similarly, many of the methacrylate RDFs do not show much dependence on alkyl side chain length; as examples the CH3(alkyl)-CH3(alkyl) and CH3(methyl)-CH3(methyl) (where “alkyl” indicates the R group connected to the ether linker and “methyl” indicates the methyl group connected to the C(sp2) atom, see Figure 1) are shown in Figure 7. In agreement with the observation for the acrylate monomers, a difference in peak heights for the CH3(alkyl)-O(carbonyl) and CHx(R-ether)O(carbonyl) RDFs is evident for methacrylate monomers depending on the length of the alkyl side chain (see Figure 7). Again, this is the result of the partial charges and steric effects. We also examined the RDFs involving the hydroxyl hydrogen and the various oxygen atoms for 2-hydroxyethyl acrylate (see Figure 8). The peak at r ≈ 1.8 Å is indicative of the formation


Figure 7. Radial distribution functions for methyl methacrylate (black), ethyl methacrylate (red), and n-butyl methacrylate (green).

Figure 8. Comparison of radial distribution functions for hydrogen bonding in 2-hydroxyethyl acrylate. The O(hydroxyl)-H(hydroxyl), O(carbonyl)-H(hydroxyl), and O(ether)-H(hydroxyl) are shown in black, red, and green, respectively.

of a hydrogen bond. Here we notice a significantly larger peak for the hydroxyl oxygen than for either the carbonyl or ether oxygen. We also obtain a small peak for the carbonyl oxygen, which indicates that if a hydrogen bond forms with one of the acrylate oxygen atoms, it is with the carbonyl oxygen rather than the ether oxygen, which does not appear to act as hydrogen bond acceptor. This can be explained by two factors: the greater accessibility of the carbonyl oxygen and the smaller charge of the ether oxygen, which is reduced by a factor of 2 for the ether oxygen in an acrylate compared to the ether oxygen in a dialkyl ether.17 4.2. Binary Mixtures. GEMC simulations in the NpT ensemble at 1 atm were performed in order to compute the temperature-composition VLCC for binary mixtures of methyl acrylate with polar and nonpolar molecules. The binary VLCC and separation factor (that indicates the difference in the compositions of the liquid and vapor phases) for the mixtures methyl acrylate/1-butanol and methyl acrylate/n-decane (see

6422


Figure 9. Temperature-composition phase diagram (top) and separation factor (bottom) for the binary mixture methyl acrylate/1-butanol. x and y indicate the compositions of the liquid and vapor phases, respectively. Experiment81 (squares) and TraPPE (circles). Dotted black lines are only included to guide the eye and do not represent experimental values.

Figures 9 and 10) agree well with experiment, with the majority of the error in the VLCC coming from the underprediction of the boiling point of 1-butanol and n-decane by the TraPPE-UA alcohol16 and alkane12 force fields. The experimental separation factor for the mixture with the polar alcohol molecule is better reproduced than for the mixture with the nonpolar alkane molecule. The larger deviation for the latter mixture indicates that the unlike interactions are somewhat too favorable, presumably because Lennard-Jones well depths for the new acrylate pseudoatoms may be too deep because they compensate for a lack of explicit partial charges in 1,3-butadiene and isoprene and for the partial charges in methyl acrylate and methyl methacrylate determined from the continuum solvation model that may be slightly too small in magnitude.24,80 The structural properties of these mixtures were also examined. To examine the microheterogeneity, the local mole fraction enhancement (LMFE) was determined as a function of distance from a chosen molecule. The local mole fraction is the average number of molecules of a given type up to a distance r from the chosen molecule divided by the total number of molecules of any type, where the number of molecules of each type is determined from the center of mass number integral (NI). To obtain the local mole fraction enhancement, the local mole

Maerzke et al.

Figure 10. Temperature-composition phase diagram (top) and separation factor (bottom) for the binary mixture methyl acrylate/n-decane. x and y indicate the compositions of the liquid and vapor phases, respectively. Experiment82 (squares) and TraPPE (circles). Dotted black lines are only included to guide the eye and do not represent experimental values.

fraction is divided by the bulk mole fraction; that is, for a binary system with molecule types a and b

LMFEaa(r) )

[

NIaa(r) 1 xa NIaa(r) + NIab(r)

]

(12)

For the methyl acrylate/1-butanol mixtures, we see a slight enhancement of the 1-butanol-1-butanol local mole fraction for low temperature/low concentrations of 1-butanol (see Figure 11). The largest LMFE occurs for T ) 360 K, xbutanol ) 0.47, and at T ) 380 K, xbutanol ) 0.94 there is a neglible LMFE as the system consists almost entirely of 1-butanol. Likewise, we see a slight enhancement of the methyl acrylate-methyl acrylate local mole fraction, which decreases with increasing methyl acrylate concentration. At the lowest concentration of methyl acrylate (T ) 380 K, xacrylate ) 0.06), the LMFE does not converge to unity. However, for this system the total number of methyl acrylate molecules is quite small (N ≈ 10), so the fact that the number integral goes to N - 1 at large distances instead of N causes a noticeable difference in the LMFE at large distances. We can use the RDFs to further probe the microstructure of these systems. In Figure 12 we show the O-H RDFs for the


Figure 11. Comparison of local mole fraction enhancements (LMFE) of 1-butanol-1-butanol (top) and methyl acrylate-methyl acrylate (bottom) for various temperatures/concentrations of the methyl-acrylate/ 1-butanol mixture. T ) 360 K, xbutanol ) 0.47 (black), T ) 365 K, xbutanol ) 0.71 (red), T ) 370 K, xbutanol ) 0.81 (green), T ) 375 K, xbutanol ) 0.88 (blue), and T ) 380 K, xbutanol ) 0.94 (magenta). The corresponding number integrals are shown in the Supporting Information.


Figure 13. Comparison of local mole fraction enhancements (LMFE) of n-decane-n-decane (top) and of methyl acrylate-methyl acrylate (bottom) for various temperatures/concentrations of the methyl acrylate/ n-decane mixture. T ) 370 K, xdecane ) 0.43 (black), T ) 385 K, xdecane ) 0.67 (red), T ) 400 K, xdecane ) 0.81 (green), T ) 415 K, xdecane ) 0.90 (blue), and T ) 425 K, xdecane ) 0.95 (magenta).

Figure 14. Comparison of CH3(acrylate)-CH3(other) radial distribution functions for the methyl acrylate/1-butanol and methyl acrylate/ n-decane mixtures. Colors as in Figures 11 and 13.

Figure 12. Comparison of radial distribution functions for hydrogen bonding in mixtures of methyl acrylate and 1-butanol. Notice the difference in scale on the y-axis. Colors as in Figure 11.

hydroxyl, carbonyl, and ether oxygens. Clearly, 1-butanol prefers to hydrogen bond with itself rather than with either of the methyl acrylate oxygens. This is most pronounced at the lowest temperature/lowest concentration of 1-butanol, where the 1-butanol-1-butanol LMFE is largest. However, if 1-butanol does form a hydrogen bond with a methyl acrylate molecule, that hydrogen bond is with the carbonyl, not the ether, oxygen. This is again likely due to the larger accessible surface area of the carbonyl oxygen.50 For the methyl acrylate/n-decane system, we see slight depletion of the n-decane-n-decane local mole fraction (see Figure 13). Based on the phase diagram, one might expect to see a enhancement instead of a depletion. However, methyl acrylate is a much more compact molecule than n-decane, so it is able to fill the spaces in between n-decane molecules, resulting

in a depletion of the n-decane-n-decane LMFE. The effect is the largest at the lowest temperature/lowest concentration of decane (T ) 370 K, xdecane ) 0.43), and diminishes as the concentration of decane increases. We also see a slight enhancement of the methyl acrylate-methyl acrylate local mole fraction. The effect increases as the concentration of methyl acrylate decreases, with the smallest enhancement occurring at the lowest temperature/highest concentration of methyl acrylate. Again, at the lowest concentration of methyl acrylate (T ) 425 K, xacrylate ) 0.05) the LMFE does not normalize to one as the total number of methyl acrylate molecules in this system is quite small. It is interesting to note that the LMFE for the methyl acrylate/ndecane mixtures persists over larger distances than the corresponding LMFE in the methyl acrylate/1-butanol mixtures; i.e. a signal for longer-range composition microheterogeneity for the mixture with the nonpolar species. Upon examining the CH3(acrylate)-CH3(other) RDFs in Figure 14, we notice that these RDFs are strikingly similar for both mixtures (i.e., do not show a significant dependence on the polarity of the second compound in the mixture) and that

6424


there is no significant temperature/concentration dependence for either mixture. 5. Conclusions The TraPPE-UA force field has been successfully extended to include acrylate and methacrylate monomers through the introduction of two new interaction sites and new torsional potentials. Excellent agreement with experiment was obtained, the average errors in the liquid density and normal boiling point of around 1%. The binary vapor-liquid equilibria for the systems methyl acrylate/1-butanol and methyl acrylate/n-decane show good agreement with the experimental separation factors; however, the underprediction of the 1-butanol and n-decane boiling points results in a slightly shifted coexistence curve. Moderate composition microheterogeneity is observed for both mixtures. Acknowledgment. Financial support from the National Science Foundation (CBET-0553911 and CBET-756641) and the 3M Company is gratefully acknowledged. Part of the computer resources were provided by the Minnesota Supercomputing Institute. Supporting Information Available: Tables listing the numerical data for all VLCC of all neat compounds and for the two binary mixtures and the extrapolated critical temperatures and densities. Figures showing the VLCCs and ClausiusClapeyron plots for all acrylates and methacrylates studied here. Figures showing the number integrals for the binary mixtures. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Beuermann, S.; Buback, M.; Drache, M.; Nelke, D.; Schmidt-Naake, G. e-Polym. 2004, 2, 1–7. (2) Allinger, N. L.; Yuh, Y. H.; Lii, J.-H. J. Am. Chem. Soc. 1989, 111, 8551–8566. (3) Korolev, G. V.; Il’in, A. A.; Solov’ev, M. E.; Mogilevich, M. M.; Evplonova, E. S. Polym. Sci. Ser. A 2001, 43, 1055–1059. (4) Korolev, G. V.; Il’in, A. A.; Solov’ev, M. E.; Srybnyi, A. V.; Mogilevich, M. M.; Evplonova, E. S. Polym. Sci. Ser. A 2002, 44, 1155– 1161. (5) Allinger, N. L. J. Am. Chem. Soc. 1977, 99, 8127–8134. (6) Burkert, U.; Allinger, N. L. Molecular Mechanics; American Chemical Society: Washington, DC, 1982. (7) Curco, D.; Aleman, C. J. Chem. Phys. 2004, 121, 9744–9752. (8) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765– 784. (9) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230–252. (10) Han, X.; Chai, L.; Shanks, R.; Pavel, D. Polym. Int. 2006, 55, 1323– 1329. (11) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897–8909. (12) Martin, M. G.; Siepmann, J. I. J. Chem. Phys. 1998, 102, 2569– 2577. (13) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4508– 4517. (14) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 5370–5379. (15) Wick, C. D.; Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 8008–8016. (16) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 3093–3104. (17) Stubbs, J. M.; Pottoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2004, 108, 17596–17605. (18) Kamath, G.; Cao, F.; Potoff, J. J. J. Phys. Chem. B 2004, 108, 14130–14136. (19) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 18974–18982. (20) Lubna, N.; Kamath, G.; Potoff, J. J.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 24100–24107.

Maerzke et al. (21) Kamath, G.; Robinson, J.; Potoff, J. J. Fluid Phase Equilib. 2006, 240, 46–55. (22) Zhao, X. S.; Chen, B.; Karaborni, S.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 5368–5374. (23) Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2007, 111, 10790–10799. (24) Rai, N.; Siepmann, J. I. Unpublished work. (25) Chen, B.; Xing, J.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 2391–4201. (26) Potoff, J. J.; Siepmann, J. I. AIChE J. 2001, 47, 1676–1682. (27) Zhang, L.; Siepmann, J. I. Theor. Chem. Acc. 2006, 115, 391–397. (28) Zhang, L.; Siepmann, J. I. Unpublished work. (29) Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Rousseau, B.; Fuchs, A. H. J. Chem. Phys. 2000, 112, 5499–5510. (30) Bourasseau, E.; Haboudou, M.; Boutin, A.; Fuchs, A. H.; Ungerer, P. J. Chem. Phys. 2003, 118, 3020–3034. (31) Perez-Pellitero, J.; Bourasseau, E.; Demachy, I.; Ridard, J.; Ungerer, P.; Mackie, A. D. J. Phys. Chem. B 2008, 112, 9853–9863. (32) Unlu, O.; Gray, N. H.; Gerek, Z. N.; Elliott, J. R. Ind. Eng. Chem. Res. 2004, 43, 1788–1793. (33) Baskaya, F. S.; Gray, N. H.; Gerek, Z. N.; Elliott, J. R. Fluid Phase Equilib. 2005, 236, 42–52. (34) Vahid, A.; Sans, A. D.; Elliott, J. R. Ind. Eng. Chem. Res. 2008, 47, 7955–7964. (35) Kelkar, M. S.; Rafferty, J. L.; Siepmann, J. I.; Maginn, E. J. Fluid Phase Equilib. 2007, 260, 218–231. (36) Rai, N.; Wagner, A. J.; Ross, R. B.; Siepmann, J. I. J. Chem. Theor. Comput. 2008, 4, 136–144. (37) Chen, B.; Siepmann, J. I. J. Am. Chem. Soc. 2000, 122, 6464– 6467. (38) Sun, L.; Siepmann, J. I.; Klotz, W. L.; Schure, M. R. J. Chromatogr. A 2006, 1126, 373–380. (39) Wick, C. D.; Siepmann, J. I.; Schure, M. R. Anal. Chem. 2002, 74, 3518–3524. (40) Wick, C. D.; Siepmann, J. I.; Klotz, W. L.; Schure, M. R. J. Chromatogr. A 2002, 954, 181–190. (41) Rafferty, J. L.; Zhang, L.; Siepmann, J. I.; Schure, M. R. Anal. Chem. 2007, 105, 1411–1417. (42) Rafferty, J. L.; Siepmann, J. I.; Schure, M. R. Anal. Chem. 2008, 80, 6214–6221. (43) Rafferty, J. L.; Siepmann, J. I.; Schure, M. R. J. Chromatogr. A 2008, 1204, 20–27. (44) Chen, B.; Siepmann, J. I.; Klein, M. J. Am. Chem. Soc. 2002, 124, 12232–12237. (45) Bresme, F.; Chacon, E.; Tarazona, P.; Tay, K. Phys. ReV. Lett. 2008, 101, 056102. (46) Smit, B. Chem. ReV. 2008, 108, 4125–4184. (47) Yazaydin, A. O.; Thompson, R. W. J. Phys. Chem. B 2006, 110, 14458–14462. (48) Dubbeldam, D.; Galvin, C. J.; Walton, K. S.; Ellis, D. E.; Snurr, R. Q. J. Am. Chem. Soc. 2008, 130, 10884–10885. (49) Wick, C. D.; Theodorou, D. N. Macromolecules 2004, 37, 7026– 7033. (50) Wick, C. D.; Siepmann, J. I.; Theodorou, D. N. J. Am. Chem. Soc. 2005, 127, 12338–12342. (51) Rissanou, A. N.; Peristeras, L. D.; Economou, I. G. Polymer 2007, 48, 3883–3892. (52) Li, H.; Curro, J. G.; Wu, D. T.; Habenschuss, A. Macromolecules 2008, 41, 2694–2700. (53) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, A. Intermolecular Forces: Their Origin and Determination; Pergamon Press: Oxford, UK, 1987. (54) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Theor. Comput. Chem. 2005, 1, 1133–1152. (55) Gaussian 03, ReVision C.01; Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian Inc.: Wallingford, CT, 2003. (56) Chamberlin, A. C.; Kelly, C. P.; Thompson, J. D.; Xidos, J. D.; Li, J.; Hawkins, G. D.; Winget, P. D.; Zhu, T.; Rinaldi, D.; Liotard, D. A.;

TraPPE-UA Force Field for Acrylates Cramer, C. J.; Truhlar, D. G.; Frisch, M. J. MN-GSM, version 6.0, University of Minnesota, Minneapolis, MN 55455-0431, 2006. (57) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813–826. (58) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527–545. (59) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59–70. (60) de Pablo, J. J.; Laso, M.; Siepmann, J. I.; Suter, U. W. Mol. Phys. 1993, 80, 55–63. (61) Vlugt, T. J. H.; Martin, M. G.; Smit, B.; Siepmann, J. I.; Krishna, R. Mol. Phys. 1998, 94, 727–733. (62) Wood, W. W.; Parker, F. R. J. Chem. Phys. 1957, 27, 720–733. (63) Ewald, P. Ann. Phys. 1921, 64, 253–287. (64) Madelung, E. Phys. Z. 1918, 19, 524–532. (65) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, UK, 1987. (66) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Oxford University Press: New York, 1989. (67) Rowlinson, J. S.; Swinton, F. L. Liquids and Liquid Mixtures; Butterworth: London, 1982. (68) Atkins, P. W. Physical Chemistry; Oxford University Press: Oxford, UK, 2002. (69) Lide, D. R., Baysinger, G., Berger, L. I., Goldberg, R. N., Kehiaian, H. V., Kuchitsu, K., Rosenblatt, G., Roth, D. L., Zwillinger, D., Eds.; CRC Handbook of Chemistry and Physics, 88th ed.; CRC Press: Boca Raton, FL, 2007. (70) NIST Chemistry WebBook, NIST Standard Reference Database Number 69, Linstrom P. J., Mallard, W. J., Eds.; National Institute of

J. Phys. Chem. B, Vol. 113, No. 18, 2009 6425 Standards and Technology: Gaithersburg MD, June 2005 (http:// webbook.nist.gov). (71) Errington, J. R.; Panagiotopoulos, A. Z. J. Phys. Chem. B 1999, 103, 6314–6322. (72) Hoy, K. L. J. Paint Technol. 1970, 42, 76–118. (73) Belmares, M.; Blanco, M.; Goddard, W. A., III; Ross, R. B.; Caldwell, G.; Chou, S.-H.; Pham, J.; Olofson, P. M.; Thomas, C. J. Comput. Chem. 2004, 25, 1814–1826. (74) 2-Hydroxyethyl acrylate, BASF Technical Information TI/ED 1632e, December 2004. (75) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2007, 120, 215–241. (76) Tsonopoulos, C.; Ambrose, D. J. Chem. Eng. Data 1996, 41, 646– 656. (77) Steele, W. V.; Chirico, R. D.; Knipmeyer, S. E.; Nguyen, A.; Smith, N. K. J. Chem. Eng. Data 1996, 41, 1285–1302. (78) Brandrup, J., Immergut, E. H., Grulke, E. A., Eds.; Polymer Handbook, 4th ed.; Wiley: New York, 1999. (79) Askadskii, A. Physical Properties of Polymers-Prediction and Control; Gordon and Breach: London, 1996. (80) Stubbs, J. M.; Siepmann, J. I. J. Phys. Chem. B 2002, 106, 3968– 3978. (81) Chubarov, G. A.; Danov, S. M.; Brovkina, G. V. J. Appl. Chem. U.S.S.R. 1976, 49, 1648–1649. (82) Chubarov, G. A.; Danov, S. M.; Logutov, V. I.; Tomashchuk, V. I. J. Appl. Chem. U.S.S.R. 1984, 57, 1891–1893.

JP810558V

TraPPE-UA Force Field for Acrylates and Monte Carlo Simulations for

Recommend Documents