[gmx-users] deltaG from PMF

Justin Lemkul jalemkul at vt.edu
Tue Aug 21 18:48:45 CEST 2012



On 8/21/12 12:46 PM, Steven Neumann wrote:
> Thanks Thomas.
> Justin, could you please comment on this?
>

I agree with everything Thomas has said.  There's not really anything to say.

-Justin

> Steven
>
> On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesier <schlesi at uni-mainz.de> wrote:
>> Am 21.08.2012 18:22, schrieb gmx-users-request at gromacs.org:
>>>
>>> On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesier<schlesi at uni-mainz.de>
>>> wrote:
>>>>
>>>>> Since your simulations of the individual windows are about 50 ns, you
>>>>> could
>>>>> first dismiss the first 10 ns for equilibration, and then perform the
>>>>> WHAM
>>>>> analysis for 10-30 ns and 30-50 ns. If everything is fine, you should
>>>>> see no
>>>>> drift.
>>>>> If you want to have more data for the analysis you could also use 5ns ;
>>>>> 5-27.5ns and 27.5-50ns.
>>>>>
>>>>>  From the PMF it seems that the equilibrium state should be around 0.6
>>>>> nm. To
>>>>> be sure, you can perform a normal simulation (without any restraints)
>>>>> from
>>>>> you initial starting window (~0.4nm) and a window near the minima
>>>>> (~0.6nm).
>>>>> Then after the equilibration phase, look at the distribution of the
>>>>> distance
>>>>> along the reaction coordinate. If in both cases the maximum is at
>>>>> ~0.6nm,
>>>>> this should be the 'true' equilibrium state of the system (instead of
>>>>> the
>>>>> first window of the PMF calculation) and i would measure \Delta G from
>>>>> this
>>>>> point.
>>>>>
>>>>> Greetings
>>>>> Thomas
>>>
>>>
>>>
>>> Thanks Thomas for this but finally I realised that my first
>>> configuration corresponds to 0.6 nm which is the minima so I take the
>>> free energy difference based on this value and plateau.
>>>
>>> I want also to calculate error bars. Would you do this:
>>>
>>> Final PMF curve for 10-50 ns
>>>
>>> Error bars from:
>>> g_wham -b 30000 -e 40000
>>>
>>> g_wham -b 50000 -e 60000
>>>
>>
>> Think this approach would be good to see if you have any drifts.
>> But for error bars there is something implemented in 'g_wham'. But i never
>> used it, since for my system umbrella sampling is not really applicable,
>> only TI. So i can't comment on it, if there is anything one should be aware
>> of, or similar. But 'g_wham -h' prints some info about how to use the error
>> analysis
>> Greetings
>> Thomas
>>
>>
>>
>>>
>>>
>>> Steven
>>>
>>>
>>>>>
>>>>>
>>>>> Am 21.08.2012 17:25, schriebgmx-users-request at gromacs.org:
>>>>>
>>>>>
>>>>>>> On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkul<jalemkul at vt.edu>
>>>>>>> wrote:
>>>>>>
>>>>>>>>>
>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> 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
>>>>>
>>>>>>>
>>>>>>> Thanks for this.
>>>>>>>
>>>>>>> Steven
>>>>>>>
>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> --
>>>>
>>>>>
>>>>>
>>>>> --
>>>>> gmx-users mailing listgmx-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 togmx-users-request at gromacs.org.
>>>>> * Can't post? Readhttp://www.gromacs.org/Support/Mailing_Lists
>>
>>
>> --
>> 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

-- 
========================================

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