Date: Sun, 20 Jun 2010 18:54:13 -0400
From: chris.neale at utoronto.ca
To: gmx-users at gromacs.org
Subject: [gmx-users] pull_pbcatom reporting strangely when set to 0
Dear gromacs users:
I am confused (again) about pull_pbcatom0. From this analysis, I am
questioning if the pbcatom reported from grompp and in the top of
the .log file is numerically defined *within* a group or in the
context of the entire .gro file.
Either way, I think that I have conclusively shown below that
setting pull_pbcatomN to 0 is not reported properly. My hunch is
that it is just a reporting problem and is working alright. For the
reporting by grompp and at the beginning of mdrun, I think that the
pbcatom is defined *within* a group, but that the selection of the
pbcatom when setting pbcatom=0 is incorrect (sometimes drastically).
I'm posting to the mailing list in the hopes that somebody can
simply take a look at the relevant code and see where the problem is.
Here, I have a .gro file that has protein, then lipid, then water. I
have the lipid as pull group 0 and the protein as pull group 1.
The protein has 14 atoms and the lipid group has 3456 atoms.
These tests have all been run with 4.0.7 (and also 4.0.5 gives
identical results).
If I specify pull_pbcatom0=10 in my .mdp, then:
grompp says:
Pull group natoms pbc atom distance at start reference at t=0
0 3456 10
1 14 7 -0.000 0.000 -0.001 0.000 0.000 0.000
and the .log file says:
pull_group 0:
atom (3456):
atom[0,...,3455] = {14,...,3469}
weight: not available
pbcatom = 9
This makes me think that at this stage the pull_pbcatom0 is defined
within the group such that the accessible numebrs are [0,...,3455]
-- If this was not the case, then pbcatom0 should be reported as
9+14 protein atoms = 25.
####################################################################################
Ok, so now I define pull_pbcatom0=0 in my .mdp, then:
grompp says:
Pull group natoms pbc atom distance at start reference at t=0
0 3456 1742
1 14 7 -0.000 0.000 -0.001 0.000 0.000 0.000
and the .log file says:
pull_group 0:
atom (3456):
atom[0,...,3455] = {14,...,3469}
weight: not available
pbcatom = 1741
Which surprised me because 3456/2=1728 and 1728 -1 (use zero as
first #) +14 (protein atoms) = 1741
However, this suggests that at this stage the pull_pbcatom0 is
defined in the context of the entire .gro file, not only within the
group. This contradicts what I saw above.
####################################################################################
Next, I added a second protein (there are now 28 atoms of protein
before the lipid coordinates start) and left pull_pbcatom0=0, now:
grompp says:
Pull group natoms pbc atom distance at start reference at t=0
0 3456 1756
1 28 14 -0.000 0.000 -0.001 0.000 0.000 0.000
and the .log file says:
pull_group 0:
atom (3456):
atom[0,...,3455] = {28,...,3483}
weight: not available
pbcatom = 1755
Where 1755 = 1741+14 so again it looks as if the pull_pbcatom0 is
defined (or at least decided) in the context of the entire .gro file.
####################################################################################
####################################################################################
####################################################################################
Now I reverse the order of the protein and the lipid, with only one
copy of the protein again.
Here, I specify pull_pbcatom0=10 (pull group 0 remains the lipid):
grompp says:
Pull group natoms pbc atom distance at start reference at t=0
0 3456 10
1 14 3463 -0.000 0.000 -0.001 0.000 0.000 0.000
and the .log file says:
pull_group 0:
atom (3456):
atom[0,...,3455] = {0,...,3455}
weight: not available
pbcatom = 9
####################################################################################
Now I specify pull_pbcatom0=0
grompp says:
Pull group natoms pbc atom distance at start reference at t=0
0 3456 1728
1 14 3463 -0.000 0.000 -0.001 0.000 0.000 0.000
and the .log file says:
pull_group 0:
atom (3456):
atom[0,...,3455] = {0,...,3455}
weight: not available
pbcatom = 1727
...
pull_group 1:
atom (14):
atom[0,...,13] = {3456,...,3469}
weight: not available
pbcatom = 3462
And from these last 2 I conclude that pull_pbcatom0 had better be
defined in the context of the entire file, or else the 14 atom
peptide having a pull_pbcatom of 3463 is going to be a real problem.
#####################################################################################
To test that, I removed all of the water so that the protein pbcatom
would go out of bounds if it was really #3462 as defined within the
group:
gromp says (the exact same as before):
Pull group natoms pbc atom distance at start reference at t=0
0 3456 1728
1 14 3463 -0.000 0.000 -0.001 0.000 0.000 0.000
and the .log file says:
pull_group 0:
atom (3456):
atom[0,...,3455] = {0,...,3455}
weight: not available
pbcatom = 1727
pull_group 1:
atom (14):
atom[0,...,13] = {3456,...,3469}
weight: not available
pbcatom = 3462
And it appears to be stable over 1000 steps of MD, so perhaps it's
not going out of bounds.
Thank you very much,
Chris.
