Expanded Ensemble Simulations Research Articles

The SAMPL6 SAMPLing challenge: assessing the reliability and efficiency of binding free energy calculations.

Approaches for computing small molecule binding free energies based on molecular simulations are now regularly being employed by academic and industry practitioners to study receptor-ligand systems and prioritize the synthesis of small molecules for ligand design. Given the variety of methods and implementations available, it is natural to ask how the convergence rates and final predictions of these methods compare. In this study, we describe the concept and results for the SAMPL6 SAMPLing challenge, the first challenge from the SAMPL series focusing on the assessment of convergence properties and reproducibility of binding free energy methodologies. We provided parameter files, partial charges, and multiple initial geometries for two octa-acid (OA) and one cucurbit[8]uril (CB8) host-guest systems. Participants submitted binding free energy predictions as a function of the number of force and energy evaluations for seven different alchemical and physical-pathway (i.e., potential of mean force and weighted ensemble of trajectories) methodologies implemented with the GROMACS, AMBER, NAMD, or OpenMM simulation engines. To rank the methods, we developed an efficiency statistic based on bias and variance of the free energy estimates. For the two small OA binders, the free energy estimates computed with alchemical and potential of mean force approaches show relatively similar variance and bias as a function of the number of energy/force evaluations, with the attach-pull-release (APR), GROMACS expanded ensemble, and NAMD double decoupling submissions obtaining the greatest efficiency. The differences between the methods increase when analyzing the CB8-quinine system, where both the guest size and correlation times for system dynamics are greater. For this system, nonequilibrium switching (GROMACS/NS-DS/SB) obtained the overall highest efficiency. Surprisingly, the results suggest that specifying force field parameters and partial charges is insufficient to generally ensure reproducibility, and we observe differences between seemingly converged predictions ranging approximately from 0.3 to 1.0 kcal/mol, even with almost identical simulations parameters and system setup (e.g., Lennard-Jones cutoff, ionic composition). Further work will be required to completely identify the exact source of these discrepancies. Among the conclusions emerging from the data, we found that Hamiltonian replica exchange-while displaying very small variance-can be affected by a slowly-decaying bias that depends on the initial population of the replicas, that bidirectional estimators are significantly more efficient than unidirectional estimators for nonequilibrium free energy calculations for systems considered, and that the Berendsen barostat introduces non-negligible artifacts in expanded ensemble simulations.

On Obtaining Boltzmann-Distributed Configurational Ensembles from Expanded Ensemble Simulations with Fast State Mixing.

In Expanded Ensemble (EXE) or Simulated Tempering simulations, the system's (effective) temperature is frequently updated to enhance configurational sampling. We investigated how short the EXE state update interval τ can become before too frequent updates impede Boltzmann sampling. Simulating alanine dipeptide in explicit water, we show that a hybrid MC/MD integrator reliably yields Boltzmann-distributed configurations regardless of τ. However, in MD-driven EXE simulations with short τ, configurational ensembles depend on the thermostat settings.

Vapor–Liquid Equilibria of Mixtures of Molecular Fluids Using the Activity Fraction Expanded Ensemble Simulation Method

The activity fraction expanded ensemble (AFEE) simulation technique has been extended to predict the vapor–liquid coexistence behavior of binary mixtures of molecular fluids. The systems of molecular mixtures studied are thiophene–n-hexane (a mixture of an aromatic compound with a n-alkane) and methane–n-hexane (mixture of n-alkanes with significant difference in size). Comparison with results from the grand canonical transition matrix Monte Carlo method and experimental data available in the literature has also been provided. Further, we have investigated the computational efficiency of the AFEE technique at temperatures where both compounds exhibit subcritical behavior (thiophene–n-hexane at 338.15 K) and where one of the components shows supercritical behavior (methane–n-hexane at 373 K). Canonical Monte Carlo simulations are also performed to obtain the radial distribution functions which have been used to study the structure of the coexisting liquid phases to gain insights into the molecular origins ...

A generalized linear response framework for expanded ensemble and replica exchange simulations.

Expanded ensemble simulation is a powerful technique for enhancing sampling over a range of thermodynamic parameters. However, although the premise is relatively simple, running successful simulations in practice still presents something of an ad hoc challenge. Three main difficulties exist: (1) the selection of the thermodynamic states, (2) the selection of the sampling weights, and (3) efficient sampling of the expanded parameter space. Here we consider these problems in the context of a pairwise linear response approach to the work fluctuation theorem. The approach offers comprehensive tactics for addressing the three difficulties and can be used as either an alternative or a complement to replica exchange simulations. Importantly, the results are trivially implemented for multi-dimensional parameter spaces and they recover results from the literature aimed at the special cases of simulated/parallel tempering and replica exchange umbrella sampling. Illustrative examples are shown using the NAMD simulation engine.

Estimating thermodynamic expectations and free energies in expanded ensemble simulations: Systematic variance reduction through conditioning.

Markov chain Monte Carlo methods are primarily used for sampling from a given probability distribution and estimating multi-dimensional integrals based on the information contained in the generated samples. Whenever it is possible, more accurate estimates are obtained by combining Monte Carlo integration and integration by numerical quadrature along particular coordinates. We show that this variance reduction technique, referred to as conditioning in probability theory, can be advantageously implemented in expanded ensemble simulations. These simulations aim at estimating thermodynamic expectations as a function of an external parameter that is sampled like an additional coordinate. Conditioning therein entails integrating along the external coordinate by numerical quadrature. We prove variance reduction with respect to alternative standard estimators and demonstrate the practical efficiency of the technique by estimating free energies and characterizing a structural phase transition between two solid phases.

Saturation properties of 1-alkyl-3-methylimidazolium based ionic liquids.

We study the liquid-vapor saturation properties of room temperature ionic liquids (RTILs) belonging to the homologous series 1-alkyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide ([Cnmim][NTf2]) using Monte Carlo simulation. We examine the effect of temperature and cation alkyl chain length n on the saturated densities, vapor pressures, and enthalpies of vaporization. These properties are explicitly calculated for temperatures spanning from 280 to 1000 K for RTILs with n = 2, 4, 6, 8, 10, and 12. We also explore how the identity of the anion influences saturation properties. Specifically, we compare results for [C(4)mim][NTf2] with those for 1-butyl-3-methylimidazolium tetrafluoroborate ([C(4)mim][BF4]) and 1-butyl-3-methylimidazolium hexafluorophosphate ([C(4)mim][PF6]). Simulations are completed with a recently developed realistic united-atom force field. A combination of direct grand canonical and isothermal-isobaric temperature expanded ensemble simulations are used to construct phase diagrams. Our results are compared with experimental data and Gibbs ensemble simulation data. Overall, we find good agreement between our results and those measured experimentally. We find that the vapor pressures and enthalpies of vaporization show a strong dependence on the size of the alkyl chain at low temperatures, whereas no particular trend is observed at high temperatures. Finally, we also discuss the effect of temperature on liquid phase nanodomains observed in RTILs with large hydrophobic groups. We do not observe a drastic change in liquid phase structure upon variation of the temperature, which suggests there is not a sharp phase transition between a nanostructured and homogeneous liquid, as has been suggested in earlier studies.

Replica exchange and expanded ensemble simulations as Gibbs sampling: Simple improvements for enhanced mixing

The widespread popularity of replica exchange and expanded ensemble algorithms for simulating complex molecular systems in chemistry and biophysics has generated much interest in discovering new ways to enhance the phase space mixing of these protocols in order to improve sampling of uncorrelated configurations. Here, we demonstrate how both of these classes of algorithms can be considered as special cases of Gibbs sampling within a Markov chain Monte Carlo framework. Gibbs sampling is a well-studied scheme in the field of statistical inference in which different random variables are alternately updated from conditional distributions. While the update of the conformational degrees of freedom by Metropolis Monte Carlo or molecular dynamics unavoidably generates correlated samples, we show how judicious updating of the thermodynamic state indices--corresponding to thermodynamic parameters such as temperature or alchemical coupling variables--can substantially increase mixing while still sampling from the desired distributions. We show how state update methods in common use can lead to suboptimal mixing, and present some simple, inexpensive alternatives that can increase mixing of the overall Markov chain, reducing simulation times necessary to obtain estimates of the desired precision. These improved schemes are demonstrated for several common applications, including an alchemical expanded ensemble simulation, parallel tempering, and multidimensional replica exchange umbrella sampling.

Thermodynamic stability of the cubatic phase of hard cut spheres evaluated by expanded ensemble simulations

The system of hard cut spheres (disk-shaped particles formed by symmetrically truncating the end caps of a sphere) exhibits an intriguing "cubatic" phase with cubic orientational symmetry. However, it is unclear whether this phase is metastable with respect to the columnar phase. We attempt to provide an answer to this question by carrying out free energy calculations by the expanded ensemble Monte Carlo method. We conclude that there may be a very small region of cubatic stability in the vicinity of the isotropic-cubatic phase transition, but that that transition would need to be determined more accurately to obtain a definitive answer. We also comment on the efficacy of the expanded ensemble method for these kinds of calculations.

Measuring excess free energies of self-assembled membrane structures

Using computer simulation of a solvent-free, coarse-grained model for amphiphilic membranes, we study the excess free energy of hourglass-shaped connections (i.e., stalks) between two apposed bilayer membranes. In order to calculate the free energy by simulation in the canonical ensemble, we reversibly transfer two apposed bilayers into a configuration with a stalk in three steps. First, we gradually replace the intermolecular interactions by an external, ordering field. The latter is chosen such that the structure of the non-interacting system in this field closely resembles the structure of the original, interacting system in the absence of the external field. The absence of structural changes along this path suggests that it is reversible; a fact which is confirmed by expanded-ensemble simulations. Second, the external, ordering field is changed as to transform the non-interacting system from the apposed bilayer structure to two-bilayers connected by a stalk. The final external field is chosen such that the structure of the non-interacting system resembles the structure of the stalk in the interacting system without a field. On the third branch of the transformation path, we reversibly replace the external, ordering field by non-bonded interactions. Using expanded-ensemble techniques, the free energy change along this reversible path can be obtained with an accuracy of 10(-3)k(B)T per molecule in the n VT-ensemble. Calculating the chemical potential, we obtain the free energy of a stalk in the grandcanonical ensemble, and employing semi-grandcanonical techniques, we calculate the change of the excess free energy upon altering the molecular architecture. This computational strategy can be applied to compute the free energy of self-assembled phases in lipid and copolymer systems, and the excess free energy of defects or interfaces.

Computation of Surface Tensions Using Expanded Ensemble Simulations

A method for the direct simulation of the surface tension is examined. The technique is based on the thermodynamic route to the interfacial tension and makes use of the expanded ensemble simulation method for the calculation of the free energy difference between two inhomogeneous systems with the same number of particles, temperature, and volume, but different interfacial area. The method is completely general and suitable for systems with either continuous or discontinuous interactions. The adequacy of the expanded ensemble method is assessed by computing the interfacial tension of the planar vapor-liquid interface of Lennard-Jones, Lennard-Jones dimers, Gay-Berne, and square-well model fluids; in the latter, the interactions are discontinuous and the present method does not exhibit the asymmetry of other related methods, such as the test area. The expanded ensemble simulation results are compared with simulation data obtained from other techniques (mechanical and test area) with overall good agreement.

Toward Efficient Chemical Potential Calculations by Expanded Ensemble Simulations; to Make the Free Energy Pathway Fairly Level

A scheme is suggested of how to construct good bias potentials ("balancing factors") to be used in expanded ensemble (EE) calculations of chemical potentials of solutions. A combination of two strategies are used: (i) to use a pathway for particle insertions that avoids large variations in free energy and (ii) to use calculated free energy derivatives to construct a bias potential that makes the pathway fairly level. Only a few very short simulations are needed to accomplish the latter, and then, a full EE simulation is done to obtain the chemical potential. By practical calculations of the chemical potential of benzene, cyclohexane, and benzylamine in water, it is shown that this method is at least equally efficient to the recent adaptive EE (AEE) method by Aberg et al. (J. Chem. Phys. 2004, 120, 3370). Furthermore, the new method provides an alternative strategy that complements existing EE methods.

Monte Carlo simulations and dynamic field theory for suspended particles in liquid crystalline systems

Monte Carlo simulations and dynamic field theory are used to study spherical particles suspended in a nematic liquid crystal. Within these two approaches, we investigate the binding of the defects to the particles, the adsorption of a particle at a solid surface, and two particles interacting with each other. Quantitative comparisons indicate good agreement between the two approaches. A Monte Carlo method based on the combination of canonical expanded ensemble simulations with a density-of-state formalism is used to determine the potential of mean force between one particle and a hard wall. On the other hand, the potential of mean force is evaluated using a dynamic field theory, where the time-dependent evolution of the second rank tensor includes two major aspects of liquid crystalline materials, namely the excluded volume and the long-range order elasticity. The results indicate an effective repulsive force that acts between the particle and the wall. Layer formation at the surface of the hard wall gives rise to local minima in the potential of mean force. The director profile for a particle at contact with a solid surface is characterized by a disclination line distorted and attracted towards the wall. The structure of the nematic for two particles at short distances is also investigated. Our results indicate a structure where the two particles are separated by a circular disclination line. The potential of mean force associated with this configuration indicates an effective attractive interaction between the two particles.

Determination of liquid–solid transition using histogram reweighting method and expanded ensemble simulations

Using the combination of the histogram reweighting and expanded ensemble simulation methods a new generic algorithm for predicting liquid–solid equilibria is proposed. This method avoids the fitting procedures and numerical integration previously required to obtain the equation of state and free energy. Using the expanded ensemble simulation we calculate the free energy difference between hard-sphere crystal and the corresponding Einstein crystal without the need to constrain the center of mass of the system. The liquid–solid equilibria for hard-sphere systems of varying size up to 1372 particles are studied with NpT Monte Carlo simulation and histogram reweighting. We find only a slight size dependence of coexistence properties with this method, and the limit for the infinitely large system is easily and accurately determined.

Potential of mean force between a spherical particle suspended in a nematic liquid crystal and a substrate

We consider a system where a spherical particle is suspended in a nematic liquid crystal confined between two walls. We calculate the liquid-crystal-mediated potential of mean force between the sphere and a substrate by means of Monte Carlo simulations. Three methods are used: a traditional Monte Carlo approach, umbrella sampling, and a novel technique that combines canonical expanded ensemble simulations with a recently proposed density-of-states formalism. The latter method offers advantages in that it facilitates good sampling of phase space without prior knowledge of the energy landscape of the system. The resulting potential of mean force, computed as a function of the normal distance between the sphere and a surface, suggests that the sphere is attracted to the surface, even in the absence of attractive molecular interactions.

Measurement of chemical potentials of systems with strong excluded volume interactions by computing the density of states

At high densities and low temperatures, the conventional Widom test particle method to compute the chemical potential of a system of particles with excluded volume interactions fails owing to bad statistics. A way to circumvent this problem is the use of expanded ensemble simulation techniques or thermodynamic integration. In this article, we will describe an alternative method to compute the chemical potential which is conceptually much easier, by computing the density of states of systems of N and N + 1 particles directly; and by performing a test particle simulation at very high temperature. The advantage of our technique is that the densities of states of the N and N + 1 particle system are computed in an ensemble in which particles can pass each other, resulting in a more efficient sampling. We will demonstrate our method not only for single particles but also for chain molecules with intramolecular interactions. By using an infinite temperature expansion and an extension of the density of states to very high energies, we will show that it is also possible to compute the chemical potential without having to compute the density of states for the N + 1 particle system.

Relationship between polymer chain conformation and phase boundaries in a supercritical fluid

We investigate the solvent density driven changes in polymer conformation and phase behavior that occur in a supercritical fluid, with a particular emphasis on conditions near the lower critical solution temperature (LCST) phase boundary. Using continuous space Monte Carlo simulations, the mean square end-to-end distance (R) and radius of gyration (Rg) are calculated for a single chain with 20 Lennard-Jones segments in a monomeric solvent over a broad range of densities and temperatures. The chains collapse as temperature increases at constant pressure, or as density decreases at constant temperature. A minimum in R and Rg occurs at a temperature slightly above the coil-to-globule transition temperature (C-GTT), where the chain adopts a quasi-ideal conformation, defined by the balance of binary attractive and repulsive interactions. Expanded ensemble simulations of finite-concentration polymer–solvent mixtures reveal that the LCST phase boundary correlates well with the single chain C-GTT. At temperatures well above the LCST, the chain expands again suggesting an upper critical solution temperature (UCST) phase boundary above the LCST.