[gmx-users] pull_pbcatom reporting strangely when set to 0

[chris.neale at utoronto.ca]

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.