On 8/21/12 12:46 PM, Steven Neumann wrote:
Thanks Thomas.
Justin, could you please comment on this?
I agree with everything Thomas has said. There's not really anything to say.
-Justin
Steven
On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesier <[email protected]> wrote:
Am 21.08.2012 18:22, schrieb [email protected]:
On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesier<[email protected]>
wrote:
Since your simulations of the individual windows are about 50 ns, you
could
first dismiss the first 10 ns for equilibration, and then perform the
WHAM
analysis for 10-30 ns and 30-50 ns. If everything is fine, you should
see no
drift.
If you want to have more data for the analysis you could also use 5ns ;
5-27.5ns and 27.5-50ns.
From the PMF it seems that the equilibrium state should be around 0.6
nm. To
be sure, you can perform a normal simulation (without any restraints)
from
you initial starting window (~0.4nm) and a window near the minima
(~0.6nm).
Then after the equilibration phase, look at the distribution of the
distance
along the reaction coordinate. If in both cases the maximum is at
~0.6nm,
this should be the 'true' equilibrium state of the system (instead of
the
first window of the PMF calculation) and i would measure \Delta G from
this
point.
Greetings
Thomas
Thanks Thomas for this but finally I realised that my first
configuration corresponds to 0.6 nm which is the minima so I take the
free energy difference based on this value and plateau.
I want also to calculate error bars. Would you do this:
Final PMF curve for 10-50 ns
Error bars from:
g_wham -b 30000 -e 40000
g_wham -b 50000 -e 60000
Think this approach would be good to see if you have any drifts.
But for error bars there is something implemented in 'g_wham'. But i never
used it, since for my system umbrella sampling is not really applicable,
only TI. So i can't comment on it, if there is anything one should be aware
of, or similar. But 'g_wham -h' prints some info about how to use the error
analysis
Greetings
Thomas
Steven
Am 21.08.2012 17:25, [email protected]:
On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkul<[email protected]>
wrote:
On 8/21/12 11:18 AM, Steven Neumann wrote:
On Tue, Aug 21, 2012 at 4:13 PM, Justin
Lemkul<[email protected]>
wrote:
On 8/21/12 11:09 AM, Steven Neumann wrote:
On Tue, Aug 21, 2012 at 3:48 PM, Justin
Lemkul<[email protected]>
wrote:
On 8/21/12 10:42 AM, Steven Neumann wrote:
Please see the example of the plot from US
simulations and
WHAM:
http://speedy.sh/Ecr3A/PMF.JPG
First grompp of frame 0 corresponds to 0.8 nm
- this is what
was shown
by grompp at the end.
The mdp file:
; Run parameters
define = -DPOSRES_T
integrator = md ; leap-frog
integrator
nsteps = 25000000 ; 100ns
dt = 0.002 ; 2 fs
tinit = 0
nstcomm = 10
; Output control
nstxout = 0 ; save coordinates
every 100 ps
nstvout = 0 ; save velocities every
nstxtcout = 50000 ; every 10 ps
nstenergy = 1000
; Bond parameters
continuation = no ; first
dynamics run
constraint_algorithm = lincs ; holonomic
constraints
constraints = all-bonds ; all bonds
(even heavy atom-H
bonds)
constrained
; Neighborsearching
ns_type = grid ; search neighboring
grid cells
nstlist = 5 ; 10 fs
vdwtype = Switch
rvdw-switch = 1.0
rlist = 1.4 ; short-range
neighborlist cutoff (in
nm)
rcoulomb = 1.4 ; short-range
electrostatic cutoff (in
nm)
rvdw = 1.2 ; short-range van der
Waals cutoff (in
nm)
ewald_rtol = 1e-5 ; relative strenght
of the
Ewald-shifted
potential rcoulomb
; Electrostatics
coulombtype = PME ; Particle Mesh
Ewald for
long-range
electrostatics
pme_order = 4 ; cubic
interpolation
fourierspacing = 0.12 ; grid spacing
for FFT
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
optimize_fft = yes
; Temperature coupling is on
tcoupl = V-rescale ;
modified
Berendsen
thermostat
tc_grps = Protein LIG_Water_and_ions ;
two coupling
groups -
more
accurate
tau_t = 0.1 0.1 ;
time constant,
in ps
ref_t = 318 318 ;
reference
temperature,
one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ;
pressure
coupling is on
for
NPT
pcoupltype = isotropic ;
uniform scaling
of box
vectors
tau_p = 2.0 ;
time constant,
in ps
ref_p = 1.0 ;
reference
pressure, in
bar
compressibility = 4.5e-5 ;
isothermal
compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off
vdW scheme
; Velocity generation
gen_vel = yes ; assign velocities
from Maxwell
distribution
gen_temp = 318 ; temperature for
Maxwell distribution
gen_seed = -1 ; generate a random
seed
; These options remove COM motion of the
system
; Pull code
pull = umbrella
pull_geometry = distance
pull_dim = N N Y
pull_start = yes
pull_ngroups = 1
pull_group0 = Protein
pull_group1 = LIG
pull_init1 = 0
pull_rate1 = 0.0
pull_k1 = 500 ; kJ mol^-1
nm^-2
pull_nstxout = 1000 ; every 2 ps
pull_nstfout = 1000 ; every 2 ps
Based on these settings you're allowing grompp to
set the
reference
distance
to whatever it finds in the coordinate file. It
seems clear to
me that
the
sampling indicates what I said before - you have
an energy
minimum
somewhere
other than where you "started" with. What that
state
corresponds to
relative to what you think is going on is for you
to decide
based on
the
nature of your system. Whatever is occurring at
0.6 nm of COM
separation
is
of particular interest, since the energy minimum
is so distinct.
So based on this the deltaG will correspond to -5.22
as the
initial
state was taken at 0.4 nm corresponding to 0 kcal/mol
as the
moment
corresponding to the minimum is the coordinate from
SMD where last
hydrogen bond was broken. Would you agree?
Based on the very little information I have, no. It
would appear
that
the
0.4 nm separation is in fact some metastable state and
the true
energy
minimum is at 0.6 nm of COM separation. What's going on
at that
location?
My mistake. The initial grompp of 1st configuartion (where
ligand
stakced on keratin surface) corresponds to 0.6 nm where
is the minimum. Thus deltaG would be -7.22 kcal/mol. Am I
right? Or
Shall I take difference between 0 and 5.22 ?
-7.22 kcal/mol sounds much more logical to me. If your first
configuration
is at the energy minimum, that's not something you ignore. The
zero
point
can be set wherever you like with the g_wham flag -zprof0, so
it's
really
rather arbitrary. The WHAM algorithm simply sets the leftmost
window
(smallest value along the reaction coordinate) to zero to
construct the
remainder of the PMF curve.
I hope you're doing a thorough analysis of
convergence if you're
generating
velocities at the outset of each run, and
removing
unequilibrated
frames
>from your analysis.
When I use WHAM I skip first 200 ps of each window as
the
equilibration.
That seems fairly short, especially given the generation
of
velocities in
conjunction with the Parrinello-Rahman barostat, which
can be very
temperamental.
Would you suggest e.g. skip 1 ns?
I'm not going to make an arbitrary guess. It's up to you to
analyze the
timeframe required for whatever relevant observables to converge.
-Justin
Thanks for this.
Steven
--
--
gmx-users mailing [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the www
interface or send it [email protected].
* Can't post? Readhttp://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the www
interface or send it to [email protected].
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
========================================
Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
--
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Only plain text messages are allowed!
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
