[gmx-users] large error bars in PMF

[Justin A. Lemkul]

Rebeca García Fandiño wrote:
> 
> Hi,
> I am trying to calculate the binding energy of two molecules using the 
> PMF (Umbrella Sampling method) and Gromacs 4.0.
> Some weeks ago I have written to the list because changing the number of 
> windows used in the Umbrella Sampling calculations different results 
> were obtained, and I was suggested to extend my simulations since the 
> error bars associated to each windows were too high.
> I have now extended my simulations from 1 ns to 8 ns, however, I cannot 
> see much different from the shorter calculations. I send you the 
> comparison of the two PMF including the error bars (attached).
> Now I am using 50 windows, but the shorter simulations were done using 
> 100 windows, so I don't think increasing the number of windows could help.
> My system has about 29200 atoms (where 29000 are chloroform atoms). The 
> *mdp file I am using is copied below.
> Would you have any suggestion to improve the results and decrease the 
> error bars in the calculations?
> 

Neither of the PMF profiles look converged.  In the case of 8-ns sampling, you 
do not have a well-defined energy minimum and the whole plot is very rough. 
This suggests to me that you still need more sampling time.  Keep in mind that 
even 8 ns is a very tiny amount of sampling in each window.

What is it that you're trying to pull?  What are r_1 and r_2?  Are they 
components of a small peptide or something else?  It looks as if the total DG is 
going to be somewhere on the order of 6 kcal/mol, so if you want reasonable 
error bars you will not lots of well-converged data.

-Justin

> ----------------------------MDP file---------------------------
> title       = Umbrella pulling simulation
> define      =
> define      =
> ; Run parameters
> integrator  = md
> dt          = 0.002
> tinit       = 0
> nsteps      = 500000   ; 1 ns
> nstcomm     = 10
> ; Output parameters
> nstxout     = 5000     ; every 10 ps
> nstvout     = 5000
> nstfout     = 5000
> nstxtcout   = 5000      ; every 10 ps
> nstenergy   = 5000
> ; Bond parameters
> constraint_algorithm    = lincs
> constraints             = all-bonds
> continuation            = yes
> ; Single-range cutoff scheme
> nstlist     = 5
> ns_type     = grid
> rlist       = 1.4
> rcoulomb    = 1.4
> rvdw        = 1.4
> ; PME electrostatics parameters
> coulombtype     = PME
> fourierspacing  = 0.12
> fourier_nx      = 0
> fourier_ny      = 0
> fourier_nz      = 0
> pme_order       = 4
> ewald_rtol      = 1e-5
> optimize_fft    = yes
> ; Berendsen temperature coupling is on in two groups
> Tcoupl      = Nose-Hoover
> tc_grps     = ACH   CL3
> tau_t       = 0.5       0.5
> ref_t       = 300       300
> ; Pressure coupling is on
> Pcoupl          = Parrinello-Rahman
> pcoupltype      = isotropic
> tau_p           = 1.0
> compressibility = 4.5e-5
> ref_p           = 1.0
> ; Generate velocities is off
> gen_vel     = no
> ; Periodic boundary conditions are on in all directions
> pbc     = xyz
> ; Long-range dispersion correction
> DispCorr    = EnerPres
> ; Pull code
> pull            = umbrella
> pull_geometry   = distance
> pull_dim        = N N Y
> pull_start      = yes
> pull_ngroups    = 1
> pull_group0     = r_1
> pull_group1     = r_2
> pull_init1      = 0
> pull_rate1      = 0.0
> pull_k1         = 1000      ; kJ mol^-1 nm^-2
> pull_nstxout    = 1000      ; every 2 ps
> pull_nstfout    = 1000      ; every 2 ps
> 
> -----------------------------------------------
> 
> 
> Thanks a lot in advance.
> 
> Best wishes,
> 
> Dr. Rebeca Garcia
> Santiago de Compostela University
> Spain
> 
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================