[gmx-users] Re: deltaG from PMF
Thomas Schlesier
schlesi at uni-mainz.de
Tue Aug 21 17:49:51 CEST 2012
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-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
>
>> >
>> >--
More information about the gromacs.org_gmx-users
mailing list