[gmx-users] deltaG from PMF

[Steven Neumann]

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


I will appreciate your help,
Steven



On Tue, Aug 21, 2012 at 3:27 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> On 8/21/12 10:23 AM, Steven Neumann wrote:
>>
>> On Tue, Aug 21, 2012 at 3:07 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>
>>>
>>>
>>> On 8/21/12 9:36 AM, Steven Neumann wrote:
>>>>
>>>>
>>>> Dear Gmx Users,
>>>>
>>>> I got confused about reading free energy difference from PMF curve.
>>>> It is the difference between maximum value on PMF curve (plateau -
>>>> final state) and the starting point corresponding to minima.
>>>>
>>>> So e.g. my curve starts at 0 [kcal/mol] going staright to minima of -2
>>>> kcal/mol and then increase obtaining plateau of 14 kcal/mol. DeltaG =
>>>> -16 kcal/mol or -12 kcal/mol ?
>>>>
>>>
>>> It seems that your system has discovered a more energetically favorable
>>> location where the minimum is at -2 kcal/mol.  It should not be assumed
>>> that
>>> whatever arbitrary "start" point we choose is the energy minimum.  In
>>> your
>>> case, it is not.  DeltaG is the difference between final and initial
>>> states.
>>> Only you know what these states are.
>>>
>>> -Justin
>>
>>
>> Thank you. Ok - mu curve starts at 0 kcal/mol corresponding to the
>> reaction coordinate of 0.5 nm. going to -2 kcal mol at 0.6 nm then
>> increasing as I said before to plateau of 14 kcal/mol.
>>
>> My first window (shown by grompp) corresponds to 0.8 nm of the
>> reaction coordinate... So the initial state should be at 0.8 nm which
>> corresponds to 0.5 kcal/mol or always take the initial state as
>> initial of the curve?
>>
>
> That depends on how you've set up your .mdp files.  If you're telling a
> configuration in which the groups are separated by 0.8 nm that they should
> be separated by 0.5 nm instead, you need to allow for some time for that to
> adjust.  The restraint distances are what you tell them to be.  Your initial
> configurations should be as close to that COM separation as possible.
>
> So, the answer is, I don't have an answer, because you didn't provide an
> .mdp file.
>
> -Justin
>
>
> --
> ========================================
>
> 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