What should be a bond in Gromacs topologies?

[jbarnoud]

A few days ago, @mnmelo and I talked about what should be considered a bond in Gromacs topologies. Indeed, Gromacs offers various ways to define a bond, an angle, a dihedral, or an improper; some of these ways are not always used with topological bonds in mind. For instance, depending on the force field and the molecule, a constraint can be a bond or not. The same goes for virtual sites.

For now, the TPR parser has a list of interactions it considers as bonds. Constraints, for instance, are part of this list meaning that there can be to much bonds read for some molecules where constraints serve an other purpose.

With the new AtomGroup model (#363), it will be possible to attach any topological attributes to a universe. Would it make sense to have each interaction type be a separate attribute? Then lists of bonds, angles, dihedrals, and impropers could get composed based on these attributes, without information getting lost in the process.

In the long shot, this information can be used in the context of #507 to guess part of the topology based on the force field.

[richardjgowers]

Hey

Yeah definitely something we can play with now. I've been thinking how I could store information against the bonds, and have this attached to the TopologyGroups too. We could then tag harmonic and constraints differently from example.

I think ideally the extra info would be similar to Attributes again. I'm pretty busy elsewhere in the code, so if you want to tinker feel free.

On Tue, 15 Dec 2015 06:08 Jonathan Barnoud notifications@github.com wrote:

A few days ago, @mnmelo https://github.com/mnmelo and I talked about what should be considered a bond in Gromacs topologies Indeed, Gromacs offers various ways to define a bond, an angle, a dihedral, or an improper; some of these ways are not always used with topological bonds in mind For instance, depending on the force field and the molecule, a constraint can be a bond or not The same goes for virtual sites

For now, the TPR parser has a list of interactions it considers as bonds Constraints, for instance, are part of this list meaning that there can be to much bonds read for some molecules where constraints serve an other purpose

With the new AtomGroup model (#363 https://github.com/MDAnalysis/mdanalysis/issues/363), it will be possible to attach any topological attributes to a universe Would it make sense to have each interaction type be a separate attribute? Then lists of bonds, angles, dihedrals, and impropers could get composed based on these attributes, without information getting lost in the process

In the long shot, this information can be used in the context of #507 https://github.com/MDAnalysis/mdanalysis/issues/507 to guess part of the topology based on the force field

— Reply to this email directly or view it on GitHub https://github.com/MDAnalysis/mdanalysis/issues/588.

[dotsdl]

@jbarnoud I'd have a look at Bonds in core/topologyattrs.py and perhaps make a subclass for constraints. I really like the idea of individual topology formats being able to find complete representation in the Topology object that's delivered by their parser. I think this could allow us to more easily write topologies in the future, too, since the information is well-contained.

[orbeckst]

On a semantic level, I would reserve "bond" for a chemical bond and not include artificial constructs such as constraints. Bonds is also what would be used to define fragments (i.e. molecules).

[jbarnoud]

Constraints can be used to describe chemical bonds. It is often used to avoid bonds vibrations and increase the step. The issue is that the same constraints can be used for other purposes. I am playing with the new topology system to store all the interactions as topology attributes and let the user choose what interactions to keep as bonds. Le 28 déc. 2015 8:32 PM, Oliver Beckstein notifications@github.com a écrit :On a semantic level, I would reserve "bond" for a chemical bond and not include artificial constructs such as constraints. Bonds is also what would be used to define fragments (i.e. molecules).

—Reply to this email directly or view it on GitHub.

[richardjgowers]

@jbarnoud that sounds sensible. I think gromacs has two types of constraints, one of which causes nonbonded exclusions (== chemical bond) and the other doesn't. Hopefully (I've remembered that correctly and) you can distinguish these inside the TPR.

I think I'd prefer if anything that was a chemical bond appeared under the .bonds attribute, but you could tag the bonds as being constraints, so maybe ag.bonds.constraints let you choose the constraints subset.

[orbeckst]

@jbarnoud is this resolved in any way or what needs to happen to close this issue? Or can you formulate it as a proposal and we can comment/"vote" on it?

[jbarnoud]

It is not fix yet. So far I think the most sensitive things to do (and that can be done quickly) is to exclude the constraints and the tabular potentials that do not create an exclusion from the "topological bonds". This is what I understand from the gromacs manual, and from my recent experience with odd topologies.

Latter on, I still plan to give access to all the interactions, and to read a maximum of information in the topology attributes.

[orbeckst]

Just to add: it would be nice if rigid water models would show up with bonds between OW and HW1 and OW and HW2 when read from a TPR:

>>> import MDAnalysis as mda; from MDAnalysis.tests.datafiles import TPR, XTC
>>> u = mda.Universe(TPR, XTC)
>>> u.select_atoms("resname SOL").bonds
<TopologyGroup containing 0 bonds>

(Not sure if there is sufficient information to make this work automatically; there are additional distance constraints to make the water model rigid and, this being TIP4P, it even contains a dummy MW site, which should not have a bond at all.)