Manual

User Manual:

Open the PDF directly: View PDF .
Page Count: 128

Download
Open PDF In Browser	View PDF

PLUMED User’s Guide
A portable plugin for free-energy calculations
with molecular dynamics

Version 1.3.0 – Nov 2011

Contents
1 Introduction
1.1 What is PLUMED? . .
1.2 Supported codes . . .
1.3 Features . . . . . . .
1.4 New in version 1.3 .
1.5 Restrictions . . . . .
1.6 The PLUMED package
1.7 Online resources . . .
1.8 Credits . . . . . . . .
1.9 Citing PLUMED . . . .
1.10 License . . . . . . . .

.
.
.
.
.
.
.
.
.
.

2 Installation
2.1 Compiling PLUMED . . . . . . . . . . . . . .
2.1.1 Compiling the ACEMD plugin with
2.2 Including reconnaissance metadynamics . .
2.3 Testing the installation . . . . . . . . . . .
2.4 Back to the original code . . . . . . . . . .
2.5 The Python interface to PLUMED . . . . . .
3 Running free-energy simulations
3.1 How to activate PLUMED . . . . . . .
3.2 The input file . . . . . . . . . . . . .
3.3 A note on units . . . . . . . . . . . .
3.4 Metadynamics . . . . . . . . . . . . .
3.4.1 Typical output . . . . . . . .
3.4.2 Bias potential . . . . . . . . .
3.4.3 Well-tempered metadynamics
1

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

. . . . .
PLUMED
. . . . .
. . . . .
. . . . .
. . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

5
5
6
7
8
9
9
10
11
11
11

.
.
.
.
.
.

13
13
17
18
19
21
21

.
.
.
.
.
.
.

23
23
25
27
27
27
28
29

3.5
3.6

3.7
3.8
3.9
3.10

3.11
3.12
3.13

3.14

3.4.4 Restarting a metadynamics run . . . . . . . . . .
3.4.5 Using GRID . . . . . . . . . . . . . . . . . . . . .
3.4.6 Multiple walkers . . . . . . . . . . . . . . . . . .
3.4.7 Monitoring a collective variable without biasing it
3.4.8 Defining an interval . . . . . . . . . . . . . . . . .
3.4.9 Inversion condition . . . . . . . . . . . . . . . . .
Running in parallel . . . . . . . . . . . . . . . . . . . . .
Replica exchange methods . . . . . . . . . . . . . . . . .
3.6.1 Parallel tempering metadynamics . . . . . . . . .
3.6.2 Bias exchange simulations . . . . . . . . . . . . .
Umbrella sampling . . . . . . . . . . . . . . . . . . . . .
Steered MD . . . . . . . . . . . . . . . . . . . . . . . . .
3.8.1 Steerplan . . . . . . . . . . . . . . . . . . . . . .
Adiabatic Bias MD . . . . . . . . . . . . . . . . . . . . .
External potentials . . . . . . . . . . . . . . . . . . . . .
3.10.1 Walls . . . . . . . . . . . . . . . . . . . . . . . . .
3.10.2 Tabulated potentials . . . . . . . . . . . . . . . .
Commitment analysis . . . . . . . . . . . . . . . . . . . .
Projection of gradients . . . . . . . . . . . . . . . . . . .
Reconnaissance metadynamics . . . . . . . . . . . . . . .
3.13.1 Typical output . . . . . . . . . . . . . . . . . . .
3.13.2 Controlling the clustering . . . . . . . . . . . . .
3.13.3 Controlling the bias . . . . . . . . . . . . . . . . .
3.13.4 Restarting a simulation . . . . . . . . . . . . . . .
3.13.5 Using a subset of the defined cvs . . . . . . . . .
Driven Adiabatic Free Energy Dynamics (d-AFED) . . .
3.14.1 Input for d-AFED . . . . . . . . . . . . . . . . .
3.14.2 Typical output for d-AFED . . . . . . . . . . . .

4 Collective variables
4.1 Absolute position . .
4.2 Distance . . . . . . .
4.3 Minimum distance .
4.4 Angles . . . . . . . .
4.5 Torsion . . . . . . . .
4.6 Coordination number
4.7 Hydrogen bonds . . .
4.8 Interfacial water . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
2

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

30
30
35
35
36
38
40
40
41
42
45
46
46
48
50
50
50
53
53
54
54
56
57
58
59
60
61
62

.
.
.
.
.
.
.
.

64
66
67
68
69
70
70
71
73

4.9
4.10
4.11
4.12
4.13
4.14
4.15
4.16
4.17
4.18

4.19
4.20
4.21
4.22
4.23
4.24
4.25
4.26
4.27
4.28

Radius of gyration . . . . . . . . . . . . . . . . . . . . . . .
4.9.1 Gyration tensor based CVs . . . . . . . . . . . . . . .
Dipole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Dihedral correlation . . . . . . . . . . . . . . . . . . . . . . .
Alpha-beta similarity . . . . . . . . . . . . . . . . . . . . . .
Alpharmsd . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Antibetarmsd . . . . . . . . . . . . . . . . . . . . . . . . . .
Parabetarmsd . . . . . . . . . . . . . . . . . . . . . . . . . .
Electrostatic potential . . . . . . . . . . . . . . . . . . . . .
Puckering coordinates . . . . . . . . . . . . . . . . . . . . .
Path collective variables . . . . . . . . . . . . . . . . . . . .
4.18.1 Mean square deviation . . . . . . . . . . . . . . . . .
4.18.2 Distance mean square deviation . . . . . . . . . . . .
4.18.3 Contact map distances . . . . . . . . . . . . . . . . .
4.18.4 Using path variables as MSD, DMSD and CMAP and
the TARGETED statement . . . . . . . . . . . . . . . .
Contact Map . . . . . . . . . . . . . . . . . . . . . . . . . .
Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Helix loops . . . . . . . . . . . . . . . . . . . . . . . . . . . .
PCA projection . . . . . . . . . . . . . . . . . . . . . . . . .
SPRINT topological variables . . . . . . . . . . . . . . . . .
Radial distribution function . . . . . . . . . . . . . . . . . .
Angular distribution function . . . . . . . . . . . . . . . . .
Polynomial combination of CVs . . . . . . . . . . . . . . . .
Function of CVs . . . . . . . . . . . . . . . . . . . . . . . . .
A note on periodic boundary conditions . . . . . . . . . . . .

5 Postprocessing
5.1 Estimating the free energy after a metadynamics run
5.1.1 Installation instructions . . . . . . . . . . . .
5.1.2 Usage . . . . . . . . . . . . . . . . . . . . . .
5.2 Evaluating collective variables on MD trajectories . .
5.2.1 Installation instructions . . . . . . . . . . . .
5.2.2 Usage . . . . . . . . . . . . . . . . . . . . . .
5.3 Processing COLVAR files . . . . . . . . . . . . . . . . .
5.4 PLUMED as a standalone program . . . . . . . . . . . .
5.4.1 Installation instructions . . . . . . . . . . . .
5.4.2 Usage . . . . . . . . . . . . . . . . . . . . . .
3

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

74
75
76
77
77
78
79
80
81
82
83
85
87
87

.
.
.
.
.
.
.
.
.
.
.

91
92
93
93
94
96
97
98
99
100
101

.
.
.
.
.
.
.
.
.
.

103
103
103
103
106
106
106
108
109
109
110

5.5

5.6

Reweighting well-tempered metadynamics calculations
5.5.1 Installation instructions . . . . . . . . . . . . .
5.5.2 Usage . . . . . . . . . . . . . . . . . . . . . . .
Bias-exchange simulations via the linux shell . . . . . .

.
.
.
.

111
112
113
114

Chapter 1
Introduction
1.1

What is PLUMED?

PLUMED[1] is a plugin for free-energy calculations in molecular systems. It
works with some of the most popular classical molecular dynamics (MD)
codes, such as GROMACS [2], NAMD [3], DL POLY [4] , LAMMPS [5]
and the SANDER module in AMBER [6]. It also works with the very fast
CUDA/GPU MD code ACEMD [7] and, more recently, it has been extended
to work with ab initio MD codes, such as Quantum-ESPRESSO [8] and
CPMD.
Free-energy calculations can be performed as a function of many order
parameters, with a particular focus on biological problems, and using stateof-the-art methods such as metadynamics [9], umbrella sampling [10, 11, 12]
and Jarzynski-equation based steered MD [13, 14].
Here is a brief outline of this guide:
• In this chapter we give an overview of the features and restrictions of
the current release of PLUMED.
• In the second chapter we describe the procedure for installing the plugin
and testing the correct installation.
• The third chapter explains how PLUMED can be used to perform freeenergy calculations, setting up the input file and analyzing the output.
• The fourth chapter contains a list of collective variables (CVs) which
are implemented and allow a wide variety of physical and chemical
problems to be addressed.
5

• The fifth chapter is dedicated to postprocessing. It explains how to
reconstruct the free-energy profile from the output of a metadynamics
run and how to extract the CV values from MD trajectories.

1.2

Supported codes

PLUMED works as an add-on to some of the most popular MD codes: NAMD,
GROMACS, SANDER, DL POLY, Quantum-ESPRESSO, ACEMD, LAMMPS
and CPMD. These codes are not distributed with the PLUMED package, but
they must be obtained separately.
NAMD 2.6/2.7/2.8
http://www.ks.uiuc.edu/Research/namd/
GROMACS 4.0/4.5.x
http://www.gromacs.org/
AMBER 10/11
http://ambermd.org/
DL POLY 2.20
http://www.cse.scitech.ac.uk/ccg/software/DL POLY/
Quantum-ESPRESSO 4.3.2
http://www.quantum-espresso.org/
ACEMD 1.2
http://multiscalelab.org/acemd
LAMMPS 27-Oct-2011
http://lammps.sandia.gov
CPMD 3.15.1
http://www.cpmd.org
Note that at the moment of this release of PLUMED, only those specific
versions of the listed codes have been tested. Moreover, some other codes
are available through the Python interface provided by Rosa Bulo.
Amsterdam Density Functional (ADF)
http://www.scm.com/
6

and all the codes supported by
Atomic Simulation Environment (ASE)
https://wiki.fysik.dtu.dk/ase/
Additionally there exist some porting to
FHI-AIMS
https://aimsclub.fhi-berlin.mpg.de/
which is not directly mantained by the PLUMED developer team. Additional information have to be retrieved by the developers of AIMS themselves.

1.3

Features

PLUMED can perform several different types of calculation:
• Metadynamics with a large variety of CVs [9];
• Well-tempered metadynamics [15];
• Multiple walkers metadynamics [16];
• Combined parallel tempering and metadynamics [17];
• Bias-exchange metadynamics [18];
• Reconaissance metadynamics [19];
• Steered MD;
• Umbrella sampling;
• Adiabatic biased molecular dynamics [20];
• Commitment analysis.

1.4

New in version 1.3

PLUMED version 1.3 presents several new features, including new collective
variables and support to new codes. Among these:
• Reconaissance Metadynamics [19];
• Driven adiabatic free energy dynamics (d-AFED, contributed by Michel
Cuendet)
• Definition of CVs as Polynomial combination of CVs and Function on
CVs;
• New INTERVAL keyword to limit a region to be sampled using Metadynamics (contributed by Alessandro Laio)
• Better scalability with Gromacs;
• New collective variables: Projection on PCA eigenvectors (contributed
by Ludovico Sutto), SPRINT topological variable and gyration tensor
based CVs;
• New tool: reweight well-tempered metadynamics simulations [21];
• Added support to ACEMD 1.2, Lammps (27-oct-2011), Quantum Espresso
4.3.2, CPMD 3.15.1, NAMD 2.8, Gromacs 4.5.x and Amber11;
• Python interface (contributed by Rosa Bulo);
• New tool to perform bias-exchange simulations via the linux shell with
any MD code;
• Bugs have been fixed and the manual has been improved thanks to the
feedbacks of many users.
A full list including also the bug fixed in the current release can be found
in the CHANGES file distributed with the package.

1.5

Restrictions

The current release of PLUMED has a few restrictions:
• Parallel tempering plus metadynamics and bias-exchange metadynamics are available only in the GROMACS version (however a tool allows
bias-exchange simulations via linux shell with any MD engine);
• The patched version of GROMACS cannot perform hamiltonian (lambda)
replica exchange;
• The ENERGY collective variable is available only for GROMACS, AMBER and DL POLY and cannot be used with multiple time step algorithm;
• Only orthorhombic cells are supported in ACEMD and DRIVER;
• Only orthorhombic and truncated octahedron cells are supported in
AMBER;
• Only Car-Parrinello and Born-Oppenheimer simulations are supported
in CPMD;
• Only NVT simulations are supported because the contribution to the
virial is not calculated.

1.6

The PLUMED package

The plugin package has the following directory structure:
• common_files. The directory containing all the basic routines.
• tests. A variety of examples of different CVs and free-energy methods
provided with topology and input files for the supported codes. These
examples, combined with a script adapted from CP2K [22], work also
as a regtest for the plugin.
• patches. A collection of patches to interface PLUMED with different
codes.
• utilities. A few small utilities:
9

– sum_hills is a post-processing program which reads the HILLS
file produced by the plugin in a metadynamics simulation and
returns the free energy by summing the Gaussians that have been
deposited.
– driver is a tool that calculates the value of selected CVs along a
MD trajectory. It requires a PDB file, a trajectory in DCD format
and a file with the same syntax of the PLUMED input file.
– standalone is a program to run PLUMED as a standalone code.
– plumedat.sh is a shell script to extract information from PLUMED
output files.
– reweight is a tool to calculate the unbiased probability distribution of variables other than the CVs in well-tempered metadynamics simulations
– python_interface a set of tools that enables to use PLUMED with
Python (hey, a Python which is plumed with feathers is a really
weird animal!).
– bias-exchange a set of tools that allow to use PLUMED to perform
bias-exchange simulations via the linux shell using any MD engine.
• manual. This manual.
• ACEMD. It contains the files needed to compile the plugin for ACEMD.
• tutorial. It contains the resources of the PLUMED tutorial 2010 (http:
//sites.google.com/site/plumedtutorial2010/)
• recon_src. It contains the source files specific for the Reconaissance
Metadynamics.

1.7

Online resources

You can find more information on the web:
http://www.plumed-code.org
For any questions, please subscribe to our mailing list:
plumed-users@googlegroups.com
10

1.8

Credits

PLUMED has been developed by Massimiliano Bonomi, Davide Branduardi,
Giovanni Bussi, Carlo Camilloni, Davide Provasi, Paolo Raiteri, Davide
Donadio, Fabrizio Marinelli, Fabio Pietrucci, Francesco Luigi Gervasio and
others. However, this work would not have been possible without the joint
effort of many people. Among these, we would like to thank (in alphabetical
order): Alessandro Barducci, Anna Berteotti, Rosa Bulo, Matteo Ceccarelli,
Michele Ceriotti, Paolo Elvati, Antonio Fortea-Rodriguez, Alessandro Laio,
Matteo Masetti, Fawzi Mohamed, Ferenc Molnar, Gabriele Petraglio, Jim
Pfaendtner and Federica Trudu. Francesco Marini is kindly acknowledged
for his technical support and Joost VandeVondele for permission to use his
regtest script.
Some PLUMED users have also contributed to the implementation of new
features and the debugging of old bugs. Among these, we would like to thank:
Toni Giorgino, Marcello Sega, Emmanuel Autieri, Gareth Tribello, Andrea
Coletta, Ludovico Sutto, Layla Martin-Samos, Walter Rocchia, Alessio Lodola,
Luigi Capoferri, Michel Cuendet, Jiri Vymetal, Fahimeh Baftizadeh, and
Katsumasa Kamiya.

1.9

Citing PLUMED

You may wish to cite the following reference if you have utilized PLUMED in
your work:
M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri,
D. Donadio, F. Marinelli, F. Pietrucci, R.A. Broglia and M. Parrinello.
PLUMED: a portable plugin for free-energy calculations with molecular
dynamics, Comp. Phys. Comm. 2009 vol. 180 (10) pp. 1961-1972.

1.10

License

PLUMED is free software: you can redistribute it and/or modify it under the
terms of the GNU Lesser General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option)
any later version. PLUMED is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See
11

the GNU Lesser General Public License for more detail. You should have
received a copy of the GNU Lesser General Public License along with PLUMED.
If not, see http://www.gnu.org/licenses/.

Chapter 2
Installation
The plugin installation requires the molecular dynamics code to be recompiled after the calls to the plugin routines have been inserted at the appropriate places in the original program. For a list of the supported molecular
dynamics codes see Sec. 1.2.
All the basic plugin routines are contained in the common files folder of
the PLUMED distribution package. The insertions are automatically performed
by a series of scripts provided in the patches directory.
In the following we will briefly describe the procedure for applying these
patches to the different MD codes supported by PLUMED. We will refer to the
root directory of the distribution package as PLUMED root.

2.1

Compiling PLUMED

A similar procedure can be followed for all the supported codes except
ACEMD. When code-specific procedures are needed, we state it explicitely.
A different proceedure is required for ACEMD (see below).
Different patches are available for different code versions. The name of
the suitable patch is plumedpatch CODENAME CODEVERSION.sh. If the patch
corresponding to the exact CODEVERSION that you are using is not available, choose the closest match. In the following, we shall indicate with
plumedpatch.sh the proper patch script.
• Extract the source code from its archive and then move into its root
directory.

• Configure the code as usual (not necessary for DL POLY)
• NAMD and LAMMPS only: In the plumedpatch script, modify the
myarch variable to match your architecture.
• Set the environmental variable plumedir to point to PLUMED root.
• Copy or link the plumedpatch.sh file from the patches folder to the
current directory.
• Execute the script: ”./plumedpatch.sh -patch ”.
• DL POLY-only: copy the proper Makefile from build/ to srcmod/ (or
src/ in older versions)
• CPMD only: you may need to change variables cxx, cxxflag, extraflag
in the plumedpatch script depending on your computer architecture:
see the examples below.
• Compile the code as usual.
Further details can be found in the patches/README file.
Example.
This is the procedure for compiling the serial version of AMBER 10 with PLUMED using g95 in the Bourne
shell.
tar zxf AMBER10.tgz
cd amber10/
export plumedir="PLUMED root"
cp $plumedir/patches/plumedpatch sander 10.sh .
cd src/
./configure g95
cd ..
./plumedpatch sander 10.sh -patch
make

Example.
This is the procedure for compiling the serial version of GROMACS, using the GNU compilers.
tar zxf gromacs-4.0.5.tar.gz
cd gromacs-4.0.5
export plumedir="PLUMED root"
cp $plumedir/patches/plumedpatch gromacs 4.0.4.sh .
CC=gcc CXX=g++ ./configure
./plumedpatch gromacs 4.0.4.sh -patch
make
make install

Example.
This is the procedure for compiling NAMD on an Intel Mac using the GNU g++ compiler and the FFTW.
tar zxf NAMD 2.6 Source.tar.gz
cd NAMD 2.6 Source
export plumedir="PLUMED root"
cp $plumedir/patches/plumedpatch namd 2.6.sh .
./config fftw MacOSX-i686-g++
Edit ./plumedpatch namd 2.6.sh by setting the myarch variable to MacOSX-i686-g++.
./plumedpatch namd 2.6.sh -patch
cd MacOSX-i686-g++
make

Example.
This is a sample procedure for compiling the scalar version of DL POLY 2.20 with the gfortran compiler
in the Bourne Shell environment.

tar zxf dl poly 2.20.tar.gz
cd dl poly 2.20
export plumedir=PLUMED root
cp $plumedir/patches/plumedpatch dlpoly 2.20.sh .
./plumedpatch dlpoly 2.20.sh -patch
cp build/MakeSEQ srcmod/Makefile
cd srcmod
make gfortran

Example.
This is a sample procedure for compiling the openmpi version of LAMMPS with the mpic++ compiler in
the Bourne Shell environment. Please, note that the LAMMPS tarball comes with the download date but
unless major changes are done in the host code, the patching procedure will stay unchanged and the
latest available patching file has to be used.
tar xvf lammps-15Jan10.tar
cd lammps-15Jan10/
Edit ./src/MAKE/Makefile.openmpi to suit your machine.
export plumedir="PLUMED root"
cp $plumedir/patches/plumedpatch lammps 15-01-2010.sh ./
Edit ./plumedpatch lammps 15-01-2010.sh to set myarch to openmpi.
./plumedpatch lammps 15-01-2010.sh -patch
cd src
make openmpi

Example.
This is a sample procedure for compiling CPMD-3.15.1 under different computer architectures.
tar zxvf cpmd-v3 15 1.tgz
cd CPMD
Generate a valid Makefile for your machine using the script ./mk config.sh, e.g.:
./mk config.sh PC-GFORTRAN > Makefile
export plumedir="PLUMED root"
cp $plumedir/patches/plumedpatch cpmd-3.15.1.sh .
If necessary, modify the variables cxx (a valid C++ compiler), cxxflag (the corresponding flag to
avoid exception handling), and extraflag (see below) in the script plumedpatch cpmd-3.15.1.sh, e.g.
for GNU compilers:
cxx="mpiCC"
cxxflag="-fno-exceptions"
for IBM machines like AIX-power7:
cxx="xlc r -c"
cxxflag="-qnoeh"
for IBM BlueGene:
cxx="bgxlc++ r -c -O"
cxxflag="-qnoeh"
for PGI compilers:
cxx="pgCC"
cxxflag="--no exceptions"
To successfully compile on some IBM architectures you may also need to set one or both of the
following extra flags:
extraflag="-DPLUMED CPMD NOUNDERSCORE -DPLUMED AIXXLF"
Finally, patch and compile:
./plumedpatch cpmd-3.15.1.sh -patch
make

2.1.1

Compiling the ACEMD plugin with PLUMED

PLUMED works with ACEMD 1.2. As the source code for ACEMD is not
available, PLUMED will work as a plug-in. For this reason the compilation
will be slightly different from that of the other codes. Moreover, given the
unavailability of source code, this version of PLUMED will not be supported
by the PLUMED development team.
• Extract the source code from its archive.
• Set the environmental variable plumedir to point to PLUMED root.
17

• Use make to compile the plugin inside the plumedir/ACEMD folder.

Example.
This is a sample procedure for compiling the ACEMD plug-in in the Bourne Shell environment.

export plumedir=PLUMED root
cd $plumedir ACEMD
make

Once the plumed.so is compiled, copy it where you will run ACEMD. Add
the following lines to the ACEMD input file:
pluginload testplug ./plumed.so
pluginarg testplug input META_INP
pluginarg testplug boxx xx
pluginarg testplug boxy yy
pluginarg testplug boxz zz
pluginfreq 1
where META INP is the PLUMED input file and xx,yy,zz are the box dimensions. At this time the support of the ACEMD plugin will be provided
by the ACEMD developers. Please note that only orthorhombic cells are
available in ACEMD 1.2.

2.2

Including reconnaissance metadynamics

The reconnaissance metadynamics routines (see section 3.13 included in PLUMED
make use of the C++ standard library and lapack. On some machines these
libraries are not installed. Hence, PLUMED has been written so that by default
no attempt to compile reconnaissance metadynamics is made. Consequentially, if you wish to perform reconnaissance metadynamics you must patch
PLUMED slightly differently.
To include reconnaissance metadynamics in your patched MD code is
simply a matter of, during the patching procedure, telling PLUMED the location of the C standard library, the location of lapack and the location of
18

a c++ compiler. This is done by making some small modifications to the
plumedpatch file. If this file is opened using a text editor you will notice at
least one of the following two lines:
RECON_LIBS=
RECON_CPP=
If at least one of these lines is not present then reconnaissance metadynamics is not implemented for the particular code you are patching. The
first of these two lines tells PLUMED the locations of the lapack and c standard libraries. If the space after the equal sign is left blank PLUMED patches
in its default fashion and no attempt is made to compile the reconnaissance
metadynamics routines. By contrast if there is any character after the equals
sign PLUMED will endeavor to compile the reconnaissance metadynamics routines. For a correct compilation the RECON_LIBS}variable should indicate the
locations of the lapack (i.e. -llapack) and c standard libraries (i.e. -lstdc++).
The RECON_CPP= variable should be set equal to a c++ compiler (e.g. g++).
Example.
To include the reconnaissance metadynamics routines in the PLUMED compilation the use must change
the values of RECON LIBS and RECON CPP in the plumedpatch file:
RECON LIBS=”-llapack -lstdc++”
RECON CPP=”g++”

There are some differences between the procedures for patching PLUMED+
reconnaissance metadynamics in the different MD codes between the different
codes. This is because some of the codes already have lapack included or are
compiled with a c++ compiler. These differences are described in table 2.1:

2.3

Testing the installation

Once the installation process has been completed successfully, the user is
encouraged to test the chosen MD package for any problems. The tests
directory contains regtest scripts for the different MD codes. These also
serve as a regularity test in case the user implements his own modifications
and as a basic illustration of the capabilities of the plugin. The user should
19

Code
NAMD
GROMACS
AMBER
DL POLY
Q-ESPRESSO
ACEMD
LAMMPS
CPMD

Reconnaissance
implemented
yes
yes
no
yes
yes
no
no
yes

C++ compiler Lapack
required
required
no
yes
yes
no
yes
yes
yes
no
yes
no

lstc++
required
yes
yes
yes
yes
yes

Table 2.1: Details on compiling PLUMED compilation with the various codes
edit the test script, setting up the path where the test suite is found and
giving the location of the executable.
For GROMACS users only. Please note that:
• The tests for GROMACS are designed for and should be executed with
the double-precision version of the code;
• Biasxmd, ptmetad, dd and pd are designed for the parallel version of
GROMACS. The user should specify in the test script the location of
the parallel executable and the version of GROMACS used. These tests
will fail if the parallel version of GROMACS has not been compiled.

Example.
In the case of NAMD, the regtest script is called do regtest namd.sh. Here, the user should modify the
dir base and the namd prefix:
dir base=/Programs/md meta/tests/namd
namd prefix=/chicco/bin/namd plugin/namd2.
Regtest scripts for the other programs require an identical procedure.

The script executes a list of tests and then compares the results with the
outcome of previously run simulations. Please note that when the script is
run for the first time it produces the reference. Finally, the script produces
a summary of the results. The test result can be one of the following:
20

• ’OK’ if the results match those of a previous run precisely. The execution time is also given.
• ’NEW’ if they have not been executed previously. The reference result
is generated automatically in this run. Tests can also be ’NEW’ if they
have been reset, i.e. been newly added to the TEST FILES RESET
files.
• ’RUNTIME FAILURE’ if they stopped unexpectedly (e.g. core dump,
or stop).
• ’WRONG RESULT’ if they produce a result that deviates from an old
reference.
The last two options generally mean that a bug has been introduced,
which requires investigation. Since regtesting only yields information relative
to a previously known result, it is most useful to do a regtest before and after
you make changes.

2.4

Back to the original code

At any time the user may want to ”unpatch” the MD code and revert back
to the original distribution. To do so, the user should go to the directory
where the PLUMED patch has been copied and type ./plumedpatch -revert.

2.5

The Python interface to PLUMED

Starting from version 1.3 PLUMED has also a Python interface contributed by
Rosa Bulo. It is possible to find it in the utilities/python_interface
directory. Inside it you might find a README file and a Makefile . After
copying (or linking) all the *.c and *.h files from the common_files directory you can immediately do:
make gnu
if you want to compile with gnu compilers. Intel compilers are also available. This creates a file that is libplumed.so that will be loaded runtime
by Python. You need to add the directory where this file resides into the
LD LIBRARY PATH so that Python will be able to find it and load it runtime.
The extension amber is only reflecting the fact that the interface is built by
21

exploiting the simple AMBER interface. Therefore the limitiations regarding
AMBER in the functionalities apply also here.
For MacOSX users the compiler flags should be changed accordingly. For
gnu compilers it requires to change the linking flags from
-shared -Wl,-soname,libplumed.so
to
-bundle -flat namespace -undefined suppress.
Next, in the directory
utilities/python_interface/pylib
you can find the file plumedamberlib.py which is the actual Python module that will deal of loading libplumed.so and create the suitable Python
interface. This provides the init metadyn and cv calculation methods
that one can practically use through Python.
A practical use of this shown by the Python script
/utilities/python interface/test/run water amber.py
where a water molecule is simply displaced and the collective coordinates are
updated at each step. Notably, this interface is embodied in the Python tools
provided with the Amsterdam Density Functional code (also called ADF, see
http://www.scm.com/).
Additionally Rosa provided an interface with the powerful Atomic Simulation
Environment ( see https://wiki.fysik.dtu.dk/ase/ ). This enable PLUMED
to be interfaced with plenty of DFT and ab-initio codes. Some examples
are present in the directory utilities/python interface/test and are
ASEmd plumed emt water.py and ASEnvt plumed emt water.py .
Many thanks to Rosa for the contribution!

Chapter 3
Running free-energy
simulations
In this chapter we describe how to activate PLUMED and how to create the
correct PLUMED input file for a specific type of free energy calculation. The
typical output of these calculations is also explained in detail.

3.1

How to activate PLUMED

PLUMED input is contained in one single file, named plumed.dat by default,
which defines the CVs, the type of run to be performed and the parameters
for the bias potential generation.
• The users of NAMD, SANDER and DL POLY can instruct the code
to parse the PLUMED input file by setting the PLUMED variable to on (or
1 for SANDER) in the MD input file. It is also possible to change
the default name for the PLUMED input file by setting the plumedfile
variable in the MD input;
• GROMACS users should activate PLUMED on the command line by specifying the flag -plumed followed by the input file name. The extension
of such a file must be .dat;
• Quantum-ESPRESSO users should activate PLUMED on the command
line by specifying the flag -plumed. The name of the PLUMED input file
is hardcoded as plumed.dat;
23

• In LAMMPS, PLUMED is activated using the fix plumed. The input file
is specified by plumedfile, the output file by outfile;
• In CPMD, PLUMED is activated on the command line by specifying the
flag -plumed after the CPMD input file and the pseudopotential directory. The PLUMED input file name is plumed.dat.
Example.
A typical SANDER input file for a metadynamics calculation:
METADYNAMICS TEST
&cntrl
imin=0, irest=0, ntx=1, ig=71278,
nstlim=1001, dt=0.0002,
ntc=1, ntf=1,
ntt=3, gamma ln=5,
tempi=300.0, temp0=300.0,
ntpr=200, ntwx=0,
ntb=0, igb=0,
cut=999.,
plumed=1 , plumedfile=’plumed.dat’
/

Example.
The NAMD and DL POLY input file should contain, in addition to the usual keywords defining the run, the
following lines;
plumed on
plumedfile plumed.dat

Example.
To perform a free-energy calculation with GROMACS, PLUMED must be activated on the command line:
mdrun -plumed plumed.dat ...

Example.
To perform a free-energy calculation with Quantum-ESPRESSO, PLUMED must be activated on the
command line:
pw.x -plumed ...

Example.
The LAMMPS input file should contain, in addition to the usual keywords defining the run, the following
line;
fix 3 all plumed plumedfile plumed.dat outfile metaout.dat

Example.
To perform a free-energy calculation with LAMMPS the file should contain, in addition to the usual
keywords defining the run, the following command:
fix ID all plumed plumedfile plumed.inp outfile plumed.out
where ID is the user assigned fix name.

Example.
To perform a free-energy calculation with CPMD, PLUMED must be activated on the command line:
cpmd.x input .

-plumed

where input is the CPMD input file, ”.” is the current directory (if pseudopotentials are located
elsewhere, substitute with the right location), and the PLUMED input file is called plumed.dat.

3.2

The input file

In the following sections we describe the syntax used in the PLUMED input file.
The commands contained in this file can be divided in two groups according
to the functionality required. A first group of commands defines the type of
simulation (metadynamics, steering and umbrella sampling, replica exchange
methods) and the parameters needed for the chosen algorithm. These are
described in this chapter. A second part defines the degrees of freedom on
which the algorithms operate, the so-called collective variables or CVs. The
details of such commands are described in the next chapter.
As a general rule, each setting is defined by a principal keyword that
must be placed at the beginning of the line, in upper case. Additional input
pertaining to the setting can be specified on the same line, using additional
25

keywords that can be added in any order. A line can be continued to the
next one by adding a backslash or an ampersand as a last character in the
line. Three kinds of keyword may exist: the directives, which must be placed
at the beginning of a line and define the argument of the line, the keywords,
which specify the attributes of the different fields in the line and flags, which
simply turn a given option on or off.
Example.
The following is an example of a complete PLUMED input file. In this case the input defines a well-tempered
metadynamics run with two CVs; the first is the distance between two atoms, and the second a dihedral
angle.
# general options
HILLS HEIGHT 0.1 W STRIDE 100
WELLTEMPERED SIMTEMP 300 BIASFACTOR 10
# print each 50 time units and add a time offset of 20.0 to COLVAR
PRINT W STRIDE 50 T OFFSET 20.0
# definition of CVs
NOTE distance between hydrogens
DISTANCE LIST 13 46 SIGMA 0.35
NOTE torsional angle
TORSION LIST 1 4 65 344 SIGMA 0.1
# wall on the CV
UWALL CV 1 LIMIT 15.0 KAPPA 100.0 EXP 4.0 EPS 1.0 OFF 0.0
ENDMETA

The ENDMETA directive is the last line read from the PLUMED input file and
any line added after this keyword will be ignored.
The symbol # allows the user to comment any line in the input file. The
directive NOTE allows the user to place comments which are copied to the
PLUMED log file.
The PRINT allows to dump a file called COLVAR (see more in 3.4) which is
written each W STRIDE timesteps and (optional) a time offset can be specified
through the T OFFSET keyword. To append the colvar values to an existing
COLVAR file, the flag APPEND can be used.

3.3

A note on units

The values in the PLUMED input file are read in the internal units for the MD
engine. For DL POLY the energy in input can be in different user-specified
units; to be consistent, the values in the PLUMED input file must be in the
same units specified in the FIELD file.

3.4

Metadynamics

The metadynamics algorithm applies additional forces to a standard molecular dynamics simulation [9, 23, 24]. In this case, the PLUMED input file must
contain the definition of at least one CV (see chapter 4 for the required
syntax), and the HILLS keyword which defines the details of the bias potential (see section 3.4.2). In this case the biasing potential will be calculated
and applied to the microscopic degrees of freedom during the run. Before
discussing the additional optional commands (section 3.4.2), we give a brief
overview of the file produced in a typical metadynamics run.

3.4.1

Typical output

Standard metadynamics will produce, in addition to the usual output files
generated by the MD engine, a file called COLVAR that contains the following
data:
• The first column contains the time step;
• The following d columns contain the values of the CV(s);
• Two additional columns contain the value V (s, t) of the bias potential
at the given point and time, and the potential due to the walls (if
defined).
Since PLUMED 1.2, a more flexible format for COLVAR has been introduced.
Here, the first line of the COLVAR file is a header. The line begins with #,
so as to be ignored by plotting programs such as gnuplot and xmgrace. A
small script plumedat.sh to interpret this file is provided in the utilities.

Another file, called HILLS, contains the information of the biasing potential needed to estimate the free-energy, and to restart metadynamics runs.
The bias potential is given by:


d
X

(0)



(si − si (kτ ))2 
V (~s, t) =
W (kτ ) exp −
.
2σi2
i=1
kτ sw , the history dependent potential is still updated according

to Eq. 3.1, but the metadynamics force is set to zero for s < sw :
∂VG (s, t)
= 0 ∀s < sw
∂s
Notice that Gaussians are added also if s < sw , as the tails of these Gaussians
influence VG in the relevant region s > sw . In this way, the force on the system
in the region s > sw comes from both metadynamics and the force field, in the
region s < sw only from the latter. This approach allows obtaining a historydependent bias potential VG that fluctuates around a stable estimator, equal
to the negative of the free energy far enough from the boundaries [26].
• In order to obtain a parallel growing VG , one should place the interval
limit sw in a region where the free energy derivative is not large. For
example, if one uses the keyword with a coordination number CV (see
Section 4.6), that is intrinsically limited to s > 0, the correct choice for
sw is not 0, but a number of the order of the width of the Gaussians
(σ in Eq. 3.1).
• This remedy has the advantage of being robust and parameter-free, but
works only for one-dimensional biases.
• If in the region s < sw the system has a free energy minimum, the
INTERVAL keyword should be used together with a soft wall at sw
(keyword LWALL, see section 3.10), namely a time-indipendent external potential that prevents the system from going there and remaining
trapped in the minimum. An example is provided below.
Example.
The following input file contains the interval condition in presence of soft walls:
HILLS HEIGHT 0.1 W STRIDE 500
COORD LIST NN 8 MM 12 D 0 0 R 0 0.25 SIGMA 0.15
g1->
26 27
g1
29 30
g2 χ1 another Gaussian centered in −s and with the same width is
added. In this case, the height of the extra Gaussians depends on V
and is given by:
w = [2V (0, t) − V (−s, t) − V (s, t)] y(s),
h

(3.4)

where y(s) = 1/ 1 + (s/(χ2 χ1 ))10 and χ2 > χ1
The second factor in Eq. (3.4) is approximately one for | s |< χ2 χ1 and goes
to zero for | s |> χ2 χ1 . This ensures that V goes smoothly to zero out of the
boundaries. χ2 χ1 thus defines the width of the inversion interval and χ2 is
a scaling factor regulated by the keyword INVERSION. χ1 is regulated by the
keyword REFLECTION (χ1 is in σ units). The keyword MAXHEIGHT defines the
largest | w | in gaussian height units and it regulates the speed of variation
of V out of the boundaries. The INVERT keyword can be used in presence of
an external potential that limits the CVs space exploration (see section 3.10)
like in the example below;

Example.
The following input file contains the inversion condition in presence of soft walls:
HILLS HEIGHT 0.1 W STRIDE 500
COORD LIST NN 8 MM 12 D 0 0 R 0 0.25 SIGMA 0.15
g1->
26 27
g1
29 30
g2 SIGMA 0.1
CA->
20 22 26 30 32
CA<# end of the input
ENDMETA

The PTMETAD directive switches on parallel tempering. All replicas have
the same CVs, in this particular case the radius of gyration defined by the
group of atoms .
The Gaussian height set by the keyword HEIGHT is automatically rescaled
with temperature, following Wi = W0 TT0i , where i is the index of a replica and
Ti its temperature.
Similarly, the simulation temperature needed to use the well-tempered
algorithm (which, for non-parallel simulations is set with the SIMTEMP keyword) is here taken directly from the GROMACS input at each replica. As
a result, the value of ∆T for the well-tempered algorithm is rescaled across
the replicas.
The plugin will produce one COLVAR file and one HILLS file for each
replica.
41

3.6.2

Bias exchange simulations

Bias exchange simulations must be run using the BIASXMD directive. Only
the GROMACS engine currently has this feature implemented within a single
parallel run, whereas the bias-exchange tool in the directory utilities allows to perform bias-exchange simulations with any MD engine via the linux
shell (see Section 5.6).
The procedure for running bias-exchange simulations is similar to the
one described earlier for parallel tempering (i.e., add -multi nrep -replex
nexch to the mdrun command line, where nrep is the number of replicas and
nexch is the number of steps between attempting exchanges). One plugin
input and one binary topology file must be provided for each replica. These
files must be named with the replica index (e.g., plumed0.dat, plumed1.dat,
..., md0.tpr, md1.tpr, ...). Each plugin input file must contain the directive
BIASXMD and the same list of CVs, in the same order. Each replica usually
has only one (or few) of the CVs active: use a NOHILLS directive for each
of the inactive CVs (or do not specify SIGMA for them). It is also possible
to have replicas with all CVs inactive: such unbiased replicas can exchange
with biased ones, jumping in this way beyond free-energy barriers, but since
they are not subject to the hills they tend to populate the free-energy basins
according to equilibrium statistics.

Example.
The following two input files define a bias-exchange metadynamics run:
# --- input file plumed0.dat --# switching on metadynamics and Gaussian parameters
HILLS HEIGHT 0.1 W STRIDE 500
# switching on bias-exchange
BIASXMD
# instruction for CVs printout
PRINT W STRIDE 50
# instruction for labeling of COLVAR and HILLS files
HILLS LABEL A
# the CVs: dihedral angles
TORSION LIST 1 5 11 13 SIGMA 0.314
TORSION LIST 11 13 15 21 SIGMA 0.314
# in this replica we bias only CV 1:
NOHILLS CV 2
# end of the input
ENDMETA

# --- input file plumed1.dat --# switching on metadynamics and Gaussian parameters
HILLS HEIGHT 0.1 W STRIDE 500
# switching on bias-exchange
BIASXMD
# instruction for CVs printout
PRINT W STRIDE 50
# instruction for labeling of COLVAR and HILLS files
HILLS LABEL B
# the CVs: dihedral angles
TORSION LIST 1 5 11 13 SIGMA 0.314
TORSION LIST 11 13 15 21 SIGMA 0.314
# in this replica we bias only CV 2:
NOHILLS CV 1
# end of the input
ENDMETA

In bias-exchange runs, the plugin will produce one COLVAR file and one
HILLS file for each replica. What is actually exchanged between replicas are
the atomic coordinates, not the bias, so that replica 0 (printing out COLVAR0)
will always be biased by HILLS0, replica 1 (printing out COLVAR1) will always
be biased by HILLS1,) and likewise for the other replicas.
In the example above, note the use of the keyword HILLS LABEL: this will
print a comment line ”#! ACTIVE 1 1 A” in files COLVAR0 and HILLS0 and
”#! ACTIVE 1 2 B” in files COLVAR1 and HILLS1, indicating the number of
active CVs, their index, and the chosen label. This facilitates the reconstruction of the multidimensional free-energy landscape from the bias-exchange
simulation according to the weighted-histogram algorithm in Ref. [29], e.g.
43

employing the METAGUI program [30] downloadable from
http://www.plumed-code.org/contributions
See also Section 5.6 for informations on bias-exchange simulations with
MD codes different from GROMACS.

3.7

Umbrella sampling

The directive UMBRELLA allows umbrella sampling calculations to be performed on the CV specified by the parameter keyword CV. The position s0
of the umbrella restraint is determined by the keyword AT, and the spring
constant k - in internal unit of the main code - by the keyword KAPPA. Optionally, also a constant force m can be specified by the keyword SLOPE. This
turns on a potential of the following functional form:
1
Vumb (s) = k(s − s0 )2 + m(s − s0 ).
2

(3.5)

Example.
The following input file defines an umbrella potential acting on the first CV and centered in s0 = −1.0.
TORSION LIST 13 15 17 1
UMBRELLA CV 1 KAPPA 200 AT -1.0
PRINT W STRIDE 100
ENDMETA

In the case of umbrella sampling runs, the value of the CVs is printed on
the COLVAR file with a stride fixed by the directive PRINT and the keyword
W STRIDE. The file contains the time, the CV values, the metadynamics potential, the harmonic potential of umbrella sampling, the CV on which the
restraint acts and the position of the restraint. The final calculation of the
free energy as a function of this CV can be done using the weighted histogram
analysis method, choosing one of the many possible implementations, for instance the wham code by Alan Grossfield [31]. The directive UMBRELLA can
be used multiple times in the case of multidimensional umbrella sampling
calculations (one directive for each CV).
The keyword RESTART can be used when restarting an umbrella sampling
calculation to append the value of the CVs on the COLVAR file.

3.8

Steered MD

PLUMED can be used to drag a system to a target value in CV space using
an harmonic potential moving at constant speed. If the process is reversible,
i.e. for velocities tending to zero, the work done in the dragging corresponds
to the free energy. In the case of finite velocity it is possible to obtain an
estimate of the free-energy from the work distribution using Jarzynski [13]
or Crooks [14] relations.
The directive STEER activates the steering on the collective variable specified by the keyword CV. The target value is determined by the keyword TO,
the velocity, in the same unit as the corresponding CV every 1000 steps, by
VEL and the spring constant by KAPPA. An additional keyword FROM can be
used to specify the starting point of the dragging, otherwise the CVs values
at the first step are taken as the starting point. The functional form of the
dragging potential is the same as the one of formula 3.5, with the reference
position s0 moving at the specified velocity.
The keyword RESTART can be used when restarting a steered MD calculation to append the value of the CVs on the COLVAR file.
Example.
The following input file defines a steered MD on the angle CV to a target value of 3.0 rad.
ANGLE LIST 13 15 17
STEER CV 1 TO 3.0 VEL 0.5 KAPPA 500.0
PRINT W STRIDE 100
ENDMETA

3.8.1

Steerplan

PLUMED can be used to drag a system on a pathway that is composed of
successive steering runs on chosen degrees of freedom in a planned fashion.
While the directive STEER is rather easy and intuitive, STEERPLAN is more
flexible and it allows to avoid a lot of scripting whenever you plan to simulate
complex transitions by means of out-of-equilibrium runs.
The directive STEERPLAN activates this option and reads the name of a
file that contains the plan as an input.

Example.
The following input file contains many collective variables and a steerplan directive.

PRINT W STRIDE 5
#
# these CVs are put here just to show that you may use more variable than the
# ones defined in the steerplan
#
TARGETED TYPE MSD FRAMESET restrained.pdb
TARGETED TYPE MSD FRAMESET waters.pdb
TARGETED TYPE MSD FRAMESET ref H.pdb
#
# difference of distances for proton transfers
#
DISTANCE LIST 53 703 DIFFDIST 703 702
DISTANCE LIST 702 913 DIFFDIST 913 912
#
# steer plan read from file myplan
#
STEERPLAN myplan
ENDMETA

The steerplan file myplan contains the control sequence of the steerplan
and it is composed as follows. The first column contains the time (in timeunits of the program) at which a particular action is planned. The following
columns are composed in blocks of five. Each block contains: the CV keyword, the index of the CV to be used, the spring constant at a given time,
the value of the center of the spring potential and a keyword among CENTRAL
(the potential is a normal parabolic shape), POSITIVE (the potential with
parabolic shape is applied only when the CV value is higher than the center
of the spring) and NEGATIVE (the potential with parabolic shape is applied
only when the CV value is lower than the center of the spring). A Block
of four columns without the latter keyword is also accepted. In this case
CENTRAL is assumed. Wildcards (*) are also accepted for the position. Their
meaning is much clearer in the following example. Comments are allowed
(#) .

Example.
An example of steerplan file.
# Bring CV 4 from wherever it is (*) to -0.5 at 2 ps by increasing
# gradually the spring constant from 0 to 800. CV 5 does the same and goes to -0.64.
0.00
CV
4
000.0
*
POSITIVE
CV
5
000.0
*
POSITIVE
2000.00 CV 4
800.0 -0.5
POSITIVE
CV
5
800.0
-0.64
POSITIVE
# Now in 400 fs keep CV 4 at its value while releasing the other
# slowly to 0 spring constant.
2400.00 CV 4
800.0 -0.5
POSITIVE
CV
5
0.0
-0.64
POSITIVE
# Now drag back only CV 4 to 1.4 with a value that acts only
# on the negative part. The potential on CV 5 is off.
8400.00 CV 4
800.0 1.4
NEGATIVE
CV
5
0.0
-0.64
POSITIVE

Please note that the wildcards have particular meaning:
Example.
Here the dragging force starts from whatever value the system assumes at 0 fs and drags it to -0.5 at 2
ps. The spring constant is linearly increasing from 0.0 to 800.0. When the starting point is a wildcard this
takes the current value.
0.00
2000.00

CV
CV

4
4

000.0
800.0

*
-0.5

POSITIVE
POSITIVE

Here the center of the harmonic potential follows the system at each step from 0 to 2 ps. When the
position of the ending point is a wildcard then the potential will follow the coordinates (this means that the
potential is not applied as the force is zero everywhere).
0.00
2000.00

CV
CV

4
4

800.0
800.0

*
*

POSITIVE
POSITIVE

This case is considered identical as the one before

0.00
2000.00

3.9

CV
CV

4
4

800.0
800.0

0.5
*

POSITIVE
POSITIVE

Adiabatic Bias MD

PLUMED can be used to evolve a system towards a target value in CV space
using an harmonic potential moving with the thermal fluctuations of the
CV[20, 32, 33].
The directive ABMD activates the biasing on the collective variable specified
by the keyword CV. The target value is determined by the keyword TO, the
48

spring constant by KAPPA. The biasing potential is implemented as
(

V (ρ(t)) =

(ρ(t) − ρm (t))2 , ρ(t) > ρm (t)
0,
ρ(t) ≤ ρm (t),
α
2

(3.6)

where
ρ(t) = (CV (t) − T O)2

(3.7)

ρm (t) = min ρ(τ ).

(3.8)

and
0≤τ ≤t

The method is based on the introduction of a biasing potential which is
zero when the system is moving towards the desired arrival point and which
damps the fluctuations when the system attempts at moving in the opposite
direction. As in the case of the ratchet and pawl system, propelled by thermal
motion of the solvent molecules, the biasing potential does not exert work
on the system.
The keyword RESTART can be used when restarting an adiabatic bias MD
calculation to append the value of the CVs on the COLVAR file and to set the
best value previously reached.
Example.
The following input file defines an adiabatic biased MD on the angle CV to a target value of 3.0 rad.
ANGLE LIST 13 15 17
ABMD CV 1 TO 3.0 KAPPA 500.0
PRINT W STRIDE 100
ENDMETA

The corresponding COLVAR
will look like:
#! FIELDS time cv1 vbias vwall vext XX XX ABMD1
0.000 1.884703657 0.000000000 0.000000000 0.000000000
0.200 2.353724409 0.000000000 0.007812004 0.000000000
0.400 2.433164746 0.000000000 0.000000000 0.000000000
0.600 2.608463835 0.000000000 0.004202352 0.000000000
0.800 2.818935636 0.000000000 0.009463992 0.000000000
1.000 2.831631093 0.000000000 0.147609095 0.000000000
1.200 2.583559043 0.000000000 7.171877688 0.000000000
1.400 2.841442463 0.000000000 0.121111568 0.000000000
1.600 2.794514101 0.000000000 0.382087213 0.000000000
1.800 2.634418597 0.000000000 4.258829096 0.000000000
2.000 2.861303394 0.000000000 0.086113847 0.000000000
where the last column is ρm (t).

ABMD
ABMD
ABMD
ABMD
ABMD
ABMD
ABMD
ABMD
ABMD
ABMD
ABMD

1
1
1
1
1
1
1
1
1
1
1

1.243885933
0.412082147
0.321302206
0.149200641
0.026631583
0.004049192
0.004049192
0.003130352
0.003130352
0.003130352
0.000677239

3.10

External potentials

3.10.1

Walls

The UWALL and LWALL keywords define a wall for the value of the CV s which
limits the region of the phase space accessible during the simulation. The
restraining potential starts acting on the system when the value of the CV is
greater (in the case of UWALL) or lower (in the case of LWALL) than a certain
limit LIMIT minus an offset OFF.
The functional form of this potential is the following:
s − LIMIT + OFF EXP
,
Vwall (s) = KAPPA
(3.9)
EPS
where KAPPA is an energy constant in internal unit of the code, EPS a rescaling
factor and EXP the exponent determining the power law.
By default: EXP = 4, EPS = 1.0, OFF = 0.

Example.
To run a well-tempered metadynamics simulation using as CV the distance between one atom and the
center of mass of a group of atoms and limiting its value below 15 Å, you have to use the following input
file:
HILLS HEIGHT 0.1 W STRIDE 100
WELLTEMPERED SIMTEMP 300 BIASFACTOR 10
PRINT W STRIDE 50
DISTANCE LIST 13 SIGMA 0.35
g1->
17 20 22 30
g1
g1->
2 3
g1
ONIONS HEIGHT 1.0 W STRIDE 1000 WIDTH 1.5
BASINS BASIN TOL 0.2 EXPAND PARAM 0.3 INITIAL SIZE 3.0
CLUSTER RUN FREQ 500000 STORE FREQ 250
TORSION LIST 13 15 17 19
TORSION LIST 15 17 19 21
TORSION LIST 17 19 21 23
TORSION LIST 19 21 23 25
cvlist->
1 2 3
cvlist<-

3.14

Driven Adiabatic Free Energy Dynamics (d-AFED)

The driven adiabatic free energy (d-AFED) algorithm [34] is a variant of
the earlier AFED method [35, 36], which required cumbersome coordinates
transformations. In the d-AFED method, an extra dynamical variable S is
coupled to a collective variable s(r), where r represents the coordinates of
a number of atoms in the system. The coupling is mediated by a potential
energy function with harmonic constant κ,
1
V (S, s(r)) = κ (S − s(r))2 .
(3.12)
2
The dynamics of the S meta-variable is adiabatically decoupled from the
dynamics of the underlying physical system by choosing a large mass mS
m̄, where m̄ is a typical mass of the physical system. Thanks to the adiabatic
separation, a temperature TS > T can be assigned to the S meta-variable.
With this choice of ms and Ts , the physical system will evolve fast at room
temperature T around the instantaneous value of s(r) = S . On the other
hand, S will evolve slowly, but have a temperature large enough to drive the
system over high free energy barriers.
In the limit of κ → ∞, it can be shown that the free energy surface
at temperature T can be recovered from the density ρadb (S) sampled at
temperature TS during the adiabatic d-AFED simulation using

∆G(S) = −kB TS log ρadb (S) .

(3.13)

This result generalizes well to the case where more than one collective variable
is used and ∆G(S) is a multi-dimensional free energy surface. Note that
the d-AFED method is similar to the Temperature Accelerated Molecular
Dynamics (TAMD) method devised by other authors[37].
The d-AFED method requires very efficient thermostatting of the metavariable S. Therefore, in the present implementation, S is coupled to a
Generalized Gaussian Moment Thermostat (GGMT) [38]. The meta-variable
is coupled to two thermostatting variables pη and pζ , with associated masses
Qη and Qζ , respectively. Given a typical time scale τ of the thermostated
system, optimal masses are Qη = kB TS τ 2 and Qζ = 38 (kB TS )3 τ 2 . The order-2
GGMT dynamics for one degree of freedom is
pη
pζ
1 p2S
= V (S, s(r)) −
pS −
kB TS +
pS ,
Qη
Qζ
3 mS
"

ṗS

(3.14)

p2S
− kB TS ,
mS
"
#2
1 p2S
ṗζ =
− (kB TS )2 ,
3 mS
pζ
pS
pη
, η̇ =
, ζ̇ =
.
Ṡ =
mS
Qη
Qζ

ṗη =

(3.15)
(3.16)
(3.17)

If multiple reaction coordinates are used, one separate GGMT thermostat
is associated to each of them. The implemented integrator for the dynamics
above is based on a Trotter decomposition of the corresponding Liouville
operator [34]. The quality of the integration can be monitored using the
quantity HS , which would be conserved if the dynamics of S was decoupled
from the physical system,
p2η
p2ζ
p2S
HS (S, pS , η, pη , ζ, pζ ) =
+ V (S) +
+
+ kB TS (η + ζ) (3.18)
2mS
2Qη 2Qζ
The heat transfer WS from the meta-variable to the physical system can be
calculated,
WS =

Z t
0

dt0 κ (S − s(r))

pS
mS

(3.19)

The effective adiabaticity of the coupling can thus be asserted. In addition,
for each collective variable, the quantity HS −WS should be strictly conserved
and provides a quality check for the simulation.

3.14.1

Input for d-AFED

For each CV, a DAFED directive is used to define the parameters of the corresponding dynamics. On the same line, the number of the CV to which the
directive applies is specified after the keyword CV. The temperature TS in
K is given after keyword TEMPERATURE. The thermostat time constant τ is
given in ps after keyword TAUTHERMO. The mass mS and harmonic constant
κ, are given after the keywords MASS and KAPPA, respectively. The units of κ
and mS depend on the nature of the CV. They should always be such that
κS 2 and mṠ 2 are both in units of energy (kJ/mol= amu nm2 / ps2 ), see the
example below.
In addition, tow optional keywords can be used with the DAFED directive.
First, for periodic CVs such as torsion angles, the S variable should also
61

evolve on a periodic interval. This is specified by the keyword PERIODIC,
followed by two numbers for the lower and upper bounds. The numbers can
be replaced by MINUS PI, PLUS PI, or PLUS 2PI to specify −π, +π, or +2π,
respectively.
The optional keyword JACOBIAN FORCE causes a bias force F = −2kB T /S
to be applied to the dynamics of S. This is useful with distance CVs in order
to counterbalance the effect of the Jacobian factor and sample a more uniform
distribution along the CV.
Example.
The following lines couple CV 1 (a distance in nm) to a meta-variable of mass 105 amu with a harmonic
constant of 106 kJ/mol/nm2 and CV 2 (a unitless number) to a meta-variable with mass 103 amu*nm2
with a harmonic constant of 104 kJ/mol. For both CV, the d-AFED temperature is 600 K and the GGMT
thermostat time constant is 0.2 ps. See text for the optional keywords JACOBIAN FORCE and PERIODIC.
DISTANCE LIST 1 34
TORSION LIST 5 15 29 36
DAFED CV 1 TEMPERATURE 600 MASS 1e5 KAPPA 1e6 TAUTHERMO 0.2 JACOBIAN FORCE
DAFED CV 2 TEMPERATURE 600 MASS 1e3 KAPPA 1e4 TAUTHERMO 0.2 PERIODIC MINUS PI PLUS PI
DAFED CONTROL RESTART checkpoint file WRITE STATE -1 N RESPA 1
PRINT W STRIDE 100
ENDMETA

A separate DAFED CONTROL directive contains general controls for the dAFED simulation. The d-AFED dynamics, including all variables described
in Eqs. (3.14 - 3.17) can be restarted exactly from a previous run using a
checkpoint file. Following the WRITE STATE keyword appears the number of
steps after which a checkpoint file is saved. A value of −1 implies that a
checkpoint file is written only when GROMACS saves its own checkpoint
file, i.e. at regular wall clock time intervals. The checkpoint file is saved in
the current directory with default name DAFED STATE. The optional keyword
RESTART is used to specify the path to the checkpoint file from which to
restart.

3.14.2

Typical output for d-AFED

WIth the d-AFED method, the COLVAR file will contain the following data,
if d collective variables are used:
62

• time step
• value of the collective variable s1 (r)...sd (r)
Then for each of the Sj , j = 1...d, appears a set of 4 columns with :
• the meta-variable S
• the instantaneous temperature of S in K
• the conserved quantity HS , see Eq. 3.18, in kJ/mol
• the work WS from S to the physical system, see Eq. 3.19, in kJ/mol
These fields are labeled sj, T sj, E sj, and W sj, respectively, in the
COLVAR header line, j = 1...d. Additional collective variables can be monitored during a d-AFED run, in which case more columns will appear before
the first d-AFED keyword.
Note that for each extended variable Sj , the quantity HSj + WSj should
be conserved, j = 1...d. Let H be the pseudo-energy of the physical system
including the associated thermostats and barostats. Then the total energy
P
of the simulation, H + dj=1 HSj , should be conserved as well. In addition,
P
considering the physical system only, the quantity H − dj=1 WSj should be
conserved.

Chapter 4
Collective variables
PLUMED contains implementations of a large number of CVs so one can properly describe the processes involved in a wide variety of interesting problems.
In the following chapter we describe all the CVs, which are implemented in
PLUMED together with their analytic first derivative.
In general to instruct PLUMED to use a particular CV a line starting with
the keyword indicating the CV type must be included in the input file. This
keyword should then be followed by the various pieces of CV-specific information required to calculate the variable along with the keywords that tell
PLUMED how this particular CV is to be employed. (N.B. unless stated explicitly these pieces of data can be specified in any order)
For metadynamics the lines defining the variables on which the user desires
there to be a biasing potential should contain the SIGMA keyword. This
keyword should then be followed by the width of the Gaussian hills (in the
units of the CV) on that particular CV. This keyword serves two functions;
namely, it instructs PLUMED to use metadynamics and tells it the widths of
the hills. Obviously the SIGMA keyword is not only required if you are running
metadynamics and is not required with other methods.

Specifying lists of atoms
Most collective variables require the user to specify one or several groups of
atoms in their definition. Whenever a set of atom groups is required, the LIST
keyword must be used. This keyword is then followed a set of tags which
specify the number and order of the groups of atoms to be employed. The
atoms involved in each of the groups invoked must be specified somewhere in
64

the input file. These group specifications must then start with the the name
of the group followed by the -> sign and finish with the same name followed
by the <- sign. Between these two delimiters the indices of the atoms which
comprise the group must then be listed, separated by spaces or line feeds.
Example.
The following syntax instructs PLUMED to use the distance between the center of mass of atoms 6 and 10
and the center of mass of atoms 8, 15 and 21 as a CV:
DISTANCE LIST SIGMA 1.0
g1->
6 10
g1
8 15 21
g2<-

Inside the group definition blocks, one can either specify the atom numbers explicitly, or one can use the LOOP keyword to define a regular sequence
of indices with a given starting number, end number and stride.
Example.
The following two commands are equivalent definitions of the group g1:
g1->
10 12 14 16 18 20
g1
LOOP 10 20 2
g1<-

For the special, and rather common, case of a group composed of a single
atom the user can specify the number of the atom of interest rather than the
corresponding tag.

Example.
The following two commands are equivalent ways instruct PLUMED to use the distance between atom 5
and group as a CV:
DISTANCE LIST SIGMA 1.0
g1->
10 12 14 16 18 20
g1
5
g2 5 SIGMA 1.0
g1->
10 12 14 16 18 20
g1<-

4.1

Absolute position

The POSITION keyword instructs PLUMED to use the absolute position of an
atom or a group of atoms, specified by using the LIST keyword, as a CV.
This CV accepts several options that allow the user to restrict the bias to
a given direction, e.g. z, to bias the position of the particle as projected on
a selected line segment or, in analogy with the path CV, to bias the atoms’
distance from a line segment. The keyword DIR accepts as input X, Y or Z
and limits the restraint to the chosen direction.
Example.
The following line instructs PLUMED to use the y coordinate of atom 13 as a CV.
POSITION LIST 13 SIGMA 0.35 DIR Y

The keyword LINE POS instructs PLUMED to use the projection of the
atoms position on a line as a CV, while the keyword LINE DIST instructs
PLUMED to use the distance from the line as a CV. In both these cases the
line is defined by stating its start and end points. The line can then either
be in constrained to be in the XY, XZ or YZ planes (keywords (XY, XZ or
YZ respectively) or it can have an arbitrary orientation in space (keyword
XYZ). Obviously if the line is constrained to the XY, XZ or YZ planes then
the start and end points are two dimensional vectors whereas if it has an
arbitrary orientation these vectors have three components.
66

Example.
The following lines instruct PLUMED to use the projection and distance of the coordinates of atom 13 on
a generic line segment defined by the start and end points (0,0,0) and (2,3,4) respectively as the two
collective coordinates.
POSITION LIST 13 SIGMA 0.35 LINE POS XYZ 0.
POSITION LIST 13 SIGMA 0.35 LINE DIST XYZ 0.

4.2

0.
0.

2.
2.

3.
3.

4.
4.

Distance

The DISTANCE keyword instructs PLUMED to use the distance between the
center of mass of two groups of atoms as a CV. Two groups must be defined
using the LIST keyword and the syntax described in section 4.
Example.
The following lines instruct PLUMED to use the distance between atom number 13 and the center of mass
of the four atoms in list as a CV.
DISTANCE LIST 13 SIGMA 0.35
g1->
17 20 22 30
g1<-

The optional flag NOPBC can be used to calculate the distance without
applying periodic boundary conditions. This should be done only if all the
atoms in the groups are part of the same molecule. See also Sec. 4.28.
The keyword DIR can be used to calculate the component of this distance
along the cartesian axes (X, Y or Z) and or the component in the planes
(XY, XZ or YZ).
Example.
The following line instruct PLUMED to use the X-component of the distance between atom number 20 and
25 as a CV.
DISTANCE LIST 20 25 DIR X SIGMA 0.35

Whenever one wants to use the difference between two distances one may
use the keyword DIFFDIST. This may turn to be useful in bond breaking/
bond formation.
67

Example.
The following line instruct PLUMED to use the difference between the distance between 20 and 25 and the
distance between 30 and 31 as CV.
DISTANCE LIST 20 25 DIFFDIST 30 31

Groups are also accepted as input instead of two atoms.
Other two useful variants are the following: the distance of a point from
an axis and the projection of the point on the axis. The first is introduced by
the additional keyword POINT FROM AXIS followed by the atom or the group
defining the point respect to which the distance has to be calculated. The
axis is defined through the standard two groups appearing in the definition
of the distance collective variable.
Example.
The following line instruct PLUMED to use the distance between one atom, 26 and the axis defined by the
two atoms 20 and 25 as CV.
DISTANCE LIST 20 25 POINT FROM AXIS 26

In a similar way it is possible to calculate the projection of this point on
an axis to be used as a CV by using the keyword PROJ ON AXIS.
Example.
The following line instruct PLUMED to use the projection of the coordinate if atom 26 on the axis defined by
the two atoms 20 and 25 as CV.
DISTANCE LIST 20 25 PROJ ON AXIS 26

Similarly to all the other variables, these two keywords may accept groups
instead of atom indexes.

4.3

Minimum distance

The MINDIST keyword instructs PLUMED to use the minimum distance between
two groups of atoms as a CV. To ensure differentiability, this quantity is
implemented as:
β
s=
,
P
log ij exp(β/||rij ||)
68

where by default β = 500. The value of β can however be tuned if needed
by using the optional keyword BETA. Much like the distance variable when
calculating the minimum distance one must define two groups using the LIST
keyword and the syntax described in section 4.
Example.
The following lines instruct PLUMED to use the minimum distance between atom number 13 and the set of
atoms in list as a CV.
MINDIST LIST 13 SIGMA 0.35 BETA 500.
g1->
17 20 22 30
g1<-

The optional flag NOPBC can be used to calculate the distance without
applying periodic boundary conditions. This should be only be done if all
the atoms in the groups are part of the same molecule. See also Sec. 4.28.

4.4

Angles

The ANGLE keyword instructs PLUMED to use the angle defined by the centers
of mass of three groups of atoms as a CV. The compulsory LIST keyword
must be followed by three properly defined groups (see Section 4).
Example.
The following lines instruct PLUMED to use the angle defined by atom number 102 and the centers of mass
of the atoms in groups g1 and g2 as a CV.
ANGLE LIST 102 SIGMA 0.05
g1->
13 15
g1
LOOP 1000 3000 3
g2<-

It is also possible to use the sine or cosine of the angle as a collective
coordinate by including the SIN or COS keywords respectively similarly to
what is done for TORSION.
69

4.5

Torsion

The TORSION keyword instructs PLUMED to use a dihedral angle as the CV.
This angle can either be defined by four atoms or, more generally, by the
positions of the centers of mass of four groups of atoms. The compulsory
LIST keyword must be followed by four properly defined groups (see Section
4).
Example.
The following lines instruct PLUMED to use to the torsion angle about the centers of mass of the four groups
, , , as a CV.
TORSION LIST SIGMA 0.35

It is also possible to use the sine or cosine of the torsional angle as a
collective coordinate by including the SIN or COS keywords respectively.
Example.
The following line instructs PLUMED to use the cosine of the torsion angle about the centers of mass of the
four groups , , , as a CV.
TORSION LIST COS SIGMA 0.35

4.6

Coordination number

The COORD keyword instructs PLUMED to use the total number of contacts between the atoms in group G1 and those in group G2 - the coordination number
between these two groups. To ensure differentiability, this is implemented as
the sum:
X X
s=
sij ,
i∈G1 j∈G2

where this sum is extended to all pairs of atoms with i ∈ G1 and j ∈ G2 . The
individual contributions sij are defined using a switching function, which, in
the present case, is given by:
sij =







1
rij n
)
r0
rij
1−( r )m
0

1−(

for rij ≤ 0
for rij > 0

where rij = |ri − rj | − d0 . The user must supply the r0 , d0 , n and m parameters, using the additional keywords R 0, D 0, NN and MM respectively and thus
has a great deal of control over the definition of the switching function. In
general a good first guess for these parameters can be achieved by looking at
the pair distribution function and setting d0 equal to the position of the first
peak in the pair distribution function, r0 as the full width at half maximum
of the peak and n and m to force sij ' 0 at the first minimum of the pair
distribution function. However, oftentimes different choices for these parameters will lead to better results because of certain specific properties of the
system of interest. An optional keyword PAIR treats the atoms in a pairwise
fashion so that instead of simply counting the number of bonds between two
groups of atoms one can define which precise bonds between the two groups
should be monitored. In this case the groups, and must have
the same number of atoms as the switching functions are on the distances
between the ith atom of group and the ith atom of group .
Example.
The following lines instruct PLUMED to use the coordination of the atoms in group g1 – 13 and 15 – with
the atoms in group solvent as the CV.
COORD LIST NN 6 MM 12 D 0 2.5 R 0 0.5 SIGMA 0.35
g1->
13 15
g1
LOOP 1000 3000 3
solvent<-

The optional flag NOPBC can be used to calculate the distance without
applying periodic boundary conditions. This should only be done when all
the atoms are part of the same molecule. See also Sec. 4.28.

4.7

Hydrogen bonds

HBONDS is the keyword for a variable that counts the number of intra-molecular
hydrogen bonds between a group of hydrogen bond donors and a group of
hydrogen bond acceptors. This is defined as:
s=

1 − ( drij0 )n

1 − ( drij0 )m
71

where i ∈ D is the group of donors and j ∈ A are the acceptors.
The two groups must be defined using the compulsory LIST keyword followed by two groups (see Section 4). PLUMED then assumes that there is only
one donor/acceptor per residue and, that within the list, the donor/acceptor
atoms on neighboring residues are consecutive. The values of r0 , n and m
can be specified using the R 0, NN and MM keywords. If no value is given for
r0 , n and m the default values of r0 = 2.5, n = 6 and m = 12 are assumed.
The TYPE keyword selects which residues to include in the count:
• With TYPE 0, all donor-acceptor pairs are included;
• If TYPE 1 is specified only those donor-acceptor pairs separated by an
odd number of residues greater than 4 are counted. This allows one to
monitor parallel β-sheet formations.
• If TYPE 2 is specified only those donor-acceptor pairs that are separated
by exactly 4 residues are included. This allows the formation of αhelical conformations (α type) to be monitored;
• If TYPE 3 is specified, only those donor-acceptor pairs that are separated by an even number of residues greater than 4 are counted. This
allows anti-parallel β-sheet formations (β-even type) to be monitored.
• If TYPE 4 is specified, only the first donor and the first acceptor and so
on are counted. This allows to monitor a set of native hydrogen bonds.
• If TYPE 5 is specified, only hydrogen bonds between atoms belonging to
different residues are counted. Moreover, pairs with an index difference
less than 5 are also discarded, so as to avoid counting H and a O that
are in the same peptide group (NH-C=O) (contributed by M. Cuendet).
This option requires the user to specify to which residue each atom is
belonging using the RESLIST keyword (see example below).

Example.
The following lines instruct PLUMED to use the count of the total number of hydrogen bonds between the
pairs in groups H and O as a CV. The default switching function with r0 = 2.5, n = 6 and m = 12 is
implied.
HBONDS LIST TYPE 0 SIGMA 0.1
H->
6 10
H
8 12
O TYPE 0 SIGMA 0.1 NN 8 MM 20 R 0 2.5

In the following example, inter-residue hydrogen bond are counted

HBONDS LIST RESLIST TYPE 5
H-> 1 9 17 H 5 13 21 O 1 2 3 resH 1 2 3 resO<-

The optional flag NOPBC can be used to calculate the distance without
applying periodic boundary conditions. This should be done only when the
atoms are part of the same molecule. See also Sec. 4.28.

4.8

Interfacial water

WATERBRIDGE is the keyword for a CV that counts the number of interfacial
contacts. This variable does this by calculating the number of atoms from
group G0 that are simultaneously in contact with atoms from both groups
G1 and G2 . A typical application of this CV is to count the number of water
molecules at the interface of two surfaces. This is calculated using:
sWatBr =

n0
X



n1
X

j| n
)
1 − ( |rir−r
0




j| m
1 − ( |rir−r
)
0


i

n2
X

j| n
1 − ( |rir−r
)
0

j| m
1 − ( |rir−r
)
0


.

The syntax of the command requires the user to specify three groups of atoms
after the keyword LIST starting with the G1 and G2 groups and closing with
73

the G0 group. (see Section 4). The parameters of the switching function are
then defined using the usual NN, MM and R 0 keywords.
Example.
The following command instructs PLUMED to use the number of atoms in the solvent group that are
simultaneously in contact with either atom 6 or 10 and either atom 8, 15 or 21 as a CV.
WATERBRIDGE LIST NN 8 MM 12 R 0 4.0 SIGMA 0.1
type1->
6 10
type1
8 15 21
type2
LOOP 100 1000 3
solvent<-

4.9

Radius of gyration

One can employ the radius of gyration of a group of atoms defined with the
compulsory additional keyword LIST by using the RGYR directive. The LIST
keyword (see Section 4) must be followed by only one properly defined group.
This CV is calculated using :
sGyr =

2
Pn |r − r
i
COM | 1/2
i
,
Pn

mi
where the sums are over the n atoms in group G and the center of mass is
defined using:
Pn
i ri m i
.
rCOM = P
n
i mi
i

Example.
The following lines instruct PLUMED to use the radius of gyration of the group as a CV.
RGYR LIST SIGMA 0.35

N.B. The radius of gyration is calculated without applying periodic boundary conditions so the atoms in group should all be part of the same
molecule. See also Sec. 4.28.
74

4.9.1

Gyration tensor based CVs

The RGYR directive can be used also to access to a number of CVs based on
gyration tensor [39]:
S=

1
N

 P 2
xi
 P
xi yi

P

xi zi

P
xy
P i2 i
y
P i

yi zi


P
xi zi
P
yz 
,
P i 2i

(4.1)

which describes the spatial distribution of mass in a molecule or complex of molecules. Alternatively, tensor of inertia may be used for the same
purpose:
I=

 P
P
P
mi (yi2 + zi2 )
−mi xi zi
−mi xi yi
P
P
P

mi (x2i + zi2 )
−mi xi yi
−mi xi zi

P
P
P
2
2

−mi xi zi

−mi yi zi

mi (xi + yi )



.

(4.2)

The weighting of atomic contribution is controlled by keyword MASS-WEIGHTED
(default) and NO-WEIGHTED. When MASS-WEIGHTED directive is applied, the
mi N
individual atomic contribution to the gyration tensor is weighted by P
mi
and the center of mass is used as origin of the system. If NO-WEIGHTED option is chosen, coordinates are related to the geometrical center of object and
no mass-weighting is performed in calculation of CV.
Diagonalizing of gyration tensor provides the principal moments of gyration - S1 ,S2 ,S3 and three eigenvectors corresponding to the principal axes of
inertia. The individual gyration tensor based CVs are available by specifying
the directive RGYR and the keyword TYPE followed by one of the following:
• TRACE: trace of inertia tensor
• GTPC 1:
√ the largest eigenvalue of gyration tensor. The square root
0
S1 = S1 of principal moment S1 is used as CV. If we approximate
0
the shape and mass distribution of the object with an ellipsoid, S1
corresponds to the largest radius of this ellipsoid.
√
0
• GTPC 2: the middle principal moment of gyration tensor, S2 = S2 .
√
0
• GTPC 3: the smallest principal moment of gyration tensor, S3 = S3 .
• RGYR 1:
√ the largest radius of gyration around principal axes of inertia,
rg1 = S1 + S2 .
75

• RGYR 2:
√ the middle radius of gyration around principal axes of inertia,
rg2 = S1 + S3 .
• RGYR 3:
√ the smallest radius of gyration around principal axes of inertia,
rg3 = S2 + S3 .
• ASPHERICITY: the deviation of mass distribution from spherical symmetry. The
√ modified version of asphericity was implemented as a
0
CV: b = b, where original Suter’s [40] asphericity is given by b =
S1 − 21 (S2 + S3 ).
• ACYLINDRICITY: the deviation
from cylindrical symmetry. Again the
√
0
modified version c = c is used as a CV rather then original form [40]
c = S2 − S3 .
S2 +S1 S3 +S2 S3
• KAPPA2: the relative shape anisotropy[40] κ2 = 1 − 3 S1(S
re2
1 +S2 +S3 )
flects both symmetry and dimensionality of an object. Limited between
values of 0 and 1, κ2 reaches 1 for an ideal linear arrangement and drops
to zero in case of highly symmetric configurations (at least tetrahedral
symmetry).

Example.
The following lines instruct PLUMED to use the trace of the inertia tensor of the group as a CV.
RGYR TYPE INERTIA LIST SIGMA 0.35

Example.
The following lines instruct PLUMED to use the asphericity of the group as a CV.
RGYR TYPE ASPHERICITY LIST SIGMA 0.35

4.10

Dipole

DIPOLE instructs PLUMED to use the electrical dipole generated by a group of
atoms as a CV:
n
sdipole = |

ri qi |.

Only the LIST keyword followed by one properly defined group is required
to define this CV (see Section 4).
Example.
The following lines instruct PLUMED to use the dipole of the group as a CV.
DIPOLE LIST SIGMA 0.35

4.11

Dihedral correlation

DIHCOR is the keyword for a CV that measures the similarity between adjacent
dihedral angles:
ND
X
1
sDC =
(1 + cos (φi − φi−1 )) .
i=2 2
The syntax for this CV requires the user to specify the number of dihedrals
ND and, in the subsequent ND lines, the indices of the four atoms defining
each dihedral φi .
Example.
The following lines instruct PLUMED to use the dihedral correlation for the 3 dihedrals listed as a CV.
DIHCOR NDIH
168 170 172
170 172 188
172 188 190

3 SIGMA 0.1
188
190
197

4.12

Alpha-beta similarity

ALPHABETA is a keyword for a CV that measures the similarity of dihedral
angles to a reference value (see also ??). It is calculated using:
sαβ =

1X
1 + cos φi − φRef
.
i
2 i=1

The syntax for this CV requires the user to specify the number of dihedrals
ND and, in the subsequent ND lines, the indices of the four atoms defining
each dihedral followed by the value of the dihedral in the reference conformation φRef
i .
Example.
The following lines instruct PLUMED to use the Alpha-beta similarity of three dihedrals as a CV.
ALPHABETA NDIH 3 SIGMA 0.1
168 170 172 188 3.14
170 172 188 190 .56
188 190 192 230 3.14

4.13

Alpharmsd

ALPHARMSD is the keyword for a CV that counts the number of 6-residue
segments in the protein chain that resemble an ideal alpha helix (i.e. the
average experimental structure). This CV is calculated using:
sα rmsd =

X
α

n RMSD {Ri }i∈Ωα ,

where n is the coordination (switching) function defined in Sec. 4.6. The sum
over α in the above runs over all the 6-residue segments in the protein chain,
where each of the residues is defined based on the positions of the backbone
atoms N, CA, C, O and CB, Ωα . The distance used in the switching function
is then the root mean square difference between the distance matrix of the
atoms in the set Ωα with the corresponding distances between these atoms
in the ideal alpha helix {R0 }:

RMSD {Ri }i∈Ωα ,

=
78

v
u
u
t

1
Npairs

X
i,j∈Ωα

dij − d0ij .

See Ref.[41] for more details (although note that in PLUMED the RMSD between distance matrices is used rather than the cartesian RMSD. However,
these two measures are essentially equivalent).
To use this CV the syntax requires the user to specify the coordination
function parameters R 0, D 0, NN and MM (see Sec. 4.6) (we suggest values of R 0=0.8 Angstrom, NN=8, MM=12, D 0=0), and a conversion factor
ANGSTROM SCALE which converts the ideal alpha positions from Angstrom to
the length units of the MD code (e.g. in Gromacs units are nanometers
therefore ANGSTROM SCALE is 0.1). In addition a list of list of atom indices for
the N, CA, C, O and CB (in this order) for each residue must be provided
for all the consecutive residues (in ascending order) which form the chain.
For those residues, such as glycine, which do not have a CB the atom index
for the corresponding hydrogen should be used.
Important note: for those MD codes which, like Gromacs4, do not keep
the protein whole but instead split it by the PBC, ALPHARMSD requires the
additional option NOPBC, together with the ALIGN ATOMS command for all the
atoms employed in the ALPHARMSD.
Example.
The following lines define an Alpharmsd CV for all 16 consecutive residues of a protein in the MD code
Gromacs4 (which uses nanometers as the length unit).
ALPHARMSD LIST SIGMA 0.5 R 0 0.08 NN 8 MM 12 ANGSTROM SCALE 0.1 NOPBC
ncacocb->
1 5 23 24 7
25 27 42 43 29
44 54 56 57 51
58 68 70 71 65
72 74 77 78 76
79 81 101 102 83
103 105 116 117 107
118 120 138 139 122
140 142 162 163 144
164 166 179 180 168
181 183 190 191 185
192 194 214 215 196
216 218 225 226 220
227 229 236 237 231
238 240 243 244 242
245 247 267 268 249
ncacocb->

4.14

Antibetarmsd

ANTIBETARMSD counts the number of pairs of 3-residue segments in the protein chain which are similar to the ideal antiparallel beta (i.e. the average
experimental structure). See Ref.[41] for details (but remember that here the
RMSD among distance matrices is used instead of the cartesian RMSD as the
two are basically equivalent). The definition, input parameters and syntax
for this CV are the same as for the ALPHARMSD and again we suggest values
of R 0=0.8 Angstrom, NN=8, MM=12, D 0=0 for the parameters. The only
difference in the implementation is the additional option STRANDS CUTOFF
which allows the user to specify a threshold distance beyond which pairs of
3-residue segments are considered far. This option considerably speeds up
the computation (often 1 nm is good choice for this quantity).
Important note: for those MD codes which, like Gromacs4, do not keep
the protein whole but instead split it by the PBC, ANTIBETARMSD requires
the additional option NOPBC, together with the ALIGN ATOMS command for
all the atoms employed in the ANTIBETARMSD.
Example.
The following lines define an Antibetarmsd CV for all 16 consecutive residues of a protein in the MD code
Gromacs4 (which uses nanometers as the length unit).
ANTIBETARMSD LIST SIGMA 0.5 R 0 0.08 NN 8 MM 12 ANGSTROM SCALE 0.1
STRANDS CUTOFF 1. NOPBC
ncacocb->
1 5 23 24 7
25 27 42 43 29
44 54 56 57 51
58 68 70 71 65
72 74 77 78 76
79 81 101 102 83
103 105 116 117 107
118 120 138 139 122
140 142 162 163 144
164 166 179 180 168
181 183 190 191 185
192 194 214 215 196
216 218 225 226 220
227 229 236 237 231
238 240 243 244 242
245 247 267 268 249
ncacocb->

4.15

Parabetarmsd

PARABETARMSD counts the number of pairs of 3-residue segments in the protein
chain which are similar to the ideal parallel beta (i.e. the average experimental structure). See Ref.[41] for details (but remember that here the RMSD
among distance matrices is used instead of the cartesian RMSD as the two are
basically equivalent). The definition, input parameters, and syntax are the
same as for theALPHARMSD and again we suggest values of R 0=0.8 Angstrom,
NN=8, MM=12, D 0=0 for the parameters. The only difference in the implementation is the additional option STRANDS CUTOFF which allows the user to
specify a threshold distance beyond which pairs of 3-residue segments are
considered far. This option considerably speeds up the computation (often
1 nm is good choice for this quantity).
Important note: for those MD codes which, like Gromacs4, do not keep
the protein whole but instead split it by the PBC, PARABETARMSD requires
the additional option NOPBC, together with the ALIGN ATOMS command for
all the atoms employed in the PARABETARMSD.
Example.
The following lines define a Parabetarmsd CV for all 16 consecutive residues of a protein in the MD code
Gromacs4 (which uses nanometers as the length unit).
PARABETARMSD LIST SIGMA 0.5 R 0 0.08 NN 8 MM 12 ANGSTROM SCALE 0.1
STRANDS CUTOFF 1. NOPBC
ncacocb->
1 5 23 24 7
25 27 42 43 29
44 54 56 57 51
58 68 70 71 65
72 74 77 78 76
79 81 101 102 83
103 105 116 117 107
118 120 138 139 122
140 142 162 163 144
164 166 179 180 168
181 183 190 191 185
192 194 214 215 196
216 218 225 226 220
227 229 236 237 231
238 240 243 244 242
245 247 267 268 249
ncacocb->

4.16

Electrostatic potential

Using the ELSTPOT keyword one can instruct PLUMED to use the electrostatic
potential exerted by a group of atoms on the center of mass of a second group
of atoms (or single atom) as a CV. This CV is calculated using:
sELST =

NA
X
i

qi
∗ f (|ri − rCOM |, R0 , CU T )
|ri − rCOM |

where

PNB

ri mi
i
rCOM = P
.
NB
mi
i
Here the sum in the first equation above is over the NA atoms in the group
whose charges exert the electric potential, while in the second equation the
sum is over the NB atoms in the group that defines the point at which this
potential is felt. f (x) is a smoothing function defined as:





1.0
x < R0
πx
f (x, R0 , CU T ) = cos 2(CU T −R0 ) R0 ≤ x ≤ CU T


 0
x > CU T
where R0 is the onset and CU T is a cutoff distance.
Example.
The following lines instruct PLUMED to use the electrostatic potential exerted by the atoms in group2 on the
center of mass of the atoms in group1 as a CV:
ELSTPOT LIST R 0 4.0 CUT 12.0 SIGMA 0.01
group1->
LOOP 1 8 1
group1
LOOP 9 16 1
group2<-

4.17

Puckering coordinates

PUCKERING refers to the set of collective coordinates for 6-membered rings in
polar coordinates [42]. Given the coordinates zj , that represent the displacements of the j−th atom from the mean ring plane, three variables Q, θ, φ can
be obtained starting from the general definition for 6-membered rings
82

q2 cos φ2 =

6
2π
1X
zj cos
(j − 1)
3 j=1
3

6
2π
1X
zj sin
q2 sin φ2 = −
(j − 1)
3 j=1
3

q3 =
as
Q=

6
1X
(−1)j−1 zj
6 j=1

q22 + q32

θ = arctan (q2 /q3 )

≥0
∈ [0, π]

∈ [0, 2π)

φ = φ2

Each one can be selected as a TYPE of PUCKERING. The CV accept the
LIST and SIGMA keywords. As of now LIST accept only 6 atoms. The atoms
in the list have to be enumerated following the chemical sequence (although
the numbering scheme in the topology does not have to be sequential). In
order to fulfill IUPAC convention for sugar hexopyranose rings the first atom
in the list has to be the ring oxygen, followed by the anomeric carbon. Q is
in general a fast degree of freedom.
Example.
The following lines define a PUCKERING CV:
PUCKERING LIST TYPE PHI SIGMA 0.1
group1->
1 2 3 4 5 6
group1<-

4.18

Path collective variables

One can instruct PLUMED to use path collective variables [43] using the S PATH
and Z PATH keywords. In this scheme one defines a path as a set of N reference
conformations that define the path in configuration space X from some initial

state to some final state. The s variable (defined with S PATH) then measures
the position along the path, and is defined as:
s=Z

−1

N
X

ie−λd(Xi ,X(t)) ,

i=1

where X(t) is the configuration of the system at any given time, d : X ×X →
PN −λd(Xi ,X(t))
R+
is a normalization factor in
0 is a metric on X , and Z =
i=1 e
which the prefactor, λ, should be chosen so as to have λd(Xi , Xi±1 ) ' 2.3 on
average.
By contrast the z variable (defined with Z PATH) measures the position
off the path, and is defined as:
z = −λ−1 log Z.
For both S PATH and Z PATH the following keywords must be used:
• TYPE, which defines the metric used to calculate distances in configuration space. The following sections provide more details on this but,
suffice to say, currently one can use either MSD (mean square deviation
see section 4.18.1), DMSD (distance root mean squared deviation see section 4.18.2) or CMAP (the distance between contact matrices see section
4.18.3) [44].
• NFRAMES which sets the number of reference structures, N , that are
used in the definition of the path.
• LAMBDA, the prefactor λ in the exponential term of the equation that
defines both s and z.
• Optionally, you can also use NEIGHLIST which defines a neighbor list on
the closest frames to speed up the calculation followed by the number of
steps in which the list must be calculated and the number of elements
it must contain (e.g. NEIGHLIST 50 10 means that each 50 steps the
neighbour list must be calculated and it will contain the 10 closest
elements to the current molecular dynamic snapshots. All the other
are discarded up to the next neighbor list calculation.).
Other keywords are specific to the type of path variable being defined
(i.e. mean square displacement in Cartesian coordinates, MSD in distances
or contact map distance).
84

An important change occurred from the version 1.2 to 1.3 : the convention for the names has changed the former RMSD has been replaced by MSD.
Similarly DRMS has been changed into DMSD. This was done since many users
noticed the inconsistencybetween the definitions and the related keywords.
N.B. before using this feature please ensure that the parameters in metadyn.h
(common files directory) are set properly. In particular look for the section
containing:
// path dimensions
#define MAXATOMS PATH 230
#define NMAX PATH 8
#define MAXFRAMES PATH 22
#define MAXATOMS RMSD 230
#define MAXCHARS PATH 40
// cmap
#define MAXDIM CMAP 3800
#define MAXNUM GROUP 10
#define MAXATOM GROUP 30
and ensure that MAXATOMS PATH is greater than or equal to the number
of atoms per frame involved in your path and MAXFRAMES PATH is greater
than or equal to the number of frames you are using to define your path. In
addition NMAX PATH maximum number of path variables that you can use
simultaneously. If you change any of these values you must subsequently
recompile all the MD codes in which you have implemented PLUMED in order
for your changes to take effect.

4.18.1

Mean square deviation

As already briefly mentioned, when using the path collective coordinates
S PATH and Z PATH, the command MSD after the TYPE keyword instructs
PLUMED to the the mean square deviation of a subset of the atoms in the
system (the displacement set B) calculated after the system has been aligned
to another subset of the atoms (the alignment set A) as the metric used in
the definition of the path. This quantity is calculated using
d(Xj , Xi ) =

NB
X

wj
(Xa(j) − Mij ({z})Xa(i) )2 ,
w
a=1 T OT
85

where Mij ({z}) is the roto-translation matrix calculated using the Kearsley [45] algorithm and wj is the number specified (called the displacement
parameter) for the beta column in the provided input PDB. If this is one
P
then wj =1 and wT OT =N (while differently wT OT = N
i wi ). The weights for
the alignment are actually taken into account when forming the alignment
matrix that compose the Mij ({z}) and are denoted here with the set {z}
defined by the beta column in the reference PDB.
To use path CVs the user must specify the coordinates of the atoms in
each of the reference frames of the path in a set of supplementary input files.
The user should specify the basename of these files (¡basename¿) after the
keyword FRAMESET and be aware that PLUMED then will expect to find N
files named ¡basename¿.i.pdb, which contain the coordinates of the atoms
in both the displacement and alignment sets for the ith frame in the path.
The particular set/s each atom is involved in is specified using the the last
two numerical fields of the frameset file. Values of 1.0 and 0.0 indicate that
the atom is to be used in the alignment set only, while values of 0.0 and 1.0
indicates that the atom is to be used in the displacement set only. Values of
1.0 and 1.0 indicate that the atom is to be used in both the alignment and
displacement sets.
Version 1.1 and higher allow these alignment and displacement indicators to
be non-integer. This allows users to perform a weighted alignment in cases
where the alignment of one region of system is considered more important
and a weighted displacement evaluation in cases where the displacement of a
particular atom/s is though to be of particular importance. Be aware however that this feature can produce strange numbers that are not trivial to
interpret and is thus perhaps best left to experienced users.
N.B. PLUMED contains a hard coded limit on the number of atoms that can
be used in the alignment so, like in the previous section, users should ensure
that the value of this hard limit (MAXATOMS RMSD in metadyn.h) is sufficiently
large for their needs prior to compilation.
The unit of distance in the PDB files is Ångstrom. Engines like GROMACS,
whose internal units are different, will perform appropriate conversions automatically.
Clearly, the atom indices in these PDB files must be the same as the indices
of the atoms they refer to in the system topology. It is therefore likely that
the atom indices in the frameset files may be non-consecutive.
Additional keywords are supported: NO ROT and NO CENTER. These two keywords prevent the rotation and the center of mass alignment respectively
86

whenever one uses MSD.

Example.
The following command defines a path using MSD metrics. 2 frames are used to define the path, which
has λ = 9.0.
S PATH TYPE MSD FRAMESET frame NFRAMES 2 LAMBDA 9.0 SIGMA 0.1
Z PATH TYPE MSD FRAMESET frame NFRAMES 2 LAMBDA 9.0 SIGMA 0.1
Two PDB files must be provided: frame 1.pdb and frame 2.pdb. The last two columns of these files specify
which atoms are to be used for alignment and which are to be used to calculate the MSD.
ATOM
ATOM
ATOM
ATOM
ATOM
END

1 C ALA 2 -0.186 -1.490 -0.181 1.00 0.00
2 O ALA 2 -0.926 -2.447 -0.497 1.00 1.00
15 N ALA 2 0.756 0.780 -0.955 1.00 0.00
17 CA ALA 2 0.634 -0.653 -1.283 1.00 1.00
19 CB ALA 2 2.063 -1.233 -1.286 1.00 1.00

4.18.2

Distance mean square deviation

Instead of using the MSD in the definition of the metric for the path in the
S PATH and Z PATH CVs the command DMSD after the TYPE keyword allows
one to use a metric based on the mean square deviation of the distances
between a subset of the atoms in the system:
d(Xj , Xi ) =

NX
NA
A −1 X
2
(j)
(i)
(rab − rab )2 ,
NA (NA − 1) a=1 b=a+1

(j)

where rab is the distance of atoms a and b in the j-th reference frame. For this
metric the coordinates of the atoms in the reference frames of the path are
specified using the keyword FRAMESET along with a set of pdb files containing
the atom coordinates. This is the same way that they are specified when the
root mean square deviation metric is used (see section 4.18.1).

4.18.3

Contact map distances

The final choice one can employ to define the metric for the path in the
S PATH and Z PATH CVs is to use the command CMAP after the TYPE keyword.

This sets the metric for the path to be the distance between the contact
matrices for a given subset of the atoms in the system:
(j)

(i)

d(Xj , Xi ) = ||Dab − Dab ||.
Given two sets of atoms a, b ∈ J , this contact matrix Dab is calculated using:

Dab (X) = θ(cab − rab )wab

(0)

1 − (rab /rab )nab
1 − (rab /rab )m
ab

where θ(x) is a step function which vanishes if x < 0.
(0)
As with other variables, the parameters rab , nab and mab allow for a
great freedom in the definition of the switching function. What is more the
parameter cab allows one to set the values of the switching function to zero at
large separations. Finally, if the formation of particular contacts is deemed
to be of great importance the weights wab can be used to change their relative
importance.
Like the other metrics for the path collective variables to use path collective variables with a contact map metric information must be provided in
supplementary files about the frames that make up the path. The syntax
requires the user to specify a file name for the indices of the atoms and the
parameters defining the calculation of the contact matrix (after the keyword
(i)
INDEX) and another filename for the values of the reference matrices Dab ,
after the keyword MAP.
The index file, specified after the INDEX keyword, must contain one line for
each of the elements in the contact matrix Dab , which specified how that
particular contact should be calculated. Each of these lines should begin
with the CONTACT keyword followed by a numerical label for the contact.
The next two fields are the indices a, b of the two atoms that make up the
(0)
contact, which are followed by the values of rab , nab , mab , cab and wab in the
switching function.
The values file, specified after the keyword MAP, contains the values of the
(i)
reference matrices Dab used in the definition of the path. Each line in this file
corresponds to one element in the contact matrix. These lines are formatted
with the numerical label for the contact as the first field. This is followed by
the indices of the two atoms that make up the contact a, b and finally the
value of this particular switching function in the reference frame. Unlike the
MSD and DMSD metrics all the reference frames are placed in a single file
with each reference frame separated by the END keyword.
88

The optional flag NOPBC can be used to calculate the distance without
applying periodic boundary conditions. This should be done only if all the
atoms in the groups are part of the same molecule. See also Sec. 4.28.
Example.
The following command defines a path using the contact map metric. 2 frames are used to define the path
with λ set equal to 0.1.
S PATH TYPE CMAP NFRAMES 2 INDEX fr.ndx MAP fr.mps LAMBDA 0.1 SIGMA 1.

NOPBC

The fr.ndx file should contain the details on how each of the switching functions in the contact matrix
should calculated:
CONTACT
CONTACT
CONTACT
CONTACT

1
2
3
4

1 2 3.0 6 10 100.0 1
1 15 3.0 6 10 100.0 1
17 2 3.0 6 10 100.0 1
15 19 3.0 6 10 100.0 1

The fr.mps file contains the values of the reference matrices:
1 1 2 0.99491645
2 1 15 0.76586085
3 17 2 0.79183088
4 15 19 0.81924184
END
1 1 2 0.99369661
2 1 15 0.76748693
3 17 2 0.76454272
4 15 19 0.72917217
END

One can also use contact matrices that involve contacts between the centers of mass of groups of atoms. However, in this case, a further additional
file containing the definition of the groups must be provided. The name of
this file is specified in the PLUMED input file, using the GROUP keyword. Then
within the specified file each group is defined on a single line starting with
the GROUP keyword. This keyword is followed by a numerical label for the
group, the number of atoms in the group and then a list of the indices of
the various atoms that make up the group. To instruct PLUMED to use these
group contacts rather than atomic contacts one must use lines starting with
the keyword GROUP in the index file. The remainder of the format of the
lines in the index file and the format of the corresponding lines in the map
file are identical to the lines used to specify atomic contacts. However, the
indices that would have specified the indices of the atoms involved in the
atomic contact must be replaced with the indices from the group file of the
two groups that make up the contact.
89

Example.
The following command defines a path using the contact matrix metrics. 2 frames are used to define the
path, and a value of λ = 0.1 is set.
S PATH TYPE CMAP NFRAMES 2 INDEX fr.ndx

MAP fr.mps GROUP fr.grp LAMBDA 0.1 SIGMA 1.

The file fr.grp defines three groups of atoms:
GROUP 1 4 23 43 56 457
GROUP 2 5 76 47 97 322 695
GROUP 3 4 17 15 19 2
The fr.ndx file contains the parameters of the contacts between groups:
GROUP 1 1 2 3.0 6 10 100.0 1
GROUP 2 1 3 3.0 6 10 100.0 1
GROUP 3 2 3 3.0 6 10 100.0 1
Finally, the fr.mps has the values of the contact matrix in the reference positions:
1 1
2 1
3 2
END
1 1
2 1
3 2
END

2 0.9949
3 0.7658
3 0.7918
2 0.9936
3 0.6674
3 0.8645

It is possible to have maps in which there are both atomic and group
contacts. However, be aware that, in the index files in which the switching
functions are specified, the definitions of the group contacts MUST follow
the definitions of the atomic CONTACT functions.
Example.
An example of an index file containing both atom-atom contacts and group contacts:
CONTACT 1 123 545 7.0 6 10 100.0 0.50
CONTACT 2 224 244 8.5 6 10 100.0 0.50
GROUP 3 1 2 3.0 6 10 100.0 1.00

4.18.4

Using path variables as MSD, DMSD and CMAP
and the TARGETED statement

If one defines a z path collective variable with a single frame it is clear from
the definition
X
z = −λ−1 log Z = −λ−1 log
e−λdf
f

that this is equivalent to using to the squared distance of the current configuration from a reference structure in the chosen metric. This CV can thus
be used in simulations which employ the standard MSD, distance DMSD or
CMAP distance from a single frame as a collective coordinate. However, in
this case we recommend use of the alias TARGETED instead, which defines the
exact same CV but with a far simpler syntax.
Example.
The following command instructs PLUMED to do steered MD towards a target frame using MSD metrics.
PRINT W STRIDE 10
TARGETED TYPE MSD FRAMESET ref frame.pdb
STEER CV 1 TO 3.0 VEL 0.5 KAPPA 500.0
ENDMETA

A single PDB file must be provided: ref frame.pdb, in which the last two columns specify which of the
atoms are to be used for alignment and which are to be used to calculate the MSD/DMSD (see section
4.18.1).
ATOM
ATOM
ATOM
ATOM
ATOM
END

In the case of CMAP the input for PLUMED is
PRINT W STRIDE 10
TARGETED TYPE CMAP INDEX CMAPINDEX MAP CMAPVALUES
STEER CV 1 TO 3.0 VEL 0.5 KAPPA 500.0
ENDMETA

where the CMAPVALUES file contains only one map (for details on the format of the CMAPINDEX and
CMPAVALUES file see section 4.18.3).

Additionally one can specify another keyword SQRT that transform the
metrics in its square root, namely the mean square deviation would be
91

changed into the more commonly used root mean square deviation. Similarly it applies to DMSD and CMAP. It should be stressed that the use of this
keyword is particularly dangerous especially for values close to zero since the
square root has a cusp there. The former example would look like:
Example.

PRINT W STRIDE 10
TARGETED TYPE MSD SQRT FRAMESET ref frame.pdb
STEER CV 1 TO 3.0 VEL 0.5 KAPPA 500.0
ENDMETA

4.19

Contact Map

The Contact Map is defined as the sum of the contacts between a number
of atom pairs specified by the user. A contact is defined in terms of the
switching functions introduced for the coordination number CV in section
4.6. Each contact can have its own set of parameters, which are defined in
a file specified by the keyword INDEX. See the documentation in paragraph
4.18.3 for a detailed explanation of the index file syntax.
As for PCV in contact map space (see paragraph 4.18.3), one can define
contacts between the centers of mass of groups of atoms. The name of the
file in which groups are defined is specified by using the GROUP keyword. See
the documentation in paragraph 4.18.3 for additional details.
The optional flag NOPBC can be used to calculate distances without applying periodic boundary conditions.
Example.
The following command defines the CV Contact Map. Distances are calculated without periodic boundary
conditions.
CMAP INDEX fr.ndx SIGMA 1.

NOPBC

The fr.ndx file should contain the details on how each of the switching functions in the contact matrix
should calculated:
CONTACT
CONTACT
CONTACT
CONTACT

1
2
3
4

1 2 3.0 6 10 100.0 1
1 15 3.0 6 10 100.0 1
17 2 3.0 6 10 100.0 1
15 19 3.0 6 10 100.0 1

4.20

Energy

The ENERGY keyword instructs PLUMED to use the total potential energy of
the system [46, 47, 48, 49] as a CV. Currently this CV is available only in
GROMACS4, AMBER and DL POLY.
Example.
The following lines instruct PLUMED to use the potential energy of the system as a CV.
ENERGY SIGMA 100.0

4.21

Helix loops

The HELIX keyword instructs PLUMED to use the number of α-helix loops as a
CV. A helix loop is formed when the pair of dihedral angles (Φ, Ψ) for three
consecutive residues along the chain all adopt a particular pair of reference
values (Φ̄, Ψ̄). Typically in an alpha helix Φ̄ and Ψ̄ have values of -1.200 and
-0.785 respectively. Having specified reference values for the dihedrals the
total number of loops is calculated using:
s=

N
−1 i+1
X
Y
i=2

ih
i
1h
cos(Φj − Φ̄i ) + 1 cos(Ψj − Ψ̄i ) + 1 ,
j=i−1 4

(4.3)

where N is the total number of residues. This CV requires the user to specify
the number of loops with the NLOOP keyword. Then for each loop the user
should provide three sets of four atoms that define the dihedrals Φi−1 , Φi , Φi+i
and a reference value for the dihedral Φ̄i along with three sets of four atoms
that define the dihedrals Ψi−1 , Ψi , Ψi+i and a reference value the dihedral Ψ̄i .
Example.
The following lines instruct PLUMED to use the number of α-helix loops as a CV.
HELIX NLOOP 2 SIGMA 0.1
10 12 14 20 20 22 24 35 35 37 39 49 -1.200 12 14 20 22 22 24 35 37 37 39 49 51 -0.785
20 22 24 35 35 37 39 49 49 51 53 59 -1.200 22 24 35 37 37 39 49 51 51 53 59 61 -0.785

4.22

PCA projection

The PCA keyword instructs PLUMED to use as CV the projection of a set
of atoms on a previously calculated Principal Components Analysis (PCA)
eigenvector[50]. A typical application of this CV is to explore the system
along its principal directions of fluctuations [51].
The current conformation X is first aligned to a reference structure X ref
(except if the optional NOALIGN keyword is used) and then projected on the
specified eigenvector e0 . The set of atoms used for alignment must be the
same as the one the eigenvector refers to. The user should specify the filename
of the reference structure after the keyword FRAME and the filename of the
eigenvector after the keyword EIGENVEC.
If the optional DIFF keyword is used, the difference between the current
(aligned) conformation and the centered reference structure is projected on
the eigenvector.
The CV is calculated as :
sP CA (X, X ref , e0 ) =

N
X

< Rref (Xi − X CM ), e0i >

i=1

or, when using the DIFF keyword, as :
sP CA (X, X ref , e0 ) =

N
X

< Rref (Xi − X CM ) − X0ref , e0i >

i=1

where N is the number of atoms used to align and project, Rref is the 3×3
rotation matrix, calculated using the Kearsley [45] algorithm, that optimally
overlaps the current centered set of atoms (Xi − X CM ) into the centered
reference set of atoms X0ref , X CM is the current centroid of conformation X
(i.e. identical masses are assumed) and <, > denotes the usual inner product.
If the optional NOALIGN keyword is used, no alignment is performed (i.e.
the above rotation matrix R is the identity). This latter keyword is supported
mainly for debug purposes.

Example.
The following command defines a CV as projection, along the eigenvector listed in egv0.dat, of the current
conformation aligned to the reference structure specified in ref.dat:
PCA FRAME ref.dat EIGENVEC egv0.dat SIGMA 0.1
where ref.dat contains the reference structure in the format ATOMID X Y Z:
# my reference structure
2 1.324 1.045 1.550
5 1.316 1.174 1.469
6 1.345 1.174 1.350
7 1.287 1.286 1.538
and egv0.dat contains the eigenvector in the same format:
# my first PCA eigenvector
2 0.12 3.28 0.19
5 0.53 1.37 1.10
6 -0.23 1.93 0.33
7 5.32 -2.56 1.44

More eigenvectors can be used with additional PCA keywords as additional
CVs, nevertheless, the set of atoms and the reference structure must be the
same for all of them.
Example.
The following command defines two CVs as projections, along two distinct eigenvectors listed in egv0.dat
and egv1.dat, of the difference between the current aligned conformation and the reference structure
specified in ref.dat:
PCA FRAME ref.dat EIGENVEC egv0.dat DIFF SIGMA 0.1
PCA FRAME ref.dat EIGENVEC egv1.dat DIFF SIGMA 0.1

Some final notes :
• the atom and eigenvector coordinates should be in engine units (e.g. in
nm for GROMACS)
• the alignment is not mass-weighted (identical masses for all the involved
atoms are assumed), and since the routine used for the alignment is the
same as the one used to calculate the RMSD, the same cautions hold. In
particular, (i) users should ensure that the value of the hardcoded limit
(MAXATOMS RMSD in metadyn.h) is sufficiently large for their needs prior
to compilation, (ii) the use of double precision code is recommended as
well as (iii) the use of the directive ALIGN ATOMS. See also 4.28.
95

• starting the simulation from the identical conformation used for alignment should be avoided since it could generate numerical instabilities
in the calculation of the rotation matrix, possibly leading to a crash of
the simulation.
• comments line beginning with ’#’ are allowed but blank lines must be
avoided from the input files.

4.23

SPRINT topological variables

The SPRINT keyword instructs PLUMED to use as CV the “Social PeRmutation
INvarianT” coordinates described in Ref. [52]:
√
Si =

max,sorted
N λmax vi

where λmax and vimax are the largest eigenvalue and the corresponding
eigenvector of the N × N contact matrix among N atoms:
N
X

Cij vjmax = λmax vimax

j=1

and Cij is the same coordination function defined for the CV COORD in Section 4.6, with the corresponding parameters specified by the keywords R 0,
D 0 (both in units of the host MD code, e.g. Bohr for CPMD), NN, and MM.
In Si the eigenvector components are sorted from the smallest to the largest,
so that S1 < S2 < S3 < ... < SN . Which Si has to be used as CV is specified
by the keyword INDEX: e.g. INDEX 2 means that the 2nd one has to be used,
namely S2 . Note that the sorting of the eigenvector is performed only within
sets of atoms of the same element (automatically identified in PLUMED by
the atomic mass): if the list of N atoms includes e.g. two carbon and four
hydrogen atoms, the SPRINT coordinates are sorted only within carbon atoms
and within hydrogen atoms, without mixing the two elements (since atoms of
different elements are not indistinguishable). Due to this, after keyword R 0
it must be specified the number of different element pairs and the R 0-value
for each possible pair of elements. The same holds for D 0. See below for an
example.

Example.
The following command defines a SPRINT CV S2 from a set of two carbon (1,2) and four hydrogen atoms
(3,4,5,6):
SPRINT LIST NN 6 MM 12 R 0 3 5.0 4.2 4.2 INDEX 2 SIGMA 1.0
all->
1 2 3 4 5 6
all<-

In the example above R 0 3 5.0 4.2 4.2 refers to the three pairs of
elements C-C, C-H, and H-H.
Note that it is desirable to keep a tail of the coordination function long
enough so that the system appears formally as a single connected cluster of
atoms. If a part is disconnected from the rest, the Perron-Frobenius theorem does not hold anymore and several Si components will go to zero (see
Ref. [52]).

4.24

Radial distribution function

The keyword RDF instructs plumed to use the number of distances between
atoms in a given range as a collective variable. A recent paper [53] showed
how a number of such collective variables could be used to define the instantaneous radial distribution function and how this description of the structure
of small clusters could be used to enhance sampling in reconnaissance metadynamics simulations.
The number of distances within a given range is calculated using:
s=

XZ b
i,j

w(r − rij )dr

(4.4)

where rij is the distance between atoms i and j and w is a Gaussian window function with width σ. When using multiple such CVs calculating the set
of distances separately for each bead in the RDF would be computational
expensive and counterproductive. As such all RDF collective coordinates
that involve the same atoms are calculated at the same time. In order that
multiple RDF collective coordinates can be used in a single calculation RDF
beads are labeled using the RDF LABEL keyword. Obviously if all your
RDF CVs come from the same RDF the RDF LABEL should be 1 for all
your RDF coordinates. An example input for this type of CV is as follows:

Example.
The following command instructs plumed to calculate two collective coordinates the number of distances
between 3.0 and 3.5 and the number of distances between 3.5 and 4.0. The value of σ in this calculation
is 0.25.
RDF RDF LABEL 1 LIST RANGE 3.0 3.5 WIDTH 0.25
RDF RDF LABEL 1 LIST RANGE 3.5 4.0 WIDTH 0.25
all->
1 2 3 4 ...
all<-

Example.
The following command instructs plumed to calculate two collective coordinates the number of distances
between atoms in group ¡all¿ between 3.0 and 3.5 and the number of distances in group ¡all2¿ between
3.5 and 4.0. The value of σ in this calculation is 0.25.
RDF RDF
RDF RDF
all->
1 2 3 4
all
5 6 7 8
all2-<

4.25

LABEL 1 LIST RANGE 3.0 3.5 WIDTH 0.25
LABEL 2 LIST RANGE 3.5 4.0 WIDTH 0.25
...

...

Angular distribution function

Rather than using the radial distribution function as a collective variable one
can use the distribution of angles between central atoms and the atoms in
the first hydration sphere [53]. Once again angular distribution functions of
this type have been combined employed successfully with the reconnaissance
metadynamics algorithm. The number of angles within a given range is
calculated using:

w =

N X
N X
N
X

σ(rij )σ(rik )

i=1 j=1 k=1

where

σ(r) =

1−

r−d0
r0

1−

r−d0
r0

Z b
a

w(θ − θjik )dθ

(4.5)

n
m

(4.6)

where w is once again a Gaussian window function with width σ. The
two switching functions σ(rij ) and σ(rik ) are used to make sure that one only
takes angles in the first hydration sphere. A detailed description of how to set
the parameters in these functions can be found in the section of this manual
on coordination numbers. Much as was described above for the RDF the
RDF LABEL keyword ensures that the code does not waste time calculating
all the angles involved in these ADF beads multiple times. If all your ADF
CVs are based on the positions of the same atoms then RDF LABEL should
be 1 for all your ADF coordinates. In other words this keyword should be
used to differentiate between the various angular distributions functions you
are calculating in input.
Example.
The following command instructs plumed to use the number of angles in the first coordination sphere that
are between 0.5 and 1.0 radians as a CV. The value of σ in this calculation is 0.25.
ADF RDF LABEL 1 LIST RANGE 0.5 1.0 WIDTH 0.25 R 0 3.0 NN 6 MM 12 D 0
1.5
all->
1 2 3 4 ...
all<-

4.26

Polynomial combination of CVs

Polynomial combinations of collective variables, given by the functional form
k
X

ci (CVi − si )ni

i=1

can be defined using the POLY directive, followed by the TERMS keyword that
specifies the number k of terms to be included in the sum. In the next k
lines, the compulsory keyword CV defines the number CVi of the collective
variable in the i-th term; optionally the values of ci , si and ni can be specified
preceded respectively by the keywords COEFF, SHIFT and EXP. If not specified they assume the default values of ci = 1 si = 0 ni = 1. Fractional and
negative values of the exponent are allowed, but singularities should be carefully avoided by imposing restraints on the regions explored by the collective
variables involved.
99

Warning! To combine more than 20 collective variables (k > 20), the hardset allocation for intpar
and vecpar in metadyn.h should be increased accrdingly and the code recompiled.

Example.
The following command defines a new CV as CV1 − CV22 :
POLY TERMS 2 SIGMA 1.0
CV 1
CV 2 COEFF -1 EXP 2

More complicated functions might need multiple auxiliary CVs to be defined.
Example.
√
The following command defines a new CV as 3CV1 + CV2
POLY TERMS 2
CV 1 COEFF 3
CV 2
NOHILLS CV 3
POLY TERMS 1 SIGMA 1.0
CV 3 EXP 0.5

4.27

Function of CVs

General function of CVs with arbitrary form can be defined using the FUNCTION
directive. You should define this function AFTER you defined the variables
you are combining.
This function, still experimental, makes use of libmatheval library that
you should download from http://www.gnu.org/software/libmatheval/
and compile along with them. You may find useful infos on the syntax in
http://www.gnu.org/software/libmatheval/manual/. It is still experimental so if you want to compile it you should hack your Makefile after
patching. I tried in NAMD2.7. You should change a couple of lines in the
Makefile.
Example.
in the Makefile (dots represent the arguments that you already find in the Makefile)
CXXBASEFLAGS = ... $(COPTI)/my/path/to/matheval/install dir x86 64/include
and later
namd2: -L/my/path/to/matheval/install dir x86 64/include/lib -lmatheval ...

100

while in SANDER (v10) I had to change
Example.
in $AMBERHOME/src/config amber.h I added the needed stuff:
AMBERBUILDFLAGS=-DAMBER -DHAVE MATHEVAL -L/my/path/to/matheval/install dir x86 64/lib
-lmatheval -I/my/path/to/matheval/install dir x86 64/include

It is crucial to put the function of CVs after all the cvs you use to define it!!!!

Example.
The following command defines a new CV as function of the two angles
PRINT W STRIDE 1
ANGLE LIST 7 9 15
ANGLE LIST 9 15 17
FUNCTION " 2*CV 1 + (CV 2+CV 1)0.2 "

Additional note is that the function is not periodic.

4.28

A note on periodic boundary conditions

PLUMED is designed so that for the majority of the CVs implemented the
periodic boundary conditions are treated in the same manner as they would
be treated in the host code. However, there are some exceptions; namely:
• Average coordinate of a group of atoms;
• RGYR;
• DISTANCE, MINDIST, COORD, HBONDS, and path collective variables S PATH
and Z PATH with contact map metrics (TYPE CMAP), if the NOPBC flag is
used;
• Path Collective Variables S PATH and Z PATH with RMSD (TYPE RMSD);

101

• ALPHARMSD, ANTIBETARMSD, PARABETARMSD.
In all these cases, it is essential that the atoms involved in the definition
of the CV are all part of a single, unbreakable object such as a molecule.
Furthermore, it is essential that in the coordinates passed to PLUMED the
molecules are kept intact. We are aware of at least two cases where this
condition is not satisfied; namely, when using domain decomposition or the
option for periodic molecules within the host code GROMACS4.
In these cases one must use the additional directive ALIGN ATOMS, which
takes the LIST keyword. By using this command the user can define an ordered group of atoms, which are to be aligned such that the distance between
adjacent atoms in the list is minimized. In the majority of cases, the atoms
in this list should be in the same order as they appear in the pdb file as
usually these files are arranged in way that reflects how close together atoms
are in the molecular topology. As for the number of atoms that must be
specified in this list it is often sufficient to just specify those atoms involved
in the CVs. However, in cases where the atoms involved are separated by
a large distances along a chain alignment of intermediate atoms will also be
required.
Example.
Input for running a metadynamics simulation using the end-to-end distance in a protein as a CV with
GROMACS4 and domain decomposition:
PRINT W STRIDE 10
HILLS HEIGHT 2.0 W STRIDE 10
DISTANCE LIST 9 238 SIGMA 0.35 NOPBC
ALIGN ATOMS LIST
C-alpha->
9 16 31 55 69 90 102 114 124 138 160 174 194 208 224 238
C-alpha
-file HILLLS
-out fes.dat
-hills nhills
-aver naver

number of collective variables NCV
CVs for the free-energy surface
mesh dimension. DEFAULT :: 100
size of the mesh of the output free energy
optional definition of the FES domain
how often the FES is written
the hills are cut off at 1.e-6
the hills are cut off at 6.25 std dev from the center
[0; 2π] periodicity on the x CV,
if -fix is not used 2pi is used
[−π; π] periodicity on the x CV,
if -fix is not used 2pi is used
kT in the energy units
apply periodicity using degrees
writing output the bias for a well tempered mtd
input file
output file
number of Gaussians that are read
time-average the bias profile over the last naver Gaussians

The program works in the following way.
Using the -file and -out flags one tells the program the names of the
input file containing the Gaussian hills file and the name of the file in which
the free-energy will be outputted. In the absence of any instruction sum hills
assumes these files are to be called HILLS and fes.dat.
The number of CVs in the HILLS file is specified using -ndim whilst the
number of CVs the output free energy surface is to be plotted as a function
of is controlled by -ndw followed by the list of CVs in the desired order. Note
that after being read the CV are reordered in the code according to ndw. As
a direct result EVERY output uses this new order. If the number of CVs
requested using -ndw is less than the number of CVs in in HILLS file the CVs
not specified for output will be integrated out with the Boltzmann weight
Kb T specified by -kt (kT must be given in the energy units that were used
by in the code in which the simulation was performed)
The position, width and height of the Gaussians are read from the file
specified by the -file option (HILLS is the default) and the free-energy
surface is printed out on a grid in gnuplot format with a blank line added
after each block of data. Gnuplot is very handy for quick visualization of 3D
data (e.g. the FES as a function of 2 CVs).
For efficiency, the Gaussians are truncated at a certain distance from their
center before being placed on the grid. -cutoff s and -cutoff e allow one
104

to tune this cutoff distance. With -cutoff s the value read specifies this
truncation distance in terms of the number of standard deviations from the
center. For consistency this should be set equal to the DP2CUTOFF used
in the metadynamics simulation (the default value for this parameter is 6.25
both in PLUMED and in sum hills). Alternatively, the user may specify this
cutoff distance as the distance at which the energy of the Gaussian hill falls to
less than some critical value specified using the -cutoff e flag. The energy
specified after this flag must be in the units used during the metadynamics
simulation. N.B. if the cutoff used in post-processesing is different to that
used during the simulation the calculated free energy differs from the bias
which was actually applied during the metadynamics simulation.
The size of the output grid can be controlled either with -ngrid followed
by the number of grid points in each CV direction or by -dp followed by
the size of the voxel in each CV direction. If neither of these options are
specified, the grid size is assumed to be 100 in each direction and the voxel
size is calculated such that all the input Gaussians fit into the grid. The
-fix option allow one to fix the boundaries of the output free energy and
after this flag. two real numbers are expected for each CV.
-stride allows one to print out the evolution of the free-energy as a
function of time, i.e. as a function of the index of the Gaussian. -stride
expects an integer which specifies how many Gaussians are added to the
calculated free energy surface between print outs. The progressive free energy
is printed in files named fes.dat.XXX where XXX is an increasing counter.
The final free-energy is printed in fes.dat. This is useful for creating movies
of the how the free-energy wells fill as a function of time.
-hills is used to integrate only a part of the Gaussian files. It expects
an integer that specifies the maximum number of Gaussians to be read from
input.
-naver is used to plot the time average of the bias profile over the last
given number of hills (see Eq. 33 in Ref. [23])
When periodic CVs like angles and dihedrals are used as CVs, periodicity
options must be specified. The flags -pi and -2pi, which are followed by
the CV index, specify that that the CV is periodic between [-pi;pi] or [0;2pi]
respectively. The -grad flag specifies that the CV values are being given
degrees rather than radians.

105

5.2

Evaluating collective variables on MD trajectories

The program driver can be used to evaluate the value of any of the CVs
implemented in PLUMED for all the structures in a given trajectory file.
For GROMACS users: note that driver has some limitation in the choice
of the simulation box shape. With GROMACS, it could be convenient to use
the -rerun option of mdrun together with PLUMED.

5.2.1

Installation instructions

Before compiling driver, the files contained in the common files directory
must be linked to the utilities/driver directory (i.e. this can be done
using ./getlinks.sh).
To compile using g95/gcc, type:
make arch=g95
To compile with gfortran/gcc, type:
make arch=gfortran
To compile with gfortran/gcc on a 64 bit machine, type:
make arch=gfortan 64
To compile with ifort/icc Intel compilers, type:
make arch=intel
To use other compilers and/or compiler flags you will need to modify the
Makefile.

5.2.2

Usage

This program takes its input parameters from the command line. If run
without options, this brief summary of options is printed out.

106

Example.
Invoking driver without arguments prints a list of the available options:
USAGE :
driver -pdb PDB FILE -dcd DCD FILE -plumed PLUMED FILE -ncv
(-interval min1 max1 min2 max2 -out OUT FILE -nopbc -cell CELLX CELLY CELLZ)
-pdb
-dcd
-plumed
-ncv
-out
-interval
-nopbc
-cell

pdb (one frame) connected to dcd file
trajectory file
PLUMED-like input file
number of collective variables
pdb output filename for clustering format (optional)
extract frames with CV in this interval (optional)
don’t apply pbc (optional)
provide fixed box dimension in Angstrom
for orthorhombic PBC (optional)

The user must provide a structure file in PDB format (only one frame,
not a movie), within which the atoms’ masses and charges are inserted in
the occupancy field and in the B-factor field respectively. This data must be
provided as it is required for the calculation of certain collective variables,
such as the center of masses or dipoles.
The trajectory file must be a CHARMM format DCD file (also NAMD
2.1 and later). To convert other trajectory files into this format, the user can
download the utility CatDCD from:
http://www.ks.uiuc.edu/Development/MDTools/catdcd/.
Within the input file driver the user can specify a printing stride (in
number of DCD frames) using the keyword W STRIDE and the list CVs he/she
wishes to calculate. This file has the same syntax as the PLUMED input file.
The CVs value will be printed on the COLVAR file. For a complete list of the
CVs driver can calculate, please see section 4.
Example.
If you wish to evaluate the distance between atom 13 and 17 for every frame in your trajectory file, the
input file for driver would be as follows:
PRINT W STRIDE 1
DISTANCE LIST 13 17
ENDMETA

By default, driver uses orthorhombic periodic boundary conditions and
looks for cell information in the DCD file. If these are not present, you must
107

either provide the fixed dimensions of the box in Angstrom (only orthorhombic PBCs are allowed) using the keyword -cell CELLX CELLY CELLZ or
switch off the pbc using the -nopbc flag.
Driver can also be used to extract frames in a specific window of the CVs
space and write these extracted frames to a PDB file. To use this functionality use the option -out output.pdb and specify the interval with -interval
CVmin CVmax.
Example.
To extract those frames in which the distance between two specified atoms is between 10.0 and 12.0 Å
while the angle defined by three specified atoms is between 2.0 and 2.3 rad, the user should prepare the
following input file (named plumed.dat):
PRINT W STRIDE 1
DISTANCE LIST 13 17
ANGLE LIST 19 20 22
ENDMETA
and then invoke driver using the following command:
./driver -pdb dialanine.pdb -dcd dialanine.dcd
window.pdb -interval 10.0 12.0 2.0 3.0

-plumed plumed.dat -ncv 2 -out

The output file window.pdb is written in a format appropriate for use with the GROMACS tool g cluster,
which can be used to perform a cluster analysis on your set of frames.

5.3

Processing COLVAR files

Since PLUMED 1.2, a more flexible format has been adopted for COLVAR files.
These files can be parsed with the simple plumedat.sh tool. If run it without
options, this brief summary of options is printed out.

108

Example.
Invoking plumedat.sh without arguments prints a list of the available options:
syntax:
plumedat.sh f1 f2 ... < file
or
plumedat.sh -l < file
file is the name of the COLVAR file
f1, f2, ... are the names of the required fields
if a required field is not available in the COLVAR file, "NA" is written in the output
with -l, the available choices are listed
example:
plumedat.sh time temp < COLVAR
prints a two-column file, with the time in the first column and the temperature in the
second column

5.4

PLUMED as a standalone program

PLUMED can be run as a standalone program so that it can be easily integrated
in a script. This is particularly useful whenever the time between one PLUMED
call and the other is large (e.g. in ab initio programs) and one aims at minimizing the time needed for implementation. plumed standalone consists in
a simple program, very much like the driver tool, that reads the configuration and prints out the forces to be added and the other various output files
typical of PLUMED.

5.4.1

Installation instructions

Before compiling plumed standalone, the files contained in the common files
directory must be linked to the utilities/standalone directory (i.e. this
can be done using ./getlinks.sh).
To compile using g95/gcc, type:
make arch=g95
To compile with gfortran/gcc, type:
make arch=gfortran
To compile with gfortran/gcc on a 64 bit machine, type:
make arch=gfortan 64
To compile with ifort/icc Intel compilers, type:
109

make arch=intel
To use other compilers and/or compiler flags you will need to modify the
Makefile. In this way you may obtain an executable which is plumed standalone.

5.4.2

Usage

PLUMED standalone expects a standard PLUMED input, a configuration file for
your system (much in the format of xyz) and gives a screen output. The
configuration file has the following format:
Example.
Example of plumed standalone input configuration
TIME 0.100000 100
AMPLI 1.000000
BOLTZ 0.001987
BOX 1000.000000 1000.000000 1000.000000
-13.523670417
0.140088464
4.263822219
-9.821878867
-0.695171589 5.042585958
-7.193416272
-2.392065561 2.830069705
-3.844793108
-3.737129154 4.056782255
.
.
.

131.207999
101.119998
163.190996
129.197997

0.000000
0.000000
0.000000
0.000000

The configuration file must have the keyword TIME followed by the step
length (in the units of the program) and AMPLI which is the conversion factor
between Å and the actual length unit in the configuration file. This is needed
when a RMSD calculation with respect to a pdb file in Å must be performed.
In this case, if the configuration is in Bohr, the AMPLI factor must be set
equal to 1.889725989. BOLTZ is the Boltzmann factor in the chosen units and
BOX followed by three numbers specifies the dimensions of the orthorombic
simulation box.
The following lines specify the configuration. The first three columns are
the x, y and z coordinates, the fourth is the mass in units of the program
(in the case above a coarse grained program takes the mass of all the single
residues) and in the last column you may put the charge (in the code above
it was not needed). PLUMED invocation may be done with command line:

110

Example.
Example of invocation of PLUMED standalone
./plumed standalone -coord x.xyz -plumed plumed input.cfg

plumed standalone returns back a file named plumed forces.dat which
contains the additional energy and forces from the bias potential calculated
by the PLUMED module. Eventually, these must be summed to the ones of the
original program.
Example.
Example of plumed standalone output file. In the first line is the additional bias potential calculated by
PLUMED, followed by the x, y and z components of the force for each atom.

2.35680
3.335998704981455 -3.684076849212145 -16.098665997131253
-0.9064447402311846 2.913567594563722 -1.0502785066801152
-0.7948355287665761 0.9861425018822194 0.6558434539267759
-1.9437110159112139 -0.6087819098852943 2.4119663558915296
.
.
.

5.5

Reweighting well-tempered metadynamics calculations

From a converged metadynamics run we can calculate directly the canonical
probability distribution of the collective variables at a given temperature.
On the contrary, the statistics of other degrees of freedom is somehow distorted by the application, during the simulation, of a time-dependent external potential on the CVs. Different possible techniques have been proposed
to reconstruct the probability distribution of variables other than the CVs
[54, 29, 21].
In well-tempered metadynamics, the reconstruction of the distribution of
variables different from the CVs is particularly simple since for long times
the amount of bias added decreases to zero and the system becomes closer
and closer to equilibrium. The algorithm described in Ref. [21] consists of
three different steps:
111

1. Accumulate the histogram of the CVs plus the other variables of interest between two updates of the bias potential;
2. When a new Gaussian is added, evolve the histogram following:
P (R, t + ∆t) = e−β(V̇ (S(R),t)−hV̇ (S,t)i)∆t P (R, t),

(5.1)

where P (R, t) is the biased probability distribution, V̇ (S(R), t) the
time derivative of the bias potential and the average in the exponent is
calculated in the biased ensemble;
3. At the end of the simulation, the unbiased distribution PB (R) can be
recovered from the histogram collected by using a standard umbrella
sampling reweighting:
PB (R) ∝ eβV (S(R),t) · P (R, t).

5.5.1

(5.2)

Installation instructions

As reweight.f90 is a simple fortran 90 program, the installation is straightforward so long as you have a fortran compiler available on your machine.
As an example, with the gnu g95 compiler one would compile reweight.f90
using the following command:
g95 -O3 reweight.f90 -o reweight

112

5.5.2

Usage

Example.
We performed a metadynamics run with 2 CVs and we are interested in reconstructing the distribution of
a third variable.
reweight -colvar COLVAR -hills HILLS -ncv 2 -nvar 3 -stride 1 -fes 3 -temp 300
-welltemp
-hills
-colvar
-out
-ncv
-nvar
-stride
-fes
-temp
-ngrid
-nreject
-timeout
-pi
-welltemp
-kjoule

HILLS filename
COLVAR filename
FES filename
number of variables in HILLS
number of variables in COLVAR
ratio between COLVAR and HILLS stride
ID of the variables for FES in output
temperature in Kelvin
histogram grid dimension
discard initial steps
stride for FES printout
ID of the variables with [−π; π] periodicity
control for well-tempered metadynamics
energy in kjoule/mol

The code needs two files with the same format of the PLUMED HILLS and
COLVAR files. In the latter, the metadynamics CVs should appear in the first
d column followed by the variables whose distribution one wants to reweight.
The ratio between the stride in COLVAR and HILLS must be constant and
greater than 1. The more data you have for the histogram, the better.
Some important things to keep in mind:
• For the choice of the bin size, please follow the suggestions described
in section 3.4.5;
• Eq. 5.1 is exact. However, at the beginning of the simulation the
average of V̇ (S(R), t) can be calculated only approximately. Luckily,
a possible initial error is recovered for long times. Alternatively, one
could discard the first part of the trajectory using the -nreject option.
Please always check that your results are robust to a discard of initial
parts of the trajectory;
• As for the calculation of the FES with sum hills, remember to control
the convergence by plotting the reconstructed distribution at different
times by using -timeout;
113

5.6

Bias-exchange simulations via the linux
shell

Bias-exchange simulations [18] can be performed using every MD engine,
employing the tool in the directory utilities/bias-exchange. Each replica
(walker) runs independently in a different directory for a duration τexch equal
to the time between two exchanges of the bias. When all runs are finished,
a script called bias-exchange.sh takes care of attempting exchanges of the
bias and re-launching the simulation for each replica.
This procedure via the linux shell can be inefficient in terms of computer
time (due to the need to stop and restart the simulations to perform exchanges) if τexch is chosen very short, e.g. < 100 timesteps. However in the
literature times larger than 10 ps are typically employed for classical simulations of biomolecules, while for ab-initio simulations this should not be an
issue already for τexch ∼ 0.1 ps. Note that for GROMACS PLUMED provides
also bias-exchange within a single parallel run, which is computationally more
efficient (see Section 3.6.2).
A simulation proceeds as follows (see also the directory example). First,
the Fortran code exchange-tool.f90 is compiled using the command make
(modify the Makefile if necessary). The script bias-exchange.sh and the
program exchange-tool.x must be copied in a chosen directory BASE DIR,
and the following variables inside the script bias-exchange.sh must be modified:
• BASE DIR
• NWALKER is the number of replicas
• KBT is the Boltzmann factor kB T in units of the MD engine
• NSIMULATIONS is the maximum number of runs between exchanges for
each walker
Then, inside BASE DIR the user must create NWALKER directories called walker1,
walker2, etc., and in each one the following files must be present:
• an input for restarting an MD simulation of length τexch
• a PLUMED input file called plumed.dat (see example below)

114

• a script run-walker.sh (see example below)
Finally, the bias-exchange simulation is started by launching the MD runs in
each walker1, walker2, etc. directories using the command ./run-walker.sh
& (the latter script must be executable: it can be made so e.g. using command chmod a+x run-walker.sh). This script must launch the MD run in
foreground (no & character at the end of the line) and it must end with the
following lines:
touch READY
../bias-exchange.sh &
After a time τexch the runs finish and control passes to the bias-exchange.sh
script (in background), which attempts to exchange the bias (i.e. files HILLS
and plumed.dat, not atomic coordinates) between pairs of randomly chosen
replicas. Detailed output is printed in the file bias-exchange.log.
It is required to employ the keyword HILLS LABEL in each plumed.dat
input file, with a different label for each walker: in this way PLUMED includes
the information (necessary for the bias-exchange tool) about the active CVs
in the COLVAR and HILLS files.

115

Example.
Example of plumed.dat input file in directory walker1. Note the keyword HILLS LABEL and that only CV
1 is biased.
HILLS RESTART HEIGHT 0.001 W STRIDE 100
PRINT W STRIDE 10
HILLS LABEL A
COORD LIST 1 NN 6 MM 12 R 0 5.0 SIGMA 0.1
g234->
2 3 4
g234 NN 6 MM 12 R 0 5.0 SIGMA 0.1
g134->
1 3 4
g134 NN 6 MM 12 R 0 5.0 SIGMA 0.1
g124->
1 2 4
g124 NN 6 MM 12 R 0 5.0 SIGMA 0.1
g123->
1 2 3
g123> output
# this creates and empty file called READY,
# telling the bias-exchange tool that the MD run has finished:
touch READY
# note that bias-exchange.sh is launched in background:
../bias-exchange.sh &

116

At the end of the bias-exchange simulation, the multidimensional freeenergy landscape as a function of several CVs can be reconstructed according to the weighted-histogram algorithm in Ref. [29], e.g. employing the
METAGUI program [30] downloadable from
http://www.plumed-code.org/contributions
The bias-exchange simulation can be continued, after NSIMULATIONS steps
are reached and all replicas have finished their run, in the following way:
• enlarge NSIMULATIONS in the script bias-exchange.sh
• in each walker-directory, delete the file READY
• in each walker-directory, restart the run using ./run-walker.sh &

117

Bibliography
[1] M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci, R. A. Broglia, M. Parrinello,
PLUMED: a portable plugin for free energy calculations with molecular
dynamics, Comp. Phys. Comm. 180 (2009) 1961.
[2] B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation, J. Chem. Theory Comput. 4 (3) (2008) 435–447.
[3] J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa,
C. Chipot, R. D. Skeel, L. Kalé, K. Schulten, Scalable molecular dynamics with NAMD, J. Comput. Chem. 26 (16) (2005) 1781–802.
[4] W. Smith, C. Yong, P. Rodger, Mol. Simulat. 28 (2002) 385–471.
[5] S. J. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comp. Phys. 117 (1995) 1–19.
[6] D. A. Case, T. A. Darden, T. E. C. III, C. L. Simmerling, J. Wang,
R. E. Duke, R. Luo, K. M. Merz, D. A. Pearlman, M. Crowley, R. C.
Walker, W. Zhang, B. Wang, S. Hayik, A. Roitberg, G. Seabra, K. F.
Wong, F. Paesani, X. Wu, S. Brozell, V. Tsui, H. Gohlke, L. Yang,
C. Tan, J. Mongan, V. Hornak, G. Cui, P. Beroza, D. H. Mathews,
C. Schafmeister, W. S. Ross, P. Kollman, AMBER 9, University of California, San Francisco.
[7] M. Harvey, G. Giupponi, G. D. Fabritiis, ACEMD: Accelerated molecular dynamics simulations in the microseconds timescale, J. Chem. Theory and Comput. 5 (2009) 1632.

118

[8] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso,
S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann,
C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,
F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto,
C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov,
P. Umari, R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular
and open-source software project for quantum simulations of materials,
J. Phys.: Condens. Matter 21 (2009) 395502.
[9] A. Laio, M. Parrinello, Escaping free energy minima, Proc. Natl. Acad.
Sci. USA 99 (2002) 12562–12566.
[10] G. Torrie, J. Valleau, Nonphysical sampling distributions in monte carlo
free energy estimation: Umbrella sampling, J. Comput. Phys. 23 (1977)
187–199.
[11] S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, P. A. Kollman,
Multidimensional free-energy calculations using the weighted histogram
analysis method, J. Comput. Chem. 16 (1995) 1339–1350.
[12] B. Roux, The calculation of the potential of mean force using computersimulations, Comput. Phys. Comm. 91 (1995) 275–282.
[13] C. Jarzynski, Nonequilibrium equality for free energy differences, Phys.
Rev. Lett. 78 (1997) 2690–2693.
[14] G. Crooks, Nonequilibrium measurements of free energy differences for
microscopically reversible markovian systems, J. Stat. Phys. 90 (5-6)
(1998) 1481–1487.
[15] A. Barducci, G. Bussi, M. Parrinello, Well-tempered metadynamics: A
smoothly converging and tunable free-energy method, Phys. Rev. Lett.
100 (2) (2008) 020603.
[16] P. Raiteri, A. Laio, F. Gervasio, C. Micheletti, M. Parrinello, Efficient
reconstruction of complex free energy landscapes by multiple walkers
metadynamics, J. Phys. Chem. B 110 (2006) 3533–3539.

119

[17] G. Bussi, F. L. Gervasio, A. Laio, M. Parrinello, Free-energy landscape
for beta hairpin folding from combined parallel tempering and metadynamics, J. Am. Chem. Soc. 128 (41) (2006) 13435–41.
[18] S. Piana, A. Laio, A bias-exchange approach to protein folding, J. Phys.
Chem. B 111 (17) (2007) 4553–9.
[19] G. A. Tribello, M. Ceriotti, M. Parrinello, A self-learning algorithm for
biased molecular dynamics, Proc. Natl. Acad. Sci. USA 107 (41) (2010)
17509–17514.
[20] M. Marchi, P. Ballone, Adiabatic bias molecular dynamics: A method
to navigate the conformational space of complex molecular systems, J.
Chem. Phys. 110 (8) (1999) 3697–3702.
[21] M. Bonomi, A. Barducci, M. Parrinello, Reconstructing the equilibrium
Boltzmann distribution from well-tempered metadynamics, J. Comput.
Chem. 30 (11) (2009) 1615–1621.
[22] J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, Quickstep:
Fast and accurate density functional calculations using a mixed gaussian
and plane waves . . . , Comp. Phys. Comm. 167 (2005) 103–128.
[23] A. Laio, F. L. Gervasio, Metadynamics: a method to simulate rare events
and reconstruct the free energy in biophysics, chemistry and material
science, Rep. Prog. Phys. 71 (2008) 126601.
[24] A. Barducci, M. Bonomi, M. Parrinello, Metadynamics, WIREs Comput. Mol. Sci. 1 (2011) 826.
[25] V. Babin, C. Roland, C. Sagui, Adaptively biased molecular dynamics
for free energy calculations, J. Chem. Phys. 128 (2008) 134101.
[26] F. Baftizadeh, P. Cossio, F. Pietrucci, A. Laio, Protein folding and
ligand-enzyme binding from bias-exchange metadynamics simulations,
Current Physical Chemistry (in press).
[27] Y. Crespo, F. Marinelli, F. Pietrucci, A. Laio, Metadynamics convergence law in a multidimensional system, Phys. Rev. E 81 (2010)
055701(R).

120

[28] C. Camilloni, D. Provasi, G. Tiana, R. A. Broglia, Exploring the protein
G helix free-energy surface by solute tempering metadynamics, Proteins
71 (2008) 1647.
[29] F. Marinelli, F. Pietrucci, A. Laio, S. Piana, A kinetic model of trpcage folding from multiple biased molecular dynamics simulations, PLoS
Comput. Biol. 5(8) (2009) e100045.
[30] X. Biarnes, F. Pietrucci, F. Marinelli, A. Laio, Metagui. a vmd interface for analyzing metadynamics and molecular dynamics simulations,
Comp. Phys. Comm. 183 (2011) 203.
[31] A.
Grossfield,
Wham.http://membrane.urmc.rochester.edu/
Software/WHAM/WHAM.html.
[32] D. Provasi, M. Filizola, Putative active states of a prototypic g-proteincoupled receptor from biased molecular dynamics, Biophys. J. 98 (2010)
2347–2355.
[33] C. Camilloni, R. A. Broglia, G. Tiana, Hierarchy of folding and unfolding
events of protein g, ci2, and acbp from explicit-solvent simulations, J.
Chem. Phys. 134 (2011) 045105.
[34] J. B. Abrams, M. E. Tuckerman, Efficient and direct generation of multidimensional free energy surfaces via adiabatic dynamics without coordinate transformations, J. Phys. Chem. B 112 (2008) 15742–15757.
[35] L. Rosso, M. E. Tuckerman, An adiabatic molecular dynamics method
for the calculation of free energy profiles, Mol. Sim. 28 (2002) 91–112.
[36] L. Rosso, P. Minary, Z. Zhu, M. E. Tuckerman, On the use of adiabatic
molecular dynamics to calculate free energy profiles, J. Chem. Phys. 116
(2002) 4389–4402.
[37] L. Maragliano, E. Vanden-Eijnden, A temperature accelerated method
for sampling free energy and determining reaction pathways in rare
events simulations, Chem. Phys. Lett. 426 (2006) 168–175.
[38] Y. Liu, M. E. Tuckerman, Generalized gaussian moment thermostatting:
A new continuous dynamical approach to the canonical ensemble, J.
Chem. Phys. 112 (2000) 1685.
121

[39] J. Vymetal, J. Vondrasek, Gyration- and inertia-tensor-based collective
coordinates for metadynamics. application on the conformational behavior of polyalanine peptides and trp-cage folding, J. Phys. Chem. A
115 (2011) 11455–11465.
[40] D. N. Theodorou, U. W. Suter, Shape of unperturbed linear polymers:
polypropylene, Macromolecules 18 (6) (1985) 1206–1214.
[41] F. Pietrucci, A. Laio, A collective variable for the efficient exploration
of protein beta-structures with metadynamics: application to sh3 and
gb1, J. Chem. Theory Comput. 5 (9) (2009) 2197–2201.
[42] M. Sega, E. Autieri, F. Pederiva, On the calculation of puckering free
energy surfaces, J. Chem. Phys. 130 (22) (2009) 225102.
[43] D. Branduardi, F. L. Gervasio, M. Parrinello, From A to B in free energy
space, J. Chem. Phys. 126 (5) (2007) 054103.
[44] M. Bonomi, D. Branduardi, F. Gervasio, M. Parrinello, The unfolded
ensemble and folding mechanism of the C-terminal GB1 β hairpin, J.
Am. Chem. Soc. 130 (42) (2008) 13938–13944.
[45] S. K. Kearsley, On the orthogonal transformation used for structural
comparison, Acta Cryst. A 45 (1989) 208–210.
[46] H. Li, D. Min, Y. Liu, W. Yang, Essential energy space random walk via
energy space metadynamics method to accelerate molecular dynamics
simulations, J. Chem. Phys. 127 (9) (2007) 094101.
[47] C. Michel, A. Laio, A. Milet, Tracing the entropy along a reactive pathway: The energy as a generalized reaction coordinate, J. Chem. Theory
Comput. 5 (9) (2009) 2193–2196.
[48] D. Donadio, P. Raiteri, M. Parrinello, Topological defects and bulk melting of hexagonal ice, J. Phys. Chem. B 109 (12) (2005) 5421–5424.
[49] M. Bonomi, M. Parrinello, Enhanced sampling in the well-tempered
ensemble, Phys. Rev. Lett. 104 (2010) 190601.
[50] A. Amadei, A. B. M. Linssen, H. J. C. Berendsen, Essential dynamics
of proteins., Proteins Struct. Funct. Genet. 17 (1993) 412–425.
122

[51] L. Sutto, M. D. Abramo, F. L. Gervasio, Comparing the efficiency of biased and unbiased molecular dynamics in reconstructing the free energy
landscape of met-enkephalin, J. Chem. Theory Comput. 6 (12) (2010)
3640–3646.
[52] F. Pietrucci, W. Andreoni, Graph theory meets ab initio molecular dynamics: atomic structures and transformations at the nanoscale, Phys.
Rev. Lett. 107 (2011) 085504.
[53] G. A. Tribello, J. Cuny, H. Eshet, M. Parrinello, Exploring the free
energy surfaces of clusters using reconnaissance metadynamics, J. Chem.
Phys 135 (11) (2011) 114109.
[54] G. Tiana, Estimation of microscopic averages from metadynamics, Eur.
Phys. J. B 63 (2) (2008) 235–238.

123

Index
Absolute position, see Collective variables, Absolute position
ADF, see Collective variables, ADF
Adiabatic bias molecular dynamics,
48
Alpha-beta similarity, see Collective
variables, Alpha-beta similarity
Alpharmsd, see Collective variables,
Alpharmsd
Angle, see Collective variables, Angle
Antibetarmsd, see Collective variables,
Antibetarmsd
Atom lists, 64

Dipole, 76
Distance from an axis, 68
Distance, 67
Electrostatic potential, 81
FUNCTION, 100
Gyration tensor based CVs, 75
Hydrogen bonds, 71
Interfacial water, 73
MSD, 91
Minimal distance, 68
POLY, 99
Parabetarmsd, 80
Path variables, 83
Projection on an axis, 68
Puckering coordinates, 82
Bias exchange simulations, 42
RDF, 97
Bias-exchange simulations via the linux
Radius of Gyration, 74
shell, 114
SPRINT, 96
Targeted, 91
Collective variables
Torsion, 70
ADF, 98
Contact Matrix distance, see CollecAbsolute position, 66
tive variables, Contact Matrix
Alpha-beta similarity, 77
distance
Alpharmsd, 78
Coordination number, see Collective
Angle, 69
variables, Coordination numAntibetarmsd, 79
ber
Contact Matrix distance, 91
Coordination number, 70
DMSD, 91
Difference of distances, 67
Dihedral correlation, 77

Difference of distances, see Collective
variables, Difference of distances
Dihedral correlation, see Collective variables, Dihedral correlation
124

Dipole, see Collective variables, Dipole
Distance, see Collective variables, Distance
Distance from an axis, see Collective
variables, Distance from an axis
DMSD, see Collective variables, DMSD
Electrostatic potential, see Collective
variables, Electrostatic potential
FUNCTION, see Collective variables,
FUNCTION
Grid for metadynamics, 30
Gyration tensor based CVs, see Collective variables, Gyration tensor based CVs
Hydrogen bonds, see Collective variables, Hydrogen bonds
Interfacial water, see Collective variables, Interfacial water
Keywords
R 0 , 97
ABMD, 48
ACYLINDRICITY, 76
ALIGN ATOMS, 95, 102
ALPHABETA, 36, 77
ALPHARMSD, 78, 102
ANGLE, 69
ANTIBETARMSD, 79, 102
ASPHERICITY, 76
BASINS, 54, 56
BASIN TOL, 56, 57
BETA, 69
BIASXMD, 42
CLUSTER, 54, 56
125

CMAP, 92
COEFF, 99
COLVAR, 53, 54, 62, 63
COMMITMENT, 53
COORD, 36, 96, 101
COS, 69, 70
CV LIST, 59
CV, 61, 99
DAFED CONTROL, 62
DAFED STATE, 62
DAFED, 61
DIFF, 94
DIHCOR, 77
DIPOLE, 76
DIR, 67
DISTANCE, 67, 101
DMSD, 84, 92
DRMS, 84
D 0, 96
EIGENVEC, 94
ELSTPOT, 81
ENDMETA, 26
ENERGY, 9, 93
EXPAND PARAM, 57
EXP, 99
EXTERNAL, 50
E sj, 63
FILENAME, 32, 33
FRAME, 94
FUNCTION, 100
GRID, 31, 32
GTPC 1, 75
GTPC 2, 75
GTPC 3, 75
HBONDS, 71, 101
HEIGHT, 58
HELIX, 93
HILLS LABEL, 43, 115

HILLS, 27, 28, 58
INDEX 2, 96
INDEX, 92, 96
INITIAL SIZE, 57
INTERVAL, 36, 37
INVERT, 36, 38, 39
JACOBIAN FORCE, 62
KAPPA2, 76
KAPPA, 61
LIST, 54, 64, 102
LWALL, 37, 50
MASS-WEIGHTED, 75
MASS, 61
MINDIST, 68, 101
MINUS PI, 62
MM, 96
MSD, 84
NLOOP, 93
NN, 96
NO-WEIGHTED, 75
NOALIGN, 94
NOHILLS, 35, 36, 42
NOPBC, 67, 69, 71, 73, 88, 92, 101
NOSPLINE, 32
NOTE, 26
NO CENTER, 86
NO ROT, 86
ONIONS, 54, 58
PARABETARMSD, 80, 102
PCA, 94, 95
PERIODIC, 62
PLUS 2PI, 62
PLUS PI, 62
POINT FROM AXIS, 68
POLY, 99
POSITION, 66
PRINT, 26
PROJ GRAD, 53, 54
126

PROJ ON AXIS, 68
PTMETAD, 41
PUCKERING, 82
RDF, 97
READ GRID, 33
RECONNAISSANCE, 54, 58, 59
RESTART, 58, 62
RGYR 1, 75
RGYR 2, 76
RGYR 3, 76
RGYR, 74, 75, 101
RMSD, 84
RUN FREQ, 56
R 0, 96
SHIFT, 99
SIGMA, 36, 38, 42, 64
SIMTEMP, 41
SIN, 69, 70
SPRINT, 96
SQRT, 91
STEERPLAN, 46
STEER, 46
STORE FREQ, 56, 57
S PATH, 83–85, 87, 101
TARGETED, 91
TAUTHERMO, 61
TEMPERATURE, 61
TERMS, 99
TORSION, 69, 70
TRACE, 75
TYPE, 75
T sj, 63
UMBRELLA, 45
UWALL, 50
WATERBRIDGE, 73
WIDTH, 58
WRITE GRID, 32
WRITE STATE, 62

W STRIDE, 32, 58
W sj, 63
Z PATH, 83–85, 87, 101
#, 26
sj, 63

Umbrella sampling, 45
Units, 27
Walls, 50
Well-tempered metadynamics Reweight,
111

Minimal distance, see Collective variables, Minimal distance
MSD, see Collective variables, MSD
Output files, Metadynamics, 27
Parabetarmsd, see Collective variables,
Parabetarmsd
Parallel tempering metadynamics, 41
Path variables, see Collective variables,
Path variables
Periodic boundary conditions, 101
plumedat.sh, 27, 108
POLY, see Collective variables, POLY
Projection on an axis, see Collective
variables, Projection on an axis
Puckering coordinates, see Collective
variables, Puckering coordinates
Radius of Gyration, see Collective variables, Radius of Gyration
RDF, see Collective variables, RDF
SPRINT, see Collective variables, SPRINT
Standalone, 109
Steered MD, 46
Steerplan, 46
Tabulated potentials, 50
Targeted, see Collective variables, Targeted
Torsion, see Collective variables, Torsion

127

Source Exif Data:

File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.4
Linearized                      : No
Page Count                      : 128
Page Mode                       : UseOutlines
Author                          : 
Title                           : 
Subject                         : 
Creator                         : LaTeX with hyperref package
Producer                        : pdfeTeX-1.21a
Create Date                     : 2013:02:01 17:38:12+01:00
PTEX Fullbanner                 : This is pdfeTeX, Version 3.141592-1.21a-2.2 (Web2C 7.5.4) kpathsea version 3.5.4

EXIF Metadata provided by EXIF.tools

Manual

Navigation menu

Versions of this User Manual:

Views

Navigation