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

CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

Efficient Generation of Torsional Energy Profiles by Multifidelity Gaussian Processes for Hindered Rotor Corrections

Accurate thermochemistry computations often require proper treatment of torsional modes. The one-dimensional hindered rotor model has proven to be a computationally efficient solution, given a sufficiently accurate potential energy surface. Methods that provide potential energies at various compromises of uncertainty and computational time demand can be optimally combined within a multifidelity treatment. In this study, we demonstrate how multifidelity modeling leads to (1) smooth interpolation along low-fidelity scan points with uncertainty estimates, (2) inclusion of high-fidelity data that change the energetic order of conformations, and (3) predicting best next-point calculations to extend an initial coarse grid. Our diverse application set comprises molecules, clusters, and transition states of alcohols, ethers, and rings. We discuss limitations for cases in which the low-fidelity computation is highly unreliable. Different features of the potential energy curve affect different quantities. To obtain "optimal" fits, we apply strategies ranging from simple minimization of deviations to developing an acquisition function tailored for statistical thermodynamics. Bayesian prediction of best next calculations can save a substantial amount of computation time for one- and multidimensional hindered rotors.

Amber free energy tools: Interoperable software for free energy simulations using generalized quantum mechanical/molecular mechanical and machine learning potentials

We report the development and testing of new integrated cyberinfrastructure for performing free energy simulations with generalized hybrid quantum mechanical/molecular mechanical (QM/MM) and machine learning potentials (MLPs) in Amber. The Sander molecular dynamics program has been extended to leverage fast, density-functional tight-binding models implemented in the DFTB+ and xTB packages, and an interface to the DeePMD-kit software enables the use of MLPs. The software is integrated through application program interfaces that circumvent the need to perform "system calls" and enable the incorporation of long-range Ewald electrostatics into the external software's self-consistent field procedure. The infrastructure provides access to QM/MM models that may serve as the foundation for QM/MM-ΔMLP potentials, which supplement the semiempirical QM/MM model with a MLP correction trained to reproduce ab initio QM/MM energies and forces. Efficient optimization of minimum free energy pathways is enabled through a new surface-accelerated finite-temperature string method implemented in the FE-ToolKit package. Furthermore, we interfaced Sander with the i-PI software by implementing the socket communication protocol used in the i-PI client-server model. The new interface with i-PI allows for the treatment of nuclear quantum effects with semiempirical QM/MM-ΔMLP models. The modular interoperable software is demonstrated on proton transfer reactions in guanine-thymine mispairs in a B-form deoxyribonucleic acid helix. The current work represents a considerable advance in the development of modular software for performing free energy simulations of chemical reactions that are important in a wide range of applications.

Electronic and Nuclear Quantum Effects on Proton Transfer Reactions of Guanine-Thymine (G-T) Mispairs Using Combined Quantum Mechanical/Molecular Mechanical and Machine Learning Potentials

Rare tautomeric forms of nucleobases can lead to Watson-Crick-like (WC-like) mispairs in DNA, but the process of proton transfer is fast and difficult to detect experimentally. NMR studies show evidence for the existence of short-time WC-like guanine-thymine (G-T) mispairs; however, the mechanism of proton transfer and the degree to which nuclear quantum effects play a role are unclear. We use a B-DNA helix exhibiting a wGT mispair as a model system to study tautomerization reactions. We perform ab initio (PBE0/6-31G*) quantum mechanical/molecular mechanical (QM/MM) simulations to examine the free energy surface for tautomerization. We demonstrate that while the ab initio QM/MM simulations are accurate, considerable sampling is required to achieve high precision in the free energy barriers. To address this problem, we develop a QM/MM machine learning potential correction (QM/MM-ΔMLP) that is able to improve the computational efficiency, greatly extend the accessible time scales of the simulations, and enable practical application of path integral molecular dynamics to examine nuclear quantum effects. We find that the inclusion of nuclear quantum effects has only a modest effect on the mechanistic pathway but leads to a considerable lowering of the free energy barrier for the GT*⇌G*T equilibrium. Our results enable a rationalization of observed experimental data and the prediction of populations of rare tautomeric forms of nucleobases and rates of their interconversion in B-DNA.

Physical Prior Mean Function-Driven Gaussian Processes Search for Minimum-Energy Reaction Paths with a Climbing-Image Nudged Elastic Band: A General Method for Gas-Phase, Interfacial, and Bulk-Phase Reactions

The climbing-image nudged elastic band (CI-NEB) method serves as an indispensable tool for computational chemists, offering insight into minimum-energy reaction paths (MEPs) by delineating both transition states (TSs) and intermediate nonstationary structures along reaction coordinates. However, executing CI-NEB calculations for reactions with extensive reaction coordinate spans necessitates a large number of images to ensure a reliable convergence of the MEPs and TS structures, presenting a computationally demanding optimization challenge, even with mildly costly electronic-structure methods. In this study, we advocate for the utilization of physically inspired prior mean function-based Gaussian processes (GPs) to expedite MEP exploration and TS optimization via the CI-NEB method. By incorporating reliable prior physical approximations into potential energy surface (PES) modeling, we demonstrate enhanced efficiency in multidimensional CI-NEB optimization with surrogate-based optimizers. Our physically informed GP approach not only outperforms traditional nonsurrogate-based optimizers in optimization efficiency but also on-the-fly learns the reaction path valley during optimization, culminating in significant advancements. The surrogate PES derived from our optimization exhibits high accuracy compared to true PES references, aligning with our emphasis on leveraging reliable physical priors for robust and efficient posterior mean learning in GPs. Through a systematic benchmark study encompassing various reaction pathways, including gas-phase, bulk-phase, and interfacial/surface reactions, our physical GPs consistently demonstrate superior efficiency and reliability. For instance, they outperform the popular fast inertial relaxation engine optimizer by approximately a factor of 10, showcasing their versatility and efficacy in exploring reaction mechanisms and surface reaction PESs.

Training machine learning potentials for reactive systems: A Colab tutorial on basic models

In the last several years, there has been a surge in the development of machine learning potential (MLP) models for describing molecular systems. We are interested in a particular area of this field - the training of system-specific MLPs for reactive systems - with the goal of using these MLPs to accelerate free energy simulations of chemical and enzyme reactions. To help new members in our labs become familiar with the basic techniques, we have put together a self-guided Colab tutorial (https://cc-ats.github.io/mlp_tutorial/), which we expect to be also useful to other young researchers in the community. Our tutorial begins with the introduction of simple feedforward neural network (FNN) and kernel-based (using Gaussian process regression, GPR) models by fitting the two-dimensional Müller-Brown potential. Subsequently, two simple descriptors are presented for extracting features of molecular systems: symmetry functions (including the ANI variant) and embedding neural networks (such as DeepPot-SE). Lastly, these features will be fed into FNN and GPR models to reproduce the energies and forces for the molecular configurations in a Claisen rearrangement reaction.

Surface-Accelerated String Method for Locating Minimum Free Energy Paths

We present a surface-accelerated string method (SASM) to efficiently optimize low-dimensional reaction pathways from the sampling performed with expensive quantum mechanical/molecular mechanical (QM/MM) Hamiltonians. The SASM accelerates the convergence of the path using the aggregate sampling obtained from the current and previous string iterations, whereas approaches like the string method in collective variables (SMCV) or the modified string method in collective variables (MSMCV) update the path only from the sampling obtained from the current iteration. Furthermore, the SASM decouples the number of images used to perform sampling from the number of synthetic images used to represent the path. The path is optimized on the current best estimate of the free energy surface obtained from all available sampling, and the proposed set of new simulations is not restricted to being located along the optimized path. Instead, the umbrella potential placement is chosen to extend the range of the free energy surface and improve the quality of the free energy estimates near the path. In this manner, the SASM is shown to improve the exploration for a minimum free energy pathway in regions where the free energy surface is relatively flat. Furthermore, it improves the quality of the free energy profile when the string is discretized with too few images. We compare the SASM, SMCV, and MSMCV using 3 QM/MM applications: a ribozyme methyltransferase reaction using 2 reaction coordinates, the 2'-O-transphosphorylation reaction of Hammerhead ribozyme using 3 reaction coordinates, and a tautomeric reaction in B-DNA using 5 reaction coordinates. We show that SASM converges the paths using roughly 3 times less sampling than the SMCV and MSMCV methods. All three algorithms have been implemented in the FE-ToolKit package made freely available.

Perspectives on Computational Enzyme Modeling: From Mechanisms to Design and Drug Development

Understanding enzyme mechanisms is essential for unraveling the complex molecular machinery of life. In this review, we survey the field of computational enzymology, highlighting key principles governing enzyme mechanisms and discussing ongoing challenges and promising advances. Over the years, computer simulations have become indispensable in the study of enzyme mechanisms, with the integration of experimental and computational exploration now established as a holistic approach to gain deep insights into enzymatic catalysis. Numerous studies have demonstrated the power of computer simulations in characterizing reaction pathways, transition states, substrate selectivity, product distribution, and dynamic conformational changes for various enzymes. Nevertheless, significant challenges remain in investigating the mechanisms of complex multistep reactions, large-scale conformational changes, and allosteric regulation. Beyond mechanistic studies, computational enzyme modeling has emerged as an essential tool for computer-aided enzyme design and the rational discovery of covalent drugs for targeted therapies. Overall, enzyme design/engineering and covalent drug development can greatly benefit from our understanding of the detailed mechanisms of enzymes, such as protein dynamics, entropy contributions, and allostery, as revealed by computational studies. Such a convergence of different research approaches is expected to continue, creating synergies in enzyme research. This review, by outlining the ever-expanding field of enzyme research, aims to provide guidance for future research directions and facilitate new developments in this important and evolving field.

Exploring torsional conformer space with physical prior mean function-driven meta-Gaussian processes

We present a novel approach for systematically exploring the conformational space of small molecules with multiple internal torsions. Identifying unique conformers through a systematic conformational search is important for obtaining accurate thermodynamic functions (e.g., free energy), encompassing contributions from the ensemble of all local minima. Traditional geometry optimizers focus on one structure at a time, lacking transferability from the local potential-energy surface (PES) around a specific minimum to optimize other conformers. In this work, we introduce a physics-driven meta-Gaussian processes (meta-GPs) method that not only enables efficient exploration of target PES for locating local minima but, critically, incorporates physical surrogates that can be applied universally across the optimization of all conformers of the same molecule. Meta-GPs construct surrogate PESs based on the optimization history of prior conformers, dynamically selecting the most suitable prior mean function (representing prior knowledge in Bayesian learning) as a function of the optimization progress. We systematically benchmarked the performance of multiple GP variants for brute-force conformational search of amino acids. Our findings highlight the superior performance of meta-GPs in terms of efficiency, comprehensiveness of conformer discovery, and the distribution of conformers compared to conventional non-surrogate optimizers and other non-meta-GPs. Furthermore, we demonstrate that by concurrently optimizing, training GPs on the fly, and learning PESs, meta-GPs exhibit the capacity to generate high-quality PESs in the torsional space without extensive training data. This represents a promising avenue for physics-based transfer learning via meta-GPs with adaptive priors in exploring torsional conformer space.