The Melting Temperature of Liquid Water with the Effective Fragment

Aug 26, 2015 - The direct simulation of the solid–liquid water interface with the effective fragment potential (EFP) via the constant enthalpy and p...
0 downloads 6 Views 1MB Size
Download PDF
Letter pubs.acs.org/JPCL

The Melting Temperature of Liquid Water with the Effective Fragment Potential Kurt R. Brorsen,† Soohaeng Yoo Willow,‡ Sotiris S. Xantheas,§ and Mark S. Gordon*,† †

Ames Laboratory, US-DOE and Department of Chemistry, Iowa State University, Ames, Iowa 50011, United States Department of Chemistry, University of Illinois at Urbana−Champaign, 600 South Mathews Avenue, Urbana, Illinois 61801, United States § Physical Sciences Division, Pacific Northwest National Laboratory, 902 Battelle Boulevard, P.O. Box 999, MS K1-83, Richland, Washington 99352, United States ‡

ABSTRACT: The direct simulation of the solid−liquid water interface with the effective fragment potential (EFP) via the constant enthalpy and pressure (NPH) ensemble was used to estimate the melting temperature (Tm) of ice-Ih. Initial configurations and velocities, taken from equilibrated constant pressure and temperature (NPT) simulations at P = 1 atm and T = 305 K, 325 K and 399 K, respectively, yielded corresponding Tm values of 378 ± 16 K, 382 ± 14 K and 384 ± 15 K. These estimates are consistently higher than experiment, albeit to the same degree as previously reported estimates using density functional theory (DFT)-based Born−Oppenheimer simulations with the Becke-Lee−Yang−Parr functional plus dispersion corrections (BLYP-D).

C

liquid water at T = 350 K.16 Based on the evidence presented by those earlier results, previous studies2,3 have suggested that DFT-based molecular dynamics (MD) simulations performed at the standard ambient temperature can be interpreted as describing supercooled rather than the liquid state of water. The notion of a supercooled state at 298 K is consistent with the existence of an overstructured RDF for water produced by DFT MD simulations at the standard ambient temperature.11 For this reason, recent ambient temperature DFT-based simulations have been performed at elevated temperatures (330 K and above) to ensure that a liquid state is sampled.17,18 As noted earlier, estimates of Tm have been reported using classical water interaction potentials.1,19−26 A simple point charge model, such as TIP3P, predicts a melting point of 146 K.27 TIP4P/2005, improves on TIP3P, but still predicts a too low Tm of 251 K.20 The too low predicted melting point was one of the reasons that the TIP4P model has been reparametrized to more accurately predict the phase diagram of ice, resulting in TIP4P/Ice. The TIP4P/Ice model has been designed to reproduce Tm by fitting the model potential parameters and predicts a value of 270 K.20 The TIP5P potential predicts an accurate value of 272 K,26 but also predicts that ice-Ih is not thermodynamically stable at 1 bar.28 With recent improvements in force field method developments,25,29−44 classical force fields can be generated from ab initio calculations as an alternative to fitting the parameters to empirical data. Examples of force fields generated by fitting parameters to ab initio calculations include iAMOEBA,25 TTM3-F,29 WHBB,41,45 and MB-Pol.43,44 The effective frag-

hemistry at standard ambient temperature and pressure is ubiquitous. For this reason, the accurate description of the phase diagram of water is an important benchmark of any method, whether it is based on electronic structure or a classical description of the interaction between the water molecules. Given the fact that the melting temperature (Tm) of ice is just 25 K below the ambient temperature, even a small error in the description of the phase diagram of water with a particular model can result in sampling the wrong phase when performing ambient temperature simulations. Estimates of Tm have previously been reported using several classical interaction potentials (cf. ref 1 for a recent review) as well as density functional theory (DFT).2−4 Whereas the reported predictions with classical potentials are generally lower than experiment (e.g., 146 K with the TIP3P potential, an exception being the TIP4P-FQ potential which predicts 303 K), the reported predictions with DFT were notably much higher than experiment.2 Previous DFT-based simulations with the PBE and BLYP functionals5−7 predicted Tm to be 417 and 411 K, respectively.2 Dispersion corrections,8−10 when added to the BLYP functional, lowers the predicted Tm to 360 K, still too high by nearly 90 K with respect to experiment.3 The DFT and DFT-D predictions for the melting temperature broadly agree with previous studies, which suggested that DFT yields an overstructured liquid at ρ = 1 g/cm 3 and ambient conditions,11−13 that a temperature of ∼415 K was necessary for the PBE functional to obtain a value for the diffusion coefficient (D) comparable to experiment at ambient conditions (based solely on the comparison of the value of D and the corresponding radial distribution function (RDF)),14 that a distinctive transition to liquid-like diffusion occurs for the PBE functional only at an elevated temperature of 400 K,15 and that the BLYP-D functional provides a consistent description of © 2015 American Chemical Society

Received: August 5, 2015 Accepted: August 26, 2015 Published: August 26, 2015 3555

DOI: 10.1021/acs.jpclett.5b01702 J. Phys. Chem. Lett. 2015, 6, 3555−3559

Letter

The Journal of Physical Chemistry Letters

chemical potential of the ice phase will be higher than the one for the liquid phase and water molecules from the ice phase will absorb heat and melt (become liquid-like) from the ice surface. The coexisting solid−liquid phase approach has been previously used in conjunction with DFT simulations to estimate Tm for metals such as Al,52 Ta,53 Si,54 MgO,55 and water under high pressure.4 The EFP method was used to describe the waters in all MD simulations, using the electronic structure program GAMESS.56 Since no information exists about the EFP phase diagram for water, the EFP liquid-ice coexistence boundary must first be identified. Note that there is no liquid-ice coexistence boundary below the pressure at the triple point.57 To confirm the use of a pressure of 1.0 atm in the EFP simulations described below, a NVT MD simulation was performed with 216 water molecules at T = 300 K and a density of 1 g/cm3, to make sure that the EFP pressure of water is not close to the triple point. This MD simulation was performed with a step size of 0.25 fs, the velocity−Verlet algorithm, and a cubic simulation cell with edge length of 18.62 Å. After 2 ps of equilibration, the next 2 ps was used to calculate the average pressure of the system to be −200 atm. Hence, the use of a pressure of 1.0 atm is reasonable for the calculation of the melting temperature for the molecular dynamic simulations in the (NPH) ensembles as the liquid-ice coexistence should be stable at this pressure. For comparison, analogous calculations using DFT with the BLYP (PBE) functional finds 10000 atm (2500 atm).2 To calculate Tm, an ice-liquid system of 192 water molecules was prepared in a simulation cell of 27.32 × 15.61 × 14.72 Å. A representative geometry of the ice-liquid system is shown in Figure 1. For the initial configuration of the ice-liquid system, 96 of the water molecules were assumed to be ice-Ih-like (top half of Figure 1) and were arranged according to the Bernal− Fowler rules such that the ice phase has a net zero dipole moment.58 The remaining 96 waters (bottom half of Figure 1) were assumed to be liquid-like and were prepared from a

ment potential (EFP) method is also obtained directly from ab initio quantum chemistry, but the EFP method contains no empirically fitted parameters.30,31 The iAMOEBA and TTM3-F potentials predicted Tm values of 261 K25 and 248 K,3 respectively. The present study will focus on the prediction of Tm with the EFP potential. In the EFP method, each water molecule is represented as a fragment with fixed internal geometry. The EFP method contains five interaction energy terms that represent the fundamental types of intermolecular interactions: E EFP = ECoul + Epol + Eexrep + Edisp + ECT

(1)

In eq 1, ECoul is the Coulombic interaction, expressed as a distributed multipole expansion. Epol is the polarization interaction, which is expressed in terms of localized molecular orbital (LMO) polarizability tensors and is iterated to selfconsistency. The polarizability term accounts for the manybody interactions via this iterative process. Eexrep represents the exchange repulsion interaction and is expressed in a power series of the intermolecular overlap. Because of the use of localized orbitals, the series is successfully truncated at the quadratic term. Edisp is the dispersion interaction, expressed in terms of LMO frequency-dependent polarizabilities, integrated over the imaginary frequency range. ECT is the charge transfer interaction, obtained in terms of the interaction between the occupied orbitals on one fragment and the virtual orbitals on the second fragment. Details regarding these terms can be obtained from the original references. Using coupled cluster CCSD(T) calculations as a reference, it has been demonstrated that the EFP predictions of interaction energies for a wide range of interaction types, including water clusters, are comparable in accuracy to that of second order perturbation theory (MP2), while requiring orders of magnitude less computational resources.46 The melting temperature can be calculated in at least two ways: by a Gibbs−Duhem integration of the free energy47 or by direct simulation of the solid−liquid interface.48−50 At the melting point, G liq (P , T )T = Tm = Gsolid(P , T )T = Tm

(2)

so the determination of Tm by Gibbs−Duhem integration requires the nontrivial determination of the Gibbs free energy for both the ice and liquid states. Additional information about the relative benefits and drawbacks of Gibbs−Duhem integration and direct coexistence simulations in different ensembles can be found in a recent review.51 An alternative approach, adopted in the present study, is to calculate Tm by direct coexistence simulations of the ice−water interface. The direct coexistence simulations can be performed using different ensembles. The present study uses the constant pressure and enthalpy (NPH) ensemble, since this allows the temperature of the system to adjust spontaneously to Tm, where eq 2 is satisfied. Unlike the (NVE) ensemble, which also allows the temperature to spontaneously adjust, the (NPH) ensemble also allows the pressure to be maintained during the simulation. During the simulation of the ice−water coexistence system with the (NPH) ensemble, if the initial temperature of the system is lower than the melting point, T < Tm, then the chemical potential of the liquid phase will be higher than the one of ice and water molecules from the liquid phase will release heat by freezing onto the ice surface. Likewise, if the initial temperature of the system is higher than the melting point, T > Tm, then the

Figure 1. A representative configuration of the ice−water system with 96 ice-Ih-like and 96 liquid-like water molecules. 3556

DOI: 10.1021/acs.jpclett.5b01702 J. Phys. Chem. Lett. 2015, 6, 3555−3559

Letter

The Journal of Physical Chemistry Letters (NVT) equilibration at T = 300 K for 100 ps with a step size of 0.25 fs using the velocity−Verlet algorithm in a simulation cell of 13.52 × 15.61 × 14.72 Å. The ice-liquid system in the present study is the same size and prepared in a similar manner as the one used in previous DFT calculations2,3 of Tm to allow the direct comparison of Tm between the EFP method and the previous DFT simulations. The predicted Tm is apparently not converged with respect to system size for a simulation cell of 192 waters. Previous studies27 with the TIP4P potential have suggested that the calculated Tm via the (NVE) ensemble increases with system size, rising from ∼200 K with a simulation cell of 192 molecules to 229 ± 1 K with a simulation cell of 12,288 water molecules. The previous DFT studies,2,3 upon which the present work is based, estimated Tm using a simulation cell of 192 water molecules. The authors of these previous studies assumed that the change of Tm as a function of simulation cell size was qualitatively similar for TIP4P and DFT, and therefore their calculated Tm with DFT was most likely a lower bound. This assumption is made in the present work as well. Additionally, nuclear quantum effects have been shown to be important for a correct theoretical description of water at this temperature range.58 The inclusion of nuclear quantum effects in previous studies59,60 suggested a downward shift of the temperature range by about 30 K. As in the previous DFT studies2,3 on which this work is based, nuclear quantum effects are not considered in this study. Initial MD simulations for 500 fs with a 0.25 fs step size were carried out in the (NPT) ensemble at P = 1.0 atm with the velocity−Verlet algorithm at temperatures between 300 and 400 K on the merged ice-liquid systems to allow the ice-liquid system to relax. The (NPT) simulations at 300, 350, and 400 K produced configurations at 305, 325, and 399 K, which were the initial starting conditions for each of the (NPH) simulations during which the Tm was estimated. The (NPH) simulations used a 0.25 fs step size, a pressure of 1.0 atm, and were performed using the velocity−Verlet algorithm. For the (NPH) simulations with initial temperatures of 325 and 399 K, the simulations were assumed to converge after 10 ps, and the next 10 ps were used to determine the value and corresponding error bars. For the (NPH) simulation with an initial temperature of 305 K, the simulation was assumed to converge after 25 ps and the next 10 ps were used to determine the value and corresponding error bars. The Tm value was obtained as the average of the individual steps of the last 10 ps of T values after the (NPH) simulation had been assumed to converge and the error bars were calculated as the standard deviation of those individual steps in the last 10 ps. Periodic boundary conditions were imposed with the minimum image convention during all MD simulations. For the periodic boundary conditions, a switching function that gradually drops the fragment−fragment interactions from full strength to a value SWR1 to zero at a value SWR2 was used. For all MD simulations, SWR2 was half the smallest box dimension, while SWR1 was 0.8 × SWR2. The temperature versus the MD simulation time for the three (NPH) simulations starting at T = 305 K (green), 325 K (blue) and 399 K (red) is shown in Figure 2. The three simulations (starting from the above three different initial conditions) converge to the values of 378 ± 16 K, 382 ± 14 K, and 384 ± 15 K, respectively. The Tm estimates from the three simulations lie within the error bars of one another.

Figure 2. Temperature of the liquid−ice (NPH) simulations starting from initial temperatures of T = 305 K (green), 325 K (blue), and 399 K (red). The horizontal bars at t = 0 indicate the starting temperatures of the three simulations.

A summary of the Tm’s predicted by the EFP method compared to other predictions reported with DFT and classical potentials is presented in Table 1. The estimate of Tm with the Table 1. Melting Temperature of Ice-Ih with Various Methods for Describing the Water−Water Interactions potential

melting temperature (K)

EFP BLYP2 PBE2 BLYP-D3 TIP3P19 TIP5P20 TTM3-F3 iAMOEBA25 TIP4P/ICE20 TIP4P/200520 POL31

381 417 411 360 146 272 248 261 270 251 180

EFP method is substantially higher than experiment, smaller than the predictions with the BLYP and PBE functionals (>400 K) and similar to the one predicted by BLYP-D. The agreement between the EFP and DFT methods (especially BLYP-D) is not entirely surprising, since the EFP method is also derived from first-principles and includes a dispersion term. In contrast, most classical potentials seem to produce values that are smaller than experiment. Based on the current results and the ones with previous DFT-based simulations, it appears that electronic state based methods systematically overestimate Tm by ∼100 K, though it should be noted that the previous DFT studies have been performed with generalized gradient approximation functionals and not with more modern meta-hybrid functionals. Additionally, the inclusion of nuclear quantum effects is expected to decrease this difference by ∼30 K or more. However, this decrease maybe offset in part by the expected increase in the melting temperature due to the size effect of the simulation box. The fact that the EFP method currently relies on the use of rigid monomers should not be a significant source of error in the estimation of Tm (even though the liquid water and ice internal geometries are different).61 Issues related to the rigid monomer geometry and its influence on the estimate of Tm will be addressed in future studies using the effective fragment molecular orbital (EFMO) method,62,63 which is based on flexible fragments. It is possible that the EFP method is 3557

DOI: 10.1021/acs.jpclett.5b01702 J. Phys. Chem. Lett. 2015, 6, 3555−3559

Letter

The Journal of Physical Chemistry Letters

(13) Kuo, I. F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I.; VandeVondele, J.; Sprik, M.; Hutter, J.; Chen, B.; Klein, M. L.; Mohamed, F.; Krack, M.; Parrinello, M. Liquid Water from First Principles: Investigation of Different Sampling Approaches. J. Phys. Chem. B 2004, 108, 12990−12998. (14) Schwegler, E.; Grossman, J. C.; Gygi, F.; Galli, G. Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water. II. J. Chem. Phys. 2004, 121, 5400. (15) Sit, P. H.-L.; Marzari, N. Static and Dynamical Properties of Heavy Water at Ambient Conditions from First-Principles Molecular Dynamics. J. Chem. Phys. 2005, 122, 204510. (16) Weber, V.; Asthagiri, D. Communication: Thermodynamics of Water Modeled using Ab Initio Simulations. J. Chem. Phys. 2010, 133, 141101. (17) Zhang, C.; Wu, J.; Galli, G.; Gygi, F. Structural and Vibrational Properties of Liquid Water from van der Waals Density Functionals. J. Chem. Theory Comput. 2011, 7, 3054−3061. (18) Baer, M. D.; Mundy, C. J.; McGrath, M. J.; Kuo, I.-F. W.; Siepmann, J. I.; Tobias, D. J. Re-examining the Properties of the Aqueous Vapor−Liquid Interface using Dispersion Corrected Density Functional Theory. J. Chem. Phys. 2011, 135, 124712. (19) Vega, C.; Sanz, E.; Abascal, J. The Melting Temperature of the Most Common Models of Water. J. Chem. Phys. 2005, 122, 114507. (20) García Fernández, R.; Abascal, J. L. F.; Vega, C. The Melting Point of Ice Ih for Common Water Models Calculated from Direct Coexistence of the Solid-Liquid Interface. J. Chem. Phys. 2006, 124, 144506. (21) Vega, C.; Abascal, J. L.; Conde, M.; Aragones, J. What Ice Can Teach Us About Water Interactions: A Critical Comparison of the Performance of Different Water Models. Faraday Discuss. 2009, 141, 251−276. (22) Paesani, F.; Xantheas, S. S.; Voth, G. A. Infrared spectroscopy and Hydrogen-Bond Dynamics of Liquid Water from Centroid Molecular Dynamics with an Ab Initio-Based Force Field. J. Phys. Chem. B 2009, 113, 13118−13130. (23) Nilsson, A.; Pettersson, L. G. Perspective on the Structure of Liquid Water. Chem. Phys. 2011, 389, 1−34. (24) Vega, C.; Abascal, J. L. Simulating Water with Rigid Nonpolarizable Models: A General Perspective. Phys. Chem. Chem. Phys. 2011, 13, 19663−19688. (25) Wang, L.-P.; Head-Gordon, T.; Ponder, J. W.; Ren, P.; Chodera, J. D.; Eastman, P. K.; Martinez, T. J.; Pande, V. S. Systematic Improvement of a Classical Molecular Model of Water. J. Phys. Chem. B 2013, 117, 9956−9972. (26) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910. (27) Wang, J.; Yoo, S.; Bai, J.; Morris, J. R.; Zeng, X. C. Melting Temperature of Ice Ih Calculated from Coexisting Solid-Liquid Phases. J. Chem. Phys. 2005, 123, 036101. (28) Sanz, E.; Vega, C.; Abascal, J.; MacDowell, L. Tracing the Phase Diagram of the Four-Site Water Potential (TIP4P). J. Chem. Phys. 2004, 121, 1165. (29) Fanourgakis, G. S.; Xantheas, S. S. Development of Transferable Interaction Potentials for Water. V. Extension of the Flexible, Polarizable, Thole-Type Model Potential (TTM3-F, v. 3.0) to Describe the Vibrational Spectra of Water Clusters and Liquid Water. J. Chem. Phys. 2008, 128, 074506. (30) Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. The Effective Fragment Potential Method: a QM-Based MM Approach to Modeling Environmental Effects in Chemistry. J. Phys. Chem. A 2001, 105, 293−307. (31) Gordon, M. S.; Slipchenko, L.; Li, H.; Jensen, J. H. The Effective Fragment Potential: a General Method for Predicting Intermolecular Interactions. Annu. Rep. Comput. Chem. 2007, 3, 177−193. (32) Gordon, M. S.; Mullin, J. M.; Pruitt, S. R.; Roskop, L. B.; Slipchenko, L. V.; Boatz, J. A. Accurate Methods for Large Molecular Systems. J. Phys. Chem. B 2009, 113, 9646−9663.

mimicking what the MP2 melting point of water would be if one performed MP2 MD simulations (excluding nuclear quantum effects) with internally frozen geometries. This possibility will be analyzed once MP2 MD simulations become feasible either by conventional MP2 or via the use of fragmentation methods33 such as the fragment molecular orbital method.64

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS K.R.B. was supported by a Computational Science Graduate Fellowship from the Department of Energy. M.S.G. was supported by a U.S. National Science Foundation Software Infrastructure (SI2) grant (ACI − 1047772). S.S.X. acknowledges support from the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences. Pacific Northwest National Laboratory (PNNL) is a multiprogram national laboratory operated for DOE by Battelle.

■

REFERENCES

(1) Muchová, E.; Gladich, I.; Picaud, S.; Hoang, P. N.; Roeselová, M. The Ice− Vapor Interface and the Melting Point of Ice I h for the Polarizable POL3 Water Model. J. Phys. Chem. A 2011, 115, 5973− 5982. (2) Yoo, S.; Zeng, X. C.; Xantheas, S. S. On the Phase Diagram of Water with Density Functional Theory Potentials: The Melting Temperature of Ice Ih with the Perdew−Burke−Ernzerhof and Becke−Lee−Yang−Parr Functionals. J. Chem. Phys. 2009, 130, 221102. (3) Yoo, S.; Xantheas, S. S. Communication: The Effect of Dispersion Corrections on the Melting Temperature of Liquid Water. J. Chem. Phys. 2011, 134 (12), 121105. (4) Schwegler, E.; Sharma, M.; Gygi, F.; Galli, G. Melting of Ice under Pressure. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 14779−14783. (5) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (6) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (7) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (8) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (9) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (10) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (11) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water. J. Chem. Phys. 2004, 120, 300. (12) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. The Influence of Temperature and Density Functional Models in Ab Initio Molecular Dynamics Simulation of Liquid Water. J. Chem. Phys. 2005, 122, 014515. 3558

DOI: 10.1021/acs.jpclett.5b01702 J. Phys. Chem. Lett. 2015, 6, 3555−3559

Letter

The Journal of Physical Chemistry Letters (33) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: a Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632−672. (34) Pruitt, S. R.; Bertoni, C.; Brorsen, K. R.; Gordon, M. S. Efficient and Accurate Fragmentation Methods. Acc. Chem. Res. 2014, 47, 2786−2794. (35) Fanourgakis, G. S.; Xantheas, S. S. The Flexible, Polarizable, Thole-Type Interaction Potential for Water (TTM2-F) Revisited. J. Phys. Chem. A 2006, 110, 4100−4106. (36) Burnham, C. J.; Xantheas, S. S. Development of Transferable Interaction Models for Water. IV. A Flexible, All-Atom Polarizable Potential (TTM2-F) Based on Geometry Dependent Charges Derived from an Ab Initio Monomer Dipole Moment Surface. J. Chem. Phys. 2002, 116, 5115. (37) Liu, H.; Wang, Y.; Bowman, J. M. Quantum Calculations of Intramolecular IR Spectra of Ice Models Using Ab Initio Potential and Dipole Moment Surfaces. J. Phys. Chem. Lett. 2012, 3, 3671−3676. (38) Liu, H.; Wang, Y.; Bowman, J. M. Quantum Calculations of the IR Spectrum of Liquid Water using Ab Initio and Model Potential and Dipole Moment Surfaces and Comparison with Experiment. J. Chem. Phys. 2015, 142, 194502. (39) Wang, Y.; Huang, X.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Flexible, Ab Initio Potential, and Dipole Moment Surfaces for Water. I. Tests and Applications for Clusters up to the 22-mer. J. Chem. Phys. 2011, 134, 094509. (40) Wang, Y.; Bowman, J. M. Towards an Ab Initio Flexible Potential for Water, and Post-Harmonic Quantum Vibrational Analysis of Water Clusters. Chem. Phys. Lett. 2010, 491, 1−10. (41) Wang, Y.; Bowman, J. M. Ab Initio Potential and Dipole Moment Surfaces for Water. II. Local-Monomer Calculations of the Infrared Spectra of Water Clusters. J. Chem. Phys. 2011, 134, 154510. (42) Wang, Y.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. FullDimensional, Ab Initio Potential Energy and Dipole Moment Surfaces for Water. J. Chem. Phys. 2009, 131, 054511. (43) Babin, V.; Leforestier, C.; Paesani, F. Development of a “First Principles” Water Potential with Flexible Monomers: Dimer Potential Energy Surface, VRT Spectrum, and Second Virial Coefficient. J. Chem. Theory Comput. 2013, 9, 5395−5403. (44) Babin, V.; Medders, G. R.; Paesani, F. Development of a “First Principles” Water Potential with Flexible Monomers. II: Trimer Potential Energy Surface, Third Virial Coefficient, and Small Clusters. J. Chem. Theory Comput. 2014, 10, 1599−1607. (45) Shank, A.; Wang, Y.; Kaledin, A.; Braams, B. J.; Bowman, J. M. Accurate Ab Initio and “Hybrid” Potential Energy Surfaces, Intramolecular Vibrational Energies, and Classical IR Spectrum of the Water Dimer. J. Chem. Phys. 2009, 130, 144314. (46) Flick, J. C.; Kosenkov, D.; Hohenstein, E. G.; Sherrill, C. D.; Slipchenko, L. V. Accurate Prediction of Noncovalent Interaction Energies with the Effective Fragment Potential Method: Comparison of Energy Components to Symmetry-Adapted Perturbation Theory for the S22 Test Set. J. Chem. Theory Comput. 2012, 8, 2835−2843. (47) Kofke, D. A. Gibbs-Duhem Integration: A New Method for Direct Evaluation of Phase Coexistence by Molecular Simulation. Mol. Phys. 1993, 78, 1331−1336. (48) Ladd, A.; Woodcock, L. Triple-Point Coexistence Properties of the Lennard-Jones System. Chem. Phys. Lett. 1977, 51, 155−159. (49) Ladd, A.; Woodcock, L. Interfacial and Co-existence Properties of the Lennard-Jones System at the Triple Point. Mol. Phys. 1978, 36, 611−619. (50) Yoo, S.; Zeng, X. C.; Morris, J. R. The Melting Lines of Model Silicon Calculated from Coexisting Solid−Liquid Phases. J. Chem. Phys. 2004, 120, 1654. (51) Vega, C.; Sanz, E.; Abascal, J.; Noya, E. Determination of Phase Diagrams via Computer Simulation: Methodology and Applications to Water, Electrolytes and Proteins. J. Phys.: Condens. Matter 2008, 20, 153101−153138. (52) Alfè, D. First-Principles Simulations of Direct Coexistence of Solid and Liquid Aluminum. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68, 064423.

(53) Taioli, S.; Cazorla, C.; Gillan, M. J.; Alfè, D. Melting Curve of Tantalum from First Principles. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 214103. (54) Yoo, S.; Xantheas, S. S.; Zeng, X. C. The Melting Temperature of Bulk Silicon from Ab Initio Molecular Dynamics Simulations. Chem. Phys. Lett. 2009, 481, 88−90. (55) Alfè, D. Melting Curve of MgO from First-Principles Simulations. Phys. Rev. Lett. 2005, 94, 235701. (56) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (57) Cengel, Y. A.; Turner, R. H. Fundamentals of Thermal-Fluid Sciences; McGraw-Hill: Boston, 2004. (58) Hayward, J.; Reimers, J. Unit Cells for the Simulation of Hexagonal Ice. J. Chem. Phys. 1997, 106, 1518. (59) Kuharski, R. A.; Rossky, P. J. A Quantum Mechanical Study of Structure in Liquid H2O and D2O. J. Chem. Phys. 1985, 82, 5164. (60) Paesani, F.; Iuchi, S.; Voth, G. A. Quantum Effects in Liquid Water from an Ab Initio-Based Polarizable Force Field. J. Chem. Phys. 2007, 127, 074506. (61) Fanourgakis, G. S.; Xantheas, S. S. The Bend Angle of Water in Ice Ih and Liquid Water: The Significance of Implementing the Nonlinear Monomer Dipole Moment Surface in Classical Interaction Potentials. J. Chem. Phys. 2006, 124, 174504. (62) Pruitt, S. R.; Steinmann, C.; Jensen, J. H.; Gordon, M. S. Fully Integrated Effective Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2013, 9, 2235−2249. (63) Steinmann, C.; Fedorov, D. G.; Jensen, J. H. Effective Fragment Molecular Orbital Method: A Merger of the Effective Fragment Potential and Fragment Molecular Orbital Methods. J. Phys. Chem. A 2010, 114, 8705−8712. (64) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment Molecular Orbital Method: An Approximate Computational Method for Large Molecules. Chem. Phys. Lett. 1999, 313, 701−706.

3559

DOI: 10.1021/acs.jpclett.5b01702 J. Phys. Chem. Lett. 2015, 6, 3555−3559