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 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.

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
#SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
# vim:ft=plumed
MOLINFO STRUCTUREcompulsory keyword
a file in pdb format containing a reference structure. =diala.pdb
phi: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-2
psi: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-2
hhphi: HISTOGRAM ARGthe input for this action is the scalar output from one or more other actions. =phi STRIDEcompulsory keyword ( default=1 )
the frequency with which the data should be collected and added to the quantity being
averaged =50 GRID_MINcompulsory keyword
the lower bounds for the grid =-pi GRID_MAXcompulsory keyword
the upper bounds for the grid =pi GRID_BINthe number of bins for the grid =600 BANDWIDTHcompulsory keyword
the bandwidths for kernel density estimation =0.1
hhpsi: HISTOGRAM ARGthe input for this action is the scalar output from one or more other actions. =psi STRIDEcompulsory keyword ( default=1 )
the frequency with which the data should be collected and added to the quantity being
averaged =50 GRID_MINcompulsory keyword
the lower bounds for the grid =-pi GRID_MAXcompulsory keyword
the upper bounds for the grid =pi GRID_BINthe number of bins for the grid =600 BANDWIDTHcompulsory keyword
the bandwidths for kernel density estimation =0.1
ffphi: CONVERT_TO_FES GRIDcompulsory keyword
the action that creates the input grid you would like to use =hhphi
ffpsi: CONVERT_TO_FES GRIDcompulsory keyword
the action that creates the input grid you would like to use =hhpsi
DUMPGRID GRIDcompulsory keyword
the action that creates the grid you would like to output =ffphi FILEcompulsory keyword ( default=density )
the file on which to write the grid. =fes_phi STRIDEcompulsory keyword ( default=0 )
the frequency with which the grid should be output to the file. =100000
DUMPGRID GRIDcompulsory keyword
the action that creates the grid you would like to output =ffpsi FILEcompulsory keyword ( default=density )
the file on which to write the grid. =fes_psi STRIDEcompulsory keyword ( default=0 )
the frequency with which the grid should be output to the file. =100000
PRINT ARGthe input for this action is the scalar output from one or more other actions. =phi,psi FILEthe name of the file on which to output these quantities =colvar.dat STRIDEcompulsory 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 #SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb # vim:ft=plumed MOLINFO STRUCTUREcompulsory keyword a file in pdb format containing a reference structure. =diala.pdb phi: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-2 psi: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-2 doubleg: CUSTOM ... ARGthe input for this action is the scalar output from one or more other actions. =phi FUNCcompulsory keyword the function you wish to evaluate =exp(__FILL__*cos(x-__FILL__))+exp(__FILL__*cos(x-__FILL__)) PERIODICcompulsory keyword if the output of your function is periodic then you should specify the periodicity of the function. =NO ... b: BIASVALUE ARGthe input for this action is the scalar output from one or more other actions. =doubleg hhphi: HISTOGRAM ARGthe input for this action is the scalar output from one or more other actions. =phi STRIDEcompulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged =50 GRID_MINcompulsory keyword the lower bounds for the grid =-pi GRID_MAXcompulsory keyword the upper bounds for the grid =pi GRID_BINthe number of bins for the grid =600 BANDWIDTHcompulsory keyword the bandwidths for kernel density estimation =0.1 hhpsi: HISTOGRAM ARGthe input for this action is the scalar output from one or more other actions. =psi STRIDEcompulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged =50 GRID_MINcompulsory keyword the lower bounds for the grid =-pi GRID_MAXcompulsory keyword the upper bounds for the grid =pi GRID_BINthe number of bins for the grid =600 BANDWIDTHcompulsory keyword the bandwidths for kernel density estimation =0.1 ffphi: CONVERT_TO_FES GRIDcompulsory keyword the action that creates the input grid you would like to use =hhphi ffpsi: CONVERT_TO_FES GRIDcompulsory keyword the action that creates the input grid you would like to use =hhpsi DUMPGRID GRIDcompulsory keyword the action that creates the grid you would like to output =ffphi FILEcompulsory keyword ( default=density ) the file on which to write the grid. =fes_phi STRIDEcompulsory keyword ( default=0 ) the frequency with which the grid should be output to the file. =100000 DUMPGRID GRIDcompulsory keyword the action that creates the grid you would like to output =ffpsi FILEcompulsory keyword ( default=density ) the file on which to write the grid. =fes_psi STRIDEcompulsory keyword ( default=0 ) the frequency with which the grid should be output to the file. =100000 PRINT ARGthe input for this action is the scalar output from one or more other actions. =phi,psi,b.bias FILEthe name of the file on which to output these quantities =colvar.dat STRIDEcompulsory 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 #SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb MOLINFO STRUCTUREcompulsory keyword a file in pdb format containing a reference structure. =diala.pdb phi: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-2 psi: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-2 tripleg: CUSTOM ... ARGthe input for this action is the scalar output from one or more other actions. =phi FUNCcompulsory keyword the function you wish to evaluate =exp(__FILL__*cos(x-__FILL__))+exp(__FILL__*cos(x-__FILL__))+exp(__FILL__*cos(x-__FILL__)) PERIODICcompulsory keyword if the output of your function is periodic then you should specify the periodicity of the function. =NO ... b: BIASVALUE ARGthe input for this action is the scalar output from one or more other actions. =tripleg hhphi: HISTOGRAM ARGthe input for this action is the scalar output from one or more other actions. =phi STRIDEcompulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged =50 GRID_MINcompulsory keyword the lower bounds for the grid =-pi GRID_MAXcompulsory keyword the upper bounds for the grid =pi GRID_BINthe number of bins for the grid =600 BANDWIDTHcompulsory keyword the bandwidths for kernel density estimation =0.1 hhpsi: HISTOGRAM ARGthe input for this action is the scalar output from one or more other actions. =psi STRIDEcompulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged =50 GRID_MINcompulsory keyword the lower bounds for the grid =-pi GRID_MAXcompulsory keyword the upper bounds for the grid =pi GRID_BINthe number of bins for the grid =600 BANDWIDTHcompulsory keyword the bandwidths for kernel density estimation =0.1 ffphi: CONVERT_TO_FES GRIDcompulsory keyword the action that creates the input grid you would like to use =hhphi ffpsi: CONVERT_TO_FES GRIDcompulsory keyword the action that creates the input grid you would like to use =hhpsi DUMPGRID GRIDcompulsory keyword the action that creates the grid you would like to output =ffphi FILEcompulsory keyword ( default=density ) the file on which to write the grid. =fes_phi STRIDEcompulsory keyword ( default=0 ) the frequency with which the grid should be output to the file. =100000 DUMPGRID GRIDcompulsory keyword the action that creates the grid you would like to output =ffpsi FILEcompulsory keyword ( default=density ) the file on which to write the grid. =fes_psi STRIDEcompulsory keyword ( default=0 ) the frequency with which the grid should be output to the file. =100000 PRINT ARGthe input for this action is the scalar output from one or more other actions. =phi,psi,b.bias FILEthe name of the file on which to output these quantities =colvar.dat STRIDEcompulsory 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 #SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb # vim:ft=plumed MOLINFO STRUCTUREcompulsory keyword a file in pdb format containing a reference structure. =diala.pdb phi: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-2 psi: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-2 tripleg: MATHEVAL ... ARGthe input for this action is the scalar output from one or more other actions. =phi FUNCcompulsory keyword the function you wish to evaluate =__FILL__ PERIODICcompulsory keyword if the output of your function is periodic then you should specify the periodicity of the function. =NO ... b: BIASVALUE ARGthe input for this action is the scalar output from one or more other actions. =tripleg as: REWEIGHT_BIAS ARGcompulsory keyword ( default=*.bias ) the biases that must be taken into account when reweighting =b.bias hhphi: HISTOGRAM ARGthe input for this action is the scalar output from one or more other actions. =phi STRIDEcompulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged =50 GRID_MINcompulsory keyword the lower bounds for the grid =-pi GRID_MAXcompulsory keyword the upper bounds for the grid =pi GRID_BINthe number of bins for the grid =600 BANDWIDTHcompulsory keyword the bandwidths for kernel density estimation =0.1 LOGWEIGHTSlist of actions that calculates log weights that should be used to weight configurations when calculating averages =as hhpsi: HISTOGRAM ARGthe input for this action is the scalar output from one or more other actions. =psi STRIDEcompulsory keyword ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged =50 GRID_MINcompulsory keyword the lower bounds for the grid =-pi GRID_MAXcompulsory keyword the upper bounds for the grid =pi GRID_BINthe number of bins for the grid =600 BANDWIDTHcompulsory keyword the bandwidths for kernel density estimation =0.1 LOGWEIGHTSlist of actions that calculates log weights that should be used to weight configurations when calculating averages =as ffphi: CONVERT_TO_FES GRIDcompulsory 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 GRIDcompulsory keyword the action that creates the input grid you would like to use =hhpsi STRIDE could not find this keyword =100000 DUMPGRID GRIDcompulsory keyword the action that creates the grid you would like to output =ffphi FILEcompulsory keyword ( default=density ) the file on which to write the grid. =ffphi.dat DUMPGRID GRIDcompulsory keyword the action that creates the grid you would like to output =ffpsi FILEcompulsory keyword ( default=density ) the file on which to write the grid. =ffpsi.dat PRINT ARGthe input for this action is the scalar output from one or more other actions. =phi,psi,b.bias FILEthe name of the file on which to output these quantities =colvar.dat STRIDEcompulsory 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
#SETTINGS MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
# vim:ft=plumed
MOLINFO STRUCTUREcompulsory keyword
a file in pdb format containing a reference structure. =diala.pdb
phi: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-2
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-3.00
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-2.75
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-2.50
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-2.25
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-2.00
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-1.75
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-1.50
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-1.25
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-1.00
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-0.75
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-0.50
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =-0.25
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =0.00
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =0.25
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =0.50
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =0.75
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =1.00
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =1.25
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =1.50
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =1.75
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =2.00
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =2.25
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =2.50
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =2.75
RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0 ATcompulsory keyword
the position of the restraint =3.00
PRINT ARGthe input for this action is the scalar output from one or more other actions. =*.bias FILEthe name of the file on which to output these quantities =biases.dat STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =10
PRINT ARGthe input for this action is the scalar output from one or more other actions. =phi FILEthe name of the file on which to output these quantities =allphi.dat STRIDEcompulsory 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
#SETTINGS FILENAME=plumed.dat MOLFILE=user-doc/tutorials/lugano-2/diala.pdb
# this is plumed.dat
# vim:ft=plumed
MOLINFO STRUCTUREcompulsory keyword
a file in pdb format containing a reference structure. =diala.pdb
phi: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-2
psi: TORSION ATOMSthe 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 ...
ARGthe input for this action is the scalar output from one or more other actions. =phi KAPPAcompulsory keyword ( default=0.0 )
specifies that the restraint is harmonic and what the values of the force constants
on each of the variables are =250.0
ATcompulsory 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 STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =20 ARGthe input for this action is the scalar output from one or more other actions. =phi,psi,restraint-phi.bias FILEthe 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

And again use the driver in parallel:

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