Intermolecular Interactions in Complex Liquids

Intermolecular Interactions in Complex Liquids: Effective Fragment Potential Investigation of Water−tert-Butanol Mixtures ABSTRACT: Structure and bonding patterns in tert-butanol (TBA)−water mixtures are investigated by using molecular dynamics simulations with the effective fragment potential (EFP) method. EFP is a model potential in which all parameters are obtained from a set of ab initio calculations on isolated fragment molecules. Mixed-basis EFP potentials (called “EFPm”) for water and TBA molecules were prepared and tested in this work. The accuracy of these EFP potentials is justified by comparison of structures and binding energies in water, TBA, and water−TBA dimers with MP2/6-311+ +G(d,p) data. It has been found that the discrepancies between EFP and MP2 do not exceed 0.1 Å in intermolecular distances and 1 kcal/mol in binding energies. Structures of TBA−water solutions with 0.0, 0.06, 0.11, 0.16, and 0.50 TBA mole fractions were analyzed by using radial distribution functions (RDFs) and coordination numbers. These results suggest that, at low TBA concentrations, the structure of water is enhanced and water and TBA are not homogeneously mixed at the molecular level. In the equimolar TBA−water solution, the microscopic mixing is more complete. Analysis of the energy components in TBA−water solutions shows that, while the electrostatic and exchange-repulsion terms provide the largest contributions to the total potential energy, the relative importance of the polarization and dispersion terms depends on the concentration of TBA. With an increase of TBA concentration, the fraction of the dispersion energy increases, while the fraction of polarization energy diminishes. However, both polarization and dispersion terms are essential for accurate description of these systems. ■ INTRODUCTION Broader use of bioalcohols, i.e., alcohols produced from biomass rather than petroleum sources, could result in energy security and lower emissions of green-house gases.1 Apart from the high production cost, a stumbling block for achieving the extensive use of alcohols is liquid−liquid phase separation in hydrocarbon−alcohol mixtures that becomes even more exacerbated by the presence of water. While different impurities (water, ions, organic molecules) can either induce or inhibit phase separation in hydrocarbon−alcohol systems,2−5 there is little known about the underlying mechanisms of these phenomena, since molecular level studies of these systems are scarce.6 In contrast, there is considerable interest in the nature of hydrophobic and hydrophilic hydration, with many recent controversial findings.7−12 For example, until recently, hydration of alcohols has been interpreted in terms of the Frank− Evans classical “iceberg” model.13 However, recent experimental and theoretical studies provide strong evidence of incomplete mixing at the molecular level and retention of the network structure of bulk water.14−19 The level of mixing, morphology, and other microscopic and macroscopic properties of liquids are governed by the nature of the underlying interactions between molecules in the liquid. Tertiary butanol (TBA) is considered to be a promising candidate for use as a biofuel. TBA is the largest monohydric alcohol that is fully soluble in water. TBA−water systems have been investigated in a number of studies.20−23 Neutron diffraction experiments on small concentrations of TBA20,21 in water suggest that the primary association between alcohol molecules is through interactions of the hydrophobic groups. Even at a high TBA concentration22 (0.86 mol fraction), direct polar interactions between alcohol molecules were not prevalent. At this concentration, a few waters coordinate several TBA molecules, similar to what was reported by Guo et al. for water−methanol solutions.15 Computer simulations of TBA−water mixtures have been performed by Lee and van der Vegt, with results differing from experiments.23 In particular, significant TBA−TBA hydrogen bonding for all concentrations above 0.04 mol fraction alcohol has been observed. In this study, the molecular structure of TBA−water mixtures is investigated by simulations with the general effective fragment potential (EFP),24−26 with particular attention paid to hydrogen bonding. The EFP method is a quantum mechanics based potential that is a computationally inexpensive way of modeling intermolecular interactions in noncovalently bound systems. Absence of fitted parameters and natural partitioning of the interaction energy into electrostatic, polarization, dispersion, and exchange-repulsion terms make it an attractive choice for analysis and interpretation of intermolecular forces. Previously, general EFP has been successfully applied for investigation of the noncovalent interactions in dimers of benzene27 and benzene derivatives,28,29 styrene dimers,30 small water−benzene31 and water−methanol19 complexes, as well as neutral-zwitterionic equilibrium in hydrated alanine.32,33 In the studies of small clusters, EFP provided accurate results as compared to ab initio methods. However, applications of the EFP potential to investigation of properties of bulk liquids are scarce.34,35 In this paper, we report results of EFP molecular dynamics (MD) simulations of TBA−water solutions at different concentrations of TBA, with a goal to characterize structure and hydrogen bonding in these mixtures. For validation of the EFP methodology, the EFP results are compared to classical simulations with the GROMOS96 force field and to experimental data.20,21 While our eventual goal is to understand the molecular-level interactions in complex hydrocarbon− alcohol−water mixtures of relevance to petroleum industry, the present study is an important first step in which we both validate the EFP methodology for description of polar heterogenic systems and provide insights on the structure and hydrogen bonding in intriguing TBA−water solutions of various concentrations. ■ THEORY There are four interaction terms in general EFP (general EFP potential has been originally called EFP2 to be distinguished from the water potential EFP126,36,37), each of which may be thought of as a truncated expansion: Coulomb (electrostatic), induction (polarization), exchange repulsion, and dispersion: E EFP − EFP = Ecoul + E pol + Edisp + Eexrep (1) The terms in the EFP potential may be grouped into longrange, (1/R)n distance dependent, and short-range interactions, which decay exponentially. The Coulomb, induction, and dispersion are long-range interactions, whereas the exchange repulsion and damping terms are short-range. EFP has been described in detail in several papers;24−26 therefore, only a brief overview is presented below. The Coulomb portion of the electrostatic interaction, Ecoul, is obtained using the Stone distributed multipolar analysis (DMA).38,39 This expansion is truncated at the octopole term. The atom centers and bond midpoints are used as expansion points. Classical Coulomb interactions become too repulsive at short range, when the electronic densities of the interacting fragments overlap and charge-penetration effects play a role. To correct for these quantum effects, electrostatic interactions are moderated by screening terms. Several variants of Coulomb screening have been developed;27,40,41 the overlapbased screening41 is used in this work. Induction (polarization), Epol, arises from the interaction of distributed induced dipoles on one fragment with the static multipole field and a field due to induced dipoles on the other fragments. The polarizability expansion is truncated at the first (dipole) term; the molecular polarizability tensor is expressed as a tensor sum of anisotropic localized molecular orbital (LMO) polarizabilities. Therefore, the number of polarizability points is equal to the number of bonds and lone pairs in the system. The induction term is iterated to self-consistency, so it is able to capture leading many body effects.42 Gaussian-type screening is employed to ensure correct short-range behavior and avoid polarization “collapse”.39 Dispersion interactions are expressed by an inverse R expansion, Edisp = ∑ CnR −n n (2) The first term in the expansion, n = 6, corresponds to the induced dipole−induced dipole (van der Waals) interactions. In EFP, coefficients C6 are derived from the (imaginary) frequency dependent polarizabilities integrated over the entire frequency range.43,44 Distributed (centered at LMOs) dynamic polarizability tensors are obtained using the time-dependent Hartree−Fock (HF) method. In addition, the contribution of the n = 8 term is estimated as one-third of the n = 6 term. Overlap-based damping is used to damp the dispersion interactions at short range.41 The exchange repulsion interaction between two fragments is derived as an expansion in the intermolecular overlap, truncated at the quadratic term.45−47 Kinetic and overlap one-electron integrals are calculated between each pair of fragments on-thefly. Thus, each effective fragment carries a basis set and localized wave function. The main drawback of the multipoles obtained from the DMA is their instability with respect to an increase in the size of the basis set. In a large basis set, and especially a basis with diffuse functions, there are many ways to express a given multipole distribution, and the optimal one may be determined variationally. In another basis set, with a different set of exponents, a similar electron density may be optimally represented by a different set of multipoles. This ambiguity increases with the size of the basis set, and it can result in very different and basis-set-dependent multipoles. Strength of the Coulomb interaction becomes somewhat weaker when the basis set is increased. This is consistent with general decrease of the intermolecular binding with increase in the basis. In the EFP model, this in particular results in destabilization of hydrogen bonding. To avoid instabilities in the DMA and weakening of hydrogen bonds, a modified potential with a smaller 6-31+G* basis set48−50 for the DMA part is employed in this work. Parameters for the rest of the energy terms, i.e., polarization, dispersion, and exchangerepulsion, were obtained in 6-311++G(3df,2p).51,52 Additionally, more accurate localized polarizability tensors, obtained with the PBE density functional,53 have been used. Thus, modified EFP potentials (“EFPm”) for water and TBA have been prepared and used throughout this work, along with the standard EFP potentials, in which all terms were generated with the 6-311++G(3df,2p) basis at the Hartree−Fock level (“EFP”). ■ COMPUTATIONAL DETAILS Simulations of mixtures of tertiary butanol (TBA) and water are performed using a classical force-field, GROMOS96,54 and the effective fragment potential (EFP) method. Classical simulations are performed in the GROMACS software package55 and the EFP calculations are carried out in the GAMESS electronic structure program.56,57 In classical simulations, a rigid six-site model was used for TBA and a rigid three-site SPC model was used for water. Intermolecular interactions were calculated as sums of LennardJones terms at CH3 (united methyl atom), C (TBA central carbon), OTBA (TBA oxygen), and OW (water oxygen) as well as Coulomb interactions between static partial charges on C, Figure 1. Lowest energy structures of the water dimer, TBA dimer, and two water−TBA dimers. OTBA, HTBA (TBA carbon, oxygen, and hydrogen), OW, and HW (water oxygen and hydrogen). Simulations were performed for 1 ns with 2208, 1906, 1818, and 864 total molecules for 0.06, 0.11, 0.16, and 0.50 mol fraction of TBA, respectively. The geometries of all molecules were kept rigid using the SHAKE algorithm with a relative tolerance of 10−4. Intermolecular interactions were calculated using a twin-range cutoff of 0.8 and 1.4 nm with nonbonded interactions updated every five time steps. The equations of motion were integrated using a leapfrog algorithm with a time step of 1 fs. A reaction field approximation (εRF) of 64.0, 50.7, 36.1, or 18.0 for 0.06, 0.11, 0.16, and 0.50 mol fraction of TBA was used for longrange interactions. Simulations were performed at 1 atm and 298 K, which were maintained constant using the Berendsen thermostat and barostat58 with coupling times of 0.1 and 0.5 ps, respectively. Two versions of EFP potentials for both TBA and water were prepared. In the standard EFP potentials (called “EFP”), all terms are generated with the 6-311++G(3df,2p) basis at the Hartree−Fock level. To avoid instabilities in DMA and weakening of hydrogen bonds, modified potentials (“EFPm”) with a smaller 6-31+G* basis set for the DMA part were also employed in this work. Parameters for the rest of the energy terms, i.e., polarization, dispersion, and exchange-repulsion, were obtained in 6-311++G(3df,2p). Additionally, more accurate localized polarizability tensors, obtained with the PBE density functional,53 were used for polarization. While preparing different EFP terms with different basis sets and levels of theory may seem to be inconsistent, the approximations beyond each of these terms result in significantly different correlation and basis set dependencies. For example, DMA is not very sensitive to the correlation level but extremely sensitive to the basis set. Static and dynamic polarizabilities monotonically improve with increase of the basis and correlation level. The exchange-repulsion term derived with the largest number of approximations is assumed to be obtained in an infinite (or at least very large) basis set, but its basis set dependence is nonmonotonic. The correlation dependence of the exchange-repulsion term in EFP is not well investigated. Thus, in this work, we chose a route of preparing the best and technically obtainable set of EFP parameters for the sake of the parameter set that is completely consistent in terms of basis and level of correlation. In the future, we plan to develop EFP potentials that are both more accurate and consistent, e.g., potentials obtained completely at the coupled cluster level of theory. MD simulations with the original EFP potential were run at 0.50 mol fraction TBA in water and neat water for 100 ps, with each molecule treated as an EFP fragment. EFPm-MD simulations were performed at 0.00, 0.06, 0.11, 0.16, and 0.50 mol fraction of TBA in water. Five different 25 ps MD simulations were run for each TBA concentration, and the resulting radial distribution functions (RDFs) were averaged. Single 100 ps MD calculations were also performed, and similar RDFs were obtained (data not shown). Simulations at 0.06 mol fraction TBA contained 150 total molecules, while those at 0.11 and 0.16 mol fraction TBA each contained 100 total molecules and those at 0.50 mol fraction TBA contained 98 total molecules. Simulations of neat water were performed with 64 water molecules. The equations of motion were integrated using the velocity Verlet algorithm with a time step of 0.5 fs. Simulations were performed in the NVT ensemble at 300 K; the temperature was maintained constant using the Nose− Hoover chain algorithm.59,60 Initial velocities were generated from a random Maxwell−Boltzmann ensemble. Periodic boundary conditions and switching functions to gradually reduce interactions to zero were used for all simulations. Initial geometries for the EFP simulations were generated by performing calculations with a leapfrog stochastic dynamics integrator with the GROMOS96 force field, followed by geometry optimization with EFP. It was observed that the geometries prepared in this way required shorter equilibration times. To ensure that the size of the systems used in the EFP simulations is sufficient, we performed several GROMOS96 simulations for systems of the same size as in EFP and compared radial distribution functions in smaller and larger systems. No quantitative differences were observed. Table 1. Interaction Energies (kcal/mol) and Representative Distances (Å) in the Water Dimer (2W), TBA Dimer (2TBA), and Water−TBA Dimers (1W1TBA) (see Figure 1 for Structures)a interaction energy a energy geometry 2W EFP EFPm MP2 EFP EFPm EFP MP2 MP2 −4.52 −5.88 −5.85 (−4.50) −5.18 −6.21 (−4.49) −5.35 2TBA −6.47 −7.90 −7.98 (−5.68) −6.83 −8.65 (−6.13) −7.39 O−O distance 1W1TBA A 1W1TBA B −6.10 −7.41 −7.08 (−5.39) −6.24 −7.53 (−5.63) −6.58 −4.07 −5.42 −5.97 (−4.31) −5.14 −6.51 (−4.33) −5.42 geometry 2W 2TBA 1W1TBA A 1W1TBA B EFP EFPm MP2 3.02 2.96 2.92 3.00 2.91 2.86 2.95 2.87 2.85 3.10 3.03 2.94 BSSE corrected energies are reported in parentheses; averages between BSSE-corrected and non-corrected values are given in italics. Table 2. EFP and EFPm Energy Components as Total Values (kcal/mol) and as a Fraction of the Total EFP/EFPm Energy for the Water Dimer (2W), TBA Dimer (2TBA), and Water−TBA Dimers (1W1TBA) EFP energy components EFP 2W electrostatic polarization dispersion exchange-repulsion total −6.59 −0.96 −1.17 4.20 −4.52 1W1TBA A EFPm 2W 1W1TBA A electrostatic polarization dispersion exchange-repulsion total −8.51 −1.16 −1.34 5.14 −5.88 −10.05 −1.94 −3.47 8.05 −7.41 1W1TBA B 2TBA 2W −7.43 −1.24 −4.55 6.75 −6.47 1.46 0.21 0.26 −0.93 1.00 1W1TBA B 2TBA −7.17 −1.11 −2.13 4.99 −5.41 −8.90 −1.79 −5.94 8.74 −7.89 −8.18 −5.39 −1.42 −0.78 −2.88 −1.72 6.38 3.82 −6.10 −4.07 EFPm energy components Dimer Calculations. Additionally, water, TBA, and mixed water−TBA dimers were optimized using the temperature annealing Monte Carlo protocol with the EFP method. Moller−Plesset second order perturbation theory (MP2) calculations in the 6-311++G(d,p) basis were performed at the lowest found structures. ■ fraction of total EFP energy RESULTS AND DISCUSSION Dimer Calculations. Two goals are pursued by investigating the structures and interaction energies of water, TBA, and mixed water−TBA dimers. First of all, one can compare relative energetics in the dimers with preferred interactions in the bulk solutions (based on their structural characteristics). This comparison provides a way to characterize the importance of many-body cooperative effects in the bulk. An additional purpose of the dimer calculations is to benchmark EFP against MP2 that is known to accurately describe hydrogen-bonded complexes.61−66 The lowest energy structures of water and TBA dimers and the two lowest structures of water−TBA dimers are shown in Figure 1. The structures were obtained with temperatureannealing Monte Carlo simulations using the EFP potential. Further, these structures were reoptimized with the modified EFP potential (EFPm) or MP2. For consistency with structural data of simulations in liquid, oxygen and hydrogen in water are called OW and HW and oxygen and hydroxyl-hydrogen in TBA 1W1TBA A 1W1TBA B 2TBA 1.34 1.32 0.23 0.19 0.47 0.42 −1.05 −0.94 1.00 1.00 fraction of total EFPm energy 1.15 0.19 0.70 −1.04 1.00 2W 1W1TBA A 1W1TBA B 2TBA 1.45 0.20 0.23 −0.87 1.00 1.36 0.26 0.47 −1.09 1.00 1.32 0.20 0.39 −0.92 1.00 1.13 0.23 0.75 −1.11 1.00 are called OTBA and HTBA. Then, each type of dimer can be characterized by the formed hydrogen bond, such as OW−HW (hydrogen bond in water dimer), OTBA−HTBA (hydrogen bond in TBA dimer), etc. Interaction energies and representative distances in the dimers are summarized in Table 1. Both types of EFP potentials and MP2 place the TBA dimer (OTBA−HTBA) as the lowest in energy, followed by the TBA−water dimer in which TBA serves as a proton acceptor (OTBA−HW). The interaction energies in the water dimer (OW−HW) and the second water−TBA dimer (with TBA being a proton donor, OW−HTBA) are very similar. At the EFP geometries, both EFP and MP2 predict lower energy in the water dimer, while the water−TBA dimer is more stable by a few tenths of kcal/mol when the MP2 geometries are used. MP2 interaction energies for the dimers are not completely converged in the 6-311++G(d,p) basis set. Therefore, basis set superposition error (BSSE) corrections were employed (shown in Table 1 in parentheses) as well as average values between BSSE-corrected and noncorrected energies were computed. We believe the averages between corrected and noncorrected values represent the most accurate estimates of the interaction energies in the dimers;31 we use these values as a reference in the following discussion. Compared to MP2, the original EFP potential generally overestimates intermolecular separations by 0.1−0.2 Å and Table 3. Comparison of EFP and EFPm Energy Components with the RVS Energy Components (Energies Are in kcal/mol) EFP energy components a EFP geometry 2W CEX (Elec + Exrep) Pol (+ CTa) CEX + Pol (+ CTa) dispersion total EFP −2.39 −0.96 −3.35 −1.17 −4.52 1W1TBA A 1W1TBA B −1.80 −1.57 −1.42 −0.78 −3.22 −2.35 −2.88 −1.72 −6.10 −4.07 EFPm energy components RVS energy components 2TBA 2W 1W1TBA A 1W1TBA B 2TBA ΔERVS‑EFPb −0.68 −1.24 −1.92 −4.55 −6.47 −3.18 −1.23 −4.41 −2.17 −2.12 −4.29 −2.66 −1.17 −3.83 −1.29 −2.13 −3.42 −0.71 −0.56 −1.28 RVS energy components EFPm geometry 2W 1W1TBA A 1W1TBA B 2TBA 2W 1W1TBA A 1W1TBA B 2TBA ΔERVS‑EFPmb CEX (Elec + Exrep) Pol (+ CTa) CEX + Pol (+ CTa) dispersion total EFPm −3.38 −1.16 −4.54 −1.34 −5.88 −2.00 −1.94 −3.94 −3.47 −7.41 −2.18 −1.11 −3.28 −2.13 −5.41 −0.16 −1.79 −1.95 −5.94 −7.89 −2.95 −1.46 −4.41 −1.78 −2.53 −4.31 −2.36 −1.40 −3.76 −0.68 −2.58 −3.26 −0.01 −0.49 −0.51 Charge-transfer (CT) term is included in RVS but not in EFP. bAverage energy difference between RVS and EFP or EFPm energy components. Figure 2. Oxygen−oxygen (OW−OW) RDFs of liquid water by EFP (black), EFPm (blue), GROMOS96 (green), and X-ray scattering experiment73 (red). underestimates the interaction energies by 0.5−1.3 kcal/mol. This observation was one of the reasons for modifying the EFP potential and strengthening its Coulomb and polarization interactions. Indeed, the modified EFP potential produces much better results. For example, intermolecular distances in EFPm and MP2 differ by less than 0.1 Å and the stabilities of the dimers by EFPm are within 0.8 kcal/mol of the MP2 stabilities, with the average differences not exceeding 0.5 kcal/ mol. Predictive description of complex liquids requires accurate relative energetics between different types of intermolecular interactions. Analysis of interaction energies in the dimers reveals that, compared to MP2, both original and modified EFPs overestimate the stability of the first (OTBA−HW) with respect to the second (OW−HTBA) water−TBA dimer (∼2.0 and ∼1.2 kcal/mol energy differences by EFP and MP2, respectively). As mentioned earlier, underestimation of the binding energy in the OW−HTBA water−TBA dimer by EFP also affected the relative energetics of water−water and water− TBA dimers. Compared to MP2, EFP potentials also underestimate the stability of the first (OTBA−HW) water− TBA dimer with respect to the TBA−TBA (OTBA−HTBA) dimer (∼0.4−0.5 kcal/mol by EFP versus ∼0.8 kcal/mol by MP2). Other differences in relative interaction energies between MP2 and EFP do not exceed ∼0.3−0.4 kcal/mol. A key feature of the EFP method is a decomposition of the total intermolecular energy into the components, shown in Table 2. While the magnitudes of the energy components differ between EFP and EFPm (partly because of slightly different geometries of EFP and EFPm), their fractions with respect to the total energy are very similar between the two potentials. As expected, electrostatic and repulsion energies contribute the most to the total hydrogen bonding energy for all of the dimers. The electrostatic energy fraction decreases in the order of water Figure 3. Hydrogen bonding RDFs (OHW, OH, OWHW, and HOW) by GROMOS96 at the following molar fractions of TBA in water: (a) XTBA = 0.06, (b) XTBA = 0.11, (c) XTBA = 0.16, and (d) XTBA = 0.50. O is TBA oxygen, H is TBA hydroxyl hydrogen, OW is water oxygen, and HW is water hydrogen. dimer, water−TBA, and TBA dimer, which is compensated by an increase of dispersive interactions in TBA containing dimers. The polarization energy fraction remains fairly constant among all of the dimers. Comparison of the EFP and EFPm energy components with the interaction energy decomposition by the reduced variational space method (RVS)67 is shown in Table 3. The original EFP underestimates the magnitude of the electrostatic term, as follows from comparison to the CEX (electrostatic plus exchange-repulsion) energies by RVS. The discrepancies between EFP and RVS in different dimers range from 0.4 to 1.1 kcal/mol, with RVS values always being more attractive. A combination of polarization and charge-transfer terms in RVS is also on average ∼0.5 kcal/mol more binding than the EFP polarization term. As a result, a sum of EFP electrostatic, exchange-repulsion, and polarization energies underestimates the RVS total energies for all dimers, with the average discrepancy being 1.3 kcal/mol. EFPm provides systematic improvement in description of the electrostatic term and a slight improvement in accuracy of the polarization term. As a result, there is a very good agreement between CEX components of EFPm and RVS, with specific errors not exceeding 0.5 kcal/mol. The agreement between the total RVS energies and corresponding EFPm components is much improved, with RVS energies being lower by 0.5 kcal/ mol on average. This analysis is consistent with the results of our investigation of accuracy of EFP on the S22 data set, where it has been shown that the modified EFP potentials provide an overall accurate description of intermolecular interactions of different types.68 On the basis of the total EFP energy of the dimers, one may expect that, in an equal mixture of TBA and water, the hydrogen bonding would form more preferably between TBA molecules than between water−TBA or water−water. In more detail, the EFP data for dimers suggest that the priority of hydrogen bonds (from highest to lowest) should be OTBA− HTBA, OTBA−HW, OW−HW, HTBA−OW, while the MP2 data suggest the order of OTBA−HTBA, OTBA−HW, HTBA−OW, OW− HW. In the next section, these predictions will be compared to the results of bulk simulations of water−TBA mixtures. However, it should be noted that the dimer energies cannot be used as the only guide in predicting the bulk structures, since many-body effects, hydrophobic interactions between TBA methyl groups, and entropic factors may and do change the picture significantly. Bulk Water. Before proceeding with water−TBA bulk simulations, we first analyze the quality of the EFP potentials in predicting the structure of bulk water. Figure 2 shows a Figure 4. OWO, CC, OWOW, and COW RDFs by GROMOS96 at the following molar fractions of TBA in water: (a) XTBA = 0.06, (b) XTBA = 0.11, (c) XTBA = 0.16, and (d) XTBA = 0.50. O is TBA oxygen, C is TBA central carbon, and OW is water oxygen. comparison of the oxygen−oxygen radial distribution functions (RDF) for liquid water by the original EFP and modified EFPm potentials. While EFPm showed only modest improvements when comparing H-bonds in dimers, it surpasses the original EFP potential in reproducing the structure of bulk water. The first peak in the OW−OW RDF by the original EFP is shifted by 0.2 Å, and the structure beyond the first shell is completely missing. This might be explained by too weak hydrogen bonding in EFP water. The modified EFPm potential provides a reasonable agreement with experimental RDFs, slightly overestimating the height of the first peak (i.e., the EFPm water is slightly overstructured) and shifting its position to the longer distances (from the experimental maximum at 2.75 to 2.79 Ǻ ). On the basis of the performance of the original and modified EFP potentials for bulk water, we exclusively use the modified EFPm potential for analysis of water−TBA mixtures. Figures 3−6 show GROMOS96 and EFPm RDFs in solutions with low concentration of TBA in water (0.06, 0.11, and 0.16 molar fractions of TBA (XTBA)) and in the equimolar mixture. Bowron, Finney, and Soper report RDFs for the low concentrations of TBA from neutron diffraction experiments.20 Overall, RDF peaks obtained from EFPm are shifted to longer distances compared to those from the classical force field and experiment. In particular, optimal hydrogen bonding distances for TBA−TBA interactions, i.e., OTBA−HTBA, are ∼0.4 Å longer in EFPm than in GROMOS96. This discrepancy occurs because GROMOS96 predicts the first maxima in the OTBA− HTBA RDF to be ∼0.2 Å too short compared to experiment, while EFPm overestimates the optimal O−H length by the same amount. The discrepancy in OTBA−HTBA distances is somewhat consistent with EFPm slightly overestimating intermolecular separations in the dimers. GROMOS96 simulations show an OTBA−HTBA peak corresponding to a Hbonded TBA−TBA complex at all concentrations of TBA. However, EFPm RDFs at the lowest concentration of TBA (0.06 mol fraction) do not have a corresponding peak. The first peak of the C−C RDF (C is the central hydrogen on TBA) is also absent at this concentration, consistent with no OTBA− HTBA hydrogen bonding being observed. While the absence of this peak may be attributed to a lack of sampling in EFPm simulations at low TBA concentrations, it is consistent with experimental observations, suggesting little TBA−TBA hydrogen bonding.20,21 As with the dimer case, we obtained the energy decomposition of the bulk systems from the EFP method (Table 4). In agreement with the dimers, the electrostatic and repulsion energies are the largest energy terms in the bulk. However, a sum of electrostatic and exchange-repulsion energies provides only a small percent to the total interaction energy, leaving polarization and dispersion to dominate the Figure 5. Hydrogen bonding RDFs (OHW, OH, OWHW, and HOW) by EFPm at the following molar fractions of TBA in water: (a) XTBA = 0.06, (b) XTBA = 0.11, (c) XTBA = 0.16, and (d) XTBA = 0.50. O is TBA oxygen, H is TBA hydroxyl hydrogen, OW is water oxygen, and HW is water hydrogen. bonding patterns in water−TBA solutions. With increasing TBA concentration, the electrostatic energy and polarization energy fraction decrease, while the dispersion energy fraction increases. These results are consistent with what one would expect when replacing highly polar water molecules with TBA molecules containing hydrophobic methyl groups. The repulsion energy fraction remains fairly constant at all TBA concentrations. It is noteworthy to mention the higher fraction (and importance) of polarization in the bulk compared to that in the dimers. For example, at all TBA concentrations, the polarization term is larger in magnitude than the sum of the first-order (electrostatic plus exchange) contributions. A similar situation was previously observed in water clusters.69,70 The increase of polarization in the bulk occurs due to many-body cooperative behavior in polar medium; the relative decrease of the polarization at higher concentrations of TBA is a sign of weakening these many-body effects in the presence of the TBA hydrophobic groups. It is worth mentioning that the accurate description of cooperativity is nontrivial; it is likely that standard polarizable force-fields with atom-centered scalar polarizabilities underestimate the amplitude of cooperativity and the magnitude of polarization energy.71,72 In addition to comparing average intermolecular distances and energy components, we calculated the coordination and excess coordination numbers for different concentrations of TBA in water. With the present notations, the coordination number (CN) is the average number of a particular type of atom j from one molecule within a specified distance of an atom i of another molecule: R Ni , j = 4πρj g (r )r 2 d r 0 i,j ∫ (3) where ρj = Nj/V is the number density of atom j, g(r) is the radial distribution function, and R is the cutoff distance. Coordination numbers depend on the order of atoms i and j; i.e., Ni,j and Nj,i are different quantities. For example, the C−OW CN determines the average number of water oxygens found in proximity to TBA central carbons. However, the OW−C CN shows the average number of TBA central carbons near water oxygen. Excess coordination numbers (ECN) represent the excess or deficiency of the number of particular atoms of one molecule around an atom of another molecule relative to that number in an ideal liquid: R Nex i , j = 4πρj [gi , j(r ) − 1]r 2 dr 0 ∫ (4) The coordination and excess coordination numbers are calculated by integrating the appropriate RDF to the first minimum except for CC, which was integrated up to 7.8 Å. Coordination numbers provide information on structural preferences; however, their analysis in complex mixtures is Figure 6. OWO, CC, OWOW, and COW RDFs by EFPm at the following molar fractions of TBA in water: (a) XTBA = 0.06, (b) XTBA = 0.11, (c) XTBA = 0.16, and (d) XTBA = 0.50. O is TBA oxygen, C is TBA central carbon, and OW is water oxygen. Table 4. EFPm Energy Components as a Fraction of Total EFPm Energy of the Bulk Water−TBA Solutions at Different TBA Mole Fractions TBA mole fraction energy component 0.00 0.06 0.11 0.16 0.50 electrostatic exchange-repulsion CEX (Elec + Exrep) polarization dispersion 1.29 −1.24 0.05 0.56 0.38 1.28 −1.24 0.04 0.50 0.45 1.26 −1.23 0.03 0.49 0.49 1.25 −1.31 −0.06 0.49 0.56 1.03 −1.20 −0.17 0.37 0.80 intricate, since sterical effects play an important role. For example, with increasing mole fraction of TBA, it is expected that the water−water OWOW and TBA−water OTBAOW CNs should decrease (since less sites are available for water molecules to surround a particular water or TBA molecule), while the water−TBA OWOTBA CN should increase. The EFPm and GROMOS96 simulations yield similar CNs and ECNs at all concentrations (see Table 5). The CNs from simulations are also in reasonable agreement with available experimental data. The largest discrepancies between EFPm and GROMOS96 are observed for the smallest TBA concentration (e.g., CC and COW coordination numbers differ by 1.2). CC and COw CNs provide information on the amount of TBA−TBA and TBA−water hydrogen-bonded pairs, respectively. As was observed from the analysis of RDFs, EFPm and GROMOS96 disagree in the character of TBA− TBA interactions in the lowest concentration solution. Differences in the CC coordination numbers confirm this conclusion. The CC coordination numbers also show the largest differences with experiment. GROMOS96 overestimates CC CNs (i.e., the amount of TBA−TBA hydrogen bonding) for all considered concentrations; EFPm is in agreement with experiment for the lowest concentration but also provides too high CC CNs for two larger concentrations. The EFPm CNs and particularly ECNs for OWHW and OWOW indicate that the structure of water in TBA solutions up to 0.16 mol fraction TBA is enhanced and strengthened. For example, the OWOW CNs remain fairly constant, while the OWOW ECNs increase from −0.6 to 1.2 when the TBA mole fraction increases from 0 to 0.16. While the CC coordination and excess coordination numbers may have large uncertainties in EFP (due to possibly insufficient sampling when only a few TBA molecules are present in a solution), they may suggest that the TBA−TBA interactions are also preferred. This is in agreement with the intuitive prediction that partly hydrophobic TBA molecules prefer to build clusters in an excess of polar water. On the other hand, interactions of water with TBA (OWC and OWOTBA) are slightly weakened at the low TBA concentrations. Very similar values of OWOTBA and OWC coordination numbers suggest that TBA and water always Table 5. Coordination and Excess Coordination Numbers for Solutions with Varying Mole Fractions of TBA (0.00, 0.06, 011, 0.16, 0.50)a coordination number 0.00 OWHW OTBAHW 2.1 HTBAOW OTBAHTBA OWOW OWOTBA OTBAOW COW OWC CC 4.4 0.06 1.9 1.8 1.4 1.0 0.8 0.0 0.1 4.7 4.3 0.2 3.1 b (1.8) (1.6) (0.9) (0.1) (4.4) (0.2) (2.7) 3.8 (2.6) 2.5 0.2 (0.2) 2.5 (3.7) 2.9 0.11 c 0.16 excess coordination number c 0.00 0.06 0.11 0.16 −2.2 −1.5 (−1.5) −1.9 (−1.8) −1.0 (−0.5) −1.8 (−0.9) −0.7 (−0.5) −1.2 (−0.9) 0.3 (0.3) 0.0 (0.1) 0.5 (0.5) −0.9 (−0.8) −0.8 (−0.5) −0.6 (−0.3) 0.1 (0.2) (0.2) 0.3 (0.5) −0.1 (0.0) −0.1 (−0.1) −0.2 (0.0) −0.2 (0.1) (4.0) 2.0 (2.5) −0.1 (0.1) 0.7 (0.7) 1.2 (1.1) 1.0 (1.6) (0.3) (2.4) 0.5 (0.4) 2.4 (2.3) 1.5 (1.5) 1.5 (1.5) −0.1 (−0.1) −1.7 (−1.8) −0.2 (−0.1) −1.5 (−1.1) −0.2 (−0.1) −1.0 (−0.7) 0.5 (0.5) 0.5 (0.5) (2.2) 2.7 (2.1) 2.0 0.5 (0.4) 7.6 (7.2) 5.8 1.6 (1.4) −4.8 (−4.7) −4.1 (−3.4) −3.3 (−2.7) −0.2 (−0.2) 1.6 (1.4) 10.7 (10.3) −0.3 (−0.3) −0.8 (0.7) −0.5 (−0.4) 0.9 (1.0) −0.6 (−0.5) 0.8 (0.8) −0.2 (−0.2) −0.1 (−0.2) 1.8 1.4 1.3 0.9 0.9 0.1 0.1 4.4 (2.1) (1.7) (4.0) 1.7 1.4 1.2 0.7 0.9 0.1 0.1 4.5 0.3 2.4 2.3 2.6 2.2 0.3 6.0 4.4 (0.8) (0.1) (0.3) (5.9) 0.50 (1.7) (1.3) 1.1 (1.0) 0.9 (0.8) (0.8) −0.6 0.50 a EFPm, GROMOS96 (in parentheses), and experimental (in italic) values are shown. Ow and Hw are water oxygen and hydrogen, respectively; OTBA and HTBA are TBA oxygen and hydroxyl hydrogen, and C is the central carbon on TBA. bExperimental values from ref 21. cExperimental values from ref 20. interact via a close contact of oxygen sites; i.e., they always form a hydrogen bond. Comparison of the OWOW and OWOTBA CNs and ECNs reveals that the water−water interactions are favored over water−TBA interactions across all concentrations. Additionally, comparison of OTBAHW and OTBAHTBA CNs and ECNs suggests that TBA−TBA hydrogen bonding is somewhat preferred over TBA−water bonding for low concentrations of TBA. This preference is more evident in GROMOS96. Comparing this pair of CNs, one needs to remember that, since water has two hydrogens HW and TBA only one hydroxylhydrogen HTBA, the OTBAHW and OTBAHTBA CNs should be scaled correspondingly. However, the latter conclusion on the preference of TBA−TBA H-bonding over TBA−water Hbonding should be considered with care due to overall small values of the OTBAHTBA coordination numbers. The above analysis suggests that, at the low concentrations of TBA (up to 0.16 mol fraction), water−water interactions and H-bonding as well as TBA−TBA interactions are enhanced, while water−TBA interactions are less preferable. This means that the solutions are microscopically heterogeneous and form TBA and water microphases. The many-body, hydrophobic, and entropic effects make the relative energetics of specific pairwise H-bonding (as was analyzed in the dimers) of little relevance. For example, the water−water H-bond was found to be one of the weakest in the dimers, but based on CNs, it is the most preferred in the low-concentrated TBA solutions. The bonding pattern is somewhat different in the equimolar TBA−water mixture: both water−TBA (OWC and OWOTBA) and TBA−water (OTBAOW and COW) excess coordination numbers increase at this concentration. In agreement with this, TBA−water H-bonding (OTBAHW) becomes more favorable than the TBA−TBA (OTBAHTBA) H-bonding in EFPm but remains less favorable in GROMOS96. In addition, TBA−TBA interactions (CC ECNs) become less preferable in the equimolar mixture. The enhancement of the water structure also becomes less obvious in the equimolar mixture. The picture emerging from this analysis is as follows: the equimolar solution experiences a more homogeneous mixing than the lowconcentrated TBA solutions, and more interactions occur between water and TBA molecules. ■ CONCLUSIONS The water−TBA clusters and bulk solutions are investigated in this paper with the ab initio-based EFP and classical GROMAS96 potentials. We prepared the modified EFP potentials for TBA and water that are better suited for description of H-bonding and showed that they accurately reproduce the experimental structure of pure water and water− TBA mixtures. We recommend using a protocol for preparing EFP potentials developed in this work for other H-bonding systems. Benchmarks of the modified EFP potential on water, TBA, and water−TBA dimers show that the energetics of H-bonding in these systems is accurately captured by EFPm, as compared to MP2. However, the relative strengths of H-bonds in the dimers cannot be unambiguously correlated with the preferences in H-bonding in liquids, since many-body and hydrophobic interactions play an important role there. In simulations of pure water, the EFPm potential reproduces experimental radial distribution functions very well, both in peak location and shape. Comparison of EFPm and GROMOS96 simulations for TBA solutions shows some discrepancies in the positions of the first peaks of H-bonding RDFs, since EFPm slightly overestimates the preferred length of H-bond (as compared to experimental RDFs), while GROMOS96 underestimates it. We observed that the EFP reproduces the experimental structural data for the lowest concentration of TBA (0.06 mol fraction) but overestimates TBA−TBA hydrogen bonding at higher concentrations. The classical simulations with GROMOS96 overestimate TBA− TBA hydrogen bonding at all concentrations. The observed discrepancies suggest that the hydrogen-bonding patterns in mixed water−alcohol systems are very sensitive to a model potential, with the last word yet to be said. Additionally, we found that the water−water interactions are promoted in the presence of a small amount of TBA and the solutions with small concentrations of TBA experience a high degree of heterogeneity. The equimolar solution is more uniformly mixed at the molecular level. The EFP method allows one to examine the magnitude of contributions to the total energy of a system. While both for the dimers and bulk systems, the electrostatic and repulsion components are the two largest components of the interaction energy, the relative weight of polarization increases from the dimers to the bulk. Indeed, the polarization energy contribution in the bulk is larger than a sum of the electrostatic and exchange-repulsion energies at all concentrations of TBA, suggesting importance of coopertivity in water−TBA solutions. Dispersion energy becomes significant at high concentrations of TBA, while electrostatic and polarization energies decrease with increasing TBA concentration. Our future work will focus on further improving the accuracy and efficiency of the EFP potential, in particular, its chargetransfer and short-range charge-penetration terms. We will also continue to unravel mysteries of hydrogen-bonded systems, including more complex ternary and tertiary solutions with hydrocarbons and/or ions, and relate microscopic properties such as differences in hydrogen bonding to macroscopic phenomena such as phase separation. ■ REFERENCES (1) Keasling, J. D.; Chou, H. Metabolic engineering delivers nextgeneration biofuels. Nat. Biotechnol. 2008, 26, 298−299. (2) Peschke, N.; Sandler, S. Liquid-Liquid Equilibria of Fuel Oxygenate Plus Water Plus Hydrocarbon Mixtures 0.1. J. Chem. Eng. Data 1995, 40, 315−320. (3) Momo Jeufack, H.; Suter, D. Effect of different sodium halides on the self-association of tertiary butanol in water. J. Chem. Phys. 2007, 126, 144501. (4) Bayrak, Y.; Iscan, M. Studies on the phase behavior of the system non-ionic surfactant/alcohol/alkane/H2O. Colloids Surf., A 2005, 268, 99−103. (5) Li, X.; Ueda, K.; Kumieda, H. Solubilization and phase behavior of microemulsions with mixed anionic-cationic surfactants and hexanol. Langmuir 1999, 15, 7973−7979. (6) Lee, M.-E.; Van der Vegt, N. F. A. Molecular Thermodynamics of Methane Solvation in tert-Butanol−Water Mixtures. J. Chem. Theory Comput. 2007, 3, 194−200. (7) Chandler, D. Interfaces and the driving force of hydrophobic assembly. Nature 2005, 437, 640−647. (8) Jensen, M. O.; Mouritsen, O. G.; Peters, G. H. The hydrophobic effect: Molecular dynamics simulations of water confined between extended hydrophobic and hydrophilic surfaces. J. Chem. Phys. 2004, 120, 9729. (9) Wojtkow, D.; Czarnecki, M. A. Two-Dimensional Attenuated Total Reflection Infrared and Near-Infrared Correlation Study of the Structure of Butyl Alcohol/Water Mixtures. Appl. Spectrosc. 2007, 61, 928−934. (10) Ferrari, E. S.; Burton, R. C.; Davey, R. J.; Gavezzotti, A. Simulation of phase separation in alcohol/water mixtures using twobody force field and standard molecular dynamics. J. Comput. Chem. 2006, 27, 1211−1219. (11) Sung, J.; Park, K.; Kim, D. Surfaces of Alcohol−Water Mixtures Studied by Sum-Frequency Generation Vibrational Spectroscopy. J. Phys. Chem. B 2005, 109, 18507−18514. (12) Asprion, N.; Hasse, H.; Maurer, G. Thermodynamic and IR spectroscopic studies of solutions with simultaneous association and solvation. Fluid Phase Equilib. 2003, 208, 23−51. (13) Frank, H. S.; Evans, M. W. Free Volume and Entropy in Condensed Systems III. Entropy in Binary Liquid Mixtures; Partial Molal Entropy in Dilute Solutions; Structure and Thermodynamics in Aqueous Electrolytes. J. Chem. Phys. 1945, 13, 507−532. (14) Dixit, S.; Crain, J.; Poon, W. C. K.; Finney, J. L.; Soper, A. K. Molecular segregation observed in a concentrated alcohol-water solution. Nature 2002, 416, 829−832. (15) Guo, J.−H.; et al. Molecular Structure of Alcohol-Water Mixtures. Phys. Rev. Lett. 2003, 91, 157401. (16) Dougan, L.; et al. Methanol-water solutions: A bi-percolating liquid mixture. J. Chem. Phys. 2004, 121, 6456. (17) Noskov, S. Y.; Lamoureux, G.; Roux, B. Molecular Dynamics Study of Hydration in Ethanol−Water Mixtures Using a Polarizable Force Field. J. Phys. Chem. B 2005, 109, 6705−6713. (18) Allison, S. K.; Fox, J. P.; Hargreaves, R.; Bates, S. P. Clustering and microimmiscibility in alcohol-water mixtures: Evidence from molecular-dynamics simulations. Phys. Rev. B 2005, 71, 024201. (19) Adamovic, I.; Gordon, M. S. Methanol−Water Mixtures: A Microsolvation Study Using the Effective Fragment Potential Method. J. Phys. Chem. A 2006, 110, 10267−10273. (20) Bowron, D. T.; Finney, J. L.; Soper, A. K. Structural Investigation of Solute−Solute Interactions in Aqueous Solutions of Tertiary Butanol. J. Phys. Chem. B 1998, 102, 3551−3563. (21) Bowron, D. T.; Soper, A. K.; Finney, J. L. Temperature dependence of the structure of a 0.06 mol fraction tertiary butanolwater solution. J. Chem. Phys. 2001, 114, 6203. (22) Bowron, D. T.; Moreno, S. D. The structure of a concentrated aqueous solution of tertiary butanol: Water pockets and resulting perturbations. J. Chem. Phys. 2002, 117, 3753. (23) Lee, M. E.; van der Vegt, N. F. A. A new force field for atomistic simulations of aqueous tertiary butanol solutions. J. Chem. Phys. 2005, 122, 114509. (24) Gordon, M. S.; Slipchenko, L. V.; Li, H.; Jensen, J. H. The Effective Fragment Potential: A General Method for Predicting Intermolecular Interactions. Ann. Rep. Comp. Chem. 2007, 3, 177−193. (25) Ghosh, D.; et al. Noncovalent Interactions in Extended Systems Described by the Effective Fragment Potential Method: Theory and Application to Nucleobase Oligomers. J. Phys. Chem. A 2010, 114, 12739−12754. (26) Gordon, M. S.; et al. The Effective Fragment Potential Method: A QM-Based MM Approach to Modeling Environmental Effects in Chemistry. J. Phys. Chem. A 2001, 105, 293−307. (27) Slipchenko, L. V.; Gordon, M. S. Electrostatic energy in the effective fragment potential method: Theory and application to benzene dimer. J. Comput. Chem. 2007, 28, 276−291. (28) Smith, Q. A.; Gordon, M. S.; Slipchenko, L. V. Benzene− Pyridine Interactions Predicted by the Effective Fragment Potential Method. J. Phys. Chem. A 2011, 115, 4598−4609. (29) Smith, T.; Slipchenko, L. V.; Gordon, M. S. Modeling π−π Interactions with the Effective Fragment Potential Method: The Benzene Dimer and Substituents. J. Phys. Chem. A 2008, 112, 5286− 5294. (30) Adamovic, I.; Li, H.; Lamm, M. H.; Gordon, M. S. Modeling Styrene−Styrene Interactions. J. Phys. Chem. A 2006, 110, 519−525. (31) Slipchenko, L. V.; Gordon, M. S. Water−Benzene Interactions: An Effective Fragment Potential and Correlated Quantum Chemistry Study. J. Phys. Chem. A 2009, 113, 2092−2102. (32) Mullin, J. M.; Gordon, M. S. Alanine: Then There Was Water. J. Phys. Chem. B 2009, 113, 8657−8669. (33) Mullin, J. M.; Gordon, M. S. Water and Alanine: From Puddles(32) to Ponds(49). J. Phys. Chem. B 2009, 113, 14413−14420. (34) Netzloff, H. M.; Gordon, M. S. The effective fragment potential: Small clusters and radial distribution functions. J. Chem. Phys. 2004, 121, 2711. (35) Pranami, G.; Slipchenko, L.; Lamm, M. H.; Gordon, M. S. Coarse-Grained Intermolecular Potentials Derived From The Effective Fragment Potential: Application To Water, Benzene, and Carbon Tetrachloride. Multi-Scale Quantum Models for Biocatalysis; Springer: The Netherlands, 2009; pp 197−218. (36) Adamovic, I.; Freitag, M. A.; Gordon, M. S. Density functional theory based effective fragment potential method. J. Chem. Phys. 2003, 118, 6725. (37) Day, P. N.; et al. An effective fragment method for modeling solvent effects in quantum mechanical calculations. J. Chem. Phys. 1996, 105, 1968. (38) Stone, A. J. Distributed multipole analysis, or how to describe a molecular charge distribution. Chem. Phys. Lett. 1981, 83, 233−239. (39) Stone, A. J. The Theory of Intermolecular Forces; Oxford University Press: Oxford, U.K., 1996. (40) Freitag, M. A.; Gordon, M. S.; Jensen, J. H.; Stevens, W. J. Evaluation of charge penetration between distributed multipolar expansions. J. Chem. Phys. 2000, 112, 7300. (41) Slipchenko, L. V.; Gordon, M. S. Damping functions in the effective fragment potential method. Mol. Phys. 2009, 107, 999. (42) Li, H.; Netzloff, H. M.; Gordon, M. S. Gradients of the polarization energy in the effective fragment potential method. J. Chem. Phys. 2006, 125, 194103. (43) Adamovic, I.; Gordon, M. S. Dynamic polarizability, dispersion coefficient C6 and dispersion energy in the effective fragment potential method. Mol. Phys. 2005, 103, 379. (44) Amos, R. D.; Handy, N. C.; Knowles, P. J.; Rice, J. E.; Stone, A. J. AB-initio prediction of properties of carbon dioxide, ammonia, and carbon dioxide...ammonia. J. Phys. Chem. 1985, 89, 2186−2192. (45) Jensen, J. H.; Gordon, M. S. An approximate formula for the intermolecular Pauli repulsion between closed shell molecules. Mol. Phys. 1996, 89, 1313−1325. (46) Jensen, J. H.; Gordon, M. S. An approximate formula for the intermolecular Pauli repulsion between closed shell molecules. II. Application to the effective fragment potential method. J. Chem. Phys. 1998, 108, 4772. (47) Jensen, J. H. Modeling intermolecular exchange integrals between nonorthogonal molecular orbitals. J. Chem. Phys. 1996, 104, 7795. (48) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. Efficient diffuse function-augmented basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li−F. J. Comput. Chem. 1983, 4, 294−301. (49) Hariharan, P. C.; Pople, J. A. The influence of polarization functions on molecular orbital hydrogenation energies. Theor. Chim. Acta 1973, 28, 213−222. (50) Hehre, W. J. SelfConsistent Molecular Orbital Methods. XII. Further Extensions of GaussianType Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257. (51) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-consistent molecular orbital methods 25. Supplementary functions for Gaussian basis sets. J. Chem. Phys. 1984, 80, 3265. (52) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Selfconsistent molecular orbital methods. XX. A basis set for correlated wave functions. J. Chem. Phys. 1980, 72, 650. (53) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865. (54) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656−1676. (55) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (56) Schmidt, M. W.; et al. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347−1363. (57) Gordon, M. S.; Schmidt, M. W. Advances in Electronic Structure Theory: GAMESS a Decade Later. In Theory and Applications of Computational Chemistry: The First Forty Years; Elsevier: Amsterdam, 2005. (58) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684. (59) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255. (60) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695. (61) Klopper, W.; Schütz, M.; Lüthi, H. P.; Leutwyler, S. An ab initio derived torsional potential energy surface for (H2O)3. II. Benchmark studies and interaction energies. J. Chem. Phys. 1995, 103, 1085. (62) Novoa, J. J.; Sosa, C. Evaluation of the Density Functional Approximation on the Computation of Hydrogen Bond Interactions. J. Phys. Chem. 1995, 99, 15837−15845. (63) Szczȩsń iak, M. M.; Chałasiński, G.; Cybulski, S. M.; Cieplak, P. Ab initio study of the potential energy surface of CH4−H2O. J. Chem. Phys. 1993, 98, 3078. (64) Sponer, J.; Hobza, P. MP2 and CCSD(T) study on hydrogen bonding, aromatic stacking and nonaromatic stacking. Chem. Phys. Lett. 1997, 267, 263−270. (65) van Duijneveldt-van de Rijdt, J. G. C. M.; van Duijneveldt, F. B. Convergence to the basis-set limit in ab initio calculations at the correlated level on the water dimer. J. Chem. Phys. 1992, 97, 5019. (66) Jurečka, P.; Šponer, J.; Č erný, J.; Hobza, P. Benchmark database of accurate (MP2 and CCSD(T) complete basis set limit) interaction energies of small model complexes, DNA base pairs, and amino acid pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985. (67) Stevens, W. J.; Fink, W. H. Frozen fragment reduced variational space analysis of hydrogen bonding interactions. Application to the water dimer. Chem. Phys. Lett. 1987, 139, 15−22. (68) Flick, J. C.; Kosenkov, D.; Hohenstein, E. G.; Sherrill, C. D.; Slipchenko, L. V. Accurate Prediction of Non-covalent Interaction Energies with the Effective Fragment Potential Method: Comparison of Energy Components to Symmetry-Adapted Perturbation Theory for the S22 Test Set. Unpublished. (69) Piquemal, J.-P.; Chevreau, H.; Gresh, N. Toward a Separate Reproduction of the Contributions to the Hartree−Fock and DFT Intermolecular Interaction Energies by Polarizable Molecular Mechanics with the SIBFA Potential†. J. Chem. Theory Comput. 2007, 3, 824−837. (70) Gresh, N. Model, Multiply Hydrogen-Bonded Water Oligomers (N = 3−20). How Closely Can a Separable, ab Initio-Grounded Molecular Mechanics Procedure Reproduce the Results of Supermolecule Quantum Chemical Computations? J. Phys. Chem. A 1997, 101, 8680−8694. (71) Piquemal, J.-P.; Chelli, R.; Procacci, P.; Gresh, N. Key Role of the Polarization Anisotropy of Water in Modeling Classical Polarizable Force Fields. J. Phys. Chem. A 2007, 111, 8170−8176. (72) York, D. M.; Yang, W. A chemical potential equalization method for molecular simulations. J. Chem. Phys. 1996, 104, 159−172. (73) Sorenson, J. M.; Hura, G.; Glaeser, R. M.; Head-Gordon, T. What can x-ray scattering tell us about the radial distribution functions of water? J. Chem. Phys. 2000, 113, 9149−9161.