On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesier <schl...@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 Steven > > > Am 21.08.2012 17:25, schrieb gmx-users-requ...@gromacs.org: > >> 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 >> >>> > >>> >-- > > > -- > 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
