Lugano tutorial: Using restraints

Aims

The aim of this tutorial is to introduce the users to the use of constant biases in PLUMED and their application to perform a simple Umbrella Sampling WHAM simulation.

Objectives

  • Apply a restraint on a simulations over one or more collective variables
  • Understand the effect of a restraint on the acquired statistics
  • Perform a simple un-biasing of a restrained simulation
  • Add an external potential in the form of an analytical or numerical function
  • Perform an Umbrella Sampling WHAM simulation

Resources

The TARBALL for this tutorial contains the following files:

  • diala.pdb: a PDB file for alanine dipeptide in vacuo
  • topol.tpr: a GROMACS run file to perform MD of alanine dipeptide
  • run_us.sh: a script to perform a serial US wham simulation
  • wham.py: a script to perform the WHAM analysis
  • do_fes.py: a script to generate fes from CV and weights

Also notice that the .solutions direction of the tarball contains correct input files for the exercises. This tutorial has been tested with version 2.5.

Introduction

PLUMED can calculate conformational properties of a system a posteriori as well as on-the-fly. This information can be use to manipulate a simulation on-the-fly. This means adding energy terms in addition to those of the original Hamiltonian. This additional energy terms are usually referred as Bias. In the following we will see how to apply a constant bias potential with PLUMED. It is preferable to run each exercise in a separate folder.

To learn more: Summary of theory

We will make use as a toy-model of alanine dipeptide: we will see how we can use an iterative approach to build a constant bias to speed up the sampling.

Note
Create a folder for each exercise and use sub-folders if you want to run the same simulation with multiple choices for the parameters

Exercise 1: Preliminary run with alanine dipeptide

Alanine dipeptide is characterized by multiple minima separated by relatively high free energy barriers. Here we will explore the conformational space of alanine dipeptide using a standard MD simulation, then instead of using the free energy as an external potential we will try to fit the potential using gnuplot and add a bias using an analytical function of a collective variable with CUSTOM and BIASVALUE .

As a first test lets run an MD and generate on-the-fly the free energy as a function of the phi and psi collective variables separately.

This is an example input file to calculate the phi and psi angles on the fly and accumulate two 1D histograms from which calculating the free energy.

Click on the labels of the actions for more information on what each action computes
tested on v2.9
#SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
# vim:ft=plumed
MOLINFO 
STRUCTURE
compulsory keyword a file in pdb format containing a reference structure.
=diala.pdb phi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@phi-2 psi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@psi-2 hhphi: HISTOGRAM
ARG
the input for this action is the scalar output from one or more other actions.
=phi
STRIDE
compulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged
=50
GRID_MIN
compulsory keyword the lower bounds for the grid
=-pi
GRID_MAX
compulsory keyword the upper bounds for the grid
=pi
GRID_BIN
the number of bins for the grid
=600
BANDWIDTH
compulsory keyword the bandwidths for kernel density estimation
=0.1 hhpsi: HISTOGRAM
ARG
the input for this action is the scalar output from one or more other actions.
=psi
STRIDE
compulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged
=50
GRID_MIN
compulsory keyword the lower bounds for the grid
=-pi
GRID_MAX
compulsory keyword the upper bounds for the grid
=pi
GRID_BIN
the number of bins for the grid
=600
BANDWIDTH
compulsory keyword the bandwidths for kernel density estimation
=0.1 ffphi: CONVERT_TO_FES
GRID
compulsory keyword the action that creates the input grid you would like to use
=hhphi ffpsi: CONVERT_TO_FES
GRID
compulsory keyword the action that creates the input grid you would like to use
=hhpsi DUMPGRID
GRID
compulsory keyword the action that creates the grid you would like to output
=ffphi
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=fes_phi
STRIDE
compulsory keyword ( default=0 ) the frequency with which the grid should be output to the file.
=100000 DUMPGRID
GRID
compulsory keyword the action that creates the grid you would like to output
=ffpsi
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=fes_psi
STRIDE
compulsory keyword ( default=0 ) the frequency with which the grid should be output to the file.
=100000 PRINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi,psi
FILE
the name of the file on which to output these quantities
=colvar.dat
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=50

run it with gromacs as

gmx mdrun -s topol -plumed plumed.dat -nb cpu -v

from the colvar file it is clear that we can quickly explore two minima but that the region for positive phi is not accessible. Ideally we would like to speed up the sampling of regions that are not visited spontaneously by MD. We have multiple possibilities. One option could be to use as a bias the opposite of the accumulated free-energy using EXTERNAL . Another option can be to fit the FES and use the fit. This is what we will do, but first of all take a look at the fes accumulated in time.

>gnuplot
plot for [i=0:9] 'analysis.'.i.'.fes_phi' u 1:2 w l t''.i
rep 'fes_phi' u 1:2 w l t'final'
plot for [i=0:9] 'analysis.'.i.'.fes_psi' u 1:2 w l t''.i
rep 'fes_psi' u 1:2 w l t'final'

So first we need to fit the opposite of the free energy as a function of phi in the region explored with a periodic function, because of the gaussian like look of the minima we can fit it using the von Mises distribution. In gnuplot

>gnuplot
gnuplot>plot 'fes_phi' u 1:(-$2) w l

Now find a value such as the fes is always positive, e.g. ~38

gnuplot>plot 'fes_phi' u 1:(-$2+38) w l
gnuplot>f(x)=exp(k1*cos(x-a1))+exp(k2*cos(x-a2))
gnuplot>k1=2
gnuplot>k2=2
gnuplot>fit [-2.9:-0.7] f(x) 'fes_phi' u 1:(-$2+38) via k1,a1,k2,a2
gnuplot>rep f(x)

The function and the resulting parameters can be used to run a new biased simulation:

Exercise 2: First biased run with alanine dipeptide

To the above file we add a few lines to define using CUSTOM a function of the angle phi.

Click on the labels of the actions for more information on what each action computes
tested on v2.9
#SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
# vim:ft=plumed
MOLINFO 
STRUCTURE
compulsory keyword a file in pdb format containing a reference structure.
=diala.pdb phi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@phi-2 psi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@psi-2 doubleg: CUSTOM ...
ARG
the input for this action is the scalar output from one or more other actions.
=phi
FUNC
compulsory keyword the function you wish to evaluate
=exp(__FILL__*cos(x-__FILL__))+exp(__FILL__*cos(x-__FILL__))
PERIODIC
compulsory keyword if the output of your function is periodic then you should specify the periodicity of the function.
=NO ... b: BIASVALUE
ARG
the input for this action is the scalar output from one or more other actions.
=doubleg hhphi: HISTOGRAM
ARG
the input for this action is the scalar output from one or more other actions.
=phi
STRIDE
compulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged
=50
GRID_MIN
compulsory keyword the lower bounds for the grid
=-pi
GRID_MAX
compulsory keyword the upper bounds for the grid
=pi
GRID_BIN
the number of bins for the grid
=600
BANDWIDTH
compulsory keyword the bandwidths for kernel density estimation
=0.1 hhpsi: HISTOGRAM
ARG
the input for this action is the scalar output from one or more other actions.
=psi
STRIDE
compulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged
=50
GRID_MIN
compulsory keyword the lower bounds for the grid
=-pi
GRID_MAX
compulsory keyword the upper bounds for the grid
=pi
GRID_BIN
the number of bins for the grid
=600
BANDWIDTH
compulsory keyword the bandwidths for kernel density estimation
=0.1 ffphi: CONVERT_TO_FES
GRID
compulsory keyword the action that creates the input grid you would like to use
=hhphi ffpsi: CONVERT_TO_FES
GRID
compulsory keyword the action that creates the input grid you would like to use
=hhpsi DUMPGRID
GRID
compulsory keyword the action that creates the grid you would like to output
=ffphi
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=fes_phi
STRIDE
compulsory keyword ( default=0 ) the frequency with which the grid should be output to the file.
=100000 DUMPGRID
GRID
compulsory keyword the action that creates the grid you would like to output
=ffpsi
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=fes_psi
STRIDE
compulsory keyword ( default=0 ) the frequency with which the grid should be output to the file.
=100000 PRINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi,psi,b.bias
FILE
the name of the file on which to output these quantities
=colvar.dat
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=50

It is now possible to run a second simulation and observe the new behavior. The system quickly explores a new minimum. While a quantitative estimate of the free energy difference of the old and new regions is out of the scope of the current exercise what we can do is to add a new von Mises function centered in the new minimum with a comparable height, in this way we can hope to facilitate a back and forth transition along the phi collective variable. Look at the old and new free energy and add a third exponential function to CUSTOM centered in the new minimum.

gnuplot> plot 'fes_phi' u 1:(-$2+38) w l
gnuplot> f(x)=exp(k3*cos(x-a3))
gnuplot>k3=2
gnuplot> fit [0.3:1.8] f(x) 'fes_phi' u 1:(-$2+38) via k3,a3

We can now run a third simulation where both regions are biased.

Exercise 3: Second biased run with alanine dipeptide

Click on the labels of the actions for more information on what each action computes
tested on v2.9
#SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
MOLINFO 
STRUCTURE
compulsory keyword a file in pdb format containing a reference structure.
=diala.pdb phi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@phi-2 psi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@psi-2 tripleg: CUSTOM ...
ARG
the input for this action is the scalar output from one or more other actions.
=phi
FUNC
compulsory keyword the function you wish to evaluate
=exp(__FILL__*cos(x-__FILL__))+exp(__FILL__*cos(x-__FILL__))+exp(__FILL__*cos(x-__FILL__))
PERIODIC
compulsory keyword if the output of your function is periodic then you should specify the periodicity of the function.
=NO ... b: BIASVALUE
ARG
the input for this action is the scalar output from one or more other actions.
=tripleg hhphi: HISTOGRAM
ARG
the input for this action is the scalar output from one or more other actions.
=phi
STRIDE
compulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged
=50
GRID_MIN
compulsory keyword the lower bounds for the grid
=-pi
GRID_MAX
compulsory keyword the upper bounds for the grid
=pi
GRID_BIN
the number of bins for the grid
=600
BANDWIDTH
compulsory keyword the bandwidths for kernel density estimation
=0.1 hhpsi: HISTOGRAM
ARG
the input for this action is the scalar output from one or more other actions.
=psi
STRIDE
compulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged
=50
GRID_MIN
compulsory keyword the lower bounds for the grid
=-pi
GRID_MAX
compulsory keyword the upper bounds for the grid
=pi
GRID_BIN
the number of bins for the grid
=600
BANDWIDTH
compulsory keyword the bandwidths for kernel density estimation
=0.1 ffphi: CONVERT_TO_FES
GRID
compulsory keyword the action that creates the input grid you would like to use
=hhphi ffpsi: CONVERT_TO_FES
GRID
compulsory keyword the action that creates the input grid you would like to use
=hhpsi DUMPGRID
GRID
compulsory keyword the action that creates the grid you would like to output
=ffphi
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=fes_phi
STRIDE
compulsory keyword ( default=0 ) the frequency with which the grid should be output to the file.
=100000 DUMPGRID
GRID
compulsory keyword the action that creates the grid you would like to output
=ffpsi
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=fes_psi
STRIDE
compulsory keyword ( default=0 ) the frequency with which the grid should be output to the file.
=100000 PRINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi,psi,b.bias
FILE
the name of the file on which to output these quantities
=colvar.dat
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=50

With this third simulation it should be possible to visit both regions as a function on the phi torsion. The resulting free energy is now reporting about the biased simulation is flatter than the former even if not flat everywhere. Now it is possible to reweight the sampling and obtain a better free energy estimate along phi.

Click on the labels of the actions for more information on what each action computes
tested on v2.9
#SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
# vim:ft=plumed
MOLINFO 
STRUCTURE
compulsory keyword a file in pdb format containing a reference structure.
=diala.pdb phi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@phi-2 psi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@psi-2 tripleg: MATHEVAL ...
ARG
the input for this action is the scalar output from one or more other actions.
=phi
FUNC
compulsory keyword the function you wish to evaluate
=__FILL__
PERIODIC
compulsory keyword if the output of your function is periodic then you should specify the periodicity of the function.
=NO ... b: BIASVALUE
ARG
the input for this action is the scalar output from one or more other actions.
=tripleg as: REWEIGHT_BIAS
ARG
compulsory keyword ( default=*.bias ) the biases that must be taken into account when reweighting
=b.bias hhphi: HISTOGRAM
ARG
the input for this action is the scalar output from one or more other actions.
=phi
STRIDE
compulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged
=50
GRID_MIN
compulsory keyword the lower bounds for the grid
=-pi
GRID_MAX
compulsory keyword the upper bounds for the grid
=pi
GRID_BIN
the number of bins for the grid
=600
BANDWIDTH
compulsory keyword the bandwidths for kernel density estimation
=0.1
LOGWEIGHTS
list of actions that calculates log weights that should be used to weight configurations when calculating averages
=as hhpsi: HISTOGRAM
ARG
the input for this action is the scalar output from one or more other actions.
=psi
STRIDE
compulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged
=50
GRID_MIN
compulsory keyword the lower bounds for the grid
=-pi
GRID_MAX
compulsory keyword the upper bounds for the grid
=pi
GRID_BIN
the number of bins for the grid
=600
BANDWIDTH
compulsory keyword the bandwidths for kernel density estimation
=0.1
LOGWEIGHTS
list of actions that calculates log weights that should be used to weight configurations when calculating averages
=as ffphi: CONVERT_TO_FES
GRID
compulsory keyword the action that creates the input grid you would like to use
=hhphi
STRIDE
could not find this keyword
=100000 ffpsi: CONVERT_TO_FES
GRID
compulsory keyword the action that creates the input grid you would like to use
=hhpsi
STRIDE
could not find this keyword
=100000 DUMPGRID
GRID
compulsory keyword the action that creates the grid you would like to output
=ffphi
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=ffphi.dat DUMPGRID
GRID
compulsory keyword the action that creates the grid you would like to output
=ffpsi
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=ffpsi.dat PRINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi,psi,b.bias
FILE
the name of the file on which to output these quantities
=colvar.dat
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=50

Of notice that the reweighting is also applied to the psi collective variable. If you have performed your simulations in different folder you can compare the effect of the bias on phi on the free energy of psi. For a single simulation with a constant bias the reweighting is simple, the weight of each frame is exp(+bias(cv(t))/kt). So it is possible to perform the reweighting by hand at any time.

Now you have performed an original Umbrella Sampling calculation. This is not particularly easy to setup nor robust, even if from a modern perspective it is a very rough implementation of METAD

In the next exercise we will perform a WHAM Umbrella Sampling simulation.

Exercise 4: WHAM Umbrella Sampling

In this case we will run many simulations with a strong harmonic restraint centered around specific values of phi in such a way to cover all possible values, keep each simulation close to its specific value, allow for overlap between neighbor simulations, i.e. simulations centered around consecutive phi values. The simulation can be either performed in parallel by preparing starting configurations close to each value or sequentially, extracting a good starting conformation from the former simulations. In the specific case of alanine dipeptide we can even just start always from the same configuration and let the bias quickly move it close to the target values.

To run the simulation in scalar you can make use of the provided bash script that is:

for AT in -3.00  -2.75 -2.50 -2.25 -2.00 -1.75 -1.50 -1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00
do

cat >plumed.dat << EOF
# vim:ft=plumed
MOLINFO STRUCTURE=diala.pdb
phi: TORSION ATOMS=@phi-2
psi: TORSION ATOMS=@psi-2
#
# Impose an umbrella potential on CV 1
# with a spring constant of 250 kjoule/mol
# at fixed points along phi 
#
restraint-phi: RESTRAINT ARG=phi KAPPA=250.0 AT=$AT
# monitor the two variables and the bias potential from the two restraints
PRINT STRIDE=20 ARG=phi,psi,restraint-phi.bias FILE=COLVAR$AT
EOF

gmx mdrun -plumed plumed.dat -nsteps 100000 -x traj$AT.xtc -c cout$AT.gro -nb cpu

done

you can run it using

./run_us.sh

Plotting the phi collective variable for all replica you will see that each simulation has explored a well defined region of the conformation space as defined by phi. To perform the WHAM merging of the windows we need to 1) collect all the frames

gmx trjcat -f traj*.xtc -cat -o concatenated.xtc

2) calculate the values for all employed biases applied on each frame for this we can write a plumed-wham.dat file including all the biases used in the former simulations:

Click on the labels of the actions for more information on what each action computes
tested on v2.9
#SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
# vim:ft=plumed
MOLINFO 
STRUCTURE
compulsory keyword a file in pdb format containing a reference structure.
=diala.pdb phi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@phi-2 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-3.00 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-2.75 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-2.50 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-2.25 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-2.00 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-1.75 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-1.50 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-1.25 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-1.00 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-0.75 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-0.50 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=-0.25 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=0.00 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=0.25 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=0.50 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=0.75 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=1.00 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=1.25 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=1.50 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=1.75 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=2.00 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=2.25 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=2.50 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=2.75 RESTRAINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=3.00 PRINT
ARG
the input for this action is the scalar output from one or more other actions.
=*.bias
FILE
the name of the file on which to output these quantities
=biases.dat
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=10 PRINT
ARG
the input for this action is the scalar output from one or more other actions.
=phi
FILE
the name of the file on which to output these quantities
=allphi.dat
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=10
plumed driver --mf_xtc concatenated.xtc --plumed plumed-wham.dat

3) run the iterative WHAM optimization and get a weight per frame

python wham.py biases.dat 25 2.49

where 25 is the number of windows and 2.49 is the temperature in energy unit. After some time the result is a file weight.dat with one weight per frame that can be used to calculate any possible property of the system. For example the free energy profile along phi.

To do so edit the weight.dat file to add 3 blank lines and then

paste allphi.dat weights.dat | grep -v \# > allphi-w.dat
python do_fes.py allphi-w.dat 1 -3.1415 3.1415 50 2.49 fes.dat

the resulting profile will be disappointing, error estimate and convergence will be discussed in the following tutorials, but clearly simulations are too short. A more advanced approach would be to use the configurations obtained from the former simulations to generate multiple replicas and then perform the US again for longer time and possible in parallel. The syntax is presented in the following but the exercise is possible only if plumed is compiled with mpi

Exercise 5: WHAM Umbrella Sampling in parallel (optional)

Here we use the "replica" syntax of plumed to write a single plumed input file for all the windows:

Click on the labels of the actions for more information on what each action computes
tested on v2.9
#SETTINGS FILENAME=plumed.dat MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
# this is plumed.dat
# vim:ft=plumed
MOLINFO 
STRUCTURE
compulsory keyword a file in pdb format containing a reference structure.
=diala.pdb phi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@phi-2 psi: TORSION
ATOMS
the four atoms involved in the torsional angle
=@psi-2 # # Impose an umbrella potential on CV 1 # with a spring constant of 250 kjoule/mol # at fixed points along phi # restraint-phi: RESTRAINT ...
ARG
the input for this action is the scalar output from one or more other actions.
=phi
KAPPA
compulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are
=250.0
AT
compulsory keyword the position of the restraint
=@replicas:{ -3.00 -2.75 -2.50 -2.25 -2.00 -1.75 -1.50 -1.25 -1.00 -0.75 -0.50 -0.25 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 } ... # monitor the two variables and the bias potential from the restraint PRINT
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=20
ARG
the input for this action is the scalar output from one or more other actions.
=phi,psi,restraint-phi.bias
FILE
the name of the file on which to output these quantities
=COLVAR

The @replicas syntax allow to define different values for a variable for the different replicas.

mpiexec -np 25 gmx_mpi -s topol -plumed plumed.dat -multi 25 -replex 100 -nb cpu -nsteps 100000 

In this case we run 25 parallel simulations and we also try to perform replica-exchange between neighbor replicas.

Once the simulation is finished the trajectories can be concatenated and analyzed with WHAM making use of the plumed native implementation:

gmx_mpi trjcat -f traj*.xtc -o concatenated.xtc -cat

Write a new plumed-wham.dat

Click on the labels of the actions for more information on what each action computes
tested on v2.9




And again use the driver in parallel:

mpiexec -np 25 plumed driver --mf_xtc concatenated.xtc --plumed plumed-wham.dat --multi 25