Anisotropic United-Atoms (AUA) Potential for Alcohols - The Journal of

J. Phys. Chem. B 2008, 112, 9853–9863

9853

Anisotropic United-Atoms (AUA) Potential for Alcohols Javier Pe´rez-Pellitero,*,† Emeric Bourasseau,‡ Isabelle Demachy,§ Jacqueline Ridard,§ Philippe Ungerer,| and Allan D. Mackie† Departament d’Enginyeria Quı´mica, ETSEQ, UniVersitat RoVira i Virgili, AVinguda dels Paı¨sos Catalans, 26 Campus Sescelades, 43007 Tarragona, Spain, De´partement de Physique Theórique et Appliqueé, CEA, BP 12, F-91680 Bruyeres Le Chatel, France, Laboratoire de Chimie-Physique, UniVersite´ Paris-Sud 11, 91405 Orsay, France, and Institut Français du Pe´trole, 1-4 AVeneda de Bois Preau, 92852 Rueil Malmaison, Cedex, France ReceiVed: March 16, 2008; ReVised Manuscript ReceiVed: May 10, 2008

An anisotropic united-atom (AUA4) intermolecular potential has been derived for the family of alkanols by first optimizing a set of charges to reproduce the electrostatic potential of the isolated molecules of methanol and ethanol and then by adjusting the parameters of the OH group to fit selected equilibrium properties. In particular, the proposed potential includes additional extra-atomic charges in order to improve the matching to the electrostatic field. Gibbs ensemble Monte Carlo simulations were performed to determine the phase equilibria, while the critical region was explored by means of grand canonical Monte Carlo simulations combined with histogram reweighting techniques. In order to increase the transferability of the model, only the parameters of the Lennard-Jones OH group have been fitted, the parameters of the other AUA groups are taken from previous works. Nevertheless, a good level of agreement was obtained for all compounds considered in this work. In particular, excellent results were obtained for the Henry constants calculation of different gases in alkanols. 1. Introduction The increasing popularity of the use of oxygenated compounds in the field of biofuels has given rise to increased interest in the family of alcohols in the petroleum industry. In particular, the simplest and most commonly used alcohols in the elaboration of biodiesels are methanol (common name methyl alcohol) and ethanol (ethyl alcohol). Alcohols are polar molecules due to the presence of the hydroxyl group in the “head” of the alkyl flexible chain which gives the molecule a nonpolar character. The strong polarity of this group is responsible for the high critical point of short alkanols like methanol and ethanol when compared with hydrocarbons of similar molecular weight. Since alcohols are also widely used as solvents, the thermodynamic properties of this family and its mixtures has been widely investigated experimentally, theoretically, or by means of molecular simulation.1,2 Molecular simulation is rapidly becoming an alternative way to predict equilibrium thermodynamic properties of pure fluids and mixtures. Of course, the accuracy of the predictions is intimately related with the intermolecular potentials used to reproduce the interactions between molecules. In the past decade, potential models, also referred to as force fields, have been developed for a wide range of compounds.3–5 One of the most common approaches for modeling hydrocarbons and other flexible molecules is the use of the united-atoms model scheme, where each chemical group is represented by one Lennard-Jones center, avoiding the explicit representation of the hydrogens in the different groups. This scheme results in a significant reduction of the computational time as compared to all-atoms * To whom correspondence should be addressed. † Universitat Rovira i Virgili. ‡ CEA. § Universite ´ Paris-Sud 11. | Institut Franc ¸ ais du Pe´trole.

models since the number of pair interactions goes as the square of the number of sites. Improvements on the standard unitedatoms model, where typically a 6-12 Lennard-Jones center of force is placed on top of the most significant atom, have been proposed. For instance, Errington et al. used a Buckingham exp-6 potential,6 where three parameters are involved in each center of force and the repulsive part of the LJ interaction is replaced by the more correct exponential term. Chen et al. introduced additional sites on the centers of the C-H bonds.7 Also, Ungerer et al.8 obtained the AUA4 model by reparameterizing the initial anisotropic united-atom (AUA) model proposed by Toxvaerd.9,10 The AUA model consists of a displacement of the Lennard-Jones centers of force toward the hydrogen atoms, converting the distance of displacement into a third adjustable parameter. In the case of the shorter alcohols, as in the case of methanol, several models have been proposed to investigate its behavior11–13 mainly at ambient temperature and density. The OPLS (optimized potential for liquid simulation) force field is found to be unreliable at high temperatures and longer chains. Vapor-liquid coexistence curves from methanol to hexanol have been obtained14 by combining the charges of the OPLS force field with the Lennard-Jones (LJ) hydroxyl group parameters derived by Van Leeuwen in ref 15 and the alkyl groups of the Siepmann-Karaborni-Smit (SKS)16 force field. Although the model of Van Leeuwen behaves much better, it requires the methyl groups of the different alcohols to be specifically readjusted. More recently, in 2001 Chen et al.2 extended the TraPPE5 (transferable potentials for phase equilibria) force field to several primary, secondary, and tertiary alcohols trying to introduce as few as possible new “pseudo-atoms” to account for the difference in electronegativity between an alkane C-C bond or a C-O bond. In general, they obtained a good agreement for a variety of alcohols and were also able to

10.1021/jp802282p CCC: $40.75  2008 American Chemical Society Published on Web 07/23/2008

9854 J. Phys. Chem. B, Vol. 112, No. 32, 2008

Pe´rez-Pellitero et al.

TABLE 1: Lennard-Jones Parameters and Distances Between the Major Atom and the AUA Force Center group

σ (Å)

ε/k (K)

δ (Å)

CH38 CH28 CHar21 Car20 OH

3.607 3.461 3.246 3.246 3.034

120.15 86.29 89.42 37.73 85.27

0.216 0.384 0.407

TABLE 3: Bending Parameters

length (Å)

C-C Car-Car (aromatic ring) C-O O-H

1.53 1.39 1.425 0.96

angle

Σ (deg)

Car-Car-Car Car-CH2-O CHx-O-H

112.7 109.47 108.9

2. Model 2.1. Intermolecular Interactions. To describe the dispersion interactions, the different molecules are represented by a set of interacting Lennard-Jones sites for each CH3, CH2, C, aromatic C and CH, or OH and the respective Coulombic interactions. The energy between two different particles, i and j, is calculated according to

[( ) ( ) ] [ ] σij rij

12

-

σij rij

6

+

qiqj 4πε0rij

(1)

where rij is the separation between a given pair of particles i and j, εij is the LJ well depth, σij is the LJ size of particles, qi and qj are the partial charges, and ε0 is the permittivity of space. To calculate the parameters between unlike united atoms, we use the Lorentz-Berthelot combining rules

εij ) √εiiεjj σij )

σii + σjj 2

k0/k (K)

CHx-CH2-O CHx-CH2-CHy CHx-CH-CHy CHx-C-CHy

109.47 114 112 109.47

59 800 74 900 72 700 70 311

torsion

reproduce correctly the azeotropic compositions of binary n-hexane/methanol mixtures. Finally, solubility data obtained by simulation can be found for the TraPPE-UA ethanol model17,18 as well as in a preliminary study of the solubility of different gases in ethanol using the AUA4 intermolecular potential of this work.19 In this work, we use an optimized intermolecular potential for the hydroxyl group of the family of alcohols aimed at giving a quantitative description of both liquid and coexistence properties based on an extension of the AUA4 intermolecular potentials already given in previous works.8,20 In order to increase the transferability of the AUA4 model, the same set of LJ parameters for the alkyl chain previously optimized for alkanes has been used without defining new pseudo-atoms. The coulombic charges have been parameterized to reproduce the electronic structure of ethanol and methanol obtained from ab initio calculations. In a second step, the LJ parameters of the OH group have been fitted to reproduce selected equilibrium properties of methanol and ethanol, where the influence of the hydroxyl group is higher.

ULJ(rij) ) 4εij

θ0 (deg)

TABLE 4: Torsion Potential Constants

0.0951

TABLE 2: Bond Lengths and Angles bond

angle

(2)

CHx-CH2-CH2-CHy CHx-CH2-CH2-OH CHx-CH2-O-H CHx-Car-O-H

c0/k (K)

c1/k (K)

c2/k (K)

c3/k (K)

0 0 0 0

670.06 353.24 419.64 0

-136.38 -106.68 -58.34 0

1582.64 1539.86 375.86 327.12

parameters for the different potentials are shown in Table 1. In Table 2 we show the bond distances as well as the bond angles of the aromatic ring. The COH angle and lengths of the bonds C-O and O-H have been determined from ab initio calculations. The angle bending is controlled by means of the following expression

Ubend 1 ) k0[cos(θ) - cos(θ0)]2 k 2

where k0 is the bending constant and θ and θ0 are the bending and equilibrium angles, respectively. While the COH angle has been considered to be rigid, as will be explained more in detail in the following section, the values for the CCO bending angles have been taken as in the case of Chen et al. from the AMBER94 force field.22 The different values of the bending parameters are listed in Table 3. The torsional potentials are taken from the OPLS-UA force field3,23

1 1 Utors(φ) ) c0 + c1(1 + cos φ) + c2(1 - cos 2φ) + 2 2 1 ( c 1 + cos 3φ) (5) 2 3 where φ is the dihedral angle and the Fourier coefficients are listed in Table 4. 2.2. Charges Optimization. Electrostatic charges have been determined by a least-squares fit of the ab initio electrostatic potential of the isolated molecule. The fit has been performed on a grid of points located on a single shell of fused spheres centered on the atoms with radii equal to two times the van der Waals radius of the corresponding atom. The points have been generated according to the methodology proposed by Singh et al.24 The geometries have been optimized at the B3LYP/6311+G(3d,3p) level, and the electrostatic potential have then been calculated at the MP2/6-311+G(3d,3p) level, using the GAUSSIAN package.25 Once the values have been determined, we can calculate the RRMS (relative root mean square) deviation of the charge distribution which will give us information about the accuracy of the optimized set of charges

∑ ( ∑ ( N

(3)

The different alkanols considered in this work have been modeled using the AUA4 for the CH3 and CH2 groups from Ungerer et al.8 In the case of the aromatic ring of phenol we also made use of the CH intermolecular potential for polyaromatic compounds by Contreras et al.21 and the aromatic C group by Ahunbay et al.20 The hydroxyl group has then been fitted to reproduce experimental properties of methanol and ethanol. The

(4)

2 Vcalcd - Vref p p )

RRMS )

p)1

(6)

N

2 Vref p )

p)1

The lower the RRMS deviation, the higher accuracy obtained in the representation of the electrostatic potential over the surface.

Anisotropic United-Atoms Potential for Alcohols This method presents the inconvenience of optimizing the charges only with respect to a single conformation of the molecule, so we are assuming the charge distribution to be valid for the different conformations adopted by the flexible molecule. The second assumption has to do with the fact that the charges have been determined for an isolated molecule, so we are assuming that the electrostatic potential created by a molecule in a vapor or liquid phase is equivalent to that created by an isolated molecule. This approximation is equivalent to neglecting the polarization term created by the surrounding molecules. In this work, we want to obtain a single optimized set of charges which can be extended to a wide variety of alcohols, either rigid or flexible. The described methodology has been applied to different conformations of methanol and ethanol in order to find the distribution of charges which represents more accurately these two alcohols, which are considered representative of the further molecules to be studied. This assumption is made based on the fact that these two alcohols have the highest share of the electrostatic energy in the total energy, so that if the set of charges is able to reproduce correctly these target molecules, it can be also expected to reproduce correctly heavier alcohols. In Figure 1 we show the two different distributions of charges that have been investigated: the first one, referred to from now on as distribution A, has three charges, two of them placed over the hydrogen and the neighboring carbon atom and the third one on the bisector of the COH angle. The second distribution, referred to as distribution B, has four different charges, two of them placed over the hydrogen and neighboring carbon atom while the third one is placed on the bond between the carbon and the oxygen and the last one on the bisector of the COHˆ angle. The first distribution is analogous to the one commonly used for water,26 so it is interesting to test the behavior of this type of distribution for the alcohols. Optimization of both distributions has been made according to three different configurations: methanol, ethanol-trans, and ethanol-gauche where the dihedral angles are 180° and 60°, respectively. These last two configurations correspond to those found to be more stable experimentally. In Table 5 we can see the values of the RRMS deviations obtained for the six charge distributions. We can see the important effect of adding the fourth charge on the bond between the carbon and the oxygen. The average deviation obtained for distribution A is around 25%, while for distribution B this value is only about 12%. Although the three optimizations for distribution B present very similar values, we have selected the set of charges optimized for the ethanol-gauche since it provides more constant RRMS values for the three different considered conformations. The values of the charges obtained in this case, which are shown in Table 6, provide a dipolar moment of 1.71 D for the different compounds (the same geometry was used) which compares well with the experimental values 1.70 and 1.69 D of methanol and ethanol, respectively.27 3. Simulation Details To determine the coexistence properties, we used the Gibbs ensemble Monte Carlo (GEMC) method combined with a configurational bias scheme4 and an additional bias for insertion of the center of mass of the molecules.28 The selected probabilities for the different moves were generally set to 0.05 for rigid-body translation, 0.05 for rigid-body rotation, 0.395 for partial regrowths, 0.5 for transfer, and 0.005 for volume change. A total of 250 molecules have been used in all simulations using a cutoff radius equal to one-half of the box length with standard finite-size Lennard-Jones corrections,29 while the Ewald sum-

J. Phys. Chem. B, Vol. 112, No. 32, 2008 9855

Figure 1. Charge distributions A (left) and B (right). The distances O-q2 and CHx-q4 are 0.3 and 0.15 Å, respectively.

TABLE 5: RRMS Values in Percent (%) of the Different Charge Distributions distribution A optimized for methanol distribution A optimized for ethanol-trans distribution A optimized for ethanol-gauche distribution B optimized for methanol distribution B optimized for ethanol-trans distribution B optimized for ethanol-gauche

methanol

ethanol-trans

ethanol-gauche

25.9

25

24

26

24

24

27

25

23

10.9

12

13

13

10

14

11

12

12

TABLE 6: Charge Values Obtained for the Selected Optimization

distribution B optimized for ethanol-gauche

charge q1

charge q2

charge q3

charge q4

-1.99

-1.49

0.70

2.78

mation technique was used to account for the long-range electrostatic interactions. Coexistence densities, average pressure (calculated through the virial), and enthalpies of vaporization (computed as the difference in enthalpy between the two simulation boxes) were calculated from the different simulations. To investigate the critical region, we used grand canonical Monte Carlo (GCMC) simulations combined with histogram reweighting.30,31 As for the Gibbs ensemble, a configurational bias scheme4 and an additional bias for insertion of the center of mass of the molecules28 were implemented to enhance the acceptance rate of insertions. The probabilities of the different moves used in the simulations were 0.05 for translations, 0.05 for rotations, 0.3 for partial regrowths, and 0.6 for insertion or deletion. A simulation box of size L ) 32 Å was used in all cases. Systems were allowed to equilibrate for at least one million MC steps followed by the production run. In the case of a liquid phase, the length of the production run was set normally to at least 25 million MC steps. 4. Determination of the Lennard-Jones Parameters In order to optimize the Lennard-Jones parameters for the OH group, we selected ethanol and methanol experimental results as the reference data. These two compounds have been selected because of being the molecules where the effect of the hydroxyl group is expected to be highest. In addition, including both compounds in the optimization method takes into account the effect of both CH3 and CH2 as adjacent groups. To optimize the Lennard-Jones parameters, a GEMC simulation was carried



out at a reduced temperature close to the critical temperature (450 K) for each of the compounds as well as an NPT simulation at a temperature near to the normal boiling point (300 K). This set of four simulations allows to optimize the parameters with reference to a total of 10 different properties of ethanol and methanol. Once the simulations have been carried out, an optimization procedure has been applied until the error has been minimized. The optimization methodology has already been described in detail32 and is thus only briefly described here The normalized error criterion selected to drive the optimization procedure is given by the following expression mod

1 (X i n i)1 n

F)

∑

- X iexp)2 si2

(7)

where n is the total number of reference data, Xmod is the value i of the ith property calculated in the simulations, Xexp is the i experimental value of the ith property, and si is the statistical uncertainty of the computed variable Xmod (estimated from the i standard block averaging technique29). F is a function of the different yi Lennard-Jones parameters to be optimized, in this case y1 ) εS and y2 ) σS. To perform this optimization we introduced in the F function vaporization enthalpies, liquid densities, and saturation pressures of ethanol and methanol at two different temperatures, the first one at an intermediate point between the critical point and the normal boiling point and the second one at a temperature below the normal boiling point. In the latter case the saturation pressure has not been used since an NPT simulation was used.

Figure 2. Coexistence curves of ethanol and methanol. Solid and dashed lines represent experimental data37 for methanol and ethanol, respectively. Asterisks denote experimental critical points. Diamond and square symbols represent the simulation results for methanol and ethanol, respectively. Open symbols denote UA model results from ref 2 with GEMC, while filled symbols are AUA4 model + HR results.

5. Calculation of the Critical Point In order to accurately locate the critical point, we applied a recently implemented33,34 methodology based on the calculation of a fourth-order cumulant (Binder parameter) combined with the use of finite size scaling35 techniques. The Binder parameter UL is defined as follows

UL )

〈m4 〉 〈m2 〉2

(8)

where m is an appropriate order parameter, in our case density. The procedure is as follows: system size-dependent critical conditions are determined for different system sizes by locating the temperature at which the system takes on the “universal” Ising value of the system. Once these values are calculated, the critical constants are expected to change as a function of the system size according to known scaling laws.35 For instance, the temperature of the finite system is expected to vary near the critical point with system size as ( ) 〈T 〉c(L) - 〈T 〉c(∞) ≈ L- θ+1 ⁄ν

(9)

where Θ ) 0.54 and ν ) 0.629.7,36 In the same way, the critical density of the finite system is expected to vary with system size as ( ) 〈F 〉c(L) - 〈F 〉c(∞) ≈ L- d-1⁄ν

(10)

where d is the dimensionality of the system. By plotting the finite size critical values according to these functions, an extrapolation can be made to estimate the corresponding infinite system size value.

Figure 3. Coexistence curve of 1-propanol. Solid line represents experimental data37 for 1-propanol. Open symbols denote UA model results from ref 2 with GEMC, while filled symbols are AUA4 model results using GEMC.

6. Performance of the Optimized Parameters The developed force field has been tested for several alcohols including methanol and ethanol, which were involved in the fitting procedure. In Figures 2 and 3, the vapor-liquid coexistence curves are presented for the density of the three shorter alcohols considered: methanol, ethanol, and 1-propanol. Saturated vapor pressures and heats of vaporization are compared to experimental data37 in Figures 4, 5, 6, and 7, respectively. The liquid densities compare well with the experimental results, particularly in the case of methanol. The deviations are slightly higher in the cases of ethanol and 1-propanol. In comparison with the UA model, we can state that our model better reproduces the critical point while the UA model describes better the coexistence curve at low temperatures. In general, it appears that both models are unable to reproduce exactly the shape of the liquid-vapor-phase coexistence curve, failing either in the near-critical region or at low temperatures. Further improvements in the model, such as taking into account the flexibility of the COH angle, parallel optimization of the partial charges

Anisotropic United-Atoms Potential for Alcohols

Figure 4. Saturated vapor pressures of ethanol and methanol. Line styles for experimental data and symbols as in Figure 2.


Figure 7. Heats of vaporization of 1-propanol. Line styles for experimental data and symbols as in Figure 3.

TABLE 7: Critical Properties and Normal Boiling Points of the Different Alcohols molecule methanol ethanol propan-1-ol octanol phenol Figure 5. Heats of vaporization of ethanol and methanol. Line styles for experimental data and symbols as in Figure 2.

type

Tc (K)

this model TraPPE-UA experiment this model TraPPE-UA experiment this model TraPPE-UA experiment this model TraPPE-UA experiment this model experiment

505.6(1) 502(2) 512.5 512.0(5) 514(2) 514 528.5(9) 538(2) 536.8 624.2(5) 629(2) 652.5 642.8(5) 694.2

Fc (kg/m3)

Tb (K)

278(3) 285(4) 273 280 (3) 281(3) 274 275(6) 290(4) 276 252(4) 270(2) 262 342(8) 411

355 340(1) 338 367 353(1) 351 377 368(2) 370.4 468 460(3) 468.4 446.2 455

experimental results. As could be expected, this effect becomes less important when we consider longer chains. Once more, following the behavior observed in the case of the saturated liquid densities, the AUA4 model reproduces better the critical region while it presents deviations at low temperatures. In the case of the UA model, it reproduces better the low temperatures while it presents higher deviations near the critical region. The

Figure 6. Saturated vapor pressures of 1-propanol. Line styles for experimental data and symbols as in Figure 3.

and the LJ parameters,38,39 or considering polarization in the model, would most likely be necessary to reproduce more accurately the shape of the diagram. The saturation pressures of methanol and ethanol are found to be in very good agreement with the experimental results, while 1-propanol is slightly overestimated. In the case of methanol, where the influence of the hydroxyl group is highest, the heats of vaporization present significant deviations from the

Figure 8. Radial distribution functions of methanol at 300 K. Solid and dashed lines represent the results of this model and experimental results from refs 44 and 45, respectively. From left to right, the order of appearance of the different peaks corresponds to gH-H, gO-O, and gCC, respectively.


Figure 9. Binder parameter intersections for different system sizes for methanol, ethanol, and 1-propanol. Inset shows more detail in the case of methanol.

Figure 10. Finite size scaling of the intersections with the universal Ising value of the Binder parameter for methanol.

Figure 11. Coexistence curves of octanol and phenol. Solid and dashed lines represent experimental data37 for octanol and phenol, respectively. Square and circle symbols represent the simulation results for octanol and phenol, respectively. Open symbols denote UA model results from ref 2 with GEMC, while filled symbols are AUA4 model + GEMC results.

fact that both models present the same trend for the vaporization enthalpies reinforces the hypothesis that incorporation of polarization to the model should be considered to better describe this property. For instance, considering the QEE variable charges


Figure 12. Saturated vapor pressures of octanol and phenol. Line styles for experimental data and symbols as in Figure 11.

Figure 13. Vaporization enthalpies of octanol and phenol. Line styles for experimental data and symbols as in Figure 11.

method40–43 may result in a better reproduction of this property without increasing dramatically the number of parameters or the computational time. The critical points of the different alcohols predicted with the AUA4 model are presented in Table 7. The estimations of the UA model are also shown, although a direct comparison is not possible due to the fact that in the work of Chen et al.2 the scaling exponent β has been fitted in order to give a better description of the critical region. In particular, the critical exponent was fitted to experimental data over the same range of temperatures used in the simulations, incorporating in this way temperatures far away from the critical point where the universal Ising exponent may not hold. They obtain in all cases significantly smaller values than the universal Ising value of 0.325. As can be observed in Figure 9, we have not observed for the different alcohols any evidence that this model belongs to a different universality class. The intersections take place close to the universal Ising value of UL ) 1.6035 as in the case of the LJ fluid.33,34 In Figure 10 we show the results obtained for the finite size scaling of the intersections with universal values of the Binder parameter for the case of methanol. Finally, the structure of methanol at 300 K has been investigated through the radial distribution function (RDF) plotted in Figure 8 where a comparison with the experimental results of Yamaguchi et al.44,45 has been established. The results obtained for the hydrogen-hydrogen (H-H), oxygen-oxygen (O-O), and carbon-carbon (C-C) curves compare qualitatively well with the experimental data. The peak positions are found to be at 2.22, 2.6, and 3.9 Å, respectively. In the case of ethanol,



TABLE 8: Intermolecular Potential Parameters of the Five Different Gases Considered position force center/ charge CH451 C N252 N N q1 O253 O1 O2 q1 q2 q3 CO254 C O1 O2 H2S55 H2S q1 H1 H2 S

LJ parameters charge

X 0

Y

Z

ε/k (K)

σ (Å) 3.7327

0

0

149.92

0.549 -0.549 0

0 0 0

0 0 0

36.0 36.0

3.30 3.30

-0.5075 -0.5075 1.015

-0.485 0.485 0 0.2 -0.2

0 0 0 0 0

0 0 0 0 0

43.183 43.183

3.1062 3.1062

0 0 4.2 -2.1 -2.1

0 1.149 -1.149

0 0 0

0 0 0

28.129 80.507 80.507

0 0 1.149 -1.149 -1.149

0 0 0.9639 -0.9639 0

0 0 0.9308 0.9308 0

250

0

2.757 3.033 3.033

0.6512 -0.3256 -0.3256

3.73

0 0.4 0.25 0.25 -0.9

the results obtained are very similar to those obtained using the TraPPE-UA force field2 being the heights of the peaks are slightly higher in the case of the AUA4 model. In Figure 11 the coexistence liquid densities for octanol and phenol are presented. The predictions of the AUA4 model compare well with the available experimental results, particularly at low temperatures. As in the case of the UA model, the critical point of octanol is underpredicted. This behavior is accentuated as compared with the shorter alcohols. This effect can be attributed to the fact that the alcohols involved in the fitting process are the shortest ones: methanol and ethanol. The saturated vapor pressures shown in Figure 12 are in very good agreement for both compounds, presenting a slight overestimation in the case of phenol. Finally, in Figure 13 are shown the predictions obtained for the heats of vaporization. In the case of octanol we obtained excellent agreement, while in the case of phenol the AUA4 model fails to reproduce the experimental trend, overestimating the vaporization enthalpy at low temper-

Figure 15. Henry constants for different gases in ethanol as a function of temperature. Filled and shaded symbols represent results from the simulations using Kong and Lorentz-Berthelot combining rules, respectively. Open symbols represent experimental results.19,65,67,68,70

Figure 16. Henry constants for different gases in 1-propanol as a function of temperature. Filled symbols represent results from the simulations, while open symbols represent experimental results.67–71

Figure 17. Henry constants for different gases in octanol as a function of temperature. Filled symbols represent results from the simulations, while open symbols represent experimental results.66,67,72,73

Figure 14. Henry constants for different gases in methanol as a function of temperature. Filled symbols represent results from the simulations, while open symbols represent experimental results.65–69

atures. This fact could be attributed to the fact that the hydroxyl group has been fitted for short alkanols. In this case, probably an additional hydroxyl group for aromatic cycles should be optimized as well as the quadrupole moment of the aromatic ring accounted for.



TABLE 9: Simulation Results for the Equilibrium Properties of Alcohols T (K) methanol 480 450 390 360 330 300 280 ethanol 480 450 420 390 360 330 300 280 1-propanol 500 475 450 420 390 360 330 300 octanol 600 550 500 450 400 350 325 phenol 650 600 550 500 450 400 350

Fvap (kg/m3)

Fliq (kg/m3)

∆Hvap (kJ/mol)

ln Psat (kPa)

76.55 ( 9 37.3 ( 3 4.81 ( 0.9 1.97 ( 0.3 0.49 ( 0.02 0.112 ( 0.012 0.041 ( 0.0098

528.4 ( 13 599.5 ( 5.2 692.9 ( 4.3 733.0 ( 4.1 761.0 ( 4.6 794.7 ( 4.1 817.8 ( 3.9

16.79 ( 0.7 22.76 ( 0.6 33.53 ( 0.7 35.67 ( 0.6 39.21 ( 0.5 41.81 ( 0.6 43.70 ( 0.2

15.17 ( 0.05 14.55 ( 0.07 12.51 ( 0.06 11.90 ( 0.07 10.57 ( 0.07 8.98 ( 0.06 7.98 ( 0.03

70.3 ( 7 26.7 ( 1.3 11.3 ( 0.8 4.91 ( 0.53 1.44 ( 0.14 0.344 ( 0.05 0.058 ( 0.008 0.0136 ( 0.003

556.0 ( 10 613.2 ( 8 660.5 ( 5.6 706.1 ( 4.5 738.6 ( 6.6 767.0 ( 5.5 791.9 ( 2.1 808.6 ( 4.2

19.86 ( 0.8 25.39 ( 0.7 33.23 ( 0.7 35.70 ( 0.6 41.92 ( 0.6 44.08 ( 0.7 46.57 ( 0.4 48.15 ( 0.3

14.90 ( 0.06 14.12 ( 0.07 13.34 ( 0.06 12.51 ( 0.07 11.28 ( 0.05 9.89 ( 0.05 8.03 ( 0.06 6.54 ( 0.05

79.3 ( 7 39.1 ( 3.8 20.5 ( 1.4 8.79 ( 0.9 3.69 ( 0.9 0.99 ( 0.2 0.33 ( 0.07 0.048 ( 0.01

520.6 ( 8 578.0 ( 7.3 634.9 ( 5.7 678.6 ( 6.4 717.6 ( 4.5 748.4 ( 4.5 777.7 ( 3.1 803.8 ( 2.8

19.60 ( 0.8 26.37 ( 0.7 32.75 ( 0.8 38.36 ( 0.8 40.32 ( 0.8 46.96 ( 0.3 49.29 ( 0.4 51.08 ( 0.6

14.89 ( 0.09 14.25 ( 0.09 13.71 ( 0.09 12.98 ( 0.08 11.95 ( 0.09 10.78 ( 0.08 9.59 ( 0.08 7.67 ( 0.07

97.0 ( 11 34.9 ( 2 10.2 ( 1.3 2.85 ( 0.2 0.469 ( 0.01 0.0313 0.001 ( 0.003 ( 0.0004

461.3 ( 12 573.8 ( 3 642.0 ( 5.2 707.6 ( 3.2 754.0 ( 6.4 800.7 ( 3.2 816.1 ( 5.3

19.53 ( 1.0 33.63 ( 0.5 44.79 ( 0.5 55.30 ( 0.9 59.02 ( 0.8 65.98 ( 0.9 67.81 ( 0.9

14.50 ( 0.08 13.74 ( 0.07 12.58 ( 0.09 11.28 ( 0.08 9.38 ( 0.07 6.55 ( 0.08 4.20 ( 0.07

142.2 ( 10 78.6 ( 7 27.2 ( 1.0 10.3 ( 0.5 2.89 ( 0.2 0.396 ( 0.01 0.154 ( 0.02

575.4 ( 10 683.9 ( 8 776.3 ( 5.9 866.0 ( 6.0 933.1 ( 4.9 976.3 ( 7 1022 ( 3.0

17.44 ( 0.7 26.5 ( 0.7 36.98 ( 0.5 45.43 ( 0.6 53.24 ( 0.6 57.37 ( 0.8 61.0 ( 0.3

15.21 ( 0.06 14.74 ( 0.05 13.91 ( 0.08 12.94 ( 0.08 11.63 ( 0.07 9.54 ( 0.06 8.46 ( 0.06

7. Calculation of Henry Constants of Gases in Alcohols One of the advantages of the use of molecular simulation lies in the possibility to estimate Henry constants. Use of the particle insertion method, also referred to as the Widom method,46 associated with statistical bias techniques,47 allows for the simple and “smart” calculation of the chemical potential µi of a species i in either a pure fluid or a mixture. Consequently, the Henry constants of such a mixture can then be calculated. In this section, we present the application of this methodology for the calculation of Henry constants of different gases in alcohols. The Henry constant of five different gases (methane, oxygen, nitrogen, carbon dioxide, hydrogen sulfide) have been determined for four different alcohols (methanol, ethanol, 1-propanol, octanol) as a function of temperature. The Henry constant KH, which has the dimension of pressure, is related in the limit of low concentration to the fugacity according to

KHxi ) fi

(11)

where xi is the concentration of the solute in the liquid phase and fi is its fugacity in equilibrium. Use of the Widom method allows us to calculate the chemical potential µi of the species i, which is related to the fugacity. In

particular, the chemical potential can be calculated by means of the following expression

{

µi ) -kTln

}

V exp(-β∆U+) Ni + 1

(12)

where ∆U+ is the potential-energy difference caused by insertion of the test molecule, V is the volume of the system, and Ni is the number of gas molecules in the solvent. As discussed by Frenkel and Smit,46 the de Broglie wavelength does not appear due to the fact that it cancels out when the ideal gas at the same temperature is taken as the reference state, that is, when µi ) 0. According to eq 12, these conditions are satisfied for a fluid of unit density (V/(Ni) + 1 ) 1) and zero interaction energy (∆U+ ) 0), which means an ideal gas of unit density. Although this reference state has no physical meaning, it makes the relationship with gas fugacity explicit through the equation of state of a system of (N + 1) molecules of an ideal gas

P0 )

(N + 1)kT

V

) kT

(13)

Then the chemical potential defined by eq 11 can be written as



() fi

µi ) kTln

(14)

0

P

Since the implicit reference pressure P0 is not constant, it is more convenient to use a reference state with fixed pressure Pref ) 1 Pa, i.e., the unit of the international system. Once we have defined this new reference state, we can also define a chemical potential µ j i which is related to µi according to

( )

( )

µ j i ) kTln

fi P0 + µi ) kTln Pref Pref

(15)

Making use of eqs 12 and 13 we can write

( )

µ j i ) kTln

{

}

P0 V - kTln exp(-β∆U+) Pref Ni + 1

(16)

The Henry constant in international system units can be determined by applying eqs 11 and 15 and considering that xi ) (Ni + 1)/(N + 1), where N is the total number of molecules in the solvent (N + 1)exp( µ j i ⁄ kT) fi KH ) ) Pref Prefxi Ni + 1

(17)

The Henry constant can then be determined in the limit of infinite dilution by making test insertions of the different gas molecules in the solvent. It is important to remark that a significant advantage of this methodology is that different Henry constants can be determined from the same single simulation. This expression works correctly for methane, since it can be considered as a single center Lennard-Jones fluid. In the case of the other gases, which involve several force centers and electrostatic charges, as previously mentioned, we need to use statistical bias methods.48 Finally, the variation of the Henry constant with temperature is related to the infinite dilution solvation Gibbs energy, ∆G∞solv,49 according to

( ) ln KH

∞ ∆Gsolv ) RT

(18)

P0

and through its temperature derivatives to infinite dilution ∞ solvation enthalpy ∆Hsolv

(

∞ ∆Hsolv ) -RT2

d ln KH dT

)

(19)

The infinite dilution solvation quantity corresponds to the transfer of one mole of the solute from the pure ideal gas state at standard pressure P0 ) 100 kPa to a hypothetical infinitely dilute solution. This quantity can also be determined from molecular simulation50 by estimating the free energy with HR + GCMC simulations along a thermodynamic path. 7.1. Intermolecular Potentials for Gases. In Table 8 are given the potential parameters of the different models for the five different gases considered. To reproduce the methane molecule, the model proposed by Mo¨ller et al. has been selected.51 It involves a single LJ force center with the absence of electrostatic charges. In the case of nitrogen, the model proposed by Delhomelle52 has been chosen. It involves two LJ force centers separated by a fixed distance of 1.098 Å as well as two negative electrostatic charges located on the atomic centers and one positive charge on the center of mass. These charges were fitted respecting the experimental quadripolar moment of the molecule. The LJ parameters were fitted to reproduce the coexistence curve of

the nitrogen molecule, thus providing an excellent account of the liquid-vapor coexistence densities of pure nitrogen. The model selected for oxygen has been proposed by Vrabec et al.,53 who determined the parameters of several fluids with a two LJ centers plus point quadripole moment on the basis of the liquid-vapor coexistence properties. As our software does not at present allow for point quadripoles, three electrostatic charges have been introduced as in the case of nitrogen or carbon dioxide but with a very close spacing of 0.2 Å. The values of these charges were determined in such a way that the quadripolar moment indicated by Vrabec et al., i.e., 0.8081 D Å, was exactly respected. For carbon dioxide, the selected model is the EPM2 potential of Harris and Yung,54 which involves three LJ force centers and three electrostatic charges on the atomic centers. Finally, the model selected for hydrogen sulfide is that proposed by Kristof and Liszi55 involving one single LJ center of force and four electrostatic charges, three of them at the atomic centers of the different atoms of the molecule and the fourth on the bisector of the angle formed by the hydrogens and the sulfur. Although Lorentz-Berthelot combining rules were used in the development of the AUA potential, we used Kong56–58 combining rules to evaluate the cross-interactions between unlike groups. These rules are

εijσij6 ) √εiσi6εjσj6 εijσij12 )

εiσi12 213

[ ( )] 1+

εjσj12

εiσi12

(20)

1⁄13 13

(21)

Use of Kong combining rules instead of Lorentz-Berthelot rules has been previously well justified for mixtures of carbon dioxide with alkanes.59 Also, these rules have been found to perform better for binary mixtures of noble gases60 and carbon dioxide with hydrogen sulfide.61 Additional simulations were carried out for the case of ethanol with the Lorentz-Berthelot rule, and the results are given in Figure 15 as shaded symbols. As can be seen from Figure 15, although the Lorentz-Berthelot rule improves agreement with experimental results for methane, the Kong rule appears to give better results for oxygen. For the other gases it is not clear which rule gives a closer agreement with the experimental value. 7.2. Results. The results obtained for the different gases in the different solvents are shown in Figures 14, 15, 16, and 17. The Henry constants are ranked in the same order for a given temperature and solvent for all figures. The Henry constant decreases as the critical temperature of the different solutes increases. In particular, nitrogen exhibits the highest critical temperature followed by oxygen, methane, carbon dioxide, and hydrogen sulfide. The explanation for this phenomenon is that the higher the critical temperature of the solute is, the larger the attraction energy in the liquid phase and the higher the solubility in the solvent. This behavior is also in good agreement with the available reference data for ethanol.19 In the case of the temperature dependence, it differs from one component to the next. Except for octanol where all the solutes increase their solubilities as temperature is increased, in the other cases the Henry constant decreases with increasing temperatures (positive ∆H∞solv) for nitrogen, oxygen, and methane while it increases with increasing temperatures (negative ∆H∞solv) for carbon dioxide and hydrogen sulfide. Further investigation of the liquid structure would be required to understand the origin of the specific temperature dependence shown by these com-

9862 J. Phys. Chem. B, Vol. 112, No. 32, 2008 pounds. In the case of ethanol, where experimental data are available, the obtained results are in good quantitative agreement for the different solutes and the predicted variations with temperature are reproduced satisfactorily. In the case of methane and oxygen, the effect of temperature on the Henry constant is found to be insignificant if statistical uncertainties are accounted for, in agreement with experiment. Although it is generally recommended to use staged insertions for evaluating chemical potentials,62 in our case we obtained reasonable estimations of the Henry constants for the different solvents. This can be explained because the temperatures are not too low, except in the case of octanol, where the lowest temperature considered is 130 K below its normal boiling point. This means that, in general, at these conditions the liquid structure still presents sufficient “holes” to accommodate the small solutes we considered and a reasonable sampling can be obtained. The second explanation is the amount of computational time dedicated to the sampling of the solvent: at least 100 million MC moves in all cases and up to 400 million MC steps for octanol at the lower temperatures. This fact can be illustrated by the average statistical deviations obtained, whose values are around 5%. Staged insertion would probably have been required when considering significantly larger solutes or significantly lower temperatures, where the computational time necessary required to apply the Widom test at the same conditions obtained using the staged deletion method is significantly higher.63,64 8. Conclusions The AUA4 potential has been extended to alcohols by adjusting the LJ parameters of the hydroxyl group and optimizing a set of charges so that they reproduce the electrostatic distributions of methanol and ethanol. The proposed distribution is intended to better describe the totality of the phase diagram by introducing extra nucleic charges that are able to better match the electrostatic field of methanol and ethanol. Although the behavior of the model could be improved individually by optimizing different sets of LJ parameters for the carbon groups next to the hydroxyl, our goal is to produce parameters that are transferable to other molecules. In this way, we maintained the LJ parameters of the groups adjacent to the hydroxyl in order to maintain as far as possible the transferability of the AUA4 model while still reproducing the equilibrium properties with an accuracy similar to previously proposed models. While the predictions for the saturated liquid densities and the vapor pressures agree well with the experimental data, the model has more difficulties to reproduce the heats of vaporization. Contrary to the TraPPE model, our model compares better near to the critical region while presenting higher deviations at low temperatures. To improve this general behavior, increasing the level of detail of the model by introducing variable charges methodologies may be considered. In the case of phenol, we think that the critical point prediction can be improved by accounting for the quadrupole moment associated with the aromatic ring, as done in other models of aromatics.74–77 From the finite size scaling study, we have shown that the presence of the electrostatic charges in this model does not change the universality class, belonging still to the Ising class. We also calculated the Henry law constants of different gases in different alcohols by applying the Widom particle insertion method and obtained very good agreement with the available experimental data. In particular, the temperature dependence of the constants is well reproduced. Acknowledgment. This study was supported by the Spanish Government, Ministerio de Ciencia e Innovacioń (Grant No.

Pe´rez-Pellitero et al. CTQ2008-06469) as well as the Institut Français du Pe´trole. The simulations have been carried out using the GIBBS Monte Carlo code developed jointly by the Institut Français du Pe´trole, Universite´ de Paris Sud, CNRS and other partners. The authors would like to acknowledge Bernard Levy (LCP, Orsay) for his contribution to the determination of charges on methanol and ethanol. Note Added after ASAP Publication. This paper was published ASAP on July 23, 2008. Equations 20 and 21 were changed. The revised paper was reposted on August 7, 2008. References and Notes (1) Nain, A. K. J. Solut. Chem. 2007, 36, 497. (2) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 3093. (3) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638. (4) Smit, B.; Karaborni, S.; Siepmann, J. I. J. Chem. Phys. 1995, 102, 2126. (5) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569. (6) Errington, J. R.; Panagiotopoulos, A. Z. J. Phys. Chem. B 1999, 103, 6314. (7) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 5370. (8) Ungerer, P. J. Chem. Phys. 2000, 112, 5499. (9) Toxvaerd, S. J. Chem. Phys. 1997, 107, 5197. (10) Toxvaerd, S. J. Chem. Phys. 1990, 107, 4290. (11) Jorgensen, W. L. J. Am. Chem. Soc. 1980, 102, 543. (12) Jorgensen, W. L. J. Am. Chem. Soc. 1981, 103, 335. (13) Haughney, M. J. Phys. Chem. 1987, 91, 4934. (14) Van Leeuwen, M. E. Mol. Phys. 1996, 87, 87. (15) Van Leeuwen, M. E. J. Phys. Chem. 1995, 99, 1831. (16) Siepmann, J. I. Nature 1993, 365, 330. (17) Cichowski, E. C.; Schmidt, T. R. Fluid Phase Equilib. 2005, 236, 58. (18) Zhang, J.; Siepmann, J. I. Theor. Chem. Acc. 2006, 115, 391. (19) Boutard, Y.; Ungerer, Ph.; Teuler, J. M.; Ahunbay, M. G.; Sabater, S. F.; Pe´rez-Pellitero, J.; Mackie, A. D.; Bourasseau, E. Fluid Phase Equilib. 2005, 236, 25. (20) Ahunbay, M. G. J. Phys. Chem B 2005, 109, 2970. (21) Contreras-Camacho, R. O. J.Phys. Chem. B 2004, 108, 14115. (22) Cornell, W. D. J. Am. Chem. Soc. 1995, 117, 5179. (23) Jorgensen, W. L. J. Phys. Chem. 1986, 90, 1276. (24) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129. (25) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, P.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; HeadGordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, Revision A.6; Gaussian, Inc.: Pittsburgh PA, 1998. (26) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Phys. Chem. 1983, 79, 926. (27) CRC Handbook of Chemistry and Physics, 75th ed.; CRC Press, Inc.: Boca Raton, FL,1994. (28) Mackie, A. D. Mol. Simul. 1997, 19, 1. (29) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (30) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1988, 61, 2635. (31) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1989, 63, 1195. (32) Bourasseau, E. J. Chem. Phys. 2003, 118, 3020. (33) Pe´rez-Pellitero, J. J. Chem. Phys. 2006, 125, 054515. (34) Pe´rez-Pellitero, J. J. Chem. Phys. 2007, 126, 079901. (35) Wilding, N. B. Phys. ReV. E 1995, 52, 602. (36) Ferrenberg, A. M.; Landau, D. P. Phys. ReV. B 1991, 44, 5081. (37) BYU DIPPR 801,Thermophysical Properties Database Public Release, January 2005 (38) Stoll, J. Fluid Phase Equilib. 2003, 209, 29. (39) Ketko, M. B. H.; Potoff, J. J. Mol. Simul. 2007, 33, 769. (40) Mortier, W.; Gosh, S.; Shankar, S. J. Am. Chem. Soc. 1986, 108, 4315.

Anisotropic United-Atoms Potential for Alcohols (41) Rappe, A.; Goddard, W. J. Chem. Phys. 1991, 95, 3358. (42) Bourasseau, E.; Maillet, J.-B.; Mondelain, L.; Anglade, P.-M. Mol. Simul 2005, 31, 705. (43) Private communication with Desbiens, N.; Bourasseau, E.; Maillet, J. B. from CEA-DAM DPTA in France,to be published. (44) Yamaguchi, T.; Hidaka, K.; Soper, A. K. Mol. Phys. 1999, 96, 1159. (45) Yamaguchi, T.; Hidaka, K.; Soper, A. K. Mol. Phys. 1999, 97, 603. (46) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Francisco, 2001. (47) Bourasseau, E. J. Phys. Chem. B 2002, 106, 5483. (48) Ungerer, P.; Tavitian, B.; Boutin, A. Applications of Molecular Simulation. In The Oil And Gas Industry; Editions Technip: Paris, 2005. (49) Dohnal, V.; Fenclova´, D.; Vrbka, P. J. Phys. Chem. Ref. Data 2006, 35, 1621. (50) Hernańdez-Cobos, J.; Mackie, A. D.; Vega, L. F. J. Chem. Phys. 2001, 114, 7527. (51) Moller, D.; Oprzynski, J.; Muller, A.; Fischer, J. Mol. Phys. 1992, 75, 363. (52) Delhommelle, J. Etablissement de potentiels d’interaction pour la simulation molećulaire-application a la pre´diction des e´quilibres liquideVapeur de me´langes binaires alcane-molećule multiolaire, Universite´ de Paris XI Orsay, France, 2000. (53) Vrabec, J.; Stoll, J.; Hasse, H. J. Phys. Chem. B 1995, 105, 12126. (54) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (55) Kristof, T. J.Phys. Chem. B 1997, 101, 5480–5483. (56) Kong, C. L. J. Chem. Phys. 1973, 59, 1953. (57) Kong, C. L. J. Chem. Phys. 1973, 59, 2464. (58) Kong, C. L. J. Phys. Chem. 1973, 77, 2668.

J. Phys. Chem. B, Vol. 112, No. 32, 2008 9863 (59) Potoff, J. J.; Errington, J. R.; Panagiotopoulos, A. Z. Mol. Phys. 1999, 97, 1073. (60) Delhommelle, J.; Millie, P. Mol. Phys. 2001, 99, 619. (61) Ungerer, P.; Wender, A.; Demoulin, G.; Bourasseau, E.; Mougin, P. Mol. Simul. 2004, 30, 631. (62) Kofke, D. A. Mol. Phys. 2004, 102, 405. (63) Boulougouris, G. C.; Economou, I. G.; Theodorou, D. N. J. Chem. Phys. 2001, 115, 8231. (64) Boulougouris, G. C.; Economou, I. G.; Theodorou, D. N. Mol. Phys. 1999, 96, 905. (65) Ukai, T.; Kodama, D.; Miyazaki, J.; Kato, M. J. Chem. Eng. Data 2002, 47, 1320. (66) Boyer, F. L.; Bincher, L. J. J. Phys. Chem. 1960, 64, 1330. (67) Bo, S.; Battino, R.; Wilhelm, E. J. Chem. Eng. Data 1993, 38, 4. (68) Fischer, K.; Wilken, M. J. Chem. Thermodyn. 2001, 33, 1285. (69) Ohgaki, K.; Katayama, T. J. Chem. Eng. Data 1976, 21, 53. (70) Tokunaga, J. J. Chem. Eng. Data 1975, 20, 41. (71) Chang, C. M. J.; Chiu, K. L.; Day, C. Y. J. Supercrit. Fluids 1998, 12, 223. (72) Weng, W. L.; Lee, M. J. Fluid Phase Equilib. 1992, 73, 117. (73) Weng, W. L.; Lee, J. M. Fluid Phase Equilib. 1996, 122, 243. (74) Bonnaud, P.; Nieto-Draghi, C.; Ungerer, P. J. Phys. Chem. B 2007, 111, 3730. (75) Nieto-Draghi, C.; Bonnaud, P.; Ungerer, P. J. Phys. Chem. C 2007, 111, 15686. (76) Nieto-Draghi, C.; Bonnaud, P.; Ungerer, P. J. Phys. Chem. C 2007, 111, 15942. (77) Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2007, 111, 10790.

JP802282P

Anisotropic United-Atoms (AUA) Potential for Alcohols - The Journal of

Recommend Documents