Predicting aqueous free energies of solvation as functions of temperature

This work introduces a model, solvation model 6 with temperature dependence (SM6T), to predict the temperature dependence of aqueous free energies of solvation for compounds containing H, C, and O in the range 273-373 K. In particular, we extend solvation model 6 (SM6), which was previously developed (Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Chem. Theory Comput. 2005, 1, 1133) for predicting aqueous free energies of solvation at 298 K, to predict the variation of the free energy of solvation relative to 298 K. Also, we describe the database of experimental aqueous free energies of solvation for compounds containing H, C, and O that was used to parametrize and test the new model. SM6T partitions the temperature dependence of the free energy of solvation into two components: the temperature dependence of the bulk electrostatic contribution to the free energy of solvation, which is computed using the generalized Born equation, and the temperature dependence of first-solvation-shell effects which is modeled using a parametrized solvent-exposed surface-area-dependent term. We found that SM6T predicts the temperature dependence of aqueous free energies of solvation with a mean unsigned error of 0.08 kcal/mol over our entire database, whereas using the experimental value at 298 K produces a mean unsigned error of 0.53 kcal/mol.     

The ionic work function and its role in estimating absolute electrode potentials

The experimental determination of the ionic work function is briefly described. Data for the proton, alkali metal ions, and halide ions in water, originally published by Randles (Randles, J. E. B. Trans Faraday Soc. 1956, 52, 1573) are recalculated on the basis of up-to-date thermodynamic tables. These calculations are extended to data for the same ions in four nonaqueous solvents, namely, methanol, ethanol, acetonitrile, and dimethyl sulfoxide. The ionic work function data are compared with estimates of the absolute Gibbs energy of solvation obtained by an extrathermodynamic route for the same ions. The work function data for the proton are used to estimate the absolute potential of the standard hydrogen electrode in each solvent. The results obtained here are compared with those published earlier by Trasatti (Trasatti, S. Electrochim. Acta 1987, 32, 843) and more recently by Kelly et al. (Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2006, 110, 16066. Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2007, 111, 408). A comparison of the ionic work function with the absolute Gibbs solvation energy permits an estimation of the surface potential of the solvent. The results show that the surface potential of water is small and positive whereas the surface potential of the nonaqueous solvents considered is negative. The sign of the surface potential is consistent with the known structure of each solvent.     

Improved methods for semiempirical solvation models

We present improved algorithms for the SMx (x = 1, 1a, 2, 3) solvation models presented previously [see the overview in C. J. Cramer and D. G. Truhlar, J. Comp.-Aided Mol. Design, 6 , 629 (1992)]. These models estimate the free energy of solvation by augmenting a semiempirical Hartree-Fock calculation on the solute with the generalized Born (GB) model for electric polarization of the solvent and a surface tension term based on solvent-accessible surface area. This article presents three improvements in the algorithms used to carry out such calculations, namely (1) an analytical accessible surface area algorithm, (2) a more efficient radial integration scheme for the dielectric screening computation in the GB model, and (3) a damping algorithm for updating the GB contribution to the Fock update during the iterations to achieve a self-consistent field. Improvements (1) and (2) decrease the computer time, and improvement (3) leads to more stable convergence. Improvement (2) removes a small systematic numerical error that was explicitly absorbed into the parameterization in the SMx models. Therefore, we have adjusted the parameters for one of the previous models to yield essentially identical performance as was obtained originally while simultaneously taking advantage of improvement (2). The resulting model is called SM2.1. The fact that we obtain similar results after removing the systematic quadrature bias attests to the robustness of the original parameterization. © 1995 by John Wiley & Sons, Inc.     

The effect of aqueous solvation upon alpha-helix formation for polyalanines

The incremental free energies of aqueous solution for acetyl(ala)NNH2 in its extended unfolded and alpha-helical conformations are compared using the SM5.2 solvation method of Cramer and Truhlar. A combination of density functional theory (DFT) at the B3LYP/D95(d,p) and AM1 has been employed using the ONIOM method. The incremental solvation energies of alpha-helical structures are very similar for both ONIOM and AM1 optimized structures as these structures do not significantly change upon solution. However, the conformations of the unfolded peptides change from extended beta-strand to polyproline II conformations upon aqueous solution. The incremental solvation free energy per residue of the polyproline II structure is about 2 kcal/mol/residue greater than that for the alpha-helix, representing an upper limit for the difference between the solvation energies. However, most of this difference disappears when the energy required to distort the optimized gas-phase extended beta-strand structure to the optimized polyproline II solution structure is included in the analysis, leaving an estimated difference in incremental solvation free energy of 0.3-0.5 kcal/mol favoring the unfolded structure. The solution structure sacrifices the stability derived from the intramolecular C5 H-bonds for more favorable interactions with the aqueous solvent.     

New analytic approximation to the standard molecular volume definition and its application to generalized Born calculations

In a recent article (Lee, M. S.; Salsbury, F. R. Jr.; Brooks, C. L., III. J Chem Phys 2002, 116, 10606), we demonstrated that generalized Born (GB) theory provides a good approximation to Poisson electrostatic solvation energy calculations if one uses the same definitions of molecular volume for each. In this work, we present a new and improved analytic method for reproducing the Lee-Richards molecular volume, which is the most common volume definition for Poisson calculations. Overall, 1% errors are achieved for absolute solvation energies of a large set of proteins and relative solvation energies of protein conformations. We also introduce an accurate SASA approximation that uses the same machinery employed by our GB method and requires a small addition of computational cost. The combined methodology is shown to yield an efficient and accurate implicit solvent representation for simulations of biopolymers.     

Aqueous solvation free energies of ions and ion-water clusters based on an accurate value for the absolute aqueous solvation free energy of the proton

Thermochemical cycles that involve pKa, gas-phase acidities, aqueous solvation free energies of neutral species, and gas-phase clustering free energies have been used with the cluster pair approximation to determine the absolute aqueous solvation free energy of the proton. The best value obtained in this work is in good agreement with the value reported by Tissandier et al. (Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. J.; Earhart, A. D.; Coe, J. V. J. Phys. Chem. A 1998, 102, 7787), who applied the cluster pair approximation to a less diverse and smaller data set of ions. We agree with previous workers who advocated the value of -265.9 kcal/mol for the absolute aqueous solvation free energy of the proton. Considering the uncertainties associated with the experimental gas-phase free energies of ions that are required to use the cluster pair approximation as well as analyses of various subsets of data, we estimate an uncertainty for the absolute aqueous solvation free energy of the proton of no less than 2 kcal/mol. Using a value of -265.9 kcal/mol for the absolute aqueous solvation free energy of the proton, we expand and update our previous compilation of absolute aqueous solvation free energies; this new data set contains conventional and absolute aqueous solvation free energies for 121 unclustered ions (not including the proton) and 147 conventional and absolute aqueous solvation free energies for 51 clustered ions containing from 1 to 6 water molecules. When tested against the same set of ions that was recently used to develop the SM6 continuum solvation model, SM6 retains its previously determined high accuracy; indeed, in most cases the mean unsigned error improves when it is tested against the more accurate reference data.     

A universal approach for continuum solvent pK a calculations: are we there yet?

This paper reviews several pK _a calculation strategies that are commonly used in aqueous acidity predictions. Among those investigated were the direct or absolute method, the proton exchange scheme, and the hybrid cluster–continuum (Pliego and Riveros) and implicit–explicit (Kelly, Cramer and Truhlar) models. Additionally, these protocols are applied in the pK _a calculation of 55 neutral organic and inorganic acids in conjunction with various solvent models, including the CPCM-UAKS/UAHF, IPCM, SM6 and COSMO-RS, with a view to identifying a universal approach for accurate pK _a predictions. The results indicate that the direct method is unsuitable for general pK _a calculations, although moderately accurate results (MAD <3 units) are possible for certain classes of acids, depending on the choice of solvent model. The proton exchange scheme generally delivers good results (MAD <2 units), with CPCM-UAKS giving the best performance. Furthermore, the sensitivity of this approach to the choice of reference acid can be substantially lessened if the solvation energies for ionic species are calculated via the IPCM cluster–continuum approach. Reference-independent hybrid approaches that include explicit water molecules can potentially give reasonably accurate values (MAD generally ~2 units) depending on the solvent model and the number of explicit water molecules added.     

Calculating solvation energies by means of a fluctuating charge model combined with continuum solvent model

Continuum solvent models have shown to be very efficient for calculating solvation energy of biomolecules in solution. However, in order to produce accurate results, besides atomic radii or volumes, an appropriate set of partial charges of the molecule is needed. Here, a set of partial charges produced by a fluctuating charge model-the atom-bond electronegativity equalization method model (ABEEMσπ) fused into molecular mechanics is used to fit for the analytical continuum electrostatics model of generalized-Born calculations. Because the partial atomic charges provided by the ABEEMσπ model can well reflect the polarization effect of the solute induced by the continuum solvent in solution, accurate and rapid calculations of the solvation energies have been performed for series of compounds involving 105 small neutral molecules, twenty kinds of dipeptides and several protein fragments. The solvation energies of small neutral molecules computed with the combination of the GB model with the fluctuating charge protocol (ABEEMσπ∕GB) show remarkable agreement with the experimental results, with a correlation coefficient of 0.97, a slope of 0.95, and a bias of 0.34 kcal∕mol. Furthermore, for twenty kinds of dipeptides and several protein fragments, the results obtained from the analytical ABEEMσπ∕GB model calculations correlate well with those from ab initio and Poisson-Boltzmann calculations. The remarkable agreement between the solvation energies computed with the ABEEMσπ∕GB model and PB model provides strong motivation for the use of ABEEMσπ∕GB solvent model in the simulation of biochemical systems.     

Some comments on valence bond representations for the radical exchange reaction X*+R:Y-->X:R+Y*

Qualitative valence bond formulations by Hiberty and co-workers (Hiberty, P. C.; Megret, C.; Song, L.; Wu, W.; Shaik, S. J. Am. Chem. Soc. 2006, 128, 2836) of mechanisms for the radical exchange reactions H*+F:H-->H:F+H* and F*+H:F-->F:H+F* are compared to a previously published formulation of the generalized radical exchange reaction X*+R:Y-->X:R+Y*. The former formulation uses covalent-ionic VB complexes, and the latter formulation, which is more general, involves the formation of reactant-like and product-like complexes at intermediate stages along the reaction coordinate.     

Optimizing the performance of the multiconfiguration molecular mechanics method

Multiconfiguration molecular mechanics (MCMM) is a general algorithm for constructing potential energy surfaces for reactive systems (Kim, Y.; Corchado, J. C.; Villà, J.; Xing, J.; Truhlar, D. G. J. Chem. Phys. 2000, 112, 2718). This paper illustrates how the performance of the MCMM method can be improved by adopting improved molecular mechanics parameters. We carry out calculations of reaction rate constants using variational transition state theory with optimized multidimensional tunneling on the MCMM PESs for three hydrogen transfer reactions, and we compare the results to direct dynamics. We find that the MCMM method with as little as one electronic structure Hessian can describe the dynamically important regions of the ground-electronic state PES, including the corner-cutting-tunneling region of the reaction swath, with practical accuracy.     

Multiple intermolecular interaction modes of positively charged residues with adenine in ATP-binding proteins

Adenosine 5'-triphosphate (ATP) plays an essential role in all forms of life. Molecular recognition of ATP in ATP-binding proteins is a subject of great importance for understanding enzymatic mechanisms and for drug design. We have carried out a large-scale data mining of the Protein Data Bank (PDB) to analyze molecular determinants for recognition of ATP, in particular, the adenine base, by ATP-binding proteins. A novel distribution pattern of charged residues around the adenine base was discovered: lysine residues tend to occupy the major groove N7 side of the adenine base, and the arginine residues situate preferentially above or below the adenine bases. Such an arrangement is advantageous because it facilitates multiple modes of intermolecular interactions, that is, cation-pi interactions and a hydrogen bond between lysine and adenine, and cation-pi and pi-pi stacking interactions between arginine and adenine. For the two representative Lys...Adenine and Arg...Adenine interactions, intermolecular interaction energies were subsequently analyzed by means of the supermolecular approach at the MP2 level with solvation free energy correction using the SM5.42R model of Cramer and Truhlar, which gave rise to significant interaction strengths.     

Generalized born model with a simple smoothing function

Based on recent developments in generalized Born (GB) theory that employ rapid volume integration schemes (M. S. Lee, F. R. Salabury, Jr., and C. L. Brooks III, J Chem Phys 2002, 116, 10606) we have recast the calculation of the self-electrostatic solvation energy to utilize a simple smoothing function at the dielectric boundary. The present GB model is formulated in this manner to provide consistency with the Poisson-Boltzmann (PB) theory previously developed to yield numerically stable electrostatic solvation forces based on finite-difference methods (W. Im, D. Beglov, and B. Roux, Comp Phys Commun 1998, 111, 59). Our comparisons show that the present GB model is indeed an efficient and accurate approach to reproduce corresponding PB solvation energies and forces. With only two adjustable parameters--a(0) to modulate the Coulomb field term, and a(1) to include a correction term beyond Coulomb field--the PB solvation energies are reproduced within 1% error on average for a variety of proteins. Detailed analysis shows that the PB energy can be reproduced within 2% absolute error with a confidence of about 95%. In addition, the solvent-exposed surface area of a biomolecule, as commonly used in calculations of the nonpolar solvation energy, can be calculated accurately and efficiently using the simple smoothing function and the volume integration method. Our implicit solvent GB calculations are about 4.5 times slower than the corresponding vacuum calculations. Using the simple smoothing function makes the present GB model roughly three times faster than GB models, which attempt to mimic the Lee-Richards molecular volume.     

Heterolytic bond dissociation in water: why is it so easy for C4H9Cl but not for C3H9SiCl?

The recently developed (Song, L.; Wu, W.; Zhang, Q.; Shaik, S. J. Phys. Chem. A 2004, 108, 6017-6024) valence bond method coupled to a polarized continuum model (VBPCM) is used to address the long standing conundrum of the heterolytic dissociation of the C-Cl and Si-Cl bonds, respectively, in tertiary-butyl chloride and trimethylsilyl chloride in condensed phases. The method is used here to compare the bond dissociation in the gas phase and in aqueous solution. In addition to the ground state reaction profile, VB theory also provides the energies of the purely covalent and purely ionic VB structures as a function of the reaction coordinate. Accordingly, the C-Cl and Si-Cl bonds are shown to be of different natures. In the gas phase, the resonance energy arising from covalent-ionic mixing at equilibrium geometry amounts to 42 kcal/mol for tertiary-butyl chloride, whereas the same quantity for trimethylsilyl chloride is significantly higher at 62 kcal/mol. With such a high value, the root cause of the Si-Cl bonding is the covalent-ionic resonance energy, and this bond belongs to the category of charge-shift bonds (Shaik, S.; Danovich, D.; Silvi, B.; Lauvergnat, D.; Hiberty, P. C. Chem.- Eur. J. 2005, 11, 6358). This difference between the C-Cl and Si-Cl bonds carries over to the solvated phase and impacts the heterolytic cleavages of the two bonds. For both molecules, solvation lowers the ionic curve below the covalent one, and hence the bond dissociation in the solvent generates the two ions, Me3E+ Cl- (E = C, Si). In both cases, the root cause of the barrier is the loss of the covalent-ionic resonance energy. In the heterolysis reaction of Si-Cl, the covalent-ionic resonance energy remains large and fully contributes to the dissociation energy, thereby leading to a high barrier for heterolytic cleavage, and thus prohibiting the generation of ions. By contrast, the covalent-ionic resonance energy is smaller for the C-Cl bond and only partially contributes to the barrier for heterolysis, which is consequently small, leading readily to ions that are commonly observed in the classical SN1 mechanism. Thus, the reluctance of R3Si-X molecules to undergo heterolysis in condensed phases and more generally the rarity of free silicenium ions under these conditions are experimental manifestations of the charge-shift character of the Si-Cl bond.     

Solvent effects on electrophilicity

Continuum solvent effect on the electrophilicity index recently proposed by Parr and co-workers (Parr, R. G.; von Szentpaly, L.; Liu, S. J. Am. Chem. Soc. 1999, 121, 1922) is discussed in detail. Solvent effect is introduced using the self-consistent isodensity polarized continuum model (SCI-PCM). A linear relationship is found between the change in electrophilicity index and the solvation energy as represented in the frame of the reaction field theory. The effect of a polarizable environment on the global electrophilicity is examined for a series of 18 well-known electrophiles presenting a wide diversity in structure and bonding properties. It is found that solvation enhances the electrophilicity power of neutral electrophilic ligands but attenuates this power in charged and ionic electrophiles.     

Accuracy of free energies of hydration using CM1 and CM3 atomic charges

Absolute free energies of hydration (DeltaGhyd) have been computed for 25 diverse organic molecules using partial atomic charges derived from AM1 and PM3 wave functions via the CM1 and CM3 procedures of Cramer, Truhlar, and coworkers. Comparisons are made with results using charges fit to the electrostatic potential surface (EPS) from ab initio 6-31G* wave functions and from the OPLS-AA force field. OPLS Lennard-Jones parameters for the organic molecules were used together with the TIP4P water model in Monte Carlo simulations with free energy perturbation theory. Absolute free energies of hydration were computed for OPLS united-atom and all-atom methane by annihilating the solutes in water and in the gas phase, and absolute DeltaGhyd values for all other molecules were computed via transformation to one of these references. Optimal charge scaling factors were determined by minimizing the unsigned average error between experimental and calculated hydration free energies. The PM3-based charge models do not lead to lower average errors than obtained with the EPS charges for the subset of 13 molecules in the original study. However, improvement is obtained by scaling the CM1A partial charges by 1.14 and the CM3A charges by 1.15, which leads to average errors of 1.0 and 1.1 kcal/mol for the full set of 25 molecules. The scaled CM1A charges also yield the best results for the hydration of amides including the E/Z free-energy difference for N-methylacetamide in water.     

Implementation of molecular symmetry in valence bond calculation

A novel strategy for the construction of many-electron symmetry-adapted wave function is proposed for ab initio valence bond (VB) calculations and is implemented for valence bond self-consistent filed (VBSCF) and breathing orbital valence bond (BOVB) methods with various orbital optimization algorithms. Symmetry-adapted VB functions are constructed by the projection operator of symmetry group. The many-electron symmetry-adapted wave function is expressed in terms of symmetry-adapted VB functions, and thus the VB calculations can be performed with the molecular symmetry restriction. Test results show that molecular symmetry reduces the computational cost of both the iteration numbers and CPU time. Furthermore, excited states with specific symmetry can be conveniently obtained in VB calculations by using symmetry-adapted VB functions.     

VBEFP: a valence bond approach that incorporates effective fragment potential method

An ab initio explicit solvation valence bond (VB) method, called VBEFP, is presented. The VBEFP method is one type of QM/MM approach in which the QM part of system is treated within the ab initio valence bond scheme and the solvent water molecules are accounted by the effective fragment potential (EFP) method, which is a polarized force field approach developed by Gordon et al. (J. Chem. Phys. 1996, 105, 1968). This hybrid method enables one to take the first-solvation shell and heterogeneous solvation effects into account explicitly with VB wave function. Therefore, the nature of chemical bonding and the mechanism of chemical reactions with explicit solvent environments can be explored at the ab inito VB level. In this paper, the hydrated metal-ligand complexes [M(2+)L](H(2)O)(n) (M(2+): Mg(2+), Zn(2+); L: NH(3), CH(2)O) are studied by the VBEFP method. Resonance energy and bond order are computed, and the influence of the solvent coordination and hydrogen bonding to the metal-ligand bonding are explored in the paper.     

Electrostatic component of solvation: comparison of SCRF continuum models

We report a systematic comparison of the electrostatic contributions to the free energy of solvation from three different kinds of quantum mechanical self-consistent reaction field (SCRF) methods. We also compare the liquid-phase dipole moments as a measure of the solute's response to the reaction field of the solvent. In particular, we compare these quantities for the generalized Born model as implemented in the SM5.42R method, the multipolar expansion model developed at Nancy, and the MST version of the polarizable continuum model. All calculations are carried out at the HF/6-31G(d) level. The effects of various choices of solute cavities and representations of the charge density are examined. The test set consists of 18 molecules containing prototypical polar groups, and three different values of the dielectric permittivity are considered.     

Teaching Chemistry with Electron Density Models. 2. Can Atomic Charges Adequately Explain Electrostatic Potential Maps?

Electrostatic potential maps generated from quantum mechanical calculations are widely used to teach students about molecular polarity and assign atomic charges (Shusterman, G. P.; Shusterman, A. J. J. Chem. Educ. 1997, 74, 771–776; Hehre, W. J.; Shusterman, A. J.; Nelson, J. E. The Molecular Modeling Workbook for Organic Chemistry: Wavefunction: Irvine, CA, 1998). The assumption that potential equals charge is only valid, however, when comparing atoms of similar size. The proper use of potential maps requires consideration of atomic charge, atomic radius, and the electron configuration (orbital occupancy) of the atom in question. These points are illustrated through the analysis of the potential maps of various halogen-containing molecules.The first article in this series is Shusterman, G. P.; Shusterman, A. J. J. Chem. Educ. 1997, 74, 771.776.