Optimized parallel tempering simulations of proteins

We apply a recently developed adaptive algorithm that systematically improves the efficiency of parallel tempering or replica exchange methods in the numerical simulation of small proteins. Feedback iterations allow us to identify an optimal set of temperatures/replicas which are found to concentrate at the bottlenecks of the simulations. A measure of convergence for the equilibration of the parallel tempering algorithm is discussed. We test our algorithm by simulating the 36-residue villin headpiece subdomain HP-36 where we find a lowest-energy configuration with a root-mean-square deviation of less than 4 A to the experimentally determined structure.     

All-exchanges parallel tempering

An alternative exchange strategy for parallel tempering simulations is introduced. Instead of attempting to swap configurations between two randomly chosen but adjacent replicas, the acceptance probabilities of all possible swap moves are calculated a priori. One specific swap move is then selected according to its probability and enforced. The efficiency of the method is illustrated first on the case of two Lennard-Jones (LJ) clusters containing 13 and 31 atoms, respectively. The convergence of the caloric curve is seen to be at least twice as fast as in conventional parallel tempering simulations, especially for the difficult case of LJ31. Further evidence for an improved efficiency is reported on the ergodic measure introduced by Mountain and Thirumalai [J. Phys. Chem. 93, 6975 (1989)], calculated here for LJ13 close to the melting point. Finally, tests on two simple spin systems indicate that the method should be particularly useful when a limited number of replicas are available.     

Multiple‐Replica Exchange with Information Retrieval

The parallel tempering simulation method was recently extended to allow for possible exchanges between non‐adjacent replicas. We introduce a multiple‐exchange variant which naturally incorporates the information from all replicas when calculating statistical averages, building on the related virtual‐move method of Coluzza and Frenkel (ChemPhysChem 2005 , 6, 1779). The method is extensively tested on three model systems, namely, a Lennard‐Jones cluster exhibiting a finite size phase transition, the Lennard‐Jones fluid, and the 2D ferromagnetic Ising model. In all cases, the present method performs significantly better and converges faster than conventional parallel tempering Monte Carlo simulations. The standard deviations are also systematically decreased with respect to virtual moves.     

Multiple scaling replica exchange for the conformational sampling of biomolecules in explicit water

A multiple scaling replica exchange method for the efficient conformational sampling of biomolecular systems in explicit solvent is presented. The method is a combination of the replica exchange with solute tempering (REST) technique and a Tsallis biasing potential. The Tsallis biasing increases the sampling efficiency, while the REST minimizes the number of replicas needed. Unbiased statistics can be obtained by reweighting of the data using a weighted histogram analysis technique. The method is illustrated by its application to a ten residue peptide in explicit water.     

Combining the biased and unbiased sampling strategy into one convenient free energy calculation method

Constructing a free energy landscape for a large molecule is difficult. One has to use either a high temperature or a strong driving force to enhance the sampling on the free energy barriers. In this work, we propose a mixed method that combines these two kinds of acceleration strategies into one simulation. First, it applies an adaptive biasing potential to some replicas of the molecule. These replicas are particularly accelerated in a collective variable space. Second, it places some unbiased and exchangeable replicas at various temperature levels. These replicas generate unbiased sampling data in the canonical ensemble. To improve the sampling efficiency, biased replicas transfer their state variables to the unbiased replicas after equilibrium by Monte Carlo trial moves. In comparison to previous integrated methods, it is more convenient for users. It does not need an initial reference biasing potential to guide the sampling of the molecule. And it is also unnecessary to insert many replicas for the requirement of passing the free energy barriers. The free energy calculation is accomplished in a single stage. It samples the data as fast as a biased simulation and it processes the data as simple as an unbiased simulation. The method provides a minimalist approach to the construction of the free energy landscape. © 2019 Wiley Periodicals, Inc.     

Computation of free energy profiles with parallel adaptive dynamics

We propose a formulation of an adaptive computation of free energy differences, in the adaptive biasing force or nonequilibrium metadynamics spirit, using conditional distributions of samples of configurations which evolve in time. This allows us to present a truly unifying framework for these methods, and to prove convergence results for certain classes of algorithms. From a numerical viewpoint, a parallel implementation of these methods is very natural, the replicas interacting through the reconstructed free energy. We demonstrate how to improve this parallel implementation by resorting to some selection mechanism on the replicas. This is illustrated by computations on a model system of conformational changes.     

On easy implementation of a variant of the replica exchange with solute tempering in GROMACS

To reduce the number of replicas required in the conventional replica exchange method for huge systems, recently the replica exchange with solute tempering (REST) method was proposed. Here we showed that a variant of REST realized by rescaling the force-field parameters can be performed with GROMACS 4 without changing the code. We tested the variant REST for alanine dipeptide and an N-terminal peptide from p53 confirming its performance nearly equal to the original REST.     

On the efficiency of exchange in parallel tempering monte carlo simulations

We introduce the concept of effective fraction, defined as the expected probability that a configuration from the lowest index replica successfully reaches the highest index replica during a replica exchange Monte Carlo simulation. We then argue that the effective fraction represents an adequate measure of the quality of the sampling technique, as far as swapping is concerned. Under the hypothesis that the correlation between successive exchanges is negligible, we propose a technique for the computation of the effective fraction, a technique that relies solely on the values of the acceptance probabilities obtained at the end of the simulation. The effective fraction is then utilized for the study of the efficiency of a popular swapping scheme in the context of parallel tempering in the canonical ensemble. For large dimensional oscillators, we show that the swapping probability that minimizes the computational effort is 38.74%. By studying the parallel tempering swapping efficiency for a 13-atom Lennard-Jones cluster, we argue that the value of 38.74% remains roughly the optimal probability for most systems with continuous distributions that are likely to be encountered in practice.     

An improved replica-exchange sampling method: temperature intervals with global energy reassignment

In a molecular dynamics (MD) simulation, representative sampling over the entire phase space is desired to obtain an accurate canonical distribution at a given temperature. For large molecules, such as proteins, this is problematic because systems tend to become trapped in local energy minima. The extensively used replica-exchange molecular dynamics (REMD) simulation technique overcomes this kinetic-trapping problem by allowing Boltzmann-weighted configuration exchange processes to occur between numerous thermally adjacent and compositionally identical simulations that are thermostated at sequentially higher temperatures. While the REMD method provides much better sampling than conventional MD, there are two substantial difficulties that are inherent in its application: (1) the large number of replicas that must be used to span a designated temperature range and (2) the subsequent long time required for configurations sampled at high temperatures to exchange down for potential inclusion within the low-temperature ensemble of interest. In this work, a new method based on temperature intervals with global energy reassignment (TIGER) is presented that overcomes both of these problems. A TIGER simulation is conducted as a series of short heating-sampling-quenching cycles. At the end of each cycle, the potential energies of all replicas are simultaneously compared at the same temperature using a Metropolis sampling method and then globally reassigned to the designated temperature levels. TIGER is compared with regular MD and REMD methods for the alanine dipeptide in water. The results indicate that TIGER increases sampling efficiency while substantially reducing the number of central processing units required for a comparable conventional REMD simulation.     

Combining smart darting with parallel tempering using Eckart space: application to Lennard-Jones clusters

The smart-darting algorithm is a Monte Carlo based simulation method used to overcome quasiergodicity problems associated with disconnected regions of configurations space separated by high energy barriers. As originally implemented, the smart-darting method works well for clusters at low temperatures with the angular momentum restricted to zero and where there are no transitions to permutational isomers. If the rotational motion of the clusters is unrestricted or if permutational isomerization becomes important, the acceptance probability of darting moves in the original implementation of the method becomes vanishingly small. In this work the smart-darting algorithm is combined with the parallel tempering method in a manner where both rotational motion and permutational isomerization events are important. To enable the combination of parallel tempering with smart darting so that the smart-darting moves have a reasonable acceptance probability, the original algorithm is modified by using a restricted space for the smart-darting moves. The restricted space uses a body-fixed coordinate system first introduced by Eckart, and moves in this Eckart space are coupled with local moves in the full 3N-dimensional space. The modified smart-darting method is applied to the calculation of the heat capacity of a seven-atom Lennard-Jones cluster. The smart-darting moves yield significant improvement in the statistical fluctuations of the calculated heat capacity in the region of temperatures where the system isomerizes. When the modified smart-darting algorithm is combined with parallel tempering, the statistical fluctuations of the heat capacity of a seven-atom Lennard-Jones cluster using the combined method are smaller than parallel tempering when used alone.     

Walking freely in the energy and temperature space by the modified replica exchange molecular dynamics method

Replica Exchange Molecular Dynamics (REMD) method is a powerful sampling tool in molecular simulations. Recently, we made a modification to the standard REMD method. It places some inactive replicas at different temperatures as well as the active replicas. The method completely decouples the number of the active replicas and the number of the temperature levels. In this article, we make a further modification to our previous method. It uses the inactive replicas in a different way. The inactive replicas first sample in their own knowledge‐based energy databases and then participate in the replica exchange operations in the REMD simulation. In fact, this method is a hybrid between the standard REMD method and the simulated tempering method. Using different active replicas, one can freely control the calculation quantity and the convergence speed of the simulation. To illustrate the performance of the method, we apply it to some small models. The distribution functions of the replicas in the energy space and temperature space show that the modified REMD method in this work can let the replicas walk freely in both of the two spaces. With the same number of the active replicas, the free energy surface in the simulation converges faster than the standard REMD. © 2016 Wiley Periodicals, Inc.     

Replica exchange with nonequilibrium switches: enhancing equilibrium sampling by increasing replica overlap

We describe a replica exchange strategy where trial swap configurations are generated by nonequilibrium switching simulations. By devoting simulation time to the switching simulations, one can systematically increase an effective overlap between replicas, which leads to an increased exchange acceptance rate and less correlated equilibrium samples. In this paper, we derive our method for a general class of stochastic dynamics, and discuss various strategies for enhancing replica overlap through novel dynamical schemes and prudent choices of switching protocols. We then demonstrate our method on a model system of alanine dipeptide in implicit solvent, characterizing decreases in data correlations and gains in sampling efficiency.     

Free-energy landscape for beta hairpin folding from combined parallel tempering and metadynamics

We develop a new free-energy method, based on the combination of parallel tempering and metadynamics, and apply this method to the calculation of the free-energy landscape of the folding beta hairpin in explicit water. We show that the combined method greatly improves the performance of both parallel tempering and metadynamics. In particular, we are able to sample the high free-energy regions, which are not accessible with conventional parallel tempering. We use our results to calculate the difference in entropy and enthalpy between the folded and the unfolded state and to characterize the most populated configurations in the relevant free-energy basins.     

The temperature intervals with global exchange of replicas empirical accelerated sampling method: parameter sensitivity and extension to a complex molecular system

The recently developed "temperature intervals with global exchange of replicas" (TIGER2) algorithm is an efficient replica-exchange sampling algorithm that provides the freedom to specify the number of replicas and temperature levels independently of the size of the system and temperature range to be spanned, thus making it particularly well suited for sampling molecular systems that are considered to be too large to be sampled using conventional replica exchange methods. Although the TIGER2 method is empirical in nature, when appropriately applied it is able to provide sampling that satisfies the balance condition and closely approximates a Boltzmann-weighted ensemble of states. In this work, we evaluated the influence of factors such as temperature range, temperature spacing, replica number, and sampling cycle design on the accuracy of a TIGER2 simulation based on molecular dynamics simulations of alanine dipeptide in implicit solvent. The influence of these factors is further examined by calculating the properties of a complex system composed of the B1 immunoglobulin-binding domain of streptococcal protein G (protein G) in aqueous solution. The accuracy of a TIGER2 simulation is particularly sensitive to the maximum temperature level selected for the simulation. A method to determine the appropriate maximum temperature level to be used in a TIGER2 simulation is presented.     

A temperature predictor for parallel tempering simulations

An algorithm is proposed that generates a set of temperatures for use in parallel tempering simulations (also known as temperature-replica exchange molecular dynamics simulations) of proteins to obtain a desired exchange probability Pdes. The input consists of the number of protein atoms and water molecules in the system, information about the use of constraints and virtual sites and the lower temperature limits. The temperatures generated yield probabilities which are very close to Pdes (correlation 97%), independent of force field and over a wide temperature range. To facilitate its use, the algorithm has been implemented as a web server at .     

Searching for globally optimal functional forms for interatomic potentials using genetic programming with parallel tempering

Molecular dynamics and other molecular simulation methods rely on a potential energy function, based only on the relative coordinates of the atomic nuclei. Such a function, called a force field, approximately represents the electronic structure interactions of a condensed matter system. Developing such approximate functions and fitting their parameters remains an arduous, time-consuming process, relying on expert physical intuition. To address this problem, a functional programming methodology was developed that may enable automated discovery of entirely new force-field functional forms, while simultaneously fitting parameter values. The method uses a combination of genetic programming, Metropolis Monte Carlo importance sampling and parallel tempering, to efficiently search a large space of candidate functional forms and parameters. The methodology was tested using a nontrivial problem with a well-defined globally optimal solution: a small set of atomic configurations was generated and the energy of each configuration was calculated using the Lennard-Jones pair potential. Starting with a population of random functions, our fully automated, massively parallel implementation of the method reproducibly discovered the original Lennard-Jones pair potential by searching for several hours on 100 processors, sampling only a minuscule portion of the total search space. This result indicates that, with further improvement, the method may be suitable for unsupervised development of more accurate force fields with completely new functional forms.     

Hamiltonian replica‐exchange simulations with adaptive biasing of peptide backbone and side chain dihedral angles

A Hamiltonian Replica‐Exchange Molecular Dynamics (REMD) simulation method has been developed that employs a two‐dimensional backbone and one‐dimensional side chain biasing potential specifically to promote conformational transitions in peptides. To exploit the replica framework optimally, the level of the biasing potential in each replica was appropriately adapted during the simulations. This resulted in both high exchange rates between neighboring replicas and improved occupancy/flow of all conformers in each replica. The performance of the approach was tested on several peptide and protein systems and compared with regular MD simulations and previous REMD studies. Improved sampling of relevant conformational states was observed for unrestrained protein and peptide folding simulations as well as for refinement of a loop structure with restricted mobility of loop flanking protein regions. © 2013 Wiley Periodicals, Inc.     

Evaluating force field accuracy with long-time simulations of a β-hairpin tryptophan zipper peptide

We have combined graphics processing unit-accelerated all-atom molecular dynamics with parallel tempering to explore the folding properties of small peptides in implicit solvent on the time scale of microseconds. We applied this methodology to the synthetic β-hairpin, trpzip2, and one of its sequence variants, W2W9. Each simulation consisted of over 8 μs of aggregated virtual time. Several measures of folding behavior showed good convergence, allowing comparison with experimental equilibrium properties. Our simulations suggest that the intramolecular interactions of tryptophan side chains are responsible for much of the stability of the native fold. We conclude that the ff99 force field combined with ff96 φ and ψ dihedral energies and an implicit solvent can reproduce plausible folding behavior in both trpzip2 and W2W9.     

Scalable free energy calculation of proteins via multiscale essential sampling

A multiscale simulation method, "multiscale essential sampling (MSES)," is proposed for calculating free energy surface of proteins in a sizable dimensional space with good scalability. In MSES, the configurational sampling of a full-dimensional model is enhanced by coupling with the accelerated dynamics of the essential degrees of freedom. Applying the Hamiltonian exchange method to MSES can remove the biasing potential from the coupling term, deriving the free energy surface of the essential degrees of freedom. The form of the coupling term ensures good scalability in the Hamiltonian exchange. As a test application, the free energy surface of the folding process of a miniprotein, chignolin, was calculated in the continuum solvent model. Results agreed with the free energy surface derived from the multicanonical simulation. Significantly improved scalability with the MSES method was clearly shown in the free energy calculation of chignolin in explicit solvent, which was achieved without increasing the number of replicas in the Hamiltonian exchange.     

Dynamics versus replicas in the random field Ising model

In a previous article we have shown, within the replica formalism, that the conventional picture of the random field Ising model breaks down, due to the effect of singularities in the interactions between fields involving several replicas below dimension eight. In the zero-replica limit, several coupling constants have thus to be considered, instead of just one. As a result we found that there is no stable fixed point in the vicinity of dimension six. It is natural to reconsider the problem in a dynamical framework, which does not require replicas, although the equilibrium properties should be recovered in the large time limit. Singularities in the zero-replica limit are a priori not visible in a dynamical picture. In this note we show that in fact new interactions are also generated in the stochastic approach. Similarly these interactions are found to be singular below dimension eight. These critical singularities require the introduction of a time origin t₀ at which initial data are given. The dynamical properties are thus dependent upon the waiting time. It is shown here that one can indeed find a complete correspondence between the equilibrium singularities in the limit at n = 0, and the singularities in the dynamics when the initial time t₀ goes to minus infinity, with n replaced by −1/t₀. There is thus complete coherence between the two approaches.