analysis.distances.between does not return AtomGroup

30 views
Skip to first unread message

rafaelpapp

unread,
Aug 25, 2022, 1:37:23 PMAug 25
to MDnalysis discussion

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

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

unread,
Aug 25, 2022, 2:40:11 PMAug 25
to MDnalysis discussion
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))

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

unread,
Aug 25, 2022, 6:09:20 PMAug 25
to mdnalysis-...@googlegroups.com
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

unread,
Aug 26, 2022, 2:23:12 AMAug 26
to MDnalysis discussion
Thanks!
I am not familiar with the test files but I'll try.

Oliver Beckstein

unread,
Aug 26, 2022, 3:47:37 AMAug 26
to mdnalysis-...@googlegroups.com
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>:



rafaelpapp

unread,
Aug 28, 2022, 3:50:25 PMAug 28
to MDnalysis discussion
Thanks, I opened a new case at the github site.

rafaelpapp

unread,
Aug 29, 2022, 3:28:35 AMAug 29
to MDnalysis discussion
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

unread,
Aug 29, 2022, 1:16:46 PMAug 29
to mdnalysis-discussion
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

rafaelpapp

unread,
Aug 29, 2022, 1:59:33 PMAug 29
to MDnalysis discussion
But intersection returns a sorted AtomGroup, isn't it?

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

Oliver Beckstein

unread,
Aug 29, 2022, 4:22:26 PMAug 29
to mdnalysis-discussion
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


--
Oliver Beckstein (he/his/him)

GitHub: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project





Reply all
Reply to author
Forward
0 new messages