[gmx-users] deltaG from PMF
Justin Lemkul
jalemkul at vt.edu
Tue Aug 21 17:21:51 CEST 2012
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
--
========================================
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
========================================
More information about the gromacs.org_gmx-users
mailing list