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
========================================
--
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
