I guess is that g_wham takes the distance from tpr file which calculates the distance using grommp.
It seems that the distances calculated using g_dist and grompp are different, as discussed in this forum about 10 days ago. -Jianguo ----- Original Message ----- From: Raphael Alhadeff <raphael.alhad...@mail.huji.ac.il> To: Discussion list for GROMACS users <gmx-users@gromacs.org> Cc: Sent: Monday, 10 September 2012, 23:54 Subject: Re: [gmx-users] umbrella sampling (PMF) position discrepancy Hi Justin, thank you for you quick reply. On Mon, Sep 10, 2012 at 6:41 PM, Justin Lemkul <jalem...@vt.edu> wrote: > > > On 9/10/12 11:22 AM, Raphael Alhadeff wrote: > >> Dear Gromacs users, >> >> I've been trying to run an umbrella sampling (for the purpose of pmf) >> using >> Gromacs 4.5.5. >> My system consists of a membrane protein transporter and a Na ion passing >> through it, from the water bulk on one side of the membrane to the water >> bulk on the other side. I've used the protein as the reference group and >> the Na ion as the pull group 1 (I've attached the mdp file below). I have >> used position geometry because to my understanding it is good for the case >> where the pulled group crosses the reference's COM. I made 31 frames where >> the ion is <0.2 nm apart, from -2 nm to +2 (in respect to the COM of the >> protein). I confirmed these distances using g_traj and g_dist on the >> starting gro files. >> >> As I run the g_wham analysis, using -v, I see that the 'position' of each >> of my frames are highly different than the distances I've measured in the >> gro files, not only in value but also relative to each other. The same >> numbers appear in the pull_x or pull_f (after converting appropriately) >> files. >> So what I don't understand is how does Gromacs calculate the values for >> pull_x (and thus for the 'position' for g_wham). I was under the >> impression >> it uses COM distance between group0 and group1, but trying to compare >> using >> g_traj or g_dist proved me wrong. >> I have pasted 5 datapoints as an example below, giving the distance >> measured by g_dist and the distance that pull_x gives: >> >> time(ps) g_dist(z) pull_x(1dz) >> 0 0.894 -1.817 >> 10 0.857 -1.698 >> 20 0.897 -1.866 >> 30 0.890 -1.913 >> 40 0.966 -1.781 >> 50 0.819 -1.76 >> >> >> I've read countless of threads before posting this, and could not find any >> answer, and will be very appreciative for some light into this. >> >> > The result of g_dist is the positive root of the distance equation. It > also uses x, y, and z components of the distance, while in your case only > the z component may be relevant. The dist.xvg file(s) will have each > component listed after the total distance in subsequent columns. > > I understand, and that is what I think I posted in the example. I have the z component of g_dist and following that my umbrella simulation pull_x 1dz value (which according to the parameters I gave in the mdp, should to the best of my knowledge, give the same number) yet the numbers are quite different. > >> I should mention that I am not using positional restraints on the protein. >> I assume that since the protein is inside the membrane and the ion is much >> smaller than the protein, the PR is not required, and I wanted the protein >> to be able to adjust slightly to the movement of the ion (this is a >> transporter, not a channel). >> If this is the reason that is causing me trouble I will be happy to have a >> short explanation on why this makes the distances seem somewhat random. >> >> Lastly, I will use this opportunity on the forum to have 2 technical >> clarifications - >> -Using pull_init = 0 (or any other number) does not overrun pull_start, >> rather it adds up, correct? >> > > Correct. The output of grompp prints the distances that will be used, as > a way to check. > > > -What is the difference between the profile.xvg created (default name) and >> the pmfintegrated.xvg that is sometimes being created, and how does g_wham >> decide whether to create one or not? >> >> > Which flag is giving you pmfintegrated.xvg? That's not a default name, so > without your g_wham command line, we can't guess. > > -Justin > > That is what I don't understand, I gave no flag for this file, my command line is simply g_wham -ix pull_x.dat -it tpr.dat -v Thank you again, Raphael > Thank you very much for the help.. >> >> Raphael >> >> mdp file: >> >> title = pmf >> >> integrator = md >> dt = 0.002 >> >> nsteps = 5000000 ; 10 ns >> nstcomm = 10 >> >> ; Output parameters >> nstxout = 5000 ; every 10 ps >> nstvout = 5000 >> nstxtcout = 5000 ; every 1 ps >> nstenergy = 5000 >> >> ; Bond parameters >> constraint_algorithm = lincs >> constraints = all-bonds >> continuation = yes ; continuing from NPT >> >> ; Single-range cutoff scheme >> nstlist = 5 >> ns_type = grid >> rlist = 1.2 >> rcoulomb = 1.2 >> rvdw = 1.2 >> >> ; PME electrostatics parameters >> coulombtype = PME >> fourierspacing = 0.16 >> pme_order = 4 >> >> ; Berendsen temperature coupling is on in two groups >> Tcoupl = Nose-Hoover >> tc_grps = Protein POP Sol_Ions >> tau_t = 0.5 0.5 0.5 >> ref_t = 310 310 310 >> >> ; Pressure coupling is on >> Pcoupl = Parrinello-Rahman >> pcoupltype = semiisotropic >> tau_p = 2.0 2.0 >> compressibility = 4.5e-5 4.5e-5 >> ref_p = 1.0 1.0 >> >> ; Generate velocities is off >> gen_vel = yes >> gen_seed = -1 >> >> ; Periodic boundary conditions are on in all directions >> pbc = xyz >> >> ; Long-range dispersion correction >> DispCorr = EnerPres >> >> ; Pull code >> pull = umbrella >> pull_geometry = position >> pull_dim = N N Y >> pull_start = yes >> pull_ngroups = 1 >> pull_group0 = Protein >> pull_group1 = r_4608 >> pull_init1 = 0 0 0 >> pull_vec1 = 0 0 1 >> pull_rate1 = 0.0 >> pull_k1 = 1000 >> >> > -- > ==============================**========== > > 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<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<http://lists.gromacs.org/mailman/listinfo/gmx-users> > * Please search the archive at http://www.gromacs.org/** > Support/Mailing_Lists/Search<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<http://www.gromacs.org/Support/Mailing_Lists> > -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * 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 * 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
