ResidueGroup and SegmentGroup: methods that return arrays of length len(group.atoms)

[dotsdl]

Discussion of this issue originally began on #385.

ResidueGroup and SegmentGroup are subclasses of AtomGroup, and for properties including resids, resnums, and segids they each yield arrays equal in length to themselves e.g.

ag = u.atoms
rg = ag.residues

len(ag.resids) == len(ag)
len(rg.resids) == len(rg)

For properties that refer "downward" in the hierarchy, an array of size equal to the appropriate level is returned e.g.

ag = u.atoms
rg = ag.residues
sg = ag.segments

len(rg.indices) == len(rg.atoms) 
len(sg.resids) == len(sg.residues)

I think this adds some needed consistency. However, there are other properties, such as rg.charges or rg.masses that currently return arrays of len(rg.atoms) but could make sense and be useful as len(rg), where the values for the atoms in each residue are summed.

Thoughts?

[dotsdl]

@richardjgowers wrote:

It gets more interesting/confusing (depending on your taste) if you think about positions too. Should residue.position be the center of mass? Then should be residuegroup.positions return an array of the center of mass for each residue?

This would be useful for when you want to select residues based on position, I think the current best method for that is something like:
rg = u.select_atoms('prop x etc etc').residues
So you're selecting atoms then propagating upwards.

Maybe you should be able to do:
rg = u.select_residues('prop x etc etc')

If @dotsdl 's idea isn't practical, then I'd like to remove residuegroup.charges altogether (so people have to explicitly go via (residuegroup.atoms.charges). It would be very simple for every residuegroup input and output to be len(rg). So a residuegroup is useful for querying and working at the residue level only.

[orbeckst]

So a residuegroup is useful for querying and working at the residue level only.

@richardjgowers 's principle seems to be a useful starting point for re-designing the hierarchical API.
If we can get users to understand it then we could implement @dotsdl 's elegant suggestion for summing extensive (charge, mass) per-atom properties. For other properties such as position the aggregate function is not always obvious (to me) and it might be better to use an explicit centroid() or center_of_mass() (although, from a coarse-graining perspective, resgroup.position as the center of mass is appealing...). Perhaps someone can come up with a list of suggestions of what the properties/methods would do at the higher levels?

Would the above also apply to other well-defined atom containers, such as Residue and Segment?

In terms of API philosophy, we would be moving away from "everything is an AtomGroup" to more context-dependent objects and names. The typical approach would then be to habitually use something.atoms.property if one wants to have atom-based properties, something.residues.property for residue based, and something.segments.property for segment-based.

Will such an approach be preferable to the current "one value per atom (unless it is an exceptions such as ResidueGroup.resids)"?

[hainm]

@dotsdl @orbeckst

quote from @dotsdl example

ag = u.atoms
rg = ag.residues
sg = ag.segments

this seems unnatural to me if a group of atoms has residues and segments.

What's about a Topology object having atoms, residues and segments?

for atom in top.atoms: print(atoms)
for residues in top.residues: print(residues)

If I asked an old-question, it's great if you if can point to the pass discussions.

There are several packages using my example, like mdtraj, pytraj, ParmEd, chemharp, ...

http://parmed.github.io/ParmEd/html/index.html

[dotsdl]

@hainm

Getting ResidueGroups and SegmentGroups from an AtomGroup are common conveniences. They are particularly handy when you obtained a complicated selection of atoms and want the residues or segments represented among them:

ag = u.select_atoms(selection)         # get an AtomGroup of atoms matching selection
rg = ag.residues                       # get a ResidueGroup of all residues represented by at least one atom in ag
sg = ag.segments                       # get a SegmentGroup of all segments represented by at least one atom in ag

If you want to iterate through all atoms, residues, or segments in the topology, you can already do that by accessing these groupings at the top level of the Universe:

u = mda.Universe('top.psf', 'traj.dcd')

for atom in u.atoms: print(atom)
for residue in u.residues: print(residue)
for segment in u.segments: print(segment)

Keep in mind that the object model of MDAnalysis does differ from other packages, so achieving the same thing in MDAnalysis as in something else may look different. That is not a bad thing.

[hainm]

u = mda.Universe('top.psf', 'traj.dcd')

for atom in u.atoms: print(atom)
for residue in u.residues: print(residue)
for segment in u.segments: print(segment)

yes, this is what I want to see in mdanalysis.

[dotsdl]

@hainm

Good to hear. And glad it's already there. :D

[hainm]

:D I just wish to have mda compat with PY3 (soon). I am now only using PY3, especially the 3.5 version is very promising.

Just remind @richardjgowers that the summer is gone :P.

[richardjgowers]

@orbeckst the position case isn't simple but it has some nice uses. For example something simple like area per lipid for a bilayer, there'd currently be the (ugly) step of
positions = np.array([r.center_of_mass() for r in rg])

Maybe something like rg.center_of_masses() would make it clear that it's not the same as positions

[richardjgowers]

Should be fixed in the new system.