SECTION 1

Umbrella Sampling

Here we will carry out umbrella sampling at 300K for every umbrella window for about 5 ns. We will run simulations to break the reaction coordinate (in our case dihedral angle Φ) up into the series of windows and apply the restraint that acts to force the reaction coordinate to remain close to centre of window. The selection for number of umbrella windows must be such that every window must sample some space of next window. Similarly, force constant must be large enough to sample only the small subset of space we want to sample but it should not be too large that there is no overlap between the two consecutive windows.

We will apply umbrella bias along Φ from –π to π (-3.2 rad to 3.2 rad for simplicity) at an interval of 0.2 rad. Till now, we have one starting structure (say Φ= -1.2 rad). If we apply restraint on it, to sample whole space of Φ then, in some part of space there will be sudden change in conformation due to sudden change in forces. To avoid this situation, we will start simulation with structure at Φ= -1.2 rad and use the restart (.rst file) from this as the input of -1.4 rad simulation and so on, upto -3.2 rad simulation. Since the value of dihedral angles are periodic, - 3.2 and 3.2 rad have same conformation. Therefore, we will use restart file of 3.2 rad simulation as input of 3.0 rad simulation and so on, upto the -1.0 rad simulation. In this way we will sample Φ by placing 33 restraining umbrella potentials.

The following input files are required for each umbrella simulation:

  1. ala_prmtop : It describes the parameter and topology of the molecules of the system
  2. ala_prod.rst : It describes the coordinates of the molecules of the system
  3. mds.in : It describes the input settings for Amber MD engine
  4. plumed.dat : It describes the input settings for PLUMED-AMBER interface

mds.in

Alanine Di Peptide
 &cntrl
  irest=1, ntx=5,
  imin = 0, ntb = 0,
  igb = 0, ntpr = 1000, ntwx = 1000,
  ntt = 3, gamma_ln = 1.0,
  tempi = 300.0, temp0 = 300.0,
  nstlim = 5000000, dt = 0.001,
  nscm=500,
  cut=999.0, rgbmax=999.0,
plumed=1,
plumedfile='plumed.dat'
&end

plumed.dat

# set up two variables for Phi and Psi dihedral  
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
##################################################
# Impose an umbrella potential on CV 1

restraint-phi: RESTRAINT ARG=phi KAPPA=500 AT=XXXX
##################################################
# monitor the two variables and the bias potential from the two restraints
PRINT STRIDE=10 ARG=phi,psi,restraint-phi.bias FILE=COLVAR

Here, we need to replace XXXX by centre of umbrella window we want to sample. Now, we can run each simulation through sander. For example to run them using 2 processors :

mpirun -np 2 $/AMBERHOME/bin/sander.MPI -O -i mds.in -o umb.out -p ala_prmtop -c ala_prod.rst -r umb.rst -x umb.mdcrd

The output files of these simulations are:

  1. umb.out : output file
  2. umb.rst : output restart file with coordinates
  3. umb.mdcrd : Output trajectory file for MD simulation
  4. COLVAR : Collective Variable file

This will generate 33 structures relaxed at the centre of every umbrella window.

Analysis of output

We should check the COLVAR file to make sure the proper sampling in every umbrella window. We will plot dihedral angle Φ versus time by following xmgrace command:

xmgrace -block COLVAR -bxy 1:2

It is showing nice sampling for window corresponding to Φ = -1.2. Similarly, we can check the sampling for all other umbrella windows. Now, will check the overlap between two consecutive windows. It is suggested that for good overlap the position of next window to be sampled must be chosen from previous window to match their estimated half maxima.

It is showing nice overlap, so now we can proceed to patch these umbrella and generated free energy surface.