Accuracy of free energies of hydration for organic molecules from 6-31g*-derived partial charges

Absolute free energies of hydration have been computed for 13 diverse organic molecules using partial charges derived from ab initio 6-31G* wave functions. Both Mulliken charges and charges fit to the electrostatic potential surface (EPS) were considered in conjunction with OPLS Lennard–Jones parameters for the organic molecules and the TIP4P model of water. Monte Carlo simulations with statistical perturbation theory yielded relative free energies of hydration. These were converted to absolute quantities through perturbations to reference molecules for which absolute free energies of hydration had been obtained previously in TIP4P water. The average errors in the computed absolute free energies of hydration are 1.1 kcal/mol for the 6-31G* EPS charges and 4.0 kcal/mol for the Mulliken charges. For the EPS charges, the largest individual errors are under 2 kcal/mol except for acetamide, in which case the error is 3.7 kcal/mol. The hydrogen bonding between the organic solutes and water has also been characterized. © John Wiley & Sons, Inc.     

PEARLS: program for energetic analysis of receptor-ligand system

Analysis of the energetics of small molecule ligand-protein, ligand-nucleic acid, and protein-nucleic acid interactions facilitates the quantitative understanding of molecular interactions that regulate the function and conformation of proteins. It has also been extensively used for ranking potential new ligands in virtual drug screening. We developed a Web-based software, PEARLS (Program for Energetic Analysis of Ligand-Receptor Systems), for computing interaction energies of ligand-protein, ligand-nucleic acid, protein-nucleic acid, and ligand-protein-nucleic acid complexes from their 3D structures. AMBER molecular force field, Morse potential, and empirical energy functions are used to compute the van der Waals, electrostatic, hydrogen bond, metal-ligand bonding, and water-mediated hydrogen bond energies between the binding molecules. The change in the solvation free energy of molecular binding is estimated by using an empirical solvation free energy model. Contribution from ligand conformational entropy change is also estimated by a simple model. The computed free energy for a number of PDB ligand-receptor complexes were studied and compared to experimental binding affinity. A substantial degree of correlation between the computed free energy and experimental binding affinity was found, which suggests that PEARLS may be useful in facilitating energetic analysis of ligand-protein, ligand-nucleic acid, and protein-nucleic acid interactions. PEARLS can be accessed at http://ang.cz3.nus.edu.sg/cgi-bin/prog/rune.pl.     

Atomic radii: Incorporation of solvation effects

Atomic radii used to define the solute cavity in continuum-based methods are determined by reproducing the solvent-accessible surface defined as the loci of minima in a potential (solvent interaction potential) between the solute and a probe. This potential includes electrostatic interaction (ion–dipole, ion–quadrupole, and ion-induced dipole) terms as well as a Lennard–Jones energy term. The method alleviates the need to distinguish solute atoms in different chemical environments. These radii, when used in the calculation of solvation free energies, are shown to be superior to fixed atom-specific radii or to radii obtained from the electron isodensity surface from quantum-mechanical calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1482–1493, 1998     

On the nature of DNA-duplex stability

The unwinding free energy of 128 DNA octamers was correlated with the sum of interaction energies among DNA bases and their solvation energies. The former energies were determined by using the recently developed density functional theory procedure augmented by London dispersion energy (RI-DFT-D) that provides accurate hydrogen-bonding and stacking energies highly comparable with CCSD(T)/complete basis set limit benchmark data. Efficient tight-binding DFT covering dispersion energy was also used and yielded satisfactory results. The latter method can be used for extended systems. The solvation energy was determined by using a C-PCM continuum solvent at HF level calculations. Various models were adopted to correlate theoretical energies with experimental unwinding free energies. Unless all energy components (hydrogen-bonding, intra- and interstrand-stacking, and solvation energies) were included and weighted individually, no satisfactory correlation resulted. The most advanced model yielded very close correlation (RMSE=0.32 kcal mol(-1)) fully comparable with the entirely empirical correlation introduced in the original paper. Analysis of the theoretical results shows the importance of inter- and intramolecular stacking energies, and especially the latter term plays a key role in determining DNA-duplex stabilization.     

Lennard-Jones parameters for the combined QM/MM method using the B3LYP/6-31G*/AMBER potential

A combined DFT quantum mechanical and AMBER molecular mechanical potential (QM/MM) is presented for use in molecular modeling and molecular simulations of large biological systems. In our approach we evaluate Lennard-Jones parameters describing the interaction between the quantum mechanical (QM) part of a system, which is described at the B3LYP/6-31+G* level of theory, and the molecular mechanical (MM) part of the system, described by the AMBER force field. The Lennard-Jones parameters for this potential are obtained by calculating hydrogen bond energies and hydrogen bond geometries for a large set of bimolecular systems, in which one hydrogen bond monomer is described quantum mechanically and the other is treated molecular mechanically. We have investigated more than 100 different bimolecular systems, finding very good agreement between hydrogen bond energies and geometries obtained from the combined QM/MM calculations and results obtained at the QM level of theory, especially with respect to geometry. Therefore, based on the Lennard-Jones parameters obtained in our study, we anticipate that the B3LYP/6-31+G*/AMBER potential will be a precise tool to explore intermolecular interactions inside a protein environment.     

Solvation free energies calculated using the GB/SA model: Sensitivity of results on charge sets,protocols, and force fields

The sensitivity of aqueous solvation free energies (SFEs), estimated using the GB/SA continuum solvent model, on charge sets, protocols, and force fields, was studied. Simple energy calculations using the GB/SA solvent model were performed on 11 monofunctional organic compounds. Results indicate that calculated SFEs are strongly dependent on the charge sets. Charges derived from electrostatic potential fitting to high level ab initio wave functions using the CHELPG procedure and “class IV” charges from AM1/CM1a or PM3/CM1p calculations yielded better results than the corresponding Mulliken charges. Calculated SFEs were similar to MC/FEP energies obtained in the presence of explicit TIP4P water. Further improvements were obtained by using GVB/6-31G** and MP2/6-31+G** (CHELPG) charge sets that included correlation effects. SFEs calculated using charge sets assigned by the OPLSA* force field gave the best results of all standard force fields (MM2*, MM3*, MMFF, AMBER*, and OPLSA*) implemented in MacroModel. Comparison of relative and absolute SFEs computed using either the GB/SA continuum model or MC/FEP calculations in the presence of explicit TIP4P water showed that, in general, relative SFEs can be estimated with greater accuracy. A second set of 20 mono- and difunctional molecules was also studied and relative SFEs estimated using energy minimization and thermodynamic cycle perturbation (TCP) protocols. SFEs calculated from TCP calculations using the GB/SA model were sensitive to bond lengths of dummy bonds (i.e., bonds involving dummy atoms). In such cases, keeping the bond lengths of dummy bonds close to the corresponding bond lengths of the starting structures improved the agreement of TCP-calculated SFEs with energy minimization results. Overall, these results indicate that GB/SA solvation free energy estimates from simple energy minimization calculations are of similar accuracy and value to those obtained using more elaborate TCP protocols. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 769–780, 1998     

Solvatochromic Linear Solvation Energy Relationships (LSER) for Solubility of Gases in Various Solvents by Target Factor Analysis

A linear solvation free energy relationship has been conducted to study the effects of solvent and solute properties on the free energy of solvation of inert gases and normal alkanes in different solvents. Factor analysis combined with target factor analysis was used to identify and quantify the factors controlling the variation of free energies of solvation, without the need to postulate any priori hypothetical method. Factor analysis of the solvation data revealed that there are two factors affecting the solubility of both types of gases in non‐polar as well as polar solvents. Target testing of the solvent parameters indicated that the Hildebrand solubility parameter of solvents is the major factor controlling the solubility of gases. Moreover, it was found that the coefficient of the Hildebrand solubility parameter in the linear solvation free energy equations has linear correlation with energy of vaporization and Lennard‐Jones force parameter of inert gases and number of carbon atoms and energy of vaporization of normal alkanes.     

Theoretical study of environmental effects on proton transfer reaction through the peptide bond in a model system

This study considered the possibility of proton transfer reactions through the peptide bond under different environments using the dipeptide and the 12-mer polyglycine α-helix models, in which diglycine is substituted by the 12-mer polyglycine helix. Ab initio molecular orbital calculations were carried out at the B3LYP/6-31+G(d) level of theory. To evaluate the free energies in solution, calculations of the solvation energies were performed using PCM. The correction functions on the calculated solvation energies were provided to reproduce experimental pKa values. The proton transfer reactions through the peptide bond are concluded to be possible in the protein for a wide range of proton acceptors. His complex has two free energy minima along a putative proton transfer pathway in spite of one minimum in the other complexes. The α-helix is estimated to suppress the proton transfer reactions through the peptide bond at the termini of the helix, although it is possible to proceed when the proton affinity of the acceptor is low. Contribution to the Mark S. Gordon 65th Birthday Festschrift Issue.     

Single-ion solvation free energies and the normal hydrogen electrode potential in methanol, acetonitrile, and dimethyl sulfoxide

The division of thermodynamic solvation free energies of electrolytes into contributions from individual ionic constituents is conventionally accomplished by using the single-ion solvation free energy of one reference ion, conventionally the proton, to set the single-ion scales. Thus, the determination of the free energy of solvation of the proton in various solvents is a fundamental issue of central importance in solution chemistry. In the present article, relative solvation free energies of ions and ion-solvent clusters in methanol, acetonitrile, and dimethyl sulfoxide (DMSO) have been determined using a combination of experimental and theoretical gas-phase free energies of formation, solution-phase reduction potentials and acid dissociation constants, and gas-phase clustering free energies. Applying the cluster pair approximation to differences between these relative solvation free energies leads to values of -263.5, -260.2, and -273.3 kcal/mol for the absolute solvation free energy of the proton in methanol, acetonitrile, and DMSO, respectively. The final absolute proton solvation free energies are used to assign absolute values for the normal hydrogen electrode potential and the solvation free energies of other single ions in the solvents mentioned above.     

Polarization charge densities provide a predictive quantification of hydrogen bond energies

A systematic density functional theory based study of hydrogen bond energies of 2465 single hydrogen bonds has been performed. In order to be closer to liquid phase conditions, different from the usual reference state of individual donor and acceptor molecules in vacuum, the reference state of donors and acceptors embedded in a perfect conductor as simulated by the COSMO solvation model has been used for the calculation of the hydrogen bond energies. The relationship between vacuum and conductor reference hydrogen bond energies is investigated and interpreted in the light of different physical contributions, such as electrostatic energy and dispersion. A very good correlation of the DFT/COSMO hydrogen bond energies with conductor polarization charge densities of separated donor and acceptor atoms was found. This provides a method to predict hydrogen bond strength in solution with a root mean square error of 0.36 kcal mol(-1) relative to the quantum chemical dimer calculations. The observed correlation is broadly applicable and allows for a predictive quantification of hydrogen bonding, which can be of great value in many areas of computational, medicinal and physical chemistry.     

一组和AMBER力场配合的普适波恩模型参数

提出了一套和AMBER力场相匹配的普适波恩模型参数.新的参数集包含21种原子类型的初始半径和屏蔽因子.参数通过遗传算法拟合359个小分子水合自由能的实验值得到.采用新的参数，预测了44个小分子的水合自由能，预测值和实验值能很好地吻合，而且大大优于采用Jayaram参数得到的结果.此外，采用新的参数，还预测了15个蛋白质的水合自由能，预测值和PB/SA预测得到的结果能很好地吻合.     

Calculation of solvation and binding free energy differences between VX-478 and its analogs by free energy perturbation and AMSOL methods

Summary VX-478 belongs to a novel class of HIV-1 protease inhibitors that are based on N,N-disubstituted benzene sulfonamides. Force field parameters for the N,N-dialkyl benzene sulfonamide moiety have been assembled from the literature and from our own ab initio calculations. These parameters were employed to calculate solvation and binding free energy differences between VX-478 and two analogs. The free energy perturbation method has been used to determine these differences using two approaches. In the first approach, intergroup interaction terms only were included in the calculation of free energies (as in most reports of free energy calculations using AMBER). In the second approach, both the inter- and intragroup interaction terms were included. The results obtained with the two approaches are in excellent agreement with each other and are also in close agreement with the experimental results. The solvation free energies of N,N-dimethyl benzene sulfonamide derivatives (truncated models of the inhibitors), calculated using continuum solvation (AMSOL) methods, are found to be in qualitative agreement with the experimental and free energy perturbation results. The binding and solvation free energy results are discussed in the context of structure-based drug design to show how physicochemical properties (for example aqueous solubilities and bioavailabilities) of these HIV-1 protease inhibitors were improved, while maintaining their inhibitory potency.     

Fast,efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation

We present the first global parameterization and validation of a novel charge model, called AM1-BCC, which quickly and efficiently generates high-quality atomic charges for computer simulations of organic molecules in polar media. The goal of the charge model is to produce atomic charges that emulate the HF/6-31G* electrostatic potential (ESP) of a molecule. Underlying electronic structure features, including formal charge and electron delocalization, are first captured by AM1 population charges; simple additive bond charge corrections (BCCs) are then applied to these AM1 atomic charges to produce the AM1-BCC charges. The parameterization of BCCs was carried out by fitting to the HF/6-31G* ESP of a training set of >2700 molecules. Most organic functional groups and their combinations were sampled, as well as an extensive variety of cyclic and fused bicyclic heteroaryl systems. The resulting BCC parameters allow the AM1-BCC charging scheme to handle virtually all types of organic compounds listed in The Merck Index and the NCI Database. Validation of the model was done through comparisons of hydrogen-bonded dimer energies and relative free energies of solvation using AM1-BCC charges in conjunction with the 1994 Cornell et al. forcefield for AMBER.(13) Homo- and hetero-dimer hydrogen-bond energies of a diverse set of organic molecules were reproduced to within 0.95 kcal/mol RMS deviation from the ab initio values, and for DNA dimers the energies were within 0.9 kcal/mol RMS deviation from ab initio values. The calculated relative free energies of solvation for a diverse set of monofunctional isosteres were reproduced to within 0.69 kcal/mol of experiment. In all these validation tests, AMBER with the AM1-BCC charge model maintained a correlation coefficient above 0.96. Thus, the parameters presented here for use with the AM1-BCC method present a fast, accurate, and robust alternative to HF/6-31G* ESP-fit charges for general use with the AMBER force field in computer simulations involving organic small molecules.     

Microscopic models for quantum mechanical calculations of chemical processes in solutions: LD/AMPAC and SCAAS/AMPAC calculations of solvation energies

Our previously developed approaches for integrating quantum mechanical molecular orbital methods with microscopic solvent models are refined and examined. These approaches consider the nonlinear solute–solvent coupling in a self-consistent way by incorporating the potential from the solvent dipoles in the solute Hamiltonian, while considering the polarization of the solvent by the potential from the solute charges. The solvent models used include the simplified Langevin Dipoles (LD) model and the much more expensive surface constrained All Atom Solvent (SCAAS) model, which is combined with a free energy pertubation (FEP) approach. Both methods are effectively integrated with the quantum mechanical AMPAC package and can be easily combined with other quantum mechanical programs. The advantages of the present approaches and their earlier versions over macroscopic reaction field models and supermolecular approaches are considered. A LD/MNDO study of solvated organic ions demonstrates that this model can yield reliable solvation energies, provided the quantum mechanical charges are scaled to have similar magnitudes to those obtained by high level ab initio methods. The incorporation of a field-dependent hydrophobic term in the LD free energy makes the present approach capable of evaluating the free energy of transfer of polar molecules from non polar solvents to aqueous solutions. The reliability of the LD approach is examined not only by evaluating a rather standard set of solvation energies of organic ions and polar molecules, but also by considering the stringent test case of sterically hindered hydrophobic ions. In this case, we compare the LD/MNDO solvation energies to the more rigorous FEP/SCAAS/MNDO solvation energies. Both methods are found to give similar results even in this challenging test case. The FEP/SCAAS/AMPAC method is incorporated into the current version of the program ENZYMIX. This option allows one to study chemical reactions in enzymes and in solutions using the MNDO and AM1 approximations. A special procedure that uses the EVB method as a reference potential for SCF MO calculations should help in improving the reliability of such studies.     

Anisotropic solvent model of the lipid bilayer. 1. Parameterization of long-range electrostatics and first solvation shell effects

A new implicit solvation model was developed for calculating free energies of transfer of molecules from water to any solvent with defined bulk properties. The transfer energy was calculated as a sum of the first solvation shell energy and the long-range electrostatic contribution. The first term was proportional to solvent accessible surface area and solvation parameters (σ(i)) for different atom types. The electrostatic term was computed as a product of group dipole moments and dipolar solvation parameter (η) for neutral molecules or using a modified Born equation for ions. The regression coefficients in linear dependencies of solvation parameters σ(i) and η on dielectric constant, solvatochromic polarizability parameter π*, and hydrogen-bonding donor and acceptor capacities of solvents were optimized using 1269 experimental transfer energies from 19 organic solvents to water. The root-mean-square errors for neutral compounds and ions were 0.82 and 1.61 kcal/mol, respectively. Quantification of energy components demonstrates the dominant roles of hydrophobic effect for nonpolar atoms and of hydrogen-bonding for polar atoms. The estimated first solvation shell energy outweighs the long-range electrostatics for most compounds including ions. The simplicity and computational efficiency of the model allows its application for modeling of macromolecules in anisotropic environments, such as biological membranes.     

Theoretical conformational analysis for neurotransmitters in the gas phase and in aqueous solution. Norepinephrine

The natural neurotransmitter (R)-norepinephrine takes the monocationic form in 93% abundance at the physiological tissue pH of 7.4. Ab initio and DFT/B3LYP calculations were performed for 12 protonated conformers of (R)-norepinephrine in the gas phase with geometry optimizations up to the MP2/6-311++G level, and with single-point calculations up to the QCISD(T) level at the HF/6-31G-optimized geometries. Four monohydrates were studied at the MP2/6-31G//HF/6-31G level. In the gas phase, the G1 conformer is the most stable with phenyl.NH(3)(+) gauche and HO(alc).NH(3)(+) gauche arrangements. A strained intramolecular hydrogen bond was found for conformers (G1 and T) with close NH(3)(+) and OH groups. Upon rotation of the NH(3)(+) group as a whole unit about the C(beta)-C(alpha) axis, a 3-fold potential was calculated with free energies for barriers of 3-12 kcal/mol at the HF/6-31G level. Only small deviations were found in MP2/6-311++G single-point calculations. A 2-fold potential was calculated for the phenyl rotation with free energies of 11-13 kcal/mol for the barriers at T = 310 K and p = 1 atm. A molecular mechanics docking study of (R)-norepinephrine in a model binding pocket of the beta-adrenergic receptor shows that the ligand takes a conformation close to the T(3) arrangement. The effect of aqueous solvation was considered by the free energy perturbation method implemented in Monte Carlo simulations. There are 4-5 strongly bound water molecules in hydrogen bonds to the conformers. Although hydration stabilizes mostly the G2 form with gauche phenyl.NH(3)(+) arrangement and a water-exposed NH(3)(+) group, the conformer population becomes T > G1 > G2, in agreement with the PMR spectroscopy measurements by Solmajer et al. (Z. Naturforsch. 1983, 38c, 758). Solvent effects reduce the free energies for barriers to 3-6 and 9-12 kcal/mol for rotations about the C(beta)-C(alpha) and the C(1)(ring)-C(beta) axes, respectively.     

Electrostatic solvation free energy of amino acid side chain analogs: implications for the validity of electrostatic linear response in water

Electrostatic free energies of solvation for 15 neutral amino acid side chain analogs are computed. We compare three methods of varying computational complexity and accuracy for three force fields: free energy simulations, Poisson-Boltzmann (PB), and linear response approximation (LRA) using AMBER, CHARMM, and OPLS-AA force fields. We find that deviations from simulation start at low charges for solutes. The approximate PB and LRA produce an overestimation of electrostatic solvation free energies for most of molecules studied here. These deviations are remarkably systematic. The variations among force fields are almost as large as the variations found among methods. Our study confirms that success of the approximate methods for electrostatic solvation free energies comes from their ability to evaluate free energy differences accurately.     

Theoretical study on the hydrogen bonding of five-membered heteroaromatics with water

The hydrogen-bonding ability of five-membered heteroaromatic molecules containing one chalcogen and two heteroatoms with nitrogen in addition to chalcogen, respectively, have been analyzed using density functional and molecular orbital methods through adduct formation with water. The stabilization energies for all the adducts are established at B3LYP/6-31+G* and MP2/6-31+G* levels after correcting for the basis set superposition error by using the counterpoise method and also corrected for zero-point vibrational energies. A natural bond orbital analysis at B3LYP/6-31+G* level and natural energy decomposition analysis at HF/6-31+G* using MP2/6-31+G* geometries have been carried out to understand the nature of hydrogen-bonding interaction in monohydrated heterocyclic adducts. Nucleus-independent chemical shift have been evaluated to understand the correlation between hydrogen bond formation and aromaticity.     

Accurate estimation of solvation free energy using polynomial fitting techniques

This report details an approach to improve the accuracy of free energy difference estimates using thermodynamic integration data (slope of the free energy with respect to the switching variable λ) and its application to calculating solvation free energy. The central idea is to utilize polynomial fitting schemes to approximate the thermodynamic integration data to improve the accuracy of the free energy difference estimates. Previously, we introduced the use of polynomial regression technique to fit thermodynamic integration data (Shyu and Ytreberg, J Comput Chem, 2009, 30, 2297). In this report we introduce polynomial and spline interpolation techniques. Two systems with analytically solvable relative free energies are used to test the accuracy of the interpolation approach. We also use both interpolation and regression methods to determine a small molecule solvation free energy. Our simulations show that, using such polynomial techniques and nonequidistant λ values, the solvation free energy can be estimated with high accuracy without using soft‐core scaling and separate simulations for Lennard‐Jones and partial charges. The results from our study suggest that these polynomial techniques, especially with use of nonequidistant λ values, improve the accuracy for ΔF estimates without demanding additional simulations. We also provide general guidelines for use of polynomial fitting to estimate free energy. To allow researchers to immediately utilize these methods, free software and documentation is provided via http://www.phys.uidaho.edu/ytreberg/software . © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010