Trieste tutorial: Analyzing trajectories using PLUMED

Aims

The aim of this tutorial is to introduce the users to the PLUMED syntax. We will go through the writing of simple collective variable and we will use them to analyze existing trajectories.

Objectives

Once this tutorial is completed students will be able to:

  • Write a simple PLUMED input file and use it with the PLUMED driver to analyze a trajectory.
  • Use the GROUP keyword to make the input file compact and easy to read and to quickly build complex atom groups.
  • Print collective variables such as distances (DISTANCE), torsional angles (TORSION), gyration radius (GYRATION), and coordination numbers (COORDINATION) using the PRINT action.
  • Computing the geometric center of a group of atoms using CENTER.
  • Know how to take care of periodic boundary conditions within PLUMED using WHOLEMOLECULES and WRAPAROUND, and be able to verify the result with DUMPATOMS.
  • Extract from a trajectory snapshots satisfying specific conditions using UPDATE_IF.

Resources

The TARBALL for this project contains the following files:

  • ref.pdb : A PDB file with a RNA duplex solvated in a water box and a Mg ion.
  • traj-whole.xtc: A trajectory for the same system in GROMACS xtc format. To make the exercise easier, RNA duplex has been made whole already.
  • traj-broken.xtc: The same trajectory as it was originally produced by GROMACS. Here the RNA duplex is broken and should be fixed.

This tutorial has been tested on a pre-release version of version 2.4. However, it should not take advantage of 2.4-only features, thus should also work with version 2.3.

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

This tutorial asks you to compute a variety of different collective variables using PLUMED for a particular trajectory and to compare the files and graphs that you obtain with the correct ones that are shown online. Compared to some of the other tutorials that are available here this tutorial contains considerably less guidance so in doing this tutorial you will have to learn how to consult the manual. If you would like a more guided introduction to PLUMED it might be better to start with the tutorials Belfast tutorial: Analyzing CVs or MARVEL tutorial: Analyzing CVs. Also notice that, whereas this tutorial was tested using a pre-release version of PLUMED 2.4, it should be completely feasible using PLUMED 2.3.

Using PLUMED from the command line

As we will see later, PLUMED provides a library that can be combined with multiple MD codes. However, in this tutorial we will only use PLUMED to analyze trajectories that have been produced already. Once PLUMED is installed you can run a plumed executable that can be used for multiple purposes:

> plumed --help

Here we will use the driver tool, that allows you to process an already existing trajectory.

> plumed driver --help

What we will need is:

  • A trajectory to be analyzed (provided).
  • A PLUMED input file (you do it!).

The syntax of the PLUMED input file is the same that we will use later to run enhanced sampling simulations, so all the things that you will learn now will be useful later when you will run PLUMED coupled to an MD code. In the following we are going to see how to write an input file for PLUMED.

The structure of a PLUMED input file

The main goal of PLUMED is to compute collective variables, which are complex descriptors than can be used to analyze a conformational change or a chemical reaction. This can be done either on the fly, that is during molecular dynamics, or a posteriori, using PLUMED as a post-processing tool. In both cases one should create an input file with a specific PLUMED syntax. A sample input file is below:

Click on the labels of the actions for more information on what each action computes
tested on v2.9
# this is optional and tell to VIM that this is a PLUMED file
# vim: ft=plumed
# see comments just below this input file
# Compute distance between atoms 1 and 10.
# Atoms are ordered as in the trajectory files and their numbering starts from 1.
# The distance is called "d" for future reference.
d: DISTANCE 
ATOMS
the pair of atom that we are calculating the distance between. 
=1,10 
# Create a virtual atom in the center between atoms 20 and 30.
# The virtual atom only exists within PLUMED and is called "center" for future reference.
center: CENTER 
ATOMS
the list of atoms which are involved the virtual atom's definition. 
=20,30 
# Compute the torsional angle between atoms 1, 10, 20, and center.
# Notice that virtual atoms can be used as real atoms here.
# The angle is called "phi" for future reference.
phi: TORSION 
ATOMS
the four atoms involved in the torsional angle 
=1,10,20,center 
# Compute some function of previously computed variables.
# In this case we compute the cosine of angle phi and we call it "d2"
d2: MATHEVAL ...
   
ARG
the input for this action is the scalar output from one or more other actions. 
=phi 
FUNC
compulsory keyword 
the function you wish to evaluate 
=cos(x) 
   
PERIODIC
compulsory keyword 
if the output of your function is periodic then you should specify the periodicity
of the function. 
=NO 
...
# The previous command has been split in multiple lines.
# It could have been equivalently written in a single line:
# d2: MATHEVAL ARG=phi FUNC=cos(x) PERIODIC=NO
# Print d and d2 every 10 step on a file named "COLVAR1".
PRINT 
ARG
the input for this action is the scalar output from one or more other actions. 
=d,d2 
STRIDE
compulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output 
=10 
FILE
the name of the file on which to output these quantities 
=COLVAR1 
# Print phi on another file names "COLVAR2" every 100 steps.
PRINT 
ARG
the input for this action is the scalar output from one or more other actions. 
=phi 
STRIDE
compulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output 
=100 
FILE
the name of the file on which to output these quantities 
=COLVAR2
Note
If you are a VIM user, you might find convenient configuring PLUMED syntax files, see Using VIM syntax file. Syntax highlighting is particularly useful for beginners since it allows you to identify simple mistakes without the need to run PLUMED. In addition, VIM has a full dictionary of available keywords and can help you by autocomplete your commands.

In the input file above, each line defines a so-called action. An action could either compute a distance, or the center between two or more atoms, or print some value on a file. Each action supports a number of keywords, whose value is specified. Action names are highlighted in green and, clicking on them, you can go to the corresponding page in the manual that contains a detailed description for each keyword. Actions that support the keyword STRIDE are those that determine how frequently things are to be done. Notice that the default value for STRIDE is always 1. In the example above, omitting STRIDE keywords the corresponding COLVAR files would have been written for every frame of the analyzed trajectory. All the other actions in the example above do not support the STRIDE keyword and are only calculated when requested. That is, d and d2 will be computed every 10 frames, and phi every 100 frames. In short, you can think that for every snapshot in the trajectory that you are analyzing PLUMED is going to execute all the listed actions, though some of them are optimized out when STRIDE is different from 1.

Also notice that PLUMED works using kJ/nm/ps as energy/length/time units. This can be personalized using UNITS, but we will here stay with default values.

Variables should be given a name (in the example above, d, phi, and d2), which is then used to refer to these variables. Instead of a: DISTANCE ATOMS=1,2 you might equivalently use DISTANCE ATOMS=1,2 LABEL=a. Lists of atoms should be provided as comma separated numbers, with no space. Virtual atoms can be created and assigned a name for later use.

You can find more information on the PLUMED syntax at Getting Started page of the manual. The complete documentation for all the supported collective variables can be found at the Collective Variables page.

To analyze the trajectory provided here, you should:

  • Create a PLUMED input file with a text editor (let us call it plumed.dat) similar to the one above.
  • Run the command plumed driver --mf_xtc traj.xtc --plumed plumed.dat.

Here traj.xtc is the trajectory that you want to analyze. Notice that driver can read multiple file formats using embedded molfile plugins from VMD (that's where the mf letters come from).

Notice that you can also visualize trajectories with VMD directly. Trajectory traj.xtc can be visualized with the command vmd ref.pdb traj-whole.xtc.

In the following we will make practice with computing and printing collective variables.

Exercise 1: Computing and printing collective variables

Analyze the traj-whole.xtc trajectory and produce a colvar file with the following collective variables.

  • The gyration radius of the solute RNA molecule (GYRATION). Look in the ref.pdb file which are the atoms that are part of RNA (search for the first occurrence of a water molecule, residue name SOL). Remember that you don't need to list all the atoms: instead of ATOMS=1,2,3,4,5 you can write ATOMS=1-5.
  • The torsional angle (TORSION) corresponding to the glycosidic chi angle \(\chi\) of the first nucleotide. Since this nucleotide is a purine (guanine), the proper atoms to compute the torsion are O4' C1 N9 C4. Find their serial number in the ref.pdb file or learn how to select a special angle reading the MOLINFO documentation.
  • The total number of contacts (COORDINATION) between all RNA atoms and all water oxygen atoms. For COORDINATION, set reference distance R_0 to 2.5 A (be careful with units!!). Try to be smart in selecting the water oxygen atoms without listing all of them explicitly.
  • Same as before but against water hydrogen. Also in this case you should be smart to select water hydrogen atoms. Documentation of GROUP might help.
  • Distance between the Mg ion and the geometric center of the RNA duplex (use CENTER and DISTANCE).

Notice that some of the atom selections can be made in a easier manner by using the MOLINFO keyword with a proper reference PDB file. Also read carefully the Groups and Virtual Atoms page before starting. Here you can find a sample plumed.dat file that you can use as a template. Whenever you see an highlighted FILL string, this is a string that you should replace.

Click on the labels of the actions for more information on what each action computes
tested on v2.9
# First load information about the molecule.
MOLINFO __FILL__
# Notice that this is special kind of "action" ("setup action")
# that is only used during setup. It will not be re-executed at each step.
# Define some group that will make the rest of the input more readable
# Here are the atoms belonging to RNA.
rna: GROUP 
ATOMS
the numerical indexes for the set of atoms in the group. 
=1-258 
# This is the Mg ion. A group with atom is also useful!
mg: GROUP 
ATOMS
the numerical indexes for the set of atoms in the group. 
=6580 
# This group should contain all the atoms belonging to water molecules.
wat: GROUP 
ATOMS
the numerical indexes for the set of atoms in the group. 
=__FILL__ 
# Select water oxygens only:
owat: GROUP __FILL__
# Select water hydrogens only:
hwat: GROUP __FILL__
# Compute gyration radius:
r: GYRATION 
ATOMS
the group of atoms that you are calculating the Gyration Tensor for. 
=__FILL__ 
# Compute the Chi torsional angle:
c: TORSION 
ATOMS
the four atoms involved in the torsional angle 
=__FILL__ 
# Compute coordination of RNA with water oxygens
co: COORDINATION 
GROUPA
First list of atoms. 
=rna 
GROUPB
Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted). 
=owat 
R_0
 could not find this keyword 
=__FILL__ 
# Compute coordination of RNA with water hydrogens
ch: COORDINATION 
GROUPA
First list of atoms. 
=rna 
GROUPB
Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted). 
=hwat __FILL__
# Compute the geometric center of the RNA molecule:
ce: CENTER 
ATOMS
the list of atoms which are involved the virtual atom's definition. 
=__FILL__ 
# Compute the distance between the Mg ion and the RNA center:
d: DISTANCE 
ATOMS
the pair of atom that we are calculating the distance between. 
=__FILL__ 
# Print the collective variables on COLVAR file
# No STRIDE means "print for every step"
PRINT 
ARG
the input for this action is the scalar output from one or more other actions. 
=r,c,co,ch,d 
FILE
the name of the file on which to output these quantities 
=COLVAR

Once your plumed.dat file is complete, you can use it with the following command 

> plumed driver --plumed plumed.dat --mf_xtc whole.xtc

Scroll in your terminal to read the PLUMED log. As you can see, PLUMED gives a lot of feedback about the input that he is reading. There's the place where you can check if PLUMED understood correctly your input.

The command above will create a file COLVAR like this one: 

#! FIELDS time r c co ch d
#! SET min_c -pi
#! SET max_c pi
 0.000000 0.788694 -2.963150 207.795793 502.027244 0.595611
 1.000000 0.804101 -2.717302 208.021688 499.792595 0.951945
 2.000000 0.788769 -2.939333 208.347867 500.552127 1.014850
 3.000000 0.790232 -2.940726 211.274315 514.749124 1.249502
 4.000000 0.796395 3.050949 212.352810 507.892198 2.270682

Notice that the first line informs you about the content of each column and the second and third lines tell you that variable c (the \(\chi\) torsion) is defined between \(-\pi\) and \(+\pi\).

In case you obtain different numbers, check your input, you might have made some mistake!

This file can then be shown with gnuplot 

gnuplot> p "COLVAR" u 1:2, "" u 1:3

As a final note, look at what happens if you run the exercise twice. The second time, PLUMED will back up the previously produced file so as not to overwrite it. You can also concatenate your files by using the action RESTART at the beginning of your input file.

To learn more: Combining collective variables

Solving periodic-boundary conditions issues

While running PLUMED can also dump the coordinate of the internally stored atoms using DUMPATOMS. This might seem useless (coordinates are already contained in the original trajectory) but can be used in the following cases:

  • To dump coordinates of virtual atoms that only exist within PLUMED (e.g. a CENTER).
  • To dump snapshots of our molecule conditioned to some value of some collective variable (see UPDATE_IF).
  • To dump coordinates of atoms that have been moved by PLUMED.

The last point is perhaps the most surprising one. Some of the PLUMED actions can indeed move the stored atoms to positions better suitable for the calculation of collective variables.

The previous exercise was done on a trajectory where the RNA was already whole. For the next exercise you will use the traj-broken.xtc file instead, which is a real trajectory produced by GROMACS. Open it with VMD to understand what we mean with broken 

> vmd ref.pdb traj-broken.xtc

Select Graphics, then Representations, then type nucleic in the box Selected Atoms. You will see that your RNA duplex is not whole. This is not a problem during MD because of periodic boundary conditions. However, it is difficult to analyze this trajectory. In addition, some collective variables that you might want to compute could require the molecules to be whole (an example of such variables is RMSD).

You might think that there are alternative programs that can be used to reconstruct the molecules the molecules that are broken by the periodic boundary correctly in your trajectory before analyzing it. However, you should keep in mind that if you need to compute CVs on the fly to add a bias potential on those (as we will to in the next tutorials) you will have to learn how to reconstruct the molecules that are broken by the periodic boundary within PLUMED. If you know alternative tools that can reconstruct the molecules that are broken by the periodic boundary, it is a good idea to also use them and compare the result with PLUMED.

Exercise 2: Solving PBC issues and dump atomic coordinates

Analyze the provided trajectory traj-broken.xtc and use the DUMPATOMS action to produce new trajectories in gro format that contain:

  • The RNA duplex made whole (not broken by periodic boundary conditions). You should read carefully the documentation of WHOLEMOLECULES.
  • The whole RNA duplex aligned to a provided template (structure reference.pdb). See FIT_TO_TEMPLATE, using TYPE=OPTIMAL. Notice that you should provide to FIT_TO_TEMPLATE a pdb file with only the atoms that you wish to align. Use the ref.pdb file as a starting point and remove the lines non containing RNA atoms. More details on PDB files in PLUMED can be found here.
  • The whole RNA duplex and Mg ion, but only including the snapshots where Mg is at a distance equal to at most 4 A from phosphorous atom of residue 8. Search for the serial number of the proper phosphorous atom in the PDB file and use the UPDATE_IF action to select the frames.
  • The whole RNA duplex plus water molecules and mg ion wrapped around the center of the duplex. Compute first the center of the duplex with CENTER then wrap the molecules with WRAPAROUND. Make sure that individual water molecules are not broken after the move!

Here you can find a template input file to be completed by you.

Click on the labels of the actions for more information on what each action computes
tested on v2.9
# First load information about the molecule.
MOLINFO __FILL__
# Define here the groups that you need.
# Same as in the previous exercise.
rna: GROUP 
ATOMS
the numerical indexes for the set of atoms in the group. 
=__FILL__ 
mg: GROUP 
ATOMS
the numerical indexes for the set of atoms in the group. 
=__FILL__ 
wat: GROUP 
ATOMS
the numerical indexes for the set of atoms in the group. 
=__FILL__ 
# Make RNA duplex whole.
WHOLEMOLECULES __FILL__
# Dump first trajectory in gro format.
# Notice that PLUMED understands the format based on the file extension
DUMPATOMS 
ATOMS
the atom indices whose positions you would like to print out. 
=rna 
FILE
compulsory keyword 
file on which to output coordinates; extension is automatically detected 
=rna-whole.gro 
# Align RNA duplex to a reference structure
# This should not be the ref.pdb file but a new file with only RNA atoms.
FIT_TO_TEMPLATE 
REFERENCE
compulsory keyword 
a file in pdb format containing the reference structure and the atoms involved in
the CV. 
=__FILL__ 
TYPE
compulsory keyword ( default=SIMPLE )
the manner in which RMSD alignment is performed. 
=OPTIMAL 
# Notice that before using FIT_TO_TEMPLATE we used WHOLEMOLECULES to make RNA whole
# This is necessary otherwise you would align a broken molecule!
# Dump the aligned RNA on a separate file
DUMPATOMS 
ATOMS
the atom indices whose positions you would like to print out. 
=rna 
FILE
compulsory keyword 
file on which to output coordinates; extension is automatically detected 
=rna-aligned.gro 
# Compute the distance between the Mg and the Phosphorous from residue 8
d: DISTANCE 
ATOMS
the pair of atom that we are calculating the distance between. 
=mg,__FILL__ ## put the serial number of the correct phosphorous here 
# here we only dump frames conditioned to the value of d
UPDATE_IF 
ARG
the input for this action is the scalar output from one or more other actions. 
=d __FILL__
DUMPATOMS 
ATOMS
the atom indices whose positions you would like to print out. 
=rna,mg 
FILE
compulsory keyword 
file on which to output coordinates; extension is automatically detected 
=rna-select.gro 
UPDATE_IF 
ARG
the input for this action is the scalar output from one or more other actions. 
=d __FILL__# this command is required to close the UPDATE_IF above 
# compute the center of the RNA molecule
center: CENTER 
ATOMS
the list of atoms which are involved the virtual atom's definition. 
=rna 
# Wrap atoms correctly
WRAPAROUND 
ATOMS
wrapped atoms. 
=mg 
AROUND
reference atoms. 
=__FILL__ 
WRAPAROUND 
ATOMS
wrapped atoms. 
=wat 
AROUND
reference atoms. 
=center __FILL__# anything missing here? 
# Dump the last trajectory
DUMPATOMS 
ATOMS
the atom indices whose positions you would like to print out. 
=rna,wat,mg 
FILE
compulsory keyword 
file on which to output coordinates; extension is automatically detected 
=rna-wrap.gro

After you have prepared a proper plumed.dat file, you can use it with the following command 

> plumed driver --plumed plumed.dat --mf_xtc broken.xtc

Visualize the resulting trajectories using VMD. Since the gro files already contain atom names, you do not need to load the pdb file first. For instance, the first trajectory can be shown with

> vmd rna-whole.gro

TODO: I should perhaps add reference plots

To learn more: Mastering WHOLEMOLECULES

To learn more: Mastering FIT_TO_TEMPLATE

Conclusions

In summary, in this tutorial you should have learned how to use PLUMED to:

  • Manipulate atomic coordinates.
  • Compute collective variables.

All of this was done by just reading an already available trajectory. Notice that there are many alternative tools that could have been used to do the same exercise. Indeed, if you are familiar with other tools, it might be a good idea to also try them and compare the results. The special things of working with PLUMED are the following:

  • PLUMED implements a vast library of useful collective variables. Browse the manual and search for ideas that are suitable for your system.
  • PLUMED has a simple and intuitive syntax to combine collective variables ending up in descriptors capable to characterize complex conformational changes.
  • And finally, the most special thing: any collective variable that can be computed within PLUMED can also be biased while you are running your MD simulation! You will learn more later about this topic.

The last point is probably the main reason why PLUMED exists and what distinguishes it from other available software.