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

[chris.neale at utoronto.ca]

You're right, although I did not realize this at first. I got confused  
(again) by the global numbering. Regarding the off-by-one numbering, I  
can figure out what is meant there and that was not really my question.

Thanks Berk, and sorry for the confusion.
Chris.


> I coded this and I looked at the code again.
> pbcatom is always the global, system wide atom number.
> But it seems that mdrun forgets to add one to the atom number
> or maybe this is part of the mdrun output scheme of starting numbering at 0.
> We would have to check this.
>
> As far as I can see that is consistent with your results, right?
>
> Berk

-- previous message --

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