[gmx-users] ACE residues at positions other than resnr=1 possibly leading to non-sequential residue numbering by editconf and trjconv

[Christopher Neale]

Dear Users:

Summary:
I believe that there is some problem in the residue numbering routine that is used by editconf and trjconv
that exists at least in 4.5.3 and 4.5.5, but not in 4.0.7, and that this is somehow related to the ACE residue.

Details:
I have a system with two peptides, water, and ions. Residues are numbered consecutively in the initial .gro file
(verified by running editconf -resnr 1). The problem is that when I use grompp to generate a .tpr file and then
use editconf -f my.tpr -o out.gro I get a .gro file with residues that are not listed consecutively (same thing if I use
trjconv to extract a .gro from a .xtc file with this .tpr).

I tested this in gromacs 4.5.5 and 4.5.3. The forcefield is either Amber99sb or OPLS-AA/L. Each peptide is capped
(starts with ACE and ends with NME).

If I order them as protein1, protein2, water, ions in the .gro and .top then the .gro that I extract from the .tpr
has residue numbers that run from 1-14 (first peptide), then again 1-15 (second peptide), and continue as
16-1413 (water), 1414-1439 (ions).

If I reorder to put ions first, then I get residue numbers that run from 16-41 (ions), 1-14 (first peptide),
1-15 (second peptide), 42-1439 (water).

In each case, however, the output from gmxdump provides molblock fields that are in the correct order.

I wonder if there could be something about the ACE residue in the force field that is causing this
strange behaviour? I've used a single capped peptide before, but never more than one and the ACE residue
has always come first in my previous work.

Here is the molblock information from gmxdump in the setup in which I put the ions before the first peptide:

topology:
   name="MD"
   #atoms               = 4652
   molblock (0):
      moltype              = 0 "NA"
      #molecules           = 15
      #atoms_mol           = 1
   molblock (1):
      moltype              = 1 "CL"
      #molecules           = 11
      #atoms_mol           = 1
   molblock (2):
      moltype              = 2 "Protein_1"
      #molecules           = 1
      #atoms_mol           = 201
   molblock (3):
      moltype              = 3 "Protein_2"
      #molecules           = 1
      #atoms_mol           = 231
   molblock (4):
      moltype              = 4 "SOL"
      #molecules           = 1398
      #atoms_mol           = 3

#############

And here is the first 30 lines of the input .gro

RUN1
4652
    1NA      NA    1   1.111   1.195   1.100
    2NA      NA    2   0.903   0.086   0.133
...
   25CL      CL   25   2.684   2.864<tel:25%C2%A0%C2%A0%202.684%C2%A0%C2%A0%202.864>   3.234
   26CL      CL   26   0.233   0.729<tel:26%C2%A0%C2%A0%200.233%C2%A0%C2%A0%200.729>   2.132
   27ACE    CH3   27   3.073   2.376   1.686
   27ACE   HH31   28   3.163   2.414   1.665
...

###############################

And here is the first 30 lines of the output .gro after grompp and editconf -f .tpr

RUN1
4652
   16NA      NA    1   1.111   1.195   1.100
   17NA      NA    2   0.903   0.086   0.133
...
   40CL      CL   25   2.684   2.864<tel:25%C2%A0%C2%A0%202.684%C2%A0%C2%A0%202.864>   3.234
   41CL      CL   26   0.233   0.729<tel:26%C2%A0%C2%A0%200.233%C2%A0%C2%A0%200.729>   2.132
    1ACE    CH3   27   3.073   2.376   1.686
    1ACE   HH31   28   3.163   2.414   1.665
...

########################

Finally, I tested this by creating an OPLS topology for either ALA or ACE-ALA. I then put two molecules
of the selected species in a box and ran grompp; editconf -f .tpr -o .gro and found that the atom numbers
started again back at zero only for the system containing ACE. This narrows it down to show that problem is
not system specific, nor does it depend on the presence of the NME residue.

Further, I tested with version 4.0.7 (pdb2gmx, grompp, editconf) and found that this problem did not exist.
Nevertheless, if I used version 4.5.5 of editconf to extract from the .tpr generated by 4.0.7 then I see improper
residue numbering.

It seems to me that there is some problem in the residue numbering routine that is used by editconf and trjconv
that exists at least in 4.5.3 and 4.5.5 but not in 4.0.7. This can be problematic for later analysis if, for example, a .ndx file is constructed from a .gro with non-sequential residue numbers. It also seems to be leading to problems with our selections for genion.

If anybody has encountered this and solved it, I'd love to hear about it.

Thank you,
Chris.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20120524/d6fdc02c/attachment.html>