On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkul <jalem...@vt.edu> wrote: > > > On 8/21/12 11:18 AM, Steven Neumann wrote: >> >> On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkul <jalem...@vt.edu> wrote: >>> >>> >>> >>> On 8/21/12 11:09 AM, Steven Neumann wrote: >>>> >>>> >>>> On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkul <jalem...@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 > > -- > ======================================== > > 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@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-requ...@gromacs.org. > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing list gmx-users@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-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
