Am 21.08.2012 18:22, schrieb gmx-users-requ...@gromacs.org:
On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesier<schl...@uni-mainz.de> 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, schriebgmx-users-requ...@gromacs.org:
>
>>On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkul<jalem...@vt.edu> wrote:
>>>
>>> >
>>> >
>>> >On 8/21/12 11:18 AM, Steven Neumann wrote:
>>>>
>>>> >>
>>>> >>On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkul<jalem...@vt.edu>
>>>> >>wrote:
>>>>>
>>>>> >>>
>>>>> >>>
>>>>> >>>
>>>>> >>>On 8/21/12 11:09 AM, Steven Neumann wrote:
>>>>>>
>>>>>> >>>>
>>>>>> >>>>
>>>>>> >>>>On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkul<jalem...@vt.edu>
>>>>>> >>>>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 listgmx-us...@gromacs.org
>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 togmx-users-requ...@gromacs.org.
>* Can't post? Readhttp://www.gromacs.org/Support/Mailing_Lists
--
gmx-users mailing list gmx-users@gromacs.org
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 gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists