Aims

The aim of this tutorial is to show how to use PLUMED to run simulations with multiple replicas and how to analyze them. In particular, we will focus on cases where the replicas feel different biasing potentials.

Objectives

Once this tutorial is completed students will be able to:

Write plumed input files suitable for multi-replica simulations.
Run replica exchange simulations with gromacs and plumed using different bias potentials in each replica.
Analyze replica exchange simulations using WHAM so as to obtain the weight of each snapshot.

Resources

The TARBALL for this project contains the following files:

topol0.tpr, topol1.tpr, topol2.tpr, and topol3.tpr, gromacs input files for analine dipeptide. Notice that two of them (0 and 2) are initialized in one free-energy well and two of them (1 and 3) in the other free-energy well.
A directory SCRIPT that contains a wham.py script to perform WHAM. Notice that this required python 3 (does not work with python 2).
A directory SETUP with the gromacs input file that can be used to generate new tpr files.

This tutorial has been tested on a pre-release version of version 2.4. In particular, it takes advantage of a special syntax for setting up multi-replica simulations that is only available since version 2.4. Exercise could be done also with version 2.3 but using a different syntax with respect to the one suggested.

Also notice that in the .solutions directory of the tarball you will find correct input files. Please only look at these files after you have tried to solve the problems yourself.

Introduction

So far we always used PLUMED to run one simulation at a time. However, PLUMED can also be used in multi-replica algorithms. When coupled with GROMACS (at least) it is also possible to run replica exchange simulations, where coordinates are swapped from time to time. In this case, PLUMED is going to take into account the different bias potentials applied to different replicas so as to compute correctly the acceptance.

Similarly to what we did before, we will first use the driver to understand how to prepare multi-replica input files. However, the very same holds when you run multi-replica MD simulations with MD codes that support them. For instance, in GROMACS that would be using the -multi option of mdrun.

Notice that this tutorial was tested using a pre-release version of PLUMED 2.4. In particular, we will use a special syntax for multi-replica input that is only available starting with PLUMED 2.4.

Multi replica input files

Imagine that you are in a directory with these files

traj.0.xtc
traj.1.xtc
traj.2.xtc
plumed.dat

That is: three trajectories and a PLUMED input files. Let's say that the content of plumed.dat is

Click on the labels of the actions for more information on what each action computes

d: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label d calculates a single scalar value
PRINT ARGthe input for this action is the scalar output from one or more other actions. 
=d FILEthe name of the file on which to output these quantities 
=COLVAR The PRINT action with label

You can use the following command to process the three trajectories simultaneously:

mpirun -np 3 plumed driver --multi 3 --plumed plumed.dat --mf_xtc traj.xtc

This command will produce three COLVAR files, namely COLVAR.0, COLVAR.1, and COLVAR.2.

How does this work?

You are here running three replicas of the plumed executable. When PLUMED opens a file for writing it adds a a suffix corresponding to the replica number. This is done for all the possible files written by PLUMED and allows you to distinguish the output of different replicas redirecting it to different files.

Attention: This is true also for input files that are opened for reading. However, when PLUMED does not find the input file with the replica suffix, it looks for the file without the suffix. That's why here it is sufficient to provide a single plumed.dat file. If by mistake you include an additional plumed.0.dat file in the same directory, PLUMED will use that file for replica 0 (and plumed.dat for replicas 1 and 2). To be precise, this is true for all the files read by PLUMED, with the exception of PDB files where the suffix is never added.

The last comment implies that if you need to process your trajectories with different input files you might do it creating files named plumed.0.dat, plumed.1.dat, and plumed.2.dat. This is correct, and this was the only way to do multi-replica simulations with different input files up to PLUMED 2.3. However, in the following we will see how to obtain the same effect with a new special command that has been introduced in PLUMED 2.4.

Using special syntax for multiple replicas

In many cases, we need to run multiple replicas with almost identical PLUMED files. These files might be prepared with cut-and-paste, which is very error prone, or could be set up with some smart bash or python script. Additionally, one can take advantage of the INCLUDE keyword so as to have a shared input file with common definitions and specific input files with replica-dependent keywords. However, as of PLUMED 2.4, we introduced a simpler manner to manipulate multiple replica inputs with tiny differences. Look at the following example:

Click on the labels of the actions for more information on what each action computes

If you prepare a single plumed.dat file like this one and feeds it to PLUMED while using 3 replicas, the 3 replicas will see the very same input except for the AT keyword, that sets the position of the restraint. Replica 0 will see a restraint centered at 1.0, replica 1 centered at 1.1, and replica 2 centered at 1.2.

The @replicas: keyword is not special for RESTRAINT or for the AT keyword. Any keyword in PLUMED can accept that syntax. For instance, the following single input file can be used to setup a bias exchange metadynamics [piana] simulations:

Click on the labels of the actions for more information on what each action computes

#SETTINGS NREPLICAS=2
# Compute distance between atoms 1 and 2
d: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label d calculates a single scalar value
# Compute a torsional angle
t: TORSION ATOMSthe four atoms involved in the torsional angle 
=30,31,32,33 The TORSION action with label t calculates a single scalar value
# Metadynamics.
METAD ...
   ARGthe input for this action is the scalar output from one or more other actions. 
=@replicas:d,t 
   HEIGHTthe heights of the Gaussian hills. =1.0 
   PACEcompulsory keyword 
the frequency for hill addition =100 
   SIGMAcompulsory keyword 
the widths of the Gaussian hills =@replicas:0.1,0.3 
   GRID_MINthe lower bounds for the grid =@replicas:0.0,-pi 
   GRID_MAXthe upper bounds for the grid =@replicas:2.0,pi 
...The METAD action with label 
# On replica 0, this means:
# METAD ARG=d HEIGHT=1.0 PACE=100 SIGMA=0.1 GRID_MIN=0.0 GRID_MAX=2.0
# On replica 1, this means:
# METAD ARG=t HEIGHT=1.0 PACE=100 SIGMA=0.3 GRID_MIN=-pi GRID_MAX=pi

This would be a typical setup for a bias exchange simulation. Notice that even though variables d and t are both read in both replicas, d is only computed on replica 0 (and t is only computed on replica 1). This is because variables that are defined but not used are never actually calculated by PLUMED.

If the value that should be provided for each replica is a vector, you should use curly braces as delimiters. For instance, if the restraint acts on two variables, you can use the following input:

Click on the labels of the actions for more information on what each action computes

#SETTINGS NREPLICAS=3
# Compute distance between atoms 1 and 2
d: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=10,20 The DISTANCE action with label d calculates a single scalar value
# Compute a torsional angle
t: TORSION ATOMSthe four atoms involved in the torsional angle 
=30,31,32,33 The TORSION action with label t calculates a single scalar value
# Apply a restraint:
RESTRAINT ...
   ARGthe input for this action is the scalar output from one or more other actions. 
=d,t 
   ATcompulsory keyword 
the position of the restraint =@replicas:{{1.0,2.0} {3.0,4.0} {5.0,6.0}}  
   KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =1.0,3.0 
...The RESTRAINT action with label 
# On replica 0 this means:
# RESTRAINT ARG=d AT=1.0,2.0 KAPPA=1.0,3.0
# On replica 1 this means:
# RESTRAINT ARG=d AT=3.0,4.0 KAPPA=1.0,3.0
# On replica 2 this means:
# RESTRAINT ARG=d AT=5.0,6.0 KAPPA=1.0,3.0

Notice the double curly braces. The outer ones are used by PLUMED to know there the argument of the AT keyword ends, whereas the inner ones are used to group the values corresponding to each replica. Also notice that the last example can be split in multiple lines exploiting the fact that within multi-line statements (enclosed by pairs of ...) newlines are replaced with simple spaces:

Click on the labels of the actions for more information on what each action computes

In short, whenever there are keywords that should vary across replicas, you should set them using the @replicas: keyword. As mentioned above, you can always use the old syntax with separate input file, and this is recommended when the number of keywords that are different is large.

Exercise 1: Running multi-replica simulations

Write a plumed file that allows you to run a multi-replica simulation of alanine dipeptide where the following four replicas are simulated:

Two replicas with no bias potential.
A replica with metadynamics on \(\phi\).
A replica with metadynamics on \(\psi\).

With PLUMED 2.3 you should use four plumed.dat files (named plumed.0.dat, plumed.1.dat, etc). In this case, the first two replicas feel no potential, so you could just use empty files. For the third and fourth replica you could recycle the files from Trieste tutorial: Metadynamics simulations with PLUMED.

However, we recommend you t