[gmx-users] Re: deltaG from PMF
Steven Neumann
s.neumann08 at gmail.com
Tue Aug 21 17:57:45 CEST 2012
On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesier <schlesi at 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
Steven
>
>
> Am 21.08.2012 17:25, schrieb gmx-users-request at gromacs.org:
>
>> On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>
>>> >
>>> >
>>> >On 8/21/12 11:18 AM, Steven Neumann wrote:
>>>>
>>>> >>
>>>> >>On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkul<jalemkul at vt.edu>
>>>> >> wrote:
>>>>>
>>>>> >>>
>>>>> >>>
>>>>> >>>
>>>>> >>>On 8/21/12 11:09 AM, Steven Neumann wrote:
>>>>>>
>>>>>> >>>>
>>>>>> >>>>
>>>>>> >>>>On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkul<jalemkul at 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 at 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-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
More information about the gromacs.org_gmx-users
mailing list