Gromacs water defined with SETTLES appear as disconnected when read from TPR

[jbarnoud]

Expected behaviour

Water molecules appear as fragments.

Actual behaviour

Water defined in Gromacs to use the SETTLES algorithm appear as disconnected when read from a TPR as the SETTLES definition is not read as bonds. As a consequence, the atoms appear disconnected.

I assume that the same system appears with water connected when read from any other topology source that will guess the bonds.

@davidercruz noticed the issue, but I think I remember @orbeckst mentioning it some time ago.

Code to reproduce the behaviour

import MDAnalysis as mda
from MDAnalysisTests.datafiles import TPR, XTC
u = mda.Universe(TPR, XTC)
water = u.select_atoms('resname SOL')
print('Number of water atoms', len(water))
print('Number of water molecules', len(water.split('molecule')))
print('Number of water fragments', len(water.fragments))
print('Number of bonds in the water', len(water.bonds))

Currently version of MDAnalysis:

(run python -c "import MDAnalysis as mda; print(mda.__version__)")

[richardjgowers]

related #588?

[jbarnoud]

It is related indeed. Though, I think there is no reason not to count SETTLES constraints as bonds while there may be reasons not to count some bonds or constraints that do not represent covalent bonds.

[richardjgowers]

@jbarnoud I'm not sure how easy it'll be to separate true chemical bonds and fudgy geometric constraints that keep the molecule sane. I'd say we can just include all constraints and the user can pick them apart themselves

[mnmelo]

In #588 @orbeckst makes the case that dummy atoms (such as in TIP4P) should not have bonds. OTOH, I recently discussed this with @davidercruz in the context of a make_whole transformation. make_whole walks over bonds, and fails for TIP4P, which is unintuitive.

[richardjgowers]

So for things like TIP4P we can probably recognise and treat it correctly (ie only pick "real" bonds), but the amount of things we can special case like this is obviously finite, which would then lead to some things getting special treatment and others getting ignored/failing....? Is the fake particle included in a tip4p trajectory? I'm guessing if it is it's useful for analysis too?

[jbarnoud]

In the case of TIP4P, the dummy particle is not affected by SETTLES. It is a virtual site that is a slightly different problem; which I have with other molecules such as Martini cholesterol.

[richardjgowers]

@jbarnoud you can open an issue for virtual site handling if it's something that can be improved