Dihedrals and selections

[mimischi]

When creating a selection, e.g., backbone = u.select_atoms('backbone') and then looking at all the dihedral angles of that selection, MDAnalysis also returns dihedrals between atoms pairs that should not be presented.

Initial code

import numpy as np
import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD

u = mda.Universe(PSF, DCD)
backbone = u.select_atoms('backbone')

backbone
<AtomGroup with 855 atoms>

Expected behaviour

Only list dihedral angles between atoms in the selection (here: backbone). Total number of dihedrals: 425.

# Extract indiced of dihedral angles phi and psi 
dihedral_indices = np.zeros([2 * (len(backbone.residues) - 1)], dtype=np.int)

j = 0
for i, dih in enumerate(backbone.dihedrals):
    if dih.atoms[0].name == 'C' and dih.atoms[-1].name == 'C' or dih.atoms[0].name == 'N' and dih.atoms[-1].name == 'N': 
        dihedral_indices[j] = i
        j += 1

for dih in backbone.dihedrals[dihedral_indices]:
    print dih.atoms

<AtomGroup [<Atom 1: N of type 56 of resname MET, resid 1 and segid 4AKE>, <Atom 5: CA of type 22 of resname MET, resid 1 and segid 4AKE>, <Atom 18: C of type 20 of resname MET, resid 1 and segid 4AKE>, <Atom 20: N of type 54 of resname ARG, resid 2 and segid 4AKE>]>
<AtomGroup [<Atom 18: C of type 20 of resname MET, resid 1 and segid 4AKE>, <Atom 20: N of type 54 of resname ARG, resid 2 and segid 4AKE>, <Atom 22: CA of type 22 of resname ARG, resid 2 and segid 4AKE>, <Atom 42: C of type 20 of resname ARG, resid 2 and segid 4AKE>]>
<AtomGroup [<Atom 20: N of type 54 of resname ARG, resid 2 and segid 4AKE>, <Atom 22: CA of type 22 of resname ARG, resid 2 and segid 4AKE>, <Atom 42: C of type 20 of resname ARG, resid 2 and segid 4AKE>, <Atom 44: N of type 54 of resname ILE, resid 3 and segid 4AKE>]>
[... truncated ...]
<AtomGroup [<Atom 3332: C of type 20 of resname LEU, resid 213 and segid 4AKE>, <Atom 3334: N of type 54 of resname GLY, resid 214 and segid 4AKE>, <Atom 3336: CA of type 23 of resname GLY, resid 214 and segid 4AKE>, <Atom 3339: C of type 32 of resname GLY, resid 214 and segid 4AKE>]>

Actual behaviour

Lists all existing dihedrals that start/end in the selection. Total number of dihedrals: 5845.

backbone.dihedrals
<TopologyGroup containing 5846 dihedrals>

for dih in backbone.dihedrals:
    print dih.atoms

<AtomGroup [<Atom 1: N of type 56 of resname MET, resid 1 and segid 4AKE>, <Atom 5: CA of type 22 of resname MET, resid 1 and segid 4AKE>, <Atom 7: CB of type 23 of resname MET, resid 1 and segid 4AKE>, <Atom 8: HB1 of type 3 of resname MET, resid 1 and segid 4AKE>]>
<AtomGroup [<Atom 1: N of type 56 of resname MET, resid 1 and segid 4AKE>, <Atom 5: CA of type 22 of resname MET, resid 1 and segid 4AKE>, <Atom 7: CB of type 23 of resname MET, resid 1 and segid 4AKE>, <Atom 9: HB2 of type 3 of resname MET, resid 1 and segid 4AKE>]>
<AtomGroup [<Atom 1: N of type 56 of resname MET, resid 1 and segid 4AKE>, <Atom 5: CA of type 22 of resname MET, resid 1 and segid 4AKE>, <Atom 7: CB of type 23 of resname MET, resid 1 and segid 4AKE>, <Atom 10: CG of type 23 of resname MET, resid 1 and segid 4AKE>]>
[... truncated ...]
<AtomGroup [<Atom 3338: HA2 of type 6 of resname GLY, resid 214 and segid 4AKE>, <Atom 3336: CA of type 23 of resname GLY, resid 214 and segid 4AKE>, <Atom 3339: C of type 32 of resname GLY, resid 214 and segid 4AKE>, <Atom 3341: OT2 of type 72 of resname GLY, resid 214 and segid 4AKE>]>

Currently version of MDAnalysis:

Both 0.15 and current 0.16-dev

[richardjgowers]

Hey @mimischi

So ag.dihedrals will return any dihedrals that any atom from ag is present in. So for example, if I had a linear chain system [1 - 2 - 3 - 4 - 5 - 6], and ag = [2 - 3 - 4 - 5], ag.dihedrals would give me [1 - 2 - 3 - 4], [2 - 3 - 4 - 5] and [3 - 4 - 5 - 6]. I'm guessing you'd just expect to get [2 - 3 - 4 - 5] from that?

[kain88-de]

Yes. The expectation is that only complete dihedral angles that are present in an AtomGroup will be returned.

The current behavior is quite a surprise for backbone selections. Because to get the dihredrals of the backbone atoms alone I would use ag.select_atoms('backbone').dihedrals. But the current implementation gives me dihedral angles of atoms that do not belong to the backbone.

[richardjgowers]

So should the shortcuts to topology groups (.bonds, .angles, .torsions and .impropers) only select bonds where all atoms are present in the AtomGroup? I'm starting to think that maybe that way is more intuitive, and can't think of a good case where I'd want the current behaviour..

It's a shame it's not a method, because then we could have had a keyword like ag.bonds(scope='any') to allow the old behaviour.

[kain88-de]

We can could try to copy the behavior of datreant. ag.bonds.any and ag.bonds

[mimischi]

It probably depends on the usage scenario, whether you want to display topology information from atoms not in the AtomGroup:

• When selecting backbone, I do not care about dihedrals formed with side chain atoms. Therefore only topology information for atoms in the selected AtomGroup should appear.
• Assuming I select a specific residue (to extract some information like dihedral angles) with u.select_atoms("resid 5"), both psi and phi contain atoms from the neighbouring residues. Thus I'd want to have topology information for atoms that are not inside the AtomGroup.

[orbeckst]

Reviving an old thread:

So should the shortcuts to topology groups (.bonds, .angles, .torsions and .impropers) only select bonds where all atoms are present in the AtomGroup.

Yes! The safer approach is to only create dihedrals fully inside the AtomGroup.

The user can always create a bigger selection if this is needed. But they have no way to not get the extended behavior.

I'd like to change this; if anyone has a different opinion now is a good time to voice it.

Otherwise the conclusion of this issue is that the current behavior should be changed in the way @richardjgowers said #1264 (comment) .

[richardjgowers]

Ok I think this needs to get added as a deprecation for 2.0 where we can switch the behaviour? I can't really see a clean way to do this except announce a break in behaviour over the 1.x / 2.x versions

[orbeckst]

Yes, I agree with you. Not sure how we will add a deprecation warning without adding nasty code that checks the result after we computed the selections. Perhaps in docs and CHANGELOG only?

[orbeckst]

What would be a postprocessing workaround, e.g.,

dih = restrict_bonded(dihedrals, ag)

so that

len(dih) <= len(dihedrals)

and dih only contains dihedrals fully inside the ag?