Dear Gromacs users, I encountered two problems in using g_wham to calculate pmf curves from a pulling simulation. My system consists of a water molecule which is pulled through a bilayer (the reference group). I used gromacs 4.5.6 with the following pull code options:
pull = umbrella pull_geometry = cylinder pull_r1 = 1.0 pull_r0 = 2.0 pull_group0 = DOPC pull_start = no pull_init1 = -1.600 pull_dim = N N Y pull_vec1 = 0 0 1 pull_nstxout = 1 pull_nstfout = 1 pull_group1 = PullWater pull_rate1 = 0.004 pull_k1 = 500 (Yes the pull is very fast, as this is only a preliminary test) So only the z-Axis is of interest. My first problem is this: Running ten simulations yields pullx.xvg files looking like this: > cat pullx.xvg @ title "Pull COM" @ xaxis label "Time (ps)" @ yaxis label "Position (nm)" @TYPE xy @ view 0.15, 0.15, 0.75, 0.85 @ legend on @ legend box on @ legend loctype view @ legend 0.78, 0.8 @ legend length 2 @ s0 legend "1 cZ" @ s1 legend "1 dZ" [...] 0.0240 3.89628 -1.59424 0.0260 3.89631 -1.59401 0.0280 3.89635 -1.59387 0.0300 3.89641 -1.59381 [...] Which contain the simulation time, reference group COM and the z-Axis distance between group0 and group1. The last column is obviously what I want to calculate my pmf from. However, when running g_wham, boundaries are automatically determined to be: Determined boundaries to 3.482960 and 4.119690 So gromacs determines boundaries and builds the histogram from the second column in the pullx files, which of course leads to non-sensical profile curves. Defining boundaries by hand does not help, as the histograms contain only data from the second data column. Question 1: How do I tell g_wham to use the third data column in my pullx-files? ------------------------------------------------------------------------------------------------------ The second question concerns the same system, but this time using the pullf.xvg files to obtain the pmf from the forces. Again, I obtain sensible looking data files containing e.g. >cat pullf.xvg @ title "Pull force" @ xaxis label "Time (ps)" @ yaxis label "Force (kJ/mol/nm)" @TYPE xy 0.0000 0.245646 0.0020 0.927539 0.0040 1.65975 [...] Visualisation in vmd and measurement with g_dist show the pulled water molecule passing the bilayer, changing the delta-Z coordinate from approx. -2 to +2 nm in each of ten simulations. But when I feed the pullf and tpr files to g_wham, again false boundaries are determined: Determined boundaries to -2.250796 and -1.270690 Apparently, calculating the position of the system along the reaction coordinate from the forces did not produce a correct result, g_wham 'sees' the system moving only along about a quarter as far along the reaction coordinate, than it actually moves. Question 2: What can cause this type of behaviour, is this problem known and how do I avoid it? __________________________________________________________________________________ I would be happy about any comments, Kind Regards, Thomas Dr. Thomas Steinbrecher Institut für Physikalische Chemie, KIT Kaiserstr. 12, 76131 Karlsruhe-- gmx-users mailing list [email protected] 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 [email protected]. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
