[gmx-users] umbrella sampling

[chris.neale at utoronto.ca]

You need to evaluate convergence yourself for any simulation. I  
suggest doing the whole thing twice (or more) with different starting  
conformations. Also, look at the PMF from block averaging (generate  
one PMF from the 0-2 ns data, another from the 2-4 ns data, and so on)  
and see if there is a systematic drift with time.

In my opinion, if you really want to converge to within 5 kcal/mol  
without using additional restraints as Justin mentioned, then you  
should be expecting to do >>100 ns per umbrella, and probably >>1000  
ns per umbrella. -- Thus additional restraints are a very good idea.  
To see about using additional restraints, check out this paper: THE  
JOURNAL OF CHEMICAL PHYSICS 125, 084902 (
2006) and this one: Journal  
of Chemical Theory and Computation 3(4):1231-1235 (2007). They are not  
US, but the idea is the same.

Chris.

-- original message --

Hi Justin,

Here is my complete code for umbrella sampling

title       = Umbrella pulling simulation
define      = -DPOSRES_2

integrator  = md
dt          = 0.002
tinit       = 0
nsteps      = 2500000   ; 5 ns
nstcomm     = 10
nstxout     = 50000     ; every 100 ps
nstvout     = 50000
nstfout     = 5000
nstxtcout   = 5000      ; every 10 ps
nstenergy   = 5000

constraint_algorithm    = lincs
constraints             = all-bonds
continuation            = yes

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      = V-rescale
tc_grps     = Protein   SOL
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     = Chain_B
pull_group1     = Chain_A
pull_init1      = 0
pull_rate1      = 0.0
pull_k1         = 1000
pull_nstxout    = 1000
pull_nstfout    = 1000


I will do this umbrella sampling for 10ns,
Do you think the COM distance can converge after long simulation time? ie
back to the restraint position.

Regards,
Vijay.




[Hide Quoted Text]
Vijayaraj wrote:
Hello,

I am doing umbrella sampling to calculate PMF curve for the detachment
of a terminal cyclic peptide with 8 aa's (CP) from the self-assembled
cyclic peptide nanotube. I extracted the reaction coordinates staring
from 5.5 ang COM distance between the terminal CP and the 2nd CP to 17.5
ang, I did umbrella sampling on 25 configurations (for 5ns) and the
window size is 0.5 ang. I used the following pulling code for umbrella
sampling,

pull            = umbrella
pull_geometry   = distance
pull_dim        = N N Y
pull_start      = yes
pull_ngroups    = 1
pull_group0     = Chain_B
pull_group1     = Chain_A
pull_init1      = 0
pull_rate1      = 0.0
pull_k1         = 1000
pull_nstxout    = 1000
pull_nstfout    = 1000


I restrained the 2nd CP unit and the pull_rate1 is 0, so the COM
distance between the chain_A (terminal) and chain_B (2nd CP) should be
restrained. after 5ns of umbrella sampling, I calculated the COM
distance between chain A and B, but it was not restrained, for the 5.5
ang COM distance restrain, the COM distance varies from 4.5 to 5.5 ang.
and also the pulling cyclic peptide undergoes large conformational
sampling. from the WHAM analysis I understood the sampling window is
poorly represented. In addition to COM distance restrain, can I restrain
the pulling CP's backbone atoms? so that the pulling groups large
conformational sampling will be reduced.
You could implement dihedral restraints to fix the backbone secondary
structure,
but I can't comment on the stability of trying to use these restraints in
addition to the pull code.  Seems like a lot going on at once, to me.

Also consider the fact that 5 ns is an extremely short timeframe to gather
meaningful data.  Perhaps you just need more time in each window to
equilibrate.
   At the shortest COM distance, your two molecules are still likely
experiencing
some interactions, and it may require a great deal of sampling in this
window to
converge the simulations.  You haven't shown the rest of your .mdp file
(always
a good idea!), so we can only guess at whether or not your other settings
should
lead to a sensible result.

-Justin