QM/MM Alchemical Free Energy Simulations: Challenges and Recent Developments

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.

Reweighting from Molecular Mechanics Force Fields to the ANI-2x Neural Network Potential

To achieve chemical accuracy in free energy calculations, it is necessary to accurately describe the system's potential energy surface and efficiently sample configurations from its Boltzmann distribution. While neural network potentials (NNPs) have shown significantly higher accuracy than classical molecular mechanics (MM) force fields, they have a limited range of applicability and are considerably slower than MM potentials, often by orders of magnitude. To address this challenge, Rufa et al. [Rufa et al. bioRxiv 2020, 10.1101/2020.07.29.227959.] suggested a two-stage approach that uses a fast and established MM alchemical energy protocol, followed by reweighting the results using NNPs, known as endstate correction or indirect free energy calculation. This study systematically investigates the accuracy and robustness of reweighting from an MM reference to a neural network target potential (ANI-2x) for an established data set in vacuum, using single-step free-energy perturbation (FEP) and nonequilibrium (NEQ) switching simulation. We assess the influence of longer switching lengths and the impact of slow degrees of freedom on outliers in the work distribution and compare the results to those of multistate equilibrium free energy simulations. Our results demonstrate that free energy calculations between NNPs and MM potentials should be preferably performed using NEQ switching simulations to obtain accurate free energy estimates. NEQ switching simulations between the MM potentials and NNPs are efficient, robust, and trivial to implement.

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.

Optimizing the Calculation of Free Energy Differences in Nonequilibrium Work SQM/MM Switching Simulations

A key step during indirect alchemical free energy simulations using quantum mechanical/molecular mechanical (QM/MM) hybrid potential energy functions is the calculation of the free energy difference ΔAlow→high between the low level (e.g., pure MM) and the high level of theory (QM/MM). A reliable approach uses nonequilibrium work (NEW) switching simulations in combination with Jarzynski's equation; however, it is computationally expensive. In this study, we investigate whether it is more efficient to use more shorter switches or fewer but longer switches. We compare results obtained with various protocols to reference free energy differences calculated with Crooks' equation. The central finding is that fewer longer switches give better converged results. As few as 200 sufficiently long switches lead to ΔAlow→high values in good agreement with the reference results. This optimized protocol reduces the computational cost by a factor of 40 compared to earlier work. We also describe two tools/ways of analyzing the raw data to detect sources of poor convergence. Specifically, we find it helpful to analyze the raw data (work values from the NEW switching simulations) in a quasi-time series-like manner. Principal component analysis helps to detect cases where one or more conformational degrees of freedom are different at the low and high level of theory.

Investigating the molecular mechanism of hydroxylated bromdiphenyl ethers to inhibit the thyroid hormone sulfotransferase SULT1A1

Hydroxylated bromodiphenyl ethers (OH-BDEs) have raised great concern due to their potential endocrine disrupting effects on humans. In vitro experiments have indicated OH-BDEs can inhibit the activity of thyroid hormone (TH) sulfotransferases (SULTs); however, the molecular mechanism has not been investigated in depth. In this work, we employed 17 OH-BDEs with five or fewer Br atoms, and performed integrated computational simulations to unravel the possible inhibition mechanism of OH-BDEs on human SULT1A1. The molecular docking results demonstrate that OH-BDEs form hydrogen bonds with residues Lys106 and His108, and the neutral OH-BDEs show comparable binding energies with their anionic counterparts. The further hybrid quantum mechanical/molecular mechanical (QM/MM) calculations unravel a metabolic mechanism of OH-BDEs comprised by proton abstraction and sulfation steps. This mechanism is involved in the SULT1A1 inhibition by some OH-BDEs comprised of three or fewer Br atoms, while other OH-BDEs likely only form ternary complexes to competitively inhibit SULT1A1 activity. Moreover, the effect of the hydroxyl group of OH-BDEs on SULT1A1 inhibition potential follows the order of ortho-OH BDE > meta-OH BDE > para-OH BDE. These results provide an insight into the inhibition mechanism of OH-BDEs to SULT1A1 at the molecular level, which are beneficial in illuminating the molecular initiating events involved in the TH disruption of OH-BDEs.

Hybrid Solvation Model with First Solvation Shell for Calculation of Solvation Free Energy

We present a hybrid solvation model with first solvation shell to calculate solvation free energies. This hybrid model combines the quantum mechanics and molecular mechanics methods with the analytical expression based on the Born solvation model to calculate solvation free energies. Based on calculated free energies of solvation and reaction profiles in gas phase, we set up a unified scheme to predict reaction profiles in solution. The predicted solvation free energies and reaction barriers are compared with experimental results for twenty bimolecular nucleophilic substitution reactions. These comparisons show that our hybrid solvation model can predict reliable solvation free energies and reaction barriers for chemical reactions of small molecules in aqueous solution.

Prediction of aqueous free energies of solvation using coupled QM and MM explicit solvent simulations

A method based on molecular dynamics simulations which employ two distinct levels of theory is proposed and tested for the prediction of Gibbs free energies of solvation for non-ionic solutes in water. The method consists of two additive contributions: (i) an evaluation of the free energy of solvation predicted by a computationally efficient molecular mechanics (MM) method; and (ii) an evaluation of the free energy difference between the potential energy surface of the MM method and that of a more computationally intensive first-principles quantum-mechanical (QM) method. The latter is computed by a thermodynamic integration method based on a series of shorter molecular dynamics simulations that employ weighted averages of the QM and MM force evaluations. The combined computational approach is tested against the experimental free energies of aqueous solvation for four solutes. For solute-solvent interactions that are found to be described qualitatively well by the MM method, the QM correction makes a modest improvement in the predicted free energy of aqueous solvation. However, for solutes that are found to not be adequately described by the MM method, the QM correction does not improve agreement with experiment. These preliminary results provide valuable insights into the novel concept of implementing thermodynamic integration between two model chemistries, suggesting that it is possible to use QM methods to improve upon the MM predictions of free energies of aqueous solvation.

Use of Interaction Energies in QM/MM Free Energy Simulations

The use of the most accurate (i.e., QM or QM/MM) levels of theory for free energy simulations (FES) is typically not possible. Primarily, this is because the computational cost associated with the extensive configurational sampling needed for converging FES is prohibitive. To ensure the feasibility of QM-based FES, the "indirect" approach is generally taken, necessitating a free energy calculation between the MM and QM/MM potential energy surfaces. Ideally, this step is performed with standard free energy perturbation (Zwanzig's equation) as it only requires simulations be carried out at the low level of theory; however, work from several groups over the past few years has conclusively shown that Zwanzig's equation is ill-suited to this task. As such, many approximations have arisen to mitigate difficulties with Zwanzig's equation. One particularly popular notion is that the convergence of Zwanzig's equation can be improved by using interaction energy differences instead of total energy differences. Although problematic numerical fluctuations (a major problem when using Zwanzig's equation) are indeed reduced, our results and analysis demonstrate that this "interaction energy approximation" (IEA) is theoretically incorrect, and the implicit approximation invoked is spurious at best. Herein, we demonstrate this via solvation free energy calculations using IEA from two different low levels of theory to the same target high level. Results from this proof-of-concept consistently yield the wrong results, deviating by ∼1.5 kcal/mol from the rigorously obtained value.

The Good, the Bad, and the Ugly: "HiPen", a New Dataset for Validating (S)QM/MM Free Energy Simulations

Indirect (S)QM/MM free energy simulations (FES) are vital to efficiently incorporating sufficient sampling and accurate (QM) energetic evaluations when estimating free energies of practical/experimental interest. Connecting between levels of theory, i.e., calculating Δ A l o w → h i g h , remains to be the most challenging step within an indirect FES protocol. To improve calculations of Δ A l o w → h i g h , we must: (1) compare the performance of all FES methods currently available; and (2) compile and maintain datasets of Δ A l o w → h i g h calculated for a wide-variety of molecules so that future practitioners may replicate or improve upon the current state-of-the-art. Towards these two aims, we introduce a new dataset, "HiPen", which tabulates Δ A g a s M M → 3 o b (the free energy associated with switching from an M M to an S C C - D F T B molecular description using the 3ob parameter set in gas phase), calculated for 22 drug-like small molecules. We compare the calculation of this value using free energy perturbation, Bennett's acceptance ratio, Jarzynski's equation, and Crooks' equation. We also predict the reliability of each calculated Δ A g a s M M → 3 o b by evaluating several convergence criteria including sample size hysteresis, overlap statistics, and bias metric ( Π ). Within the total dataset, three distinct categories of molecules emerge: the "good" molecules, for which we can obtain converged Δ A g a s M M → 3 o b using Jarzynski's equation; "bad" molecules which require Crooks' equation to obtain a converged Δ A g a s M M → 3 o b ; and "ugly" molecules for which we cannot obtain reliably converged Δ A g a s M M → 3 o b with either Jarzynski's or Crooks' equations. We discuss, in depth, results from several example molecules in each of these categories and describe how dihedral discrepancies between levels of theory cause convergence failures even for these gas phase free energy simulations.

Accelerating QM/MM Free Energy Computations via Intramolecular Force Matching

The calculation of free energy differences between levels of theory has numerous potential pitfalls. Chief among them is the lack of overlap, i.e., ensembles generated at one level of theory (e.g., "low") not being good approximations of ensembles at the other (e.g., "high"). Numerous strategies have been devised to mitigate this issue. However, the most straightforward approach is to ensure that the "low" level ensemble more closely resembles that of the "high". Ideally, this is done without increasing computational cost. Herein, we demonstrate that by reparametrizing classical intramolecular potentials to reproduce high level forces (i.e., force matching) configurational overlap between a "low" (i.e., classical) and "high" (i.e., quantum) level can be significantly improved. This procedure is validated on two test cases and results in vastly improved convergence of free energy simulations.

A Comparison of QM/MM Simulations with and without the Drude Oscillator Model Based on Hydration Free Energies of Simple Solutes

Maintaining a proper balance between specific intermolecular interactions and non-specific solvent interactions is of critical importance in molecular simulations, especially when predicting binding affinities or reaction rates in the condensed phase. The most rigorous metric for characterizing solvent affinity are solvation free energies, which correspond to a transfer from the gas phase into solution. Due to the drastic change of the electrostatic environment during this process, it is also a stringent test of polarization response in the model. Here, we employ both the CHARMM fixed charge and polarizable force fields to predict hydration free energies of twelve simple solutes. The resulting classical ensembles are then reweighted to obtain QM/MM hydration free energies using a variety of QM methods, including MP2, Hartree⁻Fock, density functional methods (BLYP, B3LYP, M06-2X) and semi-empirical methods (OM2 and AM1 ). Our simulations test the compatibility of quantum-mechanical methods with molecular-mechanical water models and solute Lennard⁻Jones parameters. In all cases, the resulting QM/MM hydration free energies were inferior to purely classical results, with the QM/MM Drude force field predictions being only marginally better than the QM/MM fixed charge results. In addition, the QM/MM results for different quantum methods are highly divergent, with almost inverted trends for polarizable and fixed charge water models. While this does not necessarily imply deficiencies in the QM models themselves, it underscores the need to develop consistent and balanced QM/MM interactions. Both the QM and the MM component of a QM/MM simulation have to match, in order to avoid artifacts due to biased solute⁻solvent interactions. Finally, we discuss strategies to improve the convergence and efficiency of multi-scale free energy simulations by automatically adapting the molecular-mechanics force field to the target quantum method.

On the convergence of multi-scale free energy simulations

In this work we employ simple model systems to evaluate the relative performance of two of the most important free energy methods: The Zwanzig equation (also known as "Free energy perturbation") and Bennett's acceptance ratio method (BAR). Although our examples should be transferable to other kinds of free energy simulations, we focus on applications of multi-scale free energy simulations. Such calculations are especially complex, since they connect two different levels of theory with very different requirements in terms of speed, accuracy, sampling and parallelizability. We try to reconcile all those different factors by developing some simple criteria to guide the early stages of the development of a free energy protocol. This is accomplished by quantifying how many λ intermediate steps and how many potential energy evaluations are necessary in order to reach a certain level of convergence.

Computing converged free energy differences between levels of theory via nonequilibrium work methods: Challenges and opportunities

We demonstrate that Jarzynski's equation can be used to reliably compute free energy differences between low and high level representations of systems. The need for such a calculation arises when employing the so-called "indirect" approach to free energy simulations with mixed quantum mechanical/molecular mechanical (QM/MM) Hamiltonians; a popular technique for circumventing extensive simulations involving quantum chemical computations. We have applied this methodology to several small and medium sized organic molecules, both in the gas phase and explicit solvent. Test cases include several systems for which the standard approach; that is, free energy perturbation between low and high level description, fails to converge. Finally, we identify three major areas in which the difference between low and high level representations make the calculation of ΔAlow→high difficult: bond stretching and angle bending, different preferred conformations, and the response of the MM region to the charge distribution of the QM region. © 2016 Wiley Periodicals, Inc.

The role of intramolecular nonbonded interaction and angle sampling in single-step free energy perturbation

Single-step free energy perturbation (sFEP) has often been proposed as an efficient tool for a quick free energy scan due to its straightforward protocol and the ability to recycle an existing molecular dynamics trajectory for free energy calculations. Although sFEP is expected to fail when the sampling of a system is inefficient, it is often expected to hold for an alchemical transformation between ligands with a moderate difference in their sizes, e.g., transforming a benzene into an ethylbenzene. Yet, exceptions were observed in calculations for anisole and methylaniline, which have similar physical sizes as ethylbenzene. In this study, we show that such exceptions arise from the sampling inefficiency on an unexpected rigid degree of freedom, namely, the bond angle θ. The distributions of θ differ dramatically between two end states of a sFEP calculation, i.e., the conformation of the ligand changes significantly during the alchemical transformation process. Our investigation also reveals the interrelation between the ligand conformation and the intramolecular nonbonded interactions. This knowledge suggests a best combination of the ghost ligand potential and the dual topology setting, which improves the accuracy in a single reference sFEP calculation by bringing down its error from around 5k_BT to k_BT.

Calculating distribution coefficients based on multi-scale free energy simulations: an evaluation of MM and QM/MM explicit solvent simulations of water-cyclohexane transfer in the SAMPL5 challenge

One of the central aspects of biomolecular recognition is the hydrophobic effect, which is experimentally evaluated by measuring the distribution coefficients of compounds between polar and apolar phases. We use our predictions of the distribution coefficients between water and cyclohexane from the SAMPL5 challenge to estimate the hydrophobicity of different explicit solvent simulation techniques. Based on molecular dynamics trajectories with the CHARMM General Force Field, we compare pure molecular mechanics (MM) with quantum-mechanical (QM) calculations based on QM/MM schemes that treat the solvent at the MM level. We perform QM/MM with both density functional theory (BLYP) and semi-empirical methods (OM1, OM2, OM3, PM3). The calculations also serve to test the sensitivity of partition coefficients to solute polarizability as well as the interplay of the quantum-mechanical region with the fixed-charge molecular mechanics environment. Our results indicate that QM/MM with both BLYP and OM2 outperforms pure MM. However, this observation is limited to a subset of cases where convergence of the free energy can be achieved.

Redox potentials and acidity constants from density functional theory based molecular dynamics

CONSPECTUS: All-atom methods treat solute and solvent at the same level of electronic structure theory and statistical mechanics. All-atom computation of acidity constants (pKa) and redox potentials is still a challenge. In this Account, we review such a method combining density functional theory based molecular dynamics (DFTMD) and free energy perturbation (FEP) methods. The key computational tool is a FEP based method for reversible insertion of a proton or electron in a periodic DFTMD model system. The free energy of insertion (work function) is computed by thermodynamic integration of vertical energy gaps obtained from total energy differences. The problem of the loss of a physical reference for ionization energies under periodic boundary conditions is solved by comparing with the proton work function computed for the same supercell. The scheme acts as a computational hydrogen electrode, and the DFTMD redox energies can be directly compared with experimental redox potentials. Consistent with the closed shell nature of acid dissociation, pKa estimates computed using the proton insertion/removal scheme are found to be significantly more accurate than the redox potential calculations. This enables us to separate the DFT error from other sources of uncertainty such as finite system size and sampling errors. Drawing an analogy with charged defects in solids, we trace the error in redox potentials back to underestimation of the energy gap of the extended states of the solvent. Accordingly the improvement in the redox potential as calculated by hybrid functionals is explained as a consequence of the opening up of the bandgap by the Hartree-Fock exchange component in hybrids. Test calculations for a number of small inorganic and organic molecules show that the hybrid functional implementation of our method can reproduce acidity constants with an uncertainty of 1-2 pKa units (0.1 eV). The error for redox potentials is in the order of 0.2 V.

Predicting hydration free energies with a hybrid QM/MM approach: an evaluation of implicit and explicit solvation models in SAMPL4

The correct representation of solute-water interactions is essential for the accurate simulation of most biological phenomena. Several highly accurate quantum methods are available to deal with solvation by using both implicit and explicit solvents. So far, however, most evaluations of those methods were based on a single conformation, which neglects solute entropy. Here, we present the first test of a novel approach to determine hydration free energies that uses molecular mechanics (MM) to sample phase space and quantum mechanics (QM) to evaluate the potential energies. Free energies are determined by using re-weighting with the Non-Boltzmann Bennett (NBB) method. In this context, the method is referred to as QM-NBB. Based on snapshots from MM sampling and accounting for their correct Boltzmann weight, it is possible to obtain hydration free energies that incorporate the effect of solute entropy. We evaluate the performance of several QM implicit solvent models, as well as explicit solvent QM/MM for the blind subset of the SAMPL4 hydration free energy challenge. While classical free energy simulations with molecular dynamics give root mean square deviations (RMSD) of 2.8 and 2.3 kcal/mol, the hybrid approach yields an improved RMSD of 1.6 kcal/mol. By selecting an appropriate functional and basis set, the RMSD can be reduced to 1 kcal/mol for calculations based on a single conformation. Results for a selected set of challenging molecules imply that this RMSD can be further reduced by using NBB to reweight MM trajectories with the SMD implicit solvent model.