[gmx-users] deltaG from PMF

[Steven Neumann]

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?

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

Steven

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