We present self-consistent reaction field (SCRF) calculations, utilizing correlated ab initio quantum mechanics, of aqueous solvation free energies for a large data base of molecular solutes. We identify a subset of chemical functional groups for which there are systematic deviations in the comparison of theory and experiment; furthermore, for one case which has been extensively investigated, methylated amines, similar deviations appear in explicit solvent free energy perturbation calculations employing several commonly used molecular mechanics potential functions. By carrying out high-level ab initio quantum chemical calculations of hydrogen-bonding energies of the solutes to a water molecule, we arrive at a coherent explanation of the disagreements between theory and experiment, namely, that hydrogen-bonding energies are in some cases poorly correlated with classical electrostatic interaction energies. We show that the deviation in hydrogen-bonding energies of a solute from a reference molecule (for which there is good agreement between the SCRF calculations and experiment) is an excellent predictor of the errors made for that solute in the SCRF calculations. A new SCRF model is developed in which short-range empirical corrections, based upon solvent accessibility, are made for these chemical functional groups; this reduces the mean error of the calculated solvation free energies for the entire data base by a factor of ∼2, to 0.37 kcal/mol. These results have significant implications for the accuracy of explicit solvent potential functions as well as dielectric continuum models. Finally, we also identify cases where the observed discrepancies in solvation free energies cannot be explained by pair hydrogen-bonding results and suggest problems here that may be specific to dielectric continuum theory.

We have compared solvation free energies obtained from a number of approximate solvation models with an accurate solution of the Poisson−Boltzmann equation for a large data set of peptide structures, ranging from a single amino acid to a peptide sequence of length nine. The models are assessed for their ability to predict relative energetics of different peptide conformations (of the same sequence) as determined from the Poisson−Boltzmann results. We find that the widely used distance dependent dielectric model yields qualitatively erroneous results; in contrast, the generalized Born model of Still and co-workers, an approximation to the Poisson−Boltzmann equation, provides reasonably good solvation free energies and performs rather well in rank ordering of conformations. A surface area based model produces results of intermediate quality. Our results suggest that the generalized Born model is presently the clearly preferred alternative if one wishes to carry out molecular dynamics simulations with a fast, approximate solvation model.

The automatic three-dimensional mesh generation system for molecular geometries developed in our laboratory is used to solve the Poisson᎐Boltzmann equation numerically using a finite element method. For a number of different systems, the results are found to be in good agreement with those obtained in finite difference calculations using the DelPhi program as well as with those from boundary element calculations using our triangulated molecular surface. The overall scaling of the method is found to be approximately linear in the number of atoms in the system. The finite element mesh structure can be exploited to compute the gradient of the polarization energy in 10᎐20% of the time required to solve the equation itself. The resulting timings for the larger systems considered indicate that energies and gradients can be obtained in about half the time required for a finite difference solution to the equation. The development of a multilevel version of the algorithm as well as future applications to structure optimization using molecular mechanics force fields are also discussed.

No abstract

The derivation of a three-dimensional integral equation for solute molecule-solvent site correlation functions is presented. The equation is obtained by averaging the Ornstein–Zernicke equation for molecular liquids over orientations of the solvent molecule consistent with one site of the solvent remaining at a fixed distance from a solute-based origin. The approach is similar to that adopted in the reduction leading to the reference interaction site model (RISM) equations but retains full three-dimensional information regarding the structure of the reference solute molecule. The proposed equation can be solved using three-dimensional HNC-like closures, of which three different forms are discussed. A formulation which allows the introduction of long range interactions through a renormalization of the equation is also presented. Applications to various molecular liquids indicate that the proposed theory provides pair correlation functions that are in better agreement with molecular dynamics simulations than those obtained using the extended RISM formulation. Furthermore, qualitative errors in the correlation functions, frequently seen in results from RISM calculations are completely eliminated through geometrical averaging of the Mayer function in the 3D HNC closure. Prospects for the development of a novel mean field theory of solvation are also discussed.