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


Am 21.08.2012 17:25, schrieb gmx-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 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

Reply via email to