GROMACS

I've attached files that should allow you to reproduce the problem.
Just put all in one dir and do:

 source source-this-rc

upon further inspection, it's a bit more complicated than I imagined.

first I have to take this back

That said, it's probably easily circumvented by modifying the residue
nrs in the final topology file working with that.
(not actually tested).

unlike what I expected the top/tpr files doesn't seem to explicitly contain residue nrs for all of the system,
so it doesn't have the information to reproduce the residue numbering in the .gro file that was used to run grompp and make the tpr file.

gmx mdrun will also have the same renumbering behaviour since it also only has a tpr to work from.

gmx dump -s ions.tpr |grep resid 
         residue (208):
            residue[0]={name="VAL", nr=17, ic=' '}
            residue[1]={name="GLU", nr=18, ic=' '}
           [..snip..]
           residue[207]={name="LYS", nr=224, ic=' '}
         residue (1):
            residue[0]={name="LIG", nr=300, ic=' '}
         residue (1):
            residue[0]={name="SOL", nr=1, ic=' '}

Here the residue nr for the protein 18-224 come from the topol.top,
the LIG residue nr 300 comes from lig_GMX_2.itp (#included in topol.top)
and the SOL resnr=1 from the amber99sb-ildn.ff/spce.itp

Since the tpr doesn't contain info about the input .gro file, for the SOL residues it obviously has to make up residue numbers for the multiple SOL residues, and simplest thing is just renumber from the last avaiable residue.
Apparently mdrun/genion gets the residue number for the first molecule from the topology and numbers on from there.
Not sure what happens with topologies with multiple protein molecules with (non)overlapping residue ranges..

I guess in theory it could use the the residue nr from the itp file (nr=300) if it is higher than the last residue but I have to admit that is a bit ugly way to do it.

I still think that conserving input residuenrs would be desirable from a user perspective but I can see it's not an easy fix

I can see a few options to make it easier to get back coordinate files with a certain residue numbering:

1. use the resnr in the itp if higher than previous (ugly)
2. include explicit residue info into the tpr file
3. add an option to genion/mdrun etc (or trjconv!) to accept a reference coordinate file from which it can take correct residue numbering assuming identical atom ordering. That would make it easy to extract pdb file with some arbitrary desirable numbering scheme after an mdrun

The 2nd option (explict resnrs in the tpr(top?) would allow general conservation of resnrs by genion/mdrun and I guess would be the least surprising behaviour for users (at least this one;-)

The 3rd option, I just realized ;-) , is already possible at least for trjconv if you give '-s' an appropriately numbered pdb or gro file !!
turns out I was doing this without realizing it ;-) https://en.wikipedia.org/wiki/Rubber_duck_debugging

Hmm, I guess this makes implementation of the second option less likely ;-(

E