An efficient method for the calculation of quantum mechanics/molecular mechanics free energies

The combination of quantum mechanics (QM) with molecular mechanics (MM) offers a route to improved accuracy in the study of biological systems, and there is now significant research effort being spent to develop QM/MM methods that can be applied to the calculation of relative free energies. Currently, the computational expense of the QM part of the calculation means that there is no single method that achieves both efficiency and rigor; either the QM/MM free energy method is rigorous and computationally expensive, or the method introduces efficiency-led assumptions that can lead to errors in the result, or a lack of generality of application. In this paper we demonstrate a combined approach to form a single, efficient, and, in principle, exact QM/MM free energy method. We demonstrate the application of this method by using it to explore the difference in hydration of water and methane. We demonstrate that it is possible to calculate highly converged QM/MM relative free energies at the MP2/aug-cc-pVDZ/OPLS level within just two days of computation, using commodity processors, and show how the method allows consistent, high-quality sampling of complex solvent configurational change, both when perturbing hydrophilic water into hydrophobic methane, and also when moving from a MM Hamiltonian to a QM/MM Hamiltonian. The results demonstrate the validity and power of this methodology, and raise important questions regarding the compatibility of MM and QM/MM forcefields, and offer a potential route to improved compatibility.     

Insights into the different dioxygen activation pathways of methane and toluene monooxygenase hydroxylases

The methane and toluene monooxygenase hydroxylases (MMOH and TMOH, respectively) have almost identical active sites, yet the physical and chemical properties of their oxygenated intermediates, designated P*, H(peroxo), Q, and Q* in MMOH and ToMOH(peroxo) in a subclass of TMOH, ToMOH, are substantially different. We review and compare the structural differences in the vicinity of the active sites of these enzymes and discuss which changes could give rise to the different behavior of H(peroxo) and Q. In particular, analysis of multiple crystal structures reveals that T213 in MMOH and the analogous T201 in TMOH, located in the immediate vicinity of the active site, have different rotatory configurations. We study the rotational energy profiles of these threonine residues with the use of molecular mechanics (MM) and quantum mechanics/molecular mechanics (QM/MM) computational methods and put forward a hypothesis according to which T213 and T201 play an important role in the formation of different types of peroxodiiron(III) species in MMOH and ToMOH. The hypothesis is indirectly supported by the QM/MM calculations of the peroxodiiron(III) models of ToMOH and the theoretically computed Mo?ssbauer spectra. It also helps explain the formation of two distinct peroxodiiron(III) species in the T201S mutant of ToMOH. Additionally, a role for the ToMOD regulatory protein, which is essential for intermediate formation and protein functioning in the ToMO system, is advanced. We find that the low quadrupole splitting parameter in the Mo?ssbauer spectrum observed for a ToMOH(peroxo) intermediate can be explained by protonation of the peroxo moiety, possibly stabilized by the T201 residue. Finally, similarities between the oxygen activation mechanisms of the monooxygenases and cytochrome P450 are discussed.     

Quantum mechanics/molecular mechanics minimum free-energy path for accurate reaction energetics in solution and enzymes: sequential sampling and optimization on the potential of mean force surface

To accurately determine the reaction path and its energetics for enzymatic and solution-phase reactions, we present a sequential sampling and optimization approach that greatly enhances the efficiency of the ab initio quantum mechanics/molecular mechanics minimum free-energy path (QM/MM-MFEP) method. In the QM/MM-MFEP method, the thermodynamics of a complex reaction system is described by the potential of mean force (PMF) surface of the quantum mechanical (QM) subsystem with a small number of degrees of freedom, somewhat like describing a reaction process in the gas phase. The main computational cost of the QM/MM-MFEP method comes from the statistical sampling of conformations of the molecular mechanical (MM) subsystem required for the calculation of the QM PMF and its gradient. In our new sequential sampling and optimization approach, we aim to reduce the amount of MM sampling while still retaining the accuracy of the results by first carrying out MM phase-space sampling and then optimizing the QM subsystem in the fixed-size ensemble of MM conformations. The resulting QM optimized structures are then used to obtain more accurate sampling of the MM subsystem. This process of sequential MM sampling and QM optimization is iterated until convergence. The use of a fixed-size, finite MM conformational ensemble enables the precise evaluation of the QM potential of mean force and its gradient within the ensemble, thus circumventing the challenges associated with statistical averaging and significantly speeding up the convergence of the optimization process. To further improve the accuracy of the QM/MM-MFEP method, the reaction path potential method developed by Lu and Yang [Z. Lu and W. Yang, J. Chem. Phys. 121, 89 (2004)] is employed to describe the QM/MM electrostatic interactions in an approximate yet accurate way with a computational cost that is comparable to classical MM simulations. The new method was successfully applied to two example reaction processes, the classical SN2 reaction of Cl-+CH3Cl in solution and the second proton transfer step of the reaction catalyzed by the enzyme 4-oxalocrotonate tautomerase. The activation free energies calculated with this new sequential sampling and optimization approach to the QM/MM-MFEP method agree well with results from other simulation approaches such as the umbrella sampling technique with direct QM/MM dynamics sampling, demonstrating the accuracy of the iterative QM/MM-MFEP method.     

Multireference ab initio quantum mechanics/molecular mechanics study on intermediates in the catalytic cycle of cytochrome P450(cam)

We have investigated the elusive reactive species of cytochrome P450(cam) (Compound I), the hydroxo complex formed during camphor hydroxylation, and the ferric hydroperoxo complex (Compound 0) by combined quantum mechanical/molecular mechanical (QM/MM) calculations, employing both density functional theory (DFT) and correlated ab initio methods. The first two intermediates appear multiconfigurational in character, especially in the doublet state and less so in the quartet state. DFT(B3LYP)/MM calculations reproduce the relative energies from correlated ab initio QM/MM treatments quite well, except for the splitting of the lowest A(1u)-A(2u) radical states. The inclusion of dynamic correlation is crucial for the proper ab initio treatment of these intermediates.     

Quantum mechanics/molecular mechanics electrostatic embedding with continuous and discrete functions

A quantum mechanics/molecular mechanics (QM/MM) implementation that uses the Gaussian electrostatic model (GEM) as the MM force field is presented. GEM relies on the reproduction of electronic density by using auxiliary basis sets to calculate each component of the intermolecular interaction. This hybrid method has been used, along with a conventional QM/MM (point charges) method, to determine the polarization on the QM subsystem by the MM environment in QM/MM calculations on 10 individual H(2)O dimers and a Mg(2+)-H(2)O dimer. We observe that GEM gives the correct polarization response in cases when the MM fragment has a small charge, while the point charges produce significant over-polarization of the QM subsystem and in several cases present an opposite sign for the polarization contribution. In the case when a large charge is located in the MM subsystem, for example, the Mg(2+) ion, the opposite is observed at small distances. However, this is overcome by the use of a damped Hermite charge, which provides the correct polarization response.     

Conical intersections in solution: formulation, algorithm, and implementation with combined quantum mechanics/molecular mechanics method

The significance of conical intersections in photophysics, photochemistry, and photodissociation of polyatomic molecules in gas phase has been demonstrated by numerous experimental and theoretical studies. Optimization of conical intersections of small- and medium-size molecules in gas phase has currently become a routine optimization process, as it has been implemented in many electronic structure packages. However, optimization of conical intersections of small- and medium-size molecules in solution or macromolecules remains inefficient, even poorly defined, due to large number of degrees of freedom and costly evaluations of gradient difference and nonadiabatic coupling vectors. In this work, based on the sequential quantum mechanics and molecular mechanics (QM/MM) and QM/MM-minimum free energy path methods, we have designed two conical intersection optimization methods for small- and medium-size molecules in solution or macromolecules. The first one is sequential QM conical intersection optimization and MM minimization for potential energy surfaces; the second one is sequential QM conical intersection optimization and MM sampling for potential of mean force surfaces, i.e., free energy surfaces. In such methods, the region where electronic structures change remarkably is placed into the QM subsystem, while the rest of the system is placed into the MM subsystem; thus, dimensionalities of gradient difference and nonadiabatic coupling vectors are decreased due to the relatively small QM subsystem. Furthermore, in comparison with the concurrent optimization scheme, sequential QM conical intersection optimization and MM minimization or sampling reduce the number of evaluations of gradient difference and nonadiabatic coupling vectors because these vectors need to be calculated only when the QM subsystem moves, independent of the MM minimization or sampling. Taken together, costly evaluations of gradient difference and nonadiabatic coupling vectors in solution or macromolecules can be reduced significantly. Test optimizations of conical intersections of cyclopropanone and acetaldehyde in aqueous solution have been carried out successfully.     

Quantum mechanics/molecular mechanics strategies for docking pose refinement: distinguishing between binders and decoys in cytochrome C peroxidase

We investigate the effect of systematically applying molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) to docked poses in an attempt to improve the correspondence between theoretical prediction and experimental observation. The proposed scheme involves running a short time scale MD simulation on a docked ligand pose (and any known structurally important crystal structure waters in the active site), followed by QM/MM minimization. Both of these steps are relatively fast for moderately sized ligands; longer time scale MD involving the protein is not found to improve the results. The final binding energy is given in terms of the QM/MM total energy, a van der Waals correction, and a term to account for desolvation effects. This methodology is first tested with a trypsin inhibitor, for which we establish the importance of running MD before reoptimizing with QM/MM. The method is then applied to cytochrome c peroxidase using a set of binders and decoys. In this example, the proposed methodology affords much better discrimination between binders and decoys than the traditional docking approach used. For both systems presented, application of this protocol results in a significantly better energetic ranking and a smaller root mean squared deviation from known crystallographic ligand poses. This work highlights the importance of including polarization effects through QM/MM and of sampling with MD to refine a set of initial docked poses.     

Fragmentation-based QM/MM simulations: length dependence of chain dynamics and hydrogen bonding of polyethylene oxide and polyethylene in aqueous solutions

Molecular fragmentation quantum mechanics (QM) calculations have been combined with molecular mechanics (MM) to construct the fragmentation QM/MM method for simulations of dilute solutions of macromolecules. We adopt the electrostatics embedding QM/MM model, where the low-cost generalized energy-based fragmentation calculations are employed for the QM part. Conformation energy calculations, geometry optimizations, and Born-Oppenheimer molecular dynamics simulations of poly(ethylene oxide), PEO(n) (n = 6-20), and polyethylene, PE(n) ( n = 9-30), in aqueous solution have been performed within the framework of both fragmentation and conventional QM/MM methods. The intermolecular hydrogen bonding and chain configurations obtained from the fragmentation QM/MM simulations are consistent with the conventional QM/MM method. The length dependence of chain conformations and dynamics of PEO and PE oligomers in aqueous solutions is also investigated through the fragmentation QM/MM molecular dynamics simulations.     

Calculation of wave‐functions with frozen orbitals in mixed quantum mechanics/molecular mechanics methods. II. Application of the local basis equation

The application of the local basis equation (Ferenczy and Adams, J. Chem. Phys. 2009 , 130, 134108) in mixed quantum mechanics/molecular mechanics (QM/MM) and quantum mechanics/quantum mechanics (QM/QM) methods is investigated. This equation is suitable to derive local basis nonorthogonal orbitals that minimize the energy of the system and it exhibits good convergence properties in a self‐consistent field solution. These features make the equation appropriate to be used in mixed QM/MM and QM/QM methods to optimize orbitals in the field of frozen localized orbitals connecting the subsystems. Calculations performed for several properties in divers systems show that the method is robust with various choices of the frozen orbitals and frontier atom properties. With appropriate basis set assignment, it gives results equivalent with those of a related approach [G. G. Ferenczy previous paper in this issue] using the Huzinaga equation. Thus, the local basis equation can be used in mixed QM/MM methods with small size quantum subsystems to calculate properties in good agreement with reference Hartree–Fock–Roothaan results. It is shown that bond charges are not necessary when the local basis equation is applied, although they are required for the self‐consistent field solution of the Huzinaga equation based method. Conversely, the deformation of the wave‐function near to the boundary is observed without bond charges and this has a significant effect on deprotonation energies but a less pronounced effect when the total charge of the system is conserved. The local basis equation can also be used to define a two layer quantum system with nonorthogonal localized orbitals surrounding the central delocalized quantum subsystem. © 2013 Wiley Periodicals, Inc.     

Protein/ligand binding free energies calculated with quantum mechanics/molecular mechanics

The calculation of binding affinities for flexible ligands has hitherto required the availability of reliable molecular mechanics parameters for the ligands, a restriction that can in principle be lifted by using a mixed quantum mechanics/molecular mechanics (QM/MM) representation in which the ligand is treated quantum mechanically. The feasibility of this approach is evaluated here, combining QM/MM with the Poisson-Boltzmann/surface area model of continuum solvation and testing the method on a set of 47 benzamidine derivatives binding to trypsin. The experimental range of the absolute binding energy (DeltaG = -3.9 to -7.6 kcal/mol) is reproduced well, with a root-mean-square (RMS) error of 1.2 kcal/mol. When QM/MM is applied without reoptimization to the very different ligands of FK506 binding protein the RMS error is only 0.7 kcal/mol. The results show that QM/MM is a promising new avenue for automated docking and scoring of flexible ligands. Suggestions are made for further improvements in accuracy.     

Converging ligand‐binding free energies obtained with free‐energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free‐energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa‐acid deep‐cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158–224 atoms). We use single‐step exponential averaging (ssEA) and the non‐Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi‐empirical PM6‐DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free‐energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Reliable treatment of electrostatics in combined QM/MM simulation of macromolecules

A robust approach for dealing with electrostatic interactions for spherical boundary conditions has been implemented in the QM/MM framework. The development was based on the generalized solvent boundary potential (GSBP) method proposed by Im et al. [J. Chem. Phys. 114, 2924 (2001)], and the specific implementation was applied to the self-consistent-charge density-functional tight-binding approach as the quantum mechanics (QM) level, although extension to other QM methods is straightforward. Compared to the popular stochastic boundary-condition scheme, the new protocol offers a balanced treatment between quantum mechanics/molecular mechanics (QM/MM) and MM/MM interactions; it also includes the effect of the bulk solvent and macromolecule atoms outside of the microscopic region at the Poisson-Boltzmann level. The new method was illustrated with application to the enzyme human carbonic anhydrase II and compared to stochastic boundary-condition simulations using different electrostatic treatments. The GSBP-based QM/MM simulations were most consistent with available experimental data, while conventional stochastic boundary simulations yielded various artifacts depending on different electrostatic models. The results highlight the importance of carefully treating electrostatics in QM/MM simulations of biomolecules and suggest that the commonly used truncation schemes should be avoided in QM/MM simulations, especially in simulations that involve extensive conformational samplings. The development of the GSBP-based QM/MM protocol has opened up the exciting possibility of studying chemical events in very complex biomolecular systems in a multiscale framework.     

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.     

Toward accurate barriers for enzymatic reactions: QM/MM case study on p-hydroxybenzoate hydroxylase

The hydroxylation reaction catalyzed by p-hydroxybenzoate hydroxylase has been investigated by quantum mechanical/molecular mechanical (QM/MM) calculations at different levels of QM theory. The solvated enzyme was modeled (approximately 23,000 atoms in total, 49 QM atoms). The geometries of reactant and transition state were optimized for ten representative pathways using semiempirical (AM1) and density functional (B3LYP) methods as QM components. Single-point calculations at B3LYP/MM optimized geometries were performed with local correlation methods [LMP2, LCCSD(T0)] and augmented triple-zeta basis sets. A careful validation of the latter approach with regard to all computational parameters indicates convergence of the QM contribution to the computed barriers to within approximately 1 kcal mol(-1). Comparison with the available experimental data supports this assessment.