[gmx-users] Pulling and g_wham - Two analysis problems

[Justin Lemkul]

On 2/8/13 6:00 AM, Steinbrecher, Thomas (IBG) wrote:
> 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,

Both of the issues sound like bugs.  Please file a report on redmine.gromacs.org 
with all files necessary to reproduce the problem.

-Justin

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

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

========================================