Issues with trjorder using TIP5P water model

User discussions
1

GROMACS version: 2020.4
GROMACS modification: No

Hi all,

I’m encountering an issue with trjorder when trying to order TIP5P water molecules based on their distance from a polymer.

I’ve never had any issue when using trjorder to order SPC/E water, here’s an example of the ordered PDB file:

ATOM  49740  OW  SOL  6484      66.620  68.720   5.780  0.00  0.30           O
ATOM  49741  HW1 SOL  6484      67.530  69.050   6.030  0.00  0.30           H
ATOM  49742  HW2 SOL  6484      65.930  69.350   6.110  0.00  0.30           H
ATOM   8022  OW  SOL  2578      63.370  68.640  13.180  0.00  0.31           O
ATOM   8023  HW1 SOL  2578      64.240  69.140  13.270  0.00  0.31           H
ATOM   8024  HW2 SOL  2578      63.490  67.710  13.520  0.00  0.31           H

However, when using the TIP5P water model, I suspect the two dummy atoms for the lone electron pairs (LP1 and LP2) are causing issue:

ATOM   8550  HW1 SOL  1659      80.520  69.850  40.270  0.00  0.25           H
ATOM   8551  HW2 SOL  1659      80.830  69.080  39.000  0.00  0.25           H   
ATOM   8552  LP1 SOL  1659      80.870  70.440  39.090  0.00  0.25
ATOM   8907  LP1 SOL  1730      82.400  68.810  44.120  0.00  0.25
ATOM   8908  LP2 SOL  1730      82.410  68.500  45.220  0.00  0.25
ATOM   8909  OW  SOL  1731      58.440  65.660  29.770  0.00  0.25           O
ATOM  14493  LP2 SOL  2847      74.260  53.470  22.730  0.00  0.26
ATOM  14494  OW  SOL  2848      77.200  45.790   7.490  0.00  0.26           O
ATOM  14495  HW1 SOL  2848      77.850  45.550   6.840  0.00  0.26           H

In both cases, I’ve specified that trjorder used the COM as the coordinate to calculate the distance. So one would expect that you’d see five consecutive entries, e.g. OW, HW1, HW2, LP1, LP2 all from the same molecule. However, as you can see above, the first three entries are from molecule 1659, then two entries from molecule 1730, etc. It seems that trjorder is writing out groups of three atoms with consecutive atomic indices (i.e. 8907, 8908, 8909) which would work for a three-point model such as SPC/E but not a five-point model such as TIP5P. I’m aware that you can specify the number of atoms per molecule using the -na flag, however specifying -na 5 for TIP5P yeilds an error that the number of atoms in SOL isn’t a multiple of 5 (the number specified is a multiple of 3 however).

Any suggestions? Perhaps this is a bug?

2

This definitely looks buggy to me, though I can’t point to a specific piece of code yet that would trigger this. Please report this to the GROMACS GitLab site.