[gmx-users] deltaG from PMF

Steven Neumann s.neumann08 at gmail.com
Tue Aug 21 17:24:13 CEST 2012


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

>
> --
> ========================================
>
> Justin A. Lemkul, Ph.D.
> Research Scientist
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> ========================================
> --
> 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