Temperature replica exchange

From CHARMM tutorial
Jump to: navigation, search

Contents

Introduction

Temperature replica exchange molecular dynamics (T-REM), also known as parallel tempering, is an enhanced sampling method that works by simulating (via MD molecular dynamics replicas of a system at a range of different temperatures and periodically exchanging between them, as is shown in the figure below.

Note: figure will go here

There are two key points to understand about this method. Firstly, if implemented correctly (which CHARMM does) and a good thermostat is used, it preserves the canonical ensemble. Secondly, replicas run at a higher temperature will explore conformational space more quickly than those which run at physiological temperature (around 300K). If you need to know how to use tempoerature control in CHARMM, please refer to the MD page of this tutorial.

The method is fully described by Sugits and Okamoto. Chem. Phys. Lett. 1999, pp. 141-151.

How temperature replica exchange works

Temperature replica exchange is based on sound principles of statistical mechanics that are explained thoroughly in the Sugita and Okamoto paper referenced above (as well as many other sources). The following discussion is rather informal,

One key point is that statistical ensembles can be combined to form a generalized ensemble. For example, consider the case where we have N non-interacting ensembles 1,...,n that are all canonical ensembles at temperatures T1,...,Tn. We know that the probablity of a conformation X existing in ensemble y is given by:

P(X,y) = \frac{1}{Z}e^{-\beta_{y}E(X)}

where

βy = (kBTy) − 1, kB is Boltzmann's constant, and Z is the canonical partition function.

The question now can be posed: if we have a generalized ensemble of n replicas as described above, what is the probability of any particular set of coordination Xg = {X1,...,Xn} existing in that ensemble?

The answer is straightforward; since replicas are not interacting, then the probability is simply the product of the probabilities of the constituent ensemble probabilities:

P(X_{g}) = \prod_{i=1}^{N}P(X_{i},y_{i}) = \prod_{i=1}^{N}\frac{1}{Z_{i}}e^{-\beta{i}E(X_{i})} \propto \exp\{-\sum_{i=1}^{N}\beta_{i}E(X_{i})\}

where i is an index running over all N replicas.

Preserving detailed balance in the generalized ensemble

Now that we've defined our generalized ensemble, we can ask if there's a way to exchange coordinates between two replicas (ensembles) at different temperatures. To put it more formally, suppose we have two replicas: j and k. Is it possible, to exchange their coordinates X_{j} \leftrightarrow X_{k} while still preserving the canonical ensemble? Note that exchanging coordinates is equivalent to exchanging temperatures.

It turns out that the answer to the question in the previous paragraph is yes, so long as the exchange is done in accordance to a particular criteria. We cannot make these exchanges as we wish, but instead we must exchange based on some exchange probability that preserves the statistical ensemble. In particular, detailed balance must be preserved. Detailed balance means that the probability of accepting a proposed trial move (exchange) must be the same as the probability of accepting the reverse move. In the replica exchange case, that if we decide we can swap replicas with some probability p, we must be able to swap them back with the same probability.

We also have to consider what happens with the velocities of the particles when swapping. We can choose a rule such that we scale the velocities of the replicas by the square root of the ratio of the temperatures if we accept the swap. In other words v_{j} = v_{j}\sqrt{T_{k}/T_{j}} and v_{k} = v_{k}\sqrt{T_{j}/T_{k}}. For more details of why we use this rule, please see the Sugita and Okamoto paper referenced above. Exchange criteria

From the previous section, we know that the probability we give to performing an exchange must be equal to the probability of performing the reverse exchange. Let's write that out more formally. We'll define the state where the replica with coordinates and velocities j is at temperature Tj and the replica with coordinates and velocities k is at temperature Tk state X and the state where the temperatures (equivalently coordinates and velocities) have been swapped as state X'.

We can write the ratio of the probability weights of the two possible swaps (X \rightarrow X' and X' \rightarrow X) using the probabilities we defined above:

<math<\frac{w(X \rightarrow X')}{w(X' \rightarrow X)} = \exp\{(-\beta_{T_{j}}E(k)-\beta_{T_{k}}E(j))-(-\beta_{T_{j}}E(j)-\beta_{T_{k}}E(k))\} = \exp\{-(\beta_{T_{j}}-\beta_{T_{i}})[E(k)-E(j)]\}</math>

Therefore, to preserve detailed balance, the exchange criteria is:

w(X \rightarrow X') = \min\{1.0,\exp\{-(\beta_{T_{j}}-\beta_{T_{i}})[E(k)-E(j)]\}\}

The complex expression \exp\{-(\beta_{T_{j}}-\beta_{T_{i}})[E(k)-E(j)]\} is commonly written as simply Δ.

In practice, code implementing replica exchange usually only considers exchanges between neighboring replicas (although the method is not limited in this way), and operates by simply calculating Δ, choosing a uniformly-distributed random number between 0 and 1, and proceeding with the exchange if the random number is less than or equal to Δ.

CHARMM example

Replica exchange is performed in CHARMM via the REPD command. In order to use this functionality, a parallel version of CHARMM must have been built with the REPDSTR key-word in pref.dat.

When the REPD command is invoked, several pieces of information must be given. The type of exchange must be specified (T-REM is given by the EXCH key-word; other types will be discussed briefly below). Also the number of desired replicas (NREP) and exchange frequency (FREQ) should be given. Optional, but highly recommended, is the ability to specify a unit number to write the exchange log for each replica. After these options are given, the range of temperatures may be specified by multiple invocations of the TEMP subcommand. A valid invocation of the REPD command might look like:

REPD EXCH NREP 4 UNIT 20 FREQ 1000 -
   TEMP 300 TEMP 330 TEMP 370 TEMP 420

This would tell CHARMM to simulate four replicas at 300, 330, 370, and 420 kelvins, while attempting exchanges between them every 1000 molecular dynamics steps.

After this command is processed, CHARMM would run 4 separate simulations, dividing up the number of total processors evenly between them. Note: the number of replicas must evenly divide the number of processors.

I/O changes with replica exchange

Because after the REPD command CHARMM is running many simulations instead of just one, there are some necessary changes to the way I/O is handled. Each replica is automatically given its own input and output file set, which can lead to some conclusion. For example, if, after issuing the REPD command given above, you were to issue the command:

OPEN UNIT 50 WRITE CARD NAME MYFILE.DAT

You would actually get 4 new files: myfile.dat_0, myfile.dat_1, myfile.dat_2, and myfile.dat_3. Each of these files would have been opened by the replica whose number is appended to it after the underscore (note: CHARMM numbers replicas starting from 0 instead of 1).

The same process holds with reading. If you were to issue:

OPEN UNIT 51 READ CARD NAME BLAH.DAT

then replica 0 would try to open blah.dat_0, replica 1 would try to open blah.dat_1, etc. If these files don't all exist, CHARMM will die. If you need all replicas to read a particular data file, it is common to read it in before the REPD command is issued or, if that is impossible, make as many copies (with the correct suffixes) as needed that are all symbolic linkls to the data file.

Note that CHARMM does not automatically split up its output by replica, but you can use these I/O features along with the output redirection command (OUTUnit) to accomplish this yourself as shown in the following example:

REPD EXCH NREP 4 UNIT 20 FREQ 1000 -
   TEMP 300 TEMP 330 TEMP 370 TEMP 420

OPEN UNIT 30 WRITE CARD NAME REPLICA.OUT
OUTU 30

OPEN UNIT 20 WRITE CARD NAME REPLICA.EXCH

In this case, the output from each replica will go to a file called replica.out_X and its exchange log will be placed in replica.exch_X (where X is the number of the replica).

Complete example input script

With all this in mind, let's look at a complete example CHARMM script that does temperature replica exchange:

read rtf card name top_all36_prot.rtf
read parameter card name par_all36_prot.prm

read psf card name struct.psf
read coor card name struct.crd

repd nrep 4 exch unit 40 freq 500 -
  temp 300 temp 330 temp 370 temp 420

open write unit 41 card name replica.out
outu 41

open write unit 40 card name replica.exch     ! exchange log file
open write unit 31 card name replica.rst      ! restart file for each replica
open write unit 32 unfo name replica.dcd      ! trajectory file for each replica

energy

scalar fbeta set 0.5 

dyna leap lang start nstep 10000 timestep 0.001 nprint 1000 -
  nsavc 1000 iuncrd 32 iunwri 31 iprfrq 1000 -
  firstt 300. finalt 300. tstruct 300. tbath 300. -
  iasors 1 iasvel 1 iscvel 0 ichecw 0 

stop

There are a couple of notes about this:

  • Before the REPD command, the topology, parameter, PSF, and coordinates are read as usual.
  • The output file is opened after the REPD command and each replica immediately switches to it via the OUTU command.
  • Separate restart and trajectory files are produced for each replica.
  • It is unnecessary to adjust the FIRSTT, FINALT, TSTRUCT, or TBATH parameters to the DYNAmics command; CHARMM does this automatically.
  • For more information on the dynamics set-up, please refer to the MD and LD sections of this tutorial.

Interpreting the exchange log

Let's take a short look at a section of the exchange log produced by CHARMM. Remember, if requested, one such log file is produced for each replica in the simulation. Below is an excerpt of such a long from a real simulation:

------------- Replica Exchange ------------
REX>EXCHANGE =          8  Step =      4000
REX>REPL     =     2  Temp = 340.000  Epot =         186.96850203
REX>NEIGHBOR =     3  Temp = 360.000  Epot =         197.43771759
REX>PROB     =  0.42281 Rand =  0.61274 Tscale =  1.0000 Success = F
------------- Replica Exchange End --------
------------- Replica Exchange ------------
REX>EXCHANGE =          9  Step =      4500
REX>REPL     =     2  Temp = 340.000  Epot =         169.55711821
REX>NEIGHBOR =     1  Temp = 320.000  Epot =         169.65014854
REX>PROB     =  1.00000 Rand =  0.50920 Tscale =  1.0308 Success = T
------------- Replica Exchange End --------

This log file is for replica 2 (the 3rd replica in the system). It has a temperature of 340 K. On the 8th exchange (which occurred after the 4000th dynamics step), it attempted an exchange with replica 3 (NEIGHBOR), which had a temperature 20 K higher - 360 K. The potential energies of each replica (~ 186.97 for replica 2 and ~ 197.44 for replica 3). Based on the exchange criterion given above, an exchange should occur about 42% of the time (PROB). In this case, however, the random number (Rand) was ~ 0.61, so the exchange failed.

The 9th exchange is interesting because the lower temperature replica had the higher potential energy. When this happens, the exchange will always succeed (as an exercise, prove that this is so). Therefore, the probability is 1.0.

Other types of replica exchange

Replica exchange can occur in principle along any variable, as long as a well-defined statistical ensemble can be calculated based on its changes. The most recent versions of CHARMM also support Hamiltonian, SGLD, and pH replica exchange. Interested readers are encouraged to consult repdstr.doc in their version of CHARMM to learn about what is available. Various papers have been published on these methods and references are available in repdstr.doc and by performing a literature search.

Personal tools
Navigation