Comparison of Implicit and Explicit Solvent Models for the Calculation

Article pubs.acs.org/JCTC

Comparison of Implicit and Explicit Solvent Models for the Calculation of Solvation Free Energy in Organic Solvents Jin Zhang,† Haiyang Zhang,‡ Tao Wu,† Qi Wang,*,† and David van der Spoel*,§ †

Department of Chemistry and Soft Matter Research Center, Zhejiang University, Hangzhou 310027, China Department of Biological Science and Engineering, School of Chemistry and Biological Engineering, University of Science and Technology Beijing, 100083 Beijing, China § Uppsala Center for Computational Chemistry, Science for Life Laboratory, Department of Cell and Molecular Biology, Uppsala University, Husargatan 3, Box 596, SE-75124 Uppsala, Sweden ‡

S Supporting Information *

ABSTRACT: Quantitative prediction of physical properties of liquids is important for many applications. Computational methods based on either explicit or implicit solvent models can be used to approximate thermodynamics properties of liquids. Here, we evaluate the predictive power of implicit solvent models for solvation free energy of organic molecules in organic solvents. We compared the results calculated with four generalized Born (GB) models (GBStill, GBHCT, GBOBCI, and GBOBCII), the Poisson−Boltzmann (PB) model, and the density-based solvent model SMD with previous solvation free energy calculations (Zhang et al. J. Chem. Inf. Model. 2015, 55, 1192−1201) and experimental data. The comparison indicates that both PB and GB give poor agreement with explicit solvent calculations and even worse agreement with experiments (rootmean-square deviation ≈ 15 kJ/mol). The main problem seems to be the prediction of the apolar contribution, which should include the solvent entropy. The quantum mechanical-based SMD model gives significantly better agreement with experimental data than do PB or GB, but it is not as good as explicit solvent calculation results. The dielectric constant ε of the solvent is found to be a powerful predictor for the polar contribution to the free energy in implicit models; however, the Onsager relation may not hold for realistic solvent, as suggested by explicit solvent and SMD calculations. From the comparison, we also find that with an optimization of the apolar contribution, the PB model gives slightly better agreement with experiments than the SMD model, whereas the correlation between the optimized GB models and experiments remains poor. Further optimization of the apolar contribution is needed for GB models to be able to treat solvents other than water.

1. INTRODUCTION Quantitative prediction of thermodynamics properties of solute molecules requires an accurate description of the solvent. Distinct solvent models may refer to either explicit solvent molecules or an implicit (analytical) description of the solvent environment. The evaluation of properties like the density, enthalpy of vaporization, heat capacity, surface tension, compressibility, expansion coefficient, and dielectric constant of pure liquids1,2 and solvation free energy (ΔGsolv) of organic molecules in organic solvents3−6 with explicit solvent models provided insights into the strengths and weakness of these organic solvent models. Following efforts to introduce longrange Lennard-Jones (LJ) interactions by implementing LJPME (particle mesh Ewald),7,8 the test set taken from the work of Caleman et al.1 was re-evaluated,9 leading to an increase in the predicted surface tension value of liquids, which was previously found to be systematically underestimated with truncated LJ interactions.1,10 The introduction of LJ-PME improved the reproduction of some liquid properties when compared to experiment, but it also highlighted a potential systematic error in force fields, an overestimated dispersion © 2017 American Chemical Society

coefficient (C6 in the Lennard-Jones potential) that is yet to be addressed. Although the accuracy of force fields is being improved systematically,11 the relatively high computational cost of explicit solvent models has resulted in implicit solvent models remaining popular.12 In most implicit solvent models, the solvent is treated as a structureless continuum with certain dielectric and interfacial properties.13−16 As a result, the number of interacting particles and the number of degrees of freedom of a system are significantly reduced, considering that explicit solvent molecules can contribute over 90% of atoms in a simulated system. Leaving out the explicit atomistic description of solvent comes at the cost of lacking hydrogen bonds with solvent, overstabilized salt bridges and hydrogen bonds within the solute, incorrect ion distribution,17 and unphysical sampling.18,19 On the other hand, the conformational search of the solute is faster due to the absence of viscosity.12,18 This is practical for timeconsuming studies like protein folding simulations,20 predicting Received: February 16, 2017 Published: February 28, 2017 1034

DOI: 10.1021/acs.jctc.7b00169 J. Chem. Theory Comput. 2017, 13, 1034−1043

Article

Journal of Chemical Theory and Computation

Because of the inconsistent results obtained using implicit solvent models, it has been suggested that more systematic benchmarks may be instructive.57 Most available benchmarks focus on hydration free energy of drug-like ligands.36,58,59 In order to faithfully describe the partition properties between different media, which is important in the association of biomolecules (drug design), protein folding, formation of micelles and membranes, protein−membrane interactions, and so forth, the accuracy of solvation free energy in solvents other than water should be considered as well.5,6 The emergence of implicit membrane models60−63 and growing interests in membrane protein simulation further highlight the importance of considering solvation free energy in more general solvents than water. A recent paper compared solvation models used to describe the dimerization free energy of β-cyclodextrins (β-CD) in varying organic solvents.19 Large discrepancies were found because hydrogen bonds between solvent and solute molecules were not explicitly taken into consideration, which resulted in strengthened H-bonds between solute molecules and hence biased sampling of solutes. Interestingly, the dimerization free energies calculated with implicit solvents were found to correlate strongly with 2(ε − 1)/(2ε + 1), obeying the Onsager relation,64,65 where ε is the dielectric constant of the solvent, whereas those calculated with explicit solvents did not.19 As suggested in that work, the experimental data ΔGsolv of our tested systems5,6 correlated poorly with 2(ε − 1)/(2ε + 1) as well.19 It should be noted that the Onsager relation was originally used to describe the relationship between ε and the change in the electrostatic potential, and the dimerization free energy clearly depends on more than just electrostatics. To clarify this issue and to evaluate ΔGsolv as predicted by implicit solvent models, we applied GB/PB models to systems employed in our previous benchmark of solvation free energy in organic solvents. The results are compared to calculations with explicit solvents and experimental data,5,6 as well as to a quantum chemistry-based implicit solvent model, SMD,66 which was developed as a universal solvent model. The GB models considered here include Still,25 HCT,26 and OBC (OBC-I and OBC-II).28

pH-dependent conformational changes, or in binding free energy calculations for drug design.12 All such applications rely on implicit solvent models being capable of accurately calculating ΔGsolv. Implicit solvent models typically uncouple polar (or electrostatic) interactions and nonpolar (or apolar) interactions (which is a crude approximation in itself) ΔGsolv = ΔGpolar + ΔGapolar

(1)

Because of the important role of the polar contribution (in high dielectric solvents such as water), electrostatic interaction theory is well-established. The polar contribution is generally approximated through solving the Poisson−Boltzmann (PB) equation21−23 or further simplified by using generalized Born (GB) model.24−28 The GB model was established to approach the results of PB with lower computational cost, and it has been claimed that, with “perfect” atomic radii, the GB model should be sufficiently accurate to reproduce the PB model.29 Indeed, development of GB models is still an active area.30 On the other hand, the nonpolar contribution is most often estimated by the solvent-accessible surface area (SASA), an argument loosely based on scaled particle theory (SPT)31,32

ΔGapolar = γA + c

(2)

where γ is a surface tension parameter obtained by fitting experimental hydration free energies of saturated linear hydrocarbons, A refers to SASA, and c is a constant, corresponding somehow to the free energy of the vacuum. Poor correlation has been found between nonpolar hydration free energies calculated through this approach and those using explicit solvent models. This finding applies to general organic molecules,33−36 longer alkanes (e.g., C100),37 and complex systems like protein−protein complexes,38−40 for which the nonpolar contribution is oversimplified. Moreover, the reported values of γ vary by over 1 order of magnitude,41−43 depending on the definition of solute surface area, the training set, and the underlying force fields used to describe the solute molecule. Work by Levy et al. highlighted the importance of the attractive van der Waals term.34,38,44,45 Wagoner et al. suggested that considering the solvent-accessible volume (SAV) and the attractive van der Waals interactions between solvent and solute in addition to SASA provides a more comprehensive picture of nonpolar interactions.45 Despite this, the SASA approach remains the main approach of estimating ΔGapolar in many applications such as high-throughput ligand screening.46,47 However, the GB/SASA approach may be unfit for ligand screening since it was found that GB electrostatics overestimates short-range interactions while simultaneously underestimating long-range screening effects.48 The SMD method used in quantum-chemical studies15 is somewhat more refined in the sense that there are separate parametrizations for different solvents described by a number of parameters (e.g., the Gaussian software49 supports over 180 different solvents). This is similar to the methodology used in the conductor-like screening models for realistic solvation (COSMO-RS) model,50−52 which provides good estimates for solvation free energies for the systems tested here.5,6 In contrast, the methods used in force field-based implicit solvents typically employ a distance-dependent dielectric constant for the electrostatic part53−55 combined with a set of γ values developed for water (except for a study by Bordner et al.56 that parametrized n-octanol and n-hexadecane solvent in GB).

2. METHODS 2.1. Test Compound Setup. Solvation free energies were calculated for 228 combinations of neutral organic solvents and solutes (molecule names, formulas, and CAS numbers can be found in Table S1), and both solvents and solutes were taken from our previous work.5,6 Organic solvents with experimental dielectric constants (ε) ranging from 2.4 (toluene) to 181.6 (Nmethylformamide) were examined. The distribution of ε is shown in Figure 1A, and ε for each solvent is given in Table S2. As shown in Figure 1B, both polar and nonpolar interactions are important in our test set. The electrostatics component (ΔGelec) and ΔGsolv were obtained from thermodynamic integration (TI) calculations with explicit solvents.5,6 The electrostatic contribution is defined here as the free energy change from state A, where the solute−solvent Coulombic interaction is turned off, to state B, where the solute fully interacts with solvent molecules. We note that this definition is well-defined technically but that the free energy at this intermediate point between states A and B does not have a direct physical meaning. 2.2. Conformation Sampling for PB and GB Calculations. It is common practice to calculate the solvation free energy using a single conformation of the solute when using 1035


Article


vacuum simulations in both GB (the ε value used for each solvent is shown in Table S2) and vacuum, and the potential energy differences (including polar and apolar components) needed for BAR were obtained to evaluate ΔGsolv. ΔGpolar was calculated in a similar way but without nonpolar calculations. For PB, ΔGpolar was calculated by solving the linearized form of the PB equation. ε = 2 was used for the solute region, and ε for solvent was set to the experimental value given in Table S2. The grid spacing was set to be 0.0208 nm. In APBS, ΔGapolar of a solute in conformation x ⃗ is calculated according to45 ΔGapolar(x ⃗) = γA(x ⃗) + pV (x ⃗)

Figure 1. (A) Dielectric constants of solvents studied. (B) Distribution of the ratio of ΔGelec and ΔGsolv from TI calculations.

+ ρ̅

∑∫ i

implicit methods. This assumes that the conformation distribution of the solute is similar in vacuum and solvated states, thus ignoring the conformational entropic and enthalpic changes of solute58 and the responses of the solvent correlated to this. Here, solute conformations of solvated and vacuum states were sampled from explicit solvent simulation and vacuum simulation at 300 K using the GROMACS software suite (version 4.5.5).67−71 The force field models of solute and solvent molecules were taken from the Web server at http:// virtualchemistry.org.1,9,72,73 The initial coordinates and velocities of the solvated state of each system were taken from the last frame of our previous simulation at λ0, where the solute molecule fully interacts with solvent.5,6 The Nosé−Hoover thermostat74,75 was used with a coupling constant of 0.5 ps for temperature coupling. Vacuum simulations were performed without periodic boundary conditions. The coordinates and velocities of solutes were taken from our previous simulations at λ1,5,6 where the solute does not interact with solvent. The vacuum simulations were performed without cutoffs for electrostatics and van der Waals interactions, and no pressure scaling was applied. The bond length of solute molecules was constrained with the LINCS algorithm.76 Each of these simulations were 5 ns long using an integration time step of 2 fs, and the snapshots were saved every 2 ps (2501 snapshots) for reprocessing using PB or GB algorithms. 2.3. Computation of the Free Energy Difference. For PB and GB, the Bennett acceptance ratio (BAR)77 method was used to estimate the free energy following a paper by Mobley et al.58 Alchemical free energy calculations in explicit solvent typically require multiple intermediate alchemical states to transition from the initial to the final state, and the potential energies from one state have to be evaluated in the other state. Mobley et al. found that the overlap between solvated and vacuum ensembles of the solute was sufficient to obtain reasonable statistical uncertainties 58 in implicit solvent calculations. Thus, no intermediate states were added here. We reprocessed the stored conformations with three types of GB model that are implemented in GROMACS,18 namely, Still,25 HCT,26 and OBC (OBC-I and OBC-II),28 and with the PB model using the adaptive Poisson−Boltzmann solver (APBS) package, version 1.3.23,45 For GB, the simulation setup is consistent with the vacuum simulation. The polar contribution in eq 1 was derived from the GB equation, and the apolar contribution was calculated from the Born radius of each atom using the ACE type approximation.78 This requires only one nonpolar solvation constant regardless of atom type, but it differs between GB models because of the different Born radii among these GB models.18 We reprocessed our solvated simulations5 and in

Ω

uiatt(x ⃗ , y ⃗ ) θ(xi⃗ , y ⃗ ) dy ⃗

(3)

where γ is a solvent surface tension parameter, A is SASA,79 p is a solvent pressure parameter, V is SAV,80 ρ is the bulk solvent number density, Ω represents the solvent accessible region outside the solute, uatt i (x ⃗ , y ⃗ ) is the attractive van der Waals potential in a Weeks−Chandler−Andersen (WCA) decomposition scheme81 for atom i at position y ⃗ , and θ(x ⃗ , y ⃗ ) is defined as the product of per-atom characteristic functions θ(x ⃗ , y ⃗ ) = ∏i θi(x ⃗ , y ⃗ ). For a particular solvent, γ was calculated by applying γ0·γsolv/γw, where γ0 is the surface tension parameter (0.0085 kJ mol−1 Å−2)45 optimized for a water probe radius of 0.65 Å, γsolv and γw are the experimental surface tensions of the solvent (given in Table S2) and water, respectively, and p was calculated by applying p0·pw/psolv, where p0 is the pressure parameter (0.2394 kJ mol−1 Å−3) optimized for the water probe radius mentioned above and psolv and pw are the experimental isothermal compressibilities of the solvent (given in Table S2) and water, respectively. As discussed below, we replaced ρ by the product of ρ and the number of heavy atom per solvent molecule to calculate ΔGapolar in our tests. The Shrake−Rupley numerical approximation was used to determine SASA.82 The SMD model66 was also considered in this work using the approach by Martins and Sousa.83 Calculations were performed with the Gaussian 09 suite of programs,49 and the default parameters were used. Following a gas-phase optimization at the MP2/6-31+G(d,p) level of theory, the solvation free energy was calculated from the single-point free energy difference between the gas phase and the liquid phase at four different levels of theory: HF/6-31+G(d,p),84 B3LYP85−87/6-31+G(d,p), M06-2X88/6-31+G(d,p), and M06-2X/cc-pVTZ.89 2.4. Statistics. For PB and GB calculations, the statistical uncertainties of implicit solvent calculations reported here were estimated from BAR (not shown in related plots for clarity). Bootstrapping analyses with 5000 iterations were applied for uncertainty estimation of the root-mean-square deviation (RMSD), mean signed error (MSE), and linear regression parameters: correlation coefficient (R), slope, and intercept. In the linear regression analysis, the uncertainties of implicit solvent calculations were ignored since they are considerablely smaller (mostly below 0.10 kJ/mol) than those of explicit solvent calculations. When comparing implicit and explicit solvent calculations, the data from explicit solvent calculations were taken as x and weighted by their uncertainties. To estimate ΔGpolar per solvent and solute, the Boltzmannenhanced discrimination of receiver-operating characteristic (BEDROC) metric90 was computed using the CROC package91 with 400 iterations in the R package.92 The BEDROC 1036


Article


(Table 1) than either GB or PB. PB calculations give significantly better agreement with both ΔGExpl and experiments (ΔGExper) than do GB calculations. However, plots of ΔGImpl versus ΔGExpl and ΔGImpl versus ΔGExper reveal rather weak correlation between these quantities. The correlation between experiments and ΔGImpl is weak for PB and GB but quite good for SMD models (Figure 3).

value is higher for solvents/solutes that appear often with high unsigned errors than for those with a uniform distribution (BEDROCunif). Student’s t test was also performed to compare the unsigned error for solvents/solutes with that for the entire test set. This helps to show whether the performence of a solvent/solute is significantly different from the whole set and is a supplement to the BEDROC approach.5,6,35

3. RESULTS 3.1. Overall Correlation. The correlation between calculated solvation free energies obtained using an explicit solvent model (ΔGExpl)5,6 and that obtained using an implicit solvent model (ΔGImpl) is given in Figure 2. Among these

Figure 3. Correlation between experimental solvation free energies (ΔGExper) and those computed from an explicit solvent model (TI, blue), PB (black), GB (red), and SMD (green) models (ΔGImpl): (A) PB, Still, and B3LYP/6-31+G(d,p), (B) HCT and HF/6-31+G(d,p), (C) OBC-I and M06-2X/6-31+G(d,p), and (D) OBC-II and M062X/cc-pVTZ. Figure 2. Correlation between calculated solvation free energies obtained using an explicit solvent model (ΔGExpl) and that obtained using PB (black), GB (red), and SMD (green) models (ΔGImpl): (A) Still and B3LYP/6-31+G(d,p), (B) HCT and HF/6-31+G(d,p), (C) OBC-I and M06-2X/6-31+G(d,p), and (D) OBC-II and M06-2X/ccpVTZ.

A breakdown of the RMSD and MSE per solvent/solute is given in Tables S3−S20 for the implicit solvent models tested. The total RMSD and MSE of all data points is given at the bottom of each related table (note that water as a solute is removed from this and the following analysis due to the discrepancy being over 400 kJ/mol, which turned out to be the result of using the rigid TIP3P model;93 the removal does not influence the main results). We find that the RMSD from

implicit solvent models, SMD at the levels of theory tested gives significantly better agreement with experimental data

Table 1. Correlation (R) and Regression to a Straight Line y = mx + b of Implicit Solvation Free Energies with Explicit Ones (ΔGExpl) and with Experimental Numbers (ΔGExper)a ΔGExpl R PB Still HCT OBC-I OBC-II HF/6-31+G(d,p) B3LYP/6-31+G(d,p) M06-2X/6-31+G(d,p) M06-2X/cc-pVTZ explicit a

0.69 0.45 0.45 0.44 0.45 0.73 0.68 0.76 0.73

± ± ± ± ± ± ± ± ±

ΔGExper

m 0.04 0.06 0.07 0.07 0.07 0.04 0.05 0.04 0.04

0.52 0.25 0.26 0.24 0.26 0.63 0.66 0.71 0.71

± ± ± ± ± ± ± ± ±

b 0.04 0.03 0.04 0.04 0.04 0.05 0.05 0.05 0.06

−2.9 −16.2 −16.5 −16.7 −16.5 −5.0 −5.8 −3.7 −5.0

R ± ± ± ± ± ± ± ± ±

1.4 0.6 0.6 0.6 0.6 1.1 1.1 1.1 1.2

0.51 0.23 0.24 0.24 0.25 0.76 0.74 0.77 0.77 0.84

± ± ± ± ± ± ± ± ±

m 0.05 0.06 0.06 0.06 0.06 0.03 0.03 0.03 0.03

0.84 0.53 0.51 0.54 0.53 1.09 0.93 0.96 0.89 1.04

± ± ± ± ± ± ± ± ±

b 0.08 0.15 0.13 0.14 0.14 0.09 0.07 0.07 0.05

−15.8 −3.4 −1.9 −1.5 −1.6 −4.2 −4.5 −5.4 −4.8

± ± ± ± ± ± ± ± ±

1.6 2.9 2.7 2.9 2.8 1.8 1.4 1.4 1.1

Explicit data are from refs 5 and 6. 1037


Article

Journal of Chemical Theory and Computation calculations with explicit solvent models lies between 4.8 ± 0.2 kJ/mol (M06-2X/cc-pVTZ) and 14.0 ± 0.5 kJ/mol (OBC-I) and the RMSD from experiments lies between 4.8 ± 0.2 kJ/mol (M06-2X/cc-pVTZ) and 14.8 ± 0.5 kJ/mol (PB). Overall, for calculations with PB and GB, the RMSD is significantly larger than the RMSD from experiments for calculations with an explicit solvent model in our previous work.5,6 For SMD, the different levels of theory sorted in order of increasing accuracy are HF/6-31+G(d,p), M06-2X/6-31+G(d,p), and B3LYP/631+G(d,p), with the most accurate one being M06-2X/ccpVTZ (Tables S8−S11). 3.2. Evaluation of ΔGpolar. For a better understanding of the discrepancy between results calculated with implicit and explicit solvent models, we further decompose ΔGsolv into ΔGpolar and ΔGapolar. As has been pointed out by Marenich et al., continuum solvation models differ in the way the polar contribution is defined;66 hence, a direct comparison of the individual part calculated by SMD with that from an explicit solvent model may not make sense. Hence, for SMD, we examined whether Onsager’s relation holds. The correlation between ΔGpolar calculated with explicit solvent (ΔGpolar‑Expl) and with implicit solvent (ΔGpolar‑Impl) is plotted in Figure 4. The correlation is strong (Table 2), and a

S33. Compared with GB calculations, PB has a stronger correlation for all solvents except for chloroform and has a larger slope (closer to 1) of the linear fits of ΔGpolar‑Expl with respect to ΔGpolar‑Impl. For PB, 1-octanol has the largest RMSD at 6.2 ± 0.7 kJ/mol and MSE at −5.4 ± 0.6 kJ/mol, followed by acetone (RMSD = 6.0 ± 1.1 kJ/mol and MSE at −4.9 ± 1.0 kJ/ mol), among the eight solvents (Table S3). The BEDROC values of 1-octanol and acetone are slightly larger than BEDROCunif (Table S23). In addition, Student’s t test (Table S27) suggests that 1-octanol is systematically off but acetone is not. For GB, 1-octanol and acetone also have larger RMSD and MSE values than the other six solvents (Tables S4−S7), and the BEDROC value is slightly larger than BEDROCunif (Table S23). However, Student’s t test suggests that the difference is not significant (Table S27). ΔGpolar can be regarded as the change in the electrostatic energy involved in transferring a molecule from vacuum to media with a dielectric constant ε, which is the only parameter that is used to represent the electrostatic response due to solvent in PB/GB. These implicit solvent models have been developed on the basis of early work by Born,94 Bell,64 Kirkwood,95 and Onsager.65 For a dipole molecule, ΔGpolar can be calculated by64 ΔGpolar = −

μ2 2(ε − 1) · 6a3 2ε + 1

(4)

where μ is the dipole moment of solute and a is the sphere radius corresponding to the volume of the compound. Thus, ΔGpolar‑Impl is expected to correlate with ε, whereas it is less clear whether ΔGpolar‑Expl correlates with ε. Figure 5 shows ΔGpolar‑Impl (GB and PB models) and ΔGpolar‑Expl (TI5,6) for different solutes as a function of solvent 2(ε − 1)/(2ε + 1), and Figure 6 does likewise for SMD at different levels of theory. Correlation between these properties for each solute with N (number of solvents) ≥ 10 is given in Figure 3 and Tables S34−S43. For PB and GB, ΔGpolar‑Impl decreases as ε increases (Figure 5) and ΔGpolar‑Impl correlates strongly with 2(ε − 1)/ (2ε + 1): R is approximately 1.0 for all solutes except for toluene (R = 0.92 ± 0.10, Still). Interestingly, similar results were found in SMD calculations, but not for some solutes, including acetone, and especially ethanol, methanol, and nitromethane, where poor correlations were found. After carefully checking the data points, we found that a poor correlation is found when a solute is capable of forming Hbonds with some of the tested solvents. Moreover, ΔGpolar‑Expl is not a monotonic function of ε. The correlation between 2(ε − 1)/(2ε + 1) and ΔGpolar‑Expl is significantly weaker: acetone has the strongest correlation (N = 17, R = 0.87 ± 0.07), and nitromethane has the weakest correlation (N = 12, R = 0.31 ± 0.21). Among the eight solutes, nitromethane has the largest RMSD and MSE values of ΔGpolar‑Impl from ΔGpolar‑Expl: RMSD = 7.3 ± 0.9 kJ/mol and MSE = −6.5 ± 0.9 kJ/mol for PB and RMSD at 24−28 kJ/mol and MSE ≈ −24 kJ/mol for GB. As disscussed in our previous paper,5 charges on the nitro group are likely too high, which may lead to an overestimated ΔGsolv in the explicit solvent calculation. ΔGpolar‑Impl is more negative than ΔGpolar‑Expl, meaning that the overestimation is magnified in PB, especially in GB. The BEDROC value and the p value (Table S24) of the Student’s t test (Table S28) suggest that the unsigned error for nitromethane as a solute is systematic. For acetone and acetonitrile, the BEDROC values and p values of

Figure 4. Correlation between the polar contribution to solvation free energy calculated with explicit solvent (ΔGpolar‑Expl) and that calculated with implicit solvent (ΔGpolar‑Impl): PB (black), (A) Still (red), (B) HCT (green), (C) OBC-I (blue), and (D) OBC-II (orange).

comparation of correlations indicates that the difference between PB and GB is significant, whereas the difference between GB models is not significant. PB gives results closest to explicit solvent calculations with an RMSD of 4.3 ± 0.2 kJ/ mol and MSE of −3.1 ± 0.2 kJ/mol (Table S3), whereas GB models have an RMSD of 12.4−14.5 kJ/mol and MSE < −10 kJ/mol (Tables S4−S7). In addition, the RMSD and MSE of ΔGpolar per solvent given in Tables S3−S7 suggest that PB gives smaller RMSD and MSE values for each solvent as well. The correlation between ΔGpolar‑Expl and ΔGpolar‑Impl per solvent with N (number of solutes) ≥ 10 is shown in Figure S1, and the linear fit parameters are given in Tables 2 and S29− 1038


Article


Table 2. Correlation Coefficients Rpolar between ΔGpolar‑Expl and ΔGpolar‑Impl Per Solvent and Rapolar between ΔGapolar‑Expl and ΔGapolar‑Impl Per Solventa Rpolar

a

solvent

N

ε

toluene chloroform ethyl acetate 1-octanol acetone ethanol methanol nitromethane

11 28 13 29 12 13 20 10

2.37 4.71 5.99 9.86 20.49 25.00 31.49 36.56

PB 0.97 0.77 0.90 0.83 0.98 0.96 0.97 0.98

± ± ± ± ± ± ± ±

Rapolar Still

0.04 0.08 0.05 0.04 0.01 0.02 0.01 0.02

0.86 0.85 0.85 0.73 0.96 0.92 0.94 0.98

± ± ± ± ± ± ± ±

0.15 0.05 0.08 0.05 0.02 0.05 0.02 0.02

PB 1.00 1.00 0.97 0.98 0.99 0.99 0.99 0.92

± ± ± ± ± ± ± ±

Still 0.02 0.00 0.05 0.01 0.03 0.01 0.01 0.06

0.91 0.91 0.67 0.81 0.79 0.78 0.83 0.57

± ± ± ± ± ± ± ±

0.13 0.06 0.23 0.08 0.17 0.15 0.11 0.28

N is the number of solutes tested. ε is the experimental dielectric constant.

not consistently well correlated with ΔGapolar‑Expl (Tables 2 and S46−48). Meanwhile, the correlation per solute for both PB and GB given in Figure 8 is much lower (Table 3). The results suggest that reconsidering the solvent-dependent parameters could improve the predictive power of ΔGapolar using PB/GB implicit solvent models. All of the solvent molecules tested in this work are larger than water; hence, a loss of accuracy may be caused by not adjusting the probe radii used to augment the vdW surface for construction of the molecular surface.79,80 The SASA, SAV, and ΔGatt terms were found to be (nearly perfectly) linearly dependent on the probe radius. Since it is not straightforward to determine the probe radii to be used for nonspherical solvent molecules, we instead optimized the surface tension parameter for the SASA model (the apolar part in GB models) and applied a scaling factor for the PB model (the results are given in section 9 of the Supporting Information). The agreement between the optimized solvation free energy calculated with implicit solvent and experimental data (and explicit solvent data) improved significantly because of this (Tables S49 and S50). PB now yields slightly better results than SMD at four levels of theory. For GB, the RMSD from experimental data is comparable with SMD calculation at the HF/6-31+G(d,p) level of theory, and the correlation is significantly higher. It should be noted that for GB, the optimized γ is negative for all solvents tested in this work; thus, simply optimizing the probe radius may not be sufficient to help improve GB calculations. Since a constant c is often used to model ΔGapolar (eq 2), we used that here too, although with the reservation that c has no physical meaning since it essentially represents thesolvent dependentΔGapolar of the vacuum.

Student’s t test highlight that the unsigned errors are systematic in GB but not in PB. 3.3. Evaluation of ΔGapolar. For GB, ΔGapolar‑Impl is calculated by applying eq 2. Since both γ and SASA are positive values, ΔGapolar‑Impl is positive as well. However, ΔGapolar calculated with explicit solvent (ΔGapolar‑Expl) has the opposite sign in most cases; hence, the slope of the linear fits in Figure 7 is negative. As noted in the Methods section, ΔGapolar‑Expl is not well-defined physically, and this might be the cause of the lacking correlation. On the other hand, there may be an apolar favorable free energy contribution due to an increased entropy that is modeled by the explicit but not by the implicit solvent description. For PB, ΔGapolar‑Impl is calculated based on eq 3, where both SASA and SAV are included to reproduce the full SPT model; in particular, the attractive van der Waals interaction term (ΔGatt) is included.45 In eq 3, the SASA and SAV terms are aimed at modeling the repulsive component of the van der Waals interaction, where ΔGrep is positive and ΔGatt is negative. The GB results stress the importance of including an attractive component like ΔGatt. However, using eq 3, ΔGapolar‑Impl for virtually all tested systems remains positive (data not shown) and correlates poorly with ΔGapolar‑Expl (R = 0.25 ± 0.06). The results suggest that either ΔGatt is underestimated or ΔGrep is overestimated. It is possible that ΔGatt is underestimated since ρ in eq 3 represents the solvent number density for water, which contains one (oxygen) atom contributeing to ΔGatt. In this work, each solvent molecule has more than one heavy atom contributing to ΔGatt. Thus, we replaced ρ by the product of the solvent number density and the number of heavy atoms per solvent molecule as a crude correction. Indeed, the correlation is significantly improved (R = 0.69 ± 0.03), and the slope becomes positive (Figure 7). The RMSD is 10.6 ± 0.4 kJ/mol, and MSE is −8.5 ± 0.4 kJ/mol (Table S3), suggesting that ΔGapolar is an important source of error in PB calculations. For a solute in conformation x ⃗ , eq 5 gives a more realistic value of ΔGatt, where uatt ij (x ⃗ , y ⃗ ) is the attractive van der Waals potential in a WCA decomposition scheme for solute atom i and solvent atom j at position y ⃗ . ΔGatt(x ⃗) = ρ ̅

∑∑∫ i

j

Ω

4. DISCUSSION The predictive power of implicit solvent models, including PB and GB, has been extensively tested in aqueous media,36,58,59 and reasonable agreement with explicit solvent simulations as well as experiments has been reported. However, the accuracy of these methods for calculating binding free energies is not convincing.19,40,48,96 Here, we calculated the solvation free energy for 226 combinations of organic solvents and solutes using PB and several types of widely used GB models (Still, HCT, OBC-I, and OBC-II). We compared the results to our previous simulations with explicit solvent models5,6 as well as experimental data. The comparison indicates poor correlation (R ≈ 0.7 for PB and 0.4 for GB; Figure 2 and Table 1) between PB/GB and explicit solvent calculations and poorer correlation (R ≈ 0.5 for PB and 0.2 for GB; Figure 3 and Table 1) when comparing to experimental data. The correlation between

uijatt(xi⃗ , y ⃗ ) θ(x ⃗ , y ⃗ ) dy ⃗ (5)

The correlation of ΔGapolar per implicit solvent with the explicit case where N ≥ 10 is given in Figure S2 (including the correction for the number of heavy atoms). For PB, the correlation coefficient is nearly 1.0 (Table 2), except for nitromethane (N = 10, R = 0.92 ± 0.06). For GB, ΔGapolar‑Impl is 1039


Article


Figure 5. Correlation between 2(ε − 1)/(2ε + 1) and the polar contribution to solvation free energy (ΔGpolar) calculated with PB (black), Still (red), HCT (green), OBC-I (blue), OBC-II (orange), and explicit solvent (TI) (purple) models per solute: (A) acetone, (B) acetonitrile, (C) chloroform, (D) dichloromethane, (E) ethanol, (F) methanol, (G) nitromethane, and (H) toluene.

Figure 6. Correlation between 2(ε − 1)/(2ε + 1) and the polar contribution to solvation free energy (ΔGpolar) calculated with explicit solvent (TI) (black) and the SMD model with HF/6-31+G(d,p) (red), B3LYP/6-31+G(d,p) (green), M06-2X/6-31+G(d,p) (blue), and M06-2X/cc-pVTZ (orange) levels of theory per solute: (A) acetone, (B) acetonitrile, (C) chloroform, (D) dichloromethane, (E) ethanol, (F) methanol, (G) nitromethane, and (H) toluene.

experiment and explicit solvent was previously found to be 0.84.5,6 Despite PB yielding statistically somewhat stronger correlation to explicit solvent than GB, all force field-based implicit solvent models give poor agreement with experimental data (RMSD is 14−15 kJ/mol). The SMD model, which was developed as a universal implicit solvent model based on quantum calculations,15 was tested as well. As reported in the literature,97 the performance of SMD tested here (R = 0.7−0.8; Table 1) is significantly better than other implicit solvent models at all four levels of theory, but it is not as good as the explicit solvent model. The SMD model differs from the PB and GB models in that solvent polarization is included and that

specific solvent parameters are used, beyond the dielectric constant ε. However, the SMD model still lacks explicit hydrogen bonding and most likely overestimates the strength of salt bridges, if present. In addition, the entropy of the solvent is difficult to model implicitly.98,99 As noted in our earlier work,5,6 more elaborate implicit models like COSMO-RS50−52 do in fact provide even better correlation than TI in explicit solvent when the GAFF force field100 employed. A decomposition of ΔG solv suggests that the large discrepancy between PB/GB calculations and explicit solvent calculations can be mainly ascribed to ΔGapolar (Tables S3−S7). 1040


Article


Table 3. Correlation Rpolar between 2(ε − 1)/(2ε + 1) and ΔGpolar Per Solute and Correlation Rapolar between ΔGapolar‑Expl and ΔGapolar‑Impl Per Solute Rpolar solute

N

acetone acetonitrile chloroform dichloromethane ethanol methanol nitromethane toluene

17 12 15 14 19 17 12 17

TI 0.87 0.67 0.39 0.63 0.77 0.77 0.31 0.67

± ± ± ± ± ± ± ±

PB 0.07 0.17 0.22 0.28 0.15 0.16 0.21 0.16

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

± ± ± ± ± ± ± ±

Rapolar Still

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.92

± ± ± ± ± ± ± ±

0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.10

B3LYP/6-31+G(d,p) 0.82 1.00 1.00 1.00 0.60 0.57 0.60 1.00

± ± ± ± ± ± ± ±

0.08 0.00 0.00 0.00 0.12 0.15 0.20 0.00

PB 0.27 0.09 0.32 0.55 0.32 0.22 0.28 0.28

± ± ± ± ± ± ± ±

Still 0.21 0.19 0.22 0.25 0.21 0.21 0.23 0.21

0.64 0.55 0.40 0.15 0.50 0.21 0.26 0.13

± ± ± ± ± ± ± ±

0.15 0.16 0.14 0.12 0.15 0.16 0.18 0.14

performance of the optimized GB model is significantly improved (RMSD from experimental data is 8−9 kJ/mol), the correlation remains poor (R = 0.5−0.6). Further corrections may be required where the dependence of parameters on the solvent dielectric constant should be taken into consideration in GB models.53 For PB, a scale factor was used to optimize the apolar contribution per solvent. The performance of the optimized PB model is slightly better (R = 0.77 ± 0.04, RMSD = 4.5 ± 0.24 kJ/mol) than the SMD model at the four levels of theory tested. In addition, the functional form of ΔGatt in eq 3 is not sufficient for solvents that have more than one heavy atom contributing to ΔGatt, and the solvent atom type may need to be taken into consideration for parametrization. As expected, for PB and GB, ΔGpolar‑Impl correlates strongly with 2(ε − 1)/(2ε + 1) (R ≥ 0.95; Table 3 and Figure 5) and with ΔGpolar‑Expl (R ≈ 0.86 for GB and 0.92 for PB; Figure 4), confirming that ε is a powerful predictor. ΔGpolar‑Impl for PB gives better agreement with ΔGpolar‑Expl (RMSD ∼ 4 kJ/mol; Table S3) than GB (RMSD > 12 kJ/mol; Tables S4−S7), which is not surprising since the GB model is an approximation of PB. Systematic errors are found for solutes including acetone, acetonitrile, and nitromethane in GB. The accuracy of the GB model crucially depends on the effective Born radii used.29 However, the dependence of effective Born radii on the solvent dielectric constant is neglected for the standard GB model, which is designed and parametrized for aqueous media of ε ∼ 80.12 Alternative GB models have been proposed to reproduce PB results with high accuracy through different approaches: explicitly considering the solvent influence on effective Born radii,53 by approximating the analytical solution for arbitrary molecular shape starting from Kirkwood’s exact solution for the case of a random charge distribution in a sphere,54,55 and by directly scaling ΔGsolv.101 However, ΔGpolar‑Expl is not consistently well correlated with 2(ε − 1)/ (2ε + 1) in our tests: R = 0.87 for acetone and R < 0.4 for chloroform and nitromethane (Table 3 and Figure 5). Although it uses a different definition of the polar contribution compared with other models discussed here, the more sophisticated implicit solvent model SMD also gives relatively poor correlation for solutes that may form H-bonds: R = 0.60 for ethanol, 0.57 for methanol, and 0.60 for nitromethane (Table 3 and Figure 6). These findings indicate that, unlike PB or GB calculations, the correlation is system dependent, implying that the Onsager relation does not hold for realistic systems. Despite differences in the description of electrostatic interactions, it remains difficult to see how any of the implicit solvent models studied could quantitatively model

Figure 7. Correlation between apolar contribution to solvation free energy calculated with explicit solvent (ΔGapolar‑Expl) and that calculated with the PB (black), (A) Still (red), (B) HCT (green), (C) OBC-I (blue), and (D) OBC-II (orange) models (ΔGapolar‑Impl).

We find that ΔGapolar‑Expl is not directly proportional to SASA (in contrast to ΔGapolar‑Impl; Figure 7); on the contrary, the nonpolar interaction is more favorable to solvation as SASA increases (Figure S2). It should also be noted that the per solvent analysis indicates that SASA is not consistently well correlated with ΔGapolar‑Expl (GB models in Table 2). The introduction of the SAV and particularly ΔGatt terms in addition to SASA in the PB model reduces the discrepancy between ΔGapolar‑Impl and ΔGapolar‑Expl (Tables S3) and significantly improves the correlation per solvent (Table 2 and Figure S2). However, the correlation per solute remains rather weak, and the RMSD and negative MSE values (with an opposite sign compared with that from the SASA model) are considerably larger than ΔGpolar for PB. The apolar term in the SMD model is optimized for various solvents, and this is likely to be part of the reason why the SMD model performs significantly better than PB or GB. As proposed by Bordner et al. for the SASA approach to calculate ΔGsolv in 1-octanol and hexadecane,56 we also find that the optimized γ for the organic solvents tested in this work is negative. The result suggest that simply optimizing the solvent probe radius may not be sufficient for improving the GB model. Although the 1041



■

Article

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (Q.W.). *E-mail: [email protected] (D.v.d.S.). ORCID

David van der Spoel: 0000-0002-7659-8526 Funding

The National Natural Science Foundation of China is acknowledged for financial support to J.Z., Q.W. (21673206), and H.Z. (21606016). H.Z. was also supported in part by the China Postdoctoral Science Foundation (2015M580993). J.Z. thanks Lu Tan from Qi Wang’s group for technical support. A grant from the Swedish Research Council (2013-5947) to D.v.d.S. is acknowledged as well. Notes

The authors declare no competing financial interest.

■

(1) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. J. Chem. Theory Comput. 2012, 8, 61−74. (2) Hub, J. S.; Caleman, C.; van der Spoel, D. Phys. Chem. Chem. Phys. 2012, 14, 9537−9545. (3) Duffy, E.; Jorgensen, W. J. Am. Chem. Soc. 2000, 122, 2878−2888. (4) Chebil, L.; Chipot, C.; Archambault, F.; Humeau, C.; Engasser, J. M.; Ghoul, M.; Dehez, F. J. Phys. Chem. B 2010, 114, 12308−12313. (5) Zhang, J.; Tuguldur, B.; van der Spoel, D. J. Chem. Inf. Model. 2015, 55, 1192−1201. (6) Zhang, J.; Tuguldur, B.; van der Spoel, D. J. Chem. Inf. Model. 2016, 56, 819−820. (7) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8592. (8) Wennberg, C. L.; Murtola, T.; Hess, B.; Lindahl, E. J. Chem. Theory Comput. 2013, 9, 3527−3537. (9) Fischer, N. M.; van Maaren, P. J.; Ditz, J. C.; Yildirim, A.; van der Spoel, D. J. Chem. Theory Comput. 2015, 11, 2938−2944. (10) Zubillaga, R. A.; Labastida, A.; Cruz, B.; Martínez, J. C.; Sánchez, E.; Alejandre, J. J. Chem. Theory Comput. 2013, 9, 1611−1615. (11) Wang, L.-P.; Martinez, T. J.; Pande, V. S. J. Phys. Chem. Lett. 2014, 5, 1885−1891. (12) Chen, J.; Brooks, C. L.; Khandogin, J. Curr. Opin. Struct. Biol. 2008, 18, 140−148. (13) Tomasi, J.; Persico, M. Chem. Rev. 1994, 94, 2027−2094. (14) Orozco, M.; Luque, F. J. Chem. Rev. 2000, 100, 4187−4226. (15) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999−3094. (16) Ren, P.; Chun, J.; Thomas, D. G.; Schnieders, M. J.; Marucho, M.; Zhang, J.; Baker, N. A. Q. Rev. Biophys. 2012, 45, 427−491. (17) Larsson, D. S. D.; van der Spoel, D. J. Chem. Theory Comput. 2012, 8, 2474−2483. (18) Larsson, P.; Lindahl, E. J. Comput. Chem. 2010, 31, 2593−2600. (19) Zhang, H.; Tan, T.; van der Spoel, D. J. Chem. Theory Comput. 2015, 11, 5103−5113. (20) Rhee, Y. M.; Sorin, E. J.; Jayachandran, G.; Lindahl, E.; Pande, V. S. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 6456−6461. (21) Davis, M. E.; McCammon, J. A. Chem. Rev. 1990, 90, 509−521. (22) Honig, B.; Nicholls, A. Science 1995, 268, 1144−1149. (23) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (24) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127−6129. (25) Qiu, D.; Shenkin, P. S.; Hollinger, F. P.; Still, W. C. J. Phys. Chem. A 1997, 101, 3005−3014. (26) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996, 100, 19824−19839. (27) Schaefer, M.; Karplus, M. J. Phys. Chem. 1996, 100, 1578−1599.

Figure 8. Correlation between apolar contributions to solvation calculated with an explicit solvent model (ΔGapolar‑Expl) and that calculated with the PB (black), Still (red), HCT (green), OBC-I (blue), and OBC-II (orange) models (ΔGapolar‑Impl) per solute: (A) acetone, (B) acetonitrile, (C) chloroform, (D) dichloromethane, (E) ethanol, (F) methanol, (G) nitromethane, and (H) toluene.

H-bonds or solvent entropy and, as a result, predict ΔGsolv for arbitrary solvents.

■

REFERENCES

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00169. Detailed tables showing the agreement between calculated solvation free energies with implicit solvent and experimental data (calculated with explicit solvent) (PDF) 1042


Article

Journal of Chemical Theory and Computation (28) Onufriev, A.; Bashford, D.; Case, D. A. Proteins: Struct., Funct., Genet. 2004, 55, 383−394. (29) Onufriev, A.; Case, D. A.; Bashford, D. J. Comput. Chem. 2002, 23, 1297−1304. (30) Aguilar, B.; Onufriev, A. V. J. Chem. Theory Comput. 2012, 8, 2404−2411. (31) Pierotti, R. A. Chem. Rev. 1976, 76, 717−726. (32) Stillinger, F. H. J. Solution Chem. 1973, 2, 141−158. (33) Gallicchio, E.; Kubo, M. M.; Levy, R. M. J. Phys. Chem. B 2000, 104, 6271−6285. (34) Tan, C.; Tan, Y.-H.; Luo, R. J. Phys. Chem. B 2007, 111, 12263− 12274. (35) Mobley, D. L.; Dill, K. A. Structure 2009, 17, 489−498. (36) Shivakumar, D.; Deng, Y.; Roux, B. J. Chem. Theory Comput. 2009, 5, 919−930. (37) Underwood, R.; Tomlinson-Phillips, J.; Ben-Amotz, D. J. Phys. Chem. B 2010, 114, 8646−8651. (38) Levy, R. M.; Zhang, L. Y.; Gallicchio, E.; Felts, A. K. J. Am. Chem. Soc. 2003, 125, 9523−9530. (39) Genheden, S.; Mikulskis, P.; Hu, L.; Kongsted, J.; Söderhjelm, P.; Ryde, U. J. Am. Chem. Soc. 2011, 133, 13081−13092. (40) Harris, R. C.; Pettitt, B. M. J. Chem. Theory Comput. 2015, 11, 4593−4600. (41) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 1978−1988. (42) Ashbaugh, H. S.; Kaler, E. W.; Paulaitis, M. E. J. Am. Chem. Soc. 1999, 121, 9243−9244. (43) Raschke, T. M.; Tsai, J.; Levitt, M. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 5965−5969. (44) Rajamani, S.; Truskett, T. M.; Garde, S. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 9475−9480. (45) Wagoner, J. A.; Baker, N. A. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 8331−8336. (46) Trott, O.; Olson, A. J. J. Comput. Chem. 2009, 31, 455−461. (47) Ruiz-Carmona, S.; Alvarez-Garcia, D.; Foloppe, N.; GarmendiaDoval, A. B.; Juhos, S.; Schmidtke, P.; Barril, X.; Hubbard, R. E.; Morley, S. D. PLoS Comput. Biol. 2014, 10, e1003571. (48) Zhang, H.; Yin, C.; Yan, H.; van der Spoel, D. J. Chem. Inf. Model. 2016, 56, 2080−2092. (49) Frisch, M. J.; et al. Gaussian 09, Revision B.01; Gaussian Inc.: Wallingford, CT, 2009. (50) Klamt, A. J. Phys. Chem. 1995, 99, 2224−2235. (51) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. J. Phys. Chem. A 1998, 102, 5074−5085. (52) Klamt, A. WIREs Comput. Mol. Sci. 2011, 1, 699−709. (53) Feig, M.; Im, W.; Brooks, C. L. J. Chem. Phys. 2004, 120, 903. (54) Sigalov, G.; Scheffel, P.; Onufriev, A. J. Chem. Phys. 2005, 122, 094511. (55) Sigalov, G.; Fenley, A.; Onufriev, A. J. Chem. Phys. 2006, 124, 124902. (56) Bordner, A. J.; Cavasotto, C. N.; Abagyan, R. A. J. Phys. Chem. B 2002, 106, 11009−11015. (57) Salari, R.; Chong, L. T. J. Phys. Chem. Lett. 2010, 1, 2844−2848. (58) Mobley, D. L.; Barber, A. E., II; Fennell, C. J.; Dill, K. A. J. Phys. Chem. B 2008, 112, 2405−2414. (59) Knight, J. L.; Brooks, C. L. J. Comput. Chem. 2011, 32, 2909− 2923. (60) Spassov, V. Z.; Yan, L.; Szalma, S. J. Phys. Chem. B 2002, 106, 8726−8738. (61) Feig, M.; Brooks, C. L. Curr. Opin. Struct. Biol. 2004, 14, 217− 224. (62) Kleinjung, J.; Fraternali, F. Curr. Opin. Struct. Biol. 2014, 25, 126−134. (63) Botello-Smith, W. M.; Luo, R. J. Chem. Inf. Model. 2015, 55, 2187−2199. (64) Bell, R. P. Trans. Faraday Soc. 1931, 27, 797. (65) Onsager, L. J. Am. Chem. Soc. 1936, 58, 1486−1493. (66) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2009, 113, 6378−6396.

(67) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43−56. (68) Lindahl, E.; Hess, B.; Van Der Spoel, D. J. Mol. Model. 2001, 7, 306−317. (69) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701−1718. (70) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435−447. (71) 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. Bioinformatics 2013, 29, 845−54. (72) van der Spoel, D.; van Maaren, P. J.; Caleman, C. Bioinformatics 2012, 28, 752−753. (73) Ghahremanpour, M.; van Maaren, P.; Ditz, J.; Lindh, R.; van der Spoel, D. J. Chem. Phys. 2016, 145, 114305. (74) Nosé, S. Mol. Phys. 1984, 52, 255−268. (75) Hoover, W. G. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (76) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463−1472. (77) Bennett, C. H. J. Comput. Phys. 1976, 22, 245−268. (78) Schaefer, M.; Bartels, C.; Karplus, M. J. Mol. Biol. 1998, 284, 835−848. (79) Lee, B.; Richards, F. J. Mol. Biol. 1971, 55, 379−400. (80) Richmond, T. J. J. Mol. Biol. 1984, 178, 63−89. (81) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237−5247. (82) Shrake, A.; Rupley, J. J. Mol. Biol. 1973, 79, 351−371. (83) Martins, S. A.; Sousa, S. F. J. Comput. Chem. 2013, 34, 1354− 1362. (84) Radom, L.; Schleyer, P. v. R.; Pople, J. A.; Hehre, W. J. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (85) Becke, A. D. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098− 3100. (86) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (87) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (88) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215−241. (89) Dunning, T. H., Jr J. Chem. Phys. 1989, 90, 1007−1023. (90) Truchon, J.-F.; Bayly, C. I. J. Chem. Inf. Model. 2007, 47, 488− 508. (91) Swamidass, S. J.; Azencott, C.-A.; Daily, K.; Baldi, P. Bioinformatics 2010, 26, 1348−1356. (92) R Core Team R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2015. (93) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (94) Born, M. Eur. Phys. J. A 1920, 1, 45−48. (95) Kirkwood, J. G. J. Chem. Phys. 1934, 2, 351. (96) Harris, R. C.; Mackoy, T.; Fenley, M. O. J. Chem. Theory Comput. 2015, 11, 705−712. (97) Kolár,̌ M.; Fanfrlík, J.; Lepšík, M.; Forti, F.; Luque, F. J.; Hobza, P. J. Phys. Chem. B 2013, 117, 5950−5962. (98) Caleman, C.; Hub, J. S.; van Maaren, P. J.; van der Spoel, D. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 6838−6842. (99) Hub, J. S.; Wolf, M. G.; Caleman, C.; van Maaren, P. J.; Groenhof, G.; van der Spoel, D. Chem. Sci. 2014, 5, 1745−1749. (100) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (101) Tjong, H.; Zhou, H.-X. J. Phys. Chem. B 2007, 111, 3055− 3061.

1043


Comparison of Implicit and Explicit Solvent Models for the Calculation

Recommend Documents