Efficiency of alchemical free energy simulations. I. A practical comparison of the exponential formula, thermodynamic integration, and Bennett's acceptance ratio method

We investigate the relative efficiency of thermodynamic integration, three variants of the exponential formula, also referred to as thermodynamic perturbation, and Bennett's acceptance ratio method to compute relative and absolute solvation free energy differences. Our primary goal is the development of efficient protocols that are robust in practice. We focus on minimizing the number of unphysical intermediate states (λ-states) required for the computation of accurate and precise free energy differences. Several indicators are presented which help decide when additional λ-states are necessary. In all tests Bennett's acceptance ratio method required the least number of λ-states, closely followed by the "double-wide" variant of the exponential formula. Use of the exponential formula in only strict "forward" or "backward" mode was not found to be competitive. Similarly, the performance of thermodynamic integration in terms of efficiency was rather poor. We show that this is caused by the use of the trapezoidal rule as method of numerical quadrature. A systematic study focusing on the optimization of thermodynamic integration is presented in a companion paper.     

Linear-scaling soft-core scheme for alchemical free energy calculations

Alchemical free energy calculations involving the removal or insertion of atoms into condensed phase systems generally make use of soft-core scaling of nonbonded interactions, designed to circumvent numerical instabilities that arise from weakly interacting "hard" atoms in close proximity. Current methods model soft-core atoms by introducing a nonlinear dependence between the shape of the interaction potential and the strength of the interaction. In this article, we propose a soft-core method that avoids introducing such a nonlinear dependence, through the application of a smooth flattening of the potential energy only in a region that is energetically accessible under normal conditions. We discuss the benefits that this entails and explore a selection of applications, including enhanced methods for the estimation of free energy differences and for the automated optimization of the placement of intermediate states in multistage alchemical calculations.     

Non-Boltzmann sampling and Bennett's acceptance ratio method: how to profit from bending the rules

The exact computation of free energy differences requires adequate sampling of all relevant low energy conformations. Especially in systems with rugged energy surfaces, adequate sampling can only be achieved by biasing the exploration process, thus yielding non-Boltzmann probability distributions. To obtain correct free energy differences from such simulations, it is necessary to account for the effects of the bias in the postproduction analysis. We demonstrate that this can be accomplished quite simply with a slight modification of Bennett's Acceptance Ratio method, referring to this technique as Non-Boltzmann Bennett. We illustrate the method by several examples and show how a creative choice of the biased state(s) used during sampling can also improve the efficiency of free energy simulations.     

Comparison of thermodynamic integration and Bennett acceptance ratio for calculating relative protein‐ligand binding free energies

The performances of Bennett's acceptance ratio method and thermodynamic integration (TI) for the calculation of free energy differences in protein simulations are compared. For the latter, the standard trapezoidal rule, Simpson's rule, and Clenshaw‐Curtis integration are used as numerical integration methods. We evaluate the influence of the number and definition of intermediate states on the precision, accuracy, and efficiency of the free energy calculations. Our results show that non‐equidistantly spaced intermediate states are in some cases beneficial for the TI methods. Using several combinations of softness parameters and the λ power dependence, it is shown that these benefits are strongly dependent on the shape of the integrand. Although TI is more user‐friendly due to its simplicity, it was found that Bennett's acceptance ratio method is the more efficient method. It is also the least dependent on the choice of the intermediate states, making it more robust than TI. © 2013 Wiley Periodicals, Inc.     

Soft-core potentials in thermodynamic integration: comparing one- and two-step transformations

Molecular dynamics-based free energy calculations allow the determination of a variety of thermodynamic quantities from computer simulations of small molecules. Thermodynamic integration (TI) calculations can suffer from instabilities during the creation or annihilation of particles. This "singularity" problem can be addressed with "soft-core" potential functions which keep pairwise interaction energies finite for all configurations and provide smooth free energy curves. "One-step" transformations, in which electrostatic and van der Waals forces are simultaneously modified, can be simpler and less expensive than "two-step" transformations in which these properties are changed in separate calculations. Here, we study solvation free energies for molecules of different hydrophobicity using both models. We provide recommended values for the two parameters α(LJ) and β(C) controlling the behavior of the soft-core Lennard-Jones and Coulomb potentials and compare one- and two-step transformations with regard to their suitability for numerical integration. For many types of transformations, the one-step procedure offers a convenient and accurate approach to free energy estimates.     

Avoiding the van der Waals endpoint problem using serial atomic insertion

In the past analyses of the so-called van der Waals end point problem focused on thermodynamic integration. Here we investigate which of the recommendations, such as the need for soft-core potentials, are still valid when Bennett's acceptance ratio method is used. We show that in combination with Bennett's acceptance ratio method intermediate states characterized by the coupling parameter λ can be replaced by intermediate states in which Lennard-Jones interactions are turned on or off on an "atom by atom" basis. By doing so, there is no necessity to use soft-core potentials. In fact, one can compute free energy differences without dedicated code, making it possible to use any molecular dynamics program to compute alchemical free energy differences. Such an approach, which we illustrate by several examples, makes it possible to exploit the tremendous computational power of the graphics processing unit.     

Reducing the bias and uncertainty of free energy estimates by using regression to fit thermodynamic integration data

This report presents the application of polynomial regression for estimating free energy differences using thermodynamic integration data, i.e., slope of free energy with respect to the switching variable λ. We employ linear regression to construct a polynomial that optimally fits the thermodynamic integration data, and thus reduces the bias and uncertainty of the resulting free energy estimate. Two test systems with analytical solutions were used to verify the accuracy and precision of the approach. Our results suggest that use of regression with high degree of polynomials provides the most accurate free energy difference estimates, but often with slightly larger uncertainty, compared to commonly used quadrature techniques. High degree polynomials possess the flexibility to closely fit the thermodynamic integration data but are often sensitive to small changes in the data points. Thus, we also used Chebyshev nodes to guide in the selection of nonequidistant λ values for use in thermodynamic integration. We conclude that polynomial regression with nonequidistant λ values delivers the most accurate and precise free energy estimates for thermodynamic integration data for the systems considered here. Software and documentation is available at http://www.phys.uidaho.edu/ytreberg/software . © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

On the direct calculation of the free energy of quantization for molecular systems in the condensed phase

Using the path integral formalism or the Feynman-Hibbs approach, various expressions for the free energy of quantization for a molecular system in the condensed phase can be derived. These lead to alternative methods to directly compute quantization free energies from molecular dynamics computer simulations, which were investigated with an eye to their practical use. For a test system of liquid neon, two methods are shown to be most efficient for a direct evaluation of the excess free energy of quantization. One of them makes use of path integral simulations in combination with a single-step free energy perturbation approach and was previously reported in the literature. The other method employs a Feynman-Hibbs effective Hamiltonian together with the thermodynamic integration formalism. However, both methods are found to give less accurate results for the excess free energy of quantization than the estimate obtained from explicit path integral calculations on the excess free energy of the neon liquid in the classical and quantum mechanical limit. Suggestions are made to make both methods more accurate.     

Calculation of sequence‐dependent free energies of hydration of dipeptides formed by alanine and glycine

The relative free energies of hydration of the dipeptides glycylalanine and alanyl‐glycine in their naturally occurring form have been calculated both for the zwitterionic and protonated species. Emphasis was laid on comparisons between the conventional cutoff method and the Particle Mesh Ewald method to account for possible differences in electrostatic contributions to the free energy. Furthermore, the convergence behavior of the total free energy and its individual contributions were examined. The results, obtained by means of the thermodynamic integration technique as implemented in the free energy module of the AMBER program suite, suggest that in aqueous solution glycylalanine is more stable than alanylglycine by 2.7 kcal/mol in the zwitterionic form and by 3.5 kcal/mol in the protonated form. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 846–860, 2001     

Evaluating thermodynamic integration performance of the new amber molecular dynamics package and assess potential halogen bonds of enoyl‐ACP reductase (FabI) benzimidazole inhibitors

Thermodynamic integration (TI) can provide accurate binding free energy insights in a lead optimization program, but its high computational expense has limited its usage. In the effort of developing an efficient and accurate TI protocol for FabI inhibitors lead optimization program, we carefully compared TI with different Amber molecular dynamics (MD) engines (sander and pmemd), MD simulation lengths, the number of intermediate states and transformation steps, and the Lennard‐Jones and Coulomb Softcore potentials parameters in the one‐step TI, using eleven benzimidazole inhibitors in complex with Francisella tularensis enoyl acyl reductase (FtFabI). To our knowledge, this is the first study to extensively test the new AMBER MD engine, pmemd, on TI and compare the parameters of the Softcore potentials in the one‐step TI in a protein‐ligand binding system. The best performing model, the one‐step pmemd TI, using 6 intermediate states and 1 ns MD simulations, provides better agreement with experimental results (RMSD = 0.52 kcal/mol) than the best performing implicit solvent method, QM/MM‐GBSA from our previous study (RMSD = 3.00 kcal/mol), while maintaining similar efficiency. Briefly, we show the optimized TI protocol to be highly accurate and affordable for the FtFabI system. This approach can be implemented in a larger scale benzimidazole scaffold lead optimization against FtFabI. Lastly, the TI results here also provide structure‐activity relationship insights, and suggest the parahalogen in benzimidazole compounds might form a weak halogen bond with FabI, which is a well‐known halogen bond favoring enzyme. © 2015 Wiley Periodicals, Inc.     

Ab initio molecular orbital study of potential energy surface for the H2NO(B1)→NO(Π)+H2 reaction

The mechanism of the H₂NO(²B₁)→NO(²Π)+H₂ reaction has been examined using ab initio molecular orbital methods. Ground-state and first-excited-state potential surfaces were plotted at the FOCI/cc-pVTZ level of theory as functions of two appropriate internal degrees of freedom. A conical intersection was found on the C_s pathway that is symmetric with respect to the plane perpendicular to the molecular plane of C_2v H₂NO(²B₁). It is therefore considered that trajectories that start from H₂NO(²B₁) towards the product region detour around the conical intersection, pass through the neighborhood of the transition state that is located at the saddle point on the C_s pathway, and finally reach the products, NO(²Π)+H₂. Thus we can explain the mechanism of the H₂NO(²B₁)→NO(²Π)+H₂ reaction, which has remained unclear to date.