diff --git a/user-doc/tutorials/a-trieste-5.txt b/user-doc/tutorials/a-trieste-5.txt new file mode 100644 index 0000000000000000000000000000000000000000..9d9f1775df5a49d70828b5a947c132f2aa84e030 --- /dev/null +++ b/user-doc/tutorials/a-trieste-5.txt @@ -0,0 +1,362 @@ +/** +\page trieste-5 Trieste tutorial: Running and analyzing multi-replica simulations. + +\section trieste-5-aims 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. + +\section trieste-5-lo Objectives + +Once this tutorial is completed students will: + +- \highlight{TODO} + +\section trieste-5-resources Resources + +The \tarball{trieste-5} for this project contains the following files: +- \highlight{TODO} + +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. Exercizes could be done also with version 2.3 but using a different syntax +with respect to the one suggested. + +\section trieste-5-introduction 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 \ref 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. + +\section trieste-5-replicas Running multi replica simulations + +Imagine that you are in a directory with these files +\verbatim +traj.0.xtc +traj.1.xtc +traj.2.xtc +plumed.dat +\endverbatim +That is: three trajectories and a PLUMED input files. +Let's say that the content of `plumed.dat` is +\plumedfile +d: DISTANCE ATOMS=1,2 +PRINT ARG=d FILE=COLVAR +\endplumedfile + +You can use the following command to process the three trajectories simultaneously: +\verbatim +mpirun -np 3 plumed driver --multi 3 --plumed plumed.dat --mf_xtc traj.xtc +\endverbatim +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 ony way to do multi-replica +simulatons 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. + +\section trieste-5-replica-special-syntax 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 \ref 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: + +\plumedfile +# Compute a distance +d: DISTANCE ATOMS=1,2 + +# Apply a restraint. +RESTRAINT ARG=d AT=@replicas:1.0,1.1,1.2 KAPPA=1.0 +# On replica 0, this means: +# RESTRAINT ARG=d AT=1.0 KAPPA=1.0 +# On replica 1, this means: +# RESTRAINT ARG=d AT=1.1 KAPPA=1.0 +# On replica 2, this means: +# RESTRAINT ARG=d AT=1.2 KAPPA=1.0 +\endplumedfile + +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 \ref 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 simulations: +\plumedfile +# Compute distance between atoms 1 and 2 +d: DISTANCE ATOMS=1,2 + +# Compute a torsional angle +t: TORSION ATOMS=30,31,32,33 + +# Metadynamics. +METAD ... + ARG=@replicas:d,t + HEIGHT=1.0 + PACE=100 + SIGMA=@replicas:0.1,0.3 + GRID_MIN=@replicas:0.0,-pi + GRID_MAX=@replicas:2.0,+pi +... +# 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 +\endplumedfile + +This would be a typical setup for a bias exchange \highlight{REF} 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: + +\plumedfile +# Compute distance between atoms 1 and 2 +d: DISTANCE ATOMS=10,20 + +# Compute a torsional angle +t: TORSION ATOMS=30,31,32,33 + +# Apply a restraint: +RESTRAINT ... + ARG=d,t + AT=@replicas:{{1.0,2.0} {3.0,4.0} {5.0,6.0}} + KAPPA=1.0,3.0 +... +# 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 +\endplumedfile + +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: +\plumedfile +d: DISTANCE ATOMS=10,20 +t: TORSION ATOMS=30,31,32,33 +RESTRAINT ... + ARG=d,t +# indentation is not required (this is not python!) +# but makes the input easier to read + AT=@replicas:{ + {1.0,2.0} + {3.0,4.0} + {5.0,6.0} + } + KAPPA=1.0 +... +\endplumedfile + +In short, whenever there are keywords that should vary across replicas, you should set them usign 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. + +\section trieste-5-ex-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 \f$\phi\f$. +- A replica with metadynamics on \f$\psi\f$. + +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 \ref trieste-4. + +However, we recommend you to take the occasion to make practice with the special replica syntax +of PLUMED 2.4. +Use the following input file as a starting point +\plumedfile +# load the pdb to better select torsions later +MOLINFO __FILL__ +# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C +phi: TORSION ATOMS=__FILL__ +# Compute the backbone dihedral angle psi, defined by atoms N-CA-C-N +psi: TORSION ATOMS=__FILL__ + +metad: METAD ... +# Activate well-tempered metadynamics +# Deposit a Gaussian every 500 time steps, with initial height equal +# to 1.2 kJoule/mol, biasfactor equal to 10.0 +# Remember: replica 0 and 1: no bias +# replica 2, bias on phi +# replica 3, bias on psi + ARG=__FILL__ + HEIGHT=__FILL__ # make sure that replicas 0 and 1 feel no potential! + PACE=500 + BIASFACTOR=10.0 + SIGMA=0.35 +# Gaussians will be written to file and also stored on grid + FILE=HILLS GRID_MIN=-pi GRID_MAX=pi +... + +# Print both collective variables and the value of the bias potential on COLVAR file +PRINT ARG=__FILL__ FILE=COLVAR STRIDE=10 +\endplumedfile + +Now check the resources provided with this tutorial. +There are 4 `tpr` files available (`topol0.tpr`, `topol1.tpr`, `topol2.tpr`, and `topol3.tpr`). +You can run a replica exchange simulation with gromacs using the following command +\verbatim +> mpirun -np 4 gmx_mpi mdrun -plumed plumed.dat -nsteps 1000000 -replex 120 +\endverbatim +Notice that coordinates swaps will be attempted between different replicas every 120 steps. +When computing the acceptance, gromacs will take into account that different replicas might +be subject to different bias potential. In this example, replicas 0 and 1 are actually subject +to the same potential, so all the exchanges will be accepted. + +\comment +The setup above is very close to the one used in bias exchange simulations. +However, in bias exchange simulations usually one would use at most one neutral (non-biased) replica. +In addition, the \ref RANDOM_EXCHANGE command is often used. + +\section trieste-5-ex-2 Analyzing a multiple-restraint simulation + +In the following we will analyze the result of the simulation above. +To this aim we will have to use the WHAM method. The WHAM procedure described here is +binless and allows one to reweight arbitrary bias potentials without the need to discretize +collective variables. As a first step it is necessary to create a concatenated trajectory +that includes all the conformations that you produced during your simulations. + +\verbatim +> gmx_mpi trjcat -f traj_comp*.xtc -cat -o fulltraj.xtc +\endverbatim + +Now we will process that directory computing the bias potentials associated with +the replicas that we simulated. Similarly to what we did in \ref trieste-4, +we will assume that reweighting is done with final bias potential \cite Branduardi:2012dl. +Notice that also here it would be possible to use \cite Tiwary_jp504920s instead. + +\plumedfile +# this will be file plumed_reweight.dat +# Remember how you did a reweighting for a single metadynamics simulation +# First put the restart keyword: +RESTART +# include here exactly the same file you used to run the simulation +# however, you should make some change to METAD: +# - replace PACE with a huge number +# - add TEMP=300 to set the temperature of the simulation +# (normally this is inherited from the MD code, but driver has no idea about it) +__FILL__ +__FILL__ +__FILL__ +__FILL__ +# you should also change the PRINT command so as not to overwrite +# the files you wrote before. +# here omit STRIDE, so that you will print +# for all the frames in your concatenated trajectory +PRINT ARG=phi,psi FILE=COLVAR_CONCAT + +# Then write the bias potential (as we did for reweighting) +# As a good practice, remember that there might be multiple bias +# potential applied (e.g. METAD and a RESTRAINT). +# to be sure that you include all of them, just combine them: +bias: COMBINE ARG=*.bias +# here omit STRIDE, so that you will print the bias +# for all the frames in your concatenated trajectory +PRINT FILE=COLVAR_WHAM ARG=bias +\endplumedfile + +Then analyze the concatenated trajectory with this command +\verbatim +> mpirun -np 4 plumed driver --mf_xtc fulltraj.xtc --plumed plumed_reweight.dat --multi 4 +\endverbatim + +Notice that since `fulltraj.xtc` the four concurrent copies of PLUMED will all analyze the same trajectory, +which is what we want to do now. +The result will be a number of COLVAR_WHAM files, one per replica, with a number +of lines corresponding to the whole concatenated trajectory. + +\verbatim +# put the bias columns in a single file ALLBIAS +> paste COLVAR_WHAM.* | grep -v \# | awk '{for(i=1;i<=NF/2;i++) printf($(i*2)" ");printf("\n");}' > ALLBIAS +\endverbatim + +Have a look at `ALLBIAS`. It should look more or less like this +\verbatim +0.000000 0.000000 75.636305 75.786147 +0.000000 0.000000 75.143215 75.127876 +0.000000 0.000000 66.071472 65.546944 +0.000000 0.000000 69.774019 69.260869 +0.000000 0.000000 64.679610 64.169994 +0.000000 0.000000 66.843788 66.312660 +0.000000 0.000000 68.703617 68.170744 +0.000000 0.000000 75.142493 75.174527 +\endverbatim +The first two columns correspond to the bias felt in replicas 0 and 1, that are unbiased, and thus should be zero. +In the other replicas on the other hand we will have large potentials accumulated by metadynamics. +The actual numbers might very depending on the length of your simulation. + +Now you can use the provided python script to analyze the `ALLBIAS` file. +`SCRIPTS/wham.py` should be provided three arguments: name of the file with bias potentials, +number of replicas (that is the number of columns) and value of kbT (2.5 in kj/mol). + +\verbatim +# use the python script to compute weights +> python SCRIPTS/wham.py ALLBIAS 4 2.5 +\endverbatim + +This will generate a file named `weights.dat` that contains the weights for each of the frame. + +\highlight{TODO: show and interpret weights} + +\highlight{TODO: compute free-energy difference} + +\highlight{TODO: repeat using only the two unbiased replicas and show that free-energy difference is wrong} + +\hidden{Analyze a solute tempering-METAD simulation} + +\highlight{TODO: prepare input files for hrex simulation} + +\endhidden + +\section trieste-5-ex-3 Computing statistical error + +\highlight{TODO: show how to use block analysis. Also discuss the problems with demuxing trajectories.} + +*/ + +link: @subpage trieste-5 + +description: This tutorial explains how to use PLUMED to analyze trajectories + +additional-files: a-trieste-5