analysis.distances.between does not return AtomGroup

[rafaelpapp]

Following the docs analysis.distances.between should return an AtomGroup.

(https://docs.mdanalysis.org/2.2.0/documentation_pages/analysis/distances.html)

MDAnalysis.analysis.distances.between(group, A, B, distance)

Return sub group of group that is within distance of both A and B

Maybe I am not using it correctly but if there are no atoms selected the type returned is int.

Is it the intedend behaviour? I was expecting a void AtomGroup.

Best regards,

Rafael

[rafaelpapp]

Here is a sample script:

import MDAnalysis as mda
from MDAnalysis.analysis import distances

u = mda.Universe("bridge.xyz")

frames = slice(0,1,1)
for ts in u.trajectory[frames]:
   oxygens = u.select_atoms('name O', updating=True)  
   reference1 = u.select_atoms('bynum 1')
   reference2 = u.select_atoms('bynum 4')
   bridge = distances.between(oxygens, reference1, reference2, 5.0)
   print(type(bridge))
   bridge = distances.between(oxygens, reference1, reference2, 4.5)
   print(type(bridge))

bridge.xyz is:

11

O      2.0230227722      1.2483216760     -1.6052231937
H      1.3837993207      0.9947679346     -0.8806594849
H      1.6115789669      1.8610493897     -2.2095671105
O      1.0869483562     -0.1508753062      4.9258681035
H      1.3197221480      0.1273109530      5.8521848683
O      0.6661784881      0.8126618640      0.5431825145
H      1.0886486325      1.1448710833      1.3557072572
H     -0.2945576571      0.7452951040      0.6992765482
O      1.8029561511      1.4030538573      2.9294587659
H      1.6929471493      0.9019294377      3.7613001396
H      2.7377147100      1.5147569534      2.8010577514

The prints show:

<class 'MDAnalysis.core.groups.AtomGroup'>
<class 'int'>

[Oliver Beckstein]

That’s a bug. Probably because we used sum() to “add” the groups and I think Python interprets sum(empty, empty) as 0 or 2… can’t test. 

Can you please raise an issue and report what you observed?

If you have a simple example with the test files then that would be even better. 

Thanks,

Oliver 

Am 8/25/22 um 12:37 schrieb rafaelpapp <rafap...@gmail.com>:



--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/bf2720a1-5547-40b4-828a-01c3f9c8871en%40googlegroups.com.

[rafaelpapp]

Thanks!

I am not familiar with the test files but I'll try.

[Oliver Beckstein]

You access test files with code like

import MDAnalysis as mda 

from MDAnalysis.tests import datafiles as data

u = mda.Universe(data.TPR, data.XTC)

and then do something with the universe. There are MANY other files in datafiles. The TPR/XTC combination is a Gromacs simulation of a protein in water with PBC splits. Another often used one is PSF/DCD, a protein without solvent and no PBC.  

Come to think, for this test case, where we need empty AtomGroups we might not really need test files (which is good between tests run faster). Just create empty AtomGroups and see if 

empty = mda.core.groups.AtomGroup()

between(empty, empty, empty, 5)

triggers the error. 

Finally, there’s also code to create an artificial test universe from scratch, which is very fast. This is somewhere in the tests but I forgot where. But might also be useful.

Best,

Oliver 

Am 8/25/22 um 23:23 schrieb rafaelpapp <rafap...@gmail.com>:



To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/0faa92cc-c10f-4bf1-ac0e-d3977805ed2fn%40googlegroups.com.

[rafaelpapp]

Thanks, I opened a new case at the github site.

[rafaelpapp]

Looking at the code in

package/MDAnalysis/analysis/distances.py

the between function reads

   ns_group = AtomNeighborSearch(group)
   resA = set(ns_group.search(A, distance))
   resB = set(ns_group.search(B, distance))
   return sum(sorted(resB.intersection(resA)))

I believe that the use of sets, sorted and sum is unnecessary, maybe something like

   ns_group = AtomNeighborSearch(group)

   resA = ns_group.search(A, distance)
   resB = ns_group.search(B, distance)
   return resB.intersection(resA)

is all we need.

Do you agree?

[Oliver Beckstein]

Thanks for raising issue https://github.com/MDAnalysis/mdanalysis/issues/3794.

Comments below:

On Aug 29, 2022, at 12:28 AM, rafaelpapp <rafap...@gmail.com> wrote:

Looking at the code in 

package/MDAnalysis/analysis/distances.py

the between function reads

   ns_group = AtomNeighborSearch(group) 
   resA = set(ns_group.search(A, distance)) 
   resB = set(ns_group.search(B, distance)) 
   return sum(sorted(resB.intersection(resA)))

I believe that the use of sets, sorted and sum is unnecessary, maybe something like

MDAnalysis selections do two things: they return unique atoms and they return them in order, sorted by index. To be consistent with how other selections work, we want both.

I think we only need to special-case the empty selection.

g = resB.intersection(resA)

return sum(sorted(g)) if g else AtomGroup()

or something similar. (I can’t remember why we need the sum() there… probably a trick to force the output of sorted back into an AtomGroup via and overloaded __add__.)

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/c9b6961a-3c9a-49b6-8004-b3879830cbcan%40googlegroups.com.

--

Oliver Beckstein (he/his/him)

orbe...@gmail.com

[rafaelpapp]

But intersection returns a sorted AtomGroup, isn't it?

At https://docs.mdanalysis.org/dev/documentation_pages/core/groups.html

"All set methods are listed in the table below. These methods treat the groups as sorted and deduplicated sets of Atoms."

[Oliver Beckstein]

If intersection already does all of this then your solution seems complete.

I’d encourage you contribute the fix yourself as a pull request (PR). We’ll review your code, give guidance on any necessary changes (primarily, adding a test or two), and you’ll be listed as one of the contributors to MDAnalysis in the AUTHORS file.

Have you got experience creating a PR for MDAnalysis (or any project)? The User Guide has a section on Contributing to the Main Code Base https://userguide.mdanalysis.org/stable/contributing_code.html .

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/c145ca74-1e6d-416d-a421-fe97d6e883f2n%40googlegroups.com.

--

Oliver Beckstein (he/his/him)

email: orbe...@mdanalysis.org
twitter: @orbeckst

GitHub: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/