Anisotropic United Atom Model Including the Electrostatic Interactions

Oct 10, 2007 - ... Amer Wafai , Simon Stephan , Maximilian Kohns , Steffen Reiser , Stephan Deublein , Martin Horsch , Hans Hasse , Jadran Vrabec...
0 downloads 0 Views 152KB Size
Download PDF
15942

J. Phys. Chem. C 2007, 111, 15942-15951

Anisotropic United Atom Model Including the Electrostatic Interactions of Methylbenzenes. II. Transport Properties† Carlos Nieto-Draghi,* Patrick Bonnaud, and Philippe Ungerer IFP, 1-4 AVenue de Bois-Pre´ au, 92852 Rueil-Malmaison, Cedex, France ReceiVed: May 15, 2007; In Final Form: August 21, 2007

A new potential model for monoaromatics based an anisotropic united atoms model and three electrostatic charges is used to simultaneously represent the dipole moments, the vapor-liquid equilibrium, the transport properties, and the structural properties while keeping computational effort minimal. The present article focuses on the transport properties and dynamic structural properties. Shear viscosities and diffusion coefficients are computed from equilibrium molecular dynamics in a large range of pressures (0.1-500 MPa) and temperatures (298-363 K). Viscosity shows a significant improvement with respect to the previous nonpolar models, as average deviations on toluene, p-xylene, m-xylene, and o-xylene stay below 9%. Self-diffusion coefficients and reorientational diffusion coefficients of toluene are also in good agreement with experimental observations. Thermal conductivity, evaluated by nonequilibrium molecular dynamics, reveals maximum deviations of 8% on benzene, toluene, and xylene isomers, which is better than those of other potentials investigated. By its lower computation time and equivalent accuracy compared with all atoms models, the proposed AUA potential is a good option for investigating complex problems such as the prediction of fuel viscosity in the highpressure conditions encountered in the injectors of diesel engines.

1. Introduction Methylbenzenes are important organic solvents and reactants used in many industrial applications. In particular, xylene isomers are important constituents of gasoline fuels due to their high octane number. The development of a new generation of low-consumption/high-efficiency motor engines,1 using high pressure injection systems, requires the knowledge of physical properties of fuels at extreme conditions (close to 250 MPa).2 At such conditions, experimental measurements are scarce, unsafe, and expensive, and classical correlative equations are not reliable to predict fluid properties. Consequently, molecular simulations techniques become one of the reliable alternatives to obtain the necessary information on the physical properties of the whole family of hydrocarbons. Though a great effort was devoted to develop intermolecular potentials able to well reproduce thermodynamic properties, there is a need for accurate intermolecular potential models that also reproduce dynamic properties with good accuracy. Several force fields have been proposed in the past decade.3-20 These models include all atoms (AA), united atoms (UA), anisotropic united atoms (AUA), and exponential 6 force fields, among others. The AUA approach has been successfully applied to different families of hydrocarbons, from linear21 and branched alkanes22 and alkenes23 and mono- or polyaromatics compounds24,25 to other sulfur26,27 compounds. Previous comparative works show that UA and AUA models without electrostatic interactions are capable of reproducing liquid-vapor equilibrium properties of alkylbenzenes quite accurately, while keeping the computational cost at a minimum.18,19 However, the price paid for this simplification is an erroneous description of the liquid structures and intermolecular †

Part of the “Keith E. Gubbins Festschrift”. * To whom correspondence should be addressed. E-mail: Carlos.nieto@ ifp.fr.

Figure 1. Local system of the reference (a) and relative angle orientation ξ of dipole moment vectors between toluene molecules (b); eˆ i⊥ is a unit vectors perpendicular to the plane of the toluene molecules i ) 1,2. (c) Schematic representation of the AUA molecules for toluene, o-xylene, m-xylene, and p-xylene (molecules are centered in their principal axes of rotation).

pair energies with respect to ab initio quantum chemistry calculations.28 In addition, these simplified potentials present problems when they are used to simulate alkene-aromatic mixtures since their interactions are overestimated.12,29 Although the electrostatic contribution in this kind of aromatics represents approximately less than 10% of the total intermolecular energy, this contribution is not required to reproduce thermodynamic

10.1021/jp073715y CCC: $37.00 © 2007 American Chemical Society Published on Web 10/10/2007

AUA Model of Methylbenzenes

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15943

TABLE 1: Different Potential Parameters and Centers of Force of Alkylbenzene Moleculesa (See Figure 1 for Details about the Position of Electrostatic Charges); All Alkylbenzene Molecules Have the Same Quadrupolar Moment as that of Benzene model ch-AUA

moleculea benzene toluene o-xylene m-xylene p-xylene model AUA nonpolar (ref 16) molecule benzene toluene o-xylene m-xylene p-xylene

group

σ (Å)

/k (K)

q (e)

δ (Å)b

Q(×10-40 C‚m2)e -20.8

3.361

75.6

0.315

3.361 3.6072

35.43 120.15

0.000 0.216

CH (rCC ) 1.4 Å) C CH3 (rC-CH3 ) 1.535 Å) +2q (y ) z ) 0.0, x ) δµ) -q (x ) y ) 0.0, z ) +0.4 and z ) -0.4 Å)c Ixxf 76.55113 76.55113 140.3724 267.8843 76.6165 group CH (rCC ) 1.4 Å)

Iyyf 76.55113 185.1308 216.3667 123.1998 331.6400 σ (Å) 3.2464

Izzf 153.1022 261.6820 356.7391 391.0841 408.2565 /k (K) 89.415

C CH3 (rC-CH3 ) 1.51 Å) Ixxf 76.55113 76.55113 140.3724 267.8843 76.6165

3.0648 3.6072 Iyyf 76.55113 185.1308 216.3667 123.1998 331.6400

42.079 120.15 Izzf 153.1022 261.6820 356.7391 391.0841 408.2565

8.13 -4.065 δµ (Å) 0.0 9.6 × 10-3 1.63 × 10-2 7.679 × 10-3 0.0 q (e)

δµ (Å)

µ (D) 0.000 0.375 0.640 0.300 0.000 δ (Å)b 0.4071

µexp (D)d 0.0 0.359 0.630 0.299 0.0 Q(×10-40 C‚m2)e 0.0

0.216 µ (D) 0.0 0.0 0.0 0.0 0.0

µexp (D) 0.0 0.359 0.630 0.299 0.0

a All models respect the planar hexagonal geometry for the aromatic ring. b Distance from the carbon atom. c Distance from the center in the axis perpendicular to the plane of the molecule. d Experimental dipole moment from the DIPPR database.42 e The experimental value44,45 range between -28.3 and -33.3 × 10-40 C‚m2. f Principal moments of inertia (in g/mol‚Å2) are computed considering the masses of the united atom groups with respect to the reference center of mass of the molecule, as shown in Figure 1c.

properties well. However, as it was demonstrated in a previous work, this is not the case for transport properties where this contribution is crucial.30 The process of the crystallization of benzene,31,32 the interaction of benzene with polar molecules,33-37 or the adsorption of benzene into charged aluminosilicate solids38-40 is a clear example where electrostatic interactions are needed. In addition, as we have shown, a simple nonelectrostatic model such as the AUA with 6 interaction sites presents an absolute average deviation of about 25% when shear viscosity is calculated. This fact is not surprising since the electrostatic interaction tends to increase the friction at high density due to slight repulsion between molecules at short distances. We intend to extend the AUA approach used to develop the charged model of benzene to describe different methylbenzene molecules such as toluene, o-xylene, m-xylene, and p-xylene. In this particular case, it is also required to include the representation of the dipolar moment on the electrostatic representation of these molecules. The objective of this work is to show how an optimal and accurate intermolecular potential model for methylbenzenes, including a minimum number of point charges to describe electrostatic interactions, is able to reproduce different transport properties well. The potential parameters of the new model are developed in an accompanying work of this series, where liquid-vapor and structural information was employed to optimize the Lennard-Jones coefficients of the aromatic carbon C, linked to the methyl group, in a standard procedure.41 The potential parameters of the CH group are based on a previous work, where we have developed a new AUA potential for benzene that accounts for electrostatic interactions of this molecule,30 and the CH3 group parameters are unchanged with respect to n-alkanes. We have computed several dynamic

properties such as reorientational and self-diffusion coefficients, shear viscosity, and thermal conductivity at different thermodynamic conditions. This work is organized as follows. In the forthcoming section 2, we present the simulation methods and algorithms used to compute the different dynamic properties. In section 3, we present the simulation details, followed by the section 4 where the main results obtained for the different properties are presented and discussed. Finally, section 5 is devoted to the main conclusions of our work. 2. Simulation Methods 2.1. Intermolecular Potential Model. The effective dispersion-repulsion interactions between two atoms or united atoms (i and j) of different molecules are represented by the LennardJones 6-12 equation

[( ) ( ) ]

ULJ ) 4ij

σij rij

12

-

σij rij

6

(1)

where σij and ij are the Lennard-Jones interaction parameters between sites i and j on different molecules and rij ≡ rj - ri is the separation distance between sites i and j. Cross interactions between centers of force of different type were treated using standard Lorentz-Berthelot combining rules. The electrostatic energy was obtained by summing the pairwise Coulombic interactions between the partial charges belonging to the different molecules

Uelec )

1 4π0

∑ i*j

qiqj rij2

(2)

15944 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Nieto-Draghi et al.

where qi and qj are the charges of different centers of interaction and 0 is the dielectric constant of vacuum. Since all molecules treated here were assumed to be rigid, no stretching, bending, or torsional energy was included. The electrostatic charges were taken from our previous work on the benzene molecule, and the dipole moment of each molecule was adjusted to the experimental value42 by displacement of the central positive charge by a distance of δµ in the plane of the molecules to match the experimental direction and value of this property. The negative charges, located at 0.4 Å above and below the ring plane on its hexagonal symmetry axis, were unchanged with respect to the previous AUA model of benzene30 (see Figure 1 for details). The equation |µ| ) q·δµ was used to determine the separation distance δµ between the projection of the coordinates of the negative charges (-q, out of the plane of the molecule) on the xy plane and the central positive charge q ) 2q in each case. In this way, we simultaneously respected the quadrupolar moment of the benzene AUA model of -20.8 × 10-40, which is slightly lower than the experimental one (i.e., -28.31 × 10-40 C‚m2 from the work of Dennis et al.43 and -33.3 × 10-40 C‚m2 from the work of Vbrancich et al.44,45), and added the necessary dipole moment for toluene, o-xylene, and m-xylene molecules. The potential parameters of the different molecules used, developed in an accompanying work, can be seen in Table 1. 2.2. Dynamic Properties and Algorithms. 2.2.1. Reorientational Dynamics. The rotational dynamics of methylbenzene molecules can be analyzed though temperature-dependent reorientational relaxation times, experimentally reported by 13C NMR dipolar spin-lattice relaxation rates, and cross correlation rates between the dipolar relaxation mechanism and the relaxation by the chemical-shift anisotropy of the molecule to compute the rotational diffusion constant around the principal axes of rotation.46,47 In addition, deuterated sample molecules are used to measure T1 relaxation times of the CH dipole moment vectors in the molecule, which permits the computation of the correlation times τ. We have computed the reorientational correlation time of several unit vectors, depending on the geometry of the molecule and the available experimental data to compare with. For the case of toluene, three vectors were defined to mimic the CH bond in the ortho, meta, and para positions with respect to the methyl group. As in a previous work, we have used a unit vector eˆ i, describing the displacement δAUA of the AUA intermolecular potential model, which is equivalent to the vector describing the geometrical position of all of the hydrogen atoms not explicitly described in this kind of united-atom-type force field. In this way, we are able to compare the reorientation of these bonds with the spin-lattice 13C-1H/2H NMR experimental data. For the case of the o-xylene and m-xylene isomers, we have monitored, in addition, one vector in the direction of the dipole moment (µ) and another vector aligned between the two methyl groups for the case of p-xylene. Finally, another vector was defined perpendicular (eˆ ⊥) to the plane of the aromatic ring and another one orthogonal to the dipole moment of the molecule (see Figure 1b and c for details of all of the vectors defined). Two orientational autocorrelation functions can be defined

Cli(t) ) 〈Pl(eˆ i(t)‚eˆ i(0))〉

(3)

where Pl is a Legendre polynomial with l ) 1 and 2 for first and second order, respectively. These functions are commonly fitted to Kohlrausch-Williams-Watts exponential functions48,49

βil

l Ci,KWW (t) ) exp[(-t/Rli) ]

(4)

for which the relaxation time can be obtained straightforwardly

τli )

R

( )

lm 1 l (t)dt ) Γ ∫0∞ Ci,KWW βlm βlm

(5)

Experimentally, one has access to correlation time τlm, which is normalized over the spherical harmonics functions Ylm. Rotational diffusion constants are computed through the integration of the angular velocity correlation function

Di )

∫0∞ 〈ωi(t)‚ωi(0)〉dt )

kBT Ii

∫0∞ Cω (t)dt i

(6)

where Di and Ii are, consequently, the rotational diffusion constant and the moment of inertia around the ith principal axis of rotation i ) x,y,z; kB is the Boltzmann constant, and Cωi(t) is the normalized angular autocorrelation function. Contrary to the benzene molecule, whose rotational diffusion constant can be related to the correlation time for the symmetric benzene diffuser, τ2m ) 1/(6D⊥ + m(D|| - D⊥)), we cannot approximate the diffusion perpendicular to the plane of the molecules with a single coefficient (D⊥) since alkylbenzene molecules present an important degree of anisotropy in their reorientation. In addition, an approach such as Di ∼ 1/6τi2 is not appropriate for the computation of the rotational diffusion constant; consequently, we have employed eq 6 to compute this property. We have compared our results with relaxation times and rotational diffusion constants obtained with 13C and 13C1H/2H spin-lattice NMR experiments for pure toluene46,50 and p-xylene51 at different temperatures and at ambient pressure. 2.2.2. Self-Diffusion. The self-diffusion coefficient of each molecule has been computed at a constant number of molecules and temperature NVT through the mean square displacement, according to the well-known Einstein relation56

D ) lim tf∞

1 〈|r(t) - r(0)|2〉 6t

(7)

where r(t) is the center of mass vector position of each molecule at time t, r(0) is the position at t ) 0, and D is the self-diffusion coefficient. 2.2.3. Shear Viscosity. The shear viscosity can be calculated for a constant number of particles and a fixed temperature, imposing either volume or pressure in thermodynamic equilibrium. The Einstein relation, which overcomes the effect of statistical noise observed in the theoretically equivalent GreenKubo formulation, was employed. The Einstein relation employed in this work is described by Smith and van Gunsteren.52 However, unlike the original work, we used all of the elements of the stress tensor to improve the convergence and statistics.53 The viscosity coefficient was then given by the expression

η)

1 V 20 kBT

lim tf∞

d dt

[

∑R 〈∆PTRR(t)2〉 + 2 R>β ∑ 〈∆PRβ(t)2〉]

(8)

Here, R and β are indices running over the three Cartesian coordinates, V is the volume (or the average volume in NpT simulations), T is the temperature, and ∆PRβ(t) denotes the displacement of the elements of the pressure tensor PRβ according to

AUA Model of Methylbenzenes

∆PRβ(t) )

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15945

∫0t 21 (PRβ(τ) + PβR(τ))dτ

(9)

and T (t) ) ∆PRβ

∫0t

(

PRR(τ) -

1 3

)

∑β Pββ(τ) dτ

(10)

The microscopic expression for the elements of the pressure tensor PRβ appearing in the integrand of eqs 9 and 10 is given by

PRβ(t) )

1 V

(

∑i

pRi(t)pβi(t) mi

+

∑ ∑ fRij(t)rβij(t) i Dm-xylene > Do-xylene

Figure 4. Variation of the self-diffusion coefficient of xylene isomers with temperature at 0.1 Mpa (a) and with the pressure at 323 K (b). Simulation results computed with the ch-AUA model are compared with the results of Rousseau et al.64 using the OPLS model.

in almost all of the thermodynamic conditions studied. Rousseau et al. have done an extensive analysis about the competition between the “hydrodynamic shape” of the isomers and the electrostatic forces to explain the complex dynamic behavior observed between the isomers.64 Since the diffusion coefficient

15948 J. Phys. Chem. C, Vol. 111, No. 43, 2007

Nieto-Draghi et al.

Figure 5. Variation of the shear viscosity of toluene with the temperature at 0.1 MPa (a) and with pressure at 333 K (b) computed with the ch-AUA model. The experimental data are taken from the work of Et-Tahir,66 Krall et al.,67 and the reference handbook for physical data.68

and shear viscosity obtained with our model and those obtained with OPLS are very close, the general conclusion raised from that work is also valid for here. It is important to remark that our model is, by far, less time consuming than the OPLS force field (which possesses 12 point charges plus 12 Lennard-Jones interaction sites). This aspect should be taken into account in order to evaluate the global performance on reproducing dynamic properties of our ch-AUA model. 4.3. Shear Viscosity Coefficient. We have analyzed the effect of electrostatic interactions in the computation of the shear viscosity reproduced by the new AUA model, including the electrostatic interactions. For the case of toluene, we have compared the temperature variation of the shear viscosity obtained in our simulations with that of Fioroni et al.65 employing an all atom, full flexible model developed using ab initio calculations in Figure 5a. Additionally, we have compared our results with the experimental data of Et-Tahir,66 Krall et al.,67 and the handbook of the physical properties of fluids.68 The agreement of the results obtained with the ch-AUA model with the different sources of experiments is good; we observe an AAD of about 9%, which is considerably less than the 29% of AAD obtained with the model of Fioroni et al. The activation energy of the viscous process reported by our model is also coherent with the experimental values, as can be seen in Table 2. In addition, we observe a general good agreement between the activation energies of the self-diffusion coefficient and shear viscosity coefficients at ambient pressure for the ch-AUA model. The knowledge of the physical properties of hydrocarbons under fuel injection conditions in the new generation of flexifuel engines is a good example of an industrial application that can be covered by molecular simulations. This new generation of injectors may reach pressures of about 250 MPa2, making the acquisition of experimental data difficult. In view of this, we have tested the viscosity reproduced by the our model under the effect of pressure from 0.1 up to 500 MPa at 333 K in Figure 5b. We compare our simulation results with the experimental data of Et-Tahir66 up to 100 MPa. In this region, our simulation shows an AAD of about 9%, which is a good performance. Over 100 MPa, our simulations are pure predictions of the behavior of the shear viscosity of toluene at extreme conditions (points at 250 and 500 MPa). The ability of a particular model to predict the behavior of the shear viscosity is essential from an industrial point of view since the correlative equations such as those of Lohrenz-Bray-Clark69 or Assael et al.70 have poor reliability in predicting the viscosity of fluids out of their range of calibration. To conclude the analysis of the shear viscosity, we show in Figure 6a and b the variation of the viscosity of xylene isomers with the temperature at 0.1 MPa and with the pressure at 323 K, respectively. We have compared our simulation results of

Figure 6. Variation of the shear viscosity coefficient of xylene isomers with temperature at 0.1 Mpa (a) and with the pressure at 323 K (b). Simulation results computed with the ch-AUA and AUA models are compared with the results of Rousseau et al.64 using the OPLS model. The experimental data are taken from Et-Tahir66 and Cadwell et al.71

the ch-AUA model with those obtained with the noncharged version of the AUA potential,24 the simulation results from the work of Rousseau et al.64 employing the OPLS model, and the experimental data of Et-Tahir66 and Cadwell et al.71,72 In Figure 6a, we can see how the OPLS model is the one that better reproduces the shear viscosity for the three isomers, followed by the ch-AUA model and the noncharged version of the AUA model. The differences between models are more pronounced for the case of o-xylene, where the OPLS presents an AAD of about 3.7% compared with the 7.3% of the charged AUA and the 13.9% for the AUA model. For the other two isomers, the new potential is better than or equal to the OPLS model, and the AAD of the noncharged model is superior to 12%. In general, our simulations respect the order in the magnitude of the shear viscosity between isomers when the temperature or the pressure is changed, that is, ηo-xylene > ηp-xylene > ηm-xylene, as it was previously confirmed by the work of Rousseau et al.64 for the case of the OPLS model. This order in the magnitude of the shear viscosity becomes more important when the pressure is increased above 100 MPa, as can be seen in Figure 6b. In

AUA Model of Methylbenzenes

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15949

Figure 7. Variation of the heat flux computed with the PeX algorithm for benzene and methylbenzene molecules (a) at ambient conditions. (b) Linear temperature gradient along the z direction of the space obtained at steady state for each system.

this pressure regime, we are able to compare our simulation results with the experimental data of Cadwell et al. for m-xylene up to 160 MPa. Here, we can observe how the simulations follow the experimental trend, and consequently, we are confident about the predictions of the viscosity of m-xylene at least at 250 MPa, with an error of about 10%. Unfortunately, we were unable to find any other source of experiments to compare with the other two isomers. Finally, according to our previous work about the development of a new AUA potential of benzene,30 we have demonstrated that the viscosities reproduced by the nonelectrostatic models of aromatics are considerably lower than experimental data (here, we obtain more than 12% of AAD for the three isomers using the noncharged AUA model). The reason exposed is still valid for the present work; only the combination of the dispersion repulsion and the electrostatic interactions is able to generate the appropriate orientation between molecules and to produce enough friction to increase the viscosity of an aromatic system. 4.4. Thermal Conductivity. As we have described in the methodological section, we have computed the thermal conductivity of the different methylbenzene molecules using the nonequilibrium method of momentum exchange (PeX).55 We can observe the time evolution of the heat flux computed for the charged AUA model for benzene, toluene, and the three xylene isomers in Figure 7a. The resulting linear temperature profile obtained for each molecule, as a consequence of the imposed heat flux, can be seen in Figure 7b. In general, the NEMD simulations generate a steady-state heat flux after a period of 200 ps for the case of the xylene isomers and after 150 ps for the case of benzene and toluene molecules. The temperature gradient of the four computed molecules with the ch-AUA model are very close; consequently, the thermal conductivities of these molecules are slightly different. The corresponding thermal conductivities obtained applying Fourier’s Law (eq 15) can be compared with the experimental data of Watanabe et al.73 at 298 K and 0.1 MPa in Table 3. In our previous work where the ch-AUA of benzene was developed,30 we did not compute the thermal conductivity of benzene. We have decided to include this molecule in our present analysis. In addition, we have also computed the thermal conductivity using the charged TraPPE model19 and the noncharged TraPPE model12 of benzene and the simulation results taken form the work of Zhang et al.74 employing the all atom OPLS model of benzene. The analysis of this coefficient for the case of benzene reveals that both AUA models present a good agreement with respect to the experimental data, with an AAD of 5.5 and 3.4% for the ch-AUA and the noncharged AUA model, respectively. These results are considerably better than those obtained with the

TABLE 3: Comparison of the Thermal Conductivity λ Computed for the Different Benzene and Methylbenzene Molecules at 298 K and 0.1MPa; The Result for OPLS Model of Benzene is taken from the Work of Zhang et al.,74 and the Percentage of Deviation with Respect to the Experimental Data is in the Last Column molecule

model

benzene

ch-AUA AUA ch-TraPPE TraPPE OPLS exp.a ch-AUA exp.a ch-AUA exp.a ch-AUA exp.a ch-AUA exp.a

toluene o-xylene m-xylene p-xylene a

heat flux 〈Jq〉z (W)

〈∇T〉z (K/Å)

λ (W‚m-1‚K-1)

dev. (%)

0.301 0.300 0.294 0.290

2.258 2.211 2.417 2.650

-5.5 -3.4 -14.1 -22.6 40.0

0.265

2.031

0.248

1.924

0.243

1.938

0.239

2.041

0.133 ( 0.002 0.136 ( 0.001 0.121 ( 0.002 0.109 ( 0.001 0.197 0.14081 0.131 ( 0.001 0.13051 0.129 ( 0.002 0.12875 0.125 ( 0.002 0.12943 0.117 ( 0.002 0.12679

0.07 0.2 -3.1 -7.7

Experimental data from Watanabe et al.73

charged and noncharged version of the TraPPE model (14.1 and 22.6% of AAD, respectively). Surprisingly, as it was remarked by Zhang et al., the OPLS model of benzene presents a poor performance, with 40% deviation. This comparison allows us to extract some conclusions; for instance, the noncharged version of the AUA model presents an AAD which is smaller than the charged AUA model. The same behavior is observed for the case of both versions of the TraPPE model. This means that the electrostatic interactions are not required to reproduce this property for the case of benzene, contrary to what is observed for the case of the shear viscosity in ref 30. This fact can be the result of the cancellation of errors, but a deeper analysis, out of the scope of this work, is required to elucidate this observation. For the rest of the compounds, we have only computed the thermal conductivity using the new ch-AUA potential, and the results obtained are in good agreement with the experimental data with a global AAD that is less than 10%. If we compare the results obtained, we can see that the accuracy of the models decreases following the reduction in the magnitude of the dipole moment of the molecules, that is, toluene > o-xylene > m-xylene > p-xylene. This fact can be the consequence of the influence of dipolar contributions in the way that heat flux is delivered inside of each system. Finally, we can say that the AUA model is better than the other models to reproduce this property for the molecules analyzed and at ambient conditions.

15950 J. Phys. Chem. C, Vol. 111, No. 43, 2007 5. Conclusions The main goal of this work is to show the performance of a new anisotropic united atom intermolecular potential of methylbenzenes that accounts for the quadrupolar and dipolar interactions of these molecules. The new potential was adjusted using experimental thermodynamic data at saturation in an accompanying work. During the process of optimization, the viscosity of the model was verified at ambient conditions to ensure coherence with dynamic properties. We have compared our potential with other potentials that include or exclude electrostatic interactions. We found the same behavior as that observed in our previous work about the development of a new AUA potential of benzene30 where we demonstrated that the viscosities reproduced by the nonelectrostatic models of aromatics were considerably lower than the experimental data (here, we obtain more than 10% of AAD for the three isomers using the noncharged AUA model). The analysis of the reorientational relaxation times of toluene showed that the ch-AUA model is in good agreement with the experimental data; in particular, our simulations reproduced the degree of anisotropy and the order between the rotational diffusion constant around the different principle axes of rotation and the relaxation time of the CH vectors at different positions (ortho, meta, and para) with respect to the methyl group. Due to the lack of experimental data, we were unable to verify the rotational diffusion constant of the xylene isomers; however, our simulation results permitted us to analyze the influence of the dipole moment and the moments of inertia on the behavior of the reorientational dynamics of these isomers. The experimental self-diffusion coefficient of toluene at ambient pressure is well reproduced with approximately 4% of AAD for the ch-AUA model. In addition, our simulations respect the activation energy of the diffusive process with an AAD of about 1.6%. Again, the lack of experimental data for the xylene isomers prevents us from confirming the accuracy of our model to reproduce this data; however, we observe a good agreement when our results are compared with those reported by the work of Rousseau et al.64 employing the OPLS model. We have analyzed the behavior of the different models to reproduce the shear viscosity at different conditions of pressure and temperature. In this case, the ch-AUA model of toluene presents an AAD of 9%, which is an excellent performance in comparison with the all atom model of Fioroni et al.65 (29% of AAD). The behavior obtained for the case of the xylene isomers is in good agreement with the experimental data with a global AAD of less than 8%. The comparison with the noncharged version of the AUA model, with an AAD of about 14%, confirms our analysis that this property requires the explicit inclusion of the electrostatic interactions to well reproduce the increment of viscosity when the pressure of the system is increased. The relative orientation between molecules and the correct description of the friction at high densities are only achieved by the models when they include electrostatic interactions to mimic the quadrupolar and dipolar moment of these molecules. The accuracy of the potential allows us to extend our simulations to predict the behavior of the shear viscosity up to 500 MPa; such conditions are out of the range of capability of the experimental devices currently available. Finally, we have computed the thermal conductivity of benzene at ambient conditions obtained with the new model and compared these results with the available experimental data and other intermolecular potentials. We observe that the simulation results using the AUA representation are better than the other intermolecular potentials (between 3.4 and 5.5% of

Nieto-Draghi et al. AAD for AUA, between 14 and 22.6% of AAD for TraPPE model, and about 40% for OPLS). For toluene, o-xylene, m-xylene, and p-xylene molecules, we have obtained an AAD of less than 8%. What is also important to remark, and one of the principal achievements of this work, is the fact that a simple model of methylbenzenes can be efficiently parametrized to reproduce several dynamic properties at the same time. Moreover, the combination of the anisotropic united atom approach and the inclusion of the electrostatic interaction allows our model to correctly reproduce the temperature and pressure dependence of almost all dynamic processes analyzed. We have shown that our new model, though having only three point charges, is able to compete in accuracy with all atoms models but with four times less computational effort. Even with the increment in computer power, computation of properties such as the shear viscosity of real fuels at extreme conditions would remain unaffordable for long simulations, if all atoms models were employed. In these particular cases, a simplified approach like the one presented in this work is a good alternative. Acknowledgment. The authors would like to thank Bernard Rousseau for the use of the Newton code. Supporting Information Available: Tables containing the data used in all figures. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Walter, B.; Gatellier, B. Oil Gas Sci. Technol. 2003, 58, 101. (2) Delacourt, E.; Demet, B.; Besson, B. Fuel 2005, 84, 859. (3) Jorgensen, W. L.; Madura, J. D. J. Am. Chem. Soc. 1984, 106, 6638. (4) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. (5) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (6) Smit, B.; Karaborni, S.; Siepmann, J. I. J. Chem. Phys. 1995, 102, 2126. (7) Sun, H. J. Phys. Chem. B 1998, 102, 7338. (8) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569. (9) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4508. (10) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 5370. (11) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 3093. (12) Wick, C. D.; Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 8008. (13) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. J. Chem. Phys. 1998, 108, 9905. (14) Nath, S. K.; Banaszak, B. J.; de Pablo, J. J. J. Chem. Phys. 2001, 114, 3612. (15) Nath, S. K.; de Pablo, J. J. J. Mol. Phys. 2000, 98, 231. (16) Errington, J. R.; Panagiotopoulos, A. Z. J. Phys. Chem. B 1999, 103, 6314. (17) Errington, J. R.; Panagiotopoulos, A. Z. J. Chem. Phys. 1999, 111, 9731. (18) Contreras-Camacho, R. O.; Ungerer, P.; Boutin, A.; Mackie, A. D. J Phys. Chem. B 2004, 108, 14109. (19) Wick, C. D.; Siepmann, J. I.; Koltz, W. L.; Schure, M. R. J. Chromatogr., A 2002, 954, 181. (20) Cacelli, I.; Cinacchi, G.; Prampolini, G.; Tani, A. J. Am. Chem. Soc. 2004, 126, 14278. (21) Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Fuchs, A. H. J. Chem. Phys. 2000, 112, 5499. (22) Bourasseau, E.; Ungerer, P.; Boutin, A.; Fuchs, A. H. Mol. Simul. 2002, 28, 317. (23) Bourasseau, E.; Haboudou, M.; Boutin, A.; Fuchs, A. H.; Ungerer, P. J. Chem. Phys. 2003, 118, 3020. (24) Ahunbay, M. G.; Perez-Pellitero, J.; Contreras-Camacho, R. O.; Teuler, J. M.; Ungerer, P.; Mackie, A. D.; Lachet, V. J. Phys. Chem. B 2005, 109, 2970. (25) Contreras-Camacho, R. O.; Ungerer, P.; Ahunbay, M. G.; Lachet, V.; Perez-Pellitero, J.; Mackie, A. D. J. Phys. Chem. B 2004, 108, 14115.

AUA Model of Methylbenzenes (26) Delhommelle, J.; Tschirwitz, C.; Ungerer, P.; Granucci, G.; Millie, P.; Pattou, D.; Fuchs, A. H. J. Phys. Chem. B 2000, 104, 4745. (27) Perez-Pellitero, J.; Ungerer, P.; Mackie, A. D. J. Phys. Chem. B 2006, 111, 4460. (28) Smith, G. D.; Jaffe, R. L. J. Phys. Chem. 1996, 100, 9624. (29) Wick, C. D.; Martin, M. G.; Siepmann, J. I.; Schure, M. R. Int. J. Thermophys. 2001, 22, 111. (30) Nieto-Draghi, C.; Ungerer, P. J. Phys. Chem. B 2007, accepted. (31) Eike, D. M.; Maginn, E. J. Chem. Phys. 2006, 124, 164503. (32) Zhao, X.; Chen, B.; Karaborni, S.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 5368. (33) Yonezawa, Y.; Nakata, K.; Takada, T.; Nakamura, H. Chem. Phys. Lett. 2006, 428, 73. (34) Laaksonen, A.; Stilbs, P.; Wasylishen, R. J. Chem. Phys. 1998, 108, 455. (35) Furukata, S.; Ikawa, S. J. Chem. Phys. 1998, 108, 5159. (36) Furukata, S.; Ikawa, S. J. Chem. Phys. 2000, 113, 1942. (37) Nieto-Draghi, C.; Avalos, J. B.; Contreras, O.; Ungerer, P.; Ridard, J. J. Chem. Phys. 2004, 121, 10566. (38) Hahn, J. R.; Jeong, H.; Jeong, S. J. Chem. Phys. 2005, 123, 244702. (39) Lachet, V.; Buttefey, S.; Boutin, A.; Fuchs, A. Phys. Chem. Chem. Phys. 2001, 3, 80. (40) Su, B. L.; Norberg, V.; Hansenne, C.; Mallmann, A. D. Adsorption 2000, 6, 61. (41) Bourasseau, E.; Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Rousseau, B.; Fuchs, A. J. Chem. Phys. 2000, 112, 5499. (42) DIPPR, Thermophysical Property Database. http://dippr.byu.edu/ database.asp. Thermophysical Properties Laboratory, version 2005. (43) Dennos, G. R.; Ritchie, G. L. D. J. Phys. Chem. 1991, 95, 656. (44) Vbrancich, J.; Ritchie, G. L. D. J. Chem. Soc., Faraday Trans. 2 1980, 76, 648. (45) Vbrancich, J.; Ritchie, G. L. D. Chem. Phys. Lett. 1983, 94, 63. (46) Sturz, L.; Do¨le, A. J. Phys. Chem. A 2001, 105, 5055. (47) Witt, R.; Sturz, L.; Do¨lle, A.; Mu¨ller-Plathe, F. J. Phys. Chem. B 2000, 104, 5716. (48) Kohlrausch, F. Poggendorff’s Ann. Phys. 1863, 119, 337. (49) Williams, G.; Watts, D. C. Trans. Faraday Soc. 1970, 66, 80. (50) Liu, G.; Mackowiak, M.; Li, Y.; Jonas, J. J. Chem. Phys. 1991, 94, 239.

J. Phys. Chem. C, Vol. 111, No. 43, 2007 15951 (51) Plomb, M.; Schreurs, M.; Bulthuis, J. J. Chem. Phys. 1988, 88, 5202. (52) Smith, P. E.; van Gunsteren, W. F. Chem. Phys. Lett. 1993, 215, 315. (53) Mondello, M.; Grest, G. S. J. Chem. Phys. 1997, 106, 9327. (54) Nieto-Draghi, C.; Avalos, J. B. Mol. Phys. 2003, 101, 2303. (55) Mu¨ller-Plathe, F. J. Chem. Phys. 1997, 106, 6082. (56) Frenkel, D.; Smit, B. Understanding Molecular Simulations: From Algorithms to Applications, 2nd ed.; Academic Press: Amsterdam: The Netherlands, 2001. (57) Evans, D. J. Mol. Phys. 1997, 34, 317. (58) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (59) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, 2nd ed.; Oxford University Press: Oxford, U.K., 1987. (60) Kru¨ger, G. J.; Weiss, R. Z. Natursforsch., A: Phys. Sci. 1970, 25, 777. (61) Trepa˜dus, V.; Raˆpeanu, S.; Pa˜dureanu, I.; Parfenov, V. A.; Navikov, A. G. J. Chem. Phys. 1974, 60, 2832. (62) Winfield, D. J. J. Chem. Phys. 1971, 54, 3643. (63) Pickup, S.; Blum, F. Macromolecules 1989, 22, 3961. (64) Rousseau, B.; Petravic, J. J. Phys. Chem. B 2002, 106, 13010. (65) Fioroni, M.; Vogt, D. J. Phys. Chem. B 2004, 108, 11774. (66) Et-Tahir, A. Ph.D. Thesis, University of Pau, 1993. (67) Krall, A.; Sengers, J. V. J. Chem. Eng. Data 1992, 37, 349. (68) Lide, D. R. Handbook of Chemistry and Physics, 73rd ed.; CRC Press: London, 1993. (69) Lohrenz, J.; Bray, B. G.; Clark, C. R. J. Pet. Technol. 1964, 16, 1171. (70) Assael, J. M.; Papadaki, M.; Wakeham, W. A. Int. J. Thermophys. 1991, 12, 449. (71) Caudwell, D.; Goodwin, A. R. H.; Trusler, J. P. M. J. Pet. Sci. Eng. 2004, 44, 333. (72) Assael, J. M.; Papadaki, M.; Wakeham, W. A. Int. J. Thermophys. 1991, 12, 449. (73) Watanabe, H.; Kato, H. J. Chem. Eng. Data 2004, 49, 809. (74) Zhang, M.; Lussetti, E.; de Souza, L. E. S.; Mu¨ller-Plathe, F. J. Phys. Chem. B 2005, 109, 15060.