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.     

Applying efficient implicit nongeometric constraints in alchemical free energy simulations

Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in molecular dynamics simulations. Advanced variations of alchemical free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ-dynamics simulation method, the conventional λ parameters that are associated with the progress variables in alchemical transformations are treated as dynamic variables and are constrained such that: 0 ≤ λ(i) ≤ 1 and Σ(i = 1)(N) λ(i) = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by the functional form: λ(i) = e(c sinθ(i))/Σ(j = 1)(N) e(c sinθ(j)) that oscillates between λ(i) = 0 and λ(i) = 1 and has relatively steep transitions between these endpoints.     

The solid surface free energy calculation. I. In defense of the multicomponent approach

The acid-base approach to the calculation of solid surface free energy and liquid-liquid interfacial tensions is a practical example of application of correlation analysis, and thus it is an approximate approach. In these limits, and provided that wide and well-obtained sets of contact angles or interfacial tension data are used for their computation, surface tension components can be considered as material properties. Although their numerical value depends on the characteristics of the chosen reference material, their chemical meaning is independent on the selected scale. Contact angles contain accessible information about intermolecular forces; using surface tension component (STC) acid-base theory, one can extract this information only making very careful use of the mathematical apparatus of correlation analysis. The specific mathematical methods used to obtain these results are illustrated by using as an example a base of data obtained by the supporters of the equation-of-state theory (EQS). The achievements are appreciably good and the agreement between STC and EQS is discussed.     

Sampling enhancement for the quantum mechanical potential based molecular dynamics simulations: a general algorithm and its extension for free energy calculation on rugged energy surface

An approach is developed in the replica exchange framework to enhance conformational sampling for the quantum mechanical (QM) potential based molecular dynamics simulations. Importantly, with our enhanced sampling treatment, a decent convergence for electronic structure self-consistent-field calculation is robustly guaranteed, which is made possible in our replica exchange design by avoiding direct structure exchanges between the QM-related replicas and the activated (scaled by low scaling parameters or treated with high "effective temperatures") molecular mechanical (MM) replicas. Although the present approach represents one of the early efforts in the enhanced sampling developments specifically for quantum mechanical potentials, the QM-based simulations treated with the present technique can possess the similar sampling efficiency to the MM based simulations treated with the Hamiltonian replica exchange method (HREM). In the present paper, by combining this sampling method with one of our recent developments (the dual-topology alchemical HREM approach), we also introduce a method for the sampling enhanced QM-based free energy calculations.     

A novel approach to the calculation of the free energy due to molecular flexibility

The contribution of the molecular flexibility to the solvation excess free energy is expressed in terms of probabilities of reaching hard limits on intramolecular coordinates in a series of calculations successively relaxing those limits. Numerical tests on the harmonic oscillator are also presented and used to make suggestion about computational issues.     

On the use of orientational restraints and symmetry corrections in alchemical free energy calculations

Alchemical free energy calculations are becoming a useful tool for calculating absolute binding free energies of small molecule ligands to proteins. Here, we find that the presence of multiple metastable ligand orientations can cause convergence problems when distance restraints alone are used. We demonstrate that the use of orientational restraints can greatly accelerate the convergence of these calculations. However, even with this acceleration, we find that sufficient sampling requires substantially longer simulations than are used in many published protocols. To further accelerate convergence, we introduce a new method of configuration space decomposition by orientation which reduces required simulation lengths by at least a factor of 5 in the cases examined. Our method is easily parallelizable, well suited for cases where a ligand cocrystal structure is not available, and can utilize initial orientations generated by docking packages.     

Efficiency of alchemical free energy simulations. II. Improvements for thermodynamic integration

We attempt to optimize the efficiency of thermodynamic integration, as defined by the minimal number of unphysical intermediate states required for the computation of accurate and precise free energy differences. The suitability of various numerical quadrature methods is tested. In particular, we compare the trapezoidal rule, Simpson's rule, Gauss-Legendre, Gauss-Kronrod-Patterson, and Clenshaw-Curtis integration, as well as integration based on a cubic spline approximation of the integrand. We find that Simpson's rule and spline integration are already significantly more efficient that the trapezoidal rule, i.e., correct free energy differences can be obtained using fewer λ-states. We demonstrate that Simpson's rule can be used advantageously with nonequidistant values of the abscissa, which increases the flexibility of the method. Efficiency is enhanced even further if higher order methods, such as Gauss-Legendre, Gauss-Kronrod-Patterson, or Clenshaw-Curtis integration, are used; no more than seven λ-states, which in the case of Clenshaw-Curtis integration include the physical end states, were required for accurate results in all test problems studied. Thus, the performance of thermodynamic integration can equal that of Bennett's acceptance ratio method. We also show, however, that the high efficiency found here relies on the particular functional form of the soft-core potential used; overall, thermodynamic integration is more susceptible to the details of the hybrid Hamiltonian used than Bennett's acceptance ratio method. Therefore, we recommend Bennett's acceptance ratio method as the most robust method to compute alchemical free energy differences; nevertheless, scenarios when thermodynamic integration may be preferable are discussed.     

Distribution-function approach to free energy computation

Connections are explored between the free energy difference of two systems and the microscopic distribution functions of the energy difference. On the basis of a rigorous relationship between the energy distribution functions and the free energy, the scheme of error minimization is introduced to derive accurate and simple methods of free energy computation. A set of distribution-function approaches are then examined against model systems, and the newly derived methods exhibit state-of-art performance. It is shown that the notion of error minimization is powerful to improve the free energy calculation using distribution functions.     

The solid surface free energy calculation; II. The limits of the Zisman and of the "equation-of-state" approaches

This paper follows the "defense" of the Good-van Oss-Chaudhury (GvOC) acid-base approach made in Part I and carries out a detailed analysis of the Zisman critical surface energy and, mainly, of the Neumann equation-of-state (EQS) theory. The analysis is made on both a "practical" and a theoretical basis, trying to highlight the acceptable fitting results of axisymmetric drop shape analysis (ADSA) methods and their independence of the assumed thermodynamic foundations of EQS. Some new and original criticisms of the EQS approach are raised and it is shown that other purely semiempirical models, represented by different fitting equations with the same number of parameters, can represent the data measured by ADSA method with the same goodness as EQS. The equation of state appears as one of many semiempirical approaches for the evaluation of surface free energy of solids. Independent of the previous analysis, the criteria used in ADSA measurements are evaluated and some comments made on them.     

Impact of domain knowledge on blinded predictions of binding energies by alchemical free energy calculations

The Drug Design Data Resource (D3R) consortium organises blinded challenges to address the latest advances in computational methods for ligand pose prediction, affinity ranking, and free energy calculations. Within the context of the second D3R Grand Challenge several blinded binding free energies predictions were made for two congeneric series of Farsenoid X Receptor (FXR) inhibitors with a semi-automated alchemical free energy calculation workflow featuring FESetup and SOMD software tools. Reasonable performance was observed in retrospective analyses of literature datasets. Nevertheless, blinded predictions on the full D3R datasets were poor due to difficulties encountered with the ranking of compounds that vary in their net-charge. Performance increased for predictions that were restricted to subsets of compounds carrying the same net-charge. Disclosure of X-ray crystallography derived binding modes maintained or improved the correlation with experiment in a subsequent rounds of predictions. The best performing protocols on D3R set1 and set2 were comparable or superior to predictions made on the basis of analysis of literature structure activity relationships (SAR)s only, and comparable or slightly inferior, to the best submissions from other groups.     

Optimal pairwise and non-pairwise alchemical pathways for free energy calculations of molecular transformation in solution phase

We estimate the global minimum variance path for computing the free energy insertion into or deletion of small molecules from a dense fluid. We perform this optimization over all pair potentials, irrespective of functional form, using functional optimization with a two-body approximation for the radial distribution function. Surprisingly, the optimal pairwise path obtained via this method is almost identical to the path obtained using a optimized generalized "soft core" potential reported by Pham and Shirts [J. Chem. Phys. 135, 034114 (2011)]. We also derive the lowest variance non-pairwise potential path for molecular insertion or deletion and compare its efficiency to the pairwise path. Under certain conditions, non-pairwise pathways can reduce the total variance by up to 60% compared to optimal pairwise pathways. However, optimal non-pairwise pathways do not appear generally feasible for practical free energy calculations because an accurate estimate of the free energy, the parameter that is itself is desired, is required for constructing this non-pairwise path. Additionally, simulations at most intermediate states of these non-pairwise paths have significantly longer correlation times, often exceeding standard simulation lengths for solvation of bulky molecules. The findings suggest that the previously obtained soft core pathway is the lowest variance pathway for molecular insertion or deletion in practice. The findings also demonstrate the utility of functional optimization for determining the efficiency of thermodynamic processes performed with molecular simulation.     

Accurate free energy calculation along optimized paths

The path‐based methods of free energy calculation, such as thermodynamic integration and free energy perturbation, are simple in theory, but difficult in practice because in most cases smooth paths do not exist, especially for large molecules. In this article, we present a novel method to build the transition path of a peptide. We use harmonic potentials to restrain its nonhydrogen atom dihedrals in the initial state and set the equilibrium angles of the potentials as those in the final state. Through a series of steps of geometrical optimization, we can construct a smooth and short path from the initial state to the final state. This path can be used to calculate free energy difference. To validate this method, we apply it to a small 10‐ALA peptide and find that the calculated free energy changes in helix‐helix and helix‐hairpin transitions are both self‐convergent and cross‐convergent. We also calculate the free energy differences between different stable states of β‐hairpin trpzip2, and the results show that this method is more efficient than the conventional molecular dynamics method in accurate free energy calculation. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Using non-Gaussian density functional fits to improve relative free energy calculations

The accurate and reliable computation of relative free energy differences remains an important long-term goal. Major stumbling blocks for achieving this goal reflect the difficulty of sampling in a known fashion along the reaction coordinate and of maximally combining information that has been collected from the simulation along the reaction coordinate. In this paper we examine the utility of a probability density functional type fit to the distribution of work events collected during a nonequilibrium sample along the reaction coordinate. This approach can readily be generalized to equilibrium sampling and has the potential to estimate the quality of a relative free energy estimate as data are being collected. The method may have the greatest utility for nonequilibrium sampling where non-Gaussian work distributions are generally present that are strongly dominated by rare event sampling in the tail region. We believe that the approach can be used to augment the design and the error analysis of relative free energy computations thus improving the ability to reliably and with known accuracy compute a relative free energy.     

A novel quantum mechanical/molecular mechanical approach to the free energy calculation for isomerization of glycine in aqueous solution

The free energy change associated with the isomerization reaction of glycine in water solution has been studied by a hybrid quantum mechanical/molecular mechanical (QM/MM) approach combined with the theory of energy representation (QM/MM-ER) recently developed. The solvation free energies for both neutral and zwitterionic form of glycine have been determined by means of the QM/MM-ER simulation. The contributions of the electronic polarization and the fluctuation of the QM solute to the solvation free energy have been investigated. It has been found that the contribution of the density fluctuation of the zwitterionic solute is estimated as -4.2 kcal/mol in the total solvation free energy of -46.1 kcal/mol, while that of the neutral form is computed as -3.0 kcal/mol in the solvation free energy of -15.6 kcal/mol. The resultant free energy change associated with the isomerization of glycine in water has been obtained as -7.8 kcal/mol, in excellent agreement with the experimental data of -7.3 or -7.7 kcal/mol, implying the accuracy of the QM/MM-ER approach. The results have also been compared with those computed by other methodologies such as the polarizable continuum model and the classical molecular simulation. The efficiency and advantage of the QM/MM-ER method has been discussed.     

Assessing the efficiency of free energy calculation methods

The efficiencies of two recently developed methods for calculating free energy changes along a generalized coordinate in a system are discussed in the context of other, related approaches. One method is based on Jarzynski's identity [Phys. Rev. Lett. 78, 2690 (1997)]. The second method relies on thermodynamic integration of the average force and is called the adaptive biasing force method [Darve and Pohorille, J. Chem. Phys. 115, 9169 (2001)]. Both methods are designed such that the system evolves along the chosen coordinate(s) without experiencing free energy barriers and they require calculating the instantaneous, unconstrained force acting on this coordinate using the formula derived by Darve and Pohorille. Efficiencies are analyzed by comparing analytical estimates of statistical errors and by considering two numerical examples-internal rotation of hydrated 1,2-dichloroethane and transfer of fluoromethane across a water-hexane interface. The efficiencies of both methods are approximately equal in the first but not in the second case. During transfer of fluoromethane the system is easily driven away from equilibrium and, therefore, the performance of the method based on Jarzynski's identity is poor.     

Efficient free energy calculation of water across lipid membranes

An efficient free energy (FE) calculation of a water molecule to go across lipid membranes is presented. Both overlapping distribution and cavity insertion Widom methods are complementarily used. The former is useful for a dense region where water molecules are found, i.e., from the interfacial to bulk water region, while the latter works well in the low density region, i.e., the hydrocarbon region. Since both methods evaluate the excess chemical potential of water, the obtained FE profile is free from the fitting problem usually arisen when a FE difference method is used. A diphytanyl phosphatidylcholine bilayer is used for our test calculations. An excellent and fast convergence of the chemical potential is obtained when each method is applied for the appropriate region. The estimated FE barrier using the Ewald method for the electrostatic interaction is approximately 7.2 kcal/mol, which is higher than that using the interaction cutoff of 20 A by about 0.9 kcal/mol.     

Chaperoned alchemical free energy simulations: a general method for QM, MM, and QM/MM potentials

A general method for alchemical free energy simulations using QM, MM, and QM/MM potential is developed by introducing "chaperones" to restrain the structures, particularly near the end points. A calculation of the free energy difference between two triazole tautomers in aqueous solution is used to illustrate the method.     

Rationalization of the behavior of solid-liquid surface free energy of water in Cassie and Wenzel wetting states on rugged solid surfaces at the nanometer scale

The present work aims to contribute to the understanding at a molecular level of the origin of the hydrophobic nature of surfaces exhibiting roughness at the nanometer scale. Graphite-based smooth and model surfaces whose roughness dimension stretches from a few angstroms to a few nanometers were used in order to generate Cassie and Wenzel wetting states of water. The corresponding solid-liquid surface free energies were computed by means of molecular dynamics simulations. The solid-liquid surface free energy of water-smooth graphite was found to be -12.7 ± 3.3 mJ/m(2), which is in reasonable agreement with a value estimated from experiments and fully consistent with the features of the employed model. All the rugged surfaces yielded higher surface free energy. In both Cassie and Wenzel states, the maximum variation of the surface free energy with respect to the smooth surface was observed to represent up to 50% of the water model surface tension. The solid-liquid surface free energy of Cassie states could be well predicted from the Cassie-Baxter equation where the surface free energies replace contact angles. The origin of the hydrophobic nature of surfaces yielding Cassie states was therefore found to be the reduction of the number of interactions between water and the solid surface where atomic defects were implemented. Wenzel's theory was found to fail to predict even qualitatively the variation of the solid-liquid surface free energy with respect to the roughness pattern. While graphite was found to be slightly hydrophilic, Wenzel states were found to be dominated by an unfavorable effect that overcame the favorable enthalpic effect induced by the implementation of roughness. From the quantitative point of view, the solid-liquid surface free energy of Wenzel states was found to vary linearly with the roughness contour length.