Identify LI as a virtual site

[shadowater]

mdtraj = 1.9.9
When loading a gro file containing LI, LI will be recognized as a virtual site. Example gro and code are provided below:

gro:

GAFF force field
 6
    1LI      LI    1   2.217   0.506   2.032
    2LI      LI    2   1.894   0.500   1.548
    3LI      LI    3   0.632   1.678   0.513
    4LI      LI    4   2.177   0.529   2.413
    5LI      LI    5   1.742   1.987   0.938
    6LI      LI    6   1.992   2.295   1.738
   2.63854   2.63854   2.63854

>>> import mdtraj as md
>>> traj = md.load("test.gro")
>>> print([traj.top.atom(i).element for i in range(traj.n_atoms)])
[(0, 'virtual_site', 'VS', 0.0, 0.0),
 (0, 'virtual_site', 'VS', 0.0, 0.0),
 (0, 'virtual_site', 'VS', 0.0, 0.0),
 (0, 'virtual_site', 'VS', 0.0, 0.0),
 (0, 'virtual_site', 'VS', 0.0, 0.0),
 (0, 'virtual_site', 'VS', 0.0, 0.0)]

seems caused by

[mattwthompson]

Hmm ... I'm not sure why the element name is stripped down from two characters to one up in L334 there. I can see why 3+ characters would need to be stripped down. It might be relevant that your residue is also named "LI", not sure.

I'm pretty sure this code is vendored from OpenMM, which may or may not have generally changed this behavior upstream in the past decade. Found it! #1767 (comment)

More generally, I'd advise against using GROMACS .gro files for storing topological (or trajectory) information. There's just not much information in them and parsing them in the different ways they might be used is hard, if possible at all.

[sukritsingh]

Since this has come up twice - would it be worth printing a warning about element assignment being based on the first letter if a .gro file is being detected as used?

[mattwthompson]

Maybe when some common footguns are hit (Lithium, Chloride, etc.); I think it'd be too noisy if a file without problems also emits the warnings

[sukritsingh]

An alternative is that we just update the documentation/docstring for md.load() to note stuff like this?

[mattwthompson]

Yep I think this behavior should be noted in the doc(string(s)) either way