[gmx-users] RE: Umbrella sampling question

[Christopher Neale]

Why not do two umbrella sampling simulations: one with initial conformations from your faster pulling and one with initial conformations from your slower pulling. Then you can run them both as regular US simulations until (a) neither US PMF is drifting systematically with increasing simulation time and (b) they converge to the same answer. This way, you might also be able to make some comment on the relative importance of doing the initial pulling slowly.

Next time, please wait until your post has appeared on the list (it is currently awaiting approval for inclusion of large images) before sending me a private email asking me to comment on-list (which is, by itself, fine with me).

Also: I can't make heads or tails of your plots. (i) I have no idea which one corresponds to which force constant; (ii) you mentioned 4 force constants but you posted 8 figures; (iii) I can not read any of the axes or labels.

Chris.
________________________________________
From: surampudi3 at gmail.com [surampudi3 at gmail.com]
Sent: 29 July 2013 19:39
To: chris.neale at utoronto.ca
Subject: Umbrella sampling question

Hi Chris,

I have seen many of your contributions to the Umbrella sampling topic in the forum. It would be of immense help to me if you can throw light on a few points.

I am running my first umbrella sampling simulation, and having hard time to decide if the PMF plot I am getting is right or wrong. I am pulling one atom from a short polymer chain from the center of a micelle and trying to obtain a PMF plot along the distance from the center of the micelle to the surface and eventually few nanometers away from the micelle. Here is the mdp file for the pulling:
title       = Umbrella pulling simulation
define      = -DPOSRES_surfpoly
; Run parameters
integrator  = md
dt          = 0.002     ; timestep (ps)
tinit       = 0
nsteps      = 250000    ;  500 ps= 0.5ns
nstcomm     = 10        ; remove COM every 10 steps
; Output parameters
nstxout     = 5000
nstvout     = 5000
nstfout     = 250
nstxtcout   = 250
nstenergy   = 500
; Bond parameters
constraint_algorithm    = lincs
constraints             = all-bonds
continuation            = yes       ; continuing from NPT
; 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     = System
tau_t       = 0.5
ref_t       = 300
; Pressure coupling is on
Pcoupl          = Parrinello-Rahman
pcoupltype      = isotropic
tau_p           = 1.0
compressibility = 4.5e-5
ref_p           = 1.0
refcoord_scaling = com
; 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  ; COM pulling using umbrella potential between reference group and other groups
;pull_geometry   = distance  ; simple distance increase
pull_geometry   = direction  ; simple distance increase
;pull_dim        = Y N N     ; pulling in x-direction
pull_vec1       = 1 0 0     ; pulling in x-direction
pull_start      = yes       ; define initial COM distance > 0
pull_ngroups    = 1         ; we are only applying a pulling force to one group
pull_group0     = surf       ; reference group for pulling
pull_group1     = poly    ; group to which pulling force is applied
pull_rate1      = 0.01      ; 0.01 nm per ps = 10 nm per ns
pull_k1         = 300       ; kJ mol^-1 nm^-2

I have repeated pull run using various force constants (pull_k1=300,500,1000,2000 kJ/mol/nm^2) to pull over 5 nm distance with 0.1nm windows and none of them is a clear pick as in the tutorial as to which one to use for umbrella sampling. From my understanding, force constant should neither be too low or too high. Is there a simpler way to decide which one is the right force constant to use from the plots? Please see the link below for the F vs time and  position vs time plots (top to bottom order --- pull_k1=2000, 1000, 500, 300 kJ/mol/nm^2).

http://gromacs.5086.x6.nabble.com/Umbrella-sampling-force-vs-time-plots-tp5009709p5010097.html

Thanks a lot,
Andy S.

_____________________________________
Sent from http://gromacs.5086.x6.nabble.com