Performance of SM8 on a test to predict small-molecule solvation free energies.

The SM8 quantum mechanical aqueous continuum solvation model is applied to a 17-molecule test set proposed by Nicholls et al. (J. Med. Chem. 2008, 51, 769) to predict free energies of solvation. With the M06-2X density functional, the 6-31G(d) basis set, and CM4M charge model, the root-mean-square error (RMSE) of SM8 is 1.08 kcal mol(-1) for aqueous geometries and 1.14 kcal mol(-1) for gas-phase geometries. These errors compare favorably with optimal explicit and continuum models reported by Nicholls et al., having RMSEs of 1.33 and 1.87 kcal mol(-1), respectively. Other models examined by these workers had RMSEs of 1.5-2.6 kcal mol(-1). We also explore the use of other density functionals and charge models with SM8 and the RMSE increases to 1.21 kcal mol(-1) for mPW1/CM4 with gas-phase geometries, to 1.50 kcal mol(-1) for M06-2X/CM4 with gas-phase geometries, and to 1.27-1.64 kcal mol(-1) with three different models at B3LYP gas-phase geometries.

Introduction

The effect of solvation on molecular structure, energetics, and dynamics can be decisive in the medicinal chemistry of drugs, druglike molecules, and biomolecules in aqueous solution as well as for numerous other technological and environmental applications. Nicholls et al.(1) recently presented a set of 17 small molecules that were proposed as a standard set against which to test computational models for the prediction of free energies of solvation. The test molecules were chosen on the basis of measured values being available and the extent to which they incorporated one or more druglike functional groups. The molecules in the test set are listed in Table 1.

Table 1

Aqueous Solvation Free Energies (kcal mol−1) from Two Continuum Solvation Models, SM8 and Poisson/ZAP-9,(1) Using Geometries Optimized in the Gas Phase and in the Solution Phase, from an Explicit-Solvent Model, and from Experiment

molecule	SM8 B3LYP-opt.a	SM8 gas-opt.b	SM8 soln-opt.c	Poissond	explicite	exptf
benzyl bromide	−2.70	−2.64	−2.67	−2.71	−1.88	−2.38
benzyl chloride	−2.53	−2.48	−2.56	−2.92	−2.09	−1.93
bis(2-chloroethyl) ether	−2.59	−2.72	−2.86	−2.28	−4.25	−4.23
1,1-diacetoxyethane	−6.13	−6.71	−7.42	−6.74	−7.38	−4.97
1,1-diethoxyethane	−1.75	−1.98	−2.12	−2.06	−3.21	−3.28
1,2-diethoxyethane	−2.47	−2.83	−3.05	−1.87	−4.74	−3.54
diethyl propanedioate	−5.69	−5.94	−6.30	−6.40	−7.69	−6.00
diethyl sulfide	−0.53	−0.47	−0.52	−1.15	−1.47	−1.43
dimethoxymethane	−2.14	−2.59	−2.81	−3.04	−4.58	−2.93
N,N-dimethyl-p-methoxybenzamide	−8.29	−8.61	−9.47	−6.96	−10.40	−11.01
1,4-dioxane	−4.96	−5.32	−5.45	−3.93	−5.79	−5.05
ethylene glycol diacetate	−7.35	−6.60	−7.21	−7.22	−8.30	−6.34
glycerol triacetate	−9.96	−9.55	−10.27	−10.70	−11.95	−8.84
imidazole	−9.16	−9.34	−9.52	−10.01	−10.18	−9.81
phenyl formate	−3.68	−4.14	−4.37	−5.93	−4.28	−3.82
m-bis(trifluoromethyl) benzene	1.11	0.90	0.78	−0.75	0.06	1.07
N,N,4-trimethylbenzamide	−7.13	−7.24	−7.97	−5.57	−8.65	−9.76
MSEg	0.49	0.35	0.03	0.24	−0.74
MUEh	0.98	0.86	0.88	1.47	1.01
RMSEi	1.25	1.14	1.08	1.88	1.33

SM8 B3LYP-opt. implies SM8/CM4M/M06-2X/6-31G(d)//B3LYP(g)/6-31G(d,p), where (g) emphasizes that the geometry was optimized in the gas phase.

SM8 gas-opt. implies SM8/CM4M/M06-2X/6-31G(d)//M06-2X(g)/6-31G(d,p).

SM8 soln-opt. implies SM8/CM4M/M06-2X/6-31G(d)//M06-2X/6-31G(d,p) with the geometry optimized in solution.

Poisson implies use of Open-Eye version 1 AM1-BCC charges and ZAP-9 radii.

Explicit implies GAFF/TIP3P calculations with Merck−Frosst version 2.6A24 AM1-BCC charges.

Experiment is from Nicholls et al.(1)

MSE is the mean signed error.

MUE is the mean unsigned error.

RMSE is the root-mean-squared error.

Nicholls et al.(1) examined two computational formalisms for predicting aqueous solvation free energies, namely, explicit[2-4] and implicit[4-6] solvation models.(7) In the former case, the free energy was computed by analysis of molecular dynamics trajectories coupled to the effective annihilation of the solute in a periodic box of explicit water molecules. In the latter case, the molecular nature of the solvent was ignored, and instead the surrounding space was represented by a continuum characterized by the appropriate dielectric constant for the solvent. Full details of the particular protocols used by Nicholls et al. are available in their original article(1) and are not repeated here. A critical aspect of the models explored by Nicholls et al., however, was that in both protocols the solutes (and molecular solvent, when explicit) were represented using molecular mechanics (MM). That is, intermolecular interactions were governed by interactions between partial atomic charges and typical MM nonbonded interactions (e.g., Lennard-Jones terms),(7) and the choice of charge model used to compute atomic partial charges had significant influence on the accuracy of the predicted intermolecular interactions. In the case of implicit solvation, their calculations were sensitive to the assignment of atomic radii from which cavities that separate the solute from the surrounding continuum were constructed. Nicholls et al.(1) reported results for various choices of charges and radii. Their optimal protocol for explicit solvation used the general Amber force field (GAFF),(8) the three-point transferable intermolecular potential (TIP3P) water model,(9) and partial atomic charges taken from version 2.6A24 of the Merck−Frosst Austin model 1 bond charge correction (AM1-BCC) model.(10) For implicit solvation, their optimal protocol employed version 1 AM1-BCC charges and determined the reaction field from numerical solution of the nonhomogeneous Poisson equation using a cavity defined by so-called ZAP 9 atomic radii (these radii having previously been optimized for the prediction of aqueous solvation free energies over a different test set of 200 small molecules). The root-mean-squared errors (RMSEs) reported by Nicholls et al.(1) for these explicit and implicit protocols were 1.33 and 1.87 kcal mol−1, respectively. Other models examined by these workers had errors of 1.53−2.57 kcal mol−1.

Agreeing wholeheartedly with the assertion of Nicholls et al.(1) that standardized test sets are useful for the benchmarking of alternative modeling protocols, we report here results for the implicit SM8 quantum mechanical continuum solvation model. Although SM8[11,12] is parametrized for both aqueous and nonaqueous solvents, the present test involves only aqueous solvation. In the spirit of viewing the test as “blind”, we report results only for previously published parameters. Our model that would be expected to be most generally reliable involves use of the M06-2X(13) density functional, the 6-31G(d)(14) basis set, and the charge model 4/M06 suite (CM4M)(11) charge model with gas-phase geometries optimized with the M06-2X(13) density functional. We recommend this protocol because it allows the user to take advantage of an accurate and reliable level of theory, M06-2X, along with a charge model, CM4M, that is specifically designed to reproduce accurate partial atomic charges using the whole M06 suite[13,15] of four density functionals. We use the 6-31G(d)(14) basis set because it produces stable partial atomic charges and reasonable gas-phase geometries. Since accurate partial atomic charges are critical for the SM8 to reliably reproduce or predict experimental free energies of solvation, we recommend the above combination with SM8 whenever it is feasible. We recognize however that it is often necessary to use other levels of theory, so we have also developed a more general alternative, CM4, which is designed to produce accurate partial atomic charges independent of the specific self-consistent-field model for a variety of basis sets. For the purpose of comparing solvation methods with the same solute geometry as was employed by the previous workers, we also use the M06-2X density functional with the same molecular geometries as those employed by Nicholls et al.(1) We briefly discuss the sensitivity of the SM8 model to choices of molecular geometry, charge model, and density functional model. We consider only the 6-31G(d) basis set because our experience has shown that this partially polarized valence double-ζ basis set is adequate to yield reasonable molecular geometries and stable partial atomic charges, whereas partial atomic charges computed by the convenient population analyses underlying the CM4 and CM4M models become more and more erratic as basis sets become larger and larger. In addition, calculations with larger basis set are more computationally demanding.

Computational Methods

The SM8 solvation model has been presented recently in full detail.(12) We note that for aqueous solution the SM8 model as very similar (not only in methods but also in numerical results) to the earlier SM6(16) model; the only significant difference between SM6 and SM8 is that SM6 is defined only for aqueous solutions whereas SM8 is defined for general solvents; this generalization results in a minor difference in the predictions for aqueous solution free energies.

The free energy of interaction between a solute and a solvent depends upon the reaction field induced in the solvent by the charge distribution of the solute. A key aspect of the SM8 model that is germane to the present application is that it does not obtain the solute reaction field by solving the nonhomogeneous Poisson equation, but instead from the generalized Born (GB)[17-23] approximation. In GB models, the bulk-electrostatic component of the free energy of solvation GP is computed according to where ε is the solvent dielectric constant, the summation runs over all atoms k and k′ having partial atomic charges q and q, and γ is an effective Coulomb integral having units of inverse distance that depends on atomic radii in much the same way as does the construction of a molecular cavity for solution of the nonhomogeneous Poisson equation. Details of the algorithm are given elsewhere.[24-26] All the free energies of solvation reported here correspond to a 1 M ideal gas and a 1 M ideal infinitely dilute solution. This convention eliminates any change in the translational entropy upon solvation, as discussed by Ben-Naim.(27)

In the SM8 solvation model, charge models 4 and 4M[11,16] (CM4 or CM4M) are used to determine partial atomic charges. The CM4 and CM4M charge models take as input partial atomic charges determined by Löwdin population analysis[28,29] of quantum mechanical density matrices and map them to improved charges. The mapping is based on parameters that were previously optimized[11,16] so that the CM4 or CM4M charges accurately reproduce experimental gas-phase dipole moments. In this respect, the models are similar in spirit to the AM1-BCC model. The CM4 model was developed to be general in the sense that it was designed for use with any “accurate” density, and the relevant parameters were optimized by considering data from calculations using a range of density functionals.(16) It represents, therefore, a compromise between, on the one hand, a flexible and convenient set of parameters that can be used with a broad choice of methods and, on the other hand, a fully optimized set of parameters. The CM4M charge model was, in contrast, optimized to give the most accurate possible partial atomic charges with our best set of density functionals, namely, the M06 suite.

In the past(12) we have not specifically recommended a particular combination of solvation model and level of theory or basis set because our methods are designed to be valid for many density functionals and basis sets. The M06-2X functional has been shown[13,15] to be the most generally reliable density functional for main-group chemistry, and the CM4M(11) charge model was developed to give the most accurate partial atomic charges when using the M06 family[13,15] of density functionals. Therefore we consider the combination of M06-2X and CM4 M to be our most reliable model. Even so, we examine both charge models, CM4 and CM4M, here. We do not consider different sets of radii but simply use the default SM8 radii which are the same as the SM6 radii for aqueous solutions. In fact, the SMx solvation models (x = 1,..., 6, 8) are designed in such a way that the nonbulk-electrostatic and nonelectrostatic contributions (which are together modeled by a term labeled CDS) are consistent with a given model (including radii) for bulk electrostatics, and one should never change the radii of a given SMx model except when creating specific reaction parameters for a specific application (which is in the category of expert usage, not routine usage).

Software

The SM8 solvation model is included in the GAMESSPLUS(30) and MN-GSM(31) software packages, and the SM6 solvation model (which is essentially identical for aqueous solution) is included in the GAMESSPLUS, MN-GSM, HONDOPLUS,(32) SMxGAUSS,(33) Jaguar,(34) and Q-Chem(35) software packages. The present calculations were carried out with MN-GSM.

Results and Discussion

In their continuum calculations, Nicholls et al.(1) considered a single geometry for each molecule; these geometries were optimized in the gas phase using the B3LYP[36-39] hybrid density functional and the 6-31G(d,p)(14) basis set. To examine the sensitivity of our results to geometry, we compare M06-2X solvation calculations with these B3LYP gas-phase structures to M06-2X solvation calculations with structures optimized, both in the gas phase and in aqeous solution, using the M06-2X hybrid density functional and the 6-31G(d) basis. These aqueous calculations employ CM4M. The results are given in Table 1 together with the values computed by Nicholls et al.(1) for their optimal implicit and explicit solvent models. The mean signed error (MSE), mean unsigned error (MUE), and RMSE values for the various models over the test set are also provided.

The RMSE of SM8 using M06-2X/6-31G(d) with B3LYP/6-31G(d,p) gas-phase optimized geometries over the test set is 1.25 kcal mol−1, which is substantially smaller than that from the Poisson approach using ZAP-9 radii, 1.88 kcal mol−1. The RMSE from explicit simulations is almost as good as that for these SM8 calculations, but the continuum calculations require several orders of magnitude less time computationally, so from a practical perspective a continuum calculation is typically a sensible first choice, especially since it is not only more economical but also more accurate (on average).

When the SM8 calculations are performed using M06-2X/6-31G(d) with M06-2X/6-31G(d,p) gas-phase or solution-phase optimized geometries, the predictions become even more accurate with a reduction of the RMSE to 1.14 and 1.08 kcal mol−1, respectively.

The largest errors from SM8 are associated with the two amides in the test set, N,N-dimethyl-p-methoxybenzamide and N,N-4-trimethylbenzamide. All of the models underestimate the solvation free energies for these two molecules, although the continuum models do considerably more poorly than the explicit solvation model, possibly reflecting specific solvation structure about the amide functional group that is not well-treated in the continuum approach. However relaxation of the amide geometry in solution causes a drastic improvement in the SM8 predictions. This observation is consistent with the observation by Nicholls et al.(26) that amides are much more polarized than other carbonyl-containing compounds. Relaxation of the solute geometry allows the amides to polarize further and significantly improves the predicted free energies of solvation.

Another functionality represented several times in the test set is the ester group, present in 1,1-diacetoxyethane, diethyl propanedioate, ethylene glycol diacetate, glycerol triacetate, and phenyl formate. The Poisson and explicit solvent models predict solvation free energies for these molecules that are too negative in every instance, with relatively large RMSEs of 2.40 and 4.48 kcal mol−1, respectively. SM8 also predicts too negative a solvation free energy in these three cases and has an RMSE over the five esters of 0.75 kcal mol−1.

It is noteworthy that all of the SM8 models using gas-phase optimized geometries have positive mean signed errors. Relaxation of molecular geometries reduces this systematic error, as shown in Table 1. (The explicit simulations of Nicholls et al. did use flexible solute molecules, so geometric relaxation in solution was considered in that case; however, Table 1 shows a large mean signed error.) While the SM8 model was originally parametrized using gas-phase geometries, it can be used with either gas-phase or solution-phase geometries. Table 1 shows that when we use geometries optimized in solution, the predicted free energy of solvation becomes more negative by 0.03−0.86 kcal mol−1. As a consequence, the MSE decreases from 0.35 to 0.03 kcal mol−1. However accounting for solute geometric relaxation does not introduce a significant improvement in the overall accuracy of the model, with the MUE increasing by 0.02 kcal mol−1 and the RMSE decreasing by only 0.06 kcal mol−1. We conclude that geometry optimization in solution is usually not worth the added expense.

As solvation free energies from continuum solvation calculations can be particularly sensitive to the treatment of electrostatics, we examined the performance of SM8 with two other hybrid density functionals, B3LYP and mPW1.(40) For these functionals, we employed the general CM4[11,16] charge model. We also considered the M06-2X functional with the CM4 (as opposed to CM4M) charge model. The results are presented in Table 2.

Table 2

Aqueous Solvation Free Energies (kcal mol−1) from SM8/CM4 with the 6-31G(d) Basis Set, B3LYP/6-31G(d,p) Gas-Phase Geometries, and Three Different Density Functionals, As Compared to Experiment

molecule	B3LYP	mPW1	M06-2X	expt
benzyl bromide	−2.32	−2.65	−2.34	−2.38
benzyl chloride	−2.14	−2.50	−2.12	−1.93
bis(2-chloroethyl) ether	−2.52	−2.61	−2.26	−4.23
1,1-diacetoxyethane	−5.76	−6.01	−5.16	−4.97
1,1-diethoxyethane	−1.50	−1.74	−1.33	−3.28
1,2-diethoxyethane	−2.23	−2.49	−1.94	−3.54
diethyl propanedioate	−5.45	−5.66	−4.77	−6.00
diethyl sulfide	−0.40	−0.52	−0.37	−1.43
dimethoxymethane	−1.94	−2.15	−1.61	−2.93
N,N-dimethyl-p-methoxybenzamide	−7.78	−8.25	−7.19	−11.01
1,4-dioxane	−4.75	−4.98	−4.34	−5.05
ethylene glycol diacetate	−6.95	−7.32	−6.33	−6.34
glycerol triacetate	−9.53	−9.93	−8.55	−8.84
imidazole	−8.46	−8.92	−8.50	−9.81
phenyl formate	−3.18	−3.55	−2.83	−3.82
m-bis(trifluoromethyl)benzene	1.26	1.14	1.68	1.07
N,N,4-trimethylbenzamide	−6.56	−7.01	−6.19	−9.76
MSE	0.83	0.54	1.18
MUE	1.10	1.00	1.23
RMSE	1.43	1.27	1.64

The computed solvation free energies using the CM4 charge model show some sensitivity to the choice of functional, with the RMSE ranging from 1.27 to 1.64 kcal mol−1. The smallest error is obtained with mPW1, consistent with this functional having been included in the original parametrization(16) of the CM4 model.

Further investigation into the sensitivity of our predictions to the choice of gas-phase geometry yields useful insight. We computed SM8 solvation free energies using mPW1/6-31G(d) and M06-2X/6-31G(d) optimized geometries. Those data are shown in Table 3. There is a consistent reduction in the errors of SM8 when the gas-phase geometries are optimized with M06-2X: compare a RMSE of 1.64 kcal mol−1 for M06-2X/CM4 at B3LYP geometries in Table 2 to 1.50 kcal mol−1 for M06-2X/CM4 using M06-2X geometries in Table 3.

Table 3

Aqueous Solvation Free Energies (kcal mol−1) from SM6 and SM8 with the 6-31G(d) Basis Set, with Two Different Charge Models, CM4 and CM4M, and with Gas-Phase Geometries Optimized with the Same Density Functional and Basis Set as the Solution Calculations, As Compared to Experiment

molecule	SM6 mPW1/CM4	SM8 mPW1/CM4	SM8 M06-2X/CM4	SM8 M06-2X/CM4M	expt
benzyl bromide	−2.44	−2.63	−2.28	−2.64	−2.38
benzyl chloride	−2.39	−2.48	−2.06	−2.48	−1.93
bis(2-chloroethyl) ether	−2.92	−2.71	−2.38	−2.72	−4.23
1,1-diacetoxyethane	−6.83	−6.74	−5.75	−6.71	−4.97
1,1-diethoxyethane	−2.17	−2.03	−1.57	−1.98	−3.28
1,2-diethoxyethane	−3.07	−2.84	−2.29	−2.83	−3.54
diethyl propanedioate	−6.03	−6.06	−5.03	−5.94	−6.00
diethyl sulfide	−0.58	−0.49	−0.31	−0.47	−1.43
dimethoxymethane	−3.01	−2.60	−2.07	−2.59	−2.93
N,N-dimethyl-p-methoxybenzamide	−8.27	−8.56	−7.50	−8.61	−11.01
1,4-dioxane	−5.94	−5.34	−4.69	−5.32	−5.05
ethylene glycol diacetate	−7.37	−7.28	−5.72	−6.60	−6.34
glycerol triacetate	−10.29	−10.17	−8.24	−9.55	−8.84
imidazole	−8.60	−9.08	−8.66	−9.34	−9.81
phenyl formate	−4.05	−4.02	−3.28	−4.14	−3.82
m-bis(trifluoromethyl)benzene	1.00	0.87	1.49	0.90	1.07
N,N,4-trimethylbenzamide	−6.65	−7.16	−6.30	−7.24	−9.76
MSE	0.27	0.29	1.04	0.35
MUE	1.00	0.95	1.14	0.86
RMSE	1.33	1.21	1.50	1.14

It is also useful to investigate the effect of choice of charge model, CM4 versus CM4M, upon the accuracy of our predictions. We find that there is a dramatic improvement when CM4M is used in combination with M06-2X as compared to using CM4; CM4M lowers the RMSE from 1.50 to 1.14 kcal mol−1. This result confirms our prior observation that there is often a tradeoff between using generally applicable parameters, as in CM4 that can be used with a variety of density functionals, and using a specifically optimized set of parameters, as in CM4M, which was optimized for the M06 suite of density functional.

Lastly it is useful to investigate how SM8 and its older and more widely available predecessor, SM6, compare on this test set. SM8 exhibits a small but definite improvement over SM6, reducing an RMSE of 1.33 kcal mol−1 for mPW1 with SM6 to 1.21 kcal mol−1 for mPW1 with SM8.

Conclusion

The most reliable implementation of the SM8 model predicts solvation free energies for a test set of 17 molecules with an RMSE of 1.08 kcal mol−1, which is as good as or better than other implicit and explicit solvation models previously reported by Nicholls et al.(1) The best error obtained by Nicholls et al.,(11) 1.33 kcal mol−1, was obtained in a competition of six solvation models. In our own trials we considered three different applications of SM8, namely, SM8 with M06-2X/6-31G(d) with B3LYP and M06-2X gas-phase optimized geometries and with solution-phase optimized geometries. The most accurate results were obtained using the solution-phase optimized geometries, which had an RMSE 1.08 kcal mol−1. The largest RMSE in the work of Nicholls et al. was 2.57 kcal mol−1, while the largest RMSE in the models considered here is 1.25 kcal mol−1. Then we further investigated the effect of the choice of geometry, the charge model, and the density functional upon the accuracy of the predictions; in these tests our best model based on gas-phase geometries had an RMSE of 1.14 kcal mol−1 using SM8 with M06-2X/6-31G(d) with CM4M and using M06-2X/6-31G(d) optimized gas-phase geometries. The worst model had an RMSE of 1.64 kcal mol−1 using SM8 with M06-2X/6-31G(d) with CM4 and using B3LYP/6-31+G(d,p) optimized gas-phase geometries. Nicholls et al.(1) emphasized in their conclusions that future improvements in current solvation modeling are critically dependent on an infusion of new experimental data for molecules characterized by high degrees of functionality (i.e., druglike). We are in complete agreement on this point and trust that the earnest pleas of the theoretical community will inspire experimental progress on this front.