Transferable Potentials for Phase Equilibria. 9 ... - ACS Publications

Aug 22, 2007 - nitrogen (two types), oxygen, and sulfur are transferable for all 13 compounds, ... TraPPE-EH force field are compared with high-level ...
15 downloads 0 Views 214KB Size
Download PDF
10790

J. Phys. Chem. B 2007, 111, 10790-10799

Transferable Potentials for Phase Equilibria. 9. Explicit Hydrogen Description of Benzene and Five-Membered and Six-Membered Heterocyclic Aromatic Compounds Neeraj Rai and J. Ilja Siepmann* Department of Chemistry and the Department of Chemical Engineering and Materials Science, UniVersity of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455 ReceiVed: May 10, 2007; In Final Form: June 26, 2007

The explicit hydrogen version of the transferable potentials for phase equilibria (TraPPE-EH) force field is extended to benzene, pyridine, pyrimidine, pyrazine, pyridazine, thiophene, furan, pyrrole, thiazole, oxazole, isoxazole, imidazole, and pyrazole. While the Lennard-Jones parameters for carbon, hydrogen (two types), nitrogen (two types), oxygen, and sulfur are transferable for all 13 compounds, the partial charges are specific for each compound. The benzene dimer energies for sandwich, T-shape, and parallel-displaced configurations obtained for the TraPPE-EH force field compare favorably with high-level electronic structure calculations. Gibbs ensemble Monte Carlo simulations were carried out to compute the single-component vapor-liquid equilibria for benzene, pyridine, three diazenes, and eight five-membered heterocycles. The agreement with experimental data is excellent with the liquid densities and vapor pressures reproduced within 1 and 5%, respectively. The critical temperatures and normal boiling points are predicted with mean deviations of 0.8 and 1.6%, respectively.

1. Introduction Aromatic molecules are an important class of organic compounds found in many biological assemblies and pharmaceutical compounds. The nonbonded interactions present in the aromatic compounds play a key role in the vertical base stacking of DNA double strands, drug-DNA binding, supramolecular chemistry, and crystal engineering applications, such as molecular devices and molecular recognition.1-12 As these interactions are relatively weak, it is often difficult to quantitatively determine them experimentally.13-15 Hence, numerous attempts using both theory4,16-22 and simulation23-29 have been made to understand and quantify aromatic interactions. High-level ab initio quantum mechanical calculations are computationally prohibitive for condensed phases and larger systems thereby limiting its use to smaller systems such as benzene dimers, trimers, and tetramers.20,21 Molecular simulations provide an alternative way of determining the intermolecular potentials. These empirical potentials (or force fields) are generally obtained by fitting Lennard-Jones parameters and, sometimes, partial atomic charges to experimental data such as the liquid density and heat of vaporization. Most of the general force fields, such as AMBER,23 OPLS,24,25 and CHARMM,26 are fitted to the experimental data at ambient conditions, while the anisotropic OPPE force field27,28 and the exp-6 force field of Errington and Panagiotopoulos29 are limited to benzene and naphthenic systems. Hence, a force field that can accurately reproduce the thermophysical properties of heterocyclic aromatic molecules over a wide range of temperatures and pressures is desirable for many applications. An ongoing effort of this research group has been to develop force fields that are simple yet accurate and transferable. These force fields, known as transferable potentials for phase equilibria * Corresponding author. Email: [email protected]

(TraPPE), have been been developed for linear and branched alkanes,30-32 alkenes,33 alcohols,34 ethers,35 aldehydes,35 ketones,35 acids,36 esters,37 amines,38 amides,38 nitroalkanes,38 sulfides,39 disulfides,39 thiols,39 and small molecules like O2,40 N2,41 CO2,41 water,42 and ammonia.43 Recently, the TraPPE force field has been used successfully to predict bubble point pressures for ethanol and 1,1,1,2,3,3,3-heptafluoropropane44 and viscosities of polyhydric alcohols45 as part of the third industrial fluid properties simulation challenge.46 The TraPPE force field has been used extensively to investigate complex chemical systems, such as retention in chromatography,47-49 octanol-water partitioning,50 adsorption in pharmaceutical solids,51 and transport properties.52-54 For alkanes and nitriles, the TraPPE force field provides a choice of using either the united atom (UA) or the explicit hydrogen (EH) representation of CHx groups. The UA version of the force field is simple and results in savings of computer time, while the EH version provides more accurate liquid and vapor densities, vapor pressures, and heats of vaporization over a wide range of temperatures and pressures at a higher computational cost. In the present work, the TraPPE-EH force field is extended to benzene, pyridine, pyrimidine, pyrazine, pyridazine, thiophene, furan, pyrrole, thiazole, oxazole, isoxazole, imidazole, and pyrazole (see Figure 1). The remainder of the article is organized as follows. The details of the force field and parametrization philosophy are discussed in section 2. The simulation details are provided in the section 3. In section 4, potential energy curves for benzene dimers and the vapor-liquid coexistence curves for benzene and aromatic heterocycles obtained with the TraPPE-EH force field are compared with high-level ab initio calculations and experimental data, respectively.

10.1021/jp073586l CCC: $37.00 © 2007 American Chemical Society Published on Web 08/22/2007

Transferable Potentials for Phase Equilibria

J. Phys. Chem. B, Vol. 111, No. 36, 2007 10791

Figure 1. Name and structures of the molecules.

TABLE 1: TraPPE-EH Lennard-Jones Parameters for Aromatic Heterocyclesa atom type

/kB (K)

σ (Å)

X(aro)-C(aro)-X(aro) H-C(aro) H-N(aro) C(aro)-N(aro)-X(aro) H-N(aro)-X2(aro) C(aro)-O(aro)-X(aro) C(aro)-S(aro)-X(aro)

30.70 25.45 12.00 57.00 141.00 70.00 158.00

3.60 2.36 0.50 3.20 3.40 2.60 3.55

a Bold and regular fonts indicate the specific atom and its bonding environment, respectively. X(aro) indicates a neighboring C(aro), N(aro), or O(aro) atom.

2. Force-Field Development The TraPPE force field uses simple Lennard-Jones (LJ) and Coulomb potentials to model nonbonded interactions

U(rij) ) 4ij

[( ) ( ) ] σij rij

12

-

σij rij

6

+

qiqj 4π0rij

(1)

where rij, σij, ij,, qi, qj, and 0 are the separation between two interacting sites, LJ size, LJ well depth, partial atomic charge on sites i and j, and the permittivity of vacuum, respectively. The Lorentz-Berthelot55,56 combining rules are used to determine LJ parameters for unlike interaction sites. Tables 1 and 2 list the LJ parameters and partial atomic charges for compounds studied in this work. Aromatic rings are modeled here as rigid entities; hence, there are no intramolecular interactions. A previous investigation for benzene has shown that allowing for angle bending and dihedral flexibility has a negligible effect on vapor-liquid equilibria.57 However, if flexibility is desirable for algorithmic reasons (e.g., to avoid the use of constraints in molecular dynamics simulations), then the harmonic force constants of other force fields can be used. The bond lengths and bond angles for the aromatic heterocyles are listed in Tables 3 and 4, respectively. For the aromatic heterocycles, the TraPPE-EH interaction sites coincide with the atomic centers for all atom types, whereas

the interaction sites are centered on the carbon atoms and on the C-H bond centers for alkanes.32 The reason for this change is that first-order electrostatic interactions are important for the aromatic heterocyles and partial charges are located on all atomic sites (whereas partial charges are not included for the alkanes). Placing partial charges on atomic sites is significantly more convenient when those charges are determined from electronic structure calculations. Additionally, when only atomic sites are used it is easier to implement the TraPPE-EH force field in molecular dynamics codes because the calculation of torques is not required. There are numerous models and algorithms available to assign partial charges to the atoms in a molecule. These algorithms generally fall into two categories: those that use a least-squares fit to the molecular electrostatic potential to determine partial charges, for example, the Merz-Singh-Kollman scheme58,59 and ChElPG method,60 and another group that uses population analysis to determine partial atomic charges, for example, Mulliken,61-63 Lowdin,64 and the charge models developed by Cramer, Truhlar, and co-workers.65-68 In the present work, we have used the CM4 charge model68 to determine the partial atomic charges for the aromatic heterocycles. This decision was influenced mainly by the need to get accurate charges in the condensed phase. Kelly et al.67 demonstrated that the CM4 model reproduces dipole moments and solvation free energies for a large set of molecules more accurately, while being less sensitive to the basis set employed in the calculations than the ChElPG method. As the CM4 model is parametrized for the use of gas-phase optimized geometries, these can be determined at a high level of theory followed by the computation of partial charges in solution. The philosophy followed in the development of the TraPPE force field is to minimize the number of different parameters required to model the homologous series of a particular functional group, hence, ensuring the ability to build larger molecules that have similar functional groups. In the present work, the same parametrization philosophy is used to determine the LJ parameters for aromatic heterocycles. However, because the charge density in aromatic compounds is much more diffuse

10792 J. Phys. Chem. B, Vol. 111, No. 36, 2007

Rai and Siepmann

TABLE 2: TraPPE-EH Partial Charges for Aromatic Heterocyclesa molecule benzene pyridine

q (e)

N1

-0.416

pyrimidine

N1, N3

-0.439

pyrazine pyridazine

N1, N4 N1, N2

-0.388 -0.290

thiophene

S1

-0.028

furan

O1

-0.192

pyrrole

N1

-0.420

thiazole

S1 N3

+0.022 -0.444

oxazole

O1 N3

-0.170 -0.445

isoxazole

O1 N2

-0.117 -0.279

imidazole

N1 N3

-0.416 -0.485

N1 N2

-0.302 -0.324

pyrazole

a

hetero

carbon

q (e)

hydrogen

q (e)

C C2, C6 C3, C5 C4 C2 C4, C6 C5 C2, C3, C5, C6 C3, C6 C4, C5 C2, C3, C4, C5

-0.95 +0.093 -0.101 -0.037 +0.284 +0.145 -0.099 +0.095 +0.119 -0.050 -0.107

C2, C5 C3, C4

+0.016 -0.135

C2, C5 C3, C4 C2 C4 C5 C2 C4 C5 C3 C4 C5

+0.019 -0.166 -0.102 +0.057 +0.123 +0.226 +0.036 +0.019 +0.090 -0.144 +0.082

C2 C4 C5

+0.224 +0.005 +0.030

C3 C4 C5

+0.060 -0.170 +0.081

H H7, H11 H8, H10 H9 H7 H8, H10 H9 H7, H8, H9, H10 H7, H10 H8, H9 H6, H9 H7, H8 H6, H9 H7, H8 H6 H7, H10 H8, H9 H6 H7 H8 H6 H7 H8 H6 H7 H8 H6 H7 H8 H9 H6 H7 H8 H9

+0.95 +0.076 +0.105 +0.107 +0.079 +0.100 +0.124 +0.099 +0.124 +0.097 +0.121 +0.107 +0.101 +0.114 +0.322 +0.098 +0.098 +0.136 +0.096 +0.114 +0.102 +0.110 +0.122 +0.119 +0.130 +0.119 +0.336 +0.097 +0.092 +0.117 +0.320 +0.098 +0.117 +0.120

Refer to Figure 1 for atom numbering.

TABLE 3: TraPPE-EH Bond Lengths (Å) for Aromatic Heterocyclesa bond type

molecule

length

bond type

molecule

length

bond type

molecule

length

C-H N-H C-C N-C C-C C-N C-C C-C C-N C-N C-C N-N

all all benzene pyridine pyridine pyrimidine pyrimidine pyrazine pyrazine pyridazine pyridazine pyridazine

1.080 1.000 1.392 1.330 1.392 1.360 1.392 1.392 1.340 1.340 1.392 1.350

S1-C2 C2-C3 C3-C4 O1-C2 C2-C3 C3-C4 N1-C2 C2-C3 C3-C4 S1-C2 S1-C5 C2-N3 N3-C4 C4-C5

thiophene thiophene thiophene furan furan furan pyrrole pyrrole pyrrole thiazole thiazole thiazole thiazole thiazole

1.71 1.37 1.42 1.36 1.36 1.43 1.37 1.38 1.42 1.72 1.72 1.28 1.38 1.34

O1-C2 C2-N3 N3-C4 C4-C5 C5-O1 O1-N2 N2-C3 C3-C4 C4-C5 C5-O1 N1-C2 N1-H6 C2-N3 N3-C4 C4-C5 C5-N1 N1-N2 N1-H5 N2-C3 C3-C4 C4-C5 C5-N1

oxazole oxazole oxazole oxazole oxazole isoxazole isoxazole isoxazole isoxazole isoxazole imidazole imidazole imidazole imidazole imidazole imidazole pyrazole pyrazole pyrazole pyrazole pyrazole pyrazole

1.35 1.28 1.38 1.35 1.36 1.38 1.30 1.42 1.35 1.33 1.36 1.00 1.31 1.37 1.37 1.37 1.34 1.00 1.33 1.41 1.38 1.35

a

Refer to Figure 1 for atom numbering.

and delocalized in comparison to aliphatic hydrocarbons, the effect of substituent heteroatoms is not confined to only the R-carbon (nearest neighbor) atoms but also affects the charge density of the entire aromatic ring. This leads to nontransferability of charges for aromatic heterocyles. In order to determine the partial charges, the CM4 charge model with the same level of theory and basis set (see below) was used for all the molecules. 1-Octanol was used as the universal continuum

solvent for determining partial charges of all aromatic compounds, as this solvent possesses both polar and nonpolar character. Once the CM4 charges were found for all atomic sites, the LJ parameters were fitted to the saturated liquid densities and vapor pressures of a few selected compounds. For determination of the LJ parameters for aromatic compounds, a stepwise procedure was used. First, the common LJ parameters for C and H were obtained by fitting to the vapor-

Transferable Potentials for Phase Equilibria

J. Phys. Chem. B, Vol. 111, No. 36, 2007 10793

TABLE 4: TraPPE-EH Bond Angles (deg) for Aromatic Heterocyclesa bond type

molecule

angle

bend type

molecule

angle

bond type

molecule

angle

C-C-C H-C-H C-C-H C2-N1-C6 N1-C2-C3 C3-C4-C5 N1-C2-H7 H8-C3-C4 H9-C4-C3 N1-C2-N3 C2-N3-C4 N3-C4-C5 C4-C5-C6 H7-C2-N3 N3-C4-H8 C4-C5-H9 C2-N1-C6 N1-C2-C3 N1-C2-H7 H7-C2-C3 N1-N2-C3 N2-C3-C4 C3-C4-C5 N2-C3-H7 C3-C4-H8

benzene benzene benzene pyridine pyridine pyridine pyridine pyridine pyridine pyrimidine pyrimidine pyrimidine pyrimidine pyrimidine pyrimidine pyrimidine pyrazine pyrazine pyrazine pyrazine pyridazine pyridazine pyridazine pyridazine pyridazine

120.0 120.0 120.0 117.3 123.6 118.5 116.1 121.3 120.7 126.9 115.8 122.3 116.8 116.6 116.3 121.6 115.3 122.3 116.7 121.0 119.0 124.1 116.9 114.5 120.9

S1-C2-C3 C2-S1-C5 C2-C3-C4 S1-C2-H6 C2-C3-H7 C2-O1-C5 O1-C2-C3 C2-C3-C4 O1-C2-H6 C2-C3-H7 C2-N1-C5 N1-C2-C3 C2-C3-C4 H6-N1-C2 N1-C2-H7 C2-C3-H8 C2-S1-C5 S1-C2-N3 C2-N3-C4 N3-C4-C5 C4-C5-S1 C4-C5-H8 S1-C2-H6 C5-C4-H7

thiophene thiophene thiophene thiophene thiophene furan furan furan furan furan pyrrole pyrrole pyrrole pyrrole pyrrole pyrrole thiazole thiazole thiazole thiazole thiazole thiazole thiazole thiazole

111.5 92.2 112.4 120.1 123.3 106.9 110.5 106.0 115.9 126.9 110.0 107.6 107.4 125.0 121.2 125.7 88.6 115.4 110.6 115.8 109.6 128.2 121.0 125.3

C5-O1-C2 O1-C2-N3 C2-N3-C4 N3-C4-C5 C4-C5-O1 O1-C2-H6 C4-C5-H8 C5-C4-H7 O1-N2-C3 N2-C3-C4 C3-C4-C5 C4-C5-O1 C5-O1-N2 C4-C3-H6 C3-C4-H7 C4-C5-H8 N1-C2-N3 C2-N3-C4 N3-C4-C5 C4-C5-N1 C5-N1-C2 C5-N1-H5 N1-C2-H7 C5-C4-H8 C4-C5-H9 N1-N2-C3 N2-C3-C4 C3-C4-C5 C4-C5-C6 C5-N1-N2 C5-N1-H6 C4-C3-H7 C3-C4-H8 C4-C5-H9

oxazole oxazole oxazole oxazole oxazole oxazole oxazole oxazole isoxazole isoxazole isoxazole isoxazole isoxazole isoxazole isoxazole isoxazole imidazole imidazole imidazole imidazole imidazole imidazole imidazole imidazole imidazole pyrazole pyrazole pyrazole pyrazole pyrazole pyrazole pyrazole pyrazole pyrazole

104.4 114.6 104.2 109.0 107.8 117.0 135.1 129.0 105.5 112.1 102.9 110.4 109.1 129.0 128.8 133.5 111.6 105.3 110.7 105.0 107.4 126.3 122.4 127.9 132.6 104.2 111.9 104.4 106.1 113.4 127.6 128.7 128.3 132.1

a

Refer to Figure 1 for atom numbering

liquid coexistence curve of neat benzene (followed by a check that this model also yields a satisfactory normal melting temperature69). These LJ parameters were then transferred to model C and H (nonpolar) in the five-membered and sixmembered rings. Thus, the only parameters to be determined in the next step were those of the heteroatoms, namely, nitrogen, oxygen, sulfur, and polar hydrogen. There are two different types of nitrogen atoms present in the aromatic systems, basic (in pyridine and diazenes) and acidic (in pyrrole). The LJ parameters for the polar hydrogen attached to an acidic nitrogen were set to σ ) 0.5 Å and /kB ) 12 K. The LJ parameters for the basic and acidic nitrogens were obtained by fitting to the saturated liquid densities and saturated vapor pressures for pyridine and pyrrole, respectively. Thiophene and furan were used to obtain the LJ parameters for sulfur and oxygen, respectively. These LJ parameter were then transferred to model diazenes, oxazoles, thiazole, imidazole, and pyrazole without the need for additional LJ parameters. 3. Simulation Details The determination of the partial atomic charges for the TraPPE-EH model follows a two-step procedure. First, the gasphase structure for each molecule was optimized using KohnSham density functional theory at the B3P86/6-311+G(3df,3pd) level of theory/basis set70,71 as implemented in the Gaussian 03 electronic structure software package.72 Second, the optimized gas-phase geometry was used in the continuum solvation calculations with 1-octanol as the universal solvent. The continuum solvation calculation employed the B3P86/6-31+G** level of theory/basis set70,71 as implemented in the MinnesotaGaussian solvation model (MN-GSM) version 6 module.73

Gibbs ensemble Monte Carlo (GEMC) simulations74,75 were carried out to compute the vapor-liquid coexistence curves (VLCC) for pure compounds. The system size consisted of 250 molecules. A spherical cutoff, rcut, at 14 Å and analytical tail corrections76 were used to compute the LJ interactions. The Ewald summation technique with tin foil boundary conditions76,77 was employed to compute the Coulomb interactions. A cutoff, rcutchg, at 14 Å was used for the real-space part, the Ewald sum convergence parameter, κ, was set to 3.2/rcutchg, and the maximum number of reciprocal space vectors included in the sum, Kmax, was set to 8. The canonical version of the GEMC approach74 employs two periodic simulation boxes in thermodynamic contact but without an interface. Molecules are swapped between the boxes with the help of the configurational-bias Monte Carlo approach78,79 to equilibrate the chemical potential, volume exchange moves are employed to equilibrate the pressure, and translational and rotational moves are used to thermally equilibrate the system. The liquid-phase boxes were larger than twice the cutoff used for LJ potential and the real-space part of the Ewald summation. The size of the vapor box was adjusted during the equilibration period so that on average approximately 20-40 molecules were in the vapor phase. The fraction of volume moves was set to have approximately 1 volume move accepted every 10 Monte Carlo (MC) cycles (one MC cycle consists of N moves, where N ) 250 is the number of molecules in the system). Similarly, the fraction of particle swap moves was adjusted to yield approximately one accepted swap move per 10-50 MC cycles (with the higher values used for lower temperatures). The remaining fraction of MC moves was equally divided between the translational and the rotational moves. At any state point,

10794 J. Phys. Chem. B, Vol. 111, No. 36, 2007

Figure 2. Schematic representations of sandwich, T-shape, and paralleldisplaced configurations for the benzene dimer.

20 000 MC cycles were used to equilibrate the system before starting the production runs. The production runs consisted of at least 40 000 MC cycles and were divided into five blocks to compute statistical uncertainties (the standard deviations). The critical temperatures, Tc, and the critical densities, Fc, were estimated using the saturated density scaling law and the law of rectilinear diameter with a scaling exponent, β*, of 0.325.80-82 The normal boiling points, Tb, were computed via the Clausius-Clapeyron equation83 for the three temperatures surrounding Tb. 4. Results and Discussion 4.1. Benzene Dimer. Although not used in the parametrization of the TraPPE force field, it is instructive to compare the minimum binding energies for benzene dimers with specific relative configurations obtained for a set of empirical force fields to high-level electronic structure calculations. Sinnokrot and Sherrill19 used the estimated coupled-cluster with singles, doubles, and perturbative triples (CCSD(T)) theory with an augcc-pVQZ* basis set (a correlation-consistent split valence basis set of contracted Gaussian functions,84 specifically, the aug-ccpVQZ basis set minus the g functions on carbons and the f functions on hydrogen atoms) and the recommended frozen monomer geometry of Gauss and Stanton85 (rC-C ) 1.3915 Å and rC-H ) 1.080 Å) to search for minima on the potential surface of the benzene dimer. These authors found configurations for three distinct minima for the benzene dimer (see Figure 2) that are called sandwich (S), T-shape (T), and paralleldisplaced (PD) configurations. For the S and T configurations, the energy depends only on the center-of-mass separation, r, whereas vertical and parallel displacements, r1 and r2, are needed to characterize the PD configuration. The binding energies and structural parameters for the benzene dimers in the three configurations computed with three different ab initio levels (MP2/aug-cc-pVQZ*,19 CCSD(T)/augcc-pVQZ*,19 and the estimated CCSD(T) complete basis set (CBS) limit18) and for six empirical potentials (TraPPE-EH, 6-site TraPPE-UA,33 9-site TraPPE-UA,86 12-site OPLS-AA,87 6-site OPPE-AUA,27 and 6-site EP/exp-629) are compared in Table 5. The binding energies for the ab initio energies were taken from the literature, whereas those for the empirical force fields were calculated for this work. The two higher-level ab initio methods predict that the PD dimer is most stable with a binding energy (for the CBS CCSD(T) method) of -11.64 kJ/ mol and is closely followed by the T dimer with a value of -11.47 kJ/mol, whereas the binding energy for the S dimer is about 4 kJ/mol smaller in magnitude. The MP2 level also yields the PD dimer as the most favorable, but in this case the binding energy is nearly -20 kJ/mol and the T and S dimers yield -15 and -14 kJ/mol, respectively. Thus, we focus on the results for the CBS CCSD(T) calculations, the best available theoretical estimate,18 for judging the performance of the empirical force fields.

Rai and Siepmann The potential energy curves obtained for the three dimer configurations with the TraPPE-EH model (this work), the OPLS-AA model (another 12-site model with partial charges at atomic sites87), the TraPPE-UA-9 model (a united atom version with three additional point charges to represent the quadrupole86), the TraPPE-UA-6 model (a six-site united atom version33), the EP/exp-6 model (a six-site model using the exp-6 potential that was developed by Errington and Panagiotopoulos29), and the OPPE-AUA model (a six-site anisotropic united atom model27) are compared numerically in Table 5 and graphically in Figures 3 and 4. The three six-sites models (TraPPE-UA-6, EP/exp-6, and OPPE-AUA) qualitatively fail to reproduce the dimer energetics; that is, they substantially overestimate Udimer for the S configuration and underestimate it for the T configuration. The two six-site models that place interaction sites at the carbon locations (TraPPE-UA-6 and EP/exp-6) yield the S dimer as the global minimum (by definition, this leads to an equivalent PD dimer with r2 ) 0 Å). Although the OPPE-AUA model yields the PD dimer as the global mimimum, this is only slightly more favorable than the S dimer. The OPPE-AUA model places the interaction sites 0.4071 Å away from the carbon atom on the carbon-hydrogen bond (i.e., its effective “bond” length between interaction sites is 1.8071 Å compared to a value of 1.40 for the other six-site models), and it consequently uses a smaller σ value of 3.2464 Å compared to about 3.70 for the other six-site models. The resulting ratio of σ/l ) 1.80 (vs 2.64) leads to a more corrugated potential of the OPPE-AUA model for parallel dimer configurations and, hence, for a slight favoring of the PD dimer for r1 < 21/6σ. Without this corrugation (and any quadrupolar interactions), the TraPPE-UA-6 and EP/exp-6 models yield r1 distances that are too large for the PD dimer, whereas the OPPE-AUA model underestimates the separation for the S dimer. Furthermore, the EP/exp-6 and OPPE-AUA models yield binding energies for the S and PD configurations that are substantially too large. The united atom model with explicit quadrupolar interactions (TraPPE-UA-9 model using six Lennard-Jones sites and three charges) yields significantly better energetics for the PD and T configurations. Its global minimum is the PD dimer but with a significantly overestimated r2 value. On the negative side, it does not yield a bound dimer for the S configuration (see Figure 3). The two 12-site models seem to perform significantly better than the six- and nine-site united atom models. Their dimer energetics are much closer to the best theoretical estimate (with mean unsigned errors of about 2 kJ/mol). However, the relative rankings are incorrect with PD > S > T (for TraPPE-EH) and T > PD > S (for OPLS-AA) versus PD > T > S (estimated for CBS CCSD(T)18). While the TraPPE-EH model yields satisfactory dimer structures, the r2 is significantly overestimated for the OPLS-AA model (see Figure 4). The relatively large deviations for the T configurations observed for the TraPPEEH model may be a consequence of using LJ potentials and ignoring induced polarization which is most important for this dimer configuration. Table 5 also lists Udimer computed for the six empirical force fields with the structures fixed at the geometric parameters obtained from the CBS CCSD(T) calculations. Such a comparison clearly penalizes any force field that overestimates the center-of-mass separation for the minimum-energy structure relative to one that underestimates the separation because of the steeply repulsive nature of intermolecular interactions at short separation (see also Figure 3). Nevertheless, the outcome

Transferable Potentials for Phase Equilibria

J. Phys. Chem. B, Vol. 111, No. 36, 2007 10795

TABLE 5: Binding Energies, Udimer (kJ/mol), and Intermolecular Distances, r, r1, and r2 (Å), for Different Configurations of the Benzene Dimera S

T

MUPEb

PD

method/force field

Udimer

r

Udimer

r

Udimer

r1

r2

/aug-cc-pVQZ*18

MP2 CCSD(T)/aug-cc-pVQZ*19 CBS CCSD(T) 18 TraPPE-EH OPLS-AA87 TraPPE-UA-986 TraPPE-UA-633 EP/exp-629 OPPE-AUA27

-14.11 -7.12 -7.58 -9.29c -7.15c

3.7 3.9 3.7 3.75 3.77

-12.71 -17.23 -16.61c

3.92 3.88 3.44

-14.82 -10.93 -11.47 -8.34 -9.02 -8.91 -6.99 -8.89 -9.60

4.9 5.0 4.9 5.20 5.10 5.03 5.07 5.06 4.80

-20.05 -10.98 -11.64 -9.84 -8.72 -10.53 -12.71 -17.23 -16.77

3.4 3.6 3.4 3.61 3.47 3.72 3.92 3.88 3.39

1.6 1.6 1.6 1.75 2.84 2.58 0.00 0.00 0.81

TraPPE-EH OPLS-AA87 TraPPE-UA-986 TraPPE-UA-633 EP/exp-629 OPPE-AUA27

-9.24 -7.05 +24.64 -11.01 -14.80 -14.90

3.7 3.7 3.7 3.7 3.7 3.7

-5.48 -7.88 -8.56 -6.51 -8.01 -9.42

4.9 4.9 4.9 4.9 4.9 4.9

-8.29 -6.78 -4.48 -2.59 +2.28 -15.72

3.4 3.4 3.4 3.4 3.4 3.4

1.6 1.6 1.6 1.6 1.6 1.6

Udimer

r

22 18

6 22

39 66 60

31 31 15

34 27 171 55 82 50

a The top part gives the energies and geometric parameters for the minimum-energy structures obtained with a given ab initio method/force field for a particular relative orientation. The global minimum is highlighted in bold font. The bottom part gives the energies computed with the geometric parameters fixed at the minimum-energy values found for CBS CCSDT. bMean unsigned percentage error. c Saddle point with respect to parallel displacement.

Figure 3. Potential energy curves for the sandwich and T-shaped configurations of the benzene dimer. The black, red, blue, green, violet, brown, and orange lines represent data obtained from the estimated CCSD(T)/aug-cc-pVQZ* level of theory19 and the TraPPE-EH, TraPPEUA-9, TraPPE-UA-6, OPLS-AA, EP/exp-6, and OPPE-AUA force fields, respectively.

of this comparison with fixed separations is the same, that is, the two 12-site models are significantly more accurate than the six-site models (and the nine-site model remains unsatisfactory for the S configuration). Mean unsigned percentage errors (relative to the CBS CCSD(T) values) for the binding energies and geometric parameters

Figure 4. Potential energy curves for the parallel-displaced configuration of the benzene dimer (top to bottom: r1 ) 3.2, 3.4, 3.6, and 3.8 Å). The line colors are the same as those used in Figure 3.

obtained for the six empirical models are also listed in Table 5. When the mean unassigned percentage errors (MUPE) are averaged for the six energies (three relative orientations at the

10796 J. Phys. Chem. B, Vol. 111, No. 36, 2007

Rai and Siepmann

TABLE 6: Comparison of the Dipole Moments (D) for Aromatic Heterocycles in the Gas Phase and in 1-Octanol Solvent gas phase molecule benzene pyridine pyrimidine pyrazine pyridazine thiophene furan pyrrole thiazole oxazole isoxazole imidazole pyrazole

solvent

CM4

expt90

CM4

2.02 2.22

2.15 2.42

2.77 3.08

3.77 0.56 0.74 1.21 1.56 1.45 2.54 3.05 1.97

4.10 0.52 0.72 1.74 1.60 1.50 2.90 3.80 2.21

5.15 0.72 0.93 1.63 2.07 1.87 3.28 3.96 2.61

minimum-energy structures for either the force field or the CBS CCSD(T) method) and four separations, the empirical models rank as TraPPE-EH (MUPE ) 19%), OPLS-AA (22%), OPPEAUA (39%), TraPPE-UA-6 (41%), and EP/exp-6 (57%). Although a failure of a given model to yield the correct energetics for the three dimer configurations can be well masked in the calculation of VLCC (see below), it is well-known that quadrupolar interactions are essential to stabilize the molecular orientations found in benzene’s low-pressure, low-temperature orthorombic phase and to yield the correct volume increment upon melting.88,89 Indeed, the new TraPPE-EH model yields an excellent normal melting temperature of 276 ( 8 K.69 4.2. Dipole Moments. The dipole moments of the aromatic heterocycles in the gas phase and in 1-octanol solvent are summarized in Table 6. Compared to the experimental gas-phase values,90 the CM4 partial atomic charges68 yield dipole moments that are mostly underpredicted except for thiophene and furan that are slightly overpredicted by 0.04 and 0.02 D, respectively. Although the computed gas-phase dipole moments for pyridine, pyrimidine, pyridazine, thiazole, oxazole, isoxazole, and pyrazole all fall within 12% of the experimental values, those for imidazole and pyrrole are underpredicted by 20 and 30%, respectively. In contrast, the computed solvent-phase dipole moments are (with the exception of pyrrole) larger than the experimental gas-phase dipole moments. Excluding imidazole and pyrrole, the solution-phase dipole moment is on average 26% larger than the experimental value (including these two compounds, the average enhancement decreases to 21%). This illustrates the importance of accounting for condensed-phase polarization effects in a systematic fashion as is done in this work through the use of 1-octanol as the universal continuum solvent for the TraPPE parametrization of partial charges. 4.3. Vapor-Liquid Equilibria. The vapor-liquid coexistence curves (VLCC) and the Clausius-Clapeyron representation of the saturated vapor pressure lines for different benzene force fields are shown in Figure 5, and the critical constants and the normal boiling point for the TraPPE-EH force field are listed in Table 7. With the exception of the OPLS-AA model, all force fields perform remarkably well for the VLCC of benzene. This should not come as a surprise because these force fields were fitted to benzene’s VLCC. For the TraPPE-EH force field, the saturated liquid densities lie within 1% of the experimental values (except for the state point closest to the critical point) and the average error for the saturated vapor pressures is less than 5%; the latter is a significant improvement over the TraPPE-UA-6 and TraPPE-UA-9 models. The critical temperature, critical density, and normal boiling point for the TraPPE-

Figure 5. Coexistence densities (left) and saturated vapor pressures (right) for benzene. The black line and star represent the experimental densities and critical point, respectively.91 Simulation data for the TraPPE-EH, TraPPE-UA-9, TraPPE-UA-6, OPLS-AA, EP/exp-6, and OPPE-AUA force fields are shown as red circles, blue squares, green diamonds, violet up triangles, brown left triangles, and orange right triangles, respectively.

TABLE 7: Comparison of the Critical Constants, Normal Boiling Points, and Heats of Vaporization Predicted by the TraPPE-EH Model with Experimental Data91 compound benzene

method

TraPPE-EH Expt. pyridine TraPPE-EH Expt. pyrimidine TraPPE-EH Expt. pyrazine TraPPE-EH Expt. pyridazine TraPPE-EH Expt. thiophene TraPPE-EH Expt. furan TraPPE-EH Expt. pyrrole TraPPE-EH Expt. thiazole TraPPE-EH Expt. oxazole TraPPE-EH Expt. isoxazole TraPPE-EH Expt. imidazole TraPPE-EH Expt. pyrazole TraPPE-EH Expt. MUPE TraPPE-EH

Tc (K) Fc (g/cm3) Tb (K) ∆Hvapa (kJ/mol) 5634 562.2 6194 620b 6322 n/a 6164 627c 7614 n/a 5874 580 4942 490 6512 640 6363 n/a 5392 n/a 5913 590 8183 n/a 7373 n/a 0.8%

0.2986 0.304 0.3194 0.32b 0.3342 n/a 0.3466 0.35c 0.3243 n/a 0.3596 0.38 0.3373 0.31 0.3184 0.33 0.3924 n/a 0.3814 n/a 0.3586 0.36 0.3293 n/a 0.3412 n/a 3.1%

3561 353.3 3913 388.6 3973 396.7 3903 388.0 5044 481.0 3591 357.3 3052 304.7 4042 402.9 3844 390.7 3252 342.7 3661 368.7 5132 530.2 4473 460.2 1.6%

35.0 33.9 43.0 40.2 44.0 49.9 39.9 60.3 53.5 35.3 34.8 27.4 27.7 43.9 45.4 41.1 32.1 32.5 38.3 37.2 57.6 n/a 50.6 n/a 4.9%

a Heats of vaporization are reported for a temperature close to 298.15 K. b Taken from ref 92 c Taken from ref 93

EH model deviate from experiments by only 0.2, 2.0, and 0.8%. The OPPE-AUA model predicts the saturated liquid densities with an error of 1% and the saturated vapor pressures within 10%.27 The EP/exp-6 model shows a slightly larger error of 2% for the saturated liquid densities but excels for the saturated vapor pressures with an average error less than 1%.29 The OPLSAA model87 performs well at at 298.15 K (where it was fitted), but its performance degrades rapidly with increasing temperature. The VLCC and saturated pressure lines for pyridine, pyrimidine, pyrazine, and pyridazine are depicted in the Figure 6. The saturated liquid densities and vapor pressures of pyridine are well reproduced by the TraPPE-EH model with mean unsigned errors less than 2% and 5%, respectively. The experimental critical temperature, critical density, and the normal boiling point of pyridine fall within the error bars of the

Transferable Potentials for Phase Equilibria

Figure 6. Coexistence densities (left) and saturated vapor pressures (right) for pyridine, pyrimidine, pyrazine, and pyridazine. The solid line, star, dashed line, and plus represent the experimental data for the saturated liquid densities of pyridine,94 the critical point for pyridine,91 the saturated vapor pressure for pyridine,94 and the critical point for pyrazine,93 respectively. The circles, triangles, diamonds, and squares represent the TraPPE-EH predictions for pyridine, pyrimidine, pyrazine, and pyridazine, respectively.

J. Phys. Chem. B, Vol. 111, No. 36, 2007 10797

Figure 8. Coexistence densities (left) and saturated vapor pressures (right) for thiazole, oxazole, isoxazole, imidazole, and pyrazole. The solid line represents experimental data for isoxazole.96 The star represents the experimental critical point for isoxazole.96 The down triangles, squares, circles, diamonds, and up triangles represent the simulation data for thiazole, oxazole, isoxazole, imidazole, and pyrazole, respectively.

Figure 7. Coexistence densities (left) and saturated vapor pressures (right) for thiophene, furan, and pyrrole. The solid, dashed, and dotted lines depict the experimental data for thiophene,95 furan,95 and pyrrole,94 respectively. The star, cross, and plus symbols represent the critical points for the thiophene, furan, and pyrrole, respectively.91 The diamonds, squares, and circles show the simulation data for the TraPPEEH force field for thiophene, furan, and pyrrole, respectively.

Figure 9. Correlation plots comparing the predictions for the TraPPEEH model with experimental data. Left: normal boiling points (circles) and critical temperatures (triangles). Right: critical densities (circles) and saturated liquid densities near the normal boiling point (triangles). Solid symbols show data for compounds that were used in the parametrization process (benzene, pyridine, thiophene, furan, and pyrrole) and open symbols represent data for the validation set. The solid line represents y ) x.

simulation results, but the TraPPE-EH force field overestimates the heat of vaporization at 298.15 K by 2.8 kJ/mol. The transferability of the basic nitrogen LJ parameters was validated through calculation of the VLCCs for the diazenes, where no additional LJ parameters are required. Unfortunatey, experimental phase equilibrium data are rather limited for these compounds. As can be seen from the data in Table 7, the predicted normal boiling points for pyrimidine and pyrazine agree to within error bars with the experimental data whereas the normal boiling point for pyridazine is overestimated by 5%. Critical data are only available for pyrazine,93 a symmetric diazene without dipole moment, and its critical temperature and density are underestimated by 1.8 and 1.2%, respectively. The results for the heats of vaporization are scattered with an underestimation by 4.9 kJ/mol for pyrimidine and an overestimation by 6.8 kJ/mol for pyridazine. Additional LJ parameters for sulfur, oxygen, and acidic nitrogen heteroatoms were obtained by fitting to the VLCC of thiophene, furan, and pyrrole, respectively. As can be seen in Figure 7, the agreement with the experimental data94,95 for the saturated liquid densities and vapor pressures of thiophene, furan, and pyrrole is excellent over the entire coexistence range with mean unsigned errors of less than 1 and 5% for the liquid densites and vapor pressures, respectively. The accuracy of the predicted vapor pressures is also reflected in small errors of 0.5, 0.1, and 0.3% for the normal boiling points of thiophene, furan, and pyrrole, respectively. Although the corresponding errors in the critical temperatures are rather small with deviations

of 1.2, 0.8, and 1.7%, those for the critical densities with deviations 5.5, 8.7, and 3.6%, respectively, are quite large. The sulfur, oxygen, and nitrogen LJ parameters were transferred to thiazole, oxazole, isoxazole, imidazole, and pyrazole. The VLCC and saturated pressure lines for these five compounds are shown together with the experimental data for isoxazole96 in Figure 8. As can be seen from the data in the Table 7, the TraPPE-EH force field performs very well for isoxazole with deviations of less than 1% for the critical temperature, critical density, and normal boiling point and a slight overprediction by 1 kJ/mol for the heat of vaporization. The normal boiling temperatures of thiazole, oxazole, imidazole, and pyrazole are all underpredicted with a mean error of about 3%, but the heat of vaporization of oxazole deviates by only 0.4 kJ/mol. Unfortunately, the critical data for the other compounds are not available. Figure 9 provides a summary of the accuracy of the TraPPEEH force field for vapor-liquid equilibrium properties of the aromatic heterocycles investigated here. As is evident, the overall agreement with experiment is excellent. The regression, R, and correlation, C, coefficients are as follows: normal boiling point, R ) 0.992 and C ) 0.985; critical temperature, R ) 0.977 and C ) 0.990; critical density, R ) 0.692 and C ) 0.848; and saturated liquid density, R ) 1.073 and C ) 0.987. As expected, the performance of the TraPPE-EH model is somewhat better for the five compounds (benzene, pyridine, thiophene, furan, and pyrrole) used in the parametrization process but the differences are not large. The range for the critical densities is

10798 J. Phys. Chem. B, Vol. 111, No. 36, 2007 relatively small, and these are more difficult to obtain experimentally. Thus, it is not surprising that R and C for the critical density are not as satisfactory as those for the other three VLCC properties. 5. Conclusions The TraPPE-EH force field is extended to benzene, pyridine, diazenes, and five-membered aromatic heterocycles with sulfur, oxygen, and nitrogen heteroatoms. As judged against high-level ab initio calculations, the TraPPE-EH model describes the potential surface for the benzene dimer fairly accurately (with a mean unsigned percentage error (MUPE) of 19% for energies and center-of-mass separations), in particular, compared to various united atom benzene force fields. The use of CM4 charges computed in a continuum 1-octanol solvent yields dipole moments that are about 25% larger (with the exception of pyrrole) than experimental gas-phase dipole moments, that is, this procedure allows for a systematic determination of partial charges that account for condensed-phase polarization effects. For benzene and the aromatic heterocyles, the TraPPE-EH force field predicts normal boiling points, critical temperatures, and critical densities with MUPE of 0.8, 1.6, and 3.1%, respectively. The MUPE in the saturated liquid densities and vapor pressures are less than 1 and 5%, respectively. The predicted heats of vaporization near the normal boiling point fall usually within 2 kJ/mol of the experiment data but those for pyridine and pyridazine are overestimated by a larger amount. Overall, the extension of the TraPPE-EH force field with five new Lennard-Jones sites is able to decribes the vaporliquid equilibrium properties of benzene and the aromatic heterocycles very well. Acknowledgment. We thank Casey Kelly for help with the MN-GSM calculations. Financial support from the National Science Foundation (CTS-0553911) is gratefully acknowledged. Part of the computer resources were provided by the Minnesota Supercomputing Institute. Supporting Information Available: Tables listing the numerical data for the TraPPE-EH force field shown in Figures 3-8. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Hunter, C. A. J. Mol. Biol. 1993, 230, 1025-1054. (2) Hunter, C. A. Chem. Soc. ReV. 1994, 23, 101-109. (3) Hunter, C. A.; Sanders, J. K. M. J. Am. Chem. Soc. 1990, 112, 5525-5534. (4) Sponer, J.; Leszczynski, J.; Hobza, P. J. Phys. Chem. A 1997, 101, 9489-9495. (5) Muller-Dethlefs, K.; Hobza, P. Chem. ReV. 2000, 100, 143-168. (6) Hunter, C. A.; Meah, M. N.; Sanders, J. K. M. J. Am. Chem. Soc. 1990, 112, 5773-5780. (7) Breault, G. A.; Hunter, C. A.; Mayers, P. C. J. Am. Chem. Soc. 1998, 120, 3402-3410. (8) Bisson, A. P.; Carver, F. J.; Hunter, C. A.; Waltho, J. P. J. Am. Chem. Soc. 1994, 116, 10292-10293. (9) Braga, D.; Grepioni, F.; Desiraju, G. R. Chem. ReV. 1998, 98, 1375-1406. (10) Thaimattam, R.; Sharma, C. V. K.; Clearfield, A.; Desiraju, G. R. Cryst. Growth Des. 2001, 1, 103-106. (11) Aitipamula, S.; Thallapally, P. K.; Thaimattam, R.; Jaskolski, M.; Desiraju, G. R. Org. Lett. 2002, 4, 921-924. (12) Desiraju, G. R.; Gopalakrishnan, B.; Jetti, R. K. R.; Nagaraju, A.; Raveendra, D.; Sarma, J. A. R. P.; Sobhia, M. E.; Thilagavathi, R. J. Med. Chem. 2002, 45, 4847-4857. (13) Felker, P. M.; Maxton, P. M.; Schaeffer, M. W. Chem. ReV. 1994, 94, 1787-1805. (14) Leopold, K. R.; Fraser, G. T.; Novick, S. E.; Klemperer, W. Chem. ReV. 1994, 94, 1807-1827.

Rai and Siepmann (15) Mueller-Dethlefs, K.; Dopfer, O.; Wright, T. G. Chem. ReV. 1994, 94, 1845-1871. (16) Hobza, P.; Selzle, H. L.; Schlag, E. W. Chem. ReV. 1994, 94, 17671785. (17) Chalasinski, G.; Szczesniak, M. M. Chem. ReV. 1994, 94, 17231765. (18) Sinnokrot, M. O.; Valeev, E. F.; Sherrill, C. D. J. Am. Chem. Soc. 2002, 124, 10887-10893. (19) Sinnokrot, M. O.; Sherrill, C. D. J. Phys. Chem. A 2004, 108, 10200-10207. (20) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2005, 109, 4209-4212. (21) Tauer, T. P.; Sherrill, C. D. J. Phys. Chem. A 2005, 109, 1047510478. (22) Sinnokrot, M. O.; Sherrill, C. D. J. Phys. Chem. A 2006, 110, 10656-10668. (23) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M. J.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (24) Jorgensen, W. L.; McDonald, N. A. J. Mol. Struct. 1998, 424, 145155. (25) Jorgensen, W. L.; McDonald, N. A. J. Phys. Chem. B 1998, 102, 8049-8059. (26) Feller, S. E.; MacKerell, A. D., Jr. J. Phys. Chem. B 2000, 104, 7510-7515. (27) Contreras-Camacho, R. O.; Ungerer, P.; Boutin, A.; Mackie, A. D. J. Phys. Chem. B 2004, 108, 14109-14114. (28) Ahunbay, M. G.; Perez-Pellitero, J.; Contreras-Camacho, R. O.; Teuler, J.; Ungerer, P.; Mackie, A. D.; Lachet, V. J. Phys. Chem. B 2005, 109, 2970-2976. (29) Errington, J. R.; Panagiotopoulos, A. Z. J. Chem. Phys. 1999, 111, 9731-9738. (30) Martin, M. G.; Siepmann, J. I. J. Chem. Phys. 1998, 102, 25692577. (31) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 45084517. (32) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 5370-5379. (33) Wick, C. D.; Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 8008-8016. (34) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 3093-3104. (35) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2004, 108, 17596-17605. (36) Kamath, G.; Cao, F.; Potoff, J. J. J. Phys. Chem. B 2004, 108, 14130-14136. (37) Kamath, G.; Robinson, J.; Potoff, J. J. Fluid Phase Equilib. 2006, 240, 46-55. (38) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 104, 18974-18982. (39) Lubna, N.; Kamath, G.; Potoff, J. J.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 24100-24107. (40) Zhang, L.; Siepmann, J. I. Theor. Chem. Acc. 2006, 115, 391397. (41) Potoff, J. J.; Siepmann, J. I. AIChE J. 2001, 47, 1676-1682. (42) Chen, B.; Xing, J.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 2391-2401. (43) Zhang, L.; Siepmann, J. I. Unpublished work. (44) Rai, N.; Rafferty, J. L.; Maiti, A.; Siepmann, J. I. Fluid Phase Equilib. 2007, doi: 10.1016/j.fluid.2007.06.034. (45) Kelkar, M. S.; Rafferty, J. L.; Siepmann, J. I.; Maginn, E. J. Fluid Phase Equilib. 2007, doi: 10.1016/j.fluid.2007.06.033. (46) Results of the Third Industrial Fluid Properties Simulation Challenge. http://fluid.properties.org/challenge/3rd/results (accessed on May 10, 2007). (47) Sun, L.; Siepmann, J. I.; Klotz, W. L.; Schure, M. R. J. Chromatogr., A 2006, 1126, 373-380. (48) Zhang, L.; Rafferty, J. L.; Siepmann, J. I.; Chen, B.; Schure, M. R. J. Chromatogr., A 2006, 1126, 219-231. (49) Sun, L.; Siepmann, J. I.; Schure, M. R. J. Phys. Chem. B 2006, 110, 10519-10525. (50) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 2006, 110, 3555-3563. (51) Wick, C. D.; Siepmann, J. I.; Sheth, A. R.; Grant, D. J. W.; Karaborni, S. Cryst. Growth Des. 2006, 6, 1318-1323. (52) Kioupis, L. I.; Maginn, E. J. J. Phys. Chem. B 1999, 103, 1078110790. (53) Kioupis, L. I.; Maginn, E. J. J. Phys. Chem. B 200, 104, 77747783. (54) Moore, J. D.; Cui, S. T.; Cochran, H. D.; Cummings, P. T. J. Chem. Phys. 2000, 113, 8833-8840. (55) Lorentz, H. A. Ann. Phys. 1881, 12, 127. (56) Berthelot, D. C. R. Hebd. Seances Acad. Sci. 1898, 126, 1703. (57) Zhao, X. S.; Bhatt, D.; Siepmann, J. I. Unpublished work. (58) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129145.

Transferable Potentials for Phase Equilibria (59) Besler, B. H.; Merz, K. M.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431-439. (60) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361373. (61) Mulliken, R. S. J. Chem. Phys. 1935, 3, 564-573. (62) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833-1840. (63) Mulliken, R. S. J. Chem. Phys. 1962, 36, 3428-3439. (64) Lowdin, P.-O. J. Chem. Phys. 1950, 18, 365-375. (65) Li, J.; Zhu, T.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 1998, 102, 1820-1831. (66) Winget, P.; Thompson, J. D.; Xidos, J. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 2002, 106, 10707-10717. (67) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Theor. Chem. Acc. 2005, 113, 133-151. (68) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Theor. Comput. Chem. 2005, 1, 1133-1152. (69) Bhatt, D.; Zhao, X. S.; Rai, N.; Siepmann, J. I. J. Chem. Phys., submitted for publication. (70) Becke, A. D. J. Chem. Phys. 1993, 98, 5648-5652. (71) Perdew, J. P. Phys. ReV. B 1986, 33, 8822-8824. (72) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.01; Gaussian, Inc.: Wallingford, CT, 2004. (73) Chamberlin, A. C.; Kelly, C. P.; Thompson, J. D.; Xidos, J. D.; Li, J.; Hawkins, G. D.; Winget, P. D.; Zhu, T.; Rinaldi, D.; Liotard, D. A.; Cramer, C. J.; Truhlar, D. G.; Frisch, M. J. MN-GSM, version 6.0; University of Minnesota: Minneapolis, MN, 2006. (74) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813-826. (75) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527-545.

J. Phys. Chem. B, Vol. 111, No. 36, 2007 10799 (76) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (77) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 1996. (78) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59-70. (79) Mooij, G. C. A. M.; Frenkel, D.; Smit, B. J. Phys.: Condens. Matter 1992, 4, L255-L259. (80) Smit, B.; de Smedt, P.; Frenkel, D. Mol. Phys. 1989, 68, 931950. (81) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Oxford University Press: New York, 1989. (82) Rowlinson, J. S.; Swinton, F. L. Liquids and Liquid Mixtures; Butterworth: London, 1982. (83) Atkins, P. W. Physical Chemistry; Oxford University Press: New York, 2002. (84) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796-6806. (85) Gauss, J.; Stanton, J. F. J. Phys. Chem. A 2000, 104, 2865-2868. (86) Wick, C. D.; Siepmann, J. I.; Klotz, W. L.; Schure, M. R. J. Chromatogr. A 2002, 954, 181-190. (87) Jorgensen, W. L.; Severance, D. L. J. Am. Chem. Soc. 1990, 112, 4768-4774. (88) Schroer, J. W.; Monson, P. A. J. Chem. Phys. 2001, 114, 41244130. (89) Zhao, X. S.; Chen, B.; Karaborni, S.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 5368-5374. (90) McClellan, A. L. Tables of Experimental Dipole Moments, Vol. 2; Rahara Enterprise: El Cerrito, 1973. (91) Lemmon, E. W.; McLinden, M. O.; Friend, D. G. Thermophysical Properties of Fluid Systems; NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom P. J., Mallard, W. G., Eds., National Institute of Standards and Technology: Gaithersburg, MD, June 2005 (http://webbook.nist.gov). (92) Chirico, R. D.; Steele, W. V.; Nguyen, A.; Klots, T. D.; Knipmeyer, S. E. J. Chem. Thermodyn. 1996, 28, 797-818. (93) Marsh, K.; Young, C.; Morton, D.; Ambrose, D.; Tsonopoulos, C. J. Chem. Eng. Data 2006, 51, 305-314. (94) Das, A.; Frenkel, M.; Gadalia, N. A. M.; Kudhchadker, S.; Marsh, K. N.; Rodgers, A. S.; Wilhoit, R. C. J. Phys. Chem. Ref. Data 1993, 22, 659-782. (95) In Perry’s Chemical Engineers’ Handbook, 7th ed.; Perry, R. H., Green D. W., Eds.; McGraw-Hill: New York,1997; Chapter 2. (96) Steele, W. V.; Chirico, R. D.; Knipmeyer, S. E.; Nguyen, A.; Smith, N. K.; Tasker, I. R. J. Chem. Eng. Data 1996, 41, 1269-1284.