Transferable Force Field for Alcohols and Polyalcohols - The Journal

Apr 3, 2009 - The Journal of Physical Chemistry B .... UMR 8000 CNRS, Université Paris-Sud. ... de Hemptinne , Rafael Lugo , Nicolas Ferrando , and J...
0 downloads 0 Views 443KB Size
Download PDF
J. Phys. Chem. B 2009, 113, 5985–5995

5985

Transferable Force Field for Alcohols and Polyalcohols Nicolas Ferrando,*,†,‡ Ve´ronique Lachet,† Jean-Marie Teuler,‡ and Anne Boutin‡ De´partement Thermodynamique et Mode´lisation Mole´culaire, IFP, 1-4, AVenue de Bois Pre´au, 92 852 Rueil-Malmaison Cedex, France, and Laboratoire de Chimie Physique, UMR 8000 CNRS, UniVersite´ Paris-Sud, baˆt. 349, 91405 Orsay Cedex, France ReceiVed: December 11, 2008; ReVised Manuscript ReceiVed: February 25, 2009

A new force field has been developed for alcohol and polyalcohol molecules. Based on the anisotropic unitedatom force field AUA4 developed for hydrocarbons, it only introduces one new anisotropic united atom corresponding to the hydroxyl group OH. In the case of polyalcohols and complex molecules, the calculation of the intramolecular electrostatic energy is revisited. These interactions are calculated between charges belonging to the different local dipoles of the molecule, one dipole being defined as a group of consecutive charges globally neutral. This new method allows avoiding the use of empirical scaling parameters commonly introduced to calculate 1-4 electrostatic interactions. The transferability of the proposed potential is demonstrated through the simulation of a wide variety of alcohol families: primary alcohols (methanol, ethanol, propan-1-ol, hexan-1-ol, octan-1-ol), secondary alcohols (propan-2-ol), tertiary alcohols (2-methylpropan-2ol), phenol, and diols (1,2-ethanediol, 1,3-propanediol, 1,5-pentanediol). Monte Carlo simulations carried out in the Gibbs ensemble lead to a good agreement between calculated and experimental data for the thermodynamic properties along the liquid/vapor saturation curve, for the critical point coordinates and for the liquid structure at room temperature. Additional simulations were performed on the methanol + n-butane system showing the capability of the proposed potential to reproduce the azeotropic behavior of such mixtures with a good agreement. 1. Introduction Alcohols play an important role in chemical industry. Due to their associating nature characterized by a strong ability to form hydrogen bonds, alcohols are for instance a widely used class of solvents in organic chemistry or in the biochemical industry. Alcohols are also involved in the production of biofuels, since they can be directly added to gasoline or gasoil blends, or be implicated in the biomass treatment processes.1 Therefore, a good knowledge of phase properties and phase equilibrium diagrams of systems containing alcohols is of first importance to design and optimize a large number of industrial processes. However, the experimental data required to model these systems are not always available and are often expensive or difficult to obtain. Predictive simulation tools then appear as an alternative solution. Molecular simulation is becoming an efficient method to accurately describe thermodynamic properties of a wide variety of pure components and mixtures.2 The accuracy of the simulations is intimately related to the ability of the force field to reproduce the intermolecular interactions. A common type of potential used in molecular simulation is based on the concept of united atoms (UA), in which a unique center of force is placed on carbon atoms that implicitly takes into account the bonded hydrogen atoms. An alternative is the anisotropic united atom (AUA) model, in which the group center is shifted from the carbon atom toward the group center of mass.3 A united atom force field for alcohols was developed by Jorgensen4 as an extension of the optimized potential for liquid * Corresponding author. Tel.: +33 147526624. Fax: +33 147527025. E-mail: [email protected]. † IFP. ‡ UMR 8000 CNRS, Universite´ Paris-Sud.

simulations (OPLS-UA) model5 in order to describe the liquidphase structure of short primary, secondary, and tertiary alcohols at room temperature. It has been shown6 that this model can also be used to predict phase properties along the liquid/vapor saturation curve for alcohols up to propanol. However, it becomes less accurate for longer alcohols. Van Leuween6 proposed an extension of the Smit-Karaboni-Siepmann (SKS) model7 to primary and secondary alcohols. Although accurate predictions of phase equilibrium were obtained for alcohols up to hexanol, this model requires a fine Lennard-Jones parameters tuning for each simulated alcohol. Hence, this approach cannot be considered as fully transferable. Chen et al.8 proposed an extension of the transferable potential for phase equilibria (TraPPE-UA) force field to primary, secondary, and tertiary alcohols, also transferable to polyalcohols. Nowadays, it is probably the most accurate transferable model for the prediction of phase properties along the saturation curve for all the alcohols investigated. This force field requires a special treatment for the simulation of polyalcohols or molecules involving both a hydroxyl group and another oxygenated function. Indeed, in such cases, an empirical scaling factor of 1/2 is introduced in the intramolecular electrostatic energy contribution between 1-4 bonded charges, and a specific repulsive potential is added between the functional hydrogen and the second oxygen of the molecule.9 More recently, a first attempt to extend the AUA4 potential10 to primary alcohols was carried out,11 but although saturated liquid densities, vaporization enthalpies, and critical properties are well predicted with this model, a significant deviation on vapor pressures is obtained. This model assumes a rigid bending angle between the carbon, oxygen, and hydrogen atoms and a four-charges distribution is introduced to account for the polarity of the molecules.

10.1021/jp810915z CCC: $40.75  2009 American Chemical Society Published on Web 04/03/2009

5986 J. Phys. Chem. B, Vol. 113, No. 17, 2009

Ferrando et al.

TABLE 1: Nonbonded Parameters for the AUA4 Force Field atom

ε (K)

σ (Å) δ (Å)

q (e)

TABLE 2: Bonded Parameters for the AUA4 Force Field

ref

120.15 3.607 0.216 +0.265 if bonded 0 elsewhere 86.29 3.461 0.384 +0.265 if bonded CH2 0 elsewhere CH 50.98 3.363 0.646 +0.265 if bonded 0 elsewhere C 15.04 2.44 0 +0.265 if bonded 0 elsewhere CHarom 89.4 3.246 0.407 +0.265 if bonded 0 elsewhere Carom 37.7 3.246 0 +0.265 if bonded 0 elsewhere O (OH) 125.01 3.081 0.010 -0.700 H (OH) 0 0 0 +0.435 CH3

[( ) ( ) ] 12

-

σij rij

ref

1.535 1.40 1.430 0.945

2 2 4 4

2

to O

2

CHx-CHy CaromdCarom CHx-O O-H

to O

2

bend

θ0 (deg)

kbend (K)

ref

to O

2

to O

2

to O

2

CHx-CH2-CHx CHx-CH-CHx CHx-C-CHx CHaromdCHaromdCHarom CHaromdCaromsO CHx-CHy-O CHx-OH

114.0 112.0 109.47 120.0 120.0 109.47 108.5

74900 72700 70311 rigid rigid 59800 61000

2 2 2 2 13 13 13

this work this work

2. Force Field Development 2.1. Molecular Model. Intermolecular Energy. The dispersive-repulsive interaction energy between two force centers i and j of different molecules is modeled through a 12-6 LennardJones potential:

σij rij

r0 (Å)

to O

In this work, we propose a new force field for primary, secondary, and tertiary alcohols and polyalcohols as an extension of the AUA4 potential. To ensure transferability, only one new anisotropic united atom is introduced corresponding to the hydroxyl group OH. All other pseudoatoms involved in the hydrocarbonated part of the molecules are directly transferred without any modification from the AUA4 original model2 with the aim to ensure transferability to a large variety of molecules. In the case of polyalcohols, a new methodology is introduced to account for the intramolecular electrostatic energy. It avoids the use of a scaling factor or an empirical additional potential to remain as transferable as possible. This paper is organized as follows: the proposed force field is presented in section 2, with a detailed description of the intramolecular electrostatic and repulsive-dispersive interactions. The simulation results of a wide variety of alcohols are discussed in section 3, highlighting the transferability of this potential and its ability to predict phase equilibrium and structure of pure alcohols as well as an alkane + alcohol mixture. Finally, section 4 gives our conclusions.

U ijLJ ) 4εij

bond length

6

(1)

Cross Lennard-Jones parameters are obtained using LorentzBerthelot combining rules:

εij ) √εiiεjj

(2)

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

(3)

All the Lennard-Jones parameters involved in this work are taken from the AUA4 potential developed for hydrocarbons2,12 and summarized in Table 1. Only one new anisotropic united atom is introduced corresponding to the hydroxyl group OH. Following the AUA principle, the Lennard-Jones sphere is not centered on the oxygen atom, but is slightly shifted on the oxygen-hydrogen bond toward the hydrogen atom to account for the hydrogen atom volume. The Lennard-Jones parameters

torsion

ai (K)

CHx-CH2-CH2-CHy a0 ) 1001.35 a2 ) -303.06 a4 ) 2226.71 a6 ) -4489.34 a8 ) 2817.37 CHx-CH2-CH2-O a0 ) 839.87 a2 ) 106.68 CHx-CH2-O-H a0 ) 339.41 a2 ) 58.34 CHx-CH-O-H a0 ) 302.29 a2 ) -62.92 CHx-C-O-H a0 ) 163.56 a2 ) 0 CHaromdCaromsOsH a0 ) 845.65 a2 ) -845.65 O-CH2CH2-O a0 ) 1530.22 a2 ) -422.78 a4 ) -772.20 a6 ) 174.93 a8 ) 534.83

a1 a3 a5 a7

ref ) ) ) )

2129.52 -3612.27 1965.93 -1736.22

2 2 2 2 2 a1 ) 2133.17 4 a3 ) -3079.72 4 a1 ) 353.97 4 a3 ) -751.72 4 a1 ) -719.09 4 a3 ) 695.68 4 a1 ) 490.68 4 a3 ) -654.24 4 a1 ) 0 14 a3 ) 0 14 a1 ) 4064.26 this work a3 ) -5770.19 this work a5 ) 1110.27 this work a7 ) 430.66 this work this work

σOH and εOH, and the AUA displacement δOH of this new united atom are thus the only three adjustable parameters. For electrostatic energy calculation, two charge distributions for united atoms force fields are described in the literature: the first one, originated from the OPLS-UA model, consists of two positive charges located on the hydrogen atom and on the R-carbon atom, and one negative charge on the oxygen atom. This distribution and the corresponding charge values lead to an overestimate of about 25% of the alcohol dipole moment in the vapor phase due to an implicit account of polarization effect.4 This approach has been adopted in the TraPPE force field.8 The second charge distribution is originated from the previous AUA4 model for alcohols.11 It is composed of four charges, not necessarily located on atoms, which have been adjusted to mimic the electrostatic field surrounding an isolated molecule of ethanol. Since one additional charge is introduced compared to the OPLS-UA force field, this distribution leads to an increase of the computation time, but it is not systematically associated with an improvement of the model accuracy. Consequently, we have decided in this work to adopt the usual charge distribution of the OPLS-UA force field (see Table 1). Intramolecular Energy. All the bonded parameters involved in the intramolecular energy calculation are given in Table 2. In the proposed force field, all bond lengths are kept fixed. The CHx-O and O-H bond lengths are taken from the OPLS-UA force field,4 and all other bonds from the AUA4 force field.2 Atoms separated by two bonds interact via a harmonic bending potential

Transferable Force Field for Alcohols

Ubend 1 ) kbend(cos θ - cos θ0)2 kB 2

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5987

(4)

where kB is the Boltzmann constant, kbend is the bending constant, and θ and θ0 are the bending angle and the equilibrium bending angle, respectively. The parameters involving atoms of the hydroxyl group are taken from the AMBER force field.13 Other bending parameters are taken from the AUA4 force field.2 For atoms separated by three bonds, a torsional potential of the following form is used:

Utors ) kB

8

∑ an(cos φ)n

(5)

n)0

where φ is the dihedral angle and ai the ith parameter. As done in the previous AUA4 force field for alcohols11 and in the TraPPE force field,8 the parameters for torsions involving the oxygen or hydrogen atoms of the hydroxyl group are taken from the OPLS-UA force field.4 For phenol, the torsion parameters of the CHarom-Carom-O-H dihedral angle are taken from the OPLS-AA force field.14 As explained further, the O-CH2CH2-O torsion encountered in 1,2-ethanediol was reparametrized in this work. All other torsion parameters are taken from the AUA4 force field.2 A 12-6 Lennard-Jones potential is used to calculate repulsive-dispersive intramolecular energy between two force centers separated by more than three bonds (eq 1). For polyalcohols where the two hydroxyl groups are separated by three bonds (1,2-ethanediol, for instance), it is well-known that difficulties related to a too strong attraction between these two hydroxyl groups occur when simulations are carried out with a united atom force field. This is due to the fact that the hydroxyl hydrogen atoms are only represented by a partial electrostatic charge and are not surrounded by a Lennard-Jones sphere. One method to overcome this difficulty is to introduce an empirical repulsive potential between the hydroxyl hydrogen and the nonbonded oxygen atom (see the TraPPE force field, for instance9). However, this approach implies the use of an additional empirical parameter. To avoid such additional parameter, we have chosen to explicitly calculate repulsivedispersive interactions between two hydroxyl groups separated by three bonds using the same Lennard-Jones potential as the one used for nonbonded interactions. Taking explicitly into account the Lennard-Jones 1-4 oxygen-oxygen interaction, as well as using a new method for intramolecular electrostatic interaction (see below), implies to reparametrize the O-CH2-CH2-O torsional potential. The corresponding parameters have been determined from ab initio calculations on an isolated 1,2-ethanediol molecule using Gaussian 03, with the MP2 method and 6-311G(d-p) basis set. The energy was obtained in the whole range of torsional angles, from 0 to 360°, with increments of 10°. For each configuration, the contributions of bending energies, of torsional energies involving other dihedral angles, and of intramolecular electrostatic and Lennard-Jones energies are subtracted from the total ab initio energy. The resulting O-CH2-CH2-O torsional energy is fitted with the analytical function given in eq 5. This torsional energy and the corresponding fit are shown in Figure 1. The O-CH2-CH2-O torsional potential used in TraPPE force field9 is also presented in this figure. The energetic barrier differences between these two potentials arise from the explicit use of 1-4 Lennard-Jones interactions between oxygen atoms

Figure 1. Torsional energy as a function of the O-CH2-CH2-O angle. An angle equal to 0 corresponds to a trans conformation, whereas an angle equal to 180° corresponds to a cis conformation. Comparison between values derived from ab initio calculations (symbols), adjusted function with eq 5 (solid line), and torsional potential from Stubbs et al.9 (dashed line).

Figure 2. Schematic representation of the 1,3-propanediol molecule. Carbon atoms are represented in light gray, oxygen atoms in dark gray, and hydrogen atoms in white. q1 to q6 represent the different partial electrostatic charges used to account for the polarity of the molecule.

in our case. In the cis conformation (dihedral angle equal to 180° in the convention used in this work), the distance between these two atoms is minimal, leading to a positive 1-4 LennardJones energy, and thus to a lower energetic barrier in the torsional potential. On the contrary, in the trans conformation, the 1-4 Lennard-Jones energy is negative and thus the torsional energy of this conformation is higher in our potential compared to the TraPPE force field. Finally, the calculation of intramolecular electrostatic energy has been revisited in this work. This energy is usually calculated using Coulombic interactions between all individual partial charges on atoms separated by more than three bonds. With such an approach, it has been shown that some corrections should be added to avoid unphysical property predictions, such as a 1/2-scaled contribution of the 1-4 electrostatic energy, as done in OPLS-AA,14 AMBER,13 or TraPPE-UA9 force fields. Since the partial electrostatic charges distribution aims to mimic the local dipoles of the molecule, we propose to calculate the intramolecular electrostatic energy between charges belonging to different local dipoles of the molecule, even if they involve bonded charged atoms, and without introducing empirical scaling parameters. One local dipole is defined as a group of consecutive charges globally neutral. For example, in the case of 1,3-propanediol, two dipoles formed by the two groups of charges {q1, q2, q3} and {q4, q5, q6}, respectively, are present, as illustrated in Figure 2. Hence, the three charges q1, q2, and q3 interact with the charges q4, q5, and q6. This means that the charges q2 and q4, q3 and q4, as well as q3 and q5 interact although belonging to bonded atoms. It is worth mentioning that the use of charge-charge formula for the calculation of dipole-dipole interactions breaks down at close proximity. This is true for both intramolecular and intermolecular energies. It

5988 J. Phys. Chem. B, Vol. 113, No. 17, 2009

Ferrando et al.

is also worth recalling here that the bonded terms have been revisited to be consistent with the electrostatic interactions computed with this new approach. Indeed, for 1,2-ethanediol, as explained before, the O-CH2-CH2-O torsion energy was determined from the total intramolecular energy from which all other bonded terms, intramolecular electrostatic and LennardJones energies have been subtracted. The H-O-CH2-CH2 torsion energy was not reparametrized since the charge distribution mimics a dipole localized close to the oxygen atom. This torsion is therefore not directly affected by the intramolecular electrostatic energy. As shown further in diols simulations, this new approach prevents any unphysical conformations and does not require any additional interaction terms. 2.2. Lennard-Jones Parameters Optimization. In the proposed force field, only three new parameters have been determined: the Lennard-Jones parameters σOH and εOH and the AUA displacement δOH of the new united atom describing the hydroxyl group. The adjustment procedure used in this work has been described elsewhere.15 It consists of minimizing the mean square relative deviation F between experimental and calculated properties Xi: calc

(X i 1 F) n i)1 n

∑

- X iexp)2 si2

(6)

where n denotes the total number of reference data and s the statistical uncertainty. The properties used in this adjustment procedure are the saturated liquid densities, vapor pressures, and vaporization enthalpies with associated statistical uncertainties of 1.5, 5, and 2%, respectively. The components taken into account in the adjustment procedure are methanol and ethanol at 325 and 450 K and 1,3-propanediol at 500 and 600 K. The final hydroxyl group parameters are given in Table 1. 3. Results 3.1. Simulation Details. All the liquid/vapor equilibrium simulations of pure components were carried out in the NVT Gibbs ensemble.17,18 A total number of 300 molecules were used for each system, except in the vicinity of the critical point where larger systems were used (up to 700 molecules). Using double size systems for different reduced temperatures, system-size effects have been estimated to be around 0.6% for liquid density, 6% for vapor pressure, and 0.6% for vaporization enthalpy. These values are within the statistical uncertainties on the studied temperature range (see Table S.1 in the Supporting Information). For monoalcohols, simulation runs lasted for 40 million steps, one step corresponding to a single Monte Carlo move. For diols, the convergence to reach equilibrium was slower, and at least 50 million steps were required. For Lennard-Jones interactions, a spherical cutoff equal to the half of the simulation box was used, and classical tail correction was employed.19 For longrange electrostatic energy, the Ewald sum technique was used, with a number of reciprocal vectors k ranging from -7 to 7 in all three directions and a Gaussian width R equal to 2 in reduced units. The different Monte Carlo moves and their attempt probabilities that were used during the simulations are: translation (20%), rigid rotation (20%), configurational-bias regrowth20 (20%), transfer (39.5%) with insertion bias,21 and volume change (0.5%). The amplitude of translations, rotations, and volume changes was adjusted during the simulation to achieve an acceptance ratio of 40% for these moves. The critical properties TC and FC of pure components were estimated by a least-squares fit of the law of rectilinear diameters:

Fl + Fv ) FC + A(T - TC) 2

(7)

where Fl and Fv are the density of the liquid and vapor phases, respectively, T is the temperature, and A is an adjustable parameter, and the critical scaling relation

Fl - Fv ) B(T - TC)β

(8)

where B is another adjustable parameter and β the universal exponent, equal to 0.325.20 This procedure works rather well for nonpolar molecules, but its application to polar systems such as alcohols is, however, questionable.22 A solution would consist of systematically adjusting the value of the universal exponent to the molecule studied.8 Following previous works, we decided in this work to keep a constant value of 0.325. Finally, the normal boiling temperatures of pure components were determined using the Clausius-Clapeyron equation. For the methanol + n-butane mixture, simulations were carried out in the NPT Gibbs ensemble.18,23 A total number of 400 molecules and 50 million steps were used for each pressure-temperature condition, except in the vicinity of the azeotrope where a larger number of molecules (up to 700) and longer simulations (80 million steps) were used. Other simulation parameters are similar to those previously described for pure component systems. The azeotropic composition and pressure were determined following a simple trial and error procedure. The azeotropic composition was first estimated from the pressure-composition diagram built from the Gibbs NPT simulations. Gibbs NVT simulations were then performed until equal compositions in the liquid and the vapor phases were obtained within statistical uncertainties. The simulations performed in this work were carried out using the in-house GIBBS Monte Carlo code jointly developed by IFP and the Laboratoire de Chimie Physique (CNRS and Universite´ Paris-Sud).2 The required CPU time for 10 million MC steps is typically 4 h using four processors (2.4 GHz). 3.2. Performance of the Force Field on Molecules Used for Adjustment. Figure 3 shows the saturated liquid densities, vapor pressures, and vaporization enthalpies of methanol, ethanol, and 1,3-propanediol. In this figure, reference data used in the adjustment procedure are plotted with filled symbols, whereas predictions obtained at other thermodynamic conditions are plotted with open symbols. Numerical values are given in the Supporting Information (Table S.1). For each of these thermodynamic properties, a good agreement is found for methanol and ethanol (deviations less than 2, 13, and 2% for saturated liquid densities, vapor pressures, and vaporization enthalpies, respectively). Vapor pressures of 1,3propanediol are remarkably well predicted (deviation about 9%), and saturated liquid densities and vaporization enthalpies are slightly underestimated for this compound (deviation about 6 and 9%, respectively). This set of parameters is the best compromise to simultaneously obtain accurate predictions on both monoalcohols and diols. The critical temperatures, critical densities, and normal boiling points of these three compounds are reported in the Supporting Information (Table S.2). A good agreement is obtained although these properties have not been used in the adjustment database. 3.3. Simulation of Primary, Secondary, and Tertiary Alcohols. Phase Equilibrium. Thermodynamic and structural properties of several primary alcohols (propan-1-ol, hexan-1ol, octan-1-ol), secondary alcohols (propan-2-ol), and tertiary

Transferable Force Field for Alcohols

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5989

Figure 3. Experimental16 (lines) and calculated (symbols) saturated liquid densities (a), vapor pressures (b), and vaporization enthalpies (c) of methanol (solid lines and squares), ethanol (dashed lines and triangles), and 1,3-propanediol (dotted lines and diamonds). Filled symbols denote data included in the adjustment database.

alcohols (2-methylpropan-2-ol) were predicted using the proposed force field. Figures 4-6 show experimental16 and calculated saturated liquid densities, vapor pressures and vaporization enthalpies of these compounds. Numerical values are reported in the Supporting Information (Table S.1). A good agreement between experimental and calculated values is achieved for saturated liquid densities. For 1-alcohols, average deviations are constant along the homologous series

and do not exceed on average 2.1% (maximum average deviation obtained for hexan-1-ol). This property is also well predicted for propan-2-ol, but slightly underestimated for 2-methylpropan-2-ol (average deviation of about 6%). The higher deviations for this property are obtained in the vicinity of the critical point. A good agreement is also obtained for vapor pressures. Here again, average deviations do not depend on the length of the alkyl chain for 1-alcohols, and do not exceed 13%

5990 J. Phys. Chem. B, Vol. 113, No. 17, 2009

Figure 4. Experimental16 (lines) and calculated (symbols) saturated liquid densities: (a) propan-1-ol (solid lines and squares), propan-2-ol (dashed lines and triangles), and octan-1-ol (dotted lines and diamonds); (b) hexan-1-ol (solid lines and squares) and 2-methylpropan-2-ol (dashed lines and triangles).

(maximum average deviation obtained for ethanol). Deviations are of the same order of magnitude as the ones obtained with the TraPPE force field8 for short alcohols, but significantly lower for longer alcohols. The higher deviations are obtained in the vicinity of the critical point for primary alcohols, but at low reduced temperatures for secondary and tertiary alcohols. Finally, the calculated vaporization enthalpies are in good agreement with experimental data. For each studied alcohol, the average deviations never exceed 11%. Experimental and calculated critical temperatures, critical densities and normal boiling temperatures are given in the Supporting Information (Table S.2). Deviations do not exceed 5% for the critical properties prediction and 2% for the normal boiling temperatures. It is particularly interesting to note that this new potential gives accurate predictions for propan-2-ol and 2-methylpropan2-ol saturated and critical properties although no additional united atom center was introduced neither for the CH group nor for the aliphatic C group located in R position to the hydroxyl group, these latter groups being directly taken from the AUA4 original model for hydrocarbons.2,24 This confirms the transferable feature of this potential. Structure of the Liquid Phase. Figure 7 compares the calculated and experimental carbon-carbon, oxygen-oxygen, oxygen-hydrogen, and carbon-oxygen radial distribution functions in methanol at 300 K. Experimental radial distribution functions of Yamaguchi et al.25 were obtained from neutron diffraction experiments. Although no structural properties were used in the adjustment of the hydroxyl group parameters, a very good agreement is found for each atom pair investigated. The first peak in the oxygen-oxygen radial distribution function is

Ferrando et al.

Figure 5. Experimental16 (lines) and calculated (symbols) vapor pressures. Same legend as in Figure 4.

located at 2.79 Å, due to intermolecular hydrogen bonds. This distance also agrees well with data obtained from X-ray diffraction experiments, which report a value of 2.7 Å for this peak.26,27 The second coordination shell peaks at about 5 Å. The angle formed by an oxygen atom and the two oxygen atoms in the first and second shells of the hydrogen bonds network was estimated to be about 125°. This value is particularly close to the experimental value of this angle in solid methanol28 (122°). The carbon-oxygen radial distribution function shows a first peak centered at about 3.53 Å in both simulated and experimental curves. The geometric criteria commonly used to define the existence of a hydrogen bond are the distance between two oxygen atoms O and O′ and the angle formed by the oxygen atom O, its bonded hydrogen H, and the second oxygen atom O′ (see Figure S.1 in the Supporting Information). In molecular simulation, the choice of the distance and angle values delimiting the existence of such a hydrogen bond is not uniquely defined.29 Figure 8 shows the number of hydrogen bonds per molecule as a function of the O-O′ distance and the O-H-O′ angle for methanol at 300 K. The hydrogen bond in methanol is strongly oriented with an angle ranging between 140 and 180°. The characterization of a hydrogen bond appears to be mostly sensitive to the distance criterion. Therefore, no angle criterion has been considered in this work. Only a distance criterion has been applied with an oxygen-oxygen distance ranging from 2.5 to 3.5 Å, which corresponds to the lower and upper limits of the first peak of the oxygen-oxygen radial distribution function (Figure 7). Using this criterion, an average hydrogen bonds number of 2.11 per methanol molecule is obtained, which is in good accordance with experiments and other UA force fields.4,8 The average number of hydrogen bonds at 300 K for the primary, secondary and tertiary alcohols studied in this work

Transferable Force Field for Alcohols

Figure 6. Experimental16 (lines) and calculated (symbols) vaporization enthalpies. Same legend as in Figure 4.

is given in the Supporting Information (Table S.3). This value is always close to 2 whatever the alcohol considered. The size of the aggregates formed by the hydrogen bond networks strongly depends on the nature of the alcohol (see Figure S.2 in the Supporting Information). The size of an aggregate is calculated by counting the number of hydrogen-bonded molecules using the previously defined distance criterion. In the case of propan-1-ol, the size distribution of aggregates is relatively uniform. On the contrary, for 2-methylpropan-2-ol, aggregates with a size of 4 molecules are predominant: about 10% of the molecules in the bulk belong to an aggregate of this size. 3.4. Simulation of Diols. In addition to the 1,3-propanediol molecule included in the adjustment database, 1,2-ethanediol and 1,5-pentanediol systems were simulated. Saturated liquid densities, vapor pressures, and vaporization enthalpies of these compounds are shown in Figure 9. Numerical values are given in the Supporting Information (Table S.1). Densities are slightly underestimated, but average deviations do not exceed 6%. Vapor pressures and vaporization enthalpies are also in good agreement with experimental data, and the higher deviations are systematically obtained in the vicinity of the critical point. The critical temperatures, critical densities, and normal boiling temperatures of these diols are given in the Supporting Information (Table S.2). Predictions are in good accordance with experimental data, although the critical temperature of 1,5-pentanediol is slightly overestimated. The overall agreement on the saturated and critical properties of diols demonstrates once again the transferability of the potential and validates the proposed method for calculating the intramolecular electrostatic energy. Histograms of the oxygen-oxygen intramolecular distances are shown in Figure 10 for the three diols at 300 K. For 1,2ethanediol, the distribution exhibits two distinct peaks: the first

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5991 one (2.9 Å) corresponds to a gauche conformation of the O-CH2-CH2-O dihedral angle, which favors the formation of an intramolecular hydrogen bond. About 22% of the molecules in the bulk are found in a gauche conformation. The second peak (3.6 Å) corresponds to a trans conformation of the molecule. For 1,3-propanediol and 1,5-pentanediol, the last peak of the distribution (4.8 and 7.4 Å, respectively) also corresponds to a geometry in which all the dihedral angles of the molecule are in trans conformation. Integration shows that in both cases about 10% of the molecules in the bulk adopt a trans conformation. Intermolecular and intramolecular hydrogen bonds in diols are determined with the geometric criterion defined above. The number of intermolecular and intramolecular hydrogen bonds per molecule at 300 K is given in the Supporting Information (Table S.3). All diols present about 4 hydrogen bonds per molecule, i.e., about 2 hydrogen bonds per hydroxyl group, which is consistent with the number of hydrogen bonds per molecule found for monoalcohols. The number of intramolecular hydrogen bonds decreases as the length of the alkyl chain increases, and no hydrogen bonds are anymore observed for the 1,5-pentanediol molecule. The effect of temperature on the average number of intermolecular and intramolecular hydrogen bonds per molecule of 1,2-ethanediol is illustrated in Table S.4 in the Supporting Information. As expected, more thermal agitation decreases the number of intermolecular hydrogen bonds. At the same time, the number of intramolecular hydrogen bonds per molecule slightly increases. These observations are consistent with the simulations performed by Stubbs et al.,9 using the TraPPE force field. 3.5. Simulation of Phenol. The transferability of the proposed force field was tested on phenol thermodynamic properties. Figure 11 shows the calculated saturated liquid densities, vapor pressures and vaporization enthalpies of phenol compared to experimental data. Numerical values are given in the Supporting Information (Table S.1). Densities and vaporization enthalpies are accurately predicted (deviations about 2%). Deviation on vapor pressures is only about 15% on average, which is particularly remarkable compared to predictions using other UA force fields where much larger deviations (about 30%) are obtained.11 It can also be noticed that these deviations are of the same order of magnitude as the deviations reported for the prediction of alkylbenzene vapor pressures with the original AUA4 potential for hydrocarbons. Experimental and calculated critical temperatures, critical densities, and normal boiling temperatures are reported in the Supporting Information (Table S.2). The critical and normal boiling temperatures are well predicted, but the critical density is slightly overestimated (deviation about 20%). Recently, it has been shown that adding some partial electrostatic charges in the original AUA4 potential to mimic the polarity induced by the aromatic ring allows to improve thermodynamic property predictions for benzene and methylbenzene.30,31 Therefore, some improvements could also be expected for phenol if such point charges are added on the aromatic ring. 3.6. Simulation of an Alkane + Alcohol Mixture. The prediction of phase equilibrium of mixtures containing both alcohols and hydrocarbons is today of particular interest in the context of biofuels. For example, ethanol can be directly added to gasoline blends in various proportions, and a good knowledge of phase properties and phase diagrams is then required for fuel formulation.32 The potential developed in this work aims to be used without any further modifications to simulate such complex mixtures. In the present work, the methanol + n-butane mixture

5992 J. Phys. Chem. B, Vol. 113, No. 17, 2009

Ferrando et al.

Figure 7. Experimental25 (dashed lines) and calculated (solid lines) carbon-carbon, oxygen-oxygen, oxygen-hydrogen, and carbon-oxygen radial distribution functions in methanol at 300 K and atmospheric pressure.

Figure 8. Number of hydrogen bonds per molecule of methanol at 300 K as a function of the oxygen (O)-oxygen (O′) distance, and the oxygen (O)-hydrogen (H)-oxygen (O′) angle.

was simulated, as a first step toward the study of these systems of industrial interest. During these simulations, the n-butane molecule was modeled using the AUA4 original force field.2 Experimental and simulated pressure-composition diagrams are shown in Figure 12. Numerical values are reported in Table 3. Experimental data are taken from the work of Leu and collaborators.33 Simulations are qualitatively and quantitatively in good agreement with experiments. The positive pressure azeotrope found in experiments is well reproduced by this force field. The azeotropic composition and pressure are accurately predicted, with a deviation of 8 and 3%, respectively. This system highlights the transferability of the force field to quantitatively predict the thermodynamics of mixtures.

4. Conclusion A new fully transferable potential for alcohols and polyalcohols has been developed, based on an extension of the AUA4 potential. Only one new pseudoatom for the hydroxyl group has been introduced. To ensure transferability to polyalcohols and complex molecules, the calculation of the intramolecular electrostatic interaction has been revisited. Instead of calculating the electrostatic interactions between charges related to nonbonded atoms, these interactions are computed between charges belonging to the different local dipoles of the molecule, one dipole being defined as a group of consecutive charges globally neutral. This new method allows avoiding the use of empirical scaling parameters commonly introduced to calculate the 1-4 electrostatic interactions. Consequently, in the proposed force

Transferable Force Field for Alcohols

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5993

Figure 9. Experimental16 (lines) and calculated (symbols) (a) saturated liquid densities, (b) vapor pressures, and (c) vaporization enthalpies of 1,2-ethanediol (solid lines and squares) and 1,5-pentanediol (dashed lines and triangles).

field, only three new parameters are required: the two LennardJones parameters and the AUA displacement of the hydroxyl group, which are adjusted on methanol, ethanol, and 1,3propanediol thermodynamic properties. The transferability of the potential is demonstrated since the introduction of these three parameters is sufficient to accurately simulate without any modifications a wide variety of alcohol families: primary alcohols (methanol, ethanol, propan-1-ol, hexan-1-ol, octan-1-

ol), secondary alcohols (propan-2-ol), tertiary alcohols (2methyl-propan-2-ol), phenol, and diols (1,2-ethanediol, 1,3propanediol, 1,5-pentanediol). Monte Carlo simulations carried out in the Gibbs ensemble have led to a good agreement between calculated and experimental data for the prediction of the thermodynamic properties along the liquid/vapor saturation curve, for the critical point coordinates, the normal boiling temperature, and the liquid structure at room temperature.

5994 J. Phys. Chem. B, Vol. 113, No. 17, 2009

Ferrando et al. TABLE 3: Calculated Pressure-Composition Data for the Methanol + n-Butane Mixture at 373.15 Ka P (MPa)

x methanol (mole fraction)

y methanol (mole fraction)

0.33417 0.806 1.000 1.216 1.500 1.630 1.76740 1.630 1.57840

1 0.9712 0.9509 0.91915 0.84120 0.55418 0.1577 0.0191 0

1 0.44420 0.34113 0.27914 0.2077 0.1876 0.1566 0.0563 0

a Italic data correspond to the azeotropic conditions. The subscript on calculated data gives the statistical uncertainty on the last decimal (for example, 0.44420 means 0.444 ( 0.020).

Figure 10. Histograms of the intramolecular oxygen-oxygen distances: 1,2-ethanediol (solid line), 1,3-propanediol (dashed line), and 1,5-pentanediol (dotted line).

Figure 12. Pressure-composition diagram for the methanol + n-butane mixture at 373.15 K. Experimental data (lines) are taken from Leu et al.33 GEMC simulations are depicted by triangles. The open square corresponds to the azeotrope calculation.

alcohol system has shown the ability of the proposed potential to reproduce the azeotropic behavior of such mixtures with a good qualitative and quantitative agreement. Work is now in progress to study the thermodynamic behavior of some more complex alcohols such as those belonging to the glycerol family. The overall transferability feature of the AUA4 potential gives the future opportunity to simulate multifunctional molecules and mixtures of industrial interest. The availability of such precise potentials allows the use of molecular simulation as a costefficient tool to generate thermodynamic data under a large range of pressure, temperature, and compositional conditions. Acknowledgment. The authors gratefully acknowledge B. Rousseau, P. Memari, and T. de Bruin for stimulating discussions. Supporting Information Available: Experimental and calculated values of saturated liquid densities, vapor pressures, and vaporization enthalpies are given in Table S.1. Critical points and normal boiling temperatures are reported in Table S.2. The number of hydrogen bonds per molecule is reported in Tables S.3 and S.4. Figure S.1 gives a schematic representation of a hydrogen bond between two molecules of methanol. Figure S.2 gives the fraction of molecules belonging to an aggregate of a given size for several alcohols. This material is available free of charge via the Internet at http://pubs.acs.org. Figure 11. Experimental16 (lines) and calculated (symbols) liquid densities, vapor pressures, and vaporization enthalpies of phenol.

Transferability is also an essential feature to simulate complex mixtures involving alcohols. The simulation of an alkane +

References and Notes (1) Huber, G. W.; Ibora, S.; Corma, A. Chem. ReV. 2006, 106, 4044– 4098. (2) Ungerer, P.; Tavitian, B.; Boutin, A. Applications of Molecular Simulation in the Oil and Gas Industry - Monte Carlo Methods; Technip: Paris, 2005.

Transferable Force Field for Alcohols (3) Toxvaerd, S. J. Chem. Phys. 1990, 93, 4290–4295. (4) Jorgensen, W. L. J. Phys. Chem. 1986, 90, 1276–1284. (5) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638–6646. (6) Van Leeuwen, M. E. Mol. Phys. 1996, 87, 87–101. (7) Smit, B.; Karaborni, S.; Siepmann, I. J. J. Chem. Phys. 1995, 105, 2126–2140. (8) Chen, B.; Potoff, J. J.; Siepmann, I. J. J. Phys. Chem. B 2001, 105, 3093–3104. (9) Stubbs, J. M.; Potoff, J. J.; Siepmann, I. J. J. Phys. Chem. B 2004, 108, 17596–17605. (10) Ungerer, P.; Beauvais, C.; Delhommelle, J.; Boutin, A.; Rousseau, B.; Fuchs, A. H. J. Chem. Phys. 2000, 112, 5499–5510. (11) Perez-Pellitero, J.; Bourasseau, E.; Demachy, I.; Ridard, J.; Ungerer, P.; Mackie, A. J. Phys. Chem. B 2008, 112, 9853–9863. (12) Ahunbay, M. G.; Ungerer, P.; Mackie, A.; Perez, J. J. Phys. Chem. B 2005, 109, 2970–2976. (13) 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–5197. (14) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. (15) Bourasseau, E.; Haboudou, M.; Boutin, A.; Fuchs, A. H.; Ungerer, P. J. Chem. Phys. 2003, 118, 3020–3034. (16) BYU DIPPR 8001, Thermophysical Properties Database Public release, January 2005. (17) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813–826. (18) Panagiotopoulos, A. Z. Mol. Sim. 1992, 9, 1–23.

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5995 (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: London, 1987. (20) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press, San Diego, CA, 1996. (21) Mackie, A. D.; Tavitian, B.; Boutin, A.; Fuchs, A. H. Mol. Simul. 1997, 19, 1–15. (22) Pitzer, K. S. J. Phys. Chem. 1995, 99, 13070–13077. (23) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527–545. (24) Bourasseau, E.; Ungerer, P.; Boutin, A.; Fuchs, A. Mol. Sim. 2002, 28, 317–336. (25) Yamaguchi, T.; Hidaka, K.; Soper, A. K. Mol. Phys. 1999, 97, 603– 605. (26) Harvey, G. G. J. Chem. Phys. 1938, 6, 111–114. (27) Wertz, D. L.; Kruh, R. K. J. Chem. Phys. 1967, 47, 388–390. (28) Jorgensen, W. L. J. Am. Chem. Soc. 1980, 102, 543–549. (29) Stubbs, J. M.; Siepmann, J. I. J. Phys. Chem. B 2002, 106, 3968– 3978. (30) Bonnaud, P.; Nieto-Draghi, C.; Ungerer, P. J. Phys. Chem. B 2007, 111, 3730–3743. (31) Nieto-Draghi, C.; Bonnaud, P.; Ungerer, P. J. Phys. Chem. C 2007, 111, 15686–15699. (32) French, R.; Malone, P. Fluid Phase Equilib. 2005, 228-229, 27– 40. (33) Leu, A. D.; Chen, C. J.; Robinson, D. B. AIChE Symp. Ser. 1989, 85, 11–16.

JP810915Z