[gmx-users] deltaG from PMF

[Justin Lemkul]

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?

>> 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.

-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

========================================