This study presents the extension of the MB-nrg (Many-Body energy) theoretical/computational framework of transferable potential energy functions (PEFs) for molecular simulations of alkali metal ion-water systems. The MB-nrg PEFs are built upon the many-body expansion of the total energy and include the explicit treatment of one-body, two-body, and three-body interactions, with all higher-order contributions described by classical induction. This study focuses on the MB-nrg two-body terms describing the full-dimensional potential energy surfaces of the M+(H2O) dimers, where M+ = Li+, Na+, K+, Rb+, and Cs+. The MB-nrg PEFs are derived entirely from “first principles” calculations carried out at the explicitly correlated coupled-cluster level including single, double, and perturbative triple excitations [CCSD(T)-F12b] for Li+ and Na+ and at the CCSD(T) level for K+, Rb+, and Cs+. The accuracy of the MB-nrg PEFs is systematically assessed through an extensive analysis of interaction energies, structures, and harmonic frequencies for all five M+(H2O) dimers. In all cases, the MB-nrg PEFs are shown to be superior to both polarizable force fields and ab initio models based on density functional theory. As previously demonstrated for halide-water dimers, the MB-nrg PEFs achieve higher accuracy by correctly describing short-range quantum-mechanical effects associated with electron density overlap as well as long-range electrostatic many-body interactions.

Electrolyte solutions are ubiquitous, occurring in cells and the environment, and mediate many industrial and technological processes.1–7 In aqueous solutions, the interactions between different ions and surrounding water molecules give rise to specific ion effects, which depend on both the physicochemical properties of the individual ions and the thermodynamic conditions of the solutions.8 While there have been several experimental and theoretical studies on ion hydration, a definitive picture of how ions are hydrated and to what extent structural differences in the organization of the hydration shells translate into differences in thermodynamic and dynamical properties of electrolyte solutions is still lacking.9–19 

Significant progress in understanding ion hydration has been made through computer simulations that rely on either force fields9,20–26 or ab initio models.8,20,27 In particular, a large number of polarizable force fields have been developed to model the behavior of both halide and alkali metal ions in solution.21,24,25,28–36 Most of these models, however, treat the water molecules as rigid, which precludes direct comparisons with experimental vibrational spectra measured from the gas and the condensed phases.37–40 A remarkable exception is the AMOEBA force field that includes both intramolecular flexibility and many-body effects through classical polarization.24 The AMOEBA force field is under continued development, with recent efforts focusing on accurate parameterizations of metal ion-water interactions, including alkali metals,41 alkaline-earth metals,42 transition metals,43,44 and lanthanides and actinides.45 

Recent years have witnessed the advent of highly accurate water potential energy functions (PEFs) rigorously derived from the many-body expansion (MBE) of the interaction energies (e.g., CC-pol,46–48 WHBB,49–51 HBB2-pol,52,53 and MB-pol54–56), enabling computer simulations of water with high accuracy.57 Among existing many-body PEFs, it has been shown that MB-pol currently provides the most accurate representation of the molecular properties of water from the gas to the condensed phase,58,59 correctly predicting the vibration-rotation tunneling spectrum of the dimer,54 the relative stability and vibrational spectra of small water clusters,60,61 the structural, thermodynamic, and dynamical properties56 as well as the IR and Raman spectra of liquid water,62 the vibrational sum-frequency generation (vSFG) spectrum of the air/water interface,63 and the energetics of different ice phases.64 

Building on the MB-pol accuracy for water, a new class of many-body PEFs (MB-nrg for “many-body energy”) has recently been introduced to describe halide-water interactions, explicitly including water flexibility.65 The MB-nrg PEFs are derived entirely from correlated electronic structure data and are expressed through the MBE as66 

EN=i=1NV1B(i)+i<jNV2B(i,j)+i<j<kNV3B(i,j,k)++VNB(1,,N),
(1)

where V1B is the water monomer energy, and VnB are the n-body (nB) interactions defined recursively as

VnB(1,,n)=En(1,,n)iV1B(i)i<jV2B(i,j)i<j<<n1V(n-1)B(i,j,,(n1)).
(2)

Since Eq. (1) converges quickly for nonmetallic systems, the MBE provides a rigorous and efficient framework for the development of full-dimensional PEFs in which the low-order interaction terms can be accurately determined using correlated electronic structure methods [e.g., coupled cluster level of theory including single, double, and perturbative triple excitations, CCSD(T), the current “gold standard” for chemical accuracy] and higher-order terms can be effectively represented by many-body induction. In this study, we introduce MB-nrg PEFs for alkali metal ion-water interactions with a specific focus on the 2B M+(H2O) term, where M+ = Li+, Na+, K+, Rb+, and Cs+. The article is organized as follows: Sec. II describes the theoretical and computational details associated with the development of the MB-nrg PEFs. The accuracy of the MB-nrg 2B term is assessed in Sec. III through the analysis of the interaction energies, structures, and harmonic frequencies of M+(H2O) dimers. A summary highlighting the main findings and discussing future directions is given in Sec. IV.

Depending on the specific alkali metal ion, all 2B energies for the M+(H2O) dimers were calculated using either coupled cluster theory with single, double, and perturbative triple excitations, CCSD(T), or explicitly correlated CCSD(T)-F12b theory.67,68 All of the reference energy calculations were performed at the complete basis set (CBS) limit, which was achieved via a two-point extrapolation69,70 between the values obtained with the correlation-consistent polarized valence triple zeta (aug-cc-pVTZ for H,O and cc-pwCVTZ for the metal ions) and quadruple zeta (aug-cc-pVQZ for H,O and cc-pwCVQZ for the metal ions) basis sets.71–74 As described in Refs. 69 and 70, the Hartree–Fock energy, the correlation energy, and the triples contribution were extrapolated separately for Li+(H2O) and Na+(H2O), where CCSD(T)-F12b. For K+(H2O), Rb+(H2O), and Cs+(H2O), the Hartree–Fock energy, the correlation energy, and the triples contribution were extrapolated separately for Li+(H2O) and Na+(H2O), for which CCSD(T)(-F12b) was used. The reference frequency calculations were obtained with CCSD(T)(-F12b) in conjunction with the aug-cc-pVTZ basis set for H and O and the cc-pwCVTZ basis set for the metal ions. The reference optimized geometries and energies were obtained with CCSD(T)(-F12b), using aug-cc-pVTZ for H and O, cc-pwCVTZ for the metal ions, and a gradient convergence of 10−6 a.u., with a subsequent single point calculation to obtain the CCSD(T)/CBS interaction energy of this structure.

Effective core potentials (ECPs) were used for K+ (ECP10MDF), Rb+ (ECP28MDF), and Cs+ (ECP46MDF).75 The dipole polarizability of each ion was computed at the CCSD(T) level of theory using the cc-pwCV5Z basis set following the methodology described in Ref. 76. The polarizability values obtained for Li+, Na+, K+, Rb+, and Cs+ are 0.0285 Å3, 0.1476 Å3, 0.8184 Å3, 1.3614 Å3, and 2.3660 Å3, respectively, in agreement with other theoretical estimates previously reported in the literature.77,78 All CCSD(T), CCSD(T)-F12b, and second order Møller-Plesset perturbation theory (MP2) calculations of interaction energies and harmonic frequencies were performed with MOLPRO, version 2015.1.79 Density functional theory (DFT) calculations of interaction energies and harmonic vibrational frequencies were carried out with Gaussian 0980 and the aug-cc-pVQZ basis set for all functionals except revPBE, revPBE0, and ωB97M-V,81 for which Q-Chem 5.082 was used. The Q-Chem calculations were performed with the aug-cc-pVQZ basis set and an integration grid consisting of 99 radial shells with 590 angular points per shell (99, 590).

M+(H2O) dimer interaction energies and harmonic vibrational frequencies were computed using previously reported classical polarizable force fields for comparison. The AMOEBA calculations were performed with version 6.3 of the TINKER package,83,84 using the AMOEBA 2009 force field, which includes the 2003 parameterizations for both water85 and alkali metal ions.41 The TTM-nrg (for “Thole-type model energy”) calculations were performed using in-house software. The original TTM-nrg parameters reported in Ref. 86 were used for Li+ and Na+. However, since higher accuracy CCSD(T)/CBS reference data for K+, Rb+, and Cs+ were obtained as part of this study using the cc-pwCVTZ and cc-pwCVQZ basis sets, improved TTM-nrg parameters were derived for these ions and used in all comparisons.

As discussed in Ref. 65, all water properties (i.e., monomer distortion, dipole moment, and polarizability, as well as water-water intermolecular interactions) are described by the MB-pol potential,54–56 since it has been shown to accurately reproduce structural, thermodynamic, dynamical, and spectroscopic properties of water from the gas to the condensed phase.58,59,62 Following MB-pol, the interactions between individual ions and water are described through the MBE in Eq. (1) and include explicit 2B (ion-water) and 3B (ion-water-water) terms, with all higher-order interactions being implicitly taken into account through many-body induction. The 2B term of the MB-nrg potentials, which is the focus of the present study, includes three contributions,

V2B(H2OM+)=Vshort2B(H2OM+)+VTTM2B(H2OM+)+Vdisp2B(H2OM+).
(3)

The electrostatic term (VTTM2B) is a modified version of the extended Thole-type model (TTM) introduced in the TTM4-F water potential,87 in which the CCSD(T) dipole polarizabilities of the alkali metal ions reported in Sec. II A are used to describe the ion-water induction energy. The dispersion energy Vdisp2B is represented as

Vdisp2B(H2OM+)=f6(D6HM,RH1M)C6HM1RH1M6f6(D6HM,RH2M)C6HM1RH2M6f6(D6OM,ROM)C6OM1ROM6,
(4)

where

fn(D6R)=1exp(D6R)k=0nD6Rk!
(5)

is the Tang-Toennies damping function.88 The C6 dispersion coefficients were calculated using the exchange-hole dipole model (XDM) developed by Becke and Johnson.89 The XDM analysis was performed with the Postg software,90,91 using the wavefunction file obtained from Gaussian 0980 calculations performed with the LC-ωPBE92–94 functional and the aug-cc-pV5Z basis set for H and O and the cc-pwCV5Z basis set for all M+ ions.

In the MB-nrg energy expression, the classical electrostatic (VTTM2B) and dispersion (Vdisp2B) contributions are supplemented by a short-range term (Vshort2B) that effectively recovers quantum-mechanical effects associated with the overlap of the monomer electron densities (e.g., charge transfer, charge penetration, and Pauli repulsion). Vshort2B is described by a permutationally invariant polynomial95 that is smoothly switched to zero for O − M+ distances larger than a predefined cutoff value

Vshort2B(H2OM+)=s2ROMRiRoRiVpoly.
(6)

Here, s2(x) is a switching function defined as

s2(x)=1if x<01+x2(2x3)if 0x<10if 1x.
(7)

The inner and outer cutoff radii (Ri and Ro, respectively) of Eq. (6) were chosen based on the one-dimensional scans of each M+(H2O) dimer. Ri corresponds to the O − M+ distance at which the total and electrostatic energy differ by 0.01 kcal/mol or less, and Ro = Ri + 1 Å. Following these criteria, the cutoffs were set to 5.0 Å and 6.0 Å for Li+, 5.5 Å and 6.5 Å for Na+, and 6.0 Å and 7.0 Å for K+, Rb+, and Cs+. The permutationally invariant polynomial (Vpoly) is a function of all distances that involve physical atoms H, O, and M+, as well as the oxygen lone-pair sites (L1 and L2) defined in MB-pol.54 All distances dn=18 entering the expression of Vpoly are listed in Table I along with the corresponding variables.

TABLE I.

Distances and the corresponding variables entering the short range part of the potential, Vpoly.

d1 H1 H2 ξ1=ekHHintrad1dHHintra 
d2 H1 ξ2=ekOHintrad2dOHintra 
d3 H2 ξ3=ekOHintrad3dOHintra 
d4 H1 ξ4=ekMHcould4dMHcould4 
d5 H2 ξ5=ekMHcould5dMHcould5 
d6 ξ6=ekMOcould6dMOcould6 
d7 L1 ξ7=ekMLd7dML 
d8 L2 ξ8=ekMLd8dML 
d1 H1 H2 ξ1=ekHHintrad1dHHintra 
d2 H1 ξ2=ekOHintrad2dOHintra 
d3 H2 ξ3=ekOHintrad3dOHintra 
d4 H1 ξ4=ekMHcould4dMHcould4 
d5 H2 ξ5=ekMHcould5dMHcould5 
d6 ξ6=ekMOcould6dMOcould6 
d7 L1 ξ7=ekMLd7dML 
d8 L2 ξ8=ekMLd8dML 

Vpoly is a function of two different types of variables, intramolecular (ξ1,ξ2,ξ3) and intermolecular (ξ4,ξ5,ξ6,ξ7,ξ8) variables, with the permutational invariance imposed with respect to permutations of equivalent sites within the same water molecule (i.e., H1, H2). A total of 429 symmetrized monomials form Vpoly: 3 first-degree monomials formed from all intermolecular variables (ξ4,ξ5,ξ6,ξ7,ξ8), 15 symmetrized second-degree monomials with at most a linear dependence on intramolecular variables, 49 symmetrized third-degree monomials, 119 symmetrized fourth-degree monomials, and 243 symmetrized fifth-degree monomials with at most a quadratic dependence on the intramolecular variables. Vpoly thus contains 429 linear fitting parameters (ci), and 10 nonlinear fitting parameters, kHHintra, kOHintra, kMHcoul, kMOcoul, kML, dHHintra, dOHintra, dMHcoul, dMOcoul, and dML. The nonlinear parameters γ and γ were set to the same values as in MB-pol, while the dispersion damping coefficients D6OX and D6HX were set at the same values as in the TTM-nrg models.86 This ensures that the underlying classical representation of alkali metal ion-water interactions is the same in both TTM-nrg and MB-nrg PEFs, with the only differences between the two models being in the description of short-range interactions.

The final training sets for the MB-nrg PEFs were generated from a uniform spherical grid of ion positions around a water molecule, including ion-water separations between 1.6 Å and 8 Å. For each position of the ion, different configurations of the water molecule were included. The distortions were obtained by uniformly elongating or shortening one or both O–H bonds, together with a change in the angle, including distortions up to 60 kcal/mol. The original grid, containing approximately 2.7 × 106 configurations, was screened using the RMSD (Root Mean Square Deviation) of the distances between all atoms and the center of mass of the M+(H2O) dimer as similarity criteria. If the RMSD between two dimers was smaller than a threshold value (TV), the two configurations were considered equivalent and the second one was removed. In the screening process, the specific value of TV varied as a function of the distance between the ion and the oxygen atom of the water molecule, according to

TV=1.01.0+e0.6(krd)+0.005,
(8)

where

d=min(dO1M1,dO2M2)
(9)

and

kr=k1ifmin(IE1,IE2)Emin2k2ifmin(IE1,IE2)Emin4k3ifmin(IE1,IE2)0.0k4ifmin(IE1,IE2)Emin+ΔEk5ifmin(IE1,IE2)>Emin+ΔE.
(10)

Here, IEi is the 2B energy of the ith M+(H2O) dimer calculated with the corresponding TTM-nrg PEF,86 and Emin is the 2B energy of the minimum energy configuration M+(H2O) dimer. After careful investigation, k1, k2, k3, k4, and k5 were set to 10, 9, 8, 7, and 6 for Li+(H2O); 11.5, 10.5, 9.5, 8, and 7 for Na+(H2O); 12, 11, 10, 9, and 8 for K+(H2O); 13, 12, 11, 10, and 8 for Rb+(H2O); and 13.5, 12.5, 11.5, 10.5, and 9 for Cs+(H2O), respectively, to obtain training sets of roughly 15 000 configurations for each alkali metal ion-water dimer. ΔE was set to the same value as in the TTM-nrg potentials.86 These initial grids were further refined to give a more uniform distribution of dimer configurations at both short and long range.

Adopting the same fitting protocol established for MB-pol54,55 and the halide-water MB-nrg PEFs,65 the linear and nonlinear parameters of the MB-nrg PEFs were optimized with singular value decomposition and the simplex algorithm, respectively, through the minimization of the regularized weighted sum of squared deviations calculated for the corresponding training set (S), commonly known as Tikhonov regularization or ridge regression,96 

χ2=nSwn[V2model(n)V2ref(n)]2+Γ2i=1429ci2,
(11)

where the weights, wn, were set to emphasize dimers with lower total energy according to

w(Ei)=ΔEEiEmin+ΔE2.
(12)

Here, ΔE is a parameter used to favorably weight low-energy configurations and was chosen such that configurations with Ei>20 kcal/mol have weights w(Ei)0.25. Specifically, ΔE was set to 50 kcal/mol for Li+(H2O), 40 kcal/mol for Na+(H2O), K+(H2O), and Rb+(H2O), and 30 kcal/mol for Cs+(H2O). The regularization parameter, Γ, was set to 5×104 in order to reduce the variation of the linear parameters without spoiling the overall accuracy of the fits. Figure 1 shows the correlation plots for the test sets, which consist of roughly 400 dimer configurations with 2B energies between −40 and 50 kcal/mol. For all five dimers, the MB-nrg PEFs accurately reproduce the CCSD(T) reference energies. The RMSDs of the test sets are 0.23 kcal/mol, 0.09 kcal/mol, 0.14 kcal/mol, 0.05 kcal/mol, and 0.05 kcal/mol for Li+(H2O), Na+(H2O), K+(H2O), Rb+(H2O), and Cs+(H2O), respectively.

FIG. 1.

Correlation plots of the test set for the 2B energy for M+(H2O) (M+ = Li+, Na+, K+, Rb+, Cs+) dimers. The CCSD(T)-F12b/CBS 2B energies are plotted on the x-axis, and the corresponding 2B energies calculated with the MB-nrg potentials are plotted on the y-axis. Each one of the five test sets contains roughly 400 configurations.

FIG. 1.

Correlation plots of the test set for the 2B energy for M+(H2O) (M+ = Li+, Na+, K+, Rb+, Cs+) dimers. The CCSD(T)-F12b/CBS 2B energies are plotted on the x-axis, and the corresponding 2B energies calculated with the MB-nrg potentials are plotted on the y-axis. Each one of the five test sets contains roughly 400 configurations.

Close modal

The MB-nrg minimum energy structures of each M+(H2O) dimer are shown in Fig. 2, and the associated 2B energies (i.e., interaction energies) are compared in Table II with the corresponding values calculated using both polarizable force fields and correlated electronic structure methods. While the MB-nrg PEFs quantitatively reproduce the CCSD(T) reference data, providing higher accuracy than MP2, both TTM-nrg and AMOEBA significantly underestimate (up to 1.5 kcal/mol) the interaction energies of all dimers.

FIG. 2.

Minimum energy geometries for the M+(H2O) dimers optimized with the MB-nrg potentials. (a) Li+(H2O), (b) Na+(H2O), (c) K+(H2O), (d) Rb+(H2O), and (e) Cs+(H2O).

FIG. 2.

Minimum energy geometries for the M+(H2O) dimers optimized with the MB-nrg potentials. (a) Li+(H2O), (b) Na+(H2O), (c) K+(H2O), (d) Rb+(H2O), and (e) Cs+(H2O).

Close modal
TABLE II.

Comparison between the 2B energies (in kcal/mol) for the optimized structures at the CCSD(T)(-F12b) level of theory for all M+(H2O) dimers, calculated using the MB-nrg potentials and the corresponding values obtained at the CCSD(T)(-F12b) and MP2 levels of theory, as well as using the AMOEBA and TTM-nrg force fields.

Method/modelLi+(H2O)Na+(H2O)K+(H2O)Rb+(H2O)Cs+(H2O)
AMOEBA24  −33.32 −23.43 −17.65 −13.95 −12.64 
TTM-nrg86  −35.29 −24.91 −17.01 −14.84 −12.75 
MP2 −34.51 −23.99 −18.09 −16.31 −14.51 
CCSD(T)-F12/CBS −34.80 −24.18 −18.03 −16.17 −14.43 
MB-nrg −34.79 −24.20 −18.03 −16.18 −14.42 
Method/modelLi+(H2O)Na+(H2O)K+(H2O)Rb+(H2O)Cs+(H2O)
AMOEBA24  −33.32 −23.43 −17.65 −13.95 −12.64 
TTM-nrg86  −35.29 −24.91 −17.01 −14.84 −12.75 
MP2 −34.51 −23.99 −18.09 −16.31 −14.51 
CCSD(T)-F12/CBS −34.80 −24.18 −18.03 −16.17 −14.43 
MB-nrg −34.79 −24.20 −18.03 −16.18 −14.42 

To assess the accuracy of the MB-nrg PEFs for molecular configurations far from the corresponding minimum-energy structures, the interaction energies calculated along one-dimensional radial scans for different orientations of the M+ ions relative to a water molecule lying on the xy-plane are shown in Fig. 3. In these calculations, the water molecule was kept in the vibrationally averaged geometry,97 with the oxygen atom defining the origin of the coordinate frame. The ion position is defined by the distance R from O and the polar (𝜃) and azimuthal (ϕ) angles. Although none of these configurations were included in the training sets, the MB-nrg PEFs (solid lines) correctly reproduce the CCSD(T)(-F12)/CBS values (dotted lines) at all ion-water separations and orientations for all M+(H2O) dimers.

FIG. 3.

Radial scans of the M+(H2O) (M+ = Li+, Na+, K+, Rb+, Cs+) dimer PESs. The distance between M+ and O is plotted on the x-axis, and the 2B energies are plotted on the y-axis. The symbols correspond to the CCSD(T)(-F12b)/CBS values, while MB-nrg energies are shown with solid lines. The different orientations, 𝜃, ϕ, of M+ relative to H2O are given within the parentheses, and the bottom right panel shows the reference frame.

FIG. 3.

Radial scans of the M+(H2O) (M+ = Li+, Na+, K+, Rb+, Cs+) dimer PESs. The distance between M+ and O is plotted on the x-axis, and the 2B energies are plotted on the y-axis. The symbols correspond to the CCSD(T)(-F12b)/CBS values, while MB-nrg energies are shown with solid lines. The different orientations, 𝜃, ϕ, of M+ relative to H2O are given within the parentheses, and the bottom right panel shows the reference frame.

Close modal

An additional test for the accuracy of the MB-nrg PEFs is represented by the analysis of the potential energy angular profiles (PEAPs) along ϕ shown in Fig. 4 for each M+(H2O) dimer, which directly probe the ability of different models to describe the anisotropy of the underlying potential energy surfaces. As for the radial scans, the water molecule was kept in the vibrationally averaged geometry and R was set to 1.83 Å, 2.21 Å, 2.59 Å, 2.75 Å, and 2.94 Å for M+ = Li+, Na+, K+, Rb+, and Cs+, respectively, which corresponds to the equilibrium distances within each M+(H2O) dimer for 𝜃=0°. Independent of M+, the MB-nrg PEFs quantitatively reproduce the CCSD(T)(-F12)/CBS data at all orientations. By contrast, the TTM-nrg PEFs clearly display some deficiencies in predicting the correct energies of the transition states at 𝜃60°, with the deviations from the reference data becoming smaller as the ion size increases, i.e., as the ion-water interactions become more classical-like. Similar behavior was exhibited by the AMOEBA force field in Ref. 86, which reinforces the notion that purely classical PEFs are not capable of correctly reproducing short-range interactions where quantum-mechanical effects (e.g., charge transfer, charge penetration, and Pauli repulsion) due to the overlap of the monomer densities are significant.

FIG. 4.

Potential energy angular profiles for M+(H2O) (M+ = Li+, Na+, K+, Rb+, Cs+) calculated as a function of ϕ for 𝜃 = 90°, as shown in the top right scheme, and R equal to the minimum energy distance of each dimer, with the water molecule fixed at its vibrationally averaged geometry. The symbols correspond to the CCSD(T)(-F12b) /CBS values. The MB-nrg energies are shown in the left panels as solid lines, and the corresponding TTM-nrg are shown as dashed lines in the right panels.

FIG. 4.

Potential energy angular profiles for M+(H2O) (M+ = Li+, Na+, K+, Rb+, Cs+) calculated as a function of ϕ for 𝜃 = 90°, as shown in the top right scheme, and R equal to the minimum energy distance of each dimer, with the water molecule fixed at its vibrationally averaged geometry. The symbols correspond to the CCSD(T)(-F12b) /CBS values. The MB-nrg energies are shown in the left panels as solid lines, and the corresponding TTM-nrg are shown as dashed lines in the right panels.

Close modal

A systematic comparison between the MB-nrg PEFs and different DFT models, with and without the D3 empirical dispersion correction,98 is presented in this section. The analysis includes results for the energetics and harmonic frequencies calculated with GGA (Generalized Gradient Approximation) functionals (BLYP,99,100 BLYP-D3, PBE,101 PBE-D3, revPBE,102 revPBE-D3), meta-GGA functionals (TPSS,103 TPSS-D3), hybrid GGA functionals (PBE0,104 PBE0-D3, revPBE0,102,104 revPBE0-D3, B3LYP,99,100,105 B3LYP-D3), and range-separated hybrid functionals (LC-ωPBE,92–94 LC-ωPBE-D3, ωB97X,106ωB97XD,107ωB97M-V81).

For each M+(H2O) dimer, the test set contains 400 CCSD(T)-F12b dimer configurations with M+ − O distances between 1.0 Å and 8.0 Å and interaction energies up to 50 kcal/mol. For these test sets, the MB-nrg RMSDs for the M+(H2O) dimers are 0.231 kcal/mol, 0.093 kcal/mol, 0.143 kcal/mol, 0.051 kcal/mol, and 0.051 kcal/mol, for M+ = Li+, Na+, K+, Rb+, and Cs+, respectively. A comparison of the RMSDs associated with the different models is shown in Fig. 5, while an analysis of the RMSDs as a function of the threshold in the interaction energies and basis sets is presented in the supplementary material.

FIG. 5.

RMSDs associated with the 2B energies of the test set with respect to CCSD(T)(-F12b)/CBS for various DFT models and the MB-nrg potentials.

FIG. 5.

RMSDs associated with the 2B energies of the test set with respect to CCSD(T)(-F12b)/CBS for various DFT models and the MB-nrg potentials.

Close modal

As expected, larger RMSD values are associated with GGA functionals than meta-GGA and hybrid functionals. Interestingly, it appears that the addition of a dispersion correction has either a negative or negligible effect on the performance of all functionals, resulting in larger RMSD values in most cases. Among all functionals considered in this study, ωB97M-V performs best on average and is the only density functional with RMSDs of less than 0.3 kcal/mol for all five dimers. PBE0 performs second best and affords RMSDs under 0.3 kcal/mol for all of the dimers except Li+(H2O). As noted in previous studies, the accuracy of PBE0 deteriorates (by nearly a factor of two) upon addition of the dispersion correction (e.g., see Ref. 57). The two range-separated functionals without dispersion corrections display opposite behavior for different ions: LC-ωPBE performs well for the smaller cations (Li+ and Na+) but not for the larger cations (K+, Rb+, and Cs+), while the reverse is true for ωB97X. It has recently been argued that revPBE and revPBE0 provide a reliable description of aqueous solutions.108–111 The present analysis of fundamental 2B interactions, however, demonstrates that both functionals actually perform worse than the original PBE and PBE0 functionals for most of the M+(H2O) dimers. Finally, the most recent ωB97M-V functional improves by more than a factor of two upon the related ωB97X-D functional for all the cations. For all M+(H2O) dimers, the MB-nrg PEFs are always associated with the lowest RMSDs, which, as mentioned above, is effectively uniform and independent of the alkali metal ion.

Two general observations can be made from the analysis of Fig. 5. First, for any given functional, the RMSD increases upon inclusion of the D3 empirical dispersion correction, and second, the performance of each functional varies in an unpredictable way with the nature of the cation, with the exception of ωB97M-V, which gives uniform RMSD values for the five alkali metal ions. Direct insights into the overall shape of each M+(H2O) PES (Potential Energy Surface) can be obtained from the analysis of the corresponding vibrational frequencies. Figure 6 shows the deviations of the harmonic frequencies of the water bend, and symmetric and asymmetric stretches calculated with both DFT models and MB-nrg PEFs relative to the CCSD(T)(-F12b) reference values. Contrary to what was found for the interaction energies, the inclusion of the D3 dispersion correction has a minimal effect on the vibrational motion of the water molecule. All functionals predict lower harmonic frequencies, with redshifts up to 60 cm−1, for the bending vibration. Different behavior is instead found for both symmetric and asymmetric stretches, with all GGA and meta-GGA functionals presented in this analysis predicting significantly lower frequencies, with redshifts up to 150 cm−1, and all hybrid functionals predicting higher harmonic frequencies, with blueshifts up to 80 cm−1, with the exception of ωB97M-V, which predicts a slightly redshifted frequency for the symmetric and asymmetric stretches in the Li+(H2O) and Na+(H2O) dimers. The primary exception to this trend is the hybrid B3LYP-D3 functional that underestimates both harmonic frequencies by 20 cm−1. While revPBE0, B3LYP, and PBE0 provide a consistently reasonable description of the vibrational frequencies of the water molecule in the five M+(H2O) dimers, ωB97M-V appears to provide, overall, the closest agreement with the CCSD(T)(-F12b) reference values among all functionals considered in this study, independent of the cation being studied. Interestingly, while ωB97X-D is one of the better functionals for reproducing the harmonic frequency of the bending vibrations, it significantly overestimates both stretching vibrations.

FIG. 6.

Histograms showing the deviations in the (a) bending, (b) symmetric stretch, and (c) asymmetric stretch of the water molecule, in the optimized geometry of each M+(H2O) dimer (M+ = Li+, Na+, K+, Rb+, Cs+), using different density functional methods and the MB-nrg potentials from the CCSD(T)(-F12b) normal modes.

FIG. 6.

Histograms showing the deviations in the (a) bending, (b) symmetric stretch, and (c) asymmetric stretch of the water molecule, in the optimized geometry of each M+(H2O) dimer (M+ = Li+, Na+, K+, Rb+, Cs+), using different density functional methods and the MB-nrg potentials from the CCSD(T)(-F12b) normal modes.

Close modal

As for the interaction energies, the analysis of the harmonic frequencies in Fig. 6 demonstrates that the performance of most functionals is particularly sensitive to the nature of the alkali metal ion. On the other hand, independent of the identity of M+, the MB-nrg PEFs accurately reproduce the coupled cluster reference values, with the deviations never being larger than 10 cm−1 for the bending and 20 cm−1 for the stretching vibrations.

In this study, we have introduced the two-body term (V2B) of the MB-nrg PEFs describing the interactions between alkali metal ions and water. By construction, V2B represents the full-dimensional potential energy surface of each M+(H2O) dimer, where M+ = Li+, Na+, K+, Rb+, and Cs+. Building upon our previous studies on many-body effects in water54–56 and halide-water systems,65,V2B is derived entirely from high-level electronic structure data and represented by combining a permutationally invariant polynomial at short range with an underlying classical description of many-body electrostatic and two-body dispersion interactions at all intermolecular separations. As described in Ref. 65, all water-water interactions in the MB-nrg PEFs are represented by the many-body MB-pol PEF that accurately predicts structural, thermodynamic, dynamical, and spectroscopic properties of water from the gas to the condensed phase.57,58

The accuracy of the 2B term of alkali metal ion-water MB-nrg PEFs has been assessed through extensive analysis of both interaction energies and harmonic frequencies for all five M+(H2O) dimers. Comparisons with different ab initio methods (i.e., coupled cluster, MP2, and several DFT models) and polarizable force fields (i.e., AMOEBA and TTM-nrg) demonstrate that the MB-nrg PEFs provide, in all cases, the highest accuracy, quantitatively reproducing the corresponding coupled cluster reference values. This analysis also shows that the addition of dispersion corrections to popular DFT models tends to have a negative effect on the energetics, while leaving the harmonic frequencies effectively unchanged. Although the magnitude of the DFT deviations from the reference values varies depending on the M+(H2O) dimer, the range-separated hybrid ωB97M-V functional overall displays the best performance among all the functionals considered in this study for both interaction energies and harmonic frequencies. On the other hand, both the AMOEBA and TTM-nrg polarizable force fields exhibit noticeable limitations, especially at short-range, as they are not able to correctly represent quantum-mechanical effects (e.g., charge transfer and penetration, and Pauli repulsion) arising from the overlap of the monomer electron densities.

See supplementary material for new TTM-nrg parameterization using CCSD(T)/CBS with correlation consistent basis sets as reference, frequency values for optimized H2O molecules using all functionals described in Sec. III B, frequency values for optimized M+(H2O) dimers, where M+ = Li+, Na+, K+, Rb+, Cs+, using all functionals described in Sec. III B, and configurations of the optimized dimers using all the above methods.

We thank Dr. Grant Hill and Dr. Kirk Peterson for providing us with the unpublished correlation consistent basis sets for potassium, rubidium, and cesium. We also thank Dr. Sandra Brown and Shelby Straight for valuable discussions. This research was supported by the National Science Foundation Center for Chemical Innovation “Center for Aerosol Impacts on Climate and the Environment” (Grant No. CHE-1305427). All expensive calculations for the training set generation (around 300 000 SUs) were performed using resources provided by the Open Science Grid,112,113 which is supported by the National Science Foundation and the U.S. Department of Energy’s Office of Science. We thank Edgar Fajardo for his help and technical support on the use of the software in the grid and the Physics Computing Facilities of UCSD for granting us access to the grid. The DFT calculations were run in the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (Grant No. ACI-1053575, Allocation No. TG-CHE110009) as well as at the Triton Shared Computing Cluster (TSCC) at the San Diego Supercomputer Center.

1.
Y.
Marèchal
,
The Hydrogen Bond and the Water Molecule: The Physics and Chemistry of Water, Aqueous and Bio-Media
(
Elsevier
,
Amsterdam
,
2006
).
2.
P.
Ball
,
Chem. Rev.
108
,
74
(
2008
).
3.
I.
Restrepo-Angulo
,
A.
De Vizcaya-Ruiz
, and
J.
Camacho
,
J. Appl. Toxicol.
30
,
497
(
2010
).
4.
P. L.
Nostro
and
B. W.
Ninham
,
Chem. Rev.
112
,
2286
(
2012
).
5.
Y.
Marcus
,
Chem. Rev.
109
,
1346
(
2009
).
6.
L.
Piatkowski
,
Z.
Zhang
,
E. H. G.
Backus
,
H. J.
Bakker
, and
M.
Bonn
,
Nat. Commun.
5
,
4083
(
2014
).
7.
V.
Migliorati
,
F.
Sessa
,
G.
Aquilanti
, and
P.
D’Angelo
,
J. Chem. Phys.
141
,
044509
(
2014
).
8.
D. J.
Tobias
,
A. C.
Stern
,
M. D.
Baer
,
Y.
Levin
, and
C. J.
Mundy
,
Annu. Rev. Phys. Chem.
64
,
339
(
2013
).
9.
D. E.
Smith
and
L. X.
Dang
,
J. Chem. Phys.
100
,
3757
(
1994
).
10.
A. W.
Omta
,
M. F.
Kropman
,
S.
Woutersen
, and
H. J.
Bakker
,
Science
301
,
347
(
2003
).
11.
D.
Spångberg
,
R.
Rey
,
J. T.
Hynes
, and
K.
Hermansson
,
J. Phys. Chem. B
107
,
4470
(
2003
).
12.
M. F.
Kropman
and
H. J.
Bakker
,
J. Am. Chem. Soc.
126
,
9135
(
2004
).
13.
J. D.
Smith
,
R. J.
Saykally
, and
P. L.
Geissler
,
J. Am. Chem. Soc.
129
,
13847
(
2007
).
14.
S.
Park
and
M.
Fayer
,
Proc. Natl. Acad. Sci. U. S. A.
104
,
16731
(
2007
).
15.
M. D.
Fayer
,
D. E.
Moilanen
,
D.
Wong
,
D. E.
Rosenfeld
,
E. E.
Fenn
, and
S.
Park
,
Acc. Chem. Res.
42
,
1210
(
2009
).
16.
D. L.
Bostick
and
C. L.
Brooks
,
J. Am. Chem. Soc.
132
,
13185
(
2010
).
17.
K. J.
Tielrooij
,
N.
Garcia-Araez
,
M.
Bonn
, and
H. J.
Bakker
,
Science
328
,
1006
(
2010
).
18.
S.
Funkner
,
G.
Niehues
,
D. A.
Schmidt
,
M.
Heyden
,
G.
Schwaab
,
K. M.
Callahan
,
D. J.
Tobias
, and
M.
Havenith
,
J. Am. Chem. Soc.
134
,
1030
(
2012
).
19.
P. L.
Geissler
,
Annu. Rev. Phys. Chem.
64
,
317
(
2013
).
20.
P.
Jungwirth
and
D. J.
Tobias
,
Chem. Rev.
106
,
1259
(
2006
).
21.
G.
Lamoureux
and
B.
Roux
,
J. Phys. Chem. B
110
,
3308
(
2006
).
22.
A.
Warshel
,
M.
Kato
, and
A. V.
Pisliakov
,
J. Chem. Theory Comput.
3
,
2034
(
2007
).
23.
T.
Peng
,
T.-M.
Chang
,
X.
Sun
,
A. V.
Nguyen
, and
L. X.
Dang
,
J. Mol. Liq.
173
,
47
(
2012
).
24.
Y.
Shi
,
Z.
Xia
,
J.
Zhang
,
R.
Best
,
C.
Wu
,
J. W.
Ponder
, and
P.
Ren
,
J. Chem. Theory Comput.
9
,
4046
(
2013
).
25.
P. T.
Kiss
and
A.
Baranyai
,
J. Chem. Phys.
141
,
114501
(
2014
).
26.
M.
Kohagen
,
E.
Pluhařová
,
P. E.
Mason
, and
P.
Jungwirth
,
J. Phys. Chem. Lett.
6
,
1563
(
2015
).
27.
A. A.
Hassanali
,
J.
Cuny
,
V.
Verdolino
, and
M.
Parrinello
,
Philos. Trans. R. Soc., A
372
,
20120482
(
2014
).
28.
C. D.
Wick
,
L. X.
Dang
, and
P.
Jungwirth
,
J. Chem. Phys.
125
,
024706
(
2006
).
29.
C. D.
Wick
,
I.-F. W.
Kuo
,
C. J.
Mundy
, and
L. X.
Dang
,
J. Chem. Theory Comput.
3
,
2002
(
2007
).
30.
A. V.
Marenich
,
R. M.
Olson
,
A. C.
Chamberlin
,
C. J.
Cramer
, and
D. G.
Truhlar
,
J. Chem. Theory Comput.
3
,
2055
(
2007
).
31.
C. D.
Wick
and
S. S.
Xantheas
,
J. Phys. Chem. B
113
,
4141
(
2009
).
32.
C. D.
Wick
,
J. Chem. Phys.
131
,
084715
(
2009
).
33.
E.
Guàrdia
,
I.
Skarmoutsos
, and
M.
Masia
,
J. Chem. Theory Comput.
5
,
1449
(
2009
).
34.
H.
Yu
,
T. W.
Whitfield
,
E.
Harder
,
G.
Lamoureux
,
I.
Vorobyov
,
V. M.
Anisimov
,
A. D.
MacKerell
, and
B.
Roux
,
J. Chem. Theory Comput.
6
,
774
(
2010
).
35.
C.
Caleman
,
J. S.
Hub
,
P. J.
van Maaren
, and
D.
van der Spoel
,
Proc. Natl. Acad. Sci. U. S. A.
108
,
6838
(
2011
).
36.
L.
Sun
,
X.
Li
,
Y.
Tu
, and
H.
Agren
,
Phys. Chem. Chem. Phys.
17
,
4311
(
2015
).
37.
J. M.
Lisy
,
Int. Rev. Phys. Chem.
16
,
267
(
1997
).
38.
J.-J.
Max
and
C.
Chapados
,
J. Chem. Phys.
115
,
2664
(
2001
).
39.
D. J.
Miller
and
J. M.
Lisy
,
J. Am. Chem. Soc.
130
,
15381
(
2008
).
40.
I. A.
Heisler
and
S. R.
Meech
,
Science
327
,
857
(
2010
).
41.
A.
Grossfield
,
P.
Ren
, and
J. W.
Ponder
,
J. Am. Chem. Soc.
125
,
15671
(
2003
).
42.
J.-P.
Piquemal
,
L.
Perera
,
G. A.
Cisneros
,
P.
Ren
,
L. G.
Pedersen
, and
T. A.
Darden
,
J. Chem. Phys.
125
,
054511
(
2006
).
43.
J. C.
Wu
,
J.-P.
Piquemal
,
R.
Chaudret
,
P.
Reinhardt
, and
P.
Ren
,
J. Chem. Theory Comput.
6
,
2059
(
2010
).
44.
J. Y.
Xiang
and
J. W.
Ponder
,
J. Comput. Chem.
34
,
739
(
2013
).
45.
A.
Marjolin
,
C.
Gourlaouen
,
C.
Clavaguéra
,
P. Y.
Ren
,
J. C.
Wu
,
N.
Gresh
,
J.-P.
Dognon
, and
J.-P.
Piquemal
,
Theor. Chem. Acc.
131
,
1198
(
2012
).
46.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
van der Avoird
,
Science
315
,
1249
(
2007
).
47.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
van der Avoird
,
J. Chem. Phys.
128
,
094313
(
2008
).
48.
R.
Bukowski
,
K.
Szalewicz
,
G. C.
Groenenboom
, and
A.
van der Avoird
,
J. Chem. Phys.
128
,
094314
(
2008
).
49.
Y.
Wang
,
B. C.
Shepler
,
B. J.
Braams
, and
J. M.
Bowman
,
J. Chem. Phys.
131
,
054511
(
2009
).
50.
Y.
Wang
,
X.
Huang
,
B. C.
Shepler
,
B. J.
Braams
, and
J. M.
Bowman
,
J. Chem. Phys.
134
,
094509
(
2011
).
51.
Y.
Wang
and
J. M.
Bowman
,
J. Chem. Phys.
134
,
154510
(
2011
).
52.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
,
J. Phys. Chem. Lett.
3
,
3765
(
2012
).
53.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
,
1103
(
2013
).
54.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
,
5395
(
2013
).
55.
V.
Babin
,
G. R.
Medders
, and
F.
Paesani
,
J. Chem. Theory Comput.
10
,
1599
(
2014
).
56.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
,
J. Chem. Theory Comput.
10
,
2906
(
2014
).
57.
G. A.
Cisneros
,
K. T.
Wikfeldt
,
L.
Ojamäe
,
J.
Lu
,
Y.
Xu
,
H.
Torabifard
,
A. P.
Bartók
,
G.
Csányi
,
V.
Molinero
, and
F.
Paesani
,
Chem. Rev.
116
,
7501
(
2016
).
58.
F.
Paesani
,
Acc. Chem. Res.
49
,
1844
(
2016
).
59.
S. K.
Reddy
,
S. C.
Straight
,
P.
Bajaj
,
C.
Huy Pham
,
M.
Riera
,
D. R.
Moberg
,
M. A.
Morales
,
C.
Knight
,
A. W.
Götz
, and
F.
Paesani
,
J. Chem. Phys.
145
,
194504
(
2016
).
60.
J. O.
Richardson
,
C.
Pérez
,
S.
Lobsiger
,
A. A.
Reid
,
B.
Temelso
,
G. C.
Shields
,
Z.
Kisiel
,
D. J.
Wales
,
B. H.
Pate
, and
S. C.
Althorpe
,
Science
351
,
1310
(
2016
).
61.
W. T. S.
Cole
,
J. D.
Farrell
,
D. J.
Wales
, and
R. J.
Saykally
,
Science
352
,
1194
(
2016
).
62.
G. R.
Medders
and
F.
Paesani
,
J. Chem. Theory Comput.
11
,
1145
(
2015
).
63.
G. R.
Medders
and
F.
Paesani
,
J. Am. Chem. Soc.
138
,
3912
(
2016
).
64.
C. H.
Pham
,
S. K.
Reddy
,
K.
Chen
,
C.
Knight
, and
F.
Paesani
,
J. Chem. Theory Comput.
13
,
1778
(
2017
).
65.
P.
Bajaj
,
A. W.
Götz
, and
F.
Paesani
,
J. Chem. Theory Comput.
12
,
2698
(
2016
).
66.
D.
Hankins
,
J.
Moskowitz
, and
F.
Stillinger
,
J. Chem. Phys.
53
,
4544
(
1970
).
67.
T. B.
Adler
,
G.
Knizia
, and
H. J.
Werner
,
J. Chem. Phys.
127
,
221106
(
2007
).
68.
G.
Knizia
,
T. B.
Adler
, and
H. J.
Werner
,
J. Chem. Phys.
130
,
054104
(
2009
).
69.
J. G.
Hill
,
K. A.
Peterson
,
G.
Knizia
, and
H.-J.
Werner
,
J. Chem. Phys.
131
,
194105
(
2009
).
70.
U.
Góra
,
R.
Podeszwa
,
W.
Cencek
, and
K.
Szalewicz
,
J. Chem. Phys.
135
,
224102
(
2011
).
71.
T. H.
Dunning
,
J. Chem. Phys.
90
,
1007
(
1989
).
72.
R. A.
Kendall
,
T. H.
Dunning
, and
R. J.
Harrison
,
J. Chem. Phys.
96
,
6796
(
1992
).
73.
D. E.
Woon
and
T. H.
Dunning
,
J. Chem. Phys.
103
,
4572
(
1995
).
74.
K.
Peterson
, private communication (2016).
75.
I. S.
Lim
,
P.
Schwerdtfeger
,
B.
Metz
, and
H.
Stoll
,
J. Chem. Phys.
122
,
104103
(
2005
).
76.
D. E.
Woon
and
T. H.
Dunning
,
J. Chem. Phys.
100
,
2975
(
1994
).
77.
I. S.
Lim
,
J. K.
Laerdahl
, and
P.
Schwerdtfeger
,
J. Chem. Phys.
116
,
172
(
2002
).
78.
Y.
Singh
,
B. K.
Sahoo
, and
B. P.
Das
,
Phys. Rev. A
88
,
062504
(
2013
).
79.
H.-J.
Werner
,
P.
Knowles
,
G.
Knizia
,
R.
Manby
,
M.
Schütz
 et al, molpro, version 2012.1, a package of ab initio programs,
2012
, see http://www.molpro.net.
80.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
B.
Mennucci
,
G. A.
Petersson
,
H.
Nakatsuji
,
M.
Caricato
,
X.
Li
,
H. P.
Hratchian
,
A. F.
Izmaylov
,
J.
Bloino
,
G.
Zheng
,
J. L.
Sonnenberg
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M.
Bearpark
,
J. J.
Heyd
,
E.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
N.
Rega
,
J. M.
Millam
,
M.
Klene
,
J. E.
Knox
,
J. B.
Cross
,
V.
Bakken
,
C.
Adamo
,
J.
Jaramillo
,
R.
Gomperts
,
R. E.
Stratmann
,
O.
Yazyev
,
A. J.
Austin
,
R.
Cammi
,
C.
Pomelli
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
V. G.
Zakrzewski
,
G. A.
Voth
,
P.
Salvador
,
J. J.
Dannenberg
,
S.
Dapprich
,
A. D.
Daniels
,
Ö.
Farkas
,
J. B.
Foresman
,
J. V.
Ortiz
,
J.
Cioslowski
, and
D. J.
Fox
, gaussian 09, Revision D.01, Gaussian, Inc., Wallingford, CT,
2009
.
81.
N.
Mardirossian
and
M.
Head-Gordon
,
J. Chem. Phys.
142
,
074111
(
2015
).
82.
Y.
Shao
,
Z.
Gan
,
E.
Epifanovsky
,
A. T. B.
Gilbert
,
M.
Wormit
,
J.
Kussmann
,
A. W.
Lange
,
A.
Behn
,
J.
Deng
,
X.
Feng
,
D.
Ghosh
,
M.
Goldey
,
P. R.
Horn
,
L. D.
Jacobson
,
I.
Kaliman
,
R. Z.
Khaliullin
,
T.
Kús
,
A.
Landau
,
J.
Liu
,
E. I.
Proynov
,
Y. M.
Rhee
,
R. M.
Richard
,
M. A.
Rohrdanz
,
R. P.
Steele
,
E. J.
Sundstrom
,
H. L.
Woodcock
 III
,
P. M.
Zimmerman
,
D.
Zuev
,
B.
Albrecht
,
E.
Alguire
,
B.
Austin
,
G. J. O.
Beran
,
Y. A.
Bernard
,
E.
Berquist
,
K.
Brandhorst
,
K. B.
Bravaya
,
S. T.
Brown
,
D.
Casanova
,
C.-M.
Chang
,
Y.
Chen
,
S. H.
Chien
,
K. D.
Closser
,
D. L.
Crittenden
,
M.
Diedenhofen
,
R. A.
DiStasio
, Jr.
,
H.
Dop
,
A. D.
Dutoi
,
R. G.
Edgar
,
S.
Fatehi
,
L.
Fusti-Molnar
,
A.
Ghysels
,
A.
Golubeva-Zadorozhnaya
,
J.
Gomes
,
M. W. D.
Hanson-Heine
,
P. H. P.
Harbach
,
A. W.
Hauser
,
E. G.
Hohenstein
,
Z. C.
Holden
,
T.-C.
Jagau
,
H.
Ji
,
B.
Kaduk
,
K.
Khistyaev
,
J.
Kim
,
J.
Kim
,
R. A.
King
,
P.
Klunzinger
,
D.
Kosenkov
,
T.
Kowalczyk
,
C. M.
Krauter
,
K. U.
Lao
,
A.
Laurent
,
K. V.
Lawler
,
S. V.
Levchenko
,
C. Y.
Lin
,
F.
Liu
,
E.
Livshits
,
R. C.
Lochan
,
A.
Luenser
,
P.
Manohar
,
S. F.
Manzer
,
S.-P.
Mao
,
N.
Mardirossian
,
A. V.
Marenich
,
S. A.
Maurer
,
N. J.
Mayhall
,
C. M.
Oana
,
R.
Olivares-Amaya
,
D. P.
O’Neill
,
J. A.
Parkhill
,
T. M.
Perrine
,
R.
Peverati
,
P. A.
Pieniazek
,
A.
Prociuk
,
D. R.
Rehn
,
E.
Rosta
,
N. J.
Russ
,
N.
Sergueev
,
S. M.
Sharada
,
S.
Sharmaa
,
D. W.
Small
,
A.
Sodt
,
T.
Stein
,
D.
Stück
,
Y.-C.
Su
,
A. J. W.
Thom
,
T.
Tsuchimochi
,
L.
Vogt
,
O.
Vydrov
,
T.
Wang
,
M. A.
Watson
,
J.
Wenzel
,
A.
White
,
C. F.
Williams
,
V.
Vanovschi
,
S.
Yeganeh
,
S. R.
Yost
,
Z.-Q.
You
,
I. Y.
Zhang
,
X.
Zhang
,
Y.
Zhou
,
B. R.
Brooks
,
G. K. L.
Chan
,
D. M.
Chipman
,
C. J.
Cramer
,
W. A.
Goddard
 III
,
M. S.
Gordon
,
W. J.
Hehre
,
A.
Klamt
,
H. F.
Schaefer
 III
,
M. W.
Schmidt
,
C. D.
Sherrill
,
D. G.
Truhlar
,
A.
Warshel
,
X.
Xua
,
A.
Aspuru-Guzik
,
R.
Baer
,
A. T.
Bell
,
N. A.
Besley
,
J.-D.
Chai
,
A.
Dreuw
,
B. D.
Dunietz
,
T. R.
Furlani
,
S. R.
Gwaltney
,
C.-P.
Hsu
,
Y.
Jung
,
J.
Kong
,
D. S.
Lambrecht
,
W.
Liang
,
C.
Ochsenfeld
,
V. A.
Rassolov
,
L. V.
Slipchenko
,
J. E.
Subotnik
,
T.
Van Voorhis
,
J. M.
Herbert
,
A. I.
Krylov
,
P. M. W.
Gill
, and
M.
Head-Gordon
,
Mol. Phys.
113
,
184
(
2015
).
83.
J. W.
Ponder
and
F. M.
Richards
,
J. Comput. Chem.
8
,
1016
(
1987
).
84.
C. E.
Kundrot
,
J. W.
Ponder
, and
F. M.
Richards
,
J. Comput. Chem.
12
,
402
(
1991
).
85.
P.
Ren
and
J. W.
Ponder
,
J. Phys. Chem. B
107
,
5933
(
2003
).
86.
M.
Riera
,
A. W.
Götz
, and
F.
Paesani
,
Phys. Chem. Chem. Phys.
18
,
30334
(
2016
).
87.
C. J.
Burnham
,
D. J.
Anick
,
P. K.
Mankoo
, and
G. F.
Reiter
,
J. Chem. Phys.
128
,
154519
(
2008
).
88.
K. T.
Tang
and
J. P.
Toennies
,
J. Chem. Phys.
80
,
3726
(
1984
).
89.
E. R.
Johnson
and
A. D.
Becke
,
J. Chem. Phys.
123
,
024101
(
2005
).
90.
F. O.
Kannemann
and
A. D.
Becke
,
J. Chem. Theory Comput.
6
,
1081
(
2010
).
91.
A.
Otero-de-la Roza
and
E. R.
Johnson
,
J. Chem. Phys.
138
,
204109
(
2013
).
92.
O. A.
Vydrov
and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
234109
(
2006
).
93.
O. A.
Vydrov
,
J.
Heyd
,
A. V.
Krukau
, and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
074106
(
2006
).
94.
O. A.
Vydrov
,
G. E.
Scuseria
, and
J. P.
Perdew
,
J. Chem. Phys.
126
,
154109
(
2007
).
95.
B. J.
Braams
and
J. M.
Bowman
,
Int. Rev. Phys. Chem.
28
,
577
(
2009
).
96.
A.
Tikhonov
,
Soviet Math. Dokl.
4
,
1035
(
1963
).
97.
E. M.
Mas
,
K.
Szalewicz
,
R.
Bukowski
, and
B.
Jeziorski
,
J. Chem. Phys.
107
,
4207
(
1997
).
98.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
99.
A. D.
Becke
,
Phys. Rev. A
38
,
3098
(
1988
).
100.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
101.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
78
,
1396
(
1996
).
102.
Y.
Zhang
and
W.
Yang
,
Phys. Rev. Lett.
80
,
890
(
1998
).
103.
J. M.
Tao
,
J. P.
Perdew
,
V. N.
Staroverov
, and
G. E.
Scuseria
,
Phys. Rev. Lett.
91
,
146401
(
2003
).
104.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
105.
A. D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
106.
J.-D.
Chai
and
M.
Head-Gordon
,
J. Chem. Phys.
128
,
084106
(
2008
).
107.
J.-D.
Chai
and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
10
,
6615
(
2008
).
108.
A.
Bankura
,
A.
Karmakar
,
V.
Carnevale
,
A.
Chandra
, and
M. L.
Klein
,
J. Phys. Chem. C
118
,
29401
(
2014
).
109.
L. B.
Skinner
,
M.
Galib
,
J. L.
Fulton
,
C. J.
Mundy
,
J. B.
Parise
,
V.-T.
Pham
,
G. K.
Schenter
, and
C. J.
Benmore
,
J. Chem. Phys.
144
,
134504
(
2016
).
110.
M.
Galib
,
M. D.
Baer
,
L. B.
Skinner
,
C. J.
Mundy
,
T.
Huthwelker
,
G. K.
Schenter
,
C. J.
Benmore
,
N.
Govind
, and
J. L.
Fulton
,
J. Chem. Phys.
146
,
084504
(
2017
).
111.
O.
Marsalek
and
T. E.
Markland
,
J. Phys. Chem. Lett.
8
,
1545
(
2017
).
112.
I.
Sfiligoi
,
D. C.
Bradley
,
B.
Holzman
,
P.
Mhashilkar
,
S.
Padhi
, and
F.
Wurthwein
,
2009 WRI World Congress on Computer Science and Information Engineering
(
IEEE
,
2009
), Vol. 2, pp.
428
432
.
113.
R.
Pordes
,
D.
Petravick
,
B.
Kramer
,
D.
Olson
,
M.
Livny
,
A.
Roy
,
P.
Avery
,
K.
Blackburn
,
T.
Wenaus
, and
F.
Würthwein
,
J. Phys.: Conf. Ser.
78
,
012057
(
2007
).

Supplementary Material