[gmx-users] deltaG from PMF
Thomas Schlesier
schlesi at uni-mainz.de
Tue Aug 21 18:37:54 CEST 2012
Am 21.08.2012 18:22, schrieb gmx-users-request at gromacs.org:
> 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
>
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-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 listgmx-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 togmx-users-request at gromacs.org.
>> >* Can't post? Readhttp://www.gromacs.org/Support/Mailing_Lists
More information about the gromacs.org_gmx-users
mailing list