Electrostatically Embedded Multiconfiguration Molecular Mechanics Based on the Combined Density Functional and Molecular Mechanical Method

Potential energy interpolation with target-customized weighting coordinates: application to excited-state dynamics of photoactive yellow protein chromophore in water

Interpolation of potential energy surfaces (PESs) can provide a practical route to performing molecular dynamics simulations with a reliability matching a high-level quantum chemical calculation. An obstacle to its widespread use is perhaps the lack of general and optimal interpolation settings that can be applied in a black-box manner for any given molecular system. How to set up the weights for interpolation is one such task, and we still need to diversify the approaches in order to treat various systems. Here, we develop a new interpolation weighting scheme, which allows us to choose the weighting coordinates in a system-specific manner, by amplifying the contribution from specific internal coordinates. The new weighting scheme with an appropriate selection of coordinates is proved to be effective in reducing the interpolation error along the reaction pathway. As a demonstration, we consider the photoactive yellow protein chromophore system, as it constitutes itself as an interesting target that bears long-standing questions related to excited-state dynamics inside protein environments. We build its two-state diabatic interpolated PES with the new weighting scheme. We indeed see the utility of our scheme by conducting nonadiabatic molecular dynamics simulations with the required semi-global PES based on a limited number of data points.

Modeling Charge Flux by Interpolating Atomic Partial Charges

Although the charge flux effect or the geometric dependence of the atomic partial charges have been known for a long time, how it can be effectively handled is not yet established. Here, we present a charge interpolation scheme as a new general tool for representing the charge flux in an analytically well-defined manner. By applying it to the anionic GFP chromophore with the diabatically represented atomic charges, we show that the charge interpolation provides a substantial improvement on the accuracy of the geometry-dependent changes in the molecular dipole moments in the gas phase. We also test the scheme toward describing the electrostatic term in the solvation energy in the aqueous environment and observe that it is also improved but that the extent of the improvement is somewhat limited. We show that the remaining errors can be largely corrected by introducing atomic polarizabilities. Overall, our results show that charge interpolation is an amenable approach for describing the charge flux effect and that its description in the condensed phase should be accompanied by proper treatments of polarization effects.

Structural properties determining low K+ affinity of the selectivity filter in the TWIK1 K+ channel

Canonical K⁺ channels are tetrameric and highly K⁺-selective, whereas two-pore-domain K⁺ (K2P) channels form dimers, but with a similar pore architecture. A two-pore-domain potassium channel TWIK1 (KCNK1 or K2P1) allows permeation of Na⁺ and other monovalent ions, resulting mainly from the presence of Thr-118 in the P1 domain. However, the mechanistic basis for this reduced selectivity is unclear. Using ion-exchange-induced difference IR spectroscopy, we analyzed WT TWIK1 and T118I (highly K⁺-selective) and L228F (substitution in the P2 domain) TWIK1 variants and found that in the presence of K⁺ ions, WT and both variants exhibit an amide-I band at 1680 cm^-1 This band corresponds to interactions of the backbone carbonyls in the selectivity filter with K⁺, a feature very similar to that of the canonical K⁺ channel KcsA. Computational analysis indicated that the relatively high frequency for the amide-I band is well explained by impairment of hydrogen bond formation with water molecules. Moreover, concentration-dependent spectral changes indicated that the K⁺ affinity of the WT selectivity filter was much lower than those of the variants. Furthermore, only the variants displayed a higher frequency shift of the 1680-cm^-1 band upon changes from K⁺ to Rb⁺ or Cs⁺ conditions. High-speed atomic force microscopy disclosed that TWIK1's surface morphology largely does not change in K⁺ and Na⁺ solutions. Our results reveal the local conformational changes of the TWIK1 selectivity filter and suggest that the amide-I bands may be useful "molecular fingerprints" for assessing the properties of other K⁺ channels.

Theoretical approaches for dynamical ordering of biomolecular systems

BACKGROUND

Living systems are characterized by the dynamic assembly and disassembly of biomolecules. The dynamical ordering mechanism of these biomolecules has been investigated both experimentally and theoretically. The main theoretical approaches include quantum mechanical (QM) calculation, all-atom (AA) modeling, and coarse-grained (CG) modeling. The selected approach depends on the size of the target system (which differs among electrons, atoms, molecules, and molecular assemblies). These hierarchal approaches can be combined with molecular dynamics (MD) simulation and/or integral equation theories for liquids, which cover all size hierarchies.

SCOPE OF REVIEW

We review the framework of quantum mechanical/molecular mechanical (QM/MM) calculations, AA MD simulations, CG modeling, and integral equation theories. Applications of these methods to the dynamical ordering of biomolecular systems are also exemplified.

MAJOR CONCLUSIONS

The QM/MM calculation enables the study of chemical reactions. The AA MD simulation, which omits the QM calculation, can follow longer time-scale phenomena. By reducing the number of degrees of freedom and the computational cost, CG modeling can follow much longer time-scale phenomena than AA modeling. Integral equation theories for liquids elucidate the liquid structure, for example, whether the liquid follows a radial distribution function.

GENERAL SIGNIFICANCE

These theoretical approaches can analyze the dynamic behaviors of biomolecular systems. They also provide useful tools for exploring the dynamic ordering systems of biomolecules, such as self-assembly. This article is part of a Special Issue entitled "Biophysical Exploration of Dynamical Ordering of Biomolecular Systems" edited by Dr. Koichi Kato.

Direct Learning Hidden Excited State Interaction Patterns from ab initio Dynamics and Its Implication as Alternative Molecular Mechanism Models

The excited states of polyatomic systems are rather complex, and often exhibit meta-stable dynamical behaviors. Static analysis of reaction pathway often fails to sufficiently characterize excited state motions due to their highly non-equilibrium nature. Here, we proposed a time series guided clustering algorithm to generate most relevant meta-stable patterns directly from ab initio dynamic trajectories. Based on the knowledge of these meta-stable patterns, we suggested an interpolation scheme with only a concrete and finite set of known patterns to accurately predict the ground and excited state properties of the entire dynamics trajectories, namely, the prediction with ensemble models (PEM). As illustrated with the example of sinapic acids, The PEM method does not require any training data beyond the clustering algorithm, and the estimation error for both ground and excited state is very close, which indicates one could predict the ground and excited state molecular properties with similar accuracy. These results may provide us some insights to construct molecular mechanism models with compatible energy terms as traditional force fields.

Perspective: Defining and quantifying the role of dynamics in enzyme catalysis

Enzymes control chemical reactions that are key to life processes, and allow them to take place on the time scale needed for synchronization between the relevant reaction cycles. In addition to general interest in their biological roles, these proteins present a fundamental scientific puzzle, since the origin of their tremendous catalytic power is still unclear. While many different hypotheses have been put forward to rationalize this, one of the proposals that has become particularly popular in recent years is the idea that dynamical effects contribute to catalysis. Here, we present a critical review of the dynamical idea, considering all reasonable definitions of what does and does not qualify as a dynamical effect. We demonstrate that no dynamical effect (according to these definitions) has ever been experimentally shown to contribute to catalysis. Furthermore, the existence of non-negligible dynamical contributions to catalysis is not supported by consistent theoretical studies. Our review is aimed, in part, at readers with a background in chemical physics and biophysics, and illustrates that despite a substantial body of experimental effort, there has not yet been any study that consistently established a connection between an enzyme's conformational dynamics and a significant increase in the catalytic contribution of the chemical step. We also make the point that the dynamical proposal is not a semantic issue but a well-defined scientific hypothesis with well-defined conclusions.

Combined Electrostatically Embedded Multiconfiguration Molecular Mechanics and Molecular Mechanical Method: Application to Molecular Dynamics Simulation of a Chemical Reaction in Aqueous Solution with Hybrid Density Functional Theory

We here combine the electrostatically embedded multiconfiguration molecular mechanics (EE-MCMM) method for generating global potential energy surfaces in the presence of an electrostatic potential with molecular mechanics (MM). The resulting EE-MCMM/MM method is illustrated by applying it to carry out a molecular dynamics simulation for the symmetric bimolecular reaction Cl(-) + CH3Cl' → ClCH3 + Cl'(-) in aqueous solution with hybrid density functional theory as the quantum mechanical level. The potential of mean force is calculated, and the free energy barrier is found to be 25.3 kcal/mol, which is in good agreement with previous work. The advantage of the combined EE-MCMM and MM method is that the number of quantum mechanical calculations required for the active subsystem is very small compared to straight direct dynamics.

Constructing polyatomic potential energy surfaces by interpolating diabatic Hamiltonian matrices with demonstration on green fluorescent protein chromophore

Simulating molecular dynamics directly on quantum chemically obtained potential energy surfaces is generally time consuming. The cost becomes overwhelming especially when excited state dynamics is aimed with multiple electronic states. The interpolated potential has been suggested as a remedy for the cost issue in various simulation settings ranging from fast gas phase reactions of small molecules to relatively slow condensed phase dynamics with complex surrounding. Here, we present a scheme for interpolating multiple electronic surfaces of a relatively large molecule, with an intention of applying it to studying nonadiabatic behaviors. The scheme starts with adiabatic potential information and its diabatic transformation, both of which can be readily obtained, in principle, with quantum chemical calculations. The adiabatic energies and their derivatives on each interpolation center are combined with the derivative coupling vectors to generate the corresponding diabatic Hamiltonian and its derivatives, and they are subsequently adopted in producing a globally defined diabatic Hamiltonian function. As a demonstration, we employ the scheme to build an interpolated Hamiltonian of a relatively large chromophore, para-hydroxybenzylidene imidazolinone, in reference to its all-atom analytical surface model. We show that the interpolation is indeed reliable enough to reproduce important features of the reference surface model, such as its adiabatic energies and derivative couplings. In addition, nonadiabatic surface hopping simulations with interpolation yield population transfer dynamics that is well in accord with the result generated with the reference analytic surface. With these, we conclude by suggesting that the interpolation of diabatic Hamiltonians will be applicable for studying nonadiabatic behaviors of sizeable molecules.

Theoretical study on excited states of bacteriochlorophyll a in solutions with density functional assessment

The excited-state properties of bacteriochlorophyll (BChl) a in triethylamine, 1-propanol, and methanol are investigated with the time-dependent density functional theory by using the quantum mechanical and molecular mechanical reweighting free energy self-consistant field method. It is found that no prevalent density functionals can reproduce the experimental excited-state properties, i.e., the absorption and reorganization energies, of BChl a in the solutions. The parameter μ in the range-separated hybrid functional is therefore optimized to reproduce the differences of the absorption energies in the solutions. We examine the origin of the differences of the absorption energies in the solutions and find that sensitive balance between contributions of structural changes and solute-solvent interactions determines the differences. The accurate description of the excitation with the density functional with the adjusted parameter is therefore essential to the understanding of the excited-state properties of BChl a in proteins and also the mechanism of the photosynthetic systems.

Reaction Path Following with Sparse Interpolation

Computing the potential energy of an N-atom molecule is an expensive optimization process of 3N - 6 molecular coordinates, so following reaction pathways as a function of all 3N - 6 coordinates is unfeasible for large molecules. In this paper, we present a method that isolates d < 3N - 6 molecular coordinates and continuously follows reaction paths on d-dimensional potential energy surfaces approximated by a Smolyak's sparse grid interpolation algorithm.1 Compared to dense grids, sparse grids efficiently improve the ratio of invested storage and computing time to approximation accuracy and thus allow one to increase the number of coordinates d in molecular reaction path following simulations. Furthermore, evaluation of the interpolant is much less expensive than the evaluation of the actual energy function, so our technique offers a computationally efficient way to simulate reaction paths on ground and excited state potential energy surfaces. To demonstrate the capabilities of our method, we present simulation results for the isomerization of 2-butene with two, three, and six degrees of freedom.

Concerted hydrogen atom and electron transfer mechanism for catalysis by lysine-specific demethylase

We calculate the free energy profile for the postulated hydride transfer reaction mechanism for the catalysis of lysine demethylation by lysine-specific demethylase LSD1. The potential energy surface is obtained by using combined electrostatically embedded multiconfiguration molecular mechanics (EE-MCMM) and single-configuration molecular mechanics (MM). We employ a constant valence bond coupling term to obtain analytical energies and gradients of the EE-MCMM subsystem, which contains 45 quantum mechanics (QM) atoms and which is parametrized with density functional calculations employing specific reaction parameters obtained by matching high-level wave function calculations. In the MM region, we employ the Amber ff03 and TIP3P force fields. The free energy of activation at 300 K is calculated by molecular dynamics (MD) umbrella sampling on a system with 102,090 atoms as the maximum of the free energy profile along the reaction coordinate as obtained by the weighted histogram analysis method with 17 umbrella sampling windows. This yields a free energy of activation of only 10 kcal/mol, showing that the previously postulated direct hydride transfer reaction mechanism is plausible, although we find that it is better interpreted as a concerted transfer of a hydrogen atom and an electron.

Multistructural variational transition state theory: kinetics of the hydrogen abstraction from carbon-2 of 2-methyl-1-propanol by hydroperoxyl radical including all structures and torsional anharmonicity

We calculated the forward and reverse rate constants of the hydrogen abstraction reaction from carbon-2 of 2-methyl-1-propanol by hydroperoxyl radical over the temperature range 250-2400 K by using multistructural canonical variational transition state theory (MS-CVT) including both multiple-structure and torsional potential anharmonicity effects by the multistructural torsional anharmonicity (MS-T) method. In these calculations, multidimensional tunneling (MT) probabilities used to compute the tunneling transmission coefficients were evaluated by the small-curvature tunneling (SCT) approximation. Comparison with the rate constants obtained by the single-structural harmonic oscillator (SS-HO) approximation shows that multistructural anharmonicity increases the forward rate constants for all temperatures, but the reverse rate constants are reduced for temperatures lower than 430 K and increased for higher temperatures. The neglect of multistructural torsional anharmonicity would lead to errors of factors of 1.5, 8.8, and 13 at 300, 1000, and 2400 K, respectively, for the forward reaction, and would lead to errors of factors of 0.76, 3.0, and 6.0, respectively, at these temperatures for the reverse reaction.

Computationally Efficient Multiconfigurational Reactive Molecular Dynamics

It is a computationally demanding task to explicitly simulate the electronic degrees of freedom in a system to observe the chemical transformations of interest, while at the same time sampling the time and length scales required to converge statistical properties and thus reduce artifacts due to initial conditions, finite-size effects, and limited sampling. One solution that significantly reduces the computational expense consists of molecular models in which effective interactions between particles govern the dynamics of the system. If the interaction potentials in these models are developed to reproduce calculated properties from electronic structure calculations and/or ab initio molecular dynamics simulations, then one can calculate accurate properties at a fraction of the computational cost. Multiconfigurational algorithms model the system as a linear combination of several chemical bonding topologies to simulate chemical reactions, also sometimes referred to as "multistate". These algorithms typically utilize energy and force calculations already found in popular molecular dynamics software packages, thus facilitating their implementation without significant changes to the structure of the code. However, the evaluation of energies and forces for several bonding topologies per simulation step can lead to poor computational efficiency if redundancy is not efficiently removed, particularly with respect to the calculation of long-ranged Coulombic interactions. This paper presents accurate approximations (effective long-range interaction and resulting hybrid methods) and multiple-program parallelization strategies for the efficient calculation of electrostatic interactions in reactive molecular simulations.

Interpolated mechanics-molecular mechanics study of internal rotation dynamics of the chromophore unit in blue fluorescent protein and its variants

The excited state dynamics of the blue fluorescent protein (BFP) and its variants, azurite, EBFP1.2, and EBFP2.0, are studied using molecular dynamics (MD) simulations on potential energy surfaces (PESs) generated with the interpolated mechanics-molecular mechanics (IM/MM) scheme. This IM/MM strategy adopts the interpolated PES for an important area of the complex and the conventional force field for the remaining part. We focus on the internal rotation dynamics of the chromophore unit, which is directly related to its fluorescence property, and analyze the time evolutions of the nonrotated chromophore fractions based on trajectories over 10 μs of aggregate simulation time. The characteristics obtained from the calculated time progresses of the nonrotated chromophore fractions in BFP and other variants agree well with experimentally observed properties. The results show that the MD simulation with an IM/MM potential is an attractive approach for studying excited state dynamics of fluorescent proteins in consideration of its efficiency and reliability. We also attempt to investigate the detailed roles that the mutated residues play in delaying the excited state chromophore twisting and thus improving the fluorescence property, and discuss the contributions by the Coulombic and the steric interactions between the chromophore and the mutated residues.

Proton Transfer Studied Using a Combined Ab Initio Reactive Potential Energy Surface with Quantum Path Integral Methodology

The rates of intramolecular proton transfer are calculated on a full-dimensional reactive electronic potential energy surface that incorporates high level ab initio calculations along the reaction path and by using classical Transition State theory, Path-Integral Quantum Transition State Theory, and the Quantum Instanton approach. The specific example problem studied is malonaldehyde. Estimates of the kinetic isotope effect using the latter two methods are found to be in reasonable agreement with each other. Improvements and extensions of this practical, yet chemically accurate framework for the calculations of quantized, reactive dynamics are also discussed.

Gradient-based multiconfiguration Shepard interpolation for generating potential energy surfaces for polyatomic reactions

This paper describes and illustrates a way to construct multidimensional representations of reactive potential energy surfaces (PESs) by a multiconfiguration Shepard interpolation (MCSI) method based only on gradient information, that is, without using any Hessian information from electronic structure calculations. MCSI, which is called multiconfiguration molecular mechanics (MCMM) in previous articles, is a semiautomated method designed for constructing full-dimensional PESs for subsequent dynamics calculations (classical trajectories, full quantum dynamics, or variational transition state theory with multidimensional tunneling). The MCSI method is based on Shepard interpolation of Taylor series expansions of the coupling term of a 2 x 2 electronically diabatic Hamiltonian matrix with the diagonal elements representing nonreactive analytical PESs for reactants and products. In contrast to the previously developed method, these expansions are truncated in the present version at the first order, and, therefore, no input of electronic structure Hessians is required. The accuracy of the interpolated energies is evaluated for two test reactions, namely, the reaction OH+H(2)-->H(2)O+H and the hydrogen atom abstraction from a model of alpha-tocopherol by methyl radical. The latter reaction involves 38 atoms and a 108-dimensional PES. The mean unsigned errors averaged over a wide range of representative nuclear configurations (corresponding to an energy range of 19.5 kcal/mol in the former case and 32 kcal/mol in the latter) are found to be within 1 kcal/mol for both reactions, based on 13 gradients in one case and 11 in the other. The gradient-based MCMM method can be applied for efficient representations of multidimensional PESs in cases where analytical electronic structure Hessians are too expensive or unavailable, and it provides new opportunities to employ high-level electronic structure calculations for dynamics at an affordable cost.

At the dawn of the 21st century: Is dynamics the missing link for understanding enzyme catalysis?

Enzymes play a key role in almost all biological processes, accelerating a variety of metabolic reactions as well as controlling energy transduction, the transcription, and translation of genetic information, and signaling. They possess the remarkable capacity to accelerate reactions by many orders of magnitude compared to their uncatalyzed counterparts, making feasible crucial processes that would otherwise not occur on biologically relevant timescales. Thus, there is broad interest in understanding the catalytic power of enzymes on a molecular level. Several proposals have been put forward to try to explain this phenomenon, and one that has rapidly gained momentum in recent years is the idea that enzyme dynamics somehow contributes to catalysis. This review examines the dynamical proposal in a critical way, considering basically all reasonable definitions, including (but not limited to) such proposed effects as "coupling between conformational and chemical motions," "landscape searches" and "entropy funnels." It is shown that none of these proposed effects have been experimentally demonstrated to contribute to catalysis, nor are they supported by consistent theoretical studies. On the other hand, it is clarified that careful simulation studies have excluded most (if not all) dynamical proposals. This review places significant emphasis on clarifying the role of logical definitions of different catalytic proposals, and on the need for a clear formulation in terms of the assumed potential surface and reaction coordinate. Finally, it is pointed out that electrostatic preorganization actually accounts for the observed catalytic effects of enzymes, through the corresponding changes in the activation free energies.

Quantum tunneling in testosterone 6beta-hydroxylation by cytochrome P450: reaction dynamics calculations employing multiconfiguration molecular-mechanical potential energy surfaces

Testosterone hydroxylation is a prototypical reaction of human cytochrome P450 3A4, which metabolizes about 50% of oral drugs on the market. Reaction dynamics calculations were carried out for the testosterone 6beta-hydrogen abstraction and the 6beta-d(1)-testosterone 6beta-duterium abstraction employing a model that consists of the substrate and the active oxidant compound I. The calculations were performed at the level of canonical variational transition state theory with multidimensional tunneling and were based on a semiglobal full-dimensional potential energy surface generated by the multiconfiguration molecular mechanics technique. The tunneling coefficients were found to be around 3, indicating substantial contributions by quantum tunneling. However, the tunneling made only modest contributions to the kinetic isotope effects. The kinetic isotope effects were computed to be about 2 in the doublet spin state and about 5 in the quartet spin state.

The EVB as a quantitative tool for formulating simulations and analyzing biological and chemical reactions

Recent years have seen dramatic improvements in computer power, allowing ever more challenging problems to be approached. In light of this, it is imperative to have a quantitative model for examining chemical reactivity, both in the condensed phase and in solution, as well as to accurately quantify physical organic chemistry (particularly as experimental approaches can often be inconclusive). Similarly, computational approaches allow for great progress in studying enzyme catalysis, as they allow for the separation of the relevant energy contributions to catalysis. Due to the complexity of the problems that need addressing, there is a need for an approach that can combine reliability with an ability to capture complex systems in order to resolve long-standing controversies in a unique way. Herein, we will demonstrate that the empirical valence bond (EVB) approach provides a powerful way to connect the classical concepts of physical organic chemistry to the actual energies of enzymatic reactions by means of computation. Additionally, we will discuss the proliferation of this approach, as well as attempts to capture its basic chemistry and repackage it under different names. We believe that the EVB approach is the most powerful tool that is currently available for studies of chemical processes in the condensed phase in general and enzymes in particular, particularly when trying to explore the different proposals about the origin of the catalytic power of enzymes.

On unjustifiably misrepresenting the EVB approach while simultaneously adopting it

In recent years, the EVB has become a widely used tool in the QM/MM modeling of reactions in condensed phases and in biological systems, with ever increasing popularity. However, despite the fact that its power and validity have been repeatedly established since 1980, a recent work (Valero, R.; et al. J. Chem. Theory Comput. 2009, 5, 1) has strongly criticized this approach, while not discussing the fact that one of the authors is effectively using it himself for both gas-phase and solution studies. Here, we have responded to the most serious unjustified assertions of that paper, covering both the more problematic aspects of that work and the more complex scientific aspects. Additionally, we have demonstrated that the poor EVB results shown in Valero et al. which where presented as verification of the unreliability of the EVB model were in fact obtained by the use of incorrect parameters, without comparing to the correct surface obtained by our program.

Theoretical studies of RNA catalysis: hybrid QM/MM methods and their comparison with MD and QM

Hybrid QM/MM methods combine the rigor of quantum mechanical (QM) calculations with the low computational cost of empirical molecular mechanical (MM) treatment allowing to capture dynamic properties to probe critical atomistic details of enzyme reactions. Catalysis by RNA enzymes (ribozymes) has only recently begun to be addressed with QM/MM approaches and is thus still a field under development. This review surveys methodology as well as recent advances in QM/MM applications to RNA mechanisms, including those of the HDV, hairpin, and hammerhead ribozymes, as well as the ribosome. We compare and correlate QM/MM results with those from QM and/or molecular dynamics (MD) simulations, and discuss scope and limitations with a critical eye on current shortcomings in available methodologies and computer resources. We thus hope to foster mutual appreciation and facilitate collaboration between experimentalists and theorists to jointly advance our understanding of RNA catalysis at an atomistic level.

Progress in ab initio QM/MM free-energy simulations of electrostatic energies in proteins: accelerated QM/MM studies of pKa, redox reactions and solvation free energies

Hybrid quantum mechanical/molecular mechanical (QM/MM) approaches have been used to provide a general scheme for chemical reactions in proteins. However, such approaches still present a major challenge to computational chemists, not only because of the need for very large computer time in order to evaluate the QM energy but also because of the need for proper computational sampling. This review focuses on the sampling issue in QM/MM evaluations of electrostatic energies in proteins. We chose this example since electrostatic energies play a major role in controlling the function of proteins and are key to the structure-function correlation of biological molecules. Thus, the correct treatment of electrostatics is essential for the accurate simulation of biological systems. Although we will be presenting different types of QM/MM calculations of electrostatic energies (and related properties) here, our focus will be on pKa calculations. This reflects the fact that pKa's of ionizable groups in proteins provide one of the most direct benchmarks for the accuracy of electrostatic models of macromolecules. While pKa calculations by semimacroscopic models have given reasonable results in many cases, existing attempts to perform pKa calculations using QM/MM-FEP have led to discrepancies between calculated and experimental values. In this work, we accelerate our QM/MM calculations using an updated mean charge distribution and a classical reference potential. We examine both a surface residue (Asp3) of the bovine pancreatic trypsin inhibitor and a residue buried in a hydrophobic pocket (Lys102) of the T4-lysozyme mutant. We demonstrate that, by using this approach, we are able to reproduce the relevant side chain pKa's with an accuracy of 3 kcal/mol. This is well within the 7 kcal/mol energy difference observed in studies of enzymatic catalysis, and is thus sufficient accuracy to determine the main contributions to the catalytic energies of enzymes. We also provide an overall perspective of the potential of QM/MM calculations in general evaluations of electrostatic free energies, pointing out that our approach should provide a very powerful and accurate tool to predict the electrostatics of not only solution but also enzymatic reactions, as well as the solvation free energies of even larger systems, such as nucleic acid bases incorporated into DNA.

Perspective on Diabatic Models of Chemical Reactivity as Illustrated by the Gas-Phase S(N)2 Reaction of Acetate Ion with 1,2-Dichloroethane

Diabatic models are widely employed for studying chemical reactivity in condensed phases and enzymes, but there has been little discussion of the pros and cons of various diabatic representations for this purpose. Here we discuss and contrast six different schemes for computing diabatic potentials for a charge rearrangement reaction. They include (i) the variational diabatic configurations (VDC) constructed by variationally optimizing individual valence bond structures and (ii) the consistent diabatic configurations (CDC) obtained by variationally optimizing the ground-state adiabatic energy, both in the nonorthogonal molecular orbital valence bond (MOVB) method, along with the orthogonalized (iii) VDC-MOVB and (iv) CDC-MOVB models. In addition, we consider (v) the fourfold way (based on diabatic molecular orbitals and configuration uniformity), and (vi) empirical valence bond (EVB) theory. To make the considerations concrete, we calculate diabatic electronic states and diabatic potential energies along the reaction path that connects the reactant and the product ion-molecule complexes of the gas-phase bimolecular nucleophilic substitution (S(N)2) reaction of 1,2-dichloethane (DCE) with acetate ion, which is a model reaction corresponding to the reaction catalyzed by haloalkane dehalogenase. We utilize ab initio block-localized molecular orbital theory to construct the MOVB diabatic states and ab initio multi-configuration quasidegenerate perturbation theory to construct the fourfold-way diabatic states; the latter are calculated at reaction path geometries obtained with the M06-2X density functional. The EVB diabatic states are computed with parameters taken from the literature. The MOVB and fourfold-way adiabatic and diabatic potential energy profiles along the reaction path are in qualitative but not quantitative agreement with each other. In order to validate that these wave-function-based diabatic states are qualitatively correct, we show that the reaction energy and barrier for the adiabatic ground state, obtained with these methods, agree reasonably well with the results of high-level calculations using the composite G3SX and G3SX(MP3) methods and the BMC-CCSD multi-coefficient correlation method. However, a comparison of the EVB gas-phase adiabatic ground-state reaction path with those obtained from MOVB and with the fourfold way reveals that the EVB reaction path geometries show a systematic shift towards the products region, and that the EVB lowest-energy path has a much lower barrier. The free energies of solvation and activation energy in water reported from dynamical calculations based on EVB also imply a low activation barrier in the gas phase. In addition, calculations of the free energy of solvation using the recently proposed SM8 continuum solvation model with CM4M partial atomic charges lead to an activation barrier in reasonable agreement with experiment only when the geometries and the gas-phase barrier are those obtained from electronic structure calculations, i.e., methods i-v. These comparisons show the danger of basing the diabatic states on molecular mechanics without the explicit calculation of electronic wave functions. Furthermore, comparison of schemes i-v with one another shows that significantly different quantitative results can be obtained by using different methods for extracting diabatic states from wave function calculations, and it is important for each user to justify the choice of diabatization method in the context of its intended use.