MDAnalysis / pmda

Parallel algorithms for MDAnalysis
https://www.mdanalysis.org/pmda/
Other
31 stars 22 forks source link

Wrong normalization for cdf calculation #117

Closed VOD555 closed 4 years ago

VOD555 commented 4 years ago

Expected behaviour

CDF represents the average number of atoms in atomgroup2 around one atom in atomgroup1.

Actual behaviour

CDF was the total number of atoms in ag2 around all atoms in ag1.

Code to reproduce the behaviour

import MDAnalysis as mda
u = mda.Universe(top, trj)

from pmda.rdf import InterRDF

ag1 = u.select_atoms('name OH2')
ag2 = u.select_atoms('name H1')
RDF = interRDF(ag1, ag2)
RDF.run()
RDF.cdf()

....
orbeckst commented 4 years ago

So cdf() is wrong when len(ag1) > 1 because it may count the same atom in ag2 at different distances from atoms in ag1?

How do you propose to change the behavior?

orbeckst commented 4 years ago

What about the RDF – is that one calculated as the average over all RDFs (ag1[i], ag2), e.g., RDF(all Na+, all OW) gives the average NA+-OW RDF?

VOD555 commented 4 years ago

Yes, RDF is calculated as the average over N1*N2 RDFs, where N1 is the number of atom in ag1 and N2 is that of ag2.

VOD555 commented 4 years ago

I think CDF should be the integrated RDF over the volume, which is the current CDF divided by N1*N2.

And we may add another attribute/function to calculate the average number of ag2 atoms within a certain range from ag1 atoms, which is the current CDF divided by N1.

orbeckst commented 4 years ago

@VOD555 was this the same issue as #120 ??? Was this fixed already or is it still open? Or is this something else?

VOD555 commented 4 years ago

Yes, this is the same problem as #120. This only mentioned CDF, but the problem is from the wrong rdf. I'll close it.

VOD555 commented 4 years ago

Fixed by PR #121