Accelerated computation of free energy profile at ab initio quantum mechanical/molecular mechanical accuracy via a semi-empirical reference potential. II. Recalibrating semi-empirical parameters with force matching

Convergence criteria for single-step free-energy calculations: the relation between the Π bias measure and the sample variance

Free energy calculations play a crucial role in simulating chemical processes, enzymatic reactions, and drug design. However, assessing the reliability and convergence of these calculations remains a challenge. This study focuses on single-step free-energy calculations using thermodynamic perturbation. It explores how the sample distributions influence the estimated results and evaluates the reliability of various convergence criteria, including Kofke's bias measure Π and the standard deviation of the energy difference ΔU, σ ΔU . The findings reveal that for Gaussian distributions, there is a straightforward relationship between Π and σ ΔU , free energies can be accurately approximated using a second-order cumulant expansion, and reliable results are attainable for σ ΔU up to 25 kcal mol-1. However, interpreting non-Gaussian distributions is more complex. If the distribution is skewed towards more positive values than a Gaussian, converging the free energy becomes easier, rendering standard convergence criteria overly stringent. Conversely, distributions that are skewed towards more negative values than a Gaussian present greater challenges in achieving convergence, making standard criteria unreliable. We propose a practical approach to assess the convergence of estimated free energies.

Free energies at QM accuracy from force fields via multimap targeted estimation

Accurate predictions of ligand binding affinities would greatly accelerate the first stages of drug discovery campaigns. However, using highly accurate interatomic potentials based on quantum mechanics (QM) in free energy methods has been so far largely unfeasible due to their prohibitive computational cost. Here, we present an efficient method to compute QM free energies from simulations using cheap reference potentials, such as force fields (FFs). This task has traditionally been out of reach due to the slow convergence of computing the correction from the FF to the QM potential. To overcome this bottleneck, we generalize targeted free energy methods to employ multiple maps-implemented with normalizing flow neural networks (NNs)-that maximize the overlap between the distributions. Critically, the method requires neither a separate expensive training phase for the NNs nor samples from the QM potential. We further propose a one-epoch learning policy to efficiently avoid overfitting, and we combine our approach with enhanced sampling strategies to overcome the pervasive problem of poor convergence due to slow degrees of freedom. On the drug-like molecules in the HiPen dataset, the method accelerates the calculation of the free energy difference of switching from an FF to a DFTB3 potential by three orders of magnitude compared to standard free energy perturbation and by a factor of eight compared to previously published nonequilibrium calculations. Our results suggest that our method, in combination with efficient QM/MM calculations, may be used in lead optimization campaigns in drug discovery and to study protein-ligand molecular recognition processes.

Bridging semiempirical and ab initio QM/MM potentials by Gaussian process regression and its sparse variants for free energy simulation

Free energy simulations that employ combined quantum mechanical and molecular mechanical (QM/MM) potentials at ab initio QM (AI) levels are computationally highly demanding. Here, we present a machine-learning-facilitated approach for obtaining AI/MM-quality free energy profiles at the cost of efficient semiempirical QM/MM (SE/MM) methods. Specifically, we use Gaussian process regression (GPR) to learn the potential energy corrections needed for an SE/MM level to match an AI/MM target along the minimum free energy path (MFEP). Force modification using gradients of the GPR potential allows us to improve configurational sampling and update the MFEP. To adaptively train our model, we further employ the sparse variational GP (SVGP) and streaming sparse GPR (SSGPR) methods, which efficiently incorporate previous sample information without significantly increasing the training data size. We applied the QM-(SS)GPR/MM method to the solution-phase SN2 Menshutkin reaction, NH3+CH3Cl→CH3NH3++Cl-, using AM1/MM and B3LYP/6-31+G(d,p)/MM as the base and target levels, respectively. For 4000 configurations sampled along the MFEP, the iteratively optimized AM1-SSGPR-4/MM model reduces the energy error in AM1/MM from 18.2 to 4.4 kcal/mol. Although not explicitly fitting forces, our method also reduces the key internal force errors from 25.5 to 11.1 kcal/mol/Å and from 30.2 to 10.3 kcal/mol/Å for the N-C and C-Cl bonds, respectively. Compared to the uncorrected simulations, the AM1-SSGPR-4/MM method lowers the predicted free energy barrier from 28.7 to 11.7 kcal/mol and decreases the reaction free energy from -12.4 to -41.9 kcal/mol, bringing these results into closer agreement with their AI/MM and experimental benchmarks.

Perspective: Reference-Potential Methods for the Study of Thermodynamic Properties in Chemical Processes: Theory, Applications, and Pitfalls

In silico investigations of enzymatic reactions and chemical reactions in condensed phases often suffer from formidable computational costs due to a large number of degrees of freedom and enormous important volume in phase space. Usually, accuracy must be compromised to trade for efficiency by lowering the reliability of the Hamiltonians employed or reducing the sampling time. Reference-potential methods (RPMs) offer an alternative approach to reaching high accuracy of simulation without much loss of efficiency. In this Perspective, we summarize the idea of RPMs and showcase some recent applications. Most importantly, the pitfalls of these methods are also discussed, and remedies to these pitfalls are presented.

Exploring Routes to Enhance the Calculation of Free Energy Differences via Non-Equilibrium Work SQM/MM Switching Simulations Using Hybrid Charge Intermediates between MM and SQM Levels of Theory or Non-Linear Switching Schemes

Non-equilibrium work switching simulations and Jarzynski's equation are a reliable method for computing free energy differences, ΔAlow→high, between two levels of theory, such as a pure molecular mechanical (MM) and a quantum mechanical/molecular mechanical (QM/MM) description of a system of interest. Despite the inherent parallelism, the computational cost of this approach can quickly become very high. This is particularly true for systems where the core region, the part of the system to be described at different levels of theory, is embedded in an environment such as explicit solvent water. We find that even for relatively simple solute-water systems, switching lengths of at least 5 ps are necessary to compute ΔAlow→high reliably. In this study, we investigate two approaches towards an affordable protocol, with an emphasis on keeping the switching length well below 5 ps. Inserting a hybrid charge intermediate state with modified partial charges, which resembles the charge distribution of the desired high level, makes it possible to obtain reliable calculations with 2 ps switches. Attempts using step-wise linear switching paths, on the other hand, did not lead to improvement, i.e., a faster convergence for all systems. To understand these findings, we analyzed the solutes' properties as a function of the partial charges used and the number of water molecules in direct contact with the solute, and studied the time needed for water molecules to reorient themselves upon a change in the solute's charge distribution.

Accelerated Quantum Mechanics/Molecular Mechanics Simulations via Neural Networks Incorporated with Mechanical Embedding Scheme

A powerful tool to study the mechanism of reactions in solutions or enzymes is to perform the ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations. However, the computational cost is too high due to the explicit electronic structure calculations at every time step of the simulation. A neural network (NN) method can accelerate the QM/MM-MD simulations, but it has long been a problem to accurately describe the QM/MM electrostatic coupling by NN in the electrostatic embedding (EE) scheme. In this work, we developed a new method to accelerate QM/MM calculations in the mechanic embedding (ME) scheme. The potentials and partial point charges of QM atoms are first learned in vacuo by the embedded atom neural networks (EANN) approach. MD simulations are then performed on this EANN/MM potential energy surface (PES) to obtain free energy (FE) profiles for reactions, in which the QM/MM electrostatic coupling is treated in the mechanic embedding (ME) scheme. Finally, a weighted thermodynamic perturbation (wTP) corrects the FE profiles in the ME scheme to the EE scheme. For two reactions in water and one in methanol, our simulations reproduced the B3LYP/MM free energy profiles within 0.5 kcal/mol with a speed-up of 30-60-fold. The results show that the strategy of combining EANN potential in the ME scheme with the wTP correction is efficient and reliable for chemical reaction simulations in liquid. Another advantage of our method is that the QM PES is independent of the MM subsystem, so it can be applied to various MM environments as demonstrated by an S_N2 reaction studied in water and methanol individually, which used the same EANN PES. The free energy profiles are in excellent accordance with the results obtained from B3LYP/MM-MD simulations. In future, this method will be applied to the reactions of enzymes and their variants.

Multireference Generalization of the Weighted Thermodynamic Perturbation Method

We describe the generalized weighted thermodynamic perturbation (gwTP) method for estimating the free energy surface of an expensive "high-level" potential energy function from the umbrella sampling performed with multiple inexpensive "low-level" reference potentials. The gwTP method is a generalization of the weighted thermodynamic perturbation (wTP) method developed by Li and co-workers [J. Chem. Theory Comput. 2018, 14, 5583-5596] that uses a single "low-level" reference potential. The gwTP method offers new possibilities in model design whereby the sampling generated from several low-level potentials may be combined (e.g., specific reaction parameter models that might have variable accuracy at different stages of a multistep reaction). The gwTP method is especially well suited for use with machine learning potentials (MLPs) that are trained against computationally expensive ab initio quantum mechanical/molecular mechanical (QM/MM) energies and forces using active learning procedures that naturally produce multiple distinct neural network potentials. Simulations can be performed with greater sampling using the fast MLPs and then corrected to the ab initio level using gwTP. The capabilities of the gwTP method are demonstrated by creating reference potentials based on the MNDO/d and DFTB2/MIO semiempirical models supplemented with the "range-corrected deep potential" (DPRc). The DPRc parameters are trained to ab initio QM/MM data, and the potentials are used to calculate the free energy surface of stepwise mechanisms for nonenzymatic RNA 2'-O-transesterification model reactions. The extended sampling made possible by the reference potentials allows one to identify unequilibrated portions of the simulations that are not always evident from the short time scale commonly used with ab initio QM/MM potentials. We show that the reference potential approach can yield more accurate ab initio free energy predictions than the wTP method or what can be reasonably afforded from explicit ab initio QM/MM sampling.

Predicting Solvent Effects on SN2 Reaction Rates: Comparison of QM/MM, Implicit, and MM Explicit Solvent Models

Solvents are one of the key variables in the optimization of a synthesis yield or properties of a synthesis product. In this paper, contemporary solvent models are applied to predict the rates of S_N2 reactions in a range of aqueous and non-aqueous solvents. High-level CCSD(T)/CBS//M06-2X/6-31+G(d) gas phase energies were combined with solvation free energies from SMD, SM12, and ADF-COSMO-RS continuum solvent models, as well as molecular mechanics (MM) explicit solvent models with different atomic charge schemes to predict the rate constants of three S_N2 reactions in eight protic and aprotic solvents. It is revealed that the prediction of rate constants in organic solvents is not necessarily less challenging than in water and popular solvent models struggle to predict their rate constants to within 3 log units of experimental values. Among the continuum solvent models, the ADF-COSMO-RS model performed the best in predicting absolute rate contants while the SM12 model was best at predicting relative rate constants with an average accuracy of about 1.5 and 0.8 log units, respectively. The use of computationally more demanding MM explicit solvent models did not translate to improvements in absolute rate constants but was quite effective at predicting relative rate constants due to systematic error cancellation. Free energy barriers obtained from umbrella sampling with explicit solvent QM/MM simulations led to excellent agreement with experimental values, provided that a validated level of theory is used to treat the QM region.

Facilitating ab initio QM/MM free energy simulations by Gaussian process regression with derivative observations

In combined quantum mechanical and molecular mechanical (QM/MM) free energy simulations, how to synthesize the accuracy of ab initio (AI) methods with the speed of semiempirical (SE) methods for a cost-effective QM treatment remains a long-standing challenge. In this work, we present a machine-learning-facilitated method for obtaining AI/MM-quality free energy profiles through efficient SE/MM simulations. In particular, we use Gaussian process regression (GPR) to learn the energy and force corrections needed for SE/MM to match with AI/MM results during molecular dynamics simulations. Force matching is enabled in our model by including energy derivatives into the observational targets through the extended-kernel formalism. We demonstrate the effectiveness of this method on the solution-phase SN2 Menshutkin reaction using AM1/MM and B3LYP/6-31+G(d,p)/MM as the base and target levels, respectively. Trained on only 80 configurations sampled along the minimum free energy path (MFEP), the resulting GPR model reduces the average energy error in AM1/MM from 18.2 to 5.8 kcal mol-1 for the 4000-sample testing set with the average force error on the QM atoms decreased from 14.6 to 3.7 kcal mol-1 Å-1. Free energy sampling with the GPR corrections applied (AM1-GPR/MM) produces a free energy barrier of 14.4 kcal mol-1 and a reaction free energy of -34.1 kcal mol-1, in closer agreement with the AI/MM benchmarks and experimental results.

Quantitative Analysis the Weak Non-Covalent Interactions of the Polymorphs of Donepezil

Donepezil has polymorphism. Different crystalline forms can exhibit different physicochemical properties and biological activities. Exploration of intermolecular interactions is essential to reveal the formation mechanism and differences in properties of polymorphs. This study explores the weak non-covalent intermolecular interactions of donepezil polymorphs through fully ab initio quantum mechanical methods, semi-empirical methods, and Hirshfeld surface analysis. The results show that the Hirshfeld surface analysis method can clearly and intuitively reveal the intermolecular interactions. Theoretical calculations using the atom-atom Coulomb-London-Pauli (AA-CLP) method were also performed to understand the interaction energies toward the total lattice energy. The value of the lattice energy was in accordance with the melting points of the donepezil polymorphs and brought to light the nature of thermal stability. In the specific energy distribution, the contribution of the dispersion force is the most prominent. Further interaction energy analysis found that within a distance of 3.8 Å from the center of the donepezil molecule, different crystalline forms of donepezil molecules have different interaction energies with surrounding molecules. The different interaction energies between polymorphs may lead to polymorphs with different physical-chemical properties.

Improving the Accuracy of Quantum Mechanics/Molecular Mechanics (QM/MM) Models with Polarized Fragment Charges

This paper introduces an economical approach for improving the accuracy and convergence of quantum mechanics/molecular mechanics (QM/MM) models. The approach is tested on a series of neutral and charged amino acids embedded in a 160-water cluster, where their intramolecular proton transfer energies (neutral amino acid → zwitterionic amino acid) were previously obtained at the ωB97X-D/6-31G(d) level of theory. When the charges on the MM atoms were replaced with those obtained at the same QM level of theory used to treat the QM atoms, this significantly improved the accuracy and convergence of the QM/MM models. In particular, the QM/MM model converged to within 1.4 kcal mol^-1 of directly calculated DFT energies for smaller (by as many as 20 waters) QM regions. The use of atomic charges obtained from the natural population analysis yielded the most significant improvement, while other charge schemes such as Mulliken, electrostatic potential, or CM5 led to poorer outcomes. It is further demonstrated that the QM atomic charges can be accurately estimated in a highly efficient manner using an iterative fragmentation approach based on the moving-domain QM/MM method. Similar observations were made when the approach was used to predict the barrier of an S_N2 reaction. Thus, the use of QM-quality atomic charges on MM atoms represents a simple and easy-to-implement strategy for improving the accuracy of QM/MM models.

Accelerating Ab Initio Quantum Mechanical and Molecular Mechanical (QM/MM) Molecular Dynamics Simulations with Multiple Time Step Integration and a Recalibrated Semiempirical QM/MM Hamiltonian

Molecular dynamics (MD) simulations employing ab initio quantum mechanical and molecular mechanical (ai-QM/MM) potentials are considered to be the state of the art, but the high computational cost associated with the ai-QM calculations remains a theoretical challenge for their routine application. Here, we present a modified protocol of the multiple time step (MTS) method for accelerating ai-QM/MM MD simulations of condensed-phase reactions. Within a previous MTS protocol [Nam J. Chem. Theory Comput. 2014, 10, 4175], reference forces are evaluated using a low-level (semiempirical QM/MM) Hamiltonian and employed at inner time steps to propagate the nuclear motions. Correction forces, which arise from the force differences between high-level (ai-QM/MM) and low-level Hamiltonians, are applied at outer time steps, where the MTS algorithm allows the time-reversible integration of the correction forces. To increase the outer step size, which is bound by the highest-frequency component in the correction forces, the semiempirical QM Hamiltonian is recalibrated in this work to minimize the magnitude of the correction forces. The remaining high-frequency modes, which are mainly bond stretches involving hydrogen atoms, are then removed from the correction forces. When combined with a Langevin or SIN(R) thermostat, the modified MTS-QM/MM scheme remains robust with an up to 8 (with Langevin) or 10 fs (with SIN(R)) outer time step (with 1 fs inner time steps) for the chorismate mutase system. This leads to an over 5-fold speedup over standard ai-QM/MM simulations, without sacrificing the accuracy in the predicted free energy profile of the reaction.

Solvent Effects on the Menshutkin Reaction

The Menshutkin reaction is a methyl transfer reaction relevant in fields ranging from biochemistry to chemical synthesis. In the present work, the energetics and solvent distributions for NH₃+MeCl and Pyr+MeBr reactions were investigated in explicit solvent (water, methanol, acetonitrile, benzene, cyclohexane) by means of reactive molecular dynamics simulations. For polar solvents (water, methanol, and acetonitrile) and benzene, strong to moderate catalytic effects for both reactions were found, whereas apolar and bulky cyclohexane interacts weakly with the solute and does not show pronounced barrier reduction. The calculated barrier heights for the Pyr+MeBr reaction in acetonitrile and cyclohexane are 23.2 and 28.1 kcal/mol compared with experimentally measured barriers of 22.5 and 27.6 kcal/mol, respectively. The solvent distributions change considerably between reactant and TS but comparatively little between TS and product conformations of the solute. As the system approaches the transition state, correlated solvent motions occur which destabilize the solvent-solvent interactions. This is required for the system to surmount the barrier. Finally, it is found that the average solvent-solvent interaction energies in the reactant, TS, and product state geometries are correlated with changes in the solvent structure around the solute.

Affordable Ab Initio Path Integral for Thermodynamic Properties via Molecular Dynamics Simulations Using Semiempirical Reference Potential

Path integral molecular dynamics (PIMD) is becoming a routinely applied method for incorporating the nuclear quantum effect in computer simulations. However, direct PIMD simulations at an ab initio level of theory are formidably expensive. Using the protonated 1,8-bis(dimethylamino)naphthalene molecule as an example, we show in this work that the computational expense for the intramolecular proton transfer between the two nitrogen atoms can be remarkably reduced by implementing the idea of reference-potential methods. The simulation time can be easily extended to a scale of nanoseconds while maintaining the accuracy on an ab initio level of theory for thermodynamic properties. In addition, postprocessing can be carried out in parallel on massive computer nodes. A 545-fold reduction in the total CPU time can be achieved in this way as compared to a direct PIMD simulation at the same ab initio level of theory.

Doubly Polarized QM/MM with Machine Learning Chaperone Polarizability

A major shortcoming of semiempirical (SE) molecular orbital methods is their severe underestimation of molecular polarizability compared with experimental and ab initio (AI) benchmark data. In a combined quantum mechanical and molecular mechanical (QM/MM) treatment of solution-phase reactions, solute described by SE methods therefore tends to generate inadequate electronic polarization response to solvent electric fields, which often leads to large errors in free energy profiles. To address this problem, here we present a hybrid framework that improves the response property of SE/MM methods through high-level molecular-polarizability fitting. Specifically, we place on QM atoms a set of corrective polarizabilities (referred to as chaperone polarizabilities), whose magnitudes are determined from machine learning (ML) to reproduce the condensed-phase AI molecular polarizability along the minimum free energy path. These chaperone polarizabilities are then used in a machinery similar to a polarizable force field calculation to compensate for the missing polarization energy in the conventional SE/MM simulations. Because QM atoms in this treatment host SE wave functions as well as classical polarizabilities, both polarized by MM electric fields, we name this method doubly polarized QM/MM (dp-QM/MM). We demonstrate the new method on the free energy simulations of the Menshutkin reaction in water. Using AM1/MM as a base method, we show that ML chaperones greatly reduce the error in the solute molecular polarizability from 6.78 to 0.03 ų with respect to the density functional theory benchmark. The chaperone correction leads to ∼10 kcal/mol of additional polarization energy in the product region, bringing the simulated free energy profiles to closer agreement with the experimental results. Furthermore, the solute-solvent radial distribution functions show that the chaperone polarizabilities modify the free energy profiles through enhanced solvation corrections when the system evolves from the charge-neutral reactant state to the charge-separated transition and product states. These results suggest that the dp-QM/MM method, enabled by ML chaperone polarizabilities, provides a very physical remedy for the underpolarization problem in SE/MM-based free energy simulations.

Targeted Free Energy Perturbation Revisited: Accurate Free Energies from Mapped Reference Potentials

We present an approach that extends the theory of targeted free energy perturbation (TFEP) to calculate free energy differences and free energy surfaces at an accurate quantum mechanical level of theory from a cheaper reference potential. The convergence is accelerated by a mapping function that increases the overlap between the target and the reference distributions. Building on recent work, we show that this map can be learned with a normalizing flow neural network, without requiring simulations with the expensive target potential but only a small number of single-point calculations, and, crucially, avoiding the systematic error that was found previously. We validate the method by numerically evaluating the free energy difference in a system with a double-well potential and by describing the free energy landscape of a simple chemical reaction in the gas phase.

Machine-Learning-Assisted Free Energy Simulation of Solution-Phase and Enzyme Reactions

Despite recent advances in the development of machine learning potentials (MLPs) for biomolecular simulations, there has been limited effort on developing stable and accurate MLPs for enzymatic reactions. Here we report a protocol for performing machine-learning-assisted free energy simulation of solution-phase and enzyme reactions at the ab initio quantum-mechanical/molecular-mechanical (ai-QM/MM) level of accuracy. Within our protocol, the MLP is built to reproduce the ai-QM/MM energy and forces on both QM (reactive) and MM (solvent/enzyme) atoms. As an alternative strategy, a delta machine learning potential (ΔMLP) is trained to reproduce the differences between the ai-QM/MM and semiempirical (se) QM/MM energies and forces. To account for the effect of the condensed-phase environment in both MLP and ΔMLP, the DeePMD representation of a molecular system is extended to incorporate the external electrostatic potential and field on each QM atom. Using the Menshutkin and chorismate mutase reactions as examples, we show that the developed MLP and ΔMLP reproduce the ai-QM/MM energy and forces with errors that on average are less than 1.0 kcal/mol and 1.0 kcal mol-1 Å-1, respectively, for representative configurations along the reaction pathway. For both reactions, MLP/ΔMLP-based simulations yielded free energy profiles that differed by less than 1.0 kcal/mol from the reference ai-QM/MM results at only a fraction of the computational cost.

Free Energy Computation for an Isomerizing Chromophore in a Molecular Cavity via the Average Solvent Electrostatic Configuration Model: Applications in Rhodopsin and Rhodopsin-Mimicking Systems

We present a novel technique for computing the free energy differences between two chromophore "isomers" hosted in a molecular environment (a generalized solvent). Such an environment may range from a relatively rigid protein cavity to a flexible solvent environment. The technique is characterized by the application of the previously reported "average electrostatic solvent configuration" method, and it is based on the idea of using the free energy perturbation theory along with a chromophore annihilation procedure in thermodynamic cycle calculations. The method is benchmarked by computing the ground-state room-temperature relative stabilities between (i) the cis and trans isomers of prototypal animal and microbial rhodopsins and (ii) the analogue isomers of a rhodopsin-like light-driven molecular switch in methanol. Furthermore, we show that the same technology can be used to estimate the activation free energy for the thermal isomerization of systems i-ii by replacing one isomer with a transition state. The results show that the computed relative stability and isomerization barrier magnitudes for the selected systems are in line with the available experimental observation in spite of their widely diverse complexity.

Reaction Path-Force Matching in Collective Variables: Determining Ab Initio QM/MM Free Energy Profiles by Fitting Mean Force

First-principles determination of free energy profiles for condensed-phase chemical reactions is hampered by the daunting costs associated with configurational sampling on ab initio quantum mechanical/molecular mechanical (AI/MM) potential energy surfaces. Here, we report a new method that enables efficient AI/MM free energy simulations through mean force fitting. In this method, a free energy path in collective variables (CVs) is first determined on an efficient reactive aiding potential. Based on the configurations sampled along the free energy path, correcting forces to reproduce the AI/MM forces on the CVs are determined through force matching. The AI/MM free energy profile is then predicted from simulations on the aiding potential in conjunction with the correcting forces. Such cycles of correction-prediction are repeated until convergence is established. As the instantaneous forces on the CVs sampled in equilibrium ensembles along the free energy path are fitted, this procedure faithfully restores the target free energy profile by reproducing the free energy mean forces. Due to its close connection with the reaction path-force matching (RP-FM) framework recently introduced by us, we designate the new method as RP-FM in collective variables (RP-FM-CV). We demonstrate the effectiveness of this method on a type-II solution-phase SN2 reaction, NH3 + CH3Cl (the Menshutkin reaction), simulated with an explicit water solvent. To obtain the AI/MM free energy profiles, we employed the semiempirical AM1/MM Hamiltonian as the base level for determining the string minimum free energy pathway, along which the free energy mean forces are fitted to various target AI/MM levels using the Hartree-Fock (HF) theory, density functional theory (DFT), and the second-order Møller-Plesset perturbation (MP2) theory as the AI method. The forces on the bond-breaking and bond-forming CVs at both the base and target levels are obtained by force transformation from Cartesian to redundant internal coordinates under the Wilson B-matrix formalism, where the linearized FM is facilitated by the use of spline functions. For the Menshutkin reaction tested, our FM treatment greatly reduces the deviations on the CV forces, originally in the range of 12-33 to ∼2 kcal/mol/Å. Comparisons with the experimental and benchmark AI/MM results, tests of the new method under a variety of simulation protocols, and analyses of the solute-solvent radial distribution functions suggest that RP-FM-CV can be used as an efficient, accurate, and robust method for simulating solution-phase chemical reactions.

The SN2 reaction and its relationship with the Walden inversion, the Finkelstein and Menshutkin reactions together with theoretical calculations for the Finkelstein reaction

This communication gives an overview of the relationships between four reactions that although related were not always perceived as such: SN2, Walden, Finkelstein, and Menshutkin. Binary interactions (SN2 & Walden, SN2 & Menshutkin, SN2 & Finkelstein, Walden & Menshutkin, Walden & Finkelstein, Menshutkin & Finkelstein) were reported. Carbon, silicon, nitrogen, and phosphorus as central atoms and fluorides, chlorides, bromides, and iodides as lateral atoms were considered. Theoretical calculations provide Gibbs free energies that were analyzed with linear models to obtain the halide contributions. The M06-2x DFT computational method and the 6-311++G(d,p) basis set have been used for all atoms except for iodine where the effective core potential def2-TZVP basis set was used. Concerning the central atom pairs, carbon/silicon vs. nitrogen/phosphorus, we reported here for the first time that the effect of valence expansion was known for Si but not for P. Concerning the lateral halogen atoms, some empirical models including the interaction between F and I as entering and leaving groups explain the Gibbs free energies.

Accelerated Computation of Free Energy Profile at Ab Initio Quantum Mechanical/Molecular Mechanics Accuracy via a Semiempirical Reference Potential. 4. Adaptive QM/MM

Although quantum mechanical/molecular mechanics (QM/MM) methods are now routinely applied to the studies of chemical reactions in condensed phases and enzymatic reactions, they may experience technical difficulties when the reactive region is varying over time. For instance, when the solvent molecules are directly participating in the reaction, the exchange of water molecules between the QM and MM regions may occur on a time scale comparable to the reaction time. To cope with this situation, several adaptive QM/MM schemes have been proposed. However, these methods either add significantly to the computational cost or introduce artificial restraints to the system. In this work, we developed a novel adaptive QM/MM scheme and applied it to the study of a nucleophilic addition reaction. In this scheme, the configuration sampling was performed with a small QM region (without solvent molecules), and the thermodynamic properties under another potential energy function with a larger QM region (with a certain number of solvent molecules and/or different levels of QM theory) are computed via extrapolation using the reference-potential method. Our simulation results show that this adaptive QM/MM scheme is numerically stable, at least for the case studied in this work. Furthermore, this method also offers an inexpensive way to examine the convergence of the QM/MM calculation with respect to the size of the QM region.

Solvation Induction of Free Energy Barriers of Decarboxylation Reactions in Aqueous Solution from Dual-Level QM/MM Simulations

Carbon dioxide capture, corresponding to the recombination process of decarboxylation reactions of organic acids, is typically barrierless in the gas phase and has a relatively low barrier in aprotic solvents. However, these processes often encounter significant solvent-reorganization-induced barriers in aqueous solution if the decarboxylation product is not immediately protonated. Both the intrinsic stereoelectronic effects and solute-solvent interactions play critical roles in determining the overall decarboxylation equilibrium and free energy barrier. An understanding of the interplay of these factors is important for designing novel materials applied to greenhouse gas capture and storage as well as for unraveling the catalytic mechanisms of a range of carboxy lyases in biological CO₂ production. A range of decarboxylation reactions of organic acids with rates spanning nearly 30 orders of magnitude have been examined through dual-level combined quantum mechanical and molecular mechanical simulations to help elucidate the origin of solvation-induced free energy barriers for decarboxylation and the reverse carboxylation reactions in water.

Computational investigations of allostery in aromatic amino acid biosynthetic enzymes

Allostery, in which binding of ligands to remote sites causes a functional change in the active sites, is a fascinating phenomenon observed in enzymes. Allostery can occur either with or without significant conformational changes in the enzymes, and the molecular basis of its mechanism can be difficult to decipher using only experimental techniques. Computational tools for analyzing enzyme sequences, structures, and dynamics can provide insights into the allosteric mechanism at the atomic level. Combining computational and experimental methods offers a powerful strategy for the study of enzyme allostery. The aromatic amino acid biosynthesis pathway is essential in microorganisms and plants. Multiple enzymes involved in this pathway are sensitive to feedback regulation by pathway end products and are known to use allostery to control their activities. To date, four enzymes in the aromatic amino acid biosynthesis pathway have been computationally investigated for their allosteric mechanisms, including 3-deoxy-d-arabino-heptulosonate 7-phosphate synthase, anthranilate synthase, chorismate mutase, and tryptophan synthase. Here we review the computational studies and findings on the allosteric mechanisms of these four enzymes. Results from these studies demonstrate the capability of computational tools and encourage future computational investigations of allostery in other enzymes of this pathway.

A simplified charge projection scheme for long-range electrostatics in ab initio QM/MM calculations

In a previous work [Pan et al., Molecules 23, 2500 (2018)], a charge projection scheme was reported, where outer molecular mechanical (MM) charges [>10 Å from the quantum mechanical (QM) region] were projected onto the electrostatic potential (ESP) grid of the QM region to accurately and efficiently capture long-range electrostatics in ab initio QM/MM calculations. Here, a further simplification to the model is proposed, where the outer MM charges are projected onto inner MM atom positions (instead of ESP grid positions). This enables a representation of the long-range MM electrostatic potential via augmentary charges (AC) on inner MM atoms. Combined with the long-range electrostatic correction function from Cisneros et al. [J. Chem. Phys. 143, 044103 (2015)] to smoothly switch between inner and outer MM regions, this new QM/MM-AC electrostatic model yields accurate and continuous ab initio QM/MM electrostatic energies with a 10 Å cutoff between inner and outer MM regions. This model enables efficient QM/MM cluster calculations with a large number of MM atoms as well as QM/MM calculations with periodic boundary conditions.

Accelerated Computation of Free Energy Profile at Ab Initio Quantum Mechanical/Molecular Mechanics Accuracy via a Semiempirical Reference Potential. 3. Gaussian Smoothing on Density-of-States

Calculations of the free energy profile, also known as potential of mean force (PMF), along a chosen collective variable (CV) are now routinely applied in the studies of chemical processes, such as enzymatic reactions and chemical reactions in condensed phases. However, if the ab initio quantum mechanical/molecular mechanics (QM/MM) level of accuracy is required for the PMF, it can be formidably demanding even with the most advanced enhanced sampling methods, such as umbrella sampling. To ameliorate this difficulty, we developed a novel method for the computation of the free energy profile based on the reference-potential method recently, in which a low-level reference Hamiltonian is employed for phase space sampling and the free energy profile can be corrected to the level of interest (the target Hamiltonian) by energy reweighting in a nonparametric way. However, when the reference Hamiltonian is very different from the target Hamiltonian, the calculated ensemble averages, including the PMF, often suffer from numerical instability, which mainly comes from the overestimation of the density-of-states (DoS) in the low-energy region. Stochastic samplings of these low-energy configurations are rare events, and some low-energy conformations may get oversampled in simulations of a finite length. In this work, an assumption of Gaussian distribution is applied to the DoS in each CV bin, and the weight of each configuration is rescaled according to the accumulated DoS. The results show that this smoothing process can remarkably reduce the ruggedness of the PMF and increase the reliability of the reference-potential method.

Predicting partition coefficients of drug-like molecules in the SAMPL6 challenge with Drude polarizable force fields

The water-octanol partition coefficient is an important physicochemical property for small molecule drug design. Here, we report our participation in the SAMPL6 logP prediction challenge with free energy perturbation (FEP) calculations in the water phase and in the 1-octanol phase using Drude polarizable force fields. Root mean square error (RMSE) and mean absolute error (MAE) of our prediction are equal to 1.85 and 1.25 logP units. The errors are not evenly distributed. Out of eleven SAMPL6 solutes, FEP/Drude performed very badly on three molecules (deviations all larger than 2 logP units) but good on the remaining eight (deviations all less than 1 logP unit). We find while FEP converges well within one nanosecond in water, simulations in 1-octanol need much longer simulation time and possibly more independent runs for sampling. We also find out that 1-octanol, albeit being a non-polar solvent, still polarizes solute molecules and forms stable hydrogen bonds with them. At the end, we attempt to reweight FEP trajectories with QM/Drude calculations and discuss possible caveats in our simulation setup.