First-Principles Many-Body Force Fields from the Gas Phase to Liquid

Article pubs.acs.org/JPCB

First-Principles Many-Body Force Fields from the Gas Phase to Liquid: A “Universal” Approach Jesse G. McDaniel and J. R. Schmidt* Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin-Madison, Madison, Wisconsin 53706, United States S Supporting Information *

ABSTRACT: We extend our previously developed approach for generating “physicallymotivated” force fields from symmetry-adapted perturbation theory by introducing explicit terms to account for nonadditive three-body exchange and dispersion interactions, yielding transferability from the gas- to condensed-phase. These Axilrod− Teller−Muto-type three-body terms require no additional parametrization and can be implemented with high computational efficiency. We demonstrate the accuracy of our force fields for a diverse set of six organic liquids/fluids, examining a wide variety of structural, thermodynamic, and dynamic properties. We find that three-body dispersion and exchange interactions make significant contributions to the internal pressure of condensed phase systems and cannot be neglected in truly ab initio force field development. These resulting force field parameters are extremely transferable over wide ranges in temperature and pressures and across chemical systems, and should be widely applicable in condensed phase simulation.

■

INTRODUCTION Molecular simulations are widely used to predict and explain the structural, thermodynamic, and dynamic properties of condensed phase chemical systems. As such, there has been considerable development of corresponding classical force fields (e.g., AMBER,1 CHARMM,2 GROMOS,3 OPLS,4 TraPPE5). These force fields are typically empirical in that they are extensively parametrized to reproduce experimental data, often providing atom-specific parameters for molecular species encompassing a wide range of chemical systems. Although such force fields have been extremely successful and widely adopted, empirically parametrized force fields have two main drawbacks that hinder their accuracy, generality, and transferability: they rely on the availability of experimental data for parametrization, which may be nonexistent for complex or novel systems; and force field parameters are usually highly coupled, and well-determined fits are possible for only a modest number of parameters, thus requiring relatively simple (but suboptimal) functional forms (e.g., Lennard-Jones). An alternative approach is to parametrize force fields from ab initio electronic structure calculations, which introduces both new opportunities and new challenges. In particular, nuclear quantum effects (at low temperatures) and nonadditive manybody interactions must be treated explicitly, rather than implicitly including them within empirically parametrized “effective” pairpotentials. In the latter case, it is generally recognized that nonadditive many-body interactions make nontrivial contributions to the energy and other properties for a wide variety of chemical systems. In the context of intermolecular perturbation theory, these many-body interactions can be separated into contributions from induction, dispersion, and exchange. Many© 2014 American Chemical Society

body induction, the only term with a corresponding classical analogue, is the most widely treated of these terms, and many force fields now employ fluctuating charge models, induced point multipoles, or Drude oscillators to account for this effect. In contrast, many-body dispersion and exchange (though well studied in many different systems6−21) have yet to be commonly incorporated into force fields in a general way. In this work, we extend our previous methodology22 for generating physically motivated, transferable force fields from symmetry-adapted perturbation theory (SAPT) to include many-body dispersion and exchange interactions. The resulting approach yields accurate, first-principle force fields for a wide variety of chemical systems, transferable from the gas-phase to condensed phase, and incorporating by construction almost all of the relevant “physics” of intermolecular interactions. Previously we obtained transferable parameters describing two-body interactions and many-body induction (in the latter case via polarizable Drude oscillators). While we were able to accurately reproduce experimental second virial coefficients and gas adsorption properties for a wide range of chemical systems, the magnitude of the error introduced by neglecting many-body dispersion and many-body exchange interactions in condensed phase systems (e.g., bulk liquids) is not a priori clear. In this work, we examine six different neat organic fluids: methane, ethane, acetone, acetaldehyde, methanol, and methyl amine. We show that inclusion of Axilrod−Teller−Muto (ATM) like three-body Special Issue: James L. Skinner Festschrift Received: January 31, 2014 Revised: March 18, 2014 Published: March 21, 2014 8042

dx.doi.org/10.1021/jp501128w | J. Phys. Chem. B 2014, 118, 8042−8053

The Journal of Physical Chemistry B

Article

force field terms to effectively model the (partially canceling) nonadditive three-body exchange and dispersion interactions yields force fields that are transferable from the gas phase to the bulk liquid, and beyond. Crucially, these ATM force field terms require no additional calculations, parametrization, or methodological changes.22

accurate two-body potential is sufficient to provide an excellent description of the structural and thermodynamic properties of condensed-phase CO2.20 Because our force fields employ Drudeoscillator polarizability, the E3b ind term is (in principle) already included and we find that the three-body portion of the Drudeoscillator energy does a good job of reproducing the three-body 20 SAPT-computed E3b Therefore, we focus on the ind energy. 3b remaining three-body exchange and dispersion, E3b exch and Edisp. In the asymptotic limit of no monomer electron overlap, threebody interactions can be written in terms of monomer properties (just as in the case of two-body interactions27). Therefore, we first examine the three-body dispersion interaction in regions 3b with no electron density overlap, Edisp−lr . The complete multipolar expansion for the three-body dispersion energy between three spherically symmetric atoms, exact in the asymptotic limit, was derived by Bell.28 The leading terms of this expansion (triple-dipole, and dipole−dipole−quadrupole) are

■

THEORY AND METHODOLOGY It is well-known that many-body nonadditive interactions make important contributions to the internal energy and pressure of condensed phase systems, and therefore must be either explicitly or implicitly accounted for in an accurate force field. The importance of different types of these nonadditive many-body interactions is system dependent; in polar systems many-body induction interactions are expected to dominate, whereas in highly nonpolar systems (e.g., noble gases) the many-body dispersion and exchange interactions dominate. For the case of noble gases,7 many-body contributions have been successfully treated with the Axilrod−Teller−Muto (ATM) potential,23,24 which is the leading order term (triple-dipole) in the multipoleexpanded three-body dispersion energy. Many-body nonadditive interactions can perhaps be best understood through the energy decomposition given by threebody SAPT10,11,25 and its computationally efficient counterpart, three-body DFT-SAPT26 (here three-body refers to three atoms or molecules). We will restrict our attention to three-body DFTSAPT as this is the method that we use for our computations. Following the notation in ref 25, the total three-body energy can be written as:

3b Edisp − lr ≈ C 9

1 + 3 cos ϕA cos ϕB cos ϕC 3 3 3 RAB RACRBC

+ PABC

⎡ (9 cos ϕC − 25 cos 3ϕC) + 6 cos(ϕA − ϕB)(3 + 5 cos 2ϕC) ⎤ ⎥ × ⎢C11ABC 3 4 4 RAB RBC RAC ⎣ ⎦ + ...

where ϕA, ϕB, and ϕC are the angles formed by the triangle connecting the three atoms, and PABC is an operator that generates all cyclic permutations of (ABC). The first term in this expansion is the well-known Axilrod−Teller−Muto triple-dipole term.6,23,24 Similar to the two-body dispersion coefficients, the three-body dispersion coefficients are given in terms of imaginary frequencydependent polarizability tensors of the atoms:10

ABC E int = E int[2, 3] + E int[3, 3]

where Eint[2,3] is the two-body contribution to the three-body energy and Eint[3,3] is the nonadditive part of the three-body energy (henceforth we will drop the “non-additive” nomenclature, so that any reference to a “three-body energy” refers to the nonadditive three-body energy). The Eint[2,3] term is already given by our previously developed force fields, and thus our focus is on constructing a methodology for describing the three-body term, Eint[3,3], in a general, transferrable manner. The three-body interactions that we consider (as in ref 20) are

C9 =

3 π

C11ABC =

∫0

∞

3 16π

αA1 (iω)αB1(iω)αC1(iω)dω

∫0

∞

αA1(iω)αB1(iω)αC2(iω) dω

where αiA(iω) is the ith rank polarizability (tensor has one unique nonzero element for spherical atoms) at the imaginary frequency iω for atom A. We have already described how to obtain transferrable, atomic parameters for these frequency-dependent polarizabilities that reflect specific molecular environments,22 allowing for the description of the multipole-expanded threebody dispersion energy with no additional parametrization. However, such a description becomes cumbersome to implement within a molecular simulation and evaluating higher-order terms in this expansion is expensive. Therefore, we also investigate fitting a scale factor to produce an effective coefficient C̃ 9, such that C̃ 9 = λ3bC9, and

SAPT(DFT) 3b 3b 3b E int [3, 3] = Eexch + Edisp + E ind , (2)

3b (1) ̃ − disp[3, 3](CKS), Eexch = Eexch [3, 3](KS) + Eexch 3b (3) Edisp = Edisp [3, 3](CKS), (2)

3b (2) ̃ − ind [3, 3](CKS) E ind = E ind [3, 3](CKS) + Eexch

where the definitions of each of these interaction terms is given in 3b ref 26. Briefly, E3b exch is the three-body exchange, Edisp is the threebody component to dispersion, and E3b is the three-body ind induction energy. The index in parentheses (i.e., “2” in (2) Eind [3,3](CKS)) denotes the perturbation order in the intermolecular interaction operator and “CKS” indicates that the response was computed using the coupled Kohn−Sham response. In addition to these terms, Bukowski and Szalewicz have also considered a fourth-order dispersion term, E(4) disp[3,3] in their study of argon,7 but the current implementation of the three-body DFT-SAPT code cannot calculate this term at the coupled KS level,26 and therefore we neglect this higher order term. We have found in our previous work that an accurate treatment of these three-body terms in conjunction with an

3b ̃ Edisp − lr ≈ C 9


where the effective C̃ 9 coefficient best reproduces the three-body dispersion energy using the truncated multipole expansion. This is similar in spirit to fitting an effective C6 coefficient for two-body dispersion interactions to account for the truncation of the twobody multipole-expanded dispersion energy.29 To reproduce the three-body dispersion energy for all configurations (including those with nonzero electron density overlap), we introduce a damping function:10,18 8043



Article

with β = 2.01, as in ref 20. Fits to the three-body dispersion and three-body exchange energies are given in the Supporting Information. The ATM C9 coefficients were generated using the frequency-dependent polarizability tensors reported in our previous work,22 and the exponents βAB used in the exchange energy fit and in the f 9 damping function were also from our previous work. A scale factor of λ3b ≈ 1.4 was required to accurately reproduce the three-body dispersion energy, clearly indicating that multipole terms of higher order than the ATM term in the three-body dispersion expansion are important and an ATM-like term can only accurately reproduce the dispersion energy using a larger effective C̃ 9 coefficient. Inspecting the three-body exchange and three-body dispersion fits, it is clear that these terms partially cancel for a majority of configurations since E3b exch is generally 3b attractive (unlike 2-body exchange) and Edisp is generally repulsive (unlike 2-body dispersion). In general, the magnitude 3b of E3b disp is greater than that of Eexch for most configurations, so that the net sum of these terms is usually repulsive. To investigate whether we can exploit this partially cancellation of terms in our force field development, we note that a truncation of the multipole expanded dispersion energy at the ATM term results in under-prediction of the repulsive E3b disp term by ∼30% (note the fitted 1.4 scale factor). Therefore, the 3b bare ATM term may reproduce the sum E3b disp + Eexch more 3b accurately than Edisp due to this partial cancellation.

3b Edisp ≈ f9 (βAB , βAC , βBC , RAB , RAC , RBC)C̃9

×


where f9 (βAB , βAC , βBC , RAB , RAC , RBC) = f3 (βAB , RAB)f3 (βAC , RAC)f3 (βBC , RBC)

is a product of three two-body Tang-Toennies damping functions, n

fn (λ , r ) = 1 − e−λr

∑ m=0

(λr )m m!

The three-body exchange energy E3b exch between three atoms is fit to 3b exch Eexch ≈ AABC exp( −βABC(RAB + RAC + RBC))

where for identical atoms A, B, C, βABC = βAB/2,20 and βAB are the exponents for the two-body exchange terms.22 Therefore in order to entirely account for the three-body interaction energies described above, two new types of parameters must be fit from ab initio data: λ3b for the threebody dispersion energies, and Aexch ABC parameters for the threebody exchange interactions. As a test system, we choose to focus on methane, for its simplicity, computational convenience, and the fact that we have previously developed an accurate two-body force field.22 Ab-Initio Investigation of Three-Body Interactions in Methane. Three-body SAPT calculations were conducted for ∼100 methane trimer configurations. These configurations were generated by sampling from a molecular dynamics (MD) simulation at a dense supercritical phase (200 K, 200 bar) employing the OPLS-AA4 force field to describe methane interactions. The trimers were generated by first randomly selecting a methane molecule, and then finding all neighboring molecules (closer than 5 Å) and picking out all unique trimers in which all three molecules were neighbors. The three-body energy component E(3) disp[3,3](CKS) was computed with aug-cc-pVDZ (AVDZ) and aug-cc-pVTZ (AVTZ) basis sets,30 and the energy (1) (2) components E exch [3,3](KS), E exch−disp [3,3](KS), and E(3) [3,3](KS) were computed with an AVDZ basis set. All disp calculations utilized the PBE functional,31,32 and were conducted using the SAPT2008 suite of codes. 2 6 The term Ẽ (2) exch−disp[3,3](CKS) was then approximated by:

Figure 1. Comparison of the Axilrod−Teller−Muto energy E3b ATM to the 3b three-body SAPT E3b disp + Eexch energies for different configurations of methane trimers.

(2)

̃ − disp[3, 3](CKS) Eexch (2) ≈ Eexch − disp[3, 3](KS)AVDZ

In Figure 1, we plot the sum of three-body SAPT terms, E3b disp + compared to the predicted energy of our damped, unscaled, ATM energy expression for methane.

(3) Edisp [3, 3](CKS)AVDZ

E3b exch,

(3) Edisp [3, 3](KS)AVDZ

3b EATM =

A complete basis set extrapolation was used to approximate E(3) disp[3,3](CKS)CBS, so that (3) Edisp [3,

≈

∑

f3 (βAB , RAB)f3 (βAC , RAC)f3 (βBC , RBC)

A,B,C = C,H

× C9ABC

3](CKS)CBS

3β (3) Edisp [3, 3](CKS)AVTZ 3β − 2β 2β (3) [3, 3](CKS)AVDZ − β Edisp 3 − 2β


where the sum is over all carbon and hydrogen atoms in each methane molecule, such that atoms A, B, C are each on separate molecules. The parameters βCC,βCH,βHH,CCCC ,CCCH ,CCHH ,CHHH 9 9 9 9 n are all either given or derived from parameters (αA(iω)) given in our previous work,22 and are listed in the Supporting 8044



Article

Information. Therefore, no new parameters were needed to generate Figure 1. The ATM term semiquantitatively reproduces the total three3b 3b body dispersion and exchange, E3b ATM ≈ Edisp + Eexch, for the vast majority of configurations − with no additional required parameters. The partial cancellation of three-body energy components has been previously observed by Bukowski and Szalewicz, in the context of noble gases.7 Those authors found that the extent of this cancellation is somewhat configuration dependent, but seems robust in typical condensed phase configurations. As such, we will examine both explicit fits to 3b 3b E3b dispand Eexch (separately), as well as EATM. Methane Third Virial Coefficient. As a test of our ability to successfully incorporate three-body effects in our force field, we compute the third virial coefficient, B3, of methane. Our previously developed force field for methane accurately reproduces the second virial coefficient,22 allowing a clean separation of two- and three-body effects. The third virial coefficient is calculated using the formula33

Figure 2. Third virial coefficient of methane as a function of temperature. The experimental data are shown as black circles and squares, while the calculations are shown as red, blue, and purple triangles. The model “2-body FF” employs no three-body terms; 3b “+E3b ATM” additionally includes ATM-type three body terms; and “+ Edisp 3b 3b + E3b exch” includes (alternatively) both the explicitly fit Edisp and Eexch three-body terms.

B3 = B3pp + B33b where Bpp 3 is the pair-potential contribution to the virial, 1 B3pp = − ⟨ 3

∫ ∫ f12 f13 f23 dr12dr13⟩Ω Ω Ω 1

2

focus to those organic systems for which we have previously derived accurate and transferable polarizable force fields,22 allowing us to isolate the role of three-body effects. In addition, we will initially limit our attention to small-molecule systems for which (we believe) a rigid-body treatment is a good approximation and thus can avoid the added uncertainties associated with intramolecular force field components. (The influence of molecular flexibility is addressed briefly in Discussion.) Therefore, we investigate six neat-organic systems: methane, ethane, acetone, acetaldehyde, methyl amine, and methanol. These systems span a range of interaction types, from nonpolar, dispersion dominated (methane, ethane), through highly polar and hydrogen bonding (methanol). In each case, we focus on bulk structural, thermodynamic, and dynamic properties for which reliable experimental reference data exists. These results should thus allow for a general analysis of the overall role of three-body interactions in condensed phases. The existing two-body force field parameters were supplemented with ATM terms (via frequency-dependent polarizabilities); in both cases parameters were taken from our previous work.22 The geometries of the organic molecules were taken from gas-phase minimum energy configurations. There were two small modifications to the force field parameters of ethane and methanol, as described in detail in the Supporting Information. Briefly, the Drude oscillators on methanol were refit to be more consistent with the other force field parameters, resulting in trivial changes in the force field and second virial coefficient. In addition, the exchange parameters on ethane were scaled slightly to better reproduce the second virial coefficient, facilitating a clean separation of errors introduced by two- and three-body terms. The modified force field parameters for ethane and methanol, as well as a comprehensive list of C9 ATM coefficients, are given in the Supporting Information. Simulation Methodology. The densities of neat fluids were calculated using Monte Carlo simulations in the NPT ensemble, while enthalpies of vaporization and structure factors were calculated in the NVT ensemble. Computations of these latter properties in the NPT ensemble lead to similar results but with errors generally slightly compounded due to small density differences; see Supporting Information. All Monte Carlo

3

with the Mayer function given by

f12 = e−E12 / kT − 1 and B3b 3 is the nonadditive three-body energy contribution to the virial, 1 B33b = − ⟨ 3

∫ ∫ (e−ΔE

123 / kT

− 1)e−(E12 + E13+ E23)/ kT dr12dr13⟩Ω1Ω2Ω3

In all cases, the average ⟨...⟩Ω1Ω2Ω3 is over molecular orientations of the three molecules. These integrals can be rewritten as (taking Bpp 3 as an example) B3pp = −

8π 2 3

∫ ∫ ∫ ⟨f12 f13 f23 ⟩Ω Ω Ω r122r132 1

2

3

× sin θ dθ dr12 dr13

where θ is the angle between the two vectors r12 and r13. The two distance integrals were numerically integrated, while Monte Carlo integration was used for all angular integrals, with 106 configurations per 2-D grid point. The calculated third virial coefficients are compared to the experimental data of Douslin et al.34 (black circles) and Pope et al.35 (black squares) in Figure 2. Both approaches, either using 3b the explicit E3b exch and Edisp fitted potentials, or using the effective 3b ATM term, EATM, do an excellent job of reproducing the experimental third virial coefficients compared to the strictly two-body methane potential. Note that, for methane, our force field does not use Drude polarization, as induction effects are very small (however, 2-body induction due to charge penetration is treated22) and therefore the three-body effects included are entirely due to exchange and dispersion. The success of the simple effective ATM term motivates us to investigate the generality of this approach.

■

RESULTS We now explore the general incorporation of three-body effects for a range of condensed-phase systems. We initially limit our 8045



Article

simulations were conducted using an in-house code.36 For each system, the number of molecules was chosen sufficiently large such that for the densest f luid a cubic box with edge length of 25− 30 Å could be used (40−45 Å, for structure factor calculations), and then this same number of molecules was used for all other lower density phase points. The particle mesh Ewald (PME) algorithm37 was used to treat electrostatics, with a Gaussian width parameter of 0.3 Å−1, 60 reciprocal space grid points, sixthorder Beta-spline interpolation, and a 10 Å real-space cutoff (0.1 Å−1 width, 60 grid points, and 25 Å cutoff, for dilute gases with large simulation boxes). All translation/rotation Monte Carlo moves employed a hybrid MD/MC algorithm,38 with trial configurations generated by numerically integrating Newton’s equations. Center of mass positions and velocities were propagated using the velocity Verlet algorithm, and Euler’s equations for rigid-body motion were propagated using a similar algorithm.39 The integration time step was dynamically updated to achieve a ∼40% acceptance rate, typically ∼10 fs. The force and energy at each configuration were calculated after optimizing the positions of all Drude oscillators using a conjugate-gradient algorithm.40 Intramolecular Drude oscillator interactions employed Thole-type screening as described previously.41 van der Waals interactions were treated using a cutoff of 12 Å (16 Å for structure factor calculations), and the standard long-range correction was added to the energy. Three-body interactions were implemented using a cell link-list,42 and were only considered if all three atomic distance separations were less than 7 Å; see Supporting Information. (Note that potentially more efficient algorithms for calculating three-body interactions are possible.43) To compute diffusion coefficients, molecular dynamics (MD) simulations were conducted using the GROMACS44,45 software package, using tabulated potentials and similar simulation parameters as those in our Monte Carlo calculations; see Supporting Information for implementation details. ATM threebody interactions as well as long-range corrections to the (cutoff) van der Waals interactions were not included in these MD simulations for practical reasons; the neglect of these terms (at constant density) should have minimal influence on the calculated diffusion coefficients. Methane. The third virial coefficient of methane has already been discussed. To test the ability of our force field(s) to model condensed phases of methane, we compute the density of methane at various points in the supercritical region of the T−P phase diagram. In Figure 3, we plot percent errors with respect to the experimental density, for calculations with and without threebody interactions. It is evident that at lower densities (P < 50 bar), all models perform well, as expected, since the second virial coefficients of all models are accurate. For higher-densities, the importance of three-body interactions becomes apparent. While the 2-body FF model exhibits errors up to 8% for the points studied, the “+E3b ATM” model reduces the maximum error (at 300 K, 200 bar) to ∼4%, while most errors are well within 2%. Inclusion of explicit and separate exchange and dispersion terms 3b 3b in the more sophisticated “+Edisp + Eexch ” model yields exceptional performance, with errors for all points within 1% or less. Note that neither model contains any adjustable parameters. The enthalpy of vaporization of methane at the experimental boiling point (111.7 K) is compared to experiment in Table 1 (the boiling point for the model has not yet been determined). Note that under these conditions there is somewhat larger error (∼6−8%) for both the predicted enthalpies of vaporization and

Figure 3. Calculated density of methane at different temperature and pressures in the supercritical regime. The definition of the models is identical to Figure 2. Percent errors relative to experiment are plotted as circles, where increasing circle size corresponds to increasing error. Solid and dashed circles correspond to positive and negative density errors, respectively. Points that appear to lie on the temperature axis correspond to pressures of 1 bar. Experimental data was taken from Lemmon et al.46

Table 1. Enthalpy of Vaporization at 111.7 K, 0.421 g/cm3, and Density at 111.7 K, 1 atm for Methanea 2-body FF NVT (0.421 g/cm3)

2.18

NPT (1 atm)

0.481

3b +E3b disp + Eexch

+ E3b ATM

ΔHvap (kcal/mol) 2.07 2.09 ρ (g/cm3) 0.453 0.453

experiment 1.95 0.421

a

The definition of the models is identical to Figure 2. Experimental data is from ref 49.

the predicted densities. We believe that this can at least be partially attributed to the neglect of nuclear quantum effects, which become important at low temperature. Calculations utilizing the Feynman−Hibbs effective potential47,48 suggest that nuclear quantum effects lead to somewhat higher internal energies (less negative) and lower densities of liquid methane at this temperature, improving the agreement with experiment; see Supporting Information for details. We next investigate the ability of our models to reproduce the liquid structure. The experimental observable is the structure factor (through the scattering function), facilitating a more direct comparison than radial distribution functions. The structure factor of methane is related to the half-Fourier sine transform of the C−C radial distribution function:50 Hd(k) ∝ 4π

∫0

∞

r 2(gCC(r ) − 1)

sin(kr ) dr kr

where Hd(k)is the distinct structure factor defined in ref 50. To account for the fact that the proportionality factor is indeterminate in the above expression, the calculated structure factor was normalized such that the amplitude of the first peak is identical to that reported by the experiment. The calculated structure factor is compared to experiment in Figure 4, and all experimentally observable peaks are extremely well reproduced. Note that all models of methane (2-body FF, 3b 3b +E3b disp + Eexch, +EATM) give essentially identical structure factors at constant density, regardless of the inclusion of three-body effects. This is in agreement with previous findings for liquid argon.7 8046



Article

The calculated third virial coefficient of ethane is compared to the experimental results of Pope et al.35 in Figure 6. The

Figure 4. Computed structure factor of methane (red) compared to experiment (black) at 92 K, ρ = 0.452 g/cm3. The experimental results were taken from Habenschuss et al.50 Simulation results for small K points are not shown due to large finite system-size effects. The +E3b ATM model was used to compute the structure factor shown, but the other models give essentially identical results.

Figure 6. Third virial coefficient of ethane as a function of temperature. The experimental data of Pope et al.35 are shown in black squares, while model results are shown as blue and purple triangles.

introduction of the three-body “+E3b ATM” terms greatly improves the agreement with experiment, and therefore these terms appear to be successfully modeling the important three-body interactions in ethane. Note that the prediction of the rapid low temperature drop in the third virial coefficient is extremely sensitive to the two-body potential; see Supporting Information. With the accurate treatment of two-body (second virial) and three-body interactions (third virial), we expect that our force field will also accurately describe the properties of condensedphase ethane. The density of ethane at a variety of temperatures and pressures spanning liquid, vapor, and supercritical regions of the phase diagram was computed and is compared to experiment in Figure 7. In the vapor region of the phase diagram, three-body effects are negligible, and the models predict very accurate densities with or without these terms. Similar to our results for methane, for dense liquid and supercritical ethane, the omission of three-body effects can regularly lead to the overprediction of density by ∼5−10%, with larger errors observed closer to the critical point. Introduction of “+E3b ATM” terms to account for

The diffusion coefficients of methane over a range of temperatures (140 K, 160 K, 223 K, 298 K) and pressures (∼5−1600 bar) are compared to experimental values in Figure 5.

Figure 5. Diffusion coefficients for methane at temperatures of 140 K (black), 160 K (red), 223 K (green), 298 K (blue). The inset shows the data at 140 and 160 K on a linear plot for greater clarity. Experimental data is shown as open circles with solid lines connecting the data points, and is taken from refs 51 and 52, while model results are shown as solid squares with error bars. Note that no three-body ATM terms were used when computing these diffusion coefficients; see Simulation Methodology section.

The diffusion coefficients were computed at the corresponding experimental fluid density. In all cases, there is excellent agreement between our force field predictions and experiment. Ethane. As mentioned above, our previous SAPT-based twobody force field for ethane was modified by introducing a small scale factor to improve the agreement with the experimental second-virial coefficient, and therefore better isolate the influence of three-body effects. All results presented in this section utilize this optimized force field; however, in the Supporting Information, we present identical results for all computed quantities using the original (unscaled) force field.

Figure 7. Calculated density of ethane at different temperature and pressures, as in Figure 3. Solid and dashed circles correspond to positive and negative density errors, respectively. Experimental data was taken from Lemmon et al.46 8047



Article

three-body effects greatly improves the agreement with experiment, reducing the errors to below 2−3% for all phase points studied. The enthalpy of vaporization of ethane at the boiling point (184.4 K) was computed and is given in Table 2. The Table 2. Enthalpy of Vaporization for Ethane at 184.4 K, 0.546 g/cm3 and Density at 184.4 K, 1 atma 2-body FF

a

NVT (0.546 g/cm3)

3.78

NPT (1 atm)

0.590

+E3b ATM ΔHvap (kcal/mol) 3.63 ρ (g/cm3) 0.556

Experiment 3.52 0.546

Experimental data is from ref 53.

introduction of three-body exchange/dispersion terms significantly reduces both the error in the computed enthalpy of vaporization (∼3% for +E3b ATM model, ∼7% for 2-body FF model) and density errors (∼2% for +E3b ATM model, ∼8% for 2-body FF model). Note that the difference in the predicted enthalpy of vaporization at NVT between the two simulations is equal to the contribution of E3b ATM to the internal energy of the liquid. The intermolecular distinct structure factor, Hd(k), of ethane (defined in ref 54) at T = 181 K, ρ=0.553 g/cm3, was computed from the intermolecular C−C radial distribution function and is compared to the experimental results of Sandler et al.54 in Figure 8. There is good agreement between our predicted and the

Figure 9. Diffusion coefficients of ethane at 250 bar and a variety of temperatures. The experimental data are the black open circles fit by a black curve, and are from ref 55. The red squares are the computed diffusion coefficients at the corresponding experimental densities. Note that no three-body ATM terms were used when computing these diffusion coefficients; see Simulation Methodology section.

body polarization effects.22 As before, we are interested in isolating the importance of the three-body dispersion and exchange contributions. Therefore two types of simulations were run, “2-body FF + Drude”, which utilize all force field 3b components except E3b ATM terms, and those called “+EATM”, which utilize all force field components including the E3b ATM energy contribution; in both cases, many-body induction is explicitly included. The density of acetone in the liquid phase at various temperatures and pressures is compared to experiment in Figure 10. The densities calculated using the (polarizable) two-body

Figure 8. Computed structure factor of ethane (red) compared to experiment (black) at 181 K, ρ = 0.553 g/cm3. The experimental results were taken from Sandler et al.54 Simulation results for small K points are not shown due large finite system-size effects. Figure 10. Calculated density of acetone at different temperature and pressures, as in Figure 3. Solid and dashed circles correspond to positive and negative density errors, respectively. Experimental data was taken from ref 56.

experimental structure factor. The structure factor that is shown was computed employing the E3b ATM energy contribution, but we find that, as before, the structure factor (in the NVT ensemble) is insensitive to the inclusion of three-body effects. Diffusion coefficients for liquid and supercritical ethane spanning ∼135−450 K at 250 bar were computed and are compared to experiment in Figure 9. All diffusion coefficients were computed at the corresponding experimental density. In general, agreement between the model predictions and experiment is excellent. Acetone. The remainder of the systems that we study have been modeled using Drude oscillators to explicitly treat many-

model systematically predict too dense liquids by 2−3%, while introduction of three-body exchange and dispersion terms systematically predicts lower densities compared to experiment by 1−2%. Therefore, the E3b ATM contribution results in a net change of ∼3−4% in the liquid densities at the temperatures and pressures studied. While there is only small net improvement in the predicted density upon inclusion of the three-body E3b ATM terms, this is likely due to the fact that the two-body potential is 8048



Article

slightly too repulsive,22 and therefore the predicted densities without the E3b ATM contribution are lower than they would be with a “perfect” two-body potential. The enthalpy of vaporization of acetone at 298 K as computed by our force field is compared to experiment in Table 3. For the

Table 5. Temperature and Pressure Dependent Density of Acetaldehydea

Table 3. Enthalpy of Vaporization of Acetone at 298 K, 0.784 g/cm3 and Density at 298 K, 1 atma 2-body FF + Drude

a

NVT (0.784 g/cm3)

7.58

NPT (1 atm)

0.806

+E3b ATM

ΔHvap (kcal/mol) 7.40 ρ (g/cm3) 0.771

temperature (K)

pressure (bar)

2-body FF + Drude

+E3b ATM

experiment

295 295

100 1000

0.869 0.921

0.833 0.884

0.786 0.855

a

All densities are reported in units of of g/cm3. Experimental results have been taken from ref 56.

experiment

interactions, etc. begin to break down. We have partially accounted for this by employing a piecewise pair-potential for the hydrogen bonding atoms (see ref 22). While empirical force fields potentially account for this hydrogen bonding in an average sense, our force field allows for a greater understanding of the necessary force field components needed to accurately model hydrogen bonding, and this is important for future force field development and refinement. In Figure 11, we compare computed and experimental densities of methanol across different points in the temper-

7.48 0.784

Experimental data is taken from ref 53.

two-body model, the enthalpy of vaporization is slightly too high (slightly too “attractive”), leading to a liquid density that is overestimated by ∼2−3%. Introduction of three-body exchange and dispersion decreases the enthalpy of vaporization and predicts an internal energy of the liquid that is slightly “too repulsive”, yielding a liquid density that is too low by ∼1−2%. These results are intuitively consistent and are entirely in agreement with the density predictions in Figure 10. Overall, the agreement of the “+ E3b ATM” model with experiment is excellent. The diffusion coefficient of acetone at 298 K, 1 atm was computed at the corresponding experimental density (0.784 g/ cm3). Our predicted value of 3.66 ± 0.20 x10−9 m2/s is in good agreement with the experimental value of 4.32 ± 0.26 x10−9 m2/s as measured by Ahn et al.57 Acetaldehyde. The liquid density and the enthalpy of vaporization at the boiling point (293.6 K) of acetaldehyde as computed by our force field are compared to experiment in Table 4. Inclusion of E3b ATM three-body terms improves the computed Table 4. Enthalpy of Vaporization of Acetaldehyde at 293.6 K, 0.778 g/cm3, and Density at 293.6 K, 1 atma 2-body FF + Drude

a

NVT (0.778 g/cm3)

6.78

NPT (1 atm)

0.867

+E3b ATM


Figure 11. Calculated density of methanol at different temperatures and pressures, as in Figure 3. Solid and dashed circles correspond to positive and negative density errors, respectively. The errors for the three gas phase points at 1 bar are essentially zero (and have been artificially increased on the plot for visibility). Experimental data was taken from Lemmon et al.46



ature−pressure phase diagram. The errors compared to experiment for dense liquid and supercritical fluids are greater than those of the previous systems that we have studied, with errors of the (polarizable) two-body model between ∼9 and 15%. Introduction of three-body terms significantly reduces these errors (between 4 and 9%), but does not entirely eliminate them. It is evident that the inclusion of E3b ATM three-body force field terms significantly enhances the agreement with experiment, thus further extending the general importance of this three-body contribution to polar, hydrogen-bonding dominated systems. The discrepancy between our force field predictions and experiment is likely due to the two-body treatment of dimer hydrogen bonding interactions; see Supporting Information. The enthalpy of vaporization of methanol is computed at 298 K and compared to experiment in Table 6. While the error of the enthalpy of vaporization predicted at constant density is on the order of the error of the predicted density at constant pressure (∼7−8% for the + E3b ATM model), we see the somewhat surprising result that the enthalpy of vaporization is too low (internal

enthalpy of vaporization by ∼3%, and improves the computed density by ∼5%. Omitting three-body terms yields a density too high by ∼11% as compared to experiment, while including E3b ATM terms reduces the error to 5%. The residual discrepancy is likely due to the fact that the two-body force field is slightly too attractive compared to experiment (see second virial coefficient in ref 22). For a more thorough analysis, we compare predicted densities to experiment at higher pressures (Table 5). In general, there is good agreement between the experimental and predicted properties of acetaldehyde with the “+E3b ATM” model. Methanol. Methanol presents a unique and challenging test because hydrogen bonding must be accurately described in order to successfully predict condensed phase properties. Hydrogen bonding is challenging to model due to the fact that there is a very short contact distance (∼1.7−2.0 Å) between the oxygen and hydrogen bonding atoms near the energy minimum. At such short distances, many approximations commonly used in force fields (including ours) such as multipole expansions, isotropic 8049



Article

Table 6. Enthalpy of Vaporization of Methanol at 298 K, 0.786 g/cm3 and Density at 298 K, 1 atma 2-body FF + Drude

a

NVT (0.786 g/cm3)

8.58

NPT (1 atm)

0.904

+E3b ATM


Table 7. Enthalpy of Vaporization of Methyl Amine at 266.6 K, 0.694 g/cm3, and Density at 266.6 K, 1 atma

experiment

2-body FF + Drude

8.95 0.786 a


energy is “too repulsive”), while the predicted density is too high. This is in contrast to the results from all previous molecular systems, in which there is a correlation between the enthalpy of vaporization and density both being predicted as either too high or too low. We again attribute this seemingly incongruous result to errors in the modeling of the two-body hydrogen bonding dimer interactions and the orientational dependence of hydrogen bonding (see Supporting Information). Yamaguchi et al.58 have extracted partial structure factors of methanol, including the hydroxyl hydrogen H−H partial structure factor, HHH(K). Because the intermolecular hydroxyl hydrogen separation depends on the length and orientation of the methanol hydrogen bond, HHH(K) provides an indirect probe of the hydrogen bond configuration. We have computed this partial structure factor for methanol at 298 K and 0.786 g/ cm3 (1 atm), and compare it to the experimental result of Yamaguchi et al. in Figure 12. We see that there is fairly good

NVT (0.694 g/cm3)

6.87

NPT (1 atm)

0.726

+E3b ATM




reduces the error in the predicted density to

First-Principles Many-Body Force Fields from the Gas Phase to Liquid

Recommend Documents