|
An introduction to replica exchange simulations: Mark Abraham, Session 1BTable of contentsOverviewThis tutorial comes in three parts. In the first, the theory behind replica exchange simulations will be briefly described. Then, we will look at how to perform temperature REMD on a small peptide in explicit solvent. Finally, we will consider the related method of Hamiltonian replica exchange, and in particular the technique of Replica-Exchange Solute Tempering. In the directory within the tarball for this tutorial, you will find two subdirectories, each with some files you can use to follow along with the tutorial. If you make a mistake, there's backups of input and output for each stage in the archive subdirectory for each stage. This tutorial assumes you are comfortable with basic GROMACS and UNIX usage, including using The text of the tutorial is available in PDF format. This gzip-ed tar archive of the input files will also be needed to work along with the tutorial. Theory of replica exchange simulationsMany molecular simulation scenarios require ergodic sampling of energy landscapes that feature many minima, and barriers between minima that can be difficult to cross at ambient temperatures over accessible simulation time scales. This means that results are confounded by the choice of initial conditions, because these determine the region of space that is explored by the simulation. Replica exchange simulations seek to enhance the sampling in such scenarios by running numerous independent replicas in slightly different ensembles, and periodically exchanging the coordinates of replicas between the ensembles. Typically, the set of replicas is constructed so that one extreme of the set of replicas is the ensemble from which sampling is wanted, and the other extreme is one where barriers can be crossed more easily. So, if
then the resulting sampling will be statistically correct and closer to the ergodic ideal. How can this work? It relies on doing Monte Carlo moves of replicas between ensembles. The probability of observing a replica in a particular state depends on the potential energy and the temperature, i.e. $P(x) \propto \exp{-\frac{U(x)}{k_B T}}$. The range of probabilities that are observed follow a normal distribution (because of the Central Limit Theorem). If two ensembles are chosen so that their distribution of states have significant overlap, then when we observe a state in one there is appreciable chance that it could have been observed in the other. ![]() Potential energy distributions of alanine dipeptide replicas. This allows us to construct a Metropolis-style Monte Carlo move that proposes that we exchange the coordinates found in two replicas based upon the probabilities that they would have been observed in the two ensembles. The exchange is accepted only if a random number has a suitable value. When done correctly, detailed balance is achieved, and the sampling in both ensembles is correct whether or not an exchange took place. What range of ensembles should be chosen? The method of REMD (also known as parallel tempering) conducts the same simulation over a range of temperatures. The highest temperature is chosen so that the barriers will be crossed. Replicas wander up and down through temperature space as exchanges are accepted. Barriers are crossed probabilistically at all of the temperatures, but at a higher rate at the higher temperatures. Those states filter down to the lower temperatures if they "belong" to the corresponding ensembles. In this way, the ergodicity of the sampling is enhanced. For more details, see (many) papers by Sugita Okamoto and Ulrich Hansmann. Stage 1Planning an REMD simulationThere are a number of inter-connected issues to consider in designing an REMD simulation:
In this tutorial, these decisions are made arbitrarily, but for real studies you will want to consult the literature for guidance here. Suggested reading includes papers by David Kofke, Ulrich Hansmann, and Mark Abraham. Temperatures should often be distributed across the replicas in a geometric progression, i.e. with the same ratio used to scale each temperature from the one below it. If the energy landscapes are sufficiently similar across the temperature space, then this will keep the overlap of the potential energy distributions constant. In turn, this will keep the probability of accepting exchanges constant. Beware, though - in real simulations you will find bottle-necks in temperature space that restrict replica flow. You should be prepared to revise your temperature scheme after some testing. As you can read in the above literature, an exchange acceptance probability around 0.2 is generally a good idea. Often you will have to experiment with the number of replicas you need to use to span the desired temperature range to achieve about that exchange probability. If you are simulating in the NPT ensemble, you might try this REMD temperature generator to help out. Here, we will use only four replicas, so that the tutorial can run reasonably. Beware that you might want many hundreds of replicas to span a decent temperature range if you were looking at something like a large protein! The temperatures we will use are 298.00, 308.00, 318.34, and 329.02. Running REMD equilibrations in GROMACSNormally, mdrun runs a single simulation. For temperature REMD we need to run a number of simulations that can communicate. This is done via the You will find directories
or if you want to check the reference temperatures for the thermostats in each, you can do that efficiently with a command like
Now we want to run grompp in each directory to build the independent
This says to loop over
you will see that there is a
to run the simultion. Use The Stage 2Running a REMD simulations in GROMACSIn this stage, we'll take the output from an equilibration run of decent length and use that as input for the REMD simulation proper. Change to the the
where you can see that we are running a continuation, rather than generating new velocities! You can also see that the only differences between the simulation
Now let's use grompp to prepare the input
This time we told
This run will take a bit longer than the equilibration run, but is still only a "toy" run. We triggered the use of REMD with the So, what happened? Let's look at one of the four
and search through it for the output peculiar to replica exchange. Fortunately, all of that is prefixed with the same string, so if you use You can see that
Here we see that at step 100 (which was the first attempt), the time was 0.2 ps. The potential energy of the state in ensemble 1 was about 3.5 kT higher than that of the state in ensemble 0, which meant that an exchange would only succeed about 3% of the time. Since this is the The next exchange attempt at step 200 in my run was a bit different (yours will be different again)
At replica-exchange steps, GROMACS alternates between two disjoint sets of replica pairs, so that each replica attempts an exchange with both of its neighbours once every two attempts. This time, replica 2 had a potential energy that was lower than replica 1, so the exchange is certain to succeed. Do read up on the Metropolis criterion if you need more information here. Down near the bottom of the The REMD statistics are accumulated only over each run, and are not stored in the This REMD simulation was also far too short to show anything useful. In the next section will take a look at what you might see in a longer run. De-multiplexing an REMD trajectory in GROMACSThere are two ways the programmer could implement REMD. At exchange attempts the replica simulations could exchange coordinates, or they could exchange temperatures. Other MD packages tend to exchange temperatures; GROMACS exchanges coordinates. This means that the trajectory written to a single file belongs to the ensemble described by the matching There are some times you want a trajectory that has continuous coordinates, despite the "jumps" in ensemble space. This means chopping up each trajectory file into pieces and pasting them back together in the right order. The Perl script ConclusionHere we've scratched the surface of the possibilities of replica-exchange. There are lots of fancy ways to change the physics of the replicas from which you don't want to sample, in order to enhance the sampling at the one you do want to sample, including actually changing the Hamiltonian! Author, feedback, and getting helpThis tutorial was prepared by Mark Abraham for the 2013 GROMACS USA Workshop. Feedback on the format and content is warmly invited. General MD and GROMACS questions arising from the above discussion should go to the GROMACS users mailing list. I am unable to give private help to individuals as they work their way through this tutorial, nor with general molecular simulation tasks. |