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

[Steinbrecher, Thomas (IBG)]

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