[gmx-users] deltaG from PMF
Steven Neumann
s.neumann08 at gmail.com
Tue Aug 21 17:18:40 CEST 2012
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 ?
>
>
>>> 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?
>
>
> -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
More information about the gromacs.org_gmx-users
mailing list