Evaluation of an ab initio quantum mechanical/molecular mechanical hybrid-potential link-atom method

 Hybrid potentials have become a common tool in the study of many condensed-phase processes and are the subject of much active research. An important aspect of the formulation of a hybrid potential concerns how to handle covalent bonds between atoms that are described with different potentials and, most notably, those at the interface of the quantum mechanical (QM) and molecular mechanical (MM) regions. Several methods have been proposed to deal with this problem, ranging from the simple link-atom method to more sophisticated hybrid-orbital techniques. Although it has been heavily criticized, the link-atom method has probably been the most widely used in applications, especially with hybrid potentials that use semiempirical QM methods. Our aim in this paper has been to evaluate the link-atom method for ab initio QM/MM hybrid potentials and to compare the results it gives with those of previously published studies. Given its simplicity and robustness, we find that the link-atom method can produce results of comparable accuracy to other methods as long as the charge distribution on the MM atoms at the interface is treated appropriately. Received: 27 September 2002 / Accepted: 21 October 2002 / Published online: 8 January 2003 Correspondence to: M. J. Field e-mail: mjfield@ibs.fr Acknowledgements. The authors thank the Institut de Biologie Structurale – Jean-Pierre Ebel, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique for support of this work.     

Exploring the quantum mechanical/molecular mechanical replica path method: a pathway optimization of the chorismate to prephenate Claisen rearrangement catalyzed by chorismate mutase

 A replica path method has been developed and extended for use in complex systems involving hybrid quantum/classical (quantum mechanical/molecular mechanical) coupled potentials. This method involves the definition of a reaction path via replication of a set of macromolecular atoms. An “important” subset of these replicated atoms is restrained with a penalty function based on weighted root-mean-square rotation/translation best-fit distances between adjacent (i±1) and next adjacent (i±2) pathway steps. An independent subset of the replicated atoms may be treated quantum mechanically using the computational engine Gamess-UK. This treatment can be performed in a highly parallel manner in which many dozens of processors can be efficiently employed. Computed forces may be projected onto a reference pathway and integrated to yield a potential of mean force (PMF). This PMF, which does not suffer from large errors associated with calculated potential-energy differences, is extremely advantageous. As an example, the QM/MM replica path method is applied to the study of the Claisen rearrangement of chorismate to prephenate which is catalyzed by the Bacillus subtilis isolated, chorismate mutase. Results of the QM/MM pathway minimizations yielded an activation enthalpy ΔH ^†† of 14.9 kcal/mol and a reaction enthalpy of −19.5 kcal/mol at the B3LYP/6-31G(d) level of theory. The resultant pathway was compared and contrasted with one obtained using a forced transition approach based on a reaction coordinate constrained repeated walk procedure (ΔH ^†† =20.1 kcal/mol, ΔH _rxn = −20.1 kcal/mol, RHF/4-31G). The optimized replica path results compare favorably to the experimental activation enthalpy of 12.7±0.4 kcal/mol. Received: 16 December 2001 / Accepted: 6 September 2002 / Published online: 8 April 2003 Contribution to the Proceedings of the Symposium on Combined QM/MM Methods at the 22nd National Meeting of the American Chemical Society, 2001. Correspondence to: H.L. Woodcock e-mail: hlwood@ccqc.uga.edu Acknowledgements. The authors thank Eric Billings, Xiongwu Wu, and Stephen Bogusz for helpful discussions and related work. The authors also show grateful appreciation to The National Institutes of Health and The National Science Foundation for support of the current research.     

The search for stationary points on a quantum mechanical/molecular mechanical potential-energy surface

 We propose a methodology to locate stationary points on a quantum mechanical/molecular mechanical potential-energy surface. This algorithm is based on a suitable approximation of an initial full Hessian matrix, either a modified Broyden–Fletcher–Goldfarg–Shanno or a Powell update formula for the location of, respectively, a minimum or a transition state, and the so-called rational function optimization. The latter avoids the Hessian matrix inversion required by a quasi-Newton–Raphson method. Some examples are presented and analyzed. Received: 16 July 2001 / Accepted: 9 October 2001 / Published online: 9 January 2002     

Generalized hybrid orbital for the treatment of boundary atoms in combined quantum mechanical and molecular mechanical calculations using the semiempirical parameterized model 3 method

The application of combined quantum mechanical (QM) and molecular mechanical methods to large molecular systems requires an adequate treatment of the boundary between the two approaches. In this article, we extend the generalized hybrid orbital (GHO) method to the semiempirical parameterized model 3 (PM3) Hamiltonian combined with the CHARMM force field. The GHO method makes use of four hybrid orbitals, one of which is included in the QM region in self-consistent field optimization and three are treated as auxiliary orbitals that do not participate in the QM optimization, but they provide an effective electric field for interactions. An important feature of the GHO method is that the semiempirical parameters for the boundary atom are transferable, and these parameters have been developed for a carbon boundary atom consistent with the PM3 model. The combined GHO-PM3/CHARMM model has been tested on molecular geometry and proton affinity for a series of organic compounds.Acknowledgement We thank the National Institutes of Health for support of this research.Contribution to the Jacopo Tomasi Honorary Issue     

Recent advances in quantum mechanical/molecular mechanical calculations of enzyme catalysis: hydrogen tunnelling in liver alcohol dehydrogenase and inhibition of elastase by α-ketoheterocycles

 Hybrid quantum mechanical (QM)/molecular mechanical (MM) calculations are used to study two aspects of enzyme catalysis, Kinetic isotope effects associated with the hydride ion transfer step in the reduction of benzyl alcohol by liver alcohol dehydrogenase are studied by employing variational transition-state theory and optimised multidimensional tunnelling. With the smaller QM region, described at the Hartree–Fock ab initio level, together with a parameterised zinc atom charge, good agreement with experiment is obtained. A comparison is made with the proton transfer in methylamine dehydrogenase. The origin of the large range in pharmacological activity shown by a series of α-ketoheterocycle inhibitors of the serine protease, elastase, is investigated by both force field and QM/MM calculations. Both models point to two different inhibition mechanisms being operative. Initial QM/MM calculations suggest that these are binding, and reaction to form a tetrahedral intermediate, the latter process occurring for only the more potent set of inhibitors. Recieved 3 October 2001 / Accepted: 6 September 2002 / Published online: 31 January 2003 Contribution to the Proceedings of the Symposium on Combined QM/MM Methods at the 222nd National Meeting of the American Chemical Society, 2001 Correspondence to: I. H. Hillier Acknowledgements. We thank EPSRC and BBSRC for support of the research and D.G. Truhlar for the use of the POLYRATE code.     

Development of a model for the rational design of molecular imprinted polymer: Computational approach for combined molecular dynamics/quantum mechanics calculations

A new rational approach for the preparation of molecularly imprinted polymer (MIP) based on the combination of molecular dynamics (MD) simulations and quantum mechanics (QM) calculations is described in this work. Before performing molecular modeling, a virtual library of functional monomers was created containing forty frequently used monomers. The MD simulations were first conducted to screen the top three monomers from virtual library in each porogen-acetonitrile, chloroform and carbon tetrachloride. QM simulations were then performed with an aim to select the optimum monomer and progen solvent in which the QM simulations were carried out; the monomers giving the highest binding energies were chosen as the candidate to prepare MIP in its corresponding solvent. The acetochlor, a widely used herbicide, was chosen as the target analyte. According to the theoretical calculation results, the MIP with acetochlor as template was prepared by emulsion polymerization method using N,N-methylene bisacrylamide (MBAAM) as functional monomer and divinylbenzene (DVB) as cross-linker in chloroform. The synthesized MIP was then tested by equilibrium-adsorption method, and the MIP demonstrated high removal efficiency to the acetochlor. Mulliken charge distribution and ¹H NMR spectroscopy of the synthesized MIP provided insight on the nature of recognition during the imprinting process probing the governing interactions for selective binding site formation at a molecular level. We think the computer simulation method first proposed in this paper is a novel and reliable method for the design and synthesis of MIP.     

Molecular recognition by β-cyclodextrin derivatives: molecular dynamics, free-energy perturbation and molecular mechanics/ Poisson–Boltzmann surface area goals and problems

 The complexation of p-tert-butylphenyl p-tert-butylbenzoate and N-(p-tert-butylphenyl)-p-tert-butylbenzamide with a β-cyclodextrin derivative formed by two cyclodextrin units linked by a disulfide bridge on one of the C6 atoms has been studied by computational methods. The better amide solubility and the better internal interactions of the ester complex explain the experimentally observed better association constant for the ester. The free-energy perturbation methodology and molecular mechanics/Poisson–Boltzmann surface area analysis have been used to explain the problem and to compare the results. Received: 14 April 2002 / Accepted: 11 August 2002 / Published online: 4 November 2002 Acknowledgements. The Kollman group at the University of California San Francisco is gratefully acknowledged for support and encouragement throughout all this study. The authors thank UAB for inland and outland fellowships to I.B.. Financial support was obtained from grant no. PPQ2000-0369 from the “Ministerio de Ciencia y Tecnologia” (Spain). Intensive computations were performed either with the computers of the Kollman group or with those of CESCA-C⁴ (Catalonia, Spain). Correspondence to: C. Jaime e-mail: carlos.jaime@uab.es     

Implementation of the IMOMM methodology for performing combined QM/MM molecular dynamics simulations and frequency calculations

A technique for implementing the integrated molecular orbital and molecular mechanics (IMOMM) methodology developed by Maseras and Morokuma that is used to perform combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations, frequency calculations and simulations of macromolecules including explicit solvent is presented. Although the IMOMM methodology is generalized to any coordinate system, the implementation first described by Maseras and Morokuma requires that the QM and MM gradients be transformed into internal coordinates before they are added together. This coordinate transformation can be cumbersome for macromolecular systems and can become ill-defined during the course of a molecular dynamics simulation. We describe an implementation of the IMOMM method in which the QM and MM gradients are combined in the cartesian coordinate system, thereby avoiding potential problems associated with using the internal coordinate system. The implementation can be used to perform combined QM/MM molecular dynamics simulations and frequency calculations within the IMOMM framework. Finally, we have examined the applicability of thermochemical data derived from IMOMM framework. Finally, we have examined the applicability of thermochemical data derived from IMOMM frequency calculations. Received: 11 May 1998 / Accepted: 14 August 1998 / Published online: 16 November 1998     

A combined quantum mechanical and molecular mechanical method using modified generalized hybrid orbitals: implementation for electronic excited states

The generalized hybrid orbital (GHO) method is implemented at the second-order approximate coupled cluster singles and doubles (CC2) level for quantum mechanical (QM)/molecular mechanical (MM) electronic excited state calculations. The linear response function of CC2 in the GHO scheme is derived and implemented. The new implementation is applied to the first singlet excited states of three aromatic amino acids, phenylalanine, tyrosine, and tryptophan, and also bacteriorhodopsin for assessment. The results obtained for aromatic amino acids agreed well with the full QM CC2 calculations, while the calculated excitation energies of bacteriorhodopsin and its chromophore, all-trans retinal, reproduced the environmental shift of the experimental data. For the bacteriorhodopsin case, the environmental shift of GHO also showed good agreements with the experimental data. The contribution of the quantum effect of certain moieties in the excited states is elucidated by changing the partitioning of QM and MM regions.     

Poisson-Boltzmann calculations versus molecular dynamics simulations for calculating the electrostatic potential of a solvated peptide

We performed a very long molecular dynamics simulation of a peptide in explicit water molecules and ions and averaged the electrostatic potential caused by peptide, water and ions at eight points in the vicinity of the peptide. These electrostatic potential values were directly compared to the potential calculated by solving the non-linear Poisson-Boltzmann equation for the system, which describes the solvent using continuum electrostatics. We analyze the contribution of dielectric constant, conformational flexibility and solvation effects on the electrostatic potential at these eight points. Received: 24 April 1998 / Accepted: 4 August 1998 / Published online: 23 November 1998     

Quantum mechanics and mixed quantum mechanics/molecular mechanics simulations of model nerve agents with acetylcholinesterase

 The accurate modeling of biological processes presents major computational difficulties owing to the inherent complexity of the macromolecular systems of interest. Simulations of biochemical reactivity tend to require highly computationally intensive quantum mechanical methods, but localized chemical effects tend to depend significantly on properties of the extended biological environment – a regime far more readily examined with lower-level classical empirical models. Mixed quantum/classical techniques are gaining in popularity as a means of bridging these competing requirements. Here we present results comparing two quantum mechanics/molecular mechanics implementations (the SIMOMM technique of Gordon et al. as implemented in GAMESS, and the ONIOM technique of Morokuma et al. found in Gaussian 98) as performed on the enzyme acetylcholinesterase and model nerve agents. This work represents part of the initial phase of a DoD HPCMP Challenge project in which we are attempting to reliably characterize the biochemical processes responsible for nerve agent activity and inhibition, thereby allowing predictions on compounds unrelated to those already studied. Received: 10 October 2001 / Accepted: 13 November 2002 / Published online: 1 April 2003 Contribution to the Proceedings of the Symposium on Combined QM/MM Methods at the 222nd National Meeting of the American Chemical Society, 2001 Correspondence to: M. M. Hurley e-mail: hurley@arl.army.mil     

The coordination of the catalytic zinc ion in alcohol dehydrogenase studied by combined quantum-chemical and molecular mechanics calculations

Summary The coordination number of the catalytic zinc ion in alcohol dehydrogenase has been studied by integrated ab initio quantum-chemical and molecular mechanics geometry optimisations involving the whole enzyme. A four-coordinate active-site zinc ion is 100–200 kJ/mol more stable than a five-coordinate one, depending on the ligands. The only stable binding site for a fifth ligand at the zinc ion is opposite to the normal substrate site, in a small cavity buried behind the zinc ion. The zinc coordination sphere has to be strongly distorted to accommodate a ligand in this site, and the ligand makes awkward contacts with surrounding atoms. Thus, the results do not support proposals attributing an important role to five-coordinate zinc complexes in the catalytic mechanism of alcohol dehydrogenase. The present approach makes it possible also to quantify the strain induced by the enzyme onto the zinc ion and its ligands; it amounts to 42–87 kJ/mol for four-coordinate active-site zinc ion complexes and 131–172 kJ/mol for five-coordinate ones. The four-coordinate structure with a water molecule bound to the zinc ion is about 20 kJ/mol less strained than the corresponding structure with a hydroxide ion, indicating that the enzyme does not speed up the reaction by forcing the zinc coordination sphere into a structure similar to the reaction intermediates.     

Aggregation mechanism of polyglutamine diseases revealed using quantum chemical calculations, fragment molecular orbital calculations, molecular dynamics simulations, and binding free energy calculations

Polyglutamine (polyQ) diseases, including Huntington’s disease (HD), are caused by expansion of polyQ-encoding repeats within otherwise unrelated gene products. The aggregation mechanism of polyQ diseases, the inhibition mechanism of Congo red, and the alleviation mechanism of trehalose were proposed here based on quantum chemical calculations and molecular dynamics simulations. The calculations and simulations revealed the following. The effective molecular bonding is between glutamine (Gln) and Gln (Gln + Gln), between Gln and Congo red (Gln + Congo red), and between Gln and trehalose (Gln + trehalose). The bonding strength is −13.1 kcal/mol for Gln + Gln, −24.4 kcal/mol for Gln + Congo red, and −12.0 kcal/mol for Gln + trehalose. In the polyQ region, both the number of intermolecular Gln + Gln formations and the total calories generated by the Gln + Gln formation are proportional to the number of repetitions of Gln. We propose an aggregation mechanism whose heat generated by the intermolecular Gln + Gln formation causes the pathogeny of polyQ disease. In our aggregation mechanism, this generated heat collapses the host protein and promotes fibrillogenesis. Without contradiction, our mechanism can explain all the experimental results reported to date. Our mechanism can also explain the inhibition mechanism by Congo red as an inhibitor of polyglutamine-induced protein aggregation and the alleviation mechanism by trehalose as an alleviator of that aggregation. The inhibition mechanism by Congo red is explained by the strong interaction with Gln and by the characteristic structure of Congo red.     

Analytic energy derivatives for the equation-of-motion coupled-cluster method: Algebraic expressions,implementation and application to theS 1 state of HFCO

Summary The general theory of analytic derivatives for the equation-of-motion coupled cluster (EOM-CC) method is reviewed. Special attention is paid to the EOM-CC singles and doubles (EOM-CCSD) approximation, which has the same computational scaling properties as the coupled-cluster singles doubles (CCSD) ground state method and is therefore applicable to a wide range of molecular systems. The detailed spin orbital equations that must be solved in EOM-CCSD gradient calculations are presented for the first time, and some guidelines are discussed regarding their computational implementation. Finally, use of the EOM-CCSD gradient method is illustrated by determining the structure, dipole moment components, harmonic frequencies and infrared intensities of formyl fluoride (HFCO) in its singlet excited (n, *) state.     

Effective modeling of intrinsic and environmental effects on the structure and electron plaramagnetic resonance parameters of nitroxides by an integrated quantum mechanical/molecular mechanics/polarizable continuum model approach

 Experimental and density functional theory geometries have been used to extend the AMBER force field to nitroxides. An optimum set of transferable atomic charges for the calculation of electrostatic interactions both in vacuo and in aqueous solution has been obtained by averaging the charges obtained by a restrained electrostatic potential fitting of representative compounds. Besides reliable structural data, our implementation allows the computation of accurate spectromagnetic properties by single-point B3LYP computations on geometries optimized at the AMBER level. Solvent shifts in aqueous solution can be reproduced quantitatively by a mixed model in which specific solvent effects are described by two water molecules strongly coordinated to the nitroxide oxygen, while bulk effects are described by the polarizable continuum model. Received: 17 September 1999 / Accepted: 3 February 2000 / Published online: 29 June 2000     

Molecular dynamics study of the tautomeric equilibrium in the 4-nitro- and 2,4,6-trichloro derivatives of 2-(N,N-dialkyloaminomethyl)phenol

 Molecular dynamics thermodynamic integration (MDTI) method and quantum chemical calculations at the density functional theory B3LYP 6-31+(d,p) level, which included the Tomasi model of the solvent reaction field, were applied to study the tautomeric equilibrium of Mannich base in methanol solution. The values obtained for the free-energy difference are in good agreement with experimental data. However, the results from quantum mechanical calculations were not as good as the results of MDTI simulations owing to inappropriate treatment of intermolecular hydrogen bonds between the solute molecule and the first shell of solvent molecules in the Tomasi model of the solvent reaction field. The radial distribution functions between solute atoms and solvent atoms confirmed the formation of hydrogen bonds between the solute molecule and surrounding methanol molecules and indicated that the zwitterionic form is associated more with an organized solvent structure at the level of the first solvation shell than is the molecular form. Received: 26 April 2002 / Accepted: 9 September 2002 / Published online: 31 March 2003     

Isotope effect of hydrogen and lithium hydride molecules. Application of the dynamic extended molecular orbital method and energy component analysis

In order to explore the isotope effect including the nuclear–electronic coupling and nuclear quantum effects under the one-particle approximation, we apply the dynamic extended molecular orbital (DEMO) method and energy component analysis to the hydrogen and lithium hydride isotope molecules. Since the DEMO method determines both electronic and nuclear wave functions simultaneously by variationally optimizing all parameters embedded in the basis sets, the virial theorem is completely satisfied and guarantees the relation of the kinetic and potential energies. We confirm the isotope effect on internuclear distances, nuclear and electronic wave functions, dipole moment, the polarizability, and each energy component. In the case of isotopic species of the hydrogen molecule, the total energy decreases from the H₂ to the T₂ molecule due to the stabilization of the nuclear–electronic potential component, as well as the nuclear kinetic one. In the case of the lithium hydride molecule, the energy lowering by replacing ⁶Li with ⁷Li is calculated to be greater than that by replacing H with D. This is mainly caused by the small destabilization of electron–electron and nuclear–nuclear repulsion in ⁷LiH compared to ⁶LiH, while the change in the repulsive components from ⁶LiH to ⁶LiD increases. Received: 24 March 1999 / Accepted: 5 August 1999 / Published online: 15 December 1999     

Investigation of the S<Subscript>0</Subscript>S<Subscript>1</Subscript> excitation in bacteriorhodopsin with the ONIOM(MO:MM) hybrid method

 We have investigated the S₀ and S₁ electronic states in bacteriorhodopsin using a variety of QM/MM levels. The decomposition of the calculated excitation energies into electronic and electrostatic components shows that the interaction of the chromophore with the protein electric field increases the excitation energy, while polarization effects are negligible. Therefore, the experimentally observed reduction in excitation energy from solution phase to protein environment (the Opsin shift) does not come from the electrostatic interaction with the protein environment, but from either the interaction ofthe chromophore with the solvent or counter ion, or structural effects. Our high-level ONIOM(TD– B3LYP:Amber) calculation predicts the excitation energy within 8 kcal/mol from experiment, the discrepancy probably being caused by the neglect of polarization of the protein environment. In addition, we have shown that the level of optimization is extremely critical for the calculation of accurate excitation energies in bacteriorhodopsin. Received: 13 October 2001 / Accepted: 6 September 2002 / Published online: 3 February 2003 Contribution to the Proceedings of the Symposium on Combined QM/MM Methods at the 222nd National Meeting of the American Chemical Society, 2001 Correspondence to: K. Morokuma e-mail: morokuma@emory.edu     

Solvent effects on the electronic absorption spectrum of formamide studied by a sequential Monte Carlo/quantum mechanical approach

Sequential Monte Carlo/quantum mechanical calculations are performed to study the solvent effects on the electronic absorption spectrum of formamide (FMA) in aqueous solution, varying from hydrogen bonds to the outer solvation shells. Full quantum-mechanical intermediate neglect of differential overlap/singly excited configuration interaction calculations are performed in the supermolecular structures generated by the Monte Carlo simulation. The largest calculation involves the ensemble average of 75 statistically uncorrelated quantum mechanical results obtained with the FMA solute surrounded by 150 water solvent molecules. We find that the n → π* transition suffers a blueshift of 1,600 cm⁻¹ upon solvation and the π → π* transition undergoes a redshift of 800 cm⁻¹. On average, 1.5 hydrogen bonds are formed between FMA and water and these contribute with about 20% and about 30% of the total solvation shifts of the n → π* and π → π* transitions, respectively. The autocorrelation function of the energy is used to sample configurations from the Monte Carlo simulation, and the solvation shifts are shown to be converged values. Received: 14 March 2002 / Accepted: 3 April 2002 / Published online: 24 June 2002