u.select_atoms(...).bonds should return only bonds with atoms *inside* of the selection

[cbouy]

No necessarily a bug, but I think it's worth mentioning:
When making a selection, should the resulting bonds only contain bonds where both atoms are in the selection, or should it return all bonds where at least one atom of the atomgroup is part of the bond ?

The current behavior is to return all bonds where at least one atom of the atomgroup is part of the bond

Code to reproduce the behavior

>>> import MDAnalysis as mda
>>> from rdkit import Chem

>>> mol = Chem.MolFromSequence("GAG")
>>> u = mda.Universe(mol)
>>> u.select_atoms("resid 2").indices
array([4, 5, 6, 7, 8])

>>> u.select_atoms("resid 2").bonds.indices
array([[2, 4],
       [4, 5],
       [5, 6],
       [5, 8],
       [6, 7],
       [6, 9]], dtype=int32)

Should it return only the green bonds, or should it return the green bonds and the yellow ones (between atoms 2-4 and 6-9) as well ?
The current behavior is green+yellow.

Current version of MDAnalysis

• Which version are you using? (run python -c "import MDAnalysis as mda; print(mda.__version__)")
  2.0.0-dev0
• Which version of Python (python -V)?
  Python 3.6.10
• Which operating system?
  Ubuntu 18.04.4 LTS through a VM (WSL2)

[lilyminium]

Try atomgroup_intersection(strict=True) -- not sure if there's a quicker way to do this. Maybe turning bonds etc into functions (ag.bonds(strict=True)) instead of properties would be easier.

ag = u.select_atoms('protein')
bonds = ag.bonds.atomgroup_intersection(ag, strict=True)

[cbouy]

Great, thank you! Yes it looks a bit heavy to be honest, but I wouldn't expect users to use this everyday so...
I'm not sure if changing to ag.bonds() is a good idea, tbh as a new user I wouldn't expect bonds to be a method, it's a bit confusing. Maybe adding a get_bonds() with parameters could work ?

[lilyminium]

I think this is what .bonds does under the hood anyway, but strict is False by default. I personally would expect .bonds to only be the ones inside the atomgroup, but it’s probably too late to change the API and not irritate people... either way an easier to find method would definitely be helpful!

…

On 3 Jul 2020, at 22:43, Cédric Bouysset ***@***.***> wrote: Great, thank you! Yes it looks a bit heavy to be honest, but I wouldn't expect users to use this everyday so... I'm not sure if changing to ag.bonds() is a good idea, tbh as a new user I wouldn't expect bonds to be a method, it's a bit confusing. Maybe adding a get_bonds() with parameters could work ? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[richardjgowers]

This is probably an API break we should make tbh

…

On Fri, Jul 3, 2020 at 14:07, Lily Wang ***@***.***> wrote: I think this is what .bonds does under the hood anyway, but strict is False by default. I personally would expect .bonds to only be the ones inside the atomgroup, but it’s probably too late to change the API and not irritate people... either way an easier to find method would definitely be helpful! > On 3 Jul 2020, at 22:43, Cédric Bouysset ***@***.***> wrote: > > Great, thank you! Yes it looks a bit heavy to be honest, but I wouldn't expect users to use this everyday so... > I'm not sure if changing to ag.bonds() is a good idea, tbh as a new user I wouldn't expect bonds to be a method, it's a bit confusing. Maybe adding a get_bonds() with parameters could work ? > > — > You are receiving this because you commented. > Reply to this email directly, view it on GitHub, or unsubscribe. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#2821 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACGSGB7XLJVWEECB2YPEXG3RZXJYXANCNFSM4OPVGHDA> .

[orbeckst]

I agree:

personally would expect .bonds to only be the ones inside the atomgroup,

Well, you can go wild for 2.0... just raise an issue to change it. Probably applies to angles, dihedrals, etc as well.

Having an way to get bonded information with the outside world when you need it would be good. If it's an expert-level thing then maybe we don't need getters for everything, just a generic ag.get_bonded_type("bonds", outside=True) or similar might do it. You can then make it work for "atoms", "angles", "dihedrals", .... It would return different types, depending on the argument, which is perhaps not desirable, but still better than a proliferation of ag.get_X() that confuses users when they also have ag.X.

[orbeckst]

just raise an issue to change it

... or we can just decide that this is the issue where we decide to change it ;-). Then it should be consistent. though, and apply to all bonded attributes.

[orbeckst]

I don't see any opposing views to the change so I am removing the proposal tag and rename the issue appropriately.