Comparative Assessment of Computational ... - ACS Publications

Article Cite This: J. Chem. Inf. Model. XXXX, XXX, XXX-XXX

pubs.acs.org/jcim

Comparative Assessment of Computational Methods for Free Energy Calculations of Ionic Hydration Haiyang Zhang,† Yang Jiang,‡ Hai Yan,*,† Ziheng Cui,‡ and Chunhua Yin*,† †

Department of Biological Science and Engineering, School of Chemistry and Biological Engineering, University of Science and Technology Beijing, 100083 Beijing, China ‡ Beijing Key Lab of Bioprocess, College of Life Science and Technology, Beijing University of Chemical Technology, Box 53, 100029 Beijing, China S Supporting Information *

ABSTRACT: Experimental observations for ionic hydration free energies are highly debated mainly due to the ambiguous absolute hydration free energy of proton, ΔGhyd * (H+). Hydration free energies (HFEs) of the 112 singly charged ions in the Minnesota solvation database were predicted by six methods with explicit and implicit solvent models, namely, thermodynamic integration (TI), energy representation module (ERmod), three-dimensional reference interaction site model (3DRISM), and continuum solvation models based on the quantum mechanical charge density (SMD) and on the Poisson−Boltzmann (PB) and generalized Born (GB) theories. Taking the solvent Galvani potential of water into account, the resulting real HFEs from TI calculations for the generalized Amber force field (GAFF) modeled ions best match the experiments based on ΔGhyd * (H+) = −262.4 kcal/mol (Randles Trans. Faraday Soc. 1956, 52, 1573−1581), in agreement with our previous work on charged amino acids (Zhang et al. J. Phys. Chem. Lett. 2017, 8, 2705−2712). The examined computational methods show an accuracy of ∼7 kcal/mol for the GAFF-modeled ions, except for SMD with a higher accuracy of ∼4 kcal/mol. A biased deficiency in modeling anionic compounds by GAFF is observed with a larger standard deviation (SD) of 9 kcal/mol than that for cations (SD ∼ 4 kcal/ mol). The relatively cheap ERmod and 3D-RISM methods reproduce TI results with good accuracy, although ERmod yields a systematic underestimation for cations by 9 kcal/mol; PB and GB generate relative (but not absolute) HFEs comparable to the TI predictions. Computational accuracy is found to be more limited by the accuracy of force fields rather than the models themselves. hydration,15−18 mainly due to the ambiguous hydration free energy (HFE) of the proton that the experiments for molecular ions depend on. Estimations of experimental HFEs are spread over a range of 10 kcal/mol. Marcus adopted a ΔG*hyd(H+) of −254.3 kcal/mol (−1056 kJ/mol minus a standard-state correction of ΔG0→* = 7.95 kJ/mol) for the experimental HFE determination of ions,15 whereas Noyes and Rosseinsky provided a different estimation16,17 for the same purpose based on Randles’s prediction19 of −262.4 kcal/mol (−260.5 kcal/ mol minus a ΔG0→* of 1.9 kcal/mol). The Marcus set was chosen by Merz and co-workers as reference for rational design of metal models,20−23 whereas Åqvist and Kamerlin used the set by Noyes16 to derive ion−water potential parameters of metal cations.24−26 Hünenberger et al. suggested a value of −1107.95 kJ/mol (−264.8 kcal/mol)27 and then adopted it as a reference to recalibrate interaction parameters of charged amino acid side chains, resulting in the GROMOS 54A8 parameter set. 8 Truhlar and co-workers followed the Tissandier’s prediction of −1112.45 kJ/mol (−265.9 kcal/ mol),28 a value often acknowledged as the best estimate of

1. INTRODUCTION Ions are ubiquitous in nature and play a fundamental role in numerous fields such as chemistry and biology. Metal ions aid, for instance, the assembly of sophisticated materials via coordinating with negatively charged linkers1 and are shown to take part in enzymatic catalysis via stabilizing reaction intermediates, activating nucleophiles, or acting as redox centers.2−4 Titratable amino acids are involved as well in promoting enzymatic catalysis via a change in protonated states such as the well-known Ser-His-Asp triad of lipases.5 The mechanisms underlying these phenomena are complex and not easy to probe experimentally. Even with computational methods, it is still a challenging endeavor to model ioncontaining systems.6 Computational models are developed typically to reproduce experimental observations for physical/chemical properties of molecules and liquids, of which the solvation free energy is an important quantity to target, usually with a given solvent of water, as did in the GROMOS 53A67 and 54A88 force field (FF) sets. If not chosen as a target for model developments, this quantity may be useful for assessing the developed models in turn.9−14 There is however no consensus to date in the literature regarding the experimental free energies of ionic © XXXX American Chemical Society

Received: August 11, 2017 Published: October 17, 2017 A

DOI: 10.1021/acs.jcim.7b00485 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling ΔG*hyd(H+),29 and this value was utilized for the Minnesota solvation database18,30 and the development of a universal solvation model (denoted as SMD).31 Using different sets of HFE references for FF development efforts of ions,20−23,25,26,32,33 as stated above, generates confusing and contradicting computational models probably, which necessitates “choosing” an experimental reference for calibration of the FF models. Another issue dealing with for instance ionic hydration is the complexity of correcting the artifacts that exist in the modeling of charged species in liquids. The artifacts have a contribution (i.e., a corrected term referred to as ΔGcor) to the ionic HFE and must be corrected properly, since the experiments were often performed at a different situation from the simulations.34 The group of Hünenberger presented a detailed analysis on the potential artifacts of ionic solvation, along with the proposed ex post correction schemes.34−36 These correction strategies lead to a methodology-independent prediction of ionic solvation from molecular simulations and have been proved useful for the calculation of single-ion solvation32,33,37 and binding free energies of charged species.38 Use of different correction schemes, however, has been reported to give different or even contradicting matches between experimental and calculated observations.8,24−26,32,33,37,39 For the ionic HFE estimations with or without corrections (ΔGcor), again, it necessitates “choosing” an experimental reference consistent with the used corrections. A very recent report on systematic HFE benchmarks of charged amino acid side chains with varied popular protein FFs indicated empirical best matches between ΔG*hyd(H+)-dependent experimental HFEs of ionic solvation and ΔGcor-dependent predictions from molecular simulations, although this does not settle up the ongoing debate of which ΔG*hyd(H+) is correct.39 Moreover, calculated solvation free energies of molecular ions depend on the techniques used in the simulations such as boundary conditions and summation schemes for potential evaluation at a particular site, and also on the cases where the correction strategies (ΔGcor) are applied or not.34,40,41 The performance of force field models and computational methods on ion-containing systems is still largely unknown therefore. Here we presented a systematic assessment of HFE calculations for a large set of 112 ions by six computational methods, including thermodynamic integration (TI), energy representation module (ERmod), three-dimensional reference interaction site model (3D-RISM), and continuum solvation models based on the quantum mechanical charge density (SMD) and on the Poisson−Boltzmann (PB) and generalized Born (GB) theories. Among these methods, TI is most accurate but with a heavy computational load; the others are relatively cheap and have shown great potential in the applications like biomolecular solvation and protein−ligand affinity predictions. Comparison of the calculated results with ΔGhyd * (H+)-dependent experimental observations was made attempting to investigate the apparent controversy in the experimental HFE of the proton. The achieved accomplishments from the state-ofart computational methods and the discrepancies/deficiencies in varied models were highlighted in detail as well. This work may be valuable for further effects of ionic model developments.

2012)18,30 were chosen as a test set including 60 anions and 52 cations. Their experimental hydration free energies (HFEs) in the database used the Tissandier et al.’s prediction of * (H+) = −265.9 kcal/mol28 as a reference by ΔGhyd * (ion) = ΔG hyd *,conv (ion) + qΔG hyd * (H+) ΔG hyd

(1)

where ΔG*hyd,conv(ion) represents the conventional hydration free energy. ΔGhyd *,conv(ion) can be determined, in principle, experimentally and is thus not as ambiguous as the absolute hydration free energy (HFE) of ΔG*hyd(ion) that depends on ΔG*hyd(H+). Table S1 gives four sets of experimental HFEs for the 112 ions, based on the generally adopted values of ΔGhyd * (H+) in the model developments of ions in computational chemistry, namely, −254.3,15 −262.4,19 −264.8,27 and −265.9 kcal/mol,28 respectively. Initial structures for these ions were also obtained from the Minnesota solvation database where these ions were optimized at the M06-2X/MG3S level of theory in the gas phase.18,30 2.2. Thermodynamic Integration. Palmer and co-workers have predicted the hydration free energies of 70 ions,37 a partial set of the 112 ions tested in this work, using TI.42 They modeled the ionic solutes using the generalized Amber force field (GAFF)43 with AM1-BCC charges44,45 in SPC/E46 water. Following a similar protocol, we here calculated the other 42 ions with a two-stage TI, i.e., decoupling Coulombic interactions first and then van der Waals interactions, in order to quantify the contributions from Coulombic (ΔGcoul) and van der Waals (ΔGvdW) interactions separately. The decoupling transitions for the Coulombic interactions were carried out along an alchemical reaction coordinate λ = 0, 0.25, 0.5, 0.75, and 1; for van der Waals interactions, λ = 0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, and 1. Each λ was subjected to a 2 ns production simulation at 298.15 K and the first 100 ps was discarded for equilibration. Each system contains only one single ion and 1024 SPC/E water molecules in a cubic box and, thus, has a unity net charge. Electrostatics interactions were computed using the particlemesh Ewald (PME) method47,48 with a real-space cutoff of 1.2 nm, and Lennard−Jones interactions were switched off from 0.9 to 1.2 nm smoothly. More details on the TI protocol have been presented in refs 10, 13, 37, and 49. All the simulations were performed with the GROMACS 5.0.4 package50−52 under full periodic boundary conditions (PBC) and the simulation setups were the same as in the work by Palmer et al.37 except that the Parrinello−Rahman algorithm53,54 was applied instead of a Berendsen barostat55 for pressure coupling at 1 bar. The Bennett acceptance ratio (BAR) method was used to compute the HFEs from the TI statistics.56 2.3. Energy Representation Module. After an energy minimization and NPT equilibration, normal MD simulations with a similar protocol to TI (i.e., free energy calculation module in GROMACS was turned off) were implemented to simulate a solution system composed of one single ion (modeled by GAFF/AM1-BCC) and 1024 SPC/E water molecules in a cubic box. A pure-solvent system of 1024 water molecules in the absence of the solute ion was simulated as well for reference. Production simulations for the solution and solvent systems were run for 1 ns and the trajectories (instantaneous configurations) was saved to disk every 50 fs (2 × 104 snapshots in total) and 200 fs (5 × 103 snapshots) for solution and solvent systems, respectively. An isolated solute was also simulated for 100 ns in vacuum with the snapshot

2. COMPUTATIONAL METHODS 2.1. Selection of Ionic Compounds. All the 112 singly charged ions in the Minnesota solvation database (version B


Article

Journal of Chemical Information and Modeling

(corresponding to water), respectively. Two methods for computing nonpolar solvation free energy were considered: one is a single term linearly proportional to the solvent accessible surface area (inp = 1 in the PBSA protocol), denoted as PB/SASA; the other is modeled by a sum of cavity (repulsive) and dispersion (attractive) terms (inp = 2),78 denoted as PB/CD here. 2.7. GBSA. Generalized Born predictions were based on the 10 ns vacuum and GB simulations with the GROMACS 5.0.4 package.50−52 Four GB models supported in GROMACS79 were examined, namely, Still,80 HCT,81 OBCI, and OBCII.82 The polar part (ΔGpol) was computed from the generalized Born equation;79,80 the nonpolar part (ΔGnp) amounts to solvent accessible surface area multiplied by a surface tension, similar to PB/SASA in AMBER,70 and was derived directly from a simple ACE-type approximation.83 van der Waals radii for evaluation of effective Born radii were taken from the TINKER84 implementation of GB implicit solvation models. The “mbondi” radii recommended for use with HCT and the “mbondi2” radii for two OBC variants were taken from the AmberTools (version 15)62 and examined as well in the GB simulations for comparison with the radii from TINKER.84 Note that for HCT and two variants of OBC models, GROMACS and AMBER share the same overlap scaling factors that were also directly extracted from TINKER.84 These radii and scaling factors are given in Table S2 in the Supporting Information (SI). Following Mobely’s work,12 the original vacuum and GB simulation trajectories were reprocessed in GB and vacuum, respectively, in order to obtain the potential energy difference needed for hydration free energy predictions with the Bennett acceptance ratio (BAR) method.56 2.8. Free Energy Corrections. In a physical process of ionic solvation, the single ion is transferring from a large distance from the solvent (i.e., vacuum), across the air/liquid surface, and then into the bulk solution.37 Ionic solvation free energy can be formally decomposed into bulk and surface contributions. Free energy calculations as done typically under full periodic boundary conditions (PBC), exempt of any boundary between air and liquid, belong to the bulk contribution and yields so-called intrinsic free energies (ΔGintr). The surface contribution (ΔGsurf) corresponding to crossing the vacuum/bulk interface is required to be added leading to the physical (real) free energy of ionic solvation (eq 2).37,85

being saved every 100 fs, and these snapshots were used as test particles to insert into a pure-solvent system at random position with random orientation. The test-particle insertion was performed 1000 times per pure-solvent configuration sampled for producing a set of solute−solvent configurations to construct energy distribution functions for ERmod analysis.57,58 Combining the solute−solvent pair distribution functions of a “real” system (i.e., solution) and that of a “pseudo” system (i.e., insertion of solute into a pure solvent), the solvation free energy of the solute ion can be determined approximately with the ERmod toolkit (version 0.3.5).57,58 A long-range correction for the Lennard-Jones interaction between the solute and solvent was considered explicitly in the ERmod calculation. Solution and solute systems with atomic charges of the solute set to zero were simulated as well to obtain the nonpolar contribution to the solvation by ERmod.57,58 2.4. Three-Dimensional Reference Interaction Site Model. A semiexplicit model of 3D-RISM with a partial series of expansion of the order 3 (PSE-3) of the hypernetted-chain (HNC) closure plus pressure correction (denoted as 3DRISM/PC+)37,59−61 was used to estimate hydration free energy of the 112 ions modeled by GAFF/AM1-BCC in SPC/E water. With the aid of AmberTools 15 package,62 water susceptibility functions were computed first using dielectrically consistent 1D-RISM,63 followed by 3D-RISM calculations with the rism3d.snglpnt program64−66 for prediction of hydration free energy. All the calculation setups followed the work of Palmer et al.,37 and these were done automatically by a Python script (https://github.com/MTS-Strathclyde/PC_plus). For comparison with AM1-BCC, restrained electrostatic potential (RESP)67 charges of the ionic compounds were also computed using the Gaussian 09 software68 at the Hartree−Fock level with the 6-31G* basis set,69 a popular scheme for the GAFF parametrization of organic molecules in the AMBER software.70 3D-RISM/PC+ calculations were carried out for the 112 GAFF/RESP ions as well in TIP3P71 water. Note that in 3DRISM/PC+,37 original water models of SPC/E46 and TIP3P71 were modified slightly via adding Lennard−Jones interactions to hydrogen for convergence improvement, denoted as coincident SPC/E (cSPC/E, σH = 1.1658 Å and εH = 0.01553 kcal/mol) and coincident TIP3P (cTIP3P, σH = 1.2363 Å and εH = 0.01520 kcal/mol), respectively.65 Such modifications were only used for solvent water models in the 3D-RISM calculations, and all the force field parameters of solute molecules were taken from GAFF.43 2.5. Universal Solvation Model (SMD). Quantum mechanical SMD predictions were computed with the Gaussian 09 software.68 The optimized structures at M06-2X/MG3S level in the gas phase30 were used as initial coordinates of the solutes. Another optimization and frequency calculations for the 112 ions were carried out in the liquid phase (i.e., with SMD water model)31 and in the gas phase separately. Three levels of theory including Hartree−Fock (HF), B3LYP,72−74 and M06-2X75 and two basis sets of 6-31+G(d,p)76 and ccpVTZ77 were examined, resulting in six combinations in total. Hydration free energy was defined as the difference in the free energy of solute calculated in the liquid phase and in the gas phase.14 2.6. PBSA. A grid-based finite difference solution to the Poisson−Boltzmann (PB) equation implemented in the AMBER 11 software70 was conducted to obtain PB-based HFE predictions. Temperature was set to 298.15 K, and interior and exterior dielectric constants were set to 1 and 78.46

ΔGreal = ΔGintr + ΔGsurf = ΔGintr + qØG

(2)

where q is the ion net charge and ØG is known as the solvent Galvani potential jump in air to liquid direction. For the intrinsic solvation free energies in classical atomistic simulation with PBC, the calculations with the particle-based summation (P-scheme) and the molecule-based summation (M-scheme) for the potential evaluation at the ionic site always differ approximately by an offset of ΔGDSC (an ad hoc term for discrete solvent correction introduced by Hünenberger et al.36), as given in eq 3. ΔGintr,P = ΔGintr,M − ΔG DSC

(3)

This difference arises because the M-summation cutoff sphere acts like an artificial vapor/liquid boundary.40 For the calculations of surface potential jump, the results between the P-scheme and M-scheme differ approximately by an offset of ΔGDSC as well (eq 4).34 qØG,P = qØG,M + ΔG DSC C

(4) DOI: 10.1021/acs.jcim.7b00485 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling Since eqs 3 and 4 have an opposite sign for ΔGDSC, one then has ΔGreal = ΔGintr,P + qØG,P = ΔGintr,M + qØG,M

without free energy corrections for all the six methods are given in Table S4 in the SI.

3. RESULTS AND DISCUSSION Solvation models examined in this study can be classified into three categories: explicit, semiexplicit, and implicit. TI and ERmod methods use an explicit description of water molecules and fall within explicit solvation models. 3D-RISM treats solvent with density correlation functions, capturing solvation effects that may be ignored by implicit solvation models, and therefore serves as a semiexplicit solvent model.37,61 SMD, PB, and GB methods use a continuous medium with certain dielectric and interfacial properties to represent the solvent implicitly and are therefore considered as implicit solvation models. In this section we compare the calculations with experimental observations and if possible with the more accurate TI results, discuss the performance for each method, and then present a systematic assessment for all the solvation models. 3.1. Explicit Solvent. As shown in section 2.1, experimental HFEs of ions depend on the HFE value of the proton, ΔGhyd * (H+), and yet there is no consensus in the literature regarding this value. The accepted and widely adopted value in the model development efforts ranges from −265 to −250 kcal/mol,23 and diverse free energy correction strategies that produce for instance P-scheme intrinsic,20−23 real,37,40 Mintrinsic 8,32,33 HFEs have been used to target varied ΔGhyd * (H+). Here we examine experimental ionic HFE references based on a large range of ΔG*hyd(H+) from −275 to −245 kcal/mol and compare them with the TI calculated results of the 112 ions (Table S3). Figure 1 shows the correlation coefficient (R) and SD from experimental references as a function of the referenced ΔGhyd * (H+). For the P-scheme intrinsic HFEs (Figure 1a), with the ΔGhyd * (H+) increasing from −275 to −247 kcal/mol, the correlation (R) between calculated and experimental observations gets stronger and R approaches to 0.96, while the SD from experiments reduces approaching the limit of 7 kcal/mol. When ΔGhyd * (H+) exceeds −247 kcal/mol, SD increases gradually and no significant changes in R are observed. For the real HFEs (Figure 1b), SD reduced to 7 kcal/mol and R increases to the limit of 0.88 as ΔGhyd * (H+) + decreases when ΔGhyd * (H ) is larger than −261 kcal/mol. For ΔG*hyd(H+) < −261 kcal/mol, SD increases and R drops sharply. For the M-scheme intrinsic HFEs (Figure 1c), the SD and R profiles display similar landscapes to that of the real HFEs; however, the smallest SD and strongest correlation (R = 0.84) take place at ΔG*hyd(H+) = −267 kcal/mol. From Figure 1, it can be seen that the accuracy of TI calculations for the 112 ions with GAFF amounts to around 7 kcal/mol, which appears not to depend on the choice of the referenced ΔGhyd * (H+) and free energy correction strategies. The partial set of 70 ions tested by Palmer and co-workers (with the same TI protocol as in this work)37 shows a higher accuracy of 3 kcal/mol compared to the full set of 112 ions (SD ∼ 7 kcal/mol), which is confirmed by a similar analysis as above (Figure S1 in the SI). However, the smallest SD and highest R occur at −250, −264, and −270 kcal/mol for P-scheme intrinsic, real, and M-scheme intrinsic HFEs, respectively (Figure S1). Palmer et al. adopted the experimental reference based on ΔG*hyd(H+) = −265.9 kcal/mol for comparison,37 while there are no significant differences in such comparisons when choosing as a reference the Randles’s ΔGhyd * (H+) of

(5)

This means that the P-scheme leads to the same real solvation free energies as the M-scheme if the bulk and surface contributions are calculated with the same summation scheme consistently. That is, the use of real solvation free energies removes the differences in the results from P- and Mschemes.34 Both P- and M-summation schemes have been used in literature for the calculation of ionic solvation free energies. For instance, Roux40,85 and Palmer37 followed the P-summation scheme, while Hünenberger and co-workers favored the Mscheme and made great efforts to seek the M-scheme intrinsic free energies and surface potential.8,34,36 The ΔGDSC correction is used to convert P-scheme results to those of M-scheme and vice visa (eqs 3 and 4) and can be computed by36 ΔG DSC = −

qγs Ns 6ε0 L3

(6)

where γs is the quadrupole-moment trace of the solvent model (0.8476 e·Å2 for SPC/E), Ns is the number of solvent molecules, ε0 is the permittivity of vacuum, and L is the length of the simulated cubic box. Note that ΔGDSC depends on the choice of the origin of the water molecule and here the oxygen atom of water is chosen as the origin.35 The raw (uncorrected) TI calculations (ΔGraw) belong to the P-scheme intrinsic HFEs (ΔGintr,P), and addition of ΔGDSC (eq 3) leads to the M-scheme intrinsic HFEs (ΔGintr,M). Taken the solvent Galvani potential (ØG,P = −13.8 kcal/mol·e)86 into account for SPC/E (eq 5), the physical real HFEs are obtained. For the ERmod and 3D-RISM calculations, only the more meaningful real HFEs are examined in the following. For consistency with the modified water models in 3D-RISM,65 surface potential (ØG,P) of −13.43 kcal/mol·e and −12.55 kcal/ mol·e are used for cSPC/E and cTIP3P water models, respectively,37 while ERmod shares the same ØG,P as in TI calculations. Finite-size effects with periodic boundary conditions and potential inaccuracies in solvent model permittivity for TI and ERmod simulations are needed to be corrected as well. The semianalytical (ANA) formula proposed by Hünenberger et al.35,36 and reformulated (eq 3 in ref 39) in our previous work was used to compute finite-size effects (denoted as ΔGANA). Potential inaccuracies in the dielectric permittivity of solvent models (ΔGD) were estimated based on the difference in the solvation free energy between PB-based calculations with a dielectric constant of 78.46 (for water in the experiment) and with 62 (for the SPC/E model)87 by the AMBER 11 package.70 We note that with the setting in this work these two corrections are quite small, approximately 0.3 and −0.1 kcal/mol for ΔGANA and ΔGD, respectively, as given in Table S3 in the SI The free energy correction/addition terms above, namely, ΔGDSC, ØG,P, ΔGANA, and ΔGD, are only applicable to the simulations with a (semi)explicit description of solvent molecules, such as TI, ERmod, and 3D-RISM/PC+. Despite the existence of correction schemes for implicit solvent models, the results for SMD, PBSA, and GBSA were not corrected here, due to the highly approximate nature and nonadditivity of the correction terms.34 For clarity, the calculated HFEs with or D


Article


(CPA) of −265.9 kcal/mol by Tissandier et al.28 Some agree that the CPA estimation by Tissandier et al.28 contains a contribution from surface potential and is thus real.37,88,89 However, a validation of the CPA protocol indicated that CPA underestimated the surface contribution by 3.7 kcal/mol, leading to a corrected (real) ΔG*hyd(H+) of −262.2 kcal/ mol, 90 which is thereby consistent with independent predictions of −261.7 kcal/mol from literature average by Hünenberger and Reif34 and of −262.4 kcal/mol from experimental measurements of surface potential by Randles.19 The TI calculations for the 112 ions in this work (Figure 1) and for the charged amino acids in our previous work39 appear to support the finding above that the real and (M-scheme) intrinsic ΔG*hyd(H+) likely amount to around −262 and −265 kcal/mol, respectively. Here the authors choose a physical real * (H+) = −262.4 kcal/mol by Randles19 as a estimation of ΔGhyd reference for all the following comparisons with the calculated real HFEs (ΔGreal) from explicit and semiexplicit solvation models. Compared to the experiments, the TI calculations with GAFF/AM1-BCC for the 112 ions yield with both root-meansquare error (RMSE) and standard deviation (SD) being ∼7 kcal/mol, a mean absolute error (MAE) of ∼5 kcal/mol, and a correlation of R ∼ 0.9, as shown in Table 1 and Figure 2a. This accuracy (∼7 kcal/mol) is much larger than the experimental uncertainty (∼3 kcal/mol) of the aqueous solvation free energy for a typical ion,91 indicating a necessity in further improvement of the force field. ERmod shows a weaker correlation (R ∼ 0.8) and a larger deviation (i.e., a worse performance) from the experiments (RMSE ∼ 11 kcal/mol and both SD and MAE being ∼9 kcal/ mol) compared to TI (Table 1 and Figure 2b), revealing a discrepancy in TI and ERmod calculations. Comparison of the raw (uncorrected) ERmod and TI predictions is presented in Figure 3 (note that these two methods share identical free energy corrections as mentioned in section 2.8). For all the ions tested, ERmod correlates strongly to TI (R = 0.99) with an RMSE of ∼7 kcal/mol (Figure 3a). ERmod reproduces the TI predictions for the anions perfectly, as indicated by a linear fit of y = 0.99x (R = 0.99); for the cations, ERmod however tends to give a systematic underestimation by about 9 kcal/mol, with a linear fit of y = x + 9.33 (R = 0.99), as shown in Figure 3a. Similar observations are detected for the polar (Coulombic) component of the total hydration free energy (Figure 3b). For the nonpolar (van der Waals) part, ERmod does not reproduce the TI results with a low correlation (R = 0.77) and displays an overall overestimation (Figure 3c and Table S5 in the SI). 3.2. Semiexplicit Solvent. The 3D-RISM/PC+ method with a so-called semiexplicit description of solvent shows a similar performance to TI, as presented in Table 1 and Figure 2c. The raw (uncorrected) predictions by 3D-RISM/PC+ have a strong correlation with the TI predictions for HFEs of the 112 ions (R = 0.99, y = x + 0.56) with a RMSE of ∼3 kcal/mol and a high correlation (R ∼ 1, Figure 4a). Compared to the TI calculation, 3D-RISM/PC+ reproduces well the polar contribution to the ionic hydration (Figure 4b) but underestimates the nonpolar component systematically, as observed by Palmer et al. for the partial set of 70 ions,37 even though the nonpolar part by 3D-RISM/PC+ correlates strongly to that by TI (R = 0.95, Figure 4c and Table S5). Modeling the ions with GAFF/ RESP in TIP3P (Figure S2) gives a very similar performance to that with GAFF/AM1-BCC by 3D-RISM/PC+ in SPC/E (Table 1).

Figure 1. Correlation coefficient (R) and standard deviation (SD, kcal/mol) from experimental observations as a function of the * (H+) for P-scheme intrinsic (a), real (b), and Mreferenced ΔGhyd scheme intrinsic (c) hydration free energies from TI calculations of the 112 ions in the Minnesota solvation database. Refer to eqs 3 and 5 and Table S4 for the derivation of P-scheme intrinsic (ΔGintr,P), real (ΔGreal), and M-scheme intrinsic (ΔGintr,M) free energies of hydration. * (H+) of −265.9 Colored bars showcase the widely adopted ΔGhyd (yellow), −264.8 (magenta), −262.4 (cyan), and −254.3 (green) kcal/ mol. The solid blue line indicates a SD of 7 kcal/mol.

−262.4 kcal/mol,19 Hünenberger’s value of −264.8 kcal/mol,27 and Tissandier’s prediction of −265.9 kcal/mol,28 as indicated in Figures 1b and S1b. In addition, the comparison indicates that, among the four widely used experimental HFE sets, the Marcus’s set based on ΔGhyd * (H+) = −254.3 kcal/mol seems more appropriate for Pscheme intrinsic HFEs from PME/PBC MD simulations, as did by Merz and co-workers for rational design of metal models.20−23 However, P-scheme intrinsic results depend on the quadrupole-moment of the water model (eqs 3 and 6),35,39 which is essentially a random number. It can be tuned entirely arbitrarily in a water model by choosing an arbitrary origin of the solvent molecule without changing its properties. Thus, the P-scheme intrinsic results are entirely meaningless in the experiments. No experiment can probe this number, since it includes a component that depends on the construction of a solvent model without affecting any of its simulated properties. The dependence of P-scheme intrinsic free energies on the solvent model were able to be eliminated to a great extent via the use of real and M-scheme intrinsic results,35,39 which means that these two kinds of HFEs likely make sense in the ionic solvation experiment physically. Hü nenberger and Reif proposed an (M-scheme) intrinsic ΔGhyd * (H+) of −264.8 kcal/mol and a real one of −261.7 kcal/mol from literature average,27,34 and they also argued that the recommended intrinsic value is consistent with the cluster pair approximation E


Article

Journal of Chemical Information and Modeling Table 1. Comparison of Calculated Hydration Free Energies with Experimental Observationsa R

method TI ERmod 3D-RISM SMD

PB GB

AM1-BCC (SPC/E) AM1-BCC (SPC/E) AM1-BCC (SPC/E) RESP (TIP3P) B3LYP/6-31+G(d,p) HF/6-31+G(d,p) M06-2X/6-31+G(d,p) B3LYP/cc-pVTZ HF/cc-pVTZ M06-2X/cc-pVTZ SASA (γA)b CD (σA + c +ΔGdisp)c Still (tinker)d HCT (tinker)d HCT (mbondi)e OBCI (tinker)d OBCI (mbondi2)e OBCII (tinker)d OBCII (mbondi2)e

0.88 0.84 0.86 0.86 0.90 0.93 0.92 0.90 0.93 0.92 0.86 0.86 0.88 0.86 0.82 0.86 0.83 0.85 0.82

RMSE

(0.02) (0.03) (0.03) (0.03) (0.02) (0.02) (0.02) (0.02) (0.01) (0.01) (0.04) (0.04) (0.02) (0.03) (0.04) (0.03) (0.05) (0.03) (0.04)

7.0 10.7 7.0 7.0 8.0 5.6 7.0 7.9 5.8 6.9 7.3 7.4 6.7 7.9 9.6 8.0 8.7 8.5 9.0

(0.5) (0.5) (0.5) (0.5) (0.4) (0.4) (0.4) (0.4) (0.4) (0.4) (0.9) (0.9) (0.4) (0.5) (1.1) (0.5) (1.0) (0.5) (1.0)

SD 6.9 8.9 6.8 6.7 4.9 4.3 4.5 4.8 4.2 4.3 7.2 7.2 6.7 7.2 9.4 7.8 8.6 7.9 8.7

(0.6) (0.7) (0.5) (0.5) (0.2) (0.3) (0.2) (0.2) (0.2) (0.2) (1.0) (0.9) (0.5) (0.5) (0.5) (0.5) (0.5) (0.5) (0.5)

MAE 5.4 9.3 5.2 5.3 6.6 4.3 5.6 6.5 4.5 5.6 5.3 5.4 5.2 6.5 7.4 6.6 6.7 7.1 7.0

(0.4) (0.5) (0.5) (0.4) (0.4) (0.4) (0.4) (0.4) (0.4) (0.4) (0.5) (0.5) (0.4) (0.4) (0.6) (0.4) (0.5) (0.4) (0.5)

Experimental observations are based on the ΔG*hyd(H+) of −262.4 kcal/mol, expect for SMD on a value of −265.9 kcal/mol because SMD was designed to target this value.31 The real free energies of hydration were considered for TI, ERmod, and 3D-RISM, whereas the raw (uncorrected) results from SMD, PB, and GB are used. Correlation coefficients (R), root-mean-square error (RMSE, kcal/mol), standard deviation (SD, kcal/mol), and mean absolute error (MAE) were calculated by the bootstrapping with 1000 iterations using the R program,92 and standard deviations for these quantities are given in parentheses. bNonpolar hydration free energies amount to the solvent accessible surface area (A) multiplied by a surface tension (γ = 0.005 kcal/mol·Å2). cNonpolar free energies are decomposed into a cavity term and a dispersion term (σ = 0.0378 kcal/mol·Å2 and c = −0.5692). dvan der Waals radii for GB models were taken from the TINKER software.84 eGB radii from the AmberTools 15 package (mbondi for HCT and mbondi2 for the two OBC variants).62 a

quantum mechanical SMD calculations with six theoretical level/basis set combinations, namely, B3LYP/6-31+G(d,p), HF/6-31+G(d,p), M06-2X/6-31+G(d,p), B3LYP/cc-pVTZ, HF/cc-pVTZ, and M06-2X/cc-pVTZ, for the 112 ions best match the experiments based on a ΔG*hyd(H+) of −265.9 kcal/ mol (Table S6). This is not surprising, since SMD models were designed to reproduce this experimental reference.31 SMD gives high correlations (R > 0.9) with a RMSE of 6−8 kcal/mol, a SD of 4−5 kcal/mol, and a MAE of 4−7 kcal/mol, depending on the used levels of theory (Table 1). When matching with other experimental references, RMSE increases with the increasing ΔGhyd * (H+) (less negative) and even rises up to 17 kcal/mol for ΔG*hyd(H+) = −254.3 kcal/mol. For the tested 112 ions, the theoretical level of HF yields with a RMSE of 6 kcal/ mol and a SD of 4 kcal/mol and appears to perform better than the B3LYP and M06-2X levels of theory, while no significant differences are observed between the basis sets of 6-31+G(d,p) and cc-pVTZ (Table 1 and Figures 2d and S2). Unlike SMD, the best match for PB predictions is the experiments with ΔGhyd * (H+) = −262.4 kcal/mol and the RMSE from the experiments is ∼7 kcal/mol with a correlation of R = 0.86 (Tables 1 and S6 and Figures 2e and S2). Considering the difference in solvation treatments of explicit and implicit models, the comparison of polar (charging) contributions to the solvation between raw TI results and PB makes no sense and only the nonpolar contribution was thus compared to then. The nonpolar term based solely on the solvent accessible surface area (SASA), denoted as PB/SASA, has almost no correlation with that by TI (R = 0.15, Figure 5a and Table S5). However, the nonpolar part modeled by two terms of cavity and dispersion contributions, denoted as PB/ CD, yield a good prediction (R = 0.9) compared to TI (Figure 5b and Table S5).

Figure 2. Comparison of the calculated hydration free energies with experimental observations of the 112 ions for representative methods of TI (a), ERmod (b), 3D-RISM with AM1-BCC charges (c), SMD with HF/6-31+G(d,p) (d), PB/SASA (e), and GB(Still) (f). Correlation coefficients (R) and the linear fits are given for the assessment between different methods. Comparisons for the other 13 methods examined are given in Figure S2. Refer to Table 1 for the comparisons in detail. Cations are shown as red circles, and anions, as blue squares.

3.3. Implicit Solvent. A continuum medium is applied to treat solvent with SMD, PBSA, and GBSA methods. All the F


Article


Figure 4. Comparison of raw hydration free energies (a) and the components of Coulombic (b) and van der Waals (c) contributions from TI and 3D-RISM/PC+ predictions with GAFF/AM1-BCC in SPC/E water. The inset box in each panel shows root-mean-square errors (RMSE), standard deviations (SD), and correlation coefficients (R) of the linear fits. Cations are shown as red circles, and anions, as blue squares.

Figure 3. Comparison of raw hydration free energies (a) and the components of Coulombic (b) and van der Waals (c) contributions from TI and ERmod predictions with GAFF/AM1-BCC in SPC/E water. The inset box in each panel shows root-mean-square errors (RMSE), standard deviations (SD), and correlation coefficients (R) of the linear fits. Hydration free energies (a) of cationic (red) and anionic (blue) compounds were fitted separately to elucidate a potential deficiency for cations (underestimation by ∼9 kcal/mol). Cations are shown as red circles, and anions, as blue squares.

3.4. Systematic Assessment. Standard deviation (SD) from the experiments is used to assess the performance of computational methods examined in this work, since this quantity does not rely on whether the ex post free energy corrections are used or not. That is, with or without the corrections, SD remains unchanged for a specific ion set with the same charge sign for TI, ERmod, and 3D-RISM/PC+. Compared to the TI calculations, ERmod predicts the ionic hydration with an SD of ∼1 and ∼2 kcal/mol for cations and anions, respectively, and 3D-RISM/PC+ gives an SD of 2−3 kcal/mol (which is larger than ERmod) for the 112 ions, as shown in Figure 6a. This indicates that ERmod (with explicit solvent) have a potential capability of predicting ionic hydration with good accuracy, whereas it gives for now a worse performance due to the observed deficiency in modeling cations (Figure 3 and Table 1), compared to 3D-RISM/PC+ (with semiexplicit solvent). PB predictions have a SD of ∼5 kcal/mol and seem slightly worse than most of GB models with a SD of 3−4 kcal/mol. Among GB models, unexpectedly, the HCT model with the recommended “mbondi” radii shows a dramatic deficiency (SD ∼ 10 kcal/mol) for the cations (Figure 6a). When compared to the experimental observations (Figure 6b), TI shows a SD of ∼9 kcal/mol for the anions, which doubles that for the cations (SD ∼ 4 kcal/mol), indicating a biased deficiency of the GAFF force field in modeling ionic compounds with negative charges (Figures 2 and 6b). This finding is observed as well for ERmod, 3D-RISM/PC+, PB and GB predictions that are all based on the GAFF-modeled ions directly (Figure 6b). The standard protocol for GAFF (i.e., with

Similar to PB, the best matches for GB models are the experiments with ΔGhyd * (H+) = −262.4 kcal/mol (Table S6). The GB(Still) model predicts hydration free energies of the 112 ions with a RMSE of ∼7 kcal/mol (R = 0.88), which is better than HCT and two variants of OBC models (Tables 1 and S6). Compared to the radii from AMBER,62 GB predictions using the radii from TINKER84 appear to give a better performance, as indicated by a smaller RMSE/SD and a higher R (Tables 1 and S6). Comparison with the experiments for the tested four GB models with the radii from TINKER84 is presented in Figures 2f and S2. As expected, the polar components of GB reproduce the results by PB well (Figure 5c and d), since GB is inherent from PB electrostatics and designed attempting to approximate PB ones. GBSA implemented in GROMACS79 and PB/SASA in AMBER70 adopted for nonpolar approximations a similar protocol that is computed from the solvent accessible surface area multiplied by a surface tension of 0.005 kcal/mol·Å2. As a result, the nonpolar term by GBSA correlates to PB/SASA strongly (R = 0.89, the HCT model in Figure 5e) but does not show any correlation with PB/CD (R = 0.16, Figure 5f). The nonpolar contributions by the other GB models of Still and OBC correlate very weakly with TI and PB/CD, but have a strong correlation with PB/ SASA (Tables S5 and S7). Nevertheless, GB results still yield a total HFE that is comparable of both PB/SASA and PB/CD methods (R > 0.95), due to the small magnitude of nonpolar contributions (Table S7). G


Article


containing elements. Mean signed errors (MSE) from the experimental free energies of ionic hydration for these categories are given in Table 2. The TI calculations yield with a MSE of 1.4 kcal/mol for the 112 ions, which is comparable of 3D-RISM/PC+ (MSE ∼ 2 kcal/mol), PB (MSE ∼ 2 kcal/mol), and GB methods (MSE = 1−3 kcal/mol). TI underestimates the hydration of the anions with alcohol, alkene, and ketone groups by 4−8 kcal/mol and also underestimates the hydration of the cations with an aromatic group and some containing chlorine and sulfur elements by 3−4 kcal/mol. And, TI overestimates the hydration of carboxylic acid derivative anions and some containing fluorine, chlorine, bromine, and sulfur elements by 4−5 kcal/mol. For the anions (Table 2), ERmod gives a similar MSE to TI, whereas it displays a large positive MSE of 10−13 kcal/mol for the cations, as observed in Figure 3a for the deficiency in underestimating the cationic hydration. 3D-RISM/PC+ appears to underestimate the hydration of all the anions by 4−9 kcal/mol, expect for the carboxylic acid derivative anions (MSE = 1−2 kcal/mol), while it reproduces the cationic hydration well with an absolute MSE less than 3 kcal/mol. For a comprehensive assessment, the corresponding mean absolute errors (MAE, Table S8) and SDs (Table S9) are given in the Supporting Information. SMD implicit solvent models underestimate the ionic hydration systematically and the Hartree−Fock (HF) theory, with a small MSE of 4 kcal/mol, appears to outperform the theoretical levels of B3LYP and M06-2X (Table 2). Of the 112 ions, the 81 unclustered ions (Table S1) were included in the training set of SMD developments targeting the experiments based on ΔG*hyd(H+) = −265.9 kcal/mol,31 and SMD gives a better performance for these ions, as indicated by a smaller MSE of 2−5 kcal/mol, than that for the other 31 ions in the test set (MSE ∼ 7−11 kcal/mol, Table S10). Despite the unclustered ions, clustered ions with a single explicit water each were included in the training set as well during the SMD design.31 The cluster-continuum model such as SMD with explicit water molecules was shown to improve the SMD performance in HFE calculations, while the prediction was shown to depend on the level of theory and the number of explicit waters used.94 This features the complexity and difficulty in the ionic model development. PB tends to give an underestimation as well expect for overestimating carboxylic acid derivative anions and the anions containing fluorine, chlorine, bromine, and sulfur elements. GB

Figure 5. Comparison of nonpolar contributions to hydration free energies between PB and TI calculations (a and b) and of polar (c and d) and nonpolar (e and f) contributions between PB and GB(HCT). Two schemes for treating nonpolar contributions in PB predictions, namely, PB/SASA (left) and PB/CD (right), are presented. Standard deviations for correlation coefficient (R) of the linear fits (black lines) were calculated by the bootstrapping with 1000 iterations using the R program.92 Cations are shown as red circles, and anions, as blue squares.

RESP charges in TIP3P) appears not to resolve such biased deficiency for anions, as indicated by 3D-RISM/PC+ predictions with GAFF/RESP in TIP3P (Figure 6b). This bias is however not observed in quantum mechanical SMD calculations, and both cations and anions are predicted with a roughly equal SD of ∼4 kcal/mol from the experiments (Figure 6b). For a deeper insight into the performance of the 6 methods assessed here, the 112 ions were classified into different categories based on the functional groups and on the

Figure 6. Standard deviations of the calculated hydration free energies for the cations (red) and anions (blue) from TI (a) and experimental (b) observations. Seven approaches including 3D-RISM with RESP/TIP3P and quantum mechanical SMD models are not available (NA) to compare with TI (AM1-BCC/SPCE) directly. All the experiments are based on ΔG*hyd(H+) = −262.4 kcal/mol, except for SMD on the value of −265.9 kcal/ mol. Refer to Table 1 for detailed comparisons with the experiments. H


Article

Journal of Chemical Information and Modeling Table 2. Mean Signed Errors (kcal/mol) From Experimental Ionic Hydration Free Energiesa 3D-RISM/PC+ solute class

N

TI

alcohol anions alkene anions amine cations aromatic compound cations carboxylic acid dev. anions heterocyclic compound cations ketone anions H, C, N, O anions F, Cl, Br, S anions all anions H, C, N, O cations Cl, S cations all cations all ions

19 6 36 17 8 10 8 43 17 60 48 4 52 112

7.6 6.3 1.4 3.2 −4.0 1.3 3.7 2.1 −4.6 0.2 2.7 3.7 2.8 1.4

alcohol anions alkene anions amine cations aromatic compound cations carboxylic acid dev. anions heterocyclic compound cations ketone anions H, C, N, O anions F, Cl, Br, S anions all anions H, C, N, O cations Cl, S cations all cations all ions

19 6 36 17 8 10 8 43 17 60 48 4 52 112

ERmod 8.1 5.7 10.5 12.4 −4.1 9.7 3.1 2.4 −4.9 0.4 11.9 12.8 12.0 5.8 PB

SPC/E

TIP3P

9.3 8.7 −0.7 1.0 0.8 −1.5 6.5 4.2 −2.4 2.3 1.0 1.7 1.1 1.8

9.2 8.7 −0.1 1.6 2.3 −0.9 7.4 4.6 −3.1 2.4 1.7 2.5 1.8 2.1

B3LYP

HF

M06-2X

B3LYP

6-31+G(d,p) 11.1 10.3 1.3 2.7 5.4 0.4 9.0 9.2 9.7 9.3 2.8 4.6 2.9 6.4

5.2 5.3 0.2 0.6 0.4 0.0 4.3 4.3 6.5 4.9 2.0 2.8 2.0 3.6

HF

M06-2X

cc-pVTZ 8.9 8.8 0.9 1.9 3.6 0.2 7.5 7.4 8.4 7.7 2.5 4.5 2.6 5.3

11.8 10.7 1.8 3.6 5.2 0.8 9.2 8.7 8.7 8.7 3.4 5.4 3.6 6.3

6.8 6.6 0.7 1.6 1.0 0.4 5.4 4.9 6.4 5.3 2.5 3.7 2.6 4.0

9.5 8.9 1.6 3.0 3.3 0.8 7.3 7.0 7.7 7.2 3.1 5.0 3.3 5.4

GB

SASA

CD

Stillb

HCTb

HCTc

OBCIb

OBCIc

OBCIIb

OBCIIc

8.3 6.6 3.5 4.0 −1.0 2.3 4.5 3.3 −4.8 1.0 2.2 0.9 2.1 1.5

8.6 6.4 3.5 4.1 −1.2 1.6 4.4 3.4 −4.3 1.2 2.3 0.9 2.2 1.7

8.4 7.9 −0.2 0.4 −3.1 0.1 5.4 4.1 −4.3 1.7 −0.3 −0.9 −0.4 0.8

11.3 9.3 3.9 5.3 −1.9 3.1 6.3 4.9 −3.1 2.7 3.8 3.9 3.8 3.2

9.8 7.2 6.4 7.6 −2.8 4.3 4.2 3.0 −5.2 0.7 3.8 3.8 3.8 2.1

9.5 7.7 3.8 5.1 −4.2 3.1 4.6 3.0 −5.5 0.6 3.5 3.5 3.5 1.9

7.2 5.3 4.7 6.3 −5.1 2.2 2.2 0.6 −8.0 −1.8 4.6 4.4 4.6 1.2

10.9 8.8 5.0 6.5 −2.8 4.1 5.8 4.2 −4.7 1.7 4.8 5.0 4.8 3.1

8.2 6.2 5.3 7.1 −4.0 2.6 3.1 1.6 −7.4 −0.9 5.3 5.3 5.3 2.0

a

Real free energies of hydration are considered for TI, ERmod, and 3D-RISM methods, whereas the raw (uncorrected) results from SMD, PB, and GB are used. The comparison with the experiments follows Table 1. All the calculations were done with GAFF/AM1-BCC, expect for SMD with B3LYP, HF, and M06-2X levels of theory. TI, ERmod, and 3D-RISM used SPC/E water. GAFF/RESP/TIP3P was tested with 3D-RISM as well. Functional groups were assigned automatically by Checkmol (version 0.5a),93 and only those occurring at least 6 times were shown here. bvan der Waals radii for GB were taken from the TINKER software.84 cWith radii from the AmberTools 15 package62 (mbondi for HCT and mbondi2 for two OBC variants).

mol, which indicates that the deficiency in the force field of ionic compounds appears to limit the accuracy of the developed methods/models,37 highlighting the necessity of further model improvements of force fields, in particular for ionic compounds with negative charges. All the tested methods seem to generate a reasonable estimation for the polar (charging) contribution to the hydration compared to the TI results, expect for the biased underestimation of cations by ∼9 kcal/mol using ERmod, while the absolute predictions on the polar component differ significantly in explicit and implicit solvents. The commonly used strategy that correlates the total nonpolar contribution with solvent accessible surface areas (SASA), as did in PB/SA and GBSA methods, fails to handle the nonpolar predictions for the 112 ions. An alternative method that breaks down the nonpolar solvation contribution into two terms of cavity and dispersion interactions by Luo and co-workers,78 as did in the PB/CD scheme, reproduces the nonpolar predictions by TI results well. More efforts are required to improve the performance of nonpolar contribution predictions with the (semi)explicit models of ERmod and 3D-RISM.

models behave similarly to PB predictions; the GB(Still) model performs best giving a MSE within 1 kcal/mol on average, followed by OBCI, and the HCT model shows a similar performance to OBCII. The GB results with HCT and two OBC models using the radii from AMBER62 yield a smaller MSE than that from TINKER,84 whereas the latter radii produce a smaller RMSE and a higher correlation with the experiments (Tables 1 and S6).

4. CONCLUDING REMARKS For the large set of ionic compounds modeled by the generalized Amber force field (GAFF),43 the methods of TI, 3D-RISM, PBSA, and GBSA provide an HFE accuracy of ∼7 kcal/mol, a value much larger that the experimental uncertainty (∼3 kcal/mol) of the aqueous solvation free energy for a typical ion.91 The GAFF models appear to show a biased deficiency in the modeling of the anionic compounds, as indicated by a large SD of 8−10 kcal/mol from the experiments; for the cationic species, however, GAFF gives a small SD of 4−5 kcal/mol. Quantum mechanical SMD models have no such bias and produce a SD of ∼4 kcal/mol. Compared to TI, ERmod (for anions only), and 3D-RISM yield with a small SD of 2−3 kcal/ I


Article

Journal of Chemical Information and Modeling Experimental thermodynamics of ion hydration depend on the experimental HFE of the proton, ΔG*hyd(H+), which differs by over 10 kcal/mol. Obviously only one of these numbers is correct but cannot be measured directly, and the choice of reference value is crucial for correct calculation of free energy of hydration for ions. As a result, the debated ΔG*hyd(H+) leads to controversial experimental thermodynamics of ionic hydration. The simplicity of the models and the approximations used in the force field calculations necessitate in addition a number of corrections (denoted as ΔGcor). These are a bit confusing and hinder both methodology and model development efforts for ions in computational chemistry. Hydration free energy (HFE) calculations with thermodynamic integration (TI) for all the 112 singly charged ions modeled by the generalized Amber force field (GAFF)43 appear * (H+)to showcase empirical best matches between ΔGhyd dependent experiments and ΔGcor-dependent predictions on the thermodynamics of ionic hydration. That is, P-scheme intrinsic free energies of hydration from for instance PME/PBC calculations as typically done in the MD simulations seems more appropriate to target the Marcus’s experimental set based on ΔG*hyd(H+) = −254.3 kcal/mol.15 Taking the air/water surface potential (ØG) into account, the resulting real hydration free energies (HFE) are more consistent with the Randles’s experimental set based on ΔG*hyd(H+) = −262.4 kcal/mol.19 For the M-scheme intrinsic free energies, the experimental set based on the Tissandier’s prediction of −265.9 kcal/mol28 appears more appropriate; note that Hünenberger’s ΔGhyd * (H+) of −264.8 kcal/mol27 is argued to be consistent with that by Tissandier et al.28 These empirical matches are in line with our previous work on the HFE calculations of charged amino acid side chains with varied protein force fields,39 and allow, in a computational sense, consistent HFE calculations with the adopted ΔG*hyd(H+) reference using the state-of-art classical bimolecular force fields. It should be noted however that the best estimation of ΔG*hyd(H+) is still an ongoing debate. Although the GAFF and most of the examined protein force field parameters39 did not rely explicitly on the ΔGhyd * (H+)-based calibration strategy, there is no reason to believe that the state-of-art classical force fields is accurate enough to deduce the ΔG*hyd(H+) they imply. The modern polarizable potentials have presented a new level of insight for free energy calculations of ions in water.95,96 Further studies on the model improvements may yield the calculations with high accuracy from which a correct ΔG*hyd(H+) can be drawn with confidence.

■

Chunhua Yin: 0000-0002-5067-2951 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Prof. David van der Spoel for fruitful comments and Prof. Tianwei Tan for a grant of computer time through ChemCloudComputing of Beijing University of Chemical Technology. This work was supported by the National Natural Science Foundation of China (21606016, 21677011, 21576017, 21436002, and 11471034), Beijing Natural Science Foundation (5174036), and the Fundamental Research Funds for the Central Universities (FRF-TP-15-018A1).

■

(1) Cook, T. R.; Zheng, Y.-R.; Stang, P. J. Metal-Organic Frameworks and Self-Assembled Supramolecular Coordination Complexes: Comparing and Contrasting the Design, Synthesis and Functionality of Metal-Organic Materials. Chem. Rev. 2013, 113, 734. (2) Vallee, B. L.; Auld, D. S. Active-Site Zinc Ligands and Activated H2O of Zinc Enzymes. Proc. Natl. Acad. Sci. U. S. A. 1990, 87, 220− 224. (3) Costas, M.; Mehn, M. P.; Jensen, M. P.; Que, L. Dioxygen Activation at Mononuclear Nonheme Iron Active Sites: Enzymes, Models, and Intermediates. Chem. Rev. 2004, 104, 939−986. (4) Herschlag, D.; Jencks, W. P. The Effect of Divalent Metal Ions on the Rate and Transition-State Structure of Phosphoryl-Transfer Reactions. J. Am. Chem. Soc. 1987, 109, 4665−4674. (5) Schrag, J. D.; Li, Y.; Wu, S.; Cygler, M. Ser-His-Glu Triad Forms the Catalytic Site of the Lipase from Geotrichum Candidum. Nature 1991, 351, 761−764. (6) Li, P.; Merz, K. M., Jr Metal Ion Modeling Using Classical Mechanics. Chem. Rev. 2017, 117, 1564−1686. (7) 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. (8) Reif, M. M.; Hünenberger, P. H.; Oostenbrink, C. New Interaction Parameters for Charged Amino Acid Side Chains in the GROMOS Force Field. J. Chem. Theory Comput. 2012, 8, 3705−3723. (9) Shirts, M. R.; Pande, V. S. Solvation Free Energies of Amino Acid Side Chain Analogs for Common Molecular Mechanics Water Models. J. Chem. Phys. 2005, 122, 134508. (10) Shirts, M. R.; Pitera, J. W.; Swope, W. C.; Pande, V. S. Extremely Precise Free Energy Calculations of Amino Acid Side Chain Analogs: Comparison of Common Molecular Mechanics Force Fields for Proteins. J. Chem. Phys. 2003, 119, 5740−5761. (11) Mobley, D. L.; Guthrie, J. P. Freesolv: A Database of Experimental and Calculated Hydration Free Energies, with Input Files. J. Comput.-Aided Mol. Des. 2014, 28, 711−720. (12) Mobley, D. L.; Dill, K. A.; Chodera, J. D. Treating Entropy and Conformational Changes in Implicit Solvent Simulations of Small Molecules. J. Phys. Chem. B 2008, 112, 938−946. (13) Zhang, J.; Tuguldur, B.; van der Spoel, D. Force Field Benchmark of Organic Liquids. 2. Gibbs Energy of Solvation. J. Chem. Inf. Model. 2015, 55, 1192−1201. (14) Martins, S. A.; Sousa, S. F. Comparative Assessment of Computational Methods for the Determination of Solvation Free Energies in Alcohol-Based Molecules. J. Comput. Chem. 2013, 34, 1354−1362. (15) Marcus, Y. Thermodynamics of Solvation of Ions. Part 5. Gibbs Free Energy of Hydration at 298.15 K. J. Chem. Soc., Faraday Trans. 1991, 87, 2995−2999. (16) Noyes, R. M. Thermodynamics of Ion Hydration as a Measure of Effective Dielectric Properties of Water. J. Am. Chem. Soc. 1962, 84, 513−522.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.7b00485. Experimental and calculated (TI) thermodynamics of ionic hydration and supplementary tables and figures (PDF)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (H.Y.). *E-mail: [email protected] (C.Y.). ORCID

Haiyang Zhang: 0000-0002-2410-7078 Yang Jiang: 0000-0003-1100-9177 J


Article


(37) Misin, M.; Fedorov, M. V.; Palmer, D. S. Hydration Free Energies of Molecular Ions from Theory and Simulation. J. Phys. Chem. B 2016, 120, 975−983. (38) Aldeghi, M.; Heifetz, A.; Bodkin, M. J.; Knapp, S.; Biggin, P. C. Accurate Calculation of the Absolute Free Energy of Binding for Drug Molecules. Chem. Sci. 2016, 7, 207−218. (39) Zhang, H.; Jiang, Y.; Yan, H.; Yin, C.; Tan, T.; van der Spoel, D. Free-Energy Calculations of Ionic Hydration Consistent with the Experimental Hydration Free Energy of the Proton. J. Phys. Chem. Lett. 2017, 8, 2705−2712. (40) Simonson, T.; Hummer, G.; Roux, B. Equivalence of M- and PSummation in Calculations of Ionic Solvation Free Energies. J. Phys. Chem. A 2017, 121, 1525−1530. (41) Duignan, T. T.; Baer, M. D.; Schenter, G. K.; Mundy, C. J. Electrostatic Solvation Free Energies of Charged Hard Spheres Using Molecular Dynamics with Density Functional Theory Interactions. ArXiv.org 2017, No. arXiv:1702.05203. (42) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3, 300−313. (43) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (44) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method. J. Comput. Chem. 2000, 21, 132−146. (45) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: II. Parameterization and Validation. J. Comput. Chem. 2002, 23, 1623−1641. (46) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (47) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (48) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N.Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (49) Zhang, H.; Lv, Y.; Tan, T.; van der Spoel, D. Atomistic Simulation of Protein Encapsulation in Metal−Organic Frameworks. J. Phys. Chem. B 2016, 120, 477−484. (50) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (51) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (52) 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. (53) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (54) Nose, S.; Klein, M. L. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055−1076. (55) 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−3690. (56) Bennett, C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976, 22, 245−268. (57) Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energetic Representation. I. Formulation. J. Chem. Phys. 2000, 113, 6070−6081. (58) Sakuraba, S.; Matubayasi, N. ERmod: Fast and Versatile Computation Software for Solvation Free Energy with Approximate Theory of Solutions. J. Comput. Chem. 2014, 35, 1592−1608. (59) Misin, M.; Fedorov, M. V.; Palmer, D. S. Communication: Accurate Hydration Free Energies at a Wide Range of Temperatures from 3D-RISM. J. Chem. Phys. 2015, 142, 091105.

(17) Rosseinsky, D. Electrode Potentials and Hydration Energies. Theories and Correlations. Chem. Rev. 1965, 65, 467−490. (18) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous Solvation Free Energies of Ions and Ion−Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton. J. Phys. Chem. B 2006, 110, 16066−16081. (19) Randles, J. The Real Hydration Energies of Ions. Trans. Faraday Soc. 1956, 52, 1573−1581. (20) Li, P.; Roberts, B. P.; Chakravorty, D. K.; Merz, K. M., Jr Rational Design of Particle Mesh Ewald Compatible Lennard-Jones Parameters for + 2 Metal Cations in Explicit Solvent. J. Chem. Theory Comput. 2013, 9, 2733−2748. (21) Li, P.; Merz, K. M., Jr Taking into Account the Ion-Induced Dipole Interaction in the Nonbonded Model of Ions. J. Chem. Theory Comput. 2014, 10, 289−297. (22) Li, P.; Song, L. F.; Merz, K. M., Jr Parameterization of Highly Charged Metal Ions Using the 12−6-4 LJ-Type Nonbonded Model in Explicit Water. J. Phys. Chem. B 2015, 119, 883−895. (23) Li, P.; Song, L. F.; Merz, K. M., Jr Systematic Parameterization of Monovalent Ions Employing the Nonbonded Model. J. Chem. Theory Comput. 2015, 11, 1645−1657. (24) Åqvist, J. Ion-Water Interaction Potentials Derived from Free Energy Perturbation Simulations. J. Phys. Chem. 1990, 94, 8021−8024. (25) Duarte, F.; Bauer, P.; Barrozo, A.; Amrein, B. A.; Purg, M.; Åqvist, J.; Kamerlin, S. C. L. Force Field Independent Metal Parameters Using a Nonbonded Dummy Model. J. Phys. Chem. B 2014, 118, 4351−4362. (26) Liao, Q.; Kamerlin, S. C. L.; Strodel, B. Development and Application of a Nonbonded Cu2+ Model That Includes the Jahn− Teller Effect. J. Phys. Chem. Lett. 2015, 6, 2657−2662. (27) Reif, M. M.; Hünenberger, P. H. Computation of MethodologyIndependent Single-Ion Solvation Properties from Molecular Simulations. IV. Optimized Lennard-Jones Interaction Parameter Sets for the Alkali and Halide Ions in Water. J. Chem. Phys. 2011, 134, 144104. (28) Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle, T. R. The Proton’s Absolute Aqueous Enthalpy and Gibbs Free Energy of Solvation from Cluster-Ion Solvation Data. J. Phys. Chem. A 1998, 102, 7787−7794. (29) Camaioni, D. M.; Schwerdtfeger, C. A. Comment on “Accurate Experimental Values for the Free Energies of Hydration of H+, OH-, and H3O+”. J. Phys. Chem. A 2005, 109, 10795−10797. (30) Marenich, A.; Kelly, C.; Thompson, J.; Hawkins, G.; Chambers, C.; Giesen, D.; Winget, P.; Cramer, C.; Truhlar, D. Minnesota Solvation Database, Version 2012; University of Minnesota: Minneapolis, 2012. (31) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (32) Jiang, Y.; Zhang, H.; Feng, W.; Tan, T. Refined Dummy Atom Model of Mg2+ by Simple Parameter Screening Strategy with Revised Experimental Solvation Free Energy. J. Chem. Inf. Model. 2015, 55, 2575−2586. (33) Jiang, Y.; Zhang, H.; Tan, T. Rational Design of MethodologyIndependent Metal Parameters Using a Nonbonded Dummy Model. J. Chem. Theory Comput. 2016, 12, 3250−3260. (34) Hünenberger, P.; Reif, M. Single-Ion Solvation: Experimental and Theoretical Approaches to Elusive Thermodynamic Quantities; Royal Society of Chemistry: Cambridge, 2011. (35) Kastenholz, M. A.; Hünenberger, P. H. Computation of Methodology-Independent Ionic Solvation Free Energies from Molecular Simulations. II. The Hydration Free Energy of the Sodium Cation. J. Chem. Phys. 2006, 124, 224501. (36) Rocklin, G. J.; Mobley, D. L.; Dill, K. A.; Hünenberger, P. H. Calculating the Binding Free Energies of Charged Species Based on Explicit-Solvent Simulations Employing Lattice-Sum Methods: An Accurate Correction Scheme for Electrostatic Finite-Size Effects. J. Chem. Phys. 2013, 139, 184103. K


Article

Journal of Chemical Information and Modeling (60) Sergiievskyi, V.; Jeanmairet, G.; Levesque, M.; Borgis, D. Solvation Free-Energy Pressure Corrections in the Three Dimensional Reference Interaction Site Model. J. Chem. Phys. 2015, 143, 184116. (61) Ratkova, E. L.; Palmer, D. S.; Fedorov, M. V. Solvation Thermodynamics of Organic Molecules by the Molecular Integral Equation Theory: Approaching Chemical Accuracy. Chem. Rev. 2015, 115, 6312−6356. (62) Case, D. A.; Berryman, J. T.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Luchko, T.; Luo, R.; Madej, B.; Merz, K. M.; Monard, G.; Needham, P.; Nguyen, H.; Nguyen, H. T.; Omelyan, I.; Onufriev, A.; Roe, D. R.; Roitberg, A.; Salomon-Ferrer, R.; Simmerling, C. L.; Smith, W.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; York, D. M.; Kollman, P. A. AmberTools 15; University of California: San Francisco, 2015. (63) Perkyns, J.; Pettitt, B. M. A Site−Site Theory for Finite Concentration Saline Solutions. J. Chem. Phys. 1992, 97, 7656−7666. (64) Kovalenko, A.; Hirata, F. Self-Consistent Description of a Metal−Water Interface by the Kohn−Sham Density Functional Theory and the Three-Dimensional Reference Interaction Site Model. J. Chem. Phys. 1999, 110, 10095−10112. (65) Luchko, T.; Gusarov, S.; Roe, D. R.; Simmerling, C.; Case, D. A.; Tuszynski, J.; Kovalenko, A. Three-Dimensional Molecular Theory of Solvation Coupled with Molecular Dynamics in Amber. J. Chem. Theory Comput. 2010, 6, 607−624. (66) Kovalenko, A.; Hirata, F. Potentials of Mean Force of Simple Ions in Ambient Aqueous Solution. I. Three-Dimensional Reference Interaction Site Model Approach. J. Chem. Phys. 2000, 112, 10391− 10402. (67) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The Resp Model. J. Phys. Chem. 1993, 97, 10269−10280. (68) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision B.01; Gaussian, Inc.: Wallingford, CT, 2009. (69) Hehre, W. J.; Ditchfield, R.; Pople, J. A. 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−2261. (70) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 11; University of California: San Francisco, CA, 2010. (71) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (72) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098.

(73) 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. (74) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (75) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (76) Hehre, W.; Radom, L.; Schleyer, P. v. R.; Pople, J. ab initio Molecular Orbital Theory; Wiley: New York, 1986. (77) Dunning, T. H., Jr Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (78) Tan, C.; Tan, Y.-H.; Luo, R. Implicit Nonpolar Solvent Models. J. Phys. Chem. B 2007, 111, 12263−12274. (79) Larsson, P.; Lindahl, E. A High-Performance ParallelGeneralized Born Implementation Enabled by Tabulated Interaction Rescaling. J. Comput. Chem. 2010, 31, 2593−2600. (80) Qiu, D.; Shenkin, P. S.; Hollinger, F. P.; Still, W. C. The GB/SA Continuum Model for Solvation. A Fast Analytical Method for the Calculation of Approximate Born Radii. J. Phys. Chem. A 1997, 101, 3005−3014. (81) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Parametrized Models of Aqueous Free Energies of Solvation Based on Pairwise Descreening of Solute Atomic Charges from a Dielectric Medium. J. Phys. Chem. 1996, 100, 19824−19839. (82) Onufriev, A.; Bashford, D.; Case, D. A. Exploring Protein Native States and Large-Scale Conformational Changes with a Modified Generalized Born Model. Proteins: Struct., Funct., Genet. 2004, 55, 383−394. (83) Schaefer, M.; Bartels, C.; Karplus, M. Solution Conformations and Thermodynamics of Structured Peptides: Molecular Dynamics Simulation with an Implicit Solvation Model. J. Mol. Biol. 1998, 284, 835−848. (84) Ponder, J. W.; Richards, F. M. An Efficient Newton-Like Method for Molecular Mechanics Energy Minimization of Large Molecules. J. Comput. Chem. 1987, 8, 1016−1024. (85) Lamoureux, G.; Roux, B. Absolute Hydration Free Energy Scale for Alkali and Halide Ions Established from Simulations with a Polarizable Force Field. J. Phys. Chem. B 2006, 110, 3308−3322. (86) Arslanargin, A.; Beck, T. L. Free Energy Partitioning Analysis of the Driving Forces That Determine Ion Density Profiles near the Water Liquid-Vapor Interface. J. Chem. Phys. 2012, 136, 104503. (87) van der Spoel, D.; van Maaren, P. J. The Origin of Layer Structure Artifacts in Simulations of Liquid Water. J. Chem. Theory Comput. 2006, 2, 1−11. (88) Beck, T. L. The Influence of Water Interfacial Potentials on Ion Hydration in Bulk Water and near Interfaces. Chem. Phys. Lett. 2013, 561, 1−13. (89) Asthagiri, D.; Pratt, L. R.; Ashbaugh, H. Absolute Hydration Free Energies of Ions, Ion−Water Clusters, and Quasichemical Theory. J. Chem. Phys. 2003, 119, 2702−2708. (90) Vlcek, L.; Chialvo, A. A.; Simonson, J. M. Correspondence between Cluster-Ion and Bulk Solution Thermodynamic Properties: On the Validity of the Cluster-Pair-Based Approximation. J. Phys. Chem. A 2013, 117, 11328−11338. (91) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. SM6: A Density Functional Theory Continuum Solvation Model for Calculating Aqueous Solvation Free Energies of Neutrals, Ions, and SoluteWater Clusters. J. Chem. Theory Comput. 2005, 1, 1133−1152. (92) R Core Team. R: A Language and Environment for Statistical Computing, version 3.2.3; R Foundation for Statistical Computing: Vienna, Austria, 2015. L


Article

Journal of Chemical Information and Modeling (93) Haider, N. Functionality Pattern Matching as an Efficient Complementary Structure/Reaction Search Tool: An Open-Source Approach. Molecules 2010, 15, 5079−5092. (94) Riccardi, D.; Guo, H.-B.; Parks, J. M.; Gu, B.; Liang, L.; Smith, J. C. Cluster-Continuum Calculations of Hydration Free Energies of Anions and Group 12 Divalent Cations. J. Chem. Theory Comput. 2013, 9, 555−569. (95) Caleman, C.; Hub, J. S.; van Maaren, P. J.; van der Spoel, D. Atomistic Simulation of Ion Solvation in Water Explains Surface Preference of Halides. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 6838− 6842. (96) Hub, J. S.; Wolf, M. G.; Caleman, C.; van Maaren, P. J.; Groenhof, G.; van der Spoel, D. Thermodynamics of Hydronium and Hydroxide Surface Solvation. Chem. Sci. 2014, 5, 1745−1749.

M


Comparative Assessment of Computational ... - ACS Publications

Recommend Documents