Comparison of RESP and IPolQ-Mod Partial Charges for Solvation

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX-XXX

pubs.acs.org/JCTC

Comparison of RESP and IPolQ-Mod Partial Charges for Solvation Free Energy Calculations of Various Solute/Solvent Pairs Andreas Mecklenfeld†,‡ and Gabriele Raabe*,† †

Institut für Thermodynamik, Technische Universität Braunschweig, Hans-Sommer-Strasse 5, 38106 Braunschweig, Germany Center of Pharmaceutical Engineering, Technische Universität Braunschweig, Franz-Liszt-Strasse 35a, 38106 Braunschweig, Germany

‡

S Supporting Information *

ABSTRACT: The calculation of solvation free energies ΔGsolv by molecular simulations is of great interest as they are linked to other physical properties such as relative solubility, partition coefficient, and activity coefficient. However, shortcomings in molecular models can lead to ΔGsolv deviations from experimental data. Various studies have demonstrated the impact of partial charges on free energy results. Consequently, calculation methods for partial charges aimed at more accurate ΔGsolv predictions are the subject of various studies in the literature. Here we compare two methods to derive partial charges for the general AMBER force field (GAFF), i.e. the default RESP as well as the physically motivated IPolQ-Mod method that implicitly accounts for polarization costs. We study 29 solutes which include characteristic functional groups of drug-like molecules, and 12 diverse solvents were examined. In total, we consider 107 solute/solvent pairs including two water models TIP3P and TIP4P/2005. Comparison with experimental results yields better agreement for TIP3P, regardless of the partial charge scheme. The overall performance of GAFF/RESP and GAFF/IPolQ-Mod is similar, though specific shortcomings in the description of ΔGsolv for both RESP and IPolQ-Mod can be identified. However, the high correlation between free energies obtained with GAFF/RESP and GAFF/IPolQ-Mod demonstrates the compatibility between the modified charges and remaining GAFF parameters.

■

for ΔGsolv calculations with these models is undefined, various studies in the literature, e.g. within the SAMPL challenges,12−16 aimed at quantifying the predictability of solvation free energies. As a consequence, several suggestions for adjustments of molecular models were made. Rizzo et al.17 computed hydration free energies for various partial charge methods. Therein, empirical and semiempirical as well as ab initio schemes were compared for Poisson−Boltzmann and Generalized Born methods. Fennell et al.18 adjusted the description of the hydroxyl group from GAFF by fitting new Lennard-Jones parameters as well as partial charges to the liquid density and enthalpy of vaporization as well as the relative permittivity of methanol. They transferred the new parameters to other alcohols and received improved hydration free energy results compared to the default GAFF. Jämbeck et al.19 compared fitting procedures, i.e. CHELP,20 CHELPG,21 MK,22,23 and RESP,24 for GAFF and calculated partial charges at the B3LYP/ 6-311+G(d,p) level of theory for vacuum as well as for condensed phases. Furthermore, different water models were tested. Semiempirical as well as ab initio calculations for a

INTRODUCTION The Gibbs free energy of solvation characterizes the change of state for transferring a molecule from a vacuum into a condensed phase at constant temperature and pressure. The corresponding free energy difference describes the strength of solute−solvent interactions and is linked to a variety of physical properties such as relative solubility, partition coefficient, or activity coefficient. Free energy differences can be calculated, e.g. with molecular dynamics simulations, and statistical uncertainties are moderate, considering a set of guidelines.1−7 When configurational space overlap between states in vacuum and condensed phase is insufficient, a linking chain of nonphysical states is defined to ensure overlap. These intermediate states are defined by scaling the solute−solvent interactions with the alchemical variable λ, ranging from λ = 0 (no interactions, vacuum phase) to λ = 1 (full interactions, condensed phase). The number and spacing of λi-states influence both the accuracy and efficiency of a ΔG solv simulation.8 This makes most free energy simulations computationally demanding. Widely used all-atom additive force fields for the description of small organic compounds such as GAFF,9 CGenFF,10 or OPLS-AA11 do not at all or only marginally consider ΔGsolv in their parametrization approaches. Since the general suitability © XXXX American Chemical Society

Received: June 30, 2017

A

DOI: 10.1021/acs.jctc.7b00692 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

So far, to our best knowledge, studies in literature comparing RESP and IPolQ-Mod only focused on hydration free energy calculations. In this work we want to analyze the performance of RESP and IPolQ-Mod charges in combination with GAFF and consider solvation free energies for 29 solutes and 12 solvents. We thereby include the water models TIP3P and TIP4P/200542 to yield a total of 107 solute/solvent pairs with ΔGsolv ranging approximately from −46 to 10 kJ/mol. Our aim is thus to analyze the performance of partial charge methods for the prediction of solvation free energies for diverse solvents with a consistent approach to allow for its application in solvent screening. The paper is structured as follows: At first, we describe the systems studied, followed by an overview of the partial charge methods used. We then summarize simulation details before discussing our results. Therein, the performance of the water models is analyzed first, followed by the evaluation of the overall results. We continue with a clustered analysis of solutes and solvents, respectively, and finally present our conclusions.

variety of quantum mechanical levels of theory both with and without implicit solvent were studied by Mobley et al.25 for GAFF with TIP3P26 and TIP4P-Ew27 water models. Ahmed and Sandler28 studied hydration free energies obtained with GAFF, CGenFF, OPLS-AA, and TraPPE29 force fields with default as well as RESP, CHELPG, and CM430 partial charges. All these publications stressed the impact of partial charges on solvation free energies but simultaneously agreed that their description needs improvement. Additive force fields such as GAFF are incapable of describing instantaneous polarization effects caused by ambient molecules. As for solvation free energy calculations, polarization costs due to the change between condensed and vacuum phase cannot be directly considered. Polarizable force fields such as AMOEBA31 or Drude-201332 are expected to overcome this issue,33 though a recent comparison of AMOEBA and GAFF for the calculation of ΔGsolv by Mohamed et al.34 demonstrated a better performance of the additive force field. The authors argued that the widespread use of GAFF led to multiple rounds of parameter refinements, resulting in a better performance compared to the more detailed but apparently less well parametrized polarizable force field. Additionally, simulations of polarizable force fields are significantly more computationally demanding. The application of additive force fields such as GAFF is therefore still relevant. GAFF RESP charges are calculated for the isolated molecule. Thereby the 6-31G* basis set is used as it overestimates the dipole moment to yield charges that approximate the effect of mean polarization in the condensed phase.35 However, the actual polarity of the solvent is not taken into account. Deriving fixed charges to capture the phase transfer of a solute with accompanying polarization costs is still the subject of research. Cerutti et al.36 presented with the IPolQ approach a new method to implicitly account for polarization costs considering an explicitly represented solvent. Ab initio MP2 calculations are performed with the cc-pV(T+d)Z basis set, and partial charges are fitted by applying the RESP method. The reaction field of the condensed phase is obtained by sampling configurations of TIP4P-Ew water molecules with adjusted partial charges around the “frozen” solute, yielding a field of point charges for which solute charges are deduced. Following the approach of Karamertzanis et al.,37 charges are calculated both for the unpolarized vacuum state and polarized condensed state. The charge sets are averaged to account for the electronic polarization of the solute molecule. The authors further implemented an iteration cycle to update Lennard-Jones parameters of certain atom types during the derivation of charges for a better description of ΔGsolv. By additionally compiling newly derived torsion parameters, Cerutti et al. published their adjusted force field as ff14ipq.38 However, due to the poor description of salt bridge interactions, parameters for partial charges, angle, and torsion potentials as well as atomic radii were again refined for the ff15ipq39 force field, but solvation free energies were no longer targeted. Parameter sets are contained in the AmberTools distribution from version 16.40 The IPolQ approach was adopted by Muddana et al., who simplified the explicit representation of the solvent by an implicit continuum model, referring to this method as IPolQMod.41 IPolQ-Mod derived charges were combined with GAFF parameters within the SAMPL4 challenge, and 47 hydration free energies were predicted with comparable accuracy as with the default RESP charges.

■

METHODS Choice of Systems. In this study we compared the performance of GAFF/RESP and GAFF/IPolQ-Mod for ΔGsolv predictions of 107 solute/solvent pairs, listed in Figure 1. The choice of solutes aims at including model compounds with prominent atom types in functional groups characteristic for potential drug candidates. We thereby covered compounds of different types, i.e. halogenated benzenes (chloro-, bromo-,

Figure 1. List of systems studied. Solvents are marked in bold, and the corresponding solutes are listed below. B


Article

Journal of Chemical Theory and Computation iodo-), alkanes (methane to n-hexane), alcohols (methanol, ethanol, phenol), ethers (methyl phenyl and ethyl phenyl), aldehydes and ketone (acetaldehyde, formaldehyde, acetone), nitriles (aceto- and benzo-), and aromatics (benzene, toluene) as well as benzamide, methylamine, and 3-methylindole as examples of amide, amine, and indole compounds. Solvents included chloroform, cyclohexane, ethanol, diethyl ether, benzene, toluene, benzonitrile, and water. Publications like those of Hess and van der Vegt,43 Shirts and Pande,44 or Jämbeck et al.19 demonstrated the effect of different water models on ΔGsolv. For simplification, however, we only studied two water models. We chose TIP3P, as it was used in the development of the AMBER35 force field and performed well in the studies previously mentioned. Furthermore, we included the more recent TIP4P/2005, which is regarded to be a general purpose model.42 For self-solvations, we additionally considered acetonitrile, ethyl phenyl ether, and methyl phenyl ether. Our choice of systems for this study is though restricted by the availability of experimental data for comparison. For most systems, experimental data were obtained from the Minnesota Solvation Database.45 For the solute 3-methylindole, experimental data were taken from the work of Villa and Mark.46 Calculation of Partial Charges. Prior to the calculation of partial charges, the molecule geometries are optimized at the HF/6-31G* level of theory. According to the RESP procedure, the electrostatic potential is calculated with the HF method and the 6-31G* basis set in vacuum. Partial charges are fitted according to the two-stage RESP scheme.24 For IPolQ-Mod charges, ab initio simulations are performed at the MP2/aug-cc-pVDZ level of theory. To account for polarization in pure fluid or solvent phase, the polarizable continuum model by Mennucci et al.47 is applied for the respective solvent (Gaussian keyword: SCRF=PCM). Solute charges for vacuum and condensed phase are fitted according to the RESP method. For an implicit representation of polarization costs, solute charges from both phases are then averaged to yield the IPolQ-Mod charges as illustrated in Figure 2. Furthermore, we employed the polarizable continuum model to account for the self-polarization in partial charges for solvent models. Merely the water models TIP3P and TIP4P/2005 remained unchanged.

In this study, all ab initio calculations were performed with Gaussian09,48 and the RESP-fitting was applied as implemented in Antechamber from AmberTools16.40 However, basis sets in Gaussian are only present for atoms not exceeding period number 4 of the periodic table, excluding the element iodine. For the parametrization of the compound iodobenzene we therefore employed a pseudopotential aug-cc-pVDZ-PP basis set49 obtained from the EMSL Basis Set Library.50,51 To our best knowledge, a comparable 6-31G* pseudopotential does not exist so that free energy calculations for iodobenzene were only performed with GAFF/IPolQ-Mod. For the comparison of averaged deviations, i.e. mean unsigned errors (MUE) and root-mean-square deviations (RMSD) between simulation results and experimental data, iodobenzene was excluded from the analysis to consider identical data sets. Simulation Details. Molecular dynamics simulations were performed using GROMACS 2016.1.52−59 GROMACS compatible GAFF topology files were created using Antechamber from AmberTools1640 and ACPYPE.60 Partial charges were calculated as previously described. All molecular dynamics simulations were performed using a Langevin dynamics integrator61 with a time step of 0.5 fs. No constraints were applied except for the rigid water models, in which case we used the SETTLE62 algorithm. During the equilibration, we employed a Berendsen63 barostat but changed to Parrinello−Rahman64 for the subsequent simulations. The pressure was set to p = 1 atm using a time constant of τp = 5 ps. The temperature was maintained by the integrator for an inverse friction constant of τT = 2.5 ps. Experimental data are tabulated at a temperature of T = 298 K in the Minnesota Solvation Database, while Villa and Mark provide results for 3methylindole at T = 293 K, and temperature target values were set accordingly. Electrostatic interactions were evaluated using Particle-Mesh Ewald (PME)65 of interpolation order 4 with a real space cutoff rcoulomb = 1.2 nm and a Fourier spacing of 0.12 nm. Lennard-Jones interactions behind a cutoff rvdw = 1.2 nm were considered with a long-range correction for energy and pressure. Solvation free energies were calculated at infinite dilution, i.e. a single solute molecule was surrounded by a sufficient number of solvent molecules while applying periodic boundary conditions. According to the good practice method proposed by Parameswaran and Mobley, a simulation box should exceed the solute size plus the doubled Lennard-Jones cutoff radius.5 We tested several box dimensions and found this criterion to be sufficient. However, for larger solvent compounds, i.e. all but acetonitrile, chloroform, ethanol, and water, less than 200 solvent molecules would be needed according to this scheme. To improve the statistics in the description of the solvent, we used at least 250 solvent molecules instead. For initial configurations, molecules were randomly distributed using PACKMOL.66 Free energies were evaluated with the Multistate Bennett Acceptance Ratio (MBAR)67 method as implemented in the ‘alchemicalanalysis.py’ tool.2 For the change of state between vacuum and condensed phase, intermediate states were defined by scaling solute−solvent interactions with the alchemical variable λ. Based on our previous study,8 we applied a soft-core potential as proposed by Beutler et al.68 for decoupling Lennard-Jones interactions according to

Figure 2. Illustration of the calculation for IPolQ-Mod partial charges with the example of chlorobenzene for solvent water. Partial charges are calculated both in vacuum and condensed phase at the MP2/augcc-pVDZ level of theory. The average of the two charge sets yields partial charges according to the IPolQ-Mod method. For a clearer representation, only four decimals of partial charge values with unit e are displayed. C


Article

Journal of Chemical Theory and Computation UijLJ(rij ,

⎡⎛

⎞12/ c 1 ⎟ λ) = 4εijλ ⎢⎜ b c⎟ α (1 − λ ) + ( r / σ ) ij ij ⎠ ⎢⎣⎝ a ⎢⎜

⎛ ⎞6/ c ⎤ 1 ⎜ ⎟ ⎥ −⎜ b c⎟ ⎥ α (1 λ ) ( r / σ ) − + ⎝ ij ij ⎠ ⎥ ⎦

(1)

with parameters a, b, c, and α set to 1, 1, 48, and 0.003, respectively. For each λi-state, an energy minimization was performed, followed by a short equilibration of 400,000 steps. Guided by the procedure presented in our recent paper,8 we ran between 3 and 5 TRIAL simulations of 800,000 steps, respectively, for an adaptive optimization of λi-state spacing, thus extending the actual equilibration phase. A suitable number of λi-states was determined by evaluating the configurational space overlap matrix dumped by the ‘alchemical-analysis.py’ script after each TRIAL run. Based on an empirically derived algorithm,8 poor overlap resulted in an increase of λi-states, while “excessive” led to its decrease. Regarding the systems studied, the number of λi-states ranges between 6 and 23. That is, after the optimization procedure, the number and spacing of λi-states were defined individually, and all systems were carefully inspected to ensure sufficient configurational space overlaps as well as a long enough equilibration phase. This was followed by five successive blocks of production runs with 4,000,000 steps each, yielding a total of 20,000,000 steps per λi-state for production. Hamiltonians were calculated between all λi-states and stored every 100 steps. The five blocks were evaluated independently, and only uncorrelated data were considered to calculate block results ΔGsolv,i. Final results ΔGsolv are given as the average of Nb = 5 block results ΔGsolv,i with the standard error of the mean SEM

Figure 3. Comparison of calculated solvation free energies with GAFF/RESP (red squares) and GAFF/IPolQ-Mod (blue circles) for water models TIP3P (filled) and TIP4P/2005 (hollow) to experimental data. For clarity, error bars are not shown. The diagonal solid line indicates perfect agreement between calculated and experimental results, and the gray area marks deviations within ±5 kJ/mol. For hydration free energies, data points are labeled by the solute name only.

while there is noticeable scattering for most of the other systems. For GAFF/IPolQ-Mod, significant deviations both to experiment and GAFF/RESP occur for the solutes acetonitrile, benzonitrile, and acetone. On the other hand, GAFF/RESP yields poor hydration free energy results for benzamide. In contrast to the other alcohols methanol and ethanol, the hydration of phenol is poorly described regardless of the partial charge method and water model used. For the analysis of deviations from experimental data, results for the solute iodobenzene as well as the self-solvations of TIP3P and TIP4P/2005 were dismissed. The mean unsigned errors and root-mean-square deviations are given in Table 1.

N

∑i =b1 (ΔGsolv, i − ΔGsolv )2

SEM =

Nb(Nb − 1)

(2)

The standard errors of the mean range from 0.02 to 0.34 kJ/ mol. All simulation results are presented in the SI. To evaluate the performance of GAFF/RESP and GAFF/IPolQ-Mod we determined the mean unsigned errors (MUE) and root-meansquare deviations (RMSD) from experimental data by MUE =

1 N

N experiment simulation − ΔGsolv, | ∑ |ΔGsolv, k k k

Table 1. Comparison of GAFF/RESP and GAFF/IPolQMod Regarding the MUE and RMSD Values in kJ/mol for Hydration Free Energies as well as the Solvation Free Energies of Toluene in Water for Models TIP3P and TIP4P/ 2005

(3)

and RMSD =

■

1 N

N experiment 2 simulation − ΔGsolv, ) ∑ (ΔGsolv, k k k

(4)

RESULTS AND DISCUSSION Comparison of Water Models. The comparison of ΔGsolv results for both water models with GAFF/RESP and GAFF/ IPolQ-Mod is shown in Figure 3. Water served as solvent for the halogenated benzenes, alkanes, alcohols, ethers, nitriles, acetone, benzamide, methylaminde, and 3-methylindole. Additionally, we analyzed water as solute in toluene and also determined the self-solvations of TIP3P and TIP4P/2005, respectively. Results for alkanes (all with positive ΔGsolv) demonstrate very little variations in ΔGsolv for the different charge schemes,

MUE RMSD

GAFF/ RESP TIP3P

GAFF/RESP TIP4P/2005

GAFF/IPolQMod TIP3P

GAFF/IPolQMod TIP4P/ 2005

2.48 3.25

3.14 3.88

2.69 3.98

3.40 4.60

The most accurate results were obtained for GAFF with RESP charges. The differences of the mean unsigned errors between RESP and IPolQ-Mod calculations are minor, while the impact of the water models on the deviation to experimental data is more relevant. The large RMSD values for both water models and GAFF/IPolQ-Mod results refer to the high deviations of acetonitrile, benzonitrile, and acetone. D


Article


solvation (= all studied solvents but cyclohexane) are remarkably worse for GAFF/IPolQ-Mod. The analysis of deviations from the experimental data is summarized in Table 2. For data set consistency, free energy results for the solute

The data reveal a better performance of TIP3P over TIP4P/ 2005 for both GAFF/RESP and GAFF/IPolQ-Mod models. Therefore, we discard the TIP4P/2005 model from further analysis. Evaluation of the Overall Performance. In this section, we evaluate the general performance of GAFF/RESP and GAFF/IPolQ-Mod, respectively, for the remaining 85 solute/ solvent pairs after the exclusion of TIP4P/2005. Simulation results with reference to experimental data are displayed in Figure 4.

Table 2. Comparison of Averaged Mean Unsigned Errors (MUE) and Root-Mean-Square Deviations (RMSD) in kJ/ mol for ΔGsolv Results with GAFF/RESP and GAFF/IPolQMod MUE RMSD

GAFF/RESP

GAFF/IPolQ-Mod

2.64 3.53

3.08 4.29

iodobenzene were not taken into account. Table 2 confirms the qualitative inspection of Figure 4. While the difference in mean unsigned error is small, the root-mean-square deviation is significantly higher for GAFF/IPolQ-Mod due to the larger number of outliers. Compared to the deviations of hydration free energy results listed in Table 1, the overall performance for ΔGsolv in different solvents is slightly worse for both GAFF/RESP and GAFF/ IPolQ-Mod. Differences in ΔGsolv by Solutes. We have identified not only acetonitrile and benzonitrile but also acetaldehyde, acetone, and formaldehyde, in combination with the more polar solvents studied here, to be poorly described with GAFF/ IPolQ-Mod. However, other solute/solvent pairs yield high ΔGsolv deviations as well, also for GAFF/RESP. We analyzed the accuracies of ΔGsolv results for one solute or a group of similar solutes in various solvents. Alkanes, acetaldehyde, acetone, and formaldehyde as well as benzene and toluene were clustered, respectively, due to their structural similarities. Furthermore, we studied the “self-solvation” when experimental data are available for comparison. For all solute/solvent pairs in a category, we determined the mean unsigned errors and calculated the sample standard deviations. The comparison of the mean unsigned errors in Figure 5 confirms our analysis regarding the description of nitriles by GAFF/IPolQ-Mod. In comparison to GAFF/RESP, the averaged MUE of acetonitrile and benzonitrile with GAFF/ IPolQ-Mod is significantly higher by 4.34 and 3.48 kJ/mol, respectively. Also, the description of the solutes acetaldehyde, acetone, and formaldehyde displayed in the category “aldehydes and ketone” is poor compared to GAFF/RESP. On the other hand, the solvation free energies of the solute benzamide are remarkably worse for GAFF/RESP than for GAFF/IPolQ-Mod with a difference in the MUE value of 3.63 kJ/mol. Furthermore, ΔGsolv results for the solute phenol show poor accuracies both with RESP and IPolQ-Mod partial charges, while calculations for methanol and ethanol are well below the total MUE values indicated by the horizontal lines in Figure 5. These compounds with MUE > 5 kJ/mol (which are not within the gray area in Figure 5) are considered as outliners. When these outliners, i.e. the solutes phenol, acetaldehyde, acetone, formaldehyde, benzamide, acetonitrile, and benzonitrile, are excluded from the analysis, the accuracies of GAFF/RESP and GAFF/IPolQ-Mod are nearly identical. MUE and RMSD values for the remaining 53 out of the former 85 solute/solvent pairs are summarized in Table 3. In Figure 5, results for benzene and toluene as well as 3methylindole are above the corresponding total MUE values as well, though the deviations are smaller. The remaining

Figure 4. Calculated solvation free energy results for GAFF/RESP (red squares) and GAFF/IPolQ-Mod (blue circles) for all compounds except TIP4P/2005 vs experimental data. For clarity, error bars are not shown. Solute/solvent systems with ΔGsolv deviations higher than ±5 kJ/mol compared to experimental data (gray area) are labeled, except for solvent TIP3P, which we previously discussed for Figure 3. The diagonal solid line indicates perfect agreement between simulation results and experimental data.

The plot shows broadly good agreement between calculated ΔGsolv results for GAFF/RESP and GAFF/IPolQ-Mod. To quantify the linear correlation between free energy results with RESP and IPolQ-Mod charges, we calculated the Pearson correlation value and obtained R = 0.9860. While a Pearson value of 1 indicates a linear relation, this still demonstrates a high correlation between free energy results obtained from GAFF/RESP and GAFF/IPolQ-Mod. Furthermore, most results agree well with experimental data, though remarkable outliers are present for GAFF/RESP as well as GAFF/IPolQ-Mod. Deviations to experimental data for both GAFF/RESP and GAFF/IPolQ-Mod occur, i.e. for the solute/ solvent pairs chlorobenzene/ethanol, phenol/TIP3P, phenol/ diethyl ether, benzamide/cyclohexane, and benzamide/diethyl ether. As discussed before, ΔGsolv of benzamid in TIP3P is poorly described by GAFF/RESP. For IPolQ-Mod, solutes acetone and acetaldehyde in solvent chloroform as well as acetone in TIP3P show significant deviations to experiment and GAFF/RESP. Additionally, the ΔGsolv calculations of both solutes acetonitrile and benzonitrile in solvents chloroform, diethyl ether, and TIP3P as well as the corresponding selfE


Article


Figure 6. Averaged mean unsigned errors for solvents in combination with various solutes. Error bars indicate sample standard deviations. Dashed red and point-dashed blue horizontal lines mark total mean unsigned errors for GAFF/RESP and GAFF/IPolQ-Mod, respectively, according to Table 2.

Figure 5. Mean unsigned errors for groups of solutes with various solvents each. Error bars indicate sample standard deviations. Dashed red and point-dashed blue horizontal lines mark total mean unsigned errors for GAFF/RESP and GAFF/IPolQ-Mod, respectively, according to Table 2. Outliers are specified as compounds with MUE > 5 kJ/ mol, thus exceeding the gray marked area.

partial charges can be mostly attributed to the description of the solutes.

■

Table 3. Comparison of Averaged Mean Unsigned Errors (MUE) and Root-Mean-Square Deviations (RMSD) in kJ/ mol for ΔGsolv Results with GAFF/RESP and GAFF/IPolQMod but with Excluded Outliersa MUE RMSD

GAFF/RESP

GAFF/IPolQ-Mod

1.77 2.12

1.73 2.17

CONCLUSIONS We compared the accuracy of solvation free energy predictions for 107 solvent/solute pairs using GAFF with partial charges obtained with the RESP and IPolQ-Mod scheme. We demonstrated that the results with TIP3P are more accurate than those with TIP4P/2005 regardless of the protocol used to calculate partial charges. Therefore, we ignored the TIP4P/ 2005 model in later analysis. Further evaluation highlighted the generally good agreement between ΔGsolv results from GAFF/ RESP and GAFF/IPolQ-Mod with a Pearson value of R = 0.9860. However, poor ΔGsolv results for certain solute/solvent pairs have been identified for both GAFF/RESP and GAFF/ IPolQ-Mod, respectively. While GAFF/RESP shows poor performance especially for the solute benzamide, GAFF/ IPolQ-Mod yields poor solvation free energy predictions for the nitrile solutes and, to a minor degree, also for the aldehydes and ketone compounds studied here. When outliers are excluded, GAFF/RESP and GAFF/IPolQ-Mod yield results of comparable accuracies, though only 53 solute/solvent pairs remain for the comparison. We then studied identical or at least similar compounds in combination with different solvents or solutes, respectively, and calculated the mean unsigned errors for these groups. This revealed specific shortcoming for both RESP and IPolQ-Mod in the description of solutes rather than solvents. That is, we have identified solute compounds for which GAFF/IPolQ-Mod performs significantly worse than GAFF/RESP but also vice versa. We want to emphasize that the number of solutes and solvents studied here is limited due to the lack of experimental data available for comparison. Considering additional functional groups, a larger number or more complex molecules might relativize or disprove our findings. While the overall accuracy of GAFF/IPolQ-Mod for the systems studied here falls behind the

a

I.e. the solutes phenol, acetaldehyde, acetone, formaldehyde, benzamide, acetonitrile, and benzonitrile.

compounds, including the self-solvations, are below the total MUE values and thus show good agreement with experimental data. The range of sample standard deviations indicates varying ΔGsolv accuracies depending on the solvents. Differences in ΔGsolv by Solvents. Besides the solutes, we furthermore studied the performance of partial charge schemes with regard to the solvents benzene, toluene, TIP3P, cyclohexane, diethyl ether, and chloroform. Mean unsigned errors for these solvents are displayed in Figure 6. Results for iodobenzene were again excluded for data set consistency. The figure shows the highest MUE value for chloroform and GAFF/IPolQ-Mod, and the deviation to GAFF/RESP is higher by 1.13 kJ/mol. The mean unsigned error for the solvent diethyl ether is also above the total MUE value (point dashed blue horizontal line) as well. For GAFF/ RESP, the highest mean unsigned error is yielded by solvent diethyl ether, while the MUE values of chloroform and cyclohexane are only slightly above the total average deviation (dashed red horizontal line). The remaining MUE values from GAFF/RESP and GAFF/IPolQ-Mod are below the total MUE values. Besides chloroform, accuracies for GAFF/RESP and GAFF/IPolQ-Mod are similar for these solvents, thus the discrepancies in ΔGsolv accuracies due to RESP and IPolQ-Mod F



■

RESP results, one should consider the following aspects to rate the results of IPolQ-Mod charges for solvation free energy calculations: 1. GAFF/RESP is a self-consistent force field and the application of IPolQ-Mod charges constitutes a perturbation of the self-consistency. 2. GAFF was developed to be compatible with HF/6-31G* RESP charges. The initial application of the 6-31G* basis set in the AMBER35 force field and the successive GAFF was motivated by the good representation of hydration free energies for various aromatics69 as well as different alchemical transformations in aqueous solutions.70 A shifting focus toward other solvents, e.g. for the prediction of relative solubilities, calls the use of this basis set into question. 3. IPolQ-Mod charges entail implicit polarization costs with respect to specific solvents and represent a physically motivated approach aimed at solvation free energy calculations. 4. Acetonitrile and benzonitrile as well as acetaldehyde, acetone, and formaldehyde in polar solvents caused the poor results for GAFF/IPolQ-Mod. The number of these solute/ solvent pairs is limited, and their exclusion yields a nearly identical performance compared to GAFF/RESP. 5. GAFF/RESP and GAFF/IPolQ-Mod showed a nearly linear correlation with a Pearson value of R = 0.9860, implying a general compatibility between IPolQ-Mod charges and remaining GAFF parameters. That is, although the general performance of GAFF is good, solvation free energies were not specifically targeted in the force field parametrization. IPolQ-Mod charges might be better suited for a broader set of solvents, though force field adjustments seem necessary, particularly for the nitrile and carbonyl groups. Overall, GAFF seems to be a sufficient basis for further force field adaption such as refitting Lennard-Jones parameters compatible with IPolQ-Mod partial charges for an improved description of solvation free energies in various solvents.

■

REFERENCES

(1) Hansen, N.; van Gunsteren, W. F. Practical Aspects of FreeEnergy Calculations: A Review. J. Chem. Theory Comput. 2014, 10, 2632−2647. (2) Klimovich, P. V.; Shirts, M. R.; Mobley, D. L. Guidelines for the analysis of free energy calculations. J. Comput.-Aided Mol. Des. 2015, 29, 397−411. (3) Bruckner, S.; Boresch, S. Efficiency of alchemical free energy simulations. I. A practical comparison of the exponential formula, thermodynamic integration, and Bennett’s acceptance ratio method. J. Comput. Chem. 2011, 32, 1303−1319. (4) Bruckner, S.; Boresch, S. Efficiency of alchemical free energy simulations. II. Improvements for thermodynamic integration. J. Comput. Chem. 2011, 32, 1320−1333. (5) Parameswaran, S.; Mobley, D. L. Box size effects are negligible for solvation free energies of neutral solutes. J. Comput.-Aided Mol. Des. 2014, 28, 825−829. (6) Naden, L. N.; Pham, T. T.; Shirts, M. R. Linear Basis Function Approach to Efficient Alchemical Free Energy Calculations. 1. Removal of Uncharged Atomic Sites. J. Chem. Theory Comput. 2014, 10, 1128−1149. (7) Naden, L. N.; Shirts, M. R. Linear basis function approach to efficient alchemical free energy calculations. 2. Inserting and deleting particles with coulombic interactions. J. Chem. Theory Comput. 2015, 11, 2536−2549. (8) Mecklenfeld, A.; Raabe, G. Efficient solvation free energy simulations: Impact of soft-core potential and a new adaptive λ-spacing method. Mol. Phys. 2017, 115, 1322−1334. (9) Wang, J.; 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. (10) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671−690. (11) Price, M. L. P.; Ostrovsky, D.; Jorgensen, W. L. Gas-phase and liquid-state properties of esters, nitriles, and nitro compounds with the OPLS-AA force field. J. Comput. Chem. 2001, 22, 1340−1352. (12) Nicholls, A.; Mobley, D. L.; Guthrie, J. P.; Chodera, J. D.; Bayly, C. I.; Cooper, M. D.; Pande, V. S. Predicting small-molecule solvation free energies: An informal blind test for computational chemistry. J. Med. Chem. 2008, 51, 769−779. (13) Guthrie, J. P. A blind challenge for computational solvation free energies: Introduction and overview. J. Phys. Chem. B 2009, 113, 4501−4507. (14) Geballe, M. T.; Skillman, A. G.; Nicholls, A.; Guthrie, J. P.; Taylor, P. J. The SAMPL2 blind prediction challenge: introduction and overview. J. Comput.-Aided Mol. Des. 2010, 24, 259−279. (15) Mobley, D. L.; Wymer, K. L.; Lim, N. M.; Guthrie, J. P. Blind prediction of solvation free energies from the SAMPL4 challenge. J. Comput.-Aided Mol. Des. 2014, 28, 135−150. (16) Muddana, H. S.; Fenley, A. T.; Mobley, D. L.; Gilson, M. K. The SAMPL4 host-guest blind prediction challenge: An overview. J. Comput.-Aided Mol. Des. 2014, 28, 305−317. (17) Rizzo, R. C.; Aynechi, T.; Case, D. A.; Kuntz, I. D. Estimation of Absolute Free Energies of Hydration Using Continuum Methods: Accuracy of Partial Charge Models and Optimization of Nonpolar Contributions. J. Chem. Theory Comput. 2006, 2, 128−139. (18) Fennell, C. J.; Wymer, K. L.; Mobley, D. L. A fixed-charge model for alcohol polarization in the condensed phase, and its role in small molecule hydration. J. Phys. Chem. B 2014, 118, 6438−6446. (19) Jambeck, J. P. M.; Mocci, F.; Lyubartsev, A. P.; Laaksonen, A. Partial atomic charges and their impact on the free energy of solvation. J. Comput. Chem. 2013, 34, 187−197. (20) Chirlian, L. E.; Francl, M. M. Atomic charges derived from electrostatic potentials: A detailed study. J. Comput. Chem. 1987, 8, 894−905.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00692.

■

Article

Calculated solvation free energy results (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Andreas Mecklenfeld: 0000-0003-0790-8057 Funding

Andreas Mecklenfeld acknowledges funding by a GeorgChristoph-Lichtenberg fellowship by the federal state of Lower Saxony, Germany. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS All simulations were performed on a cluster of North-German Supercomputing Alliance (HLRN). We greatly appreciate the support. G


Article


H.; Gotz, A. W. et al. AMBER; University of California: San Francisco, CA, 2015. (41) Muddana, H. S.; Sapra, N. V.; Fenley, A. T.; Gilson, M. K. The SAMPL4 hydration challenge: evaluation of partial charge sets with explicit-water molecular dynamics simulations. J. Comput.-Aided Mol. Des. 2014, 28, 277−287. (42) Abascal, J. L. F.; Vega, C. A general purpose model for the condensed phases of water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (43) Hess, B.; van der Vegt, N. F. A. Hydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models. J. Phys. Chem. B 2006, 110, 17616−17626. (44) 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. (45) Marenich, A. V.; Kelly, C. P.; Thompson, J. D.; Hawkins, G. D.; Chambers, C. C.; Giesen, D. J.; Winget, P.; Cramer, C. J.; Truhlar, D. G. Minnesota Solvation Database − version 2012; University of Minnesota: Minneapolis, 2012. (46) Villa, A.; Mark, A. E. Calculation of the free energy of solvation for neutral analogs of amino acid side chains. J. Comput. Chem. 2002, 23, 548−553. (47) Mennucci, B.; Cammi, R.; Tomasi, J. Excited states and solvatochromic shifts within a nonequilibrium solvation approach: A new formulation of the integral equation formalism method at the selfconsistent field, configuration interaction, and multiconfiguration selfconsistent field level. J. Chem. Phys. 1998, 109, 2798−2807. (48) 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. et al. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (49) Peterson, K. A.; Shepler, B. C.; Figgen, D.; Stoll, H. On the spectroscopic and thermochemical properties of ClO, BrO, IO, and their anions. J. Phys. Chem. A 2006, 110, 13877−13883. (50) Feller, D. The role of databases in support of computational chemistry calculations. J. Comput. Chem. 1996, 17, 1571−1586. (51) Schuchardt, K. L.; Didier, B. T.; Elsethagen, T.; Sun, L.; Gurumoorthi, V.; Chase, J.; Li, J.; Windus, T. L. Basis set exchange: A community database for computational sciences. J. Chem. Inf. Model. 2007, 47, 1045−1052. (52) Bekker, H.; Berendsen, H. J. C.; Dijkstra, E. J.; Achterop, S.; van Drunen, R.; van der Spoel, D.; Sijbers, A.; Keegstra, H.; Reitsma, B.; Renardus, M. K. R. Gromacs: A parallel computer for molecular dynamics simulations. PHYSICS COMPUTING '92; de Groot, R. A., Nadrchal, J., Eds.; World Scientific Publishing: Singapore, 1993; pp 252−256. (53) Berendsen, H.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 43−56. (54) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: A package for molecular simulation and trajectory analysis. J. Mol. Model. 2001, 7, 306−317. (55) 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. (56) 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. (57) 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.; et al. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845−854. (58) Páll, S.; Abraham, M. J.; Kutzner, C.; Hess, B.; Lindahl, E. Tackling Exascale Software Challenges in Molecular Dynamics Simulations with GROMACS. In Solving Software Challenges for Exascale; Markidis, S., Laure, E., Eds.; Lecture Notes in Computer Science; Springer International Publishing: Cham, 2015; pp 3−27, DOI: 10.1007/978-3-319-15976-8_1.

(21) Breneman, C. M.; Wiberg, K. B. Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis. J. Comput. Chem. 1990, 11, 361−373. (22) Singh, U. C.; Kollman, P. A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 1984, 5, 129− 145. (23) Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic charges derived from semiempirical methods. J. Comput. Chem. 1990, 11, 431− 439. (24) 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. (25) Mobley, D. L.; Dumont, E.; Chodera, J. D.; Dill, K. A. Comparison of charge models for fixed-charge force fields: smallmolecule hydration free energies in explicit solvent. J. Phys. Chem. B 2007, 111, 2242−2254. (26) 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. (27) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665−9678. (28) Ahmed, A.; Sandler, S. I. Hydration Free Energies of Multifunctional Nitroaromatic Compounds. J. Chem. Theory Comput. 2013, 9, 2774−2785. (29) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (30) 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. (31) Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P. The Polarizable Atomic Multipole-based AMOEBA Force Field for Proteins. J. Chem. Theory Comput. 2013, 9, 4046−4063. (32) Lopes, P. E. M.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell, A. D., JR. Force Field for Peptides and Proteins based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2013, 9, 5430−5449. (33) Leontyev, I. V.; Stuchebrukhov, A. A. Electronic Polarizability and the Effective Pair Potentials of Water. J. Chem. Theory Comput. 2010, 6, 3153−3161. (34) Mohamed, N. A.; Bradshaw, R. T.; Essex, J. W. Evaluation of solvation free energies for small molecules with the AMOEBA polarizable force field. J. Comput. Chem. 2016, 37, 2749−2758. (35) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (36) Cerutti, D. S.; Rice, J. E.; Swope, W. C.; Case, D. A. Derivation of fixed partial charges for amino acids accommodating a specific water model and implicit polarization. J. Phys. Chem. B 2013, 117, 2328− 2338. (37) Karamertzanis, P. G.; Raiteri, P.; Galindo, A. The Use of Anisotropic Potentials in Modeling Water and Free Energies of Hydration. J. Chem. Theory Comput. 2010, 6, 1590−1607. (38) Cerutti, D. S.; Swope, W. C.; Rice, J. E.; Case, D. A. ff14ipq: A Self-Consistent Force Field for Condensed-Phase Simulations of Proteins. J. Chem. Theory Comput. 2014, 10, 4515−4534. (39) Debiec, K. T.; Cerutti, D. S.; Baker, L. R.; Gronenborn, A. M.; Case, D. A.; Chong, L. T. Further along the Road Less Traveled: AMBER ff15ipq, an Original Protein Force Field Built on a SelfConsistent Physical Model. J. Chem. Theory Comput. 2016, 12, 3926. (40) 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


Article

Journal of Chemical Theory and Computation (59) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1−2, 19−25. (60) Sousa da Silva, A. W.; Vranken, W. F. ACPYPE - AnteChamber PYthon Parser interfacE. BMC Res. Notes 2012, 5, 367. (61) Goga, N.; Rzepiela, A. J.; de Vries, A. H.; Marrink, S. J.; Berendsen, H. J. C. Efficient Algorithms for Langevin and DPD Dynamics. J. Chem. Theory Comput. 2012, 8, 3637−3649. (62) Miyamoto, S.; Kollman, P. A. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 1992, 13, 952−962. (63) 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. (64) Parrinello, M.; Rahman, A. Crystal Structure and Pair Potentials: A Molecular-Dynamics Study. Phys. Rev. Lett. 1980, 45, 1196−1199. (65) 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. (66) Martinez, L.; Andrade, R.; Birgin, E. G.; Martinez, J. M. PACKMOL: a package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 2157−2164. (67) Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. J. Chem. Phys. 2008, 129, 124105. (68) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222, 529−539. (69) Kuyper, L. F.; Hunter, R. N.; Ashton, D. Free energy calculations on the relative solvation free energies of benzene, anisole, and 1,2,3-trimethoxybenzene: Theoretical and experimental analysis of aromatic methoxy solvation. J. Phys. Chem. 1991, 95, 6661−6666. (70) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollmann, P. A. Application of RESP charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation. J. Am. Chem. Soc. 1993, 115, 9620−9631.

I


Comparison of RESP and IPolQ-Mod Partial Charges for Solvation

Recommend Documents