Molecular Dynamics in the Multicanonical Ensemble: Equivalence of Wang-Landau Sampling, Statistical Temperature Molecular Dynamics, and Metadynamics

Reduced critical slowing down for statistical physics simulations

Wang-Landau simulations offer the possibility to integrate explicitly over a collective coordinate and stochastically over the remainder of configuration space. We propose to choose the so-called "slow mode," which is responsible for large autocorrelation times and thus critical slowing down, for collective integration. We study this proposal for the Ising model and the linear-log-relaxation (LLR) method as simulation algorithm. We first demonstrate supercritical slowing down in a phase with spontaneously broken symmetry and for the heat-bath algorithms, for which autocorrelation times grow exponentially with system size. By contrast, using the magnetization as collective coordinate, we present evidence that supercritical slowing down is absent. We still observe a polynomial increase of the autocorrelation time with volume (critical slowing down), which is, however, reduced by orders of magnitude when compared to local update techniques.

Mechanistic Understanding from Molecular Dynamics in Pharmaceutical Research 2: Lipid Membrane in Drug Design

We review the use of molecular dynamics (MD) simulation as a drug design tool in the context of the role that the lipid membrane can play in drug action, i.e., the interaction between candidate drug molecules and lipid membranes. In the standard "lock and key" paradigm, only the interaction between the drug and a specific active site of a specific protein is considered; the environment in which the drug acts is, from a biophysical perspective, far more complex than this. The possible mechanisms though which a drug can be designed to tinker with physiological processes are significantly broader than merely fitting to a single active site of a single protein. In this paper, we focus on the role of the lipid membrane, arguably the most important element outside the proteins themselves, as a case study. We discuss work that has been carried out, using MD simulation, concerning the transfection of drugs through membranes that act as biological barriers in the path of the drugs, the behavior of drug molecules within membranes, how their collective behavior can affect the structure and properties of the membrane and, finally, the role lipid membranes, to which the vast majority of drug target proteins are associated, can play in mediating the interaction between drug and target protein. This review paper is the second in a two-part series covering MD simulation as a tool in pharmaceutical research; both are designed as pedagogical review papers aimed at both pharmaceutical scientists interested in exploring how the tool of MD simulation can be applied to their research and computational scientists interested in exploring the possibility of a pharmaceutical context for their research.

Calculation of the residual entropy of Ice Ih by Monte Carlo simulation with the combination of the replica-exchange Wang-Landau algorithm and multicanonical replica-exchange method

We estimated the residual entropy of Ice Ih by the recently developed simulation protocol, namely, the combination of the replica-exchange Wang-Landau algorithm and multicanonical replica-exchange method. We employed a model with the nearest neighbor interactions on the three-dimensional hexagonal lattice, which satisfied the ice rules in the ground state. The results showed that our estimate of the residual entropy is in accordance with various previous results. In this article, we not only give our latest estimate of the residual entropy of Ice Ih but also discuss the importance of the uniformity of a random number generator in Monte Carlo simulations.

Flat-histogram method comparison on the two-dimensional Ising model

We compare the convergence of several flat-histogram methods applied to the two-dimensional Ising model, including the recently introduced stochastic approximation with a dynamic update factor (SAD) method. We compare this method to the Wang-Landau (WL) method, the 1/t variant of the WL method, and standard stochastic approximation Monte Carlo (SAMC). In addition, we consider a procedure WL followed by a "production run" with fixed weights that refines the estimation of the entropy. We find that WL followed by a production run does converge to the true density of states, in contrast to pure WL. Three of the methods converge robustly: SAD, 1/t-WL, and WL followed by a production run. Of these, SAD does not require a priori knowledge of the energy range. This work also shows that WL followed by a production run performs superior to other forms of WL while ensuring both ergodicity and detailed balance.

Stochastic approximation Monte Carlo with a dynamic update factor

We present a Monte Carlo algorithm based on the stochastic approximation Monte Carlo (SAMC) algorithm for directly calculating the density of states. The proposed method is stochastic approximation with a dynamic update factor (SAD), which dynamically adjusts the update factor γ_{t} during the course of the simulation. We test this method on a square-well fluid and a 31-atom Lennard-Jones cluster and compare the convergence behavior of several related Monte Carlo methods. We find that both the SAD and 1/t-Wang-Landau (1/t-WL) methods rapidly converge to the correct density of states without the need for the user to specify an arbitrary tunable parameter t_{0} as in the case of SAMC. SAD requires as input the temperature range of interest, in contrast with 1/t-WL, which requires that the user identify the interesting range of energies. The convergence of the 1/t-WL method is very sensitive to the energy range chosen for the low-temperature heat capacity of the Lennard-Jones cluster. Thus, SAD is more powerful in the common case in which the range of energies is not known in advance.

Efficient simulation protocol for determining the density of states: Combination of replica-exchange Wang-Landau method and multicanonical replica-exchange method

By combining two generalized-ensemble algorithms, the replica-exchange Wang-Landau (REWL) method and the multicanonical replica-exchange method (MUCAREM), we propose an effective simulation protocol to determine the density of states with high accuracy. The new protocol is referred to as REWL-MUCAREM, and REWL is first performed and then MUCAREM is performed. In order to verify the effectiveness of our protocol, we performed simulations of a square-lattice Ising model using the three methods, namely REWL, MUCAREM, and REWL-MUCAREM. The results showed that the density of states obtained by REWL-MUCAREM is more accurate than that is estimated by the two methods separately.

Simulation of fluid/gel phase equilibrium in lipid vesicles

Simulation of single component dipalmitoylphosphatidylcholine (DPPC) coarse-grained DRY-MARTINI lipid vesicles of diameter 10 nm (1350 lipids), 20 nm (5100 lipids) and 40 nm (17 600 lipids) is performed using statistical temperature molecular dynamics (STMD), to study finite size effects upon the order-disorder gel/fluid transition. STMD obtains enhanced sampling using a generalized ensemble, obtaining a flat energy distribution between upper and lower cutoffs, with little computational cost over canonical molecular dynamics. A single STMD trajectory of moderate length is sufficient to sample 20+ transition events, without trapping in the gel phase, and obtain well averaged properties. Phase transitions are analyzed via the energy-dependence of the statistical temperature, T_S(U). The transition temperature decreases with decreasing diameter, in agreement with experiment, and the transition changes from first order to borderline first-second order. The size- and layer-dependence of the structure of both stable phases, and of the pathway of the phase transition, are determined. It is argued that the finite size effects are primarily caused by the disruption of the gel packing by curvature. Inhomogeneous states with faceted gel patches connected by unusual fluid seams are observed at high curvature, with visually different structure in the inner and outer layers due to the different curvatures. Thus a simple physical picture describes phase transitions in nanoscale finite systems far from the thermodynamic limit.

An Overview of Recent Standard and Accelerated Molecular Dynamics Simulations of Helium Behavior in Tungsten

One of the most critical challenges for the successful adoption of nuclear fusion power corresponds to plasma-facing materials. Due to its favorable properties in this context (low sputtering yield, high thermal conductivity, high melting point, among others), tungsten is a leading candidate material. Nevertheless, tungsten is affected by the plasma and fusion byproducts. Irradiation by helium nuclei, in particular, strongly modifies the surface structure by a synergy of processes, whose origin is the nucleation and growth of helium bubbles. In this review, we present recent advances in the understanding of helium effects in tungsten from a simulational approach based on accelerated molecular dynamics, which emphasizes the use of realistic parameters, as are expected in experimental and operational fusion power conditions.

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.

Optimal updating magnitude in adaptive flat-distribution sampling

We present a study on the optimization of the updating magnitude for a class of free energy methods based on flat-distribution sampling, including the Wang-Landau (WL) algorithm and metadynamics. These methods rely on adaptive construction of a bias potential that offsets the potential of mean force by histogram-based updates. The convergence of the bias potential can be improved by decreasing the updating magnitude with an optimal schedule. We show that while the asymptotically optimal schedule for the single-bin updating scheme (commonly used in the WL algorithm) is given by the known inverse-time formula, that for the Gaussian updating scheme (commonly used in metadynamics) is often more complex. We further show that the single-bin updating scheme is optimal for very long simulations, and it can be generalized to a class of bandpass updating schemes that are similarly optimal. These bandpass updating schemes target only a few long-range distribution modes and their optimal schedule is also given by the inverse-time formula. Constructed from orthogonal polynomials, the bandpass updating schemes generalize the WL and Langfeld-Lucini-Rago algorithms as an automatic parameter tuning scheme for umbrella sampling.

Macroscopically constrained Wang-Landau method for systems with multiple order parameters and its application to drawing complex phase diagrams

A generalized approach to Wang-Landau simulations, macroscopically constrained Wang-Landau, is proposed to simulate the density of states of a system with multiple macroscopic order parameters. The method breaks a multidimensional random-walk process in phase space into many separate, one-dimensional random-walk processes in well-defined subspaces. Each of these random walks is constrained to a different set of values of the macroscopic order parameters. When the multivariable density of states is obtained for one set of values of fieldlike model parameters, the density of states for any other values of these parameters can be obtained by a simple transformation of the total system energy. All thermodynamic quantities of the system can then be rapidly calculated at any point in the phase diagram. We demonstrate how to use the multivariable density of states to draw the phase diagram, as well as order-parameter probability distributions at specific phase points, for a model spin-crossover material: an antiferromagnetic Ising model with ferromagnetic long-range interactions. The fieldlike parameters in this model are an effective magnetic field and the strength of the long-range interaction.

Ligand diffusion in proteins via enhanced sampling in molecular dynamics

Computational simulations in biophysics describe the dynamics and functions of biological macromolecules at the atomic level. Among motions particularly important for life are the transport processes in heterogeneous media. The process of ligand diffusion inside proteins is an example of a complex rare event that can be modeled using molecular dynamics simulations. The study of physical interactions between a ligand and its biological target is of paramount importance for the design of novel drugs and enzymes. Unfortunately, the process of ligand diffusion is difficult to study experimentally. The need for identifying the ligand egress pathways and understanding how ligands migrate through protein tunnels has spurred the development of several methodological approaches to this problem. The complex topology of protein channels and the transient nature of the ligand passage pose difficulties in the modeling of the ligand entry/escape pathways by canonical molecular dynamics simulations. In this review, we report a methodology involving a reconstruction of the ligand diffusion reaction coordinates and the free-energy profiles along these reaction coordinates using enhanced sampling of conformational space. We illustrate the above methods on several ligand-protein systems, including cytochromes and G-protein-coupled receptors. The methods are general and may be adopted to other transport processes in living matter.

Wang-Landau Reaction Ensemble Method: Simulation of Weak Polyelectrolytes and General Acid-Base Reactions

We present a novel method for the study of weak polyelectrolytes and general acid-base reactions in molecular dynamics and Monte Carlo simulations. The approach combines the advantages of the reaction ensemble and the Wang-Landau sampling method. Deprotonation and protonation reactions are simulated explicitly with the help of the reaction ensemble method, while the accurate sampling of the corresponding phase space is achieved by the Wang-Landau approach. The combination of both techniques provides a sufficient statistical accuracy such that meaningful estimates for the density of states and the partition sum can be obtained. With regard to these estimates, several thermodynamic observables like the heat capacity or reaction free energies can be calculated. We demonstrate that the computation times for the calculation of titration curves with a high statistical accuracy can be significantly decreased when compared to the original reaction ensemble method. The applicability of our approach is validated by the study of weak polyelectrolytes and their thermodynamic properties.

Dilute Semiflexible Polymers with Attraction: Collapse, Folding and Aggregation

We review the current state on the thermodynamic behavior and structural phases of self- and mutually-attractive dilute semiflexible polymers that undergo temperature-driven transitions. In extreme dilution, polymers may be considered isolated, and this single polymer undergoes a collapse or folding transition depending on the internal structure. This may go as far as to stable knot phases. Adding polymers results in aggregation, where structural motifs again depend on the internal structure. We discuss in detail the effect of semiflexibility on the collapse and aggregation transition and provide perspectives for interesting future investigations.

μ-tempered metadynamics: Artifact independent convergence times for wide hills

Recent analysis of well-tempered metadynamics (WTmetaD) showed that it converges without mollification artifacts in the bias potential. Here, we explore how metadynamics heals mollification artifacts, how healing impacts convergence time, and whether alternative temperings may be used to improve efficiency. We introduce "μ-tempered" metadynamics as a simple tempering scheme, inspired by a related mollified adaptive biasing potential, that results in artifact independent convergence of the free energy estimate. We use a toy model to examine the role of artifacts in WTmetaD and solvated alanine dipeptide to compare the well-tempered and μ-tempered frameworks demonstrating fast convergence for hill widths as large as 60(∘) for μTmetaD.

Molecular Dynamics and Monte Carlo simulations in the microcanonical ensemble: Quantitative comparison and reweighting techniques

Molecular Dynamics (MD) and Monte Carlo (MC) simulations are the most popular simulation techniques for many-particle systems. Although they are often applied to similar systems, it is unclear to which extent one has to expect quantitative agreement of the two simulation techniques. In this work, we present a quantitative comparison of MD and MC simulations in the microcanonical ensemble. For three test examples, we study first- and second-order phase transitions with a focus on liquid-gas like transitions. We present MD analysis techniques to compensate for conservation law effects due to linear and angular momentum conservation. Additionally, we apply the weighted histogram analysis method to microcanonical histograms reweighted from MD simulations. By this means, we are able to estimate the density of states from many microcanonical simulations at various total energies. This further allows us to compute estimates of canonical expectation values.

Towards an Optimal Flow: Density-of-States-Informed Replica-Exchange Simulations

Replica exchange (RE) is one of the most popular enhanced-sampling simulations technique in use today. Despite widespread successes, RE simulations can sometimes fail to converge in practical amounts of time, e.g., when sampling around phase transitions, or when a few hard-to-find configurations dominate the statistical averages. We introduce a generalized RE scheme, density-of-states-informed RE, that addresses some of these challenges. The key feature of our approach is to inform the simulation with readily available, but commonly unused, information on the density of states of the system as the RE simulation proceeds. This enables two improvements, namely, the introduction of resampling moves that actively move the system towards equilibrium and the continual adaptation of the optimal temperature set. As a consequence of these two innovations, we show that the configuration flow in temperature space is optimized and that the overall convergence of RE simulations can be dramatically accelerated.

Accelerated weight histogram method for exploring free energy landscapes

Calculating free energies is an important and notoriously difficult task for molecular simulations. The rapid increase in computational power has made it possible to probe increasingly complex systems, yet extracting accurate free energies from these simulations remains a major challenge. Fully exploring the free energy landscape of, say, a biological macromolecule typically requires sampling large conformational changes and slow transitions. Often, the only feasible way to study such a system is to simulate it using an enhanced sampling method. The accelerated weight histogram (AWH) method is a new, efficient extended ensemble sampling technique which adaptively biases the simulation to promote exploration of the free energy landscape. The AWH method uses a probability weight histogram which allows for efficient free energy updates and results in an easy discretization procedure. A major advantage of the method is its general formulation, making it a powerful platform for developing further extensions and analyzing its relation to already existing methods. Here, we demonstrate its efficiency and general applicability by calculating the potential of mean force along a reaction coordinate for both a single dimension and multiple dimensions. We make use of a non-uniform, free energy dependent target distribution in reaction coordinate space so that computational efforts are not wasted on physically irrelevant regions. We present numerical results for molecular dynamics simulations of lithium acetate in solution and chignolin, a 10-residue long peptide that folds into a β-hairpin. We further present practical guidelines for setting up and running an AWH simulation.

Computational Calorimetry: High-Precision Calculation of Host-Guest Binding Thermodynamics

We present a strategy for carrying out high-precision calculations of binding free energy and binding enthalpy values from molecular dynamics simulations with explicit solvent. The approach is used to calculate the thermodynamic profiles for binding of nine small molecule guests to either the cucurbit[7]uril (CB7) or β-cyclodextrin (βCD) host. For these systems, calculations using commodity hardware can yield binding free energy and binding enthalpy values with a precision of ∼0.5 kcal/mol (95% CI) in a matter of days. Crucially, the self-consistency of the approach is established by calculating the binding enthalpy directly, via end point potential energy calculations, and indirectly, via the temperature dependence of the binding free energy, i.e., by the van't Hoff equation. Excellent agreement between the direct and van't Hoff methods is demonstrated for both host-guest systems and an ion-pair model system for which particularly well-converged results are attainable. Additionally, we find that hydrogen mass repartitioning allows marked acceleration of the calculations with no discernible cost in precision or accuracy. Finally, we provide guidance for accurately assessing numerical uncertainty of the results in settings where complex correlations in the time series can pose challenges to statistical analysis. The routine nature and high precision of these binding calculations opens the possibility of including measured binding thermodynamics as target data in force field optimization so that simulations may be used to reliably interpret experimental data and guide molecular design.

Toward Improved Force-Field Accuracy through Sensitivity Analysis of Host-Guest Binding Thermodynamics

Improving the capability of atomistic computer models to predict the thermodynamics of noncovalent binding is critical for successful structure-based drug design, and the accuracy of such calculations remains limited by nonoptimal force field parameters. Ideally, one would incorporate protein-ligand affinity data into force field parametrization, but this would be inefficient and costly. We now demonstrate that sensitivity analysis can be used to efficiently tune Lennard-Jones parameters of aqueous host-guest systems for increasingly accurate calculations of binding enthalpy. These results highlight the promise of a comprehensive use of calorimetric host-guest binding data, along with existing validation data sets, to improve force field parameters for the simulation of noncovalent binding, with the ultimate goal of making protein-ligand modeling more accurate and hence speeding drug discovery.

First-order phase transition and tricritical scaling behavior of the Blume-Capel model: A Wang-Landau sampling approach

We investigate the tricritical scaling behavior of the two-dimensional spin-1 Blume-Capel model by using the Wang-Landau method of measuring the joint density of states for lattice sizes up to 48×48 sites. We find that the specific heat deep in the first-order area of the phase diagram exhibits a double-peak structure of the Schottky-like anomaly appearing with the transition peak. The first-order transition curve is systematically determined by employing the method of field mixing in conjunction with finite-size scaling, showing a significant deviation from the previous data points. At the tricritical point, we characterize the tricritical exponents through finite-size-scaling analysis including the phenomenological finite-size scaling with thermodynamic variables. Our estimation of the tricritical eigenvalue exponents, yt=1.804(5), yg=0.80(1), and yh=1.925(3), provides the first Wang-Landau verification of the conjectured exact values, demonstrating the effectiveness of the density-of-states-based approach in finite-size scaling study of multicritical phenomena.

Extended Hamiltonian approach to continuous tempering

We introduce an enhanced sampling simulation technique based on continuous tempering, i.e., on continuously varying the temperature of the system under investigation. Our approach is mathematically straightforward, being based on an extended Hamiltonian formulation in which an auxiliary degree of freedom, determining the effective temperature, is coupled to the physical system. The physical system and its temperature evolve continuously in time according to the equations of motion derived from the extended Hamiltonian. Due to the Hamiltonian structure, it is easy to show that a particular subset of the configurations of the extended system is distributed according to the canonical ensemble for the physical system at the correct physical temperature.