Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Annu Rev Phys Chem. Author manuscript; available in PMC 2013 July 30.
Published in final edited form as:
PMCID: PMC3727228



Combined QM/MM methods provide an accurate and efficient energetic description of complex chemical and biological systems, leading to significant advances in the understanding of chemical reactions in solution and in enzymes. Progress in QM/MM methodology and application will be reviewed, with a focus on ab initio QM based approaches. Ab initio QM/MM methods capitalize on the accuracy and reliability of the associated quantum mechanical approaches, however at a much higher computational cost compared with semiempirical quantum mechanical approaches. Thus reaction path and activation free energy calculations based on ab initio QM/MM methods encounter unique challenges in simulation timescales and phase space sampling. Recent developments overcoming these challenges and enabling accurate free energy determination for reaction processes in solution and enzymes will be featured in this review, along with applications.

Keywords: enzyme catalysis, solution reaction, enzyme proficiency, conformational dynamics, potential of mean force


Understanding chemical reactions in solution and in enzymes is of ultimate importance. The majority of chemical reactions for synthesis or manufacture takes place in solution. Investigating solvent effects on reaction mechanisms has long been and continues to be an area of important and challenging inquiry. In biological systems, most biological functions are accomplished through the binding of ligands with proteins and/or a series of chemical reactions catalyzed by specific enzymes with different efficiency and specificity. The structure-function relationship and the catalytic role of enzymes is thus one of the most fundamental subjects in biochemical research. It is also essential for the development of new or better inhibitors and enzymes, which have important practical applications ranging from drug design to the development of novel catalysts in industrial processes. As an increasing amount of structural information for proteins and enzymes becomes available, the structure-function relationship becomes an even more important link in biological science. Quantitative tools such as simulations will make key contributions to the investigation of this topic.

In the process of chemical reactions, only a small number of atoms directly participate in bond forming or breaking processes. Many other atoms in the system do not undergo changes in electronic structure, but instead serve as a steric and electrostatic environment to influence the properties and reactivity of the active site. The QM/MM approach, first developed by Warshel and Levitt [1], is a multiscale/multiresolution simulation approach: the interesting part of the system, such as the active site of an enzyme, is described at the electronic level with quantum mechanics, while the rest of the system is described at the atomistic level with molecular mechanics. This combined quantum mechanical/molecular mechanical methodology allows reliable electronic structure calculations for enzymatic reactions with a realistic and atomistic description of enzyme environment. Such an approach takes advantage of the applicability and accuracy of the QM methods for chemical reactions in systems of several tens of atoms and of the computational efficiency of the MM description for the rest of the enzyme and solvent, which normally consists of many thousands of atoms. Development of combined QM/MM methods has enabled simulations of complex chemical and biological processes, leading to significant advances in our understanding. In particular, simulations have generated considerable insight into chemical reaction mechanisms in solution and in enzymes, as discussed in several recent reviews [2, 3, 4, 5, 6, 7, 8].

The QM/MM approach can be classified into two types according to the level of QM theory used. One type involves the use of semiempirical methods such as MNDO, AM1, PM3, empirical valence bond (EVB), and the recently developed self-consistent charge density functional tight binding (SCC-DFTB) method [9, 10, 11, 12]. The majority of the work in the field is carried out with these methods. The computational efficiency of semiempirical QM/MM is such that direct MD sampling is readily affordable, and thus free energy and reaction dynamics calculations can be routinely performed.

The other type of QM/MM calculation [13, 14, 15, 16, 17, 18, 19, 20, 7] uses ab initio QM via wavefunction theory or density functional theory. DFT is the most popular approach because of the optimal balance of efficiency and accuracy [21, 22, 23]. Recent development in approximate functionals has resulted in even better accuracy at a similar computational cost [24]. However, because ab initio QM calculations are computationally much more demanding than semiempirical methods, rigorous statistical mechanics sampling and reaction dynamics calculations with an ab initio QM/MM method are most challenging.

In the present review, we focus on the theoretical development and application of ab initio QM/MM free energy simulation methods for chemical reactions in solution and in enzymes.

Reaction free energy and rate

The fundamental properties of a chemical reaction process include the equilibrium distribution of states, kinetics, and the reaction mechanism. Thermodynamic laws tell us that equilibrium in a reaction process is determined by the free energy difference between the reactant and product states. This free energy difference determines the equilibrium, but not the reaction rate in general.

Reaction rate can be determined rigorously using theory based on classical statistical mechanics or quantum statistical mechanics. The exact classical rate constant k(T) at temperature T is [25, 3, 2]


where γ(T) is the transmission coefficient, β = 1/kBT, kB is the Boltzmann constant, and ΔG(T) is the molar activation free energy, i.e. the free energy difference between the reactant and TS. Within the transition state approximation, γ(T) = 1 , and ΔG(T) then determines the rate.

TS is a special state of the reacting system whose population controls how rapidly the system is able to transfer from the reactant valley to the product valley. If the dynamics of a system is projected onto an order parameter, or reaction coordinate, that characterizes the reaction process, TS is the point with the highest (free) energy on this one-dimension PMF. Identification of TS is crucial not only for the rate determination, but also for the understanding of the reaction mechanism.

Modeling enzyme catalysis

For enzymatic reactions, experimental studies provide crucial and indispensable information concerning the reaction mechanism, thermodynamics and kinetics. However, experimental data is often insufficient for determining the detailed reaction mechanism [26]. Most importantly, experimental study cannot directly determine the structure of the transition state, which is crucial for biomedicinal research such as inhibitor or drug design.

Complementary to experimental study, simulations can yield atomistic or even electronic information regarding the effects of site-specific interactions on the reaction process, the reaction path, and the structure of the transition state[14, 6, 2, 3, 12, 4, 27, 16]. Simulation studies of many enzymatic processes have addressed and also raised many important issues in enzyme catalysis such as covalent catalytic mechanisms [28, 29], contribution of the pre-organized electrostatic environment of enzymes[30, 31], the effects of strain and conformational dynamics of the enzyme-substrate complex, non-equilibrium dynamical effects, and quantum tunnel ling effects[32, 33].

Despite numerous simulation studies, the origin of enzymatic proficiency remains a current topic of interest because the contributions from different sources may vary in different enzymes [26, 34]. Accurate simulation methods will definitely help the understanding of those complex interactions in the enzymatic catalysis.

Unique role of solution reactions

Solvent molecules play significant roles in chemical reactions in solution [9, 35, 36, 33]. Energetically, the complex electrostatic interactions between the solvent and solute molecules create a reaction environment drastically different from gas phase, resulted in different chemical equilibrium and reaction rate. Dynamically, fluctuation and diffusion of the solvent molecules serve an energy bath to the motions of the reaction moieties, resulted in different relaxation and barrier-crossing dynamics. Last but not least, direct participation of the solvent molecule in the reaction not only creates a concentration effect, but also may alter the reaction path.

Solution reactions are important references for enzymatic reactions. In most cases, the slow reaction rate of the non-catalyzed solution reaction is the evolutionary driving force for enzymes; the exceptional proficiencies of enzymes are usually determined by the degree of difficulty of the solution reaction [37]. From this perspective, understanding enzymatic proficiency requires determining the differences between the enzyme catalyzed reaction and its corresponding solution reaction.

Solution reactions are often simulated with varying levels of approximation, particularly for simulations employing ab initio approaches. One common practice is to simplify the description of the solvent by using a continuum representation. This reduces the number of degrees of freedom in the simulation, but the isotropic continuum-medium model cannot always correctly reproduce the anisotropic, site-specific interactions between the solute and solvent molecules. To make reliable comparison, both the solution and enzymatic reactions must be simulated at the same level of accuracy.


The first requirement for simulating a reaction process is a potential energy surface that accurately describes the chemical changes with the redistribution of electrons. To account for this effect, in general quantum mechanics becomes necessary.

With currently available computing resources, a high-level ab initio QM description of the entire system is affordable only for small systems of less than 100 atoms. For reactions in solution or in enzymes, too many electronic degrees of freedom of the system make it difficult, if not impossible, to describe the whole system quantum mechanically. Although semiempirical QM methods using a linear scaling method such as the divide-and-conquer approach have been applied to large proteins and enzymes [38, 39, 40, 41, 42, 43], the computational speed is still a factor and the results may be less reliable because of the semiempirical approximations.

An effective and attractive solution to the energetic description of macromolecules is found in the combined QM/MM method [1, 11]. By partitioning the entire system into QM and MM subsystems, the total internal energy of the QM/MM system can be written as


The first two terms on the right hand side are the QM internal energy and the electrostatic energy between the QM and MM subsystems, for which the calculation is described in next section. The remaining three terms are the van der Waals energy between the QM and MM subsystems, the covalent interaction energy between the two subsystems, and the purely MM interaction energy of the MM subsystem, respectively, which are described with a given MM force field.

Due to the diversity in their forms and parameters, simulations performed with different MM force fields may show somewhat large discrepancies in the results [44]. Instead of MM force fields, an effective fragment potential method has been developed to derive a simple potential energy function for complex systems from ab initio QM calculations [45].

QM/MM electrostatic interactions

Computing the first two terms on the right hand side of Eq. (2) is the core of the QM/MM methods. A straightforward procedure is to compute the two terms together by including the electrostatic potential from the MM atoms in the QM calculation; that is


where Heff is the effective QM Hamiltonian including the MM electrostatic potential. In this charge embedding scheme, the QM internal energy and the QM/MM electrostatic energy are computed together in a self-consistent manner. In other words, the polarization of the QM subsystem due to the presence of the MM subsystem is captured at the same level of QM theory.

In the mechanical embedding approach represented by the class of ONIOM methods[46], the systems are hierarchically split into different layers that are described at different levels of theory. The energy function of the outer layer automatically includes the energy of the inner layer but it is described with a lower-accuracy method.

QM/MM boundary

When one or more covalent bonds connect the QM and MM subsystems, how to handle the boundary between the QM and MM subsystems becomes an important question. Several approaches have been developed [11][13, 47][48, 19, 49, 20, 50]. Traditionally, “link” hydrogen atoms are added to the MM-bonded QM atoms to saturate the valence orbital of the QM subsystem[11]. Using additional hydrogen atoms to cap the QM subsystem is easy to implement, but the resulting system is thermodynamically different from the original one because the total number of atoms is different.

In contrast, two other approaches, namely the pseudobond method [13, 47] and the frozen local orbital method [48, 19, 49, 20, 50], do not bear such problems. In the pseudobond method, the MM-bonded QM atoms are assigned a special basis set and an effective core potential that is designed to mimic the correct covalent bonding involving the boundary QM atoms. By making such atoms with a free valence of 1, there is no need for additional atoms to saturate the QM/MM covalent bonds. Several methods, following the pseudobond method in using effective core potentials, have subsequently been developed, including: the quantum capping potential method [51, 52], the effective group potential technique [53], effective Hamiltonians from a minimum principle [54], variational optimization of effective core potentials for molecular properties [55] and multicentered valence-electron effective potentials [56].

In the frozen local orbital method, a set of specially designed local orbitals are assigned to the boundary QM or MM atoms to maintain closure of the QM subsystem. In different implementation schemes, the magnitude of the neighboring MM charges, the positions of the MM point charges, and the positions of the frozen orbitals can vary.

In connection with this QM/MM boundary problem, another important issue is how to compute the electrostatic interactions between the QM subsystem and the nearby MM atomic charges, in particular those MM atoms in covalent contact with the QM atoms. Direct inclusion of those MM point charges in the QM effective Hamiltonian can cause charge penetration and off-balance polarization of the QM electrons. To address this problem, rescaling of the MM charges [13] or smearing the MM charges by Gaussian distributions have been developed [57].

Long-range QM/MM electrostatic interactions

The importance of the correct description of long-range QM/MM electrostatic interactions has not been explored in most simulations, as many QM/MM simulations have been performed with stochastic boundary conditions. The inclusion of long-range QM/MM electrostatic interaction is not trivial because the well-established Ewald summation method may not be directly applicable. In fact, several technical concerns must be addressed.

The first issue concerns the treatment of a periodic QM subsystem. In biomolecular simulations, usually the simulation box must be large enough so that the interaction between the QM subsystem, which is often charged, with its images is negligible. The second concern is how to consider the periodically distributed MM charges in the QM calculation if the charge-embedding scheme is used. Even though Ewald-type methods have been developed for semiempirical methods[58, 59], a similar development for ab initio QM/MM methods is still lacking.

No doubt the long-range electrostatic interaction is important for stabilizing the enzyme structures and may also be important for enzyme functions. However, the key question in both concerns might be whether or not the long-range polarization is significant for QM subsystem. It has been observed that even for a charge-transfer reaction, the QM polarization effects are small when the MM point charges are more than 9 ~ 14 Å from the QM atoms [60]. Beyond this distance, the MM point charges merely contribute by providing a static electrostatic potential. Therefore a simple technical scheme has been adopted in our simulations with periodic boundary condition [61, 62, 60] in which a cutoff of 9 ~ 14 Å was used for selecting MM charges for the QM SCF calculation. The QM calculation yielded polarized QM ESP charges that were in turn used to represent the QM subsystem during MD simulations with the long range electrostatic interactions treated by the Particle-mesh Ewald Method [63].

To avoid the technical difficulties associated with periodic QM/MM systems, Cui and co-workers have made important progress in approximating long-range QM/MM electrostatic interactions by implementing a stochastic generalized solvent boundary with the SCC-DFTB method [64, 65]. York and co-workers have developed a variational electrostatic projection method that employs a continuum solvent model for the long-range electrostatic interactions [66].


A reaction path provides a clear picture of a chemical reaction mechanism as well as quantitative energetic information regarding the reactant, transition state, and product, which are required for reaction rate calculations. A reaction path also serves as the structural basis for the detailed analysis of site-specific interactions involved in the reaction process and can be used in the design of new inhibitors.

For reactions of small molecules in the gas phase, reaction paths can be determined as MEP on the total PES. The TS on the PES of the entire system plays the most important role for the rate. The free energies of the reactant and the TS can be determined with the harmonic approximation, based on QM calculations of frequencies.

For reactions in solution and in enzymes, the TS on the PES of the entire system only represents one of many such states, because there are many degrees of freedom (e.g. for the solvent and the protein) that are not directly related to the chemical changes but modulate the shape of PES. In such cases, the progress of a reaction can also often be marked by the structural changes of a reduced set of atoms, e.g. the solute conformation for solution reactions or the conformation of the substrate plus certain active site residues for enzymatic reactions. Quite conveniently, the reduced set of coordinates can be the same as the QM subsystems. Thus, one can construct the MFEP, a reaction path defined on the PMF of a reduced set of coordinates [61, 67, 68, 69, 70]. The contributions from the rest of the system are ensemble-averaged quantities obtained during MD simulations.

Reaction coordinate

Usually defined as a set of geometric or energetic parameters, the reaction coordinate characterize a reaction as the one-dimensional profile connecting reactants to products. The choice of the reaction coordinate is critical. It has been shown that the improper choice of a reaction coordinate can bias the simulation and yield slower convergence [71]. The problem of an inappropriately chosen reaction coordinate is more severe in simulations using coordinate driving techniques in which the choice of the reaction coordinate not only strongly influences the efficiency of sampling the phase space, but also causes technical difficulties in generating continuous reaction path. In many cases, the determination of the reaction coordinate is non-trivial, especially for many complicated reactions catalyzed by enzymes as the changes in specific geometrical quantities (e.g. interatomic distances) might be stepwise or nonlinearly correlated. Furthermore, the environmental degrees of freedom may also contribute to the reaction coordinate as dynamic effects [25].

Reaction path optimization

The mathematical methods for determining reaction path are the same for the MEP on the total PES as for the MFEP on the PMF surface of the QM subsystem. The classical method for determining a reaction path is the intrinsic reaction coordinate method in which the steepest descent path is computed from the TS by following the gradient downhill. This method is not efficient for macromolecular systems and also requires knowledge about the structure of the transition state, which is a difficult task without prior knowledge of the reaction path.

Several algorithms have been developed for the determination of reaction paths for large systems in which a chain of conformations along the reaction coordinate is simultaneously optimized. Commonly used methods include the nudged elastic band (NEB) method[72], the Ayala-Schlegel second order minimum energy path method[73], and the string methods [74]. Our experience indicates that the NEB method is simple to implement but converges slowly and has difficulty locating transition states. The Ayala-Schlegel method can converge to the correct path but its application to large molecular systems is troublesome because it requires an initial guess that is enough close to the exact path. As an alternative, a quadratic string method (QSM) was developed in our laboratory which has been shown to yield better performance than the NEB method [75]. It has also been found that for large systems it is often important to select a small subset of coordinates as the chemical metric to define the path length for properly positioning the states along the reaction path [76].

One key component of MEP determination is the minimization of the QM/MM total energy, in which the QM calculation of energy and gradient is the bottleneck. To reduce the number of QM energy and gradient evaluations, Zhang, Liu, and Yang first developed an efficient iterative optimization for such a minimization in which the QM and MM subsystem optimizations are carried out sequentially, rather than concurrently [14]. The key is to use a simplified, often QM ESP charge based, energy function during the optimization of the MM subsystem. Friesner and co-workers introduced a correction term to improve the ESP charge approximation [20], and a similar method was proposed by Thiel and co-workers [77]. To improve the ESP charge approximation, Morokuma and co-workers have employed multiple QM calculations in the MM optimization process, as well as ESP multipoles [78].

Transition path sampling

Instead of classical rate theory that depends on the determination of the TS, the reaction rate can be determined by the probability ratio of finding successful dynamic reaction paths in multi-dimensional phase space, following the transition path sampling method [79]. The simulation of a large scale conformational transition in an enzyme has also been reported [80]. The major obstacle for the application of this method to reaction processes is the high computational cost for ab initio QM/MM calculations because a large number of paths must be sampled to ensure a converged ratio of successful/unsuccessful paths.

Reaction Path Potential based on ab initio QM/MM methods

The calculation of the QM internal energy and QM/MM electrostatic interactions in Eq. (3) is the bottleneck for QM/MM calculations. In the charge embedding scheme, rigorous treatment of this term requires a SCF calculation for each different QM or MM conformation, which is quite costly for ab initio QM methods. It would be desirable if one can develop an approximate, ab initio-based QM/MM energy function for phase space sampling and long timescale simulations without the need for SCF calculations.

Lu and Yang developed the ab initio RPP method [81]. The approach involves separating the QM energy into two components: a QM internal energy and an electrostatic interaction energy between the QM and MM subsystems. Each component is then expanded analytically in terms of both the fluctuations of the QM geometry and the MM electrostatic potential. The resulting RPP provides a simple, analytic expression for the QM/MM total energy that is valid in the vicinity of the reaction path.

RPP approximates the electrostatic term as the Coulombic interactions between the QM ESP charges (and multipoles) and the MM atomic charges, i.e.


Here, Qi(rQM, rMM) is the ESP fitted charge of QM atom i, and qj is the point charge of MM atom j from the MM force field. Whereas the MM atomic charges are constant in common force fields, the QM electrostatic charges are clearly dependent on both rQM and rMM. The QM internal energy, E1(rQM, rMM), is then defined as


This QM internal energy is the energy of the QM system in the presence of the electrostatic potential of the MM atoms, minus the Coulombic interactions between the QM ESP charges and MM atomic charges. The QM energy also depends on both rQM and rMM. Obviously, the term EQMMMESP(rQM,rMM) captures the essence of the electrostatic interactions between the QM and MM subsystems. The significance of this QM ESP charge expression is that it provides us great flexibility as the interactions are now expressed in a classical pairwise MM form. To improve the accuracy of ESP fitting, we can add electrostatic multipoles to each atom or bond.

The most important advantage of this separation scheme is that we can now introduce polarization effects into both terms as high-order perturbations. In the RPP model, the QM ESP charges are assumed to respond linearly to changes in both the external electrostatic potentials at the atomic sites and the geometries of the QM system, as


where QQM,i0 is a reference ESP charge for QM atom i. The QM reference charges are determined at a given QM geometry {rQM,j0} and with a given external MM electrostatic potential {vMM(rQM,j0)} at the position of the QM atoms. After the perturbation of the QM geometry and the external MM electrostatic potential, the polarized charges are determined through two response kernels, namely [81, 82],


and [81]


The second kernel can now be computed analytically with a recently introduced ESP fitting method[83]. Just as in the polarization effects introduced in the term EQMMMESP(rQM,rMM), we also expand the term E1(rQM, rMM) in a similar manner [81, 62].

The approximate QM/MM total energy becomes


It has been shown that this linear-polarization QM ESP charge model yields quite accurate energetics for reaction systems, given a good initial reference state [81, 61]. It is even possible to further simplify this QM ESP charge model by truncating the polarization effects of the QM ESP charges [61, 62]. Applications with these simplified QM ESP charge models in the ab initio QM/MM simulation of reaction processes have been successful.


Classical alchemical free energy simulation

Many methods have been developed for computing the free energies of important chemical processes on the basis of classical statistical simulations. Most popular methods include free energy perturbation, thermodynamic integration (TI), umbrella sampling, slow growth, and a recently developed fast growth method [84]. If one is only interested in equilibrium states, a classical free energy simulation method might be employed together with a so-called “alchemical” molecular transformation technique to provide the free energy difference between stationary states of a reaction process[85].

It has been pointed out that the direct application of TI and FEP with a QM/MM energy function requires special care for handling the end states [86]. In classical MM force fields, the energy functions are continuous and (almost always) infinitely differentiable with respect to the variables. As a result, without any difficulty one can create nonphysical states of the molecule, e.g. partial existence of an atom or a group, in the simulation process. Nonetheless QM methods in general do not easily allow a direct smooth process of creating or annihilating atoms. To overcome this difficulty, many methods have been developed including mixing of energy functions [9, 86], utilization of a reference state [87, 88, 89, 90], and dual-transformation [91].

Compared to their critical roles in the MM simulations, classical alchemical free energy simulations are less popular in the QM/MM simulations of the condensed phase reactions. Two exceptions are the calculations of residue pKa of protein molecules [92] and redox potentials [86, 89, 60].

PMF calculation along a reaction coordinate

When information regarding the reaction rate or the TS structure is sought, explicit modeling of the reaction process becomes necessary. Given an appropriate reaction coordinate, sampling-based free energy calculation methods can be applied to the study of reaction processes if the QM calculations are fast enough, as in semiempirical QM methods. In such cases, the system is driven from one state to another along the reaction coordinate and the PMF is computed during this driving process.

Umbrella sampling method computes the PMF by improving the sampling of high-energy, low-probability regions of phase space. Usually the results of umbrella sampling are processed with the weighted histogram analysis method, which provides reduced statistical error. Recently, Thiel and co-workers have developed an “umbrella integration” method which combines the advantages of umbrella sampling and thermodynamic integration [93].

Parrinello and co-workers [94] have developed a meta-dynamics method in which the energy surface of the target system is gradually modulated to allow enhanced sampling of rare conformational states.

Slow growth is in fact a variant of TI. Instead of sampling a few fixed points on the reaction/transformation path and computing a converged free energy gradient for each point, slow growth slowly drives the system from one state to the other using a very small step size. The work exerted during this process is summed up to yield an upper-bound estimate of the free energy difference between the two states. The reverse process yields a lower-bound estimate. Together, the results of multiple forward and backward processes yield a good estimate of the free energy difference [95]. When the slow growth method is combined with the coordinate driving technique, the free energy change with respect to the reaction coordinate may be obtained.

In analogy to the slow growth method, a “fast growth” method [84] has been developed and has been used to simulate chemical reactions [96, 97]. Unlike the slow growth method in which each simulated work is an estimate of the true free energy, in fast growth the ensemble averaged work is the upper-bound (or lower-bound for the reverse process) estimate of the true free energy. This estimate slowly converges to the true free energy upon the convergence of the ensemble. Since each individual simulation is often only in pico-second timescale range, ab initio QM/MM simulations at this timescale might be possible. However, the main drawback of the fast growth method is that it requires many simulations to yield a good estimate of the free energy, which always possess a large amount of uncertainty [98].

It is important to realize that the phase space sampling methods based on PMF calculation along a reaction coordinate usually require converged sampling, and of course prior knowledge about the reaction coordinate. Typically, this method is affordable for semiempirical QM/MM methods. As for ab initio QM/MM, few simulations have been reported[99, 62]. With the steady improvement of computer speed, direct sampling methods may play more important roles in the reaction simulations.

Aimed to reduce the costs of ab initio QM calculations, two methods have been developed utilizing the information about the reaction coordinate and path. Jorgensen and co-workers have developed a QM-FE method, in which the reaction path optimized for gas-phase reaction process is used to carry out free energy simulation in condensed phase [35]. Warshel and co-workers have developed a QM(ai)/MM method, the key of which is that the sampling and free energies are first computed with a simplified EVB potential, and then corrected to ab initio QM level with FEP combined with the linear response approximation [87, 90].


Realizing the importance of the reaction path, Yang and co-workers have devoted much effort into developing reaction path optimization and free energy calculation techniques based on ab initio QM/MM methods [14, 61, 62]. The first method, termed QM/MM-FE (or QM/MMFEP) [14], is carried out in two stages: optimization of the reaction path and calculation of the free energies for the path.

The reaction path determination has been designed as a sequential optimization process. To exemplify the advantages of this sequential approach, let us first examine the situation where the QM and MM degrees of freedom are optimized concurrently. In this case, each change of the QM and/or MM geometries will require a new QM evaluation of the energy and gradient. While the energy and gradient calculation is computationally inexpensive for semiempirical QM methods, it is expensive for ab initio QM methods. As it usually takes hundreds to thousands of minimization steps to obtain an adequately converged structure for the whole QM/MM system, the cost for ab initio QM calculations would be prohibitive. Inspired by the fact that in certain types of optimization problems it may be desirable to break up a process into an iterative sequence of two or more steps, we perform the optimization as follows: (i) optimize a subset of the system, A, with the rest of the system, B, held constant; (ii) optimize B with A fixed according to the results obtained in step (i). Thus in the QM/MM-FE method, instead of a concurrent optimization of the QM and MM degrees of freedom, an iterative, sequential optimization protocol was developed that has proven to be effective in reducing the number of QM energy and gradient evaluations. The main idea is that, starting from a given structure of the QM/MM molecular system, the optimization is separated into two processes. One first optimizes rQM with fixed rMM, which is at an approximate minimum in the MM degrees of freedom. Afterwards, the conformation of the MM subsystem, rMM, is optimized with fixed rQM. To reduce the number of QM evaluations, an approximate QM/MM total energy, E~(rQM,rMM), is used for the MM optimization. In this approximate QM/MM total energy, the electrostatic interactions between the QM and MM atoms are approximated by the Coulombic interactions between the point charges of the MM atoms and the ESP fitted charges of the QM atoms. This process is then iterated until convergence, which is normally achieved within a few iterations (often less than 10). Expressed as an algorithm, the ab initio QM/MM-FE optimization procedure is as follows:

  1. Initiate a structure of the QM subsystem rQM(0), and set the cycle number n = 0;
  2. Increase the cycle number n = n + 1;
    1. Carry out an MM minimization with QM atoms fixed at rQM(n1):
    2. Carry out a QM optimization with MM atoms fixed at rMM(n):
  3. Go to Step (2) until converged

This process can be carried out individually for a single point on the reaction path, e.g. the reactant and product states, or simultaneously for a chain of conformations along the reaction coordinate with a chain-of-states optimization algorithm such as the NEB method, the Ayala-Schlegel method, or the superlinearly convergent QSM.

Once the reaction path is determined, a FEP simulation with the approximate QM/MM energy function E~(rQM,rMM) can be carried out for the optimized QM conformations on the reaction path, similar to the QM-FE method [35]. The validity of QM/MM-FE has been confirmed by other laboratories [100, 101, 102].

QM/MM-MFEP: Path optimization on a QM PMF surface

The development of the QM/MM-FE method has provided a viable way for computing accurate free energies of reactions in enzymes. But one limitation has hampered the application of this method to the simulation of solution reactions. In the QM/MM-FE method, the reaction path is optimized on the QM/MM PES, starting from a given initial structure. As a result, the optimized path is influenced by the choice of the initial conformation. In many enzyme-substrate complexes, the dependence of the initial conformation may not cause problems because the active site of the enzyme is usually protected from bulk solvent. When the reaction occurs in solution, or the enzyme active site is exposed to solvent, this dependence will be significant. The rapid exchange of solvent molecules can also cause difficulty for the convergence of the path optimization process. Similar observations have been made for the enzymatic reactions in which the enzyme undergoes significant conformational changes during the reaction process.

To overcome this problem, the QM/MM-MFEP method has been developed in which the reaction path is optimized on the PMF surface of the QM degrees of freedom, instead of the total energy surface [61]. Within the QM/MM context, the thermodynamics of the entire system is simplified by defining the PMF of a QM/MM system in terms of the QM conformation as


where E(rQM, rMM) is the total energy of the entire system expressed as a function of the coordinates of the QM and MM subsystems, rQM and rMM, respectively. The gradient of the PMF, also known as the free energy gradient, is then


which appears conveniently as the ensemble average of the gradient of the QM atoms, obtained from MD simulations of the MM atoms. In practice, we used the free energy perturbation method and its associated gradient, instead of Eqs. (12-13) [61, 62].

Our construction of the PMF and PMF gradient is different from other work in terms of the variables of the PMF. For reasons we have discussed in the section of reaction coordinates, we allow all of QM degrees of freedom to contribute to the reaction coordinate whereas others have often used one or a few predefined geometric terms. In the latter case, a Jacobian term may be required to correctly include the effects of geometrical constraints on the reaction coordinates.

Because the QM degrees of freedom are coupled with the MM degrees of freedom, a straightforward minimization algorithm requires each step in the optimization of the QM conformations to be associated with converged sampling of the MM ensemble. In this optimization scheme, every QM optimization step to a new conformation on the PMF surface is followed by a course of MD sampling of the MM conformations, usually with a simulation time of 100 ~ 1000 ps. Such extensive MM sampling is required so that the QM PMF and gradient obtained are sufficiently accurate for successful structure and reaction path optimization of the QM subsystem. Thus, this method is less efficient in practical simulations. First, 100 ~ 1000 ps of MD simulations on a system of ~ 10,000 atoms is expensive. Second, such MD simulations must be repeated for each step in the QM optimization process. The intrinsic fluctuations and limited simulation times for the MD sampling may also contribute to slow convergence in the path optimization if a new MD simulation is always begun immediately after every QM geometry optimization step.

Sequential sampling and optimization for QM/MM-MFEP method

As in the QM/MM-FE method, to improve the efficiency of the QM/MM-MFEP method, we can reformulate the concurrent optimization of the QM subsystem and the statistical sampling of the MM subsystem into an iterative step of sequential MD sampling of the MM system at a fixed QM structure and subsequent optimization of the QM subsystem within the fixed MM conformational ensemble [62]. The approximate energy function E~(rQM,rMM), introduced in the QM/MM-FE method and well described in the RPP model, plays a crucial role in reducing the number of QM energy and gradient evaluations required. In the MD sampling of the MM conformations, E~(rQM,rMM) acts as a reference sampling energy function to drive the motion of the MM atoms without performing a QM calculation at every MD step. In the optimization of the QM subsystem within the fixed-size MM conformational ensemble, it can be used to avoid QM calculations associated with each new conformation in the MM ensemble.

The algorithm of the QM/MM-MFEP method can be described as follows.

  1. Initiate a structure of the QM subsystem, rQM(0), and set cycle number n = 0.
  2. Increase cycle number n = n + 1;
    1. Carry out MD sampling of the MM ensemble with QM atoms fixed at rQM(n1):
      {rMM(n)(τ),τ=1,,N}MD sampling based onEref(rMM)
      where τ is the step of the MD simulation, N is the number of MD steps, the reference QM structure is derived from rQM(n1), the QM geometry from the previous iteration.
    2. Carry out a QM optimization with the MM ensemble fixed at {rMM(n)(τ)} where the object of minimization is the QM PMF (or QM free energy) in the n-th iteration given by a finite sum representation of free energy perturbation as
      and the corresponding gradient with respect to the i-th QM coordinate is also given by the finite sum
      which accounts for the fact that the samples were obtained from a fixed MD simulation of a reference state. Eref and Aref are the reference energy function and reference free energy in the free energy perturbation expression.
    3. Update the reference structure based on the minimized QM structure rQM(n).
  3. Go to Step (2) until converged.

The key feature of our new QM/MM-MFEP algorithm is the iterative QM optimization in a fixed MM ensemble, Eq. (14) above. Because the MM ensemble, {rMM(n)(τ),τ=1,,N}, is finite and remains fixed throughout the course of the QM optimization for rQM(n), one can obtain the precise PMF, Eq. (15), and its gradient, Eq. (16), defined within this ensemble. This circumvents the difficult and costly convergence problems associated with MM sampling. The optimization of the PMF can be carried out efficiently using classical numerical optimization tools. Each optimized QM structure rQM(n) in turn provides the next reference QM structure and its energy function, rQMref(n), for the next round of MD sampling of the MM conformations. Each optimized QM structure should improve on the previous one by providing a better QM geometry and corresponding ESP charges for the MM simulation in the next cycle.

The use of a finite, fixed-size ensemble of MM conformations improves the utilization of the MM conformations and avoids repetitive MD sampling at each step of the QM structure optimization. Thus, instead of performing excessive MD samplings that likely have significant overlap with each other, a few cycles of MD simulation are sufficient to yield converged results in the current method. Applications have shown that our QM/MM-MFEP method converges as efficiently as the QM/MM-FE method.


The classical dynamic effects beyond the classical TS theory is described in the transmission coefficient γ(T), the determination of which requires an MD simulation starting from the TS [25]. Main quantum mechanical effects can be included in the quantum PMF calculations. There are only several ab initio QM/MM studies reported so far. Cui calculated γ(T) for the proton transfer reaction catalyzed in TIM using variational transition state theory [103]. Wang et al. used the reaction path potential, which enables rapid evaluation of the potential energy around the reaction path, to calculate γ(T) [104] and also the quantum PMF with a centroid path integral approach [105]. When the reaction path and its TS captures the correct physics, TS theory is a good approximation, and γ(T) is nearly 1 [3, 103, 104, 105]. Wang’s calculations on TIM resulted in kinetic isotope effects in excellent agreement with experiments and revealed that the main QM effect is in the zero-point energy.


Many enzymatic reaction processes have been modeled by QM/MM methods, of which most were investigated with semiempirical QM/MM methods. As only a small number of cases have been studied using ab initio QM/MM calculations, the number of applications of ab initio QM/MM free energy simulations are even fewer because of the difficulty of combining ab initio QM with statistical sampling. A comprehensive list can be found in Ref [8]. In this section, we briefly review several enzyme systems whose mechanisms have been scrutinized by ab initio QM/MM, some in combination with free energy simulation.

The enzyme orotidine 5′-monophosphate decarboxylase (ODCase) is one of the most proficient enzymes known and catalyzes the decarboxylation of orotidine 5′-monophosphate without any cofactor or metal ions. The solution reaction has been modeled with a continuum solvent model [106] which supported a carbene intermediate state. Semiempirical simulations have supported the direct decarboxylation mechanism [107, 108]. However the reaction barriers computed in later DFT simulations, including a combined CPMD and fast growth simulation, were still higher than the experimental measurement [109, 110, 97].

The enzyme 4-oxalocrotonate tautomerase (4OT) catalyzes the conversion between 2-oxo-4-hexenedioate and 2-oxo-3-hexenedioate. Cisneros et al. have carried out extensive QM/MM-FE simulations on the reaction mechanism and identified the absence of a catalytic acid and the role of an ordered water and the protein backbone in the catalysis. [16, 111, 112, 113] Thiel and co-works have also studied the mechanism of the enzyme and its mutant [114].

One enzyme that has undergone extensive characterization is triosephosphate isomerase (TIM), which catalyzes the reversible isomerization of dihydroxyacetone phosphate (DHAP) to D-glyceraldehyde-3-phosphate. The mechanism has been examined by semiempirical QM/MM methods [115, 116]. Later, several DFT-based QM/MM studies were been carried out [14, 117, 17, 118]. With a higher resolution structure of the enzyme, Friesner’s group modeled the reaction process and obtained results in good agreement with experimental data [119]. The correlation between the reaction and the loop motion has been explored in an interesting way in which the QM/MM simulation is employed in combination with a loop structure prediction algorithm.

The reaction of catechol O-methyltransferase has been simulated by Kollman and co-workers with ab initio QM/MM and QM-FE methods [120, 121, 122]. Rod and co-workers also simulated the reaction mechanism with a recently developed quantum mechanical thermodynamic cycle perturbation method [101].

Enolase has been simulated by both semiempirical and ab initio QM/MM methods [123, 15]. Chorismate mutase catalyzes the Claisen rearrangement from chorismate to prephenate. In addition to many semiempirical QM/MM simulations, recently it was modeled by an ab initio QM/MM method with optimization of reaction path [124, 125, 126]. For its medicinal significance, β-lactamase has been explored by many groups. Different reaction steps have been examined in detail [127, 128, 129]. Acetylcholinesterase has been simulated by QM/MM-FE with an emphasis on the effects of conformational dynamics on the reaction process [130, 131]. cAMP-dependent protein kinase [132, 133], histone lysine methyltransferase SET7/9 [134, 135, 136], and bacterial peptide deformylase [137, 138] have been studied by the QM/MM-FE method. The reaction mechanism of methane monooxygenase has been explored in detail by Friesner and co-workers[139, 140], as well as the P450cam pathway [141].

A large class of protein and enzymes have been investigated with the CPMD/MM method, including ion channels [142], aspartic protease [143], HIV-1 protease [144], and caspases [145].


An important issue in modeling enzyme catalysis is the role of conformational dynamics. Special interest has been raised recently for the contribution of conformational dynamics to enzyme catalysis [27]. For the study on this topic, a great deal of attention must be focused on two considerations: the timescale of the conformational dynamics and the coupling scheme between the enzyme dynamics and the reaction process. As agreement has been reached that the transmission coefficient in the transition state rate equation is affected by the stochastic dynamics of the enzyme but its variation has a small effect on the reaction rate in general [3], the challenge is to understand if/how the collective conformational dynamics of the enzyme has an impact on the reaction rate [27].

In the current popular form of the QM/MM approach, MM force fields are not polarizable. As efforts are underway to develop polarizable MM force fields, schemes have been proposed to combine QM with polarizable MM force fields [146, 147, 148, 149, 150], but more testings might be needed before a specific scheme is employed. In two scenarios this issue will become more complicated. One is the mutual polarization at the QM/MM boundary; the other is the charge transfer effect between the QM and MM atoms. Generally, the implementation and testing of polarizable MM force fields would be one immediate important topic in the QM/MM methodology development.


With the recent development of ab initio QM/MM free energy simulation methods, realistic and accurate simulations of reactions in solution and in enzymes become feasible. Application of such methods will provide comprehensive understanding of reactions in solution and in enzymes.


We are grateful to the National Institute of Health and the National Science Foundation for support.


quantum mechanics
molecular mechanics
free energy perturbation
minimum free energy path
potential of mean force
transition state
molecular dynamics
density functional theory
reaction path potential
electrostatic potential
minimum energy path
self-consistent field
potential energy surface


[1] Warshel A, Levitt M. Theoretical studies of enzymic reactions - dielectric, electrostatic and steric stabilization of carbonium-ion in reaction of lysozyme. Journal of Molecular Biology. 1976;103(2):227–249. [PubMed]
[2] Gao JL, Truhlar DG. Quantum mechanical methods for enzyme kinetics. Annual Review of Physical Chemistry. 2002;53:467–505. [PubMed]
[3] Garcia-Viloca M, Gao J, Karplus M, Truhlar DG. How enzymes work: Analysis by modern rate theory and computer simulations. Science. 2004 Jan;303(5655):186–195. [PubMed]
[4] Friesner RA, Guallar V. Ab initio quantum chemical and mixed quantum mechanics/molecular mechanics (qm/mm) methods for studying enzymatic catalysis. Annual Review of Physical Chemistry. 2005;56:389–427. [PubMed]
[5] Mulholland AJ. Modelling enzyme reaction mechanisms, specificity and catalysis. Drug Discovery Today. 2005 Oct;10(20):1393–1402. [PubMed]
[6] Warshel A, Sharma PK, Kato M, Xiang Y, Liu HB, Olsson MHM. Electrostatic basis for enzyme catalysis. Chemical Reviews. 2006 Aug;106(8):3210–3235. [PubMed]
[7] Zhang YK. Pseudobond ab initio qm/mm approach and its applications to enzyme reactions. Theoretical Chemistry Accounts. 2006 Aug;116(1-3):43–50.
[8] Senn HM, Thiel W. Atomistic Approaches in Modern Biology: from Quantum Chemistry to Molecular Simulations, volume 268 of Topics in Current Chemistry. 2007. Qm/mm methods for biological systems; pp. 173–290. Times Cited: 1.
[9] Hwang JK, King G, Creighton S, Warshel A. Simulation of free-energy relationships and dynamics of sn2 reactions in aqueous-solution. Journal of the American Chemical Society. 1988 Aug;110(16):5297–5311.
[10] Gao JL, Xia XF. A priori evaluation of aqueous polarization effects through monte-carlo qm-mm simulations. Science. 1992 Oct;258(5082):631–635. [PubMed]
[11] Field MJ, Bash PA, Karplus M. A combined quantum-mechanical and molecular mechanical potential for molecular-dynamics simulations. Journal of Computational Chemistry. 1990 Jul;11(6):700–733.
[12] Riccardi D, Schaefer P, Yang Y, Yu HB, Ghosh N, Prat-Resina X, Konig P, Li GH, Xu DG, Guo H, Elstner M, Cui Q. Development of effective quantum mechanical/molecular mechanical (qm/mm) methods for complex biological processes. Journal of Physical Chemistry B. 2006 Apr;110(13):6458–6469. [PubMed]
[13] Zhang YK, Lee TS, Yang WT. A pseudobond approach to combining quantum mechanical and molecular mechanical methods. Journal of Chemical Physics. 1999 Jan;110(1):46–54.
[14] Zhang YK, Liu HY, Yang WT. Free energy calculation on enzyme reactions with an efficient iterative procedure to determine minimum energy paths on a combined ab initio qm/mm potential energy surface. Journal of Chemical Physics. 2000 Feb;112(8):3483–3492.
[15] Liu HY, Zhang YK, Yang WT. How is the active site of enolase organized to catalyze two different reaction steps? Journal of the American Chemical Society. 2000 Jul;122(28):6560–6570.
[16] Cisneros GA, Liu HY, Zhang YK, Yang WT. Ab initio qm/mm study shows there is no general acid in the reaction catalyzed by 4-oxalocrotonate tautornerase. Journal of the American Chemical Society. 2003 Aug;125(34):10384–10393. [PubMed]
[17] Cui Q, Karplus M. Triosephosphate isomerase: A theoretical comparison of alternative pathways. Journal of the American Chemical Society. 2001 Mar;123(10):2284–2290. [PubMed]
[18] Cui Q, Elstner M, Karplus M. A theoretical analysis of the proton and hydride transfer in liver alcohol dehydrogenase (ladh) Journal of Physical Chemistry B. 2002 Mar;106(10):2721–2740.
[19] Philipp DM, Friesner RA. Mixed ab initio qm/mm modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide. Journal of Computational Chemistry. 1999 Nov;20(14):1468–1494.
[20] Murphy RB, Philipp DM, Friesner RA. A mixed quantum mechanics/molecular mechanics (qm/mm) method for large-scale modeling of chemistry in protein environments. Journal of Computational Chemistry. 2000 Dec;21(16):1442–1457.
[21] Parr RG, Yang WT. Density-Functional Theory of Atoms and Molecules. Oxford University Press; New York: 1989.
[22] Becke AD. Density-functional thermochemistry. iii. the role of exact exchange. J. Chem. Phys. 1993;98:5648–5652.
[23] Lee C, Yang WT, Parr RG. Development of the colle-salvetti correlationenergy formula into a functional of the electron density. Phys. Rev. B. 1988;37:785–789. [PubMed]
[24] Mori-Sanchez P, Cohen AJ, Yang WT. Self-interaction-free exchange-correlation functional for thermochemistry and kinetics. Journal of Chemical Physics. 2006 Mar;124(9):091102. [PubMed]
[25] Truhlar DG, Garrett BC, Klippenstein SJ. Current status of transition-state theory. Journal of Physical Chemistry. 1996 Aug;100(31):12771–12800.
[26] Kraut DA, Carroll KS, Herschlag D. Challenges in enzyme mechanism and energetics. Annual Review of Biochemistry. 2003;72:517–571. [PubMed]
[27] Hammes-Schiffer S, Benkovic SJ. Relating protein motion to catalysis. Annual Review of Biochemistry. 2006;75:519–541. [PubMed]
[28] Zhang XY, Houk KN. Why enzymes are proficient catalysts: Beyond the pauling paradigm. Accounts of Chemical Research. 2005 May;38(5):379–385. [PubMed]
[29] Bruice TC, Bruice PY. Covalent intermediates and enzyme proficiency. Journal of the American Chemical Society. 2005 Sep;127(36):12478–12479. [PubMed]
[30] Warshel A. Electrostatic origin of the catalytic power of enzymes and the role of preorganized active sites. Journal of Biological Chemistry. 1998 Oct;273(42):27035–27038. [PubMed]
[31] Cannon WR, Benkovic SJ. Solvation, reorganization energy, and biological catalysis. Journal of Biological Chemistry. 1998 Oct;273(41):26257–26260. [PubMed]
[32] Nam K, Prat-Resina X, Garcia-Viloca M, Devi-Kesavan LS, Gao JL. Dynamics of an enzymatic substitution reaction in haloalkane dehalogenase. Journal of the American Chemical Society. 2004 Feb;126(5):1369–1376. [PubMed]
[33] Olsson MHM, Warshel A. Solute solvent dynamics and energetics in enzyme catalysis: The s(n)2 reaction of dehalogenase as a general benchmark. Journal of the American Chemical Society. 2004 Nov;126(46):15167–15179. [PubMed]
[34] Kraut DA, Sigala PA, Pybus B, Liu CW, Ringe D, Petsko GA, Herschlag D. Testing electrostatic complementarity in enzyme catalysis: Hydrogen bonding in the ketosteroid isomerase oxyanion hole. Plos Biology. 2006 Apr;4(4):501–519. [PMC free article] [PubMed]
[35] Jorgensen WL. Free-energy calculations - a breakthrough for modeling organic-chemistry in solution. Accounts of Chemical Research. 1989 May;22(5):184–189.
[36] Gao JL, Garcia-Viloca M, Poulsen TD, Mo YR. Solvent effects, reaction coordinates, and reorganization energies on nucleophilic substitution reactions in aqueous solution. Advances in Physical Organic Chemistry. 2003;Vol 38, volume 38 of Advances in Physical Organic Chemistry:161–181. Times Cited: 6.
[37] Wolfenden R, Snider MJ. The depth of chemical time and the power of enzymes as catalysts. Accounts of Chemical Research. 2001 Dec;34(12):938–945. [PubMed]
[38] Yang WT. Direct calculation of electron-density in density-functional theory. Physical Review Letters. 1991 Mar;66(11):1438–1441. [PubMed]
[39] Yang WT, Lee TS. A density-matrix divide-and-conquer approach for electronic-structure calculations of large molecules. Journal of Chemical Physics. 1995 Oct;103(13):5674–5678.
[40] Yang WT, Pérez-Jordá J. Linear scaling methods for electronic structure calculations. In: Schleyer P.v.R., editor. Encyclopedia of Computational Chemistry. John Wiley & Sons; New York: 1998. pp. 1496–1513.
[41] Lee TS, York DM, Yang WT. Linear-scaling semiempirical quantum calculations for macromolecules. Journal of Chemical Physics. 1996 Aug;105(7):2744–2750.
[42] Lewis JP, Carter CW, Hermans J, Pan W, Lee TS, Yang WT. Active species for the ground-state complex of cytidine deaminase: A linear-scaling quantum mechanical investigation. Journal of the American Chemical Society. 1998 Jun;120(22):5407–5410.
[43] Liu HY, Elstner M, Kaxiras E, Frauenheim T, Hermans J, Yang WT. Quantum mechanics simulation of protein dynamics on long timescale. Proteins-Structure Function and Genetics. 2001 Sep;44(4):484–489. [PubMed]
[44] Hu H, Elstner M, Hermans J. Comparison of a qm/mm force field and molecular mechanics force fields in simulations of alanine and glycine dipeptides (ace-ala-nme and ace-gly-nme) in water in relation to the problem of modeling the unfolded peptide backbone in solution. Proteins: Structure, Function, and Genetics. 2003;50(3):451–463. [PubMed]
[45] Jensen JH, Day PN, Gordon MS, Basch H, Cohen D, Garmer DR, Kraus M, Stevens WJ. Modeling the Hydrogen Bond. American Chemical Society; 1994. An effective fragment method for modeling intermolecular hydrogen bonding-effects on quantum mechanical calculations.
[46] Svensson M, Humbel S, Froese RDJ, Matsubara T, Sieber S, Morokuma K. Oniom: A multilayered integrated mo+mm method for geometry optimizations and single point energy predictions. a test for diels-alder reactions and pt(p(t-bu)(3))(2)+h-2 oxidative addition. Journal of Physical Chemistry. 1996 Dec;100(50):19357–19363.
[47] Zhang YK. Improved pseudobonds for combined ab initio quantum mechanical/molecular mechanical methods. Journal of Chemical Physics. 2005 Jan;122(2):024114. [PubMed]
[48] Théry V, Rinaldi D, Rivail J-L, Maigret B, Ferenczy GG. Quantum mechanical computations on very large molecular systems: the local self-consistent field method. J. Comput. Chem. 1994;15:269–282.
[49] Murphy RB, Philipp DM, Friesner RA. Frozen orbital qm/mm methods for density functional theory. Chemical Physics Letters. 2000 Apr;321(1-2):113–120.
[50] Pu JZ, Gao JL, Truhlar DG. Generalized hybrid-orbital method for combining density functional theory with molecular mechanicals. Chemphyschem. 2005 Sep;6(9):1853–1865. [PMC free article] [PubMed]
[51] DiLabio GA, Hurley MM, Christiansen PA. Simple one-electron quantum capping potentials for use in hybrid qm/mm studies of biological molecules. J. Chem. Phys. 2002;116:9578–9584.
[52] Dilabio GA, Wolkow RA, Johnson ER. Efficient silicon surface and cluster modeling using quantum capping potentials. J. Chem. Phys. 2005;122:044708. [PubMed]
[53] Poteau R, Ortega I, Alary F, Solis AR, Barthelat J-C, Daudey J-P. Effective group potentials. 1. method. J. Phys. Chem. A. 2001;105:198–205.
[54] Yasuda K, Yamaki D. Simple minimum principle to derive a quantummechanical/molecular-mechanical method. J. Chem. Phys. 2004;121:3964–3972. [PubMed]
[55] von Lilienfeld OA, Tavernelli I, Rothlisberger U, Sebastiani D. Variational optimization of effective atom centered potentials for molecular properties. J. Chem. Phys. 2005;122:014113. [PubMed]
[56] Slavicek P, Martinez TJ. Multicentered valence electron effective potentials: a solution to the link atom problem for ground and excited electronic states. J. Chem. Phys. 2006;124:084107. [PubMed]
[57] Das D, Eurenius KP, Billings EM, Sherwood P, Chatfield DC, Hodoscek M, Brooks BR. Optimization of quantum mechanical molecular mechanical partitioning schemes: Gaussian delocalization of molecular mechanical charges and the double link atom method. Journal of Chemical Physics. 2002 Dec;117(23):10534–10547.
[58] Gao JL, Alhambra C. A hybrid semiempirical quantum mechanical and lattice-sum method for electrostatic interactions in fluid simulations. Journal of Chemical Physics. 1997 Jul;107(4):1212–1217.
[59] Nam K, Gao JL, York DM. An efficient linear-scaling ewald method for long-range electrostatic interactions in combined qm/mm calculations. Journal of Chemical Theory and Computation. 2005 Jan-Feb;1(1):2–13. [PubMed]
[60] Zeng XC, Hu H, Hu XQ, Cohen AJ, Yang WT. Qm/mm calculation of electron transfer process: Fractional electrons approach. J. Chem. Phys. 2007:xxxx.
[61] Hu H, Lu ZY, Yang WT. Qm/mm minimum free-energy path: Methodology and application to triosephosphate isomerase. Journal of Chemical Theory and Computation. 2007 Mar-Apr;3(2):390–406. [PMC free article] [PubMed]
[62] Hu H, Lu ZY, Parks JM, Burger SK, Yang WT. Qm/mm minimum free energy path for accurate reaction energetics in solution and enzymes: Sequential sampling and optimization on the potential of mean force surface. J. Chem. Phys. 2007:xxxx. [PubMed]
[63] Darden T, York D, Pederson L. Particle mesh ewald: An n log(n) method for ewald sums in large systems. J. Chem. Phys. 1993;98:10089–10092.
[64] Dinner AR, Lopez X, Karplus M. A charge-scaling method to treat solvent in qm/mm simulations. Theoretical Chemistry Accounts. 2003 Apr;109(3):118–124.
[65] Schaefer P, Riccardi D, Cui Q. Reliable treatment of electrostatics in combined qm/mm simulation of macromolecules. Journal of Chemical Physics. 2005 Jul;123(1):014905. [PubMed]
[66] Gregersen BA, York DM. Variational electrostatic projection (vep) methods for efficient modeling of the macromolecular electrostatic and solvation environment in activated dynamics simulations. Journal of Physical Chemistry B. 2005 Jan;109(1):536–556. [PubMed]
[67] Nagaoka M, Okuyama-Yoshida N, Yamabe T. Origin of the transition state on the free energy surface: Intramolecular proton transfer reaction of glycine in aqueous solution. Journal of Physical Chemistry A. 1998;102(42):8202–8208.
[68] Fleurat-Lessard P, Ziegler T. Tracing the minimum-energy path on the free-energy surface. Journal of Chemical Physics. 2005;123:084101. [PubMed]
[69] Maragliano L, Fischer A, Vanden-Eijnden E, Ciccotti G. String method in collective variables: Minimum free energy paths and isocommittor surfaces. JCP. 2006;125:024106. [PubMed]
[70] Higashi M, Hayashi S, Kato S. Transition state determination of enzyme reaction on free energy surface: Application to chorismate mutase. Chemical Physics Letters. 2007 Apr;437(4-6):293–297.
[71] Bolhuis PG, Dellago C, Chandler D. Reaction coordinates of biomolecular isomerization. Proceedings of the National Academy of Sciences of the United States of America. 2000 May;97(11):5877–5882. [PubMed]
[72] JAnsson H, Mills G, Jacobsen KW. Classical and Quantum Dynamics in Condensed Phase Simulations. World Scientific; Singapore: 1998. Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions.
[73] Ayala PY, Schlegel HB. A combined method for determining reaction paths, minima, and transition state geometries. Journal of Chemical Physics. 1997;107(2):375–384.
[74] Ren WNE,WQ, Vanden-Eijnden E. String method for the study of rare events. Phys. Rev. B. 2002;66:052301.
[75] Burger SK, Yang WT. Quadratic string method for determining the minimum-energy path based on multiobjective optimization. Journal of Chemical Physics. 2006 Feb;124(5):054109. [PubMed]
[76] Xie L, Liu HY, Yang WT. Adapting the nudged elastic band method for determining minimum-energy paths of chemical reactions in enzymes. Journal of Chemical Physics. 2004 May;120(17):8039–8052. [PubMed]
[77] Kastner J, Thiel S, Senn HM, Sherwood P, Thiel W. Exploiting qm/mm capabilities in geometry optimization: A microiterative approach using electrostatic embedding. Journal of Chemical Theory and Computation. 2007 May-Jun;3(3):1064–1072. [PubMed]
[78] Vreven T, Morokuma K, Farkas O, Schlegel HB, Frisch MJ. Geometry optimization with qm/mm, oniom, and other combined methods. i. microiterations and constraints. Journal of Computational Chemistry. 2003 Apr;24(6):760–769. [PubMed]
[79] Bolhuis PG, Dellago C, Geissler PL, Chandler D. Transition path sampling: throwing ropes over mountains in the dark. Journal of Physics-Condensed Matter. 2000 Feb;12(8A):A147–A152.
[80] Radhakrishnan R, Schlick T. Orchestration of cooperative events in dna synthesis and repair mechanism unraveled by transition pathsampling of dna polymerase b’s closing. Proc. Natl. Acad. Sci. USA. 2004;101:5970–5975. [PubMed]
[81] Lu ZY, Yang WT. Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations. Journal of Chemical Physics. 2004 Jul;121(1):89–100. [PubMed]
[82] Morita A, Kato S. Ab initio molecular orbital theory on intramolecular charge polarization: Effect of hydrogen abstraction on the charge sensitivity of aromatic and nonaromatic species. J. Am. Chem. Soc. 1997;119:4021–4032.
[83] Hu H, Lu ZY, Yang WT. Fitting molecular electrostatic potentials from quantum mechanical calculations. Journal of Chemical Theory and Computation. 2007 May-Jun;3(3):1004–1013. [PubMed]
[84] Jarzynski C. Nonequilibrium equality for free energy differences. Phys. Rev. Lett. 1997;78:2690–2693.
[85] Beveridge DL, DiCapua FM. Free energy via molecular simulation: Applications to chemical and biomolecular systems. Annu. Rev. Biophys. Biophys. Chem. 1989;18:431–492. [PubMed]
[86] Li GH, Zhang XD, Cui Q. Free energy perturbation calculations with combined qm/mm potentials complications, simplifications, and applications to redox potential calculations. Journal of Physical Chemistry B. 2003 Aug;107(33):8643–8653.
[87] Bentzien J, Muller RP, Florian J, Warshel A. Hybrid ab initio quantum mechanics molecular mechanics calculations of free energy surfaces for enzymatic reactions: The nucleophilic attack in subtilisin. Journal of Physical Chemistry B. 1998 Mar;102(12):2293–2301.
[88] Wood RH, Yezdimer EM, Sakane S, Barriocanal JA, Doren DJ. Free energies of solvation with quantum mechanical interaction energies from classical mechanical simulations. Journal of Chemical Physics. 1999;110(3):1329–1337.
[89] Blumberger J, Tavernelli I, Klein ML, Sprik M. Diabatic free energy curves and coordination fluctuations for the aqueous ag+/ag2+ redox couple: A biased born-oppenheimer molecular dynamics investigation. J. Chem. Phys. 2006;124:064507. [PubMed]
[90] Rosta E, Klahn M, Warshel A. Towards accurate ab initio qm/mm calculations of free-energy profiles of enzymatic reactions. Journal of Physical Chemistry B. 2006 Feb;110(6):2934–2941. [PubMed]
[91] Hu H, Yang WT. Dual-topology/dual-coordinate free-energy simulation using qm/mm force field. Journal of Chemical Physics. 2005 Jul;123(4):041102. [PubMed]
[92] Li GH, Cui Q. pk(a) calculations with qm/mm free energy perturbations. Journal of Physical Chemistry B. 2003 Dec;107(51):14521–14528.
[93] Kastner J, Thiel W. Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: “umbrella integration” Journal of Chemical Physics. 2005 Oct;123(14):144104. [PubMed]
[94] Laio A, Parrinello M. Escaping free-energy minima. Proceedings of the National Academy of Sciences USA. 2002;99:12562–12566. [PubMed]
[95] Hermans J. A simple analysis of noise and hysteresis in free energy simulations. JPC. 1991;95:9029–9032.
[96] Crespo A, Martí MA, Estrin DA, Roitberg AE. Multiple-steering qmmm calculation of the free energy profile in chorismate mutase. J. Am. Chem. Soc. 2005;127:6940–6941. [PubMed]
[97] Raugei S, Cascella M, Carloni P. A proficient enzyme: Insights on the mechanism of orotidine monophosphate decarboxylase from computer simulations. Journal of the American Chemical Society. 2004 Dec;126(48):15730–15737. [PubMed]
[98] Hu H, Yun RH, Hermans J. Reversibility of free energy simulations: slow growth may have a unique advantage. (with a note on use of ewald summation) Mol. Sim. 2002;28(1-2):67–80.
[99] Wang SL, Hu P, Zhang YK. Ab initio quantum mechanical/molecular mechanical molecular dynamics simulation of enzyme catalysis: The case of histone lysine methyltransferase set7/9. Journal of Physical Chemistry B. 2007 Apr;111(14):3758–3764. [PMC free article] [PubMed]
[100] Ishida T, Kato S. Theoretical perspectives on the reaction mechanism of serine proteases: The reaction free energy profiles of the acylation process. Journal of the American Chemical Society. 2003 Oct;125(39):12035–12048. [PubMed]
[101] Rod TH, Ryde U. Accurate qm/mm free energy calculations of enzyme reactions: Methylation by catechol o-methyltransferase. J. Chem. Theory Comput. 2005;1:1240–1251. [PubMed]
[102] Kastner J, Senn HM, Thiel S, Otte N, Thiel W. Qm/mm free-energy perturbation compared to thermodynamic integration and umbrella sampling: Application to an enzymatic reaction. Journal of Chemical Theory and Computation. 2006 Mar-Apr;2(2):452–461. [PubMed]
[103] Cui QA, Karplus M. Promoting modes and demoting modes in enzyme-catalyzed proton transfer reactions: A study of models and realistic systems. Journal of Physical Chemistry B. 2002 Aug;106(32):7927–7947.
[104] Wang ML, Lu ZY, Yang WT. Transmission coefficient calculation for proton transfer in triosephosphate isomerase based on the reaction path potential method. Journal of Chemical Physics. 2004 Jul;121(1):101–107. [PubMed]
[105] Wang ML, Lu ZY, Yang WT. Nuclear quantum effects on an enzymecatalyzed reaction with reaction path potential: Proton transfer in triosephosphate isomerase. Journal of Chemical Physics. 2006 Mar;124(12):124516. [PubMed]
[106] Lee JK, Houk KN. A proficient enzyme revisited: The predicted mechanism for orotidine monophosphate decarboxylase. Science. 1997 May;276(5314):942–945. [PubMed]
[107] Wu N, Mo YR, Gao JL, Pai EF. Electrostatic stress in catalysis: Structure and mechanism of the enzyme orotidine monophosphate decarboxylase. Proceedings of the National Academy of Sciences of the United States of America. 2000 Feb;97(5):2017–2022. [PubMed]
[108] Warshel A, Strajbl M, Villa J, Florian J. Remarkable rate enhancement of orotidine 5 ′-monophosphate decarboxylase is due to transition-state stabilization rather than to ground-state destabilization. Biochemistry. 2000 Dec;39(48):14728–14738. [PubMed]
[109] Lee TS, Chong LT, Chodera JD, Kollman PA. An alternative explanation for the catalytic proficiency of orotidine 5 ′-phosphate decarboxylase. Journal of the American Chemical Society. 2001 Dec;123(51):12837–12848. [PubMed]
[110] Lundberg M, Blomberg MRA, Siegbahn PEM. Density functional models of the mechanism for decarboxylation in orotidine decarboxylase. Journal of Molecular Modeling. 2002;8(4):119–130. [PubMed]
[111] Cisneros GA, Wang M, Silinski P, Fitzgerald MC, Yang WT. The protein backbone makes important contributions to 4-oxalocrotonate tautomerase enzyme catalysis: Understanding from theory and experiment. Biochemistry. 2004 Jun;43(22):6885–6892. [PubMed]
[112] Cisneros GA, Liu HY, Lu ZY, Yang WT. Reaction path determination for quantum mechanical/molecular mechanical modeling of enzyme reactions by combining first order and second order “chain-of-replicas” methods. Journal of Chemical Physics. 2005 Mar;122(11):114502. [PubMed]
[113] Cisneros GA, Wang M, Silinski P, Fitzgerald MC, Yang WT. Theoretical and experimental determination on two substrates turned over by 4-oxalocrotonate tautomerase. Journal of Physical Chemistry A. 2006 Jan;110(2):700–708. [PubMed]
[114] Tuttle T, Keinan E, Thiel W. Understanding the enzymatic activity of 4-oxalocrotonate tautomerase and its mutant analogues: A computational study. Journal of Physical Chemistry B. 2006 Oct;110(39):19685–19695. [PubMed]
[115] Bash PA, Field MJ, Davenport RC, Petsko GA, Ringe D, Karplus M. Computer-simulation and analysis of the reaction pathway of triosephosphate isomerase. Biochemistry. 1991 Jun;30(24):5826–5832. [PubMed]
[116] Aqvist J, Fothergill M. Computer simulation of the triosephosphate isomerase catalyzed reaction. Journal of Biological Chemistry. 1996 Apr;271(17):10010–10016. [PubMed]
[117] Alagona G, Ghio C, Kollman PA. The intramolecular mechanism for the second proton transfer in triosephosphate isomerase (tim): A qm/fe approach. Journal of Computational Chemistry. 2003 Jan;24(1):46–56. [PubMed]
[118] Cui Q, Karplus M. Quantum mechanics/molecular mechanics studies of triosephosphate isomerase-catalyzed reactions: Effect of geometry and tunneling on proton-transfer rate constants. Journal of the American Chemical Society. 2002 Mar;124(12):3093–3124. [PubMed]
[119] Guallar V, Jacobson M, McDermott A, Friesner RA. Computational modeling of the catalytic reaction in triosephosphate isomerase. Journal of Molecular Biology. 2004 Mar;337(1):227–239. [PubMed]
[120] Lee TS, Massova I, Kuhn B, Kollman PA. Qm and qm-fe simulations on reactions of relevance to enzyme catalysis: trypsin, catechol o-methyltransferase, beta-lactamase and pseudouridine synthase. Journal of the Chemical Society-Perkin Transactions. 2000;2(3):409–415.
[121] Kuhn B, Kollman PA. Qm-fe and molecular dynamics calculations on catechol omethyltransferase: Free energy of activation in the enzyme and in aqueous solution and regioselectivity of the enzyme-catalyzed reaction. Journal of the American Chemical Society. 2000 Mar;122(11):2586–2596.
[122] Donini O, Darden T, Kollman PA. Qm-fe calculations of aliphatic hydrogen abstraction in citrate synthase and in solution: Reproduction of the effect of enzyme catalysis and demonstration that an enolate rather than an enol is formed. Journal of the American Chemical Society. 2000 Dec;122(49):12270–12280.
[123] Alhambra C, Gao JL, Corchado JC, Villa J, Truhlar DG. Quantum mechanical dynamical effects in an enzyme-catalyzed proton transfer reaction. Journal of the American Chemical Society. 1999 Mar;121(10):2253–2258.
[124] Woodcock HL, Hodoscek M, Sherwood P, Lee YS, Schaefer HF, Brooks BR. Exploring the quantum mechanical/molecular mechanical replica path method: a pathway optimization of the chorismate to prephenate claisen rearrangement catalyzed by chorismate mutase. Theoretical Chemistry Accounts. 2003 Apr;109(3):140–148.
[125] Lee YS, Worthington SE, Krauss M, Brooks BR. Reaction mechanism of chorismate mutase studied by the combined potentials of quantum mechanics and molecular mechanics. Journal of Physical Chemistry B. 2002 Nov;106(46):12059–12065.
[126] Woodcock HL, Hodoscek M, Gilbert ATB, Gill PMW, Schaefer HF, Brooks BR. Interfacing q-chem and charmm to perform qm/mm reaction path calculations. Journal of Computational Chemistry. 2007 Jul;28(9):1485–1502. [PubMed]
[127] Hermann JC, Hensen C, Ridder L, Mulholland AJ, Holtje HD. Mechanisms of antibiotic resistance: Qm/mm modeling of the acylation reaction of a class a beta-lactamase with benzylpenicillin. Journal of the American Chemical Society. 2005 Mar;127(12):4454–4465. [PubMed]
[128] Diaz N, Suarez D, Sordo TL, Merz KM. Acylation of class a beta-lactamases by penicillins: A theoretical examination of the role of serine 130 and the beta-lactam carboxylate group. Journal of Physical Chemistry B. 2001 Nov;105(45):11302–11313.
[129] Diaz N, Sordo TL, Merz KM, Suarez D. Insights into the acylation mechanism of class a beta-lactamases from molecular dynamics simulations of the tem-1 enzyme complexed with benzylpenicillin. Journal of the American Chemical Society. 2003 Jan;125(3):672–684. [PubMed]
[130] Zhang YK, Kua J, McCammon JA. Role of the catalytic triad and oxyanion hole in acetylcholinesterase catalysis: An ab initio qm/mm study. J. Am. Chem. Soc. 2002;124:10572–10577. [PubMed]
[131] Zhang YK, Kua J, McCammon JA. Influence of structural fluctuation on enzyme reaction energy barriers in combined quantum mechanical/molecular mechanical studies. J. Phys. Chem. B. 2003;107:4459–4463.
[132] Cheng YH, Zhang YK, McCammon JA. How does activation loop phosphorylation modulate catalytic activity in the camp-dependent protein kinase: A theoretical study. Protein Science. 2006 Apr;15(4):672–683. [PubMed]
[133] Cheng YH, Zhang YK, McCammon JA. How does the camp-dependent protein kinase catalyze the phosphorylation reaction: An ab initio qm/mm study. Journal of the American Chemical Society. 2005 Feb;127(5):1553–1562. [PubMed]
[134] Wang LH, Yu XY, Hu P, Broyde S, Zhang YK. A water-mediated and substrate-assisted catalytic mechanism for sulfolobus solfataricus dna polymerase iv. Journal of the American Chemical Society. 2007 Apr;129(15):4731–4737. [PMC free article] [PubMed]
[135] Corminboeut C, Hu P, Tuckerman ME, Zhang YK. Unexpected deacetylation mechanism suggested by a density functional theory qm/mm study of histonedeacetylase-like protein. Journal of the American Chemical Society. 2006 Apr;128(14):4530–4531. [PubMed]
[136] Hu P, Zhang YK. Catalytic mechanism and product specificity of the histone lysine methyltransferase set7/9: An ab initio qm/mm-fe study with multiple initial structures. Journal of the American Chemical Society. 2006 Feb;128(4):1272–1278. [PubMed]
[137] Xiao CY, Zhang YK. Catalytic mechanism and metal specificity of bacterial peptide deformylase: A density functional theory qm/mm study. Journal of Physical Chemistry B. 2007 Jun;111(22):6229–6235. [PubMed]
[138] Karambelkar VV, Xiao CY, Zhang YK, Sarjeant AAN, Goldberg DP. Geometric preferences in iron(ii) and zinc(ii) model complexes of peptide deformylase. Inorganic Chemistry. 2006 Feb;45(4):1409–1411. [PMC free article] [PubMed]
[139] Gherman BF, Baik MH, Lippard SJ, Friesner RA. Dioxygen activation in methane monooxygenase: A theoretical study. Journal of the American Chemical Society. 2004 Mar;126(9):2978–2990. [PubMed]
[140] Gherman BF, Lippard SJ, Friesner RA. Substrate hydroxylation in methane monooxygenase: Quantitative modeling via mixed quantum mechanics/molecular mechanics techniques. Journal of the American Chemical Society. 2005 Jan;127(3):1025–1037. [PubMed]
[141] Guallar V, Friesner RA. Cytochrome p450cam enzymatic catalysis cycle: A quantum mechanics molecular mechanics study. Journal of the American Chemical Society. 2004 Jul;126(27):8501–8508. [PubMed]
[142] Bucher D, Raugei S, Guidoni L, Dal Peraro M, Rothlisberger U, Carloni P, Klein ML. Polarization effects and charge transfer in the kcsa potassium channel. Biophysical Chemistry. 2006 Dec;124(3):292–301. [PubMed]
[143] Cascella M, Micheletti C, Rothlisberger U, Carloni P. Evolutionarily conserved functional mechanics across pepsin-like and retroviral aspartic proteases. Journal of the American Chemical Society. 2005 Mar;127(11):3734–3742. [PubMed]
[144] Piana S, Bucher D, Carloni P, Rothlisberger U. Reaction mechanism of hiv-1 protease by hybrid carparrinello/classical md simulations. Journal of Physical Chem-istry B. 2004 Jul;108(30):11139–11149.
[145] Sulpizi M, Rothlisberger U, Laio A, Cattaneo A, Carloni P. Reaction mechanism of caspases: Insights from qm/mm car-parrinello simulations. Biophysical Journal. 2002 Jan;82(1):359A–360A. [PubMed]
[146] Field MJ. Hybrid quantum mechanical/molecular mechanical fluctuating charge models for condensed phase simulations. Mol. Phys. 1997;91:835–835.
[147] Thompson MA. Qm/mmpol: a consistent model for solute/solvent polarization. application to the aqueous solvation and spectroscopy of formaldehyde, acetaldehyde, and acetone. J. Phys. Chem. 1996;100:14492–14507.
[148] Bakowies D, Thiel D. Hybrid models for combined quantum mechanical and molecular mechanical approaches. J. Phys. Chem. 1996;100:10580–10594.
[149] Gao JL. Energy components of aqueous solution: Insight from hybrid qm/mm simulations using a polarizable solvent model. Journal of Computational Chemistry. 1997 Jun;18(8):1061–1071.
[150] Dupuis M, Aida M, Kawashima Y, Hirao K. A polarizable mixed hamiltonian model of electronic structure for micro-solvated excited states. i. energy and gradients formulation and application to formaldehyde ((1)a(2)) J. Chem. Phys. 2002;117:1242–1255.