[gmx-users] Umbrella Sampling Results and Issues

[Agnivo Gosai]

Dear Users

I am doing an umbrella sampling MD to calculate the binding energy between
a DNA and a protein. I am following Dr. Lemkuls tutorial for the same.

Firstly I put the complex in a rectangular box and used the
AMBER99SB-Parmbsc0 forcefield for the topology generation. Then I used SD
for minimization , 1 ns restrained NVT equilibration ( V-rescale , T = 300
K) , 1 ns restrained NPT equilibration ( Berendsen couplings in both , ref
Press = 1 bar , however achieved average press = 0.4 bar with very low
error , RMSD ) , 500 ps Steered Molecular Dynamics keeping DNA fixed and
pulling Protein ( Using Nose-Hoover thermostat , Parinello - Rahman
barostat ).

I pull for a total of 5 nm.

This is the relevant mdp entry for the SMD

; Pull code
pull            = umbrella
pull_geometry   = distance  ; simple distance increase
pull_dim        = N N Y     ; Z axis reaction coordinate
pull_start      = yes       ; define initial COM distance > 0
pull_ngroups    = 1         ; applying puling force to protein only
pull_group0     = DNA       ; reference fixed group for pulling
pull_group1     = Protein   ; the molecule to which pulling force is applied
pull_rate1      = 0.01      ; 0.01 nm per ps = 10 nm per ns
pull_k1         = 1000      ; kJ mol^-1 nm^-2 force constant for pulling

After this I extracted the conf.gro files and using the perl/python scripts
from the tutorial extracted 21 windows (0.2 nm spacing) from the SMD
trajectory for the umbrella sampling. In each window I did a 100 ps short
NPT equilibration restraining only the DNA and then an 1 ns umbrella
sampling simulation restraining the DNA with the following mdp entry :

; Pull code
pull            = umbrella
pull_geometry   = distance
pull_dim        = N N Y
pull_start      = yes
pull_ngroups    = 1
pull_group0     = DNA
pull_group1     = Protein
pull_init1      = 0
pull_rate1      = 0.0
pull_k1         = 1000      ; kJ mol^-1 nm^-2
pull_nstxout    = 100      ; every 0.2 ps
pull_nstfout    = 100      ; every 0.2 ps

Then I ran the g_wham tool on the results.

I am attaching the PMF profile and the Histogram for reference :

​
 PMF
<https://docs.google.com/file/d/0B-U8uULVZjfRaHE4STc2QVo3VmM/edit?usp=drive_web>
​​
 HISTOGRAM
<https://docs.google.com/file/d/0B-U8uULVZjfRZ2IxbXZSeTc1QUE/edit?usp=drive_web>
​
As observed , the PMF is having 3 steps for particular reaction coordinates
and at those regions the Histogram is also not showing ideal overlapping
behaviour.

This is the ERROR output of g_wham :

WARNING, no data point in bin 26 (z=3.04567) !
You may not get a reasonable profile. Check your histograms!

The error file is also attached :

​
 WHAM1_ERROR.txt
<https://docs.google.com/file/d/0B-U8uULVZjfRbkllV1hZLXR6MjQ/edit?usp=drive_web>
​
So, my questions are :

Should I decrease the nm spacing to say 0.1 nm around those particular
reaction coordinates where the "no data point in bin" error is coming ? ?
Should I be going for asymmetric windows

OR , should I simultaneously increase my simulation time for umbrella
sampling to say 5 ns from 1 ns ?

Kindly suggest

------------------------------------------------------------------------


Thanks & Regards
Agnivo Gosai
Grad Student, Iowa State University.