Study of peptide conformation in terms of the ABEEM/MM method

The ABEEM/MM model (atom-bond electronegativity equalization method fused into molecular mechanics) is applied to study of the polypeptide conformations. The Lennard-Jones and torsional parameters were optimized to be consistent with the ABEEM/MM fluctuating charge electrostatic potential. The hydrogen bond was specially treated with an electrostatic fitting function. Molecular dipole moments, dimerization energies, and hydrogen bond lengths of complexes are reasonably achieved by our model, compared to ab initio results. The ABEEM/MM fluctuating charge model reproduces both the peptide conformational energies and structures with satisfactory accuracy with low computer cost. The transferability is tested by applying the parameters of our model to the tetrapeptide of alanine and another four dipeptides. The overall RMS deviations in conformational energies and key dihedral angles for four di- or tetrapeptide, is 0.39 kcal/mol and 7.7 degrees . The current results agree well with those by the accurate ab initio method, and are comparable to those from the best existing force fields. The results make us believe that our fluctuating charge model can obtain more promising results in protein and macromolecular modeling with good accuracy but less computer cost.     

Alcohols,ethers, carbohydrates,and related compounds. IV. Carbohydrates

Ab initio calculations [B3LYP/6-311++G(2d,2p)] have been carried out on 84 conformations of 12 different sugars (hexoses), in both pyranose and furanose forms, with the idea of generating a data base for carbohydrate structural energies that may be used for developing the predictive value of molecular mechanics calculations for carbohydrates. The average value for the apparent gas phase anomeric effect for a series of 31 pairs of pyranose conformations was found to be 1.83 kcal/mol (vs. 2.67 kcal/mol with a smaller basis set used in earlier calculations). In developing MM4 to reproduce these data, it was necessary first to have good energies for simple alcohols and ethers, together with an adequate treatment of hydrogen bonding, and then to include the anomeric effect, and the ethylene glycol type system, as was previously recognized. It was also found that the so-called delta-2 effect, long recognized in carbohydrates, must be explicitly included, in order to obtain acceptable results. When a force field that included all of these items as developed from the small molecules based on the MM4 hydrocarbon force field was applied without any parameter adjustment to the set of hexopyranose and furanose conformations mentioned earlier, the E(beta) - E(alpha) was found to have an average value of 1.88 kcal/mol, versus 1.74 for the quantum calculations. The signed average and RMS deviations of the MM4 from the QM results were +0.15 and 0.87 kcal/mol.     

Reactivity of isobutane on zeolites: a first principles study

In this work, ab initio and density functional theory methods are used to study isobutane protolytic cracking, primary hydrogen exchange, tertiary hydrogen exchange, and dehydrogenation reactions catalyzed by zeolites. The reactants, products, and transition-state structures are optimized at the B3LYP/6-31G* level, and the final energies are calculated using the CBS-QB3 composite energy method. The computed activation barriers are 52.3 kcal/mol for cracking, 29.4 kcal/mol for primary hydrogen exchange, 29.9 kcal/mol for tertiary hydrogen exchange, and 59.4 kcal/mol for dehydrogenation. The zeolite acidity effects on the reaction barriers are also investigated by changing the cluster terminal Si-H bond lengths. The analytical expressions between activation barriers and zeolite deprotonation energies for each reaction are proposed so that accurate activation barriers can be obtained when using different zeolites as catalysts.     

Accurate interaction energies of hydrogen-bonded nucleic acid base pairs

Hydrogen-bonded nucleic acids base pairs substantially contribute to the structure and stability of nucleic acids. The study presents reference ab initio structures and interaction energies of selected base pairs with binding energies ranging from -5 to -47 kcal/mol. The molecular structures are obtained using the RI-MP2 (resolution of identity MP2) method with extended cc-pVTZ basis set of atomic orbitals. The RI-MP2 method provides results essentially identical with the standard MP2 method. The interaction energies are calculated using the Complete Basis Set (CBS) extrapolation at the RI-MP2 level. For some base pairs, Coupled-Cluster corrections with inclusion of noniterative triple contributions (CCSD(T)) are given. The calculations are compared with selected medium quality methods. The PW91 DFT functional with the 6-31G basis set matches well the RI-MP2/CBS absolute interaction energies and reproduces the relative values of base pairing energies with a maximum relative error of 2.6 kcal/mol when applied with Becke3LYP-optimized geometries. The Becke3LYP DFT functional underestimates the interaction energies by few kcal/mol with relative error of 2.2 kcal/mol. Very good performance of nonpolarizable Cornell et al. force field is confirmed and this indirectly supports the view that H-bonded base pairs are primarily stabilized by electrostatic interactions.     

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.     

Development and testing of a general amber force field

We describe here a general Amber force field (GAFF) for organic molecules. GAFF is designed to be compatible with existing Amber force fields for proteins and nucleic acids, and has parameters for most organic and pharmaceutical molecules that are composed of H, C, N, O, S, P, and halogens. It uses a simple functional form and a limited number of atom types, but incorporates both empirical and heuristic models to estimate force constants and partial atomic charges. The performance of GAFF in test cases is encouraging. In test I, 74 crystallographic structures were compared to GAFF minimized structures, with a root-mean-square displacement of 0.26 A, which is comparable to that of the Tripos 5.2 force field (0.25 A) and better than those of MMFF 94 and CHARMm (0.47 and 0.44 A, respectively). In test II, gas phase minimizations were performed on 22 nucleic acid base pairs, and the minimized structures and intermolecular energies were compared to MP2/6-31G* results. The RMS of displacements and relative energies were 0.25 A and 1.2 kcal/mol, respectively. These data are comparable to results from Parm99/RESP (0.16 A and 1.18 kcal/mol, respectively), which were parameterized to these base pairs. Test III looked at the relative energies of 71 conformational pairs that were used in development of the Parm99 force field. The RMS error in relative energies (compared to experiment) is about 0.5 kcal/mol. GAFF can be applied to wide range of molecules in an automatic fashion, making it suitable for rational drug design and database searching.     

Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations

Combined ab initio quantum mechanical and molecular mechanical calculations have been widely used for modeling chemical reactions in complex systems such as enzymes, with most applications being based on the determination of a minimum energy path connecting the reactant through the transition state to the product in the enzyme environment. However, statistical mechanics sampling and reaction dynamics calculations with a combined ab initio quantum mechanical (QM) and molecular mechanical (MM) potential are still not feasible because of the computational costs associated mainly with the ab initio quantum mechanical calculations for the QM subsystem. To address this issue, a reaction path potential energy surface is developed here for statistical mechanics and dynamics simulation of chemical reactions in enzymes and other complex systems. The reaction path potential follows the ideas from the reaction path Hamiltonian of Miller, Handy and Adams for gas phase chemical reactions but is designed specifically for large systems that are described with combined ab initio quantum mechanical and molecular mechanical methods. The reaction path potential is an analytical energy expression of the combined quantum mechanical and molecular mechanical potential energy along the minimum energy path. An expansion around the minimum energy path is made in both the nuclear and the electronic degrees of freedom for the QM subsystem internal energy, while the energy of the subsystem described with MM remains unchanged from that in the combined quantum mechanical and molecular mechanical expression and the electrostatic interaction between the QM and MM subsystems is described as the interaction of the MM charges with the QM charges. The QM charges are polarizable in response to the changes in both the MM and the QM degrees of freedom through a new response kernel developed in the present work. The input data for constructing the reaction path potential are energies, vibrational frequencies, and electron density response properties of the QM subsystem along the minimum energy path, all of which can be obtained from the combined quantum mechanical and molecular mechanical calculations. Once constructed, it costs much less for its evaluation. Thus, the reaction path potential provides a potential energy surface for rigorous statistical mechanics and reaction dynamics calculations of complex systems. As an example, the method is applied to the statistical mechanical calculations for the potential of mean force of the chemical reaction in triosephosphate isomerase.     

Convergence in the QM‐only and QM/MM modeling of enzymatic reactions: A case study for acetylene hydratase

We report systematic quantum mechanics‐only (QM‐only) and QM/molecular mechanics (MM) calculations on an enzyme‐catalyzed reaction to assess the convergence behavior of QM‐only and QM/MM energies with respect to the size of the chosen QM region. The QM and MM parts are described by density functional theory (typically B3LYP/def2‐SVP) and the CHARMM force field, respectively. Extending our previous work on acetylene hydratase with QM regions up to 157 atoms (Liao and Thiel, J. Chem. Theory Comput. 2012, 8, 3793), we performed QM/MM geometry optimizations with a QM region M4 composed of 408 atoms, as well as further QM/MM single‐point calculations with even larger QM regions up to 657 atoms. A charge deletion analysis was conducted for the previously used QM/MM model ( M3a , with a QM region of 157 atoms) to identify all MM residues with strong electrostatic contributions to the reaction energetics (typically more than 2 kcal/mol), which were then included in M4 . QM/MM calculations with this large QM region M4 lead to the same overall mechanism as the previous QM/MM calculations with M3a , but there are some variations in the relative energies of the stationary points, with a mean absolute deviation (MAD) of 2.7 kcal/mol. The energies of the two relevant transition states are close to each other at all levels applied (typically within 2 kcal/mol), with the first (second) one being rate‐limiting in the QM/MM calculations with M3a ( M4 ). QM‐only gas‐phase calculations give a very similar energy profile for QM region M4 (MAD of 1.7 kcal/mol), contrary to the situation for M3a where we had previously found significant discrepancies between the QM‐only and QM/MM results (MAD of 7.9 kcal/mol). Extension of the QM region beyond M4 up to M7 (657 atoms) leads to only rather small variations in the relative energies from single‐point QM‐only and QM/MM calculations (MAD typically about 1–2 kcal/mol). In the case of acetylene hydratase, a model with 408 QM atoms thus seems sufficient to achieve convergence in the computed relative energies to within 1–2 kcal/mol.Copyright © 2013 Wiley Periodicals, Inc.     

Comparison of linear-scaling semiempirical methods and combined quantum mechanical/molecular mechanical methods for enzymic reactions. II. An energy decomposition analysis

QM/MM methods have been developed as a computationally feasible solution to QM simulation of chemical processes, such as enzyme-catalyzed reactions, within a more approximate MM representation of the condensed-phase environment. However, there has been no independent method for checking the quality of this representation, especially for highly nonisotropic protein environments such as those surrounding enzyme active sites. Hence, the validity of QM/MM methods is largely untested. Here we use the possibility of performing all-QM calculations at the semiempirical PM3 level with a linear-scaling method (MOZYME) to assess the performance of a QM/MM method (PM3/AMBER94 force field). Using two model pathways for the hydride-ion transfer reaction of the enzyme dihydrofolate reductase studied previously (Titmuss et al., Chem Phys Lett 2000, 320, 169-176), we have analyzed the reaction energy contributions (QM, QM/MM, and MM) from the QM/MM results and compared them with analogous-region components calculated via an energy partitioning scheme implemented into MOZYME. This analysis further divided the MOZYME components into Coulomb, resonance and exchange energy terms. For the model in which the MM coordinates are kept fixed during the reaction, we find that the MOZYME and QM/MM total energy profiles agree very well, but that there are significant differences in the energy components. Most significantly there is a large change (approximately 16 kcal/mol) in the MOZYME MM component due to polarization of the MM region surrounding the active site, and which arises mostly from MM atoms close to (<10 A) the active-site QM region, which is not modelled explicitly by our QM/MM method. However, for the model where the MM coordinates are allowed to vary during the reaction, we find large differences in the MOZYME and QM/MM total energy profiles, with a discrepancy of 52 kcal/mol between the relative reaction (product-reactant) energies. This is largely due to a difference in the MM energies of 58 kcal/mol, of which we can attribute approximately 40 kcal/mol to geometry effects in the MM region and the remainder, as before, to MM region polarization. Contrary to the fixed-geometry model, there is no correlation of the MM energy changes with distance from the QM region, nor are they contributed by only a few residues. Overall, the results suggest that merely extending the size of the QM region in the QM/MM calculation is not a universal solution to the MOZYME- and QM/MM-method differences. They also suggest that attaching physical significance to MOZYME Coulomb, resonance and exchange components is problematic. Although we conclude that it would be possible to reparameterize the QM/MM force field to reproduce MOZYME energies, a better way to account for both the effects of the protein environment and known deficiencies in semiempirical methods would be to parameterize the force field based on data from DFT or ab initio QM linear-scaling calculations. Such a force field could be used efficiently in MD simulations to calculate free energies.     

Ab initio quantum mechanics-based free energy perturbation method for calculating relative solvation free energies

A free energy perturbation (FEP) method was developed that uses ab initio quantum mechanics (QM) for treating the solute molecules and molecular mechanics (MM) for treating the surroundings. Like our earlier results using AM1 semi empirical QMs, the ab initio QM/MM-based FEP method was shown to accurately calculate relative solvation free energies for a diverse set of small molecules that differ significantly in structure, aromaticity, hydrogen bonding potential, and electron density. Accuracy was similar to or better than conventional FEP methods. The QM/MM-based methods eliminate the need for time-consuming development of MM force field parameters, which are frequently required for drug-like molecules containing structural motifs not adequately described by MM. Future automation of the method and parallelization of the code for Linux 128/256/512 clusters is expected to enhance the speed and increase its use for drug design and lead optimization.     

Mapping the interaction energy surfaces of cyclic alkanes: evaluating the transferability of an ab initio based potential model

Detailed interaction energy maps are computed for symmetric cyclopropane and tetrahedrane dimer systems using ab initio methods. Interaction energies of cubane and cyclohexane dimers are also reported. The global minimum energy structures of cyclopropane and tetrahedrane systems are both D(3d) structures with energies of -1.850 and -2.171 kcal mol(-1). The ability of NIPE potential model, based on ab initio nonbonding data of neopentane (N), isobutane (I), propane (P), ethane (E) and all their combinations to predict the pair interaction energies of these strained cyclic hydrocarbons is also investigated. The difference between the energies predicted by NIPE and those obtained from the ab initio calculations increases with ring strain In general, NIPE values are in close agreement with the ab initio results for alkane ring structures having low ring strain.     

Theoretical studies on the effects of methods and parameterization on the calculated free energy of hydration for small molecules

Free energies of hydration (FEH) have been computed for 13 neutral and nine ionic species as a difference of theoretically calculated Gibbs free energies in solution and in the gas phase. In‐solution calculations have been performed using both SCIPCM and PCM polarizable continuum models at the density functional theory (DFT)/B3LYP and ab initio Hartree–Fock levels with two basis sets (6‐31G* and 6‐311++G**). Good linear correlation has been obtained for calculated and experimental gas‐phase dipole moments, with an increase by ～30% upon solvation due to solute polarization. The geometry distortion in solution turns out to be small, whereas solute polarization energies are up to 3 kcal/mol for neutral molecules. Calculation of free energies of hydration with PCM provides a balanced set of values with 6‐31G* and 6‐311++G** basis sets for neutral molecules and ionic species, respectively. Explicit solvent calculations within Monte Carlo simulations applying free energy perturbation methods have been considered for 12 neutral molecules. Four different partial atomic charge sets have been studied, obtained by a fit to the gas‐phase and in‐solution molecular electrostatic potentials at in‐solution optimized geometries. Calculated FEH values depend on the charge set and the atom model used. Results indicate a preference for the all‐atom model and partial charges obtained by a fit to the molecular electrostatic potential of the solute computed at the SCIPCM/B3LYP/6‐31G* level. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2004     

Molecular mechanics models for tetracycline analogs

Tetracyclines (Tcs) are an important family of antibiotics that bind to the ribosome and several proteins. To model Tc interactions with protein and RNA, we have developed a molecular mechanics force field for 12 tetracyclines, consistent with the CHARMM force field. We considered each Tc variant in its zwitterionic tautomer, with and without a bound Mg(2+). We used structures from the Cambridge Crystallographic Data Base to identify the conformations likely to be present in solution and in biomolecular complexes. A conformational search by simulated annealing was undertaken, using the MM3 force field, for tetracycline, anhydrotetracycline, doxycycline, and tigecycline. Resulting, low-energy structures were optimized with an ab initio method. We found that Tc and its analogs all adopt an extended conformation in the zwitterionic tautomer and a twisted one in the neutral tautomer, and the zwitterionic-extended state is the most stable in solution. Intermolecular force field parameters were derived from a standard supermolecule approach: we considered the ab initio energies and geometries of a water molecule interacting with each Tc analog at several different positions. The final, rms deviation between the ab initio and force field energies, averaged over all forms, was 0.35 kcal/mol. Intramolecular parameters were adopted from either the standard CHARMM force field, the ab initio structure, or the earlier, plain Tc force field. The model reproduces the ab initio geometry and flexibility of each Tc. As tests, we describe MD and free energy simulations of a solvated complex between three Tcs and the Tet repressor protein.     

Novel quantum mechanical/molecular mechanical method combined with the theory of energy representation: free energy calculation for the Beckmann rearrangement promoted by proton transfers in the supercritical water

The Beckmann rearrangement of acetone oxime promoted by proton transfers in the supercritical water has been investigated by means of the hybrid quantum mechanical/molecular mechanical approach combined with the theory of energy representation (QM/MM-ER) recently developed. The transition state (TS) structures have been explored by ab initio calculations for the reaction of hydrated acetone oxime on the assumption that the reaction is catalyzed by proton transfers along the hydrogen bonds connecting the solute and the solvent water molecules. Up to two water molecules have been considered as reactants that take part in the proton transfers. As a result of the density functional theory calculations with B3LYP functional and aug-cc-pVDZ basis set, it has been found that participation of two water molecules in the reaction reduces the activation free energy by -12.3 kcal/mol. Furthermore, the QM/MM-ER simulations have revealed that the TS is more stabilized than the reactant state in the supercritical water by 2.7 kcal/mol when two water molecules are involved in the reaction. Solvation free energies of the reactant and the TS have been decomposed into terms due to the electronic polarization of the solute, electron density fluctuation, and others to elucidate the origin of the stabilization of the TS as compared with the reactant. It has been revealed that the promotion of the chemical reaction due to the hydration mainly originates from the interaction between the nonpolarized solute and the solvent water molecules at the supercritical state.     

Potential energy functions for an intramolecular proton transfer reaction in the ground and excited state

We describe the development of empirical potential functions for the study of the excited state intramolecular proton transfer reaction in 1-(trifuloroacetylamino)-naphtaquinone (TFNQ). The potential is a combination of the standard CHARMM27 force field for the backbone structure of TFNQ and an empirical valence bond formalism for the proton transfer reaction. The latter is parameterized to reproduce the potential energies both in the ground and the excited state, determined at the CASPT2 level of theory. Parameters describing intermolecular interactions are fitted to reproduce molecular dipole moments computed at the CASSCF level of theory and to reproduce ab initio hydrogen bonding energies and geometries for TFNQ-water bimolecular complexes. The utility of this potential energy function was examined by computing the potentials of mean force for the proton transfer reactions in the gas phase and in water, in both electronic states. The ground state PMF exhibits little solvent effects, whereas computed potential of mean force shows a solvent stabilization of 2.5 kcal mol⁻¹ in the product state region, suggesting proton transfer is more pronounced in polar solvents, consistent with experimental findings. Electronic Supplementary Material The online version of this article (doi:) contains supplementary material, which is available to authorized users. Contribution to the Fernando Bernardi Memorial Issue.     

Ab initio study of chiral recognition in the propylene imine.hydrogen peroxide complex

High level ab initio calculations were employed to study the chiral recognition effect in several chiral molecular pairs that consist of the propylene imine and hydrogen peroxide molecules. The potential energy surfaces for the complexes formed between S-cis-1,2-propylene imine and the two enantiomeric forms of hydrogen peroxide were constructed, using the calculated interaction energies at different separations and orientations. The energy calculations were done using the MOLPRO suite of programs with CCSD(T)/cc-pVDZ. The energies were counterpoise corrected at every point to eliminate the basis set superposition error. Complete geometry optimizations were further carried out for the molecular complexes consisting of the cis- or trans-propylene imine isomers and the two enantiomeric forms of hydrogen peroxide. The geometry optimizations were done using the Gaussian 98 and 03 suites of programs, with MP2/aug-cc-pVDZ being the highest level used. Altogether, eight stable complexes were identified, and the corresponding dissociation energies were calculated with MP2/aug-cc-pVTZ. The largest chirodiastaltic energy is found at 0.74 kcal mol(-1) for the (syn)trans-propylene imine.hydrogen peroxide complexes, where hydrogen peroxide acts as a hydrogen donor and is on the opposite side of the ring from the methyl group. The rotational constants, dipole moments, and harmonic frequencies of the complexes are presented to assist future spectroscopic investigations.     

Active site structures and the redox properties of blue copper proteins: atomic resolution structure of azurin II and electronic structure calculations of azurin, plastocyanin and stellacyanin

Understanding how the active site structures of blue copper proteins determine their redox properties is the central structure-function relationship question of this important class of protein, also referred to as cupredoxins. We here describe both experimental and computational studies of azurin, plastocyanin and stellacyanin designed to define more accurately the geometric structures of the active site of the reduced and oxidized species, and thus to understand how these structures determine the redox potentials of these proteins. To this end the crystal structure of reduced azurin II has been determined at an atomic resolution of 1.13 Angstrom and is presented here. Co-ordinates and structure factors have been deposited in the RCSB Protein Data Bank with accession codes 2ccw and r2ccwsf respectively. The improved accuracy provided by the atomic resolution for the metal stereochemistry are utilised in conjunction with the EXAFS data for theoretical calculations. Multilevel calculations involving density functional theory and molecular mechanical potentials are used to predict both the geometric and electronic structure of the active sites of azurin, plastocyanin and stellacyanin and to estimate the relative redox potentials of these three proteins. We have also compared the relative energies of the structures obtained from experiment at varying resolutions, and from the isolated and embedded cluster calculations. We find significant energy differences between low and high (atomic) resolution structures arising primarily due to inaccuracies in the Cu-ligand distances in the lower resolution structures, emphasising the importance of accurate, very high resolution structural information. QM/MM structures are only approximately 1 kcal mol(-1) lower in energy than the 1.13 Angstrom structure while the optimized gas phase structure is 13.0 kcal mol(-1) lower in energy.     

Evaluation of the charge penetration energy between non-orthogonal molecular orbitals using the Spherical Gaussian Overlap approximation

An overlap dependent formula for evaluating the charge penetration energy between non-orthogonal molecular orbitals is derived using the Spherical Gaussian Overlap approximation. When combined with an accurate multipole representation of the electrostatic energy, such as in the effective fragment potential method, ab initio electrostatic energies are generally reproduced to within 0.2 kcal/mol for a variety of molecular dimers and basis sets. The only larger error is for the DMSO dimer, where the electrostatic energy is overestimated by 0.7 kcal/mol.     

Comparison of charge models for fixed-charge force fields: small-molecule hydration free energies in explicit solvent

In molecular simulations with fixed-charge force fields, the choice of partial atomic charges influences numerous computed physical properties, including binding free energies. Many molecular mechanics force fields specify how nonbonded parameters should be determined, but various choices are often available for how these charges are to be determined for arbitrary small molecules. Here, we compute hydration free energies for a set of 44 small, neutral molecules in two different explicit water models (TIP3P and TIP4P-Ew) to examine the influence of charge model on agreement with experiment. Using the AMBER GAFF force field for nonbonded parameters, we test several different methods for obtaining partial atomic charges, including two fast methods exploiting semiempirical quantum calculations and methods deriving charges from the electrostatic potentials computed with several different levels of ab initio quantum calculations with and without a continuum reaction field treatment of solvent. We find that the best charge sets give a root-mean-square error from experiment of roughly 1 kcal/mol. Surprisingly, agreement with experimental hydration free energies does not increase substantially with increasing level of quantum theory, even when the quantum calculations are performed with a reaction field treatment to better model the aqueous phase. We also find that the semiempirical AM1-BCC method for computing charges works almost as well as any of the more computationally expensive ab initio methods and that the root-mean-square error reported here is similar to that for implicit solvent models reported in the literature. Further, we find that the discrepancy with experimental hydration free energies grows substantially with the polarity of the compound, as does its variation across theory levels.     

Ab initio force field for simulations of proteins and nucleic acids

We present a new self-consistent set of ab initio analytical pair potential to predict specific nonbonded interactions of protein with nucleic acid, of protein with protein, and of nucleic acid with nucleic acid. The purpose of this study is to represent the interaction between biological molecules with an accuracy equivalent to the ab initio molecular orbital calculations, which are used as reference data to obtain the pair potentials. Atoms in nucleic acids and proteins are classified according to their chemical environments. An “effective charge,” a modification of a charge obtained from the Mulliken population analysis, is introduced and used to represent the electrostatic energy. More than 30,000 SCF interaction energies have been calculated to provide the reference data for the fitting procedure that we have adopted in the parameterization of the potentials. The standard deviation is 1.61 kcal/mol for interaction energies spanning the range from about ?220 kcal/mol to +20 kcal/mol. Molecular dynamics simulations, using the above new set of force field, have been performed successfully for the systems where adequate treatments of specific interactions are required: The stability of α-helix of C-peptide and the interaction of spermine with oligonucleotide are examined as preliminary examples.