VB/MM--the validity of the underlying approximations

Two recently developed methods, VB/MM and density embedded VB/MM (DE-VB/MM), are compared, and their respective approximations are examined. The two methods combine valence-bond (VB) calculations with molecular mechanics (MM) and aim to allow VB analysis of reactions in large biological environments. Furthermore, the two methods utilize two major approximations regarding both the overlap and the reduced resonance integral between the various VB configurations. The difference between the two methods, however, is that VB/MM employs these approximations for the overall interaction of the reacting fragments with their surrounding, whereas DE-VB/MM employs the approximations only with regards to the van der Waals (VdW) interactions whereas the electrostatic interactions are calculated rigorously at the quantum level. The approximations that lay the grounds for the two methods involve the assumption that the overlap between the VB configurations and the respective reduced resonance integral are both invariant to the environment. Similar approximations are utilized in several other VB-based QM/MM methods. However, although extensively used, these approximations were never rigorously proved. Here, we exploit the development of the DE-VB/MM method to numerically examine the approximations by calculating the accurate as well as the approximated values of overlap and reduced resonance integral for systems where the environment involves only electrostatic interactions. The quality of the approximations is examined together with their effect on the absolute energies, the wave function, and the overall energetics. Three test cases are chosen, the dissociation of CH 3F and LiF and the identity S N2 reaction. It is shown that the approximations are usually good with the exception of cases where extreme changes are expected in the wave function. Furthermore, the impact of the approximations on the overall wave function and the overall energetics is found to be quite small. It is concluded that VB/MM, where the approximations are used more extensively, can serve as the first method of choice.     

Hybrid ab initio VB/MM method--a valence bond ride through classical landscapes

The development of a new hybrid (QM/MM) method, where the QM part is treated by ab initio valence bond (VB) theory is presented. This VB/MM method has the advantages of empirical VB (EVB) methodology but does not rely on empirical parameterization for the quantum part. The method implements embedding of the quantum region of each diabatic state separately, by treating the electrostatic interactions between QM and MM regions classically. Additionally, it assumes that changes of the off diagonal matrix element due to different environments are such that the overall resonance integral does not change. These assumptions are discussed in detail and the validity of the method is shown to be successful in three different bond dissociation processes, where bond dissociation as well as solvation energies compare very well with the experimental data.     

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.     

Sr2+和Ba2+水分子体系结构和结合能的从头算和ABEEM/MM研究

应用从头算方法和ABEEM/MM浮动电荷分子力场, 研究了水合碱土离子团簇Sr²⁺/Ba²⁺(H₂O)_n (n=1-6), 构建了离子-水相互作用的ABEEM/MM势能函数, 获得了水合离子团簇的稳定结构, 计算了结合能. 计算结果表明, ABEEM/MM方法的结果和从头算方法的结果有很好的一致性. 进一步应用ABEEM/MM对Sr²⁺和Ba²⁺水溶液进行了分子动力学模拟. 对Sr²⁺水溶液, 得到的Sr²⁺-水中氧原子的径向分布函数的第一和第二最高峰分别位于0.257和0.464 nm处, 第一和第二水合层的配位水分子数分别为9.2和11.4; 对Ba²⁺水溶液, 得到的Ba²⁺与水中氧原子的径向分布函数的第一和第二最高峰分别位于0.269和0.467 nm处, 第一和第二水合层的配位水分子数分别为9.9和12.4. 这与实验值或其它理论模拟结果有较好的一致性. 对比外层的水分子, 金属离子的极化作用使得溶液中第一水合层中水分子的O―H键长增长, HOH键角减小.     

Role of the catalytic triad and oxyanion hole in acetylcholinesterase catalysis: an ab initio QM/MM study

The initial step of the acylation reaction catalyzed by acetylcholinesterase (AChE) has been studied by a combined ab initio quantum mechanical/molecular mechanical (QM/MM) approach. The reaction proceeds through the nucleophilic addition of the Ser203 O to the carbonyl C of acetylcholine, and the reaction is facilitated by simultaneous proton transfer from Ser203 to His447. The calculated potential energy barrier at the MP2(6-31+G) QM/MM level is 10.5 kcal/mol, consistent with the experimental reaction rate. The third residue of the catalytic triad, Glu334, is found to be essential in stabilizing the transition state through electrostatic interactions. The oxyanion hole, formed by peptidic NH groups from Gly121, Gly122, and Ala204, is also found to play an important role in catalysis. Our calculations indicate that, in the AChE-ACh Michaelis complex, only two hydrogen bonds are formed between the carbonyl oxygen of ACh and the peptidic NH groups of Gly121 and Gly122. As the reaction proceeds, the distance between the carbonyl oxygen of ACh and NH group of Ala204 becomes smaller, and the third hydrogen bond is formed both in the transition state and in the tetrahedral intermediate.     

Producing DFT/MM enzyme reaction trajectories from SCC‐DFTB/MM driving forces to probe the underlying electronics of a glycosyltransferase reaction

The SCC‐DFTB/MIO/CHARMM free energy surface for a glycosyltransferase, TcTS, is benchmarked against a DFT/MM reaction trajectory using the same CHARMM MM force field ported to the NWChem package. The popular B3LYP functional, against which the MIO parameter set was parameterized is used to optimize TS structures and run DFT reaction dynamics. A novel approach was used to generate reaction forces from a SCC‐DFTB/MIO/CHARMM reaction surface to drive B3LYP/6‐31G/MM and B3LYP/6‐31G(d)/MM reaction trajectories. Although TS structures compare favorably, differences stemming primarily from a minimal basis set approximation prevented a successful 6‐31G(d) FEARCF reaction dynamics trajectory. None the less, the dynamic evolution of the B3LYP/6‐31G/MM‐computed electron density provided an opportunity to perform NBO analysis along the reaction trajectory. Here, we illustrate that a successful ab initio reaction trajectory is computationally accessible when the underlying potential energy function of the semi‐empirical method used to produce driving forces is sufficiently close to the ab initio potential. © 2017 Wiley Periodicals, Inc.     

Prediction and rationalization of protein pKa values using QM and QM/MM methods

We describe the development and application of a computational method for the prediction and rationalization of pKa values of ionizable residues in proteins, based on ab initio quantum mechanics (QM) and the effective fragment potential (EFPs) method (a hybrid QM/MM method). The theoretical developments include (1) a covalent boundary method based on frozen localized orbitals, (2) divide-and-conquer methods for the ab initio computation of protein EFPs consisting of multipoles up to octupoles and dipole polarizability tensors, (3) a method for computing vibrational free energies for a localized molecular region, and (4) solutions of the polarized continuum model of bulk solvation equations for protein-sized systems. The QM-based pKa prediction method is one of the most accurate methods currently available and can be used in cases where other pKa prediction methods fail. Preliminary analysis of the computed results indicate that many pKa values (1) are primarily determined by hydrogen bonds rather than long-range charge-charge interactions and (2) are relatively insensitive to large-scale dynamical fluctuations of the protein structure.     

Flexible effective fragment QM/MM method: validation through the challenging tests

A new version of the QM/MM method, which is based on the effective fragment potential (EFP) methodology [Gordon, M. et al., J Phys Chem A 2001, 105, 293] but allows flexible fragments, is verified through calculations of model molecular systems suggested by different authors as challenging tests for QM/MM approaches. For each example, the results of QM/MM calculations for a partitioned system are compared to the results of an all-electron ab initio quantum chemical study of the entire system. In each case we were able to achieve approximately similar or better accuracy of the QM/MM results compared to those described in original publications. In all calculations we kept the same set of parameters of our QM/MM scheme. A new test example is considered when calculating the potential of internal rotation in the histidine dipeptide around the C(alpha)bond;C(beta) side chain bond.     

Valence bond calculations of hydrogen transfer reactions: a general predictive pattern derived from theory

Hydrogen abstraction reactions of the type X(*) + H-H' --> X-H + H'(*) (X = F, Cl, Br, I) are studied by ab initio valence bond methods and the VB state correlation diagram (VBSCD) model. The reaction barriers and VB parameters of the VBSCD are computed by using the breathing orbital valence bond and valence bond configuration interaction methods. The combination of the VBSCD model and semiempirical VB theory leads to analytical expressions for the barriers and other VB quantities that match the ab initio VB calculations fairly well. The barriers are influenced by the endo- or exothermicity of the reaction, but the fundamental factor of the barrier is the average singlet-triplet gap of the bonds that are broken or formed in the reactions. Some further approximations lead to a simple formula that expresses the barrier for nonidentity and identity hydrogen abstraction reactions as a function of the bond strengths of reactants and products. The semiempirical expressions are shown to be useful not only for the model reactions that are studied in this work, but also for other nonidentity and identity hydrogen abstraction reactions that have been studied in previous articles.     

How does the cAMP-dependent protein kinase catalyze the phosphorylation reaction: an ab initio QM/MM study

We have carried out density functional theory QM/MM calculations on the catalytic subunit of cAMP-dependent protein kinase (PKA). The QM/MM calculations indicate that the phosphorylation reaction catalyzed by PKA is mainly dissociative, and Asp166 serves as the catalytic base to accept the proton delivered by the substrate peptide. Among the key interactions in the active site, the Mg(2+) ions, glycine rich loop, and Lys72 are found to stabilize the transition state through electrostatic interactions. On the other hand, Lys168, Asn171, Asp184, and the conserved waters bound to Mg(2+) ions do not directly contribute to lower the energy barrier of the phosphorylation reaction, and possible roles for these residues are proposed. The QM/MM calculations with different QM/MM partition schemes or different initial structures yield consistent results. In addition, we have carried out 12 ns molecular dynamics simulations on both wild type and K168A mutated PKA, respectively, to demonstrate that the catalytic role of Lys168 is to keep ATP and substrate peptide in the near-attack reactive conformation.     

On the construction of diabatic and adiabatic potential energy surfaces based on ab initio valence bond theory

A theoretical model is presented for deriving effective diabatic states based on ab initio valence bond self-consistent field (VBSCF) theory by reducing the multiconfigurational VB Hamiltonian into an effective two-state model. We describe two computational approaches for the optimization of the effective diabatic configurations, resulting in two ways of interpreting such effective diabatic states. In the variational diabatic configuration (VDC) method, the energies of the diabatic states are variationally minimized. In the consistent diabatic configuration (CDC) method, both the configuration coefficients and orbital coefficients are simultaneously optimized to minimize the adiabatic ground-state energy in VBSCF calculations. In addition, we describe a mixed molecular orbital and valence bond (MOVB) approach to construct the CDC diabatic and adiabatic states for a chemical reaction. Note that the VDC-MOVB method has been described previously. Employing the symmetric S(N)2 reaction between NH(3) and CH(3)NH(3)(+) as a test system, we found that the results from ab initio VBSCF and from ab initio MOVB calculations using the same basis set are in good agreement, suggesting that the computationally efficient MOVB method is a reasonable model for VB simulations of condensed phase reactions. The results indicate that CDC and VDC diabatic states converge, respectively, to covalent and ionic states as the molecular geometries are distorted from the minimum of the respective diabatic state along the reaction coordinate. Furthermore, the resonance energy that stabilizes the energy of crossing between the two diabatic states, resulting in the transition state of the adiabatic ground-state reaction, has a strong dependence on the overlap integral between the two diabatic states and is a function of both the exchange integral and the total diabatic ground-state energy.     

Molecular mechanics (MM3) calculations on lithium amide compounds

The MM3 force field has been extended to deal with the lithium amide molecules that are widely used as efficient catalysts for stereoselective asymmetric synthesis. The MM3 force field parameters have been determined on the basis of the ab initio MP2/6-31G* and/or DFT (B3LYP/6-31G*, B3-PW91/6-31G*) geometry optimization calculations. To evaluate the electronic interactions specific to the lithium amides derived from the diamine molecules properly, the Lewis bonding potential term for the interaction between the lithium atom and the nonbonded adjacent electronegative atom such as nitrogen was introduced into the MM3 force field. The bond dipoles were evaluated correctly from the electronic charges on the atoms calculated by fitting to the electrostatic potential at points selected. The MM3 results on the molecular structures, conformational energies, and vibrational spectra show good agreement with those from the quantum mechanical calculations.     

Ab initio quantum mechanical charge field (QMCF) molecular dynamics: a QM/MM – MD procedure for accurate simulations of ions and complexes

A new formalism for quantum mechanical / molecular mechanical (QM/MM) dynamics of chemical species in solution has been developed, which does not require the construction of any other potential functions except those for solvent–solvent interactions, maintains all the advantages of large simulation boxes and ensures the accuracy of ab initio quantum mechanics for all forces acting in the chemically most relevant region. Interactions between solute and more distant solvent molecules are incorporated by a dynamically adjusted force field corresponding to the actual molecular configuration of the simulated system and charges derived from the electron distribution in the solvate. The new formalism has been tested with some examples of hydrated ions, for which accurate conventional ab initio QM/MM simulations have been previously performed, and the comparison shows equivalence and in some aspects superiority of the new method. As this simulation procedure does not require any tedious construction of two-and three-body interaction potentials inherent to conventional QM/MM approaches, it opens the straightforward access to ab initio molecular dynamics simulations of any kind of solutes, such as metal complexes and other composite species in solution.     

Pseudobond ab initio QM/MM approach and its applications to enzyme reactions

This perspective article mainly focuses on the development and applications of a pseudobond ab initio QM/MM approach to study enzyme reactions. The following aspects of methodology development are discussed: the approaches for the QM/MM covalent boundary problem, an efficient iterative optimization procedure, the methods to determine enzyme reaction paths, and the approaches to calculate free energy change in enzyme reactions. Several applications are described to illustrate the capability of the methods. Finally, future directions are discussed.     

Structure of an accurate ab initio model of the aqueous Cl- ion at high temperatures

The structure of an accurate ab initio model of aqueous chloride ion was calculated at two high-temperature state points (573 K, 0.725 g/cm(3) and 723 K, 0.0098 g/cm(3)) by a two-step procedure. First, the structure of an approximate model was calculated from a molecular dynamics simulation of the model. Then the difference between the structure of the ab initio model and the approximate model was calculated by non-Boltzmann weighting of a sample of configurations taken from the approximate model simulation. Radial distribution functions, average coordination numbers, the distribution of coordination numbers, an analysis of orientations of water in the first coordination shell, and the free energy of hydration of the chloride ion are reported for both state points. The most common water structure has one hydrogen close to the chloride ion and one pointing away (46% at 573 K and 57% at 723 K). Waters in the first coordination shell that are not strongly bound to the chloride ions are common. Several variations of the method were tested. Models in which the water-water interaction is calculated with ab initio methods predict only a slightly different structure than models in which water-water interactions are determined from the approximate models. Similarly, using the approximate model for solute-water interactions when the water is far from the chloride ion did not affect the results. Uncertainties due to the limited sample of configurations are estimated and found to be small. The results are in qualitative agreement with X-ray and neutron diffraction experiments and with simulations of approximate models.