MDAnalysis.analysis.align.alignto won't reindex the atomgroup as is documented

[xiki-tempula]

Expected behavior

In the doc (https://docs.mdanalysis.org/1.0.0/documentation_pages/analysis/align.html#MDAnalysis.analysis.align.alignto), it is said:

Warning

The atom order for mobile and reference is only preserved when select is either “all” or None. In any other case, a new selection will be made that will sort the resulting AtomGroup by index and therefore destroy the correspondence between the two groups. It is safest not to mix ordered AtomGroups with selection strings.

So I would imagine that is select keyword is used, the atomgroup will be sorted by index before passing for alignment, which is not the case.

Actual behavior

The atom group won't be sorted even when the select is used.

Code to reproduce the behavior

>>> from MDAnalysisTests.datafiles import PDB_helix
>>> from MDAnalysis.analysis.align import alignto
>>> u = mda.Universe(PDB_helix)
>>> mobile = u.atoms[:4]
>>> mobile.ids
array([1, 2, 3, 4])
>>> mobile.ix
array([0, 1, 2, 3])
>>> ref = sum([u.atoms[3], u.atoms[2], u.atoms[1], u.atoms[0]])
>>> ref.ids
array([4, 3, 2, 1])
>>> ref.ix
array([3, 2, 1, 0])
>>> alignto(mobile, ref, select='bynum 1-4')
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "~/GitHub/mdanalysis/package/MDAnalysis/analysis/align.py", line 492, in alignto
    match_atoms=match_atoms)
  File "~/GitHub/mdanalysis/package/MDAnalysis/analysis/align.py", line 1367, in get_matching_atoms
    raise SelectionError(errmsg)
MDAnalysis.exceptions.SelectionError: Inconsistent selections, masses differ by more than 0.1; mis-matching atoms are shown above.
>>> alignto(mobile, ref, select='bynum 1-3')
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "~/GitHub/mdanalysis/package/MDAnalysis/analysis/align.py", line 492, in alignto
    match_atoms=match_atoms)
  File "~/GitHub/mdanalysis/package/MDAnalysis/analysis/align.py", line 1367, in get_matching_atoms
    raise SelectionError(errmsg)
MDAnalysis.exceptions.SelectionError: Inconsistent selections, masses differ by more than 0.1; mis-matching atoms are shown above.

So mobile and ref are essentially the same things with a different order. Using the select should sort them and thus, give them the same order, which is not the case. Related to #2976

Current version of MDAnalysis

• Which version are you using? 2.0.0-dev0
• Which version of Python (python -V)? 3.6
• Which operating system? osx

[orbeckst]

It is possible that the "will be resorted" comes from a time when select automatically sorted. I am not 100% sure, though.

Is it desirable that the AtomGroups get sorted automatically by index? Or is it actually better to allow users to provide custom sorted AtomGroups together with selections?

[lilyminium]

What happens if you switched the orders, ie ref is sorted? I believe select is only applied to the mobile atom group and that’s why ref is not sorted.

Is it desirable that the AtomGroups get sorted automatically by index? Or is it actually better to allow users to provide custom sorted AtomGroups together with selections?

This is probably a more general question for @coredevs. They could always sort their atomgroups before using it as an input to an analysis though.

…

On Jul 17, 2021, at 5:11 PM, Oliver Beckstein ***@***.***> wrote:  It is possible that the "will be resorted" comes from a time when select automatically sorted. I am not 100% sure, though. Is it desirable that the AtomGroups get sorted automatically by index? Or is it actually better to allow users to provide custom sorted AtomGroups together with selections? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[xiki-tempula]

For me, I think retaining the sorted atomgroup is more desirable.
One could always sort the atomgroup before passing to the alignment.
But if one wants to provide a specific mapping to the alignment, this couldn't be done in the current setting.

[lilyminium]

Could you sort the mobile and reference Atomgroups beforehand and use the default select=None?

…

On Jul 18, 2021, at 3:34 PM, Zhiyi Wu ***@***.***> wrote:  For me, I think retaining the sorted atomgroup is more desirable. One could always sort the atomgroup before passing to the alignment. But if one wants to provide a specific mapping to the alignment, this couldn't be done in the current setting. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[xiki-tempula]

@lilyminium

>>> alignto(mobile, ref, select='bynum 1-4')
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/Users/zhiyiwu/miniconda3/envs/mda_dev/lib/python3.7/site-packages/MDAnalysis/analysis/align.py", line 497, in alignto
    match_atoms=match_atoms)
  File "/Users/zhiyiwu/miniconda3/envs/mda_dev/lib/python3.7/site-packages/MDAnalysis/analysis/align.py", line 1374, in get_matching_atoms
    raise SelectionError(errmsg)
MDAnalysis.exceptions.SelectionError: Inconsistent selections, masses differ by more than 0.1; mis-matching atoms are shown above.
>>> alignto(mobile, ref, select=None)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/Users/zhiyiwu/miniconda3/envs/mda_dev/lib/python3.7/site-packages/MDAnalysis/analysis/align.py", line 497, in alignto
    match_atoms=match_atoms)
  File "/Users/zhiyiwu/miniconda3/envs/mda_dev/lib/python3.7/site-packages/MDAnalysis/analysis/align.py", line 1374, in get_matching_atoms
    raise SelectionError(errmsg)
MDAnalysis.exceptions.SelectionError: Inconsistent selections, masses differ by more than 0.1; mis-matching atoms are shown above.
>>> alignto(ref, mobile,select=None)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/Users/zhiyiwu/miniconda3/envs/mda_dev/lib/python3.7/site-packages/MDAnalysis/analysis/align.py", line 497, in alignto
    match_atoms=match_atoms)
  File "/Users/zhiyiwu/miniconda3/envs/mda_dev/lib/python3.7/site-packages/MDAnalysis/analysis/align.py", line 1374, in get_matching_atoms
    raise SelectionError(errmsg)
MDAnalysis.exceptions.SelectionError: Inconsistent selections, masses differ by more than 0.1; mis-matching atoms are shown above.
>>> alignto(ref, mobile,select='bynum 1-4')
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/Users/zhiyiwu/miniconda3/envs/mda_dev/lib/python3.7/site-packages/MDAnalysis/analysis/align.py", line 497, in alignto
    match_atoms=match_atoms)
  File "/Users/zhiyiwu/miniconda3/envs/mda_dev/lib/python3.7/site-packages/MDAnalysis/analysis/align.py", line 1374, in get_matching_atoms
    raise SelectionError(errmsg)
MDAnalysis.exceptions.SelectionError: Inconsistent selections, masses differ by more than 0.1; mis-matching atoms are shown above.

[lilyminium]

Oh, interesting, select_atoms actually doesn't reorder atoms. I always thought it did and it's documented as doing so. So this is a more general issue. I wondered if I changed this when I made the more general selection class, but it happens with resids too (which I didn't touch):

>>> ref.select_atoms("bynum 1-4").ix
array([3, 2, 1, 0])
>>> ag = u.atoms[[-1, 0, 10, -2, 21, 20, 2]]
>>> ag.resids
array([13,  1,  1, 13,  2,  2,  1])
>>> ag.select_atoms("resid 1 to 2").ix
array([ 0, 10, 21, 20,  2])

I think this comes from a change in 0.19.0, as the test for uniqueness does not check if it's sorted.

	.. versionchanged:: 0.19.0 If the :class:`AtomGroup` is already unique,
            :attr:`AtomGroup.unique` now returns the group itself instead of a
            copy.

This is a pretty fundamental change since I know I've been relying on the sortedness, so it seems important to decide what we want and fix either the code or the docs for 2.0. Personally I'm in favour of select_atoms sorting the atomgroup for me.

Edit: to break it down more, the selection classes use .unique to sort the output -- so it doesn't matter if select_atoms does not return the original atomgroup

[lilyminium]

I believe select is only applied to the mobile atom group and that’s why ref is not sorted.

And I was wrong here, select is applied to both in alignto and AlignTraj. Thanks for correcting me @xiki-tempula

[orbeckst]

Thanks for tracking down the changes in select_atoms — that's a fundamental issue and it does affect 1.1.1. Thanks for raising #3364.