Re: [mdnalysis-discussion] calculate rdf when only positions of residues available

[Oliver Beckstein]

Hello Poojitha Ramachandra,

Welcome to the MDAnalysis maling list!

On Apr 8, 2021, at 4:32 AM, Poojitha Ramachandra <poojitha.r...@gmail.com> wrote:

Hi,

I have an array of positions of molecules of certain residue group

I want to calculate the rdf among these molecules(I dont care what kind of molecules these are)  

I want to use the api this way:

  from MDAnalysis.analysis.rdf import InterRDF

    rdf = InterRDF(pos,pos,range=(0,10.05),nbins=500)

InterRDF() takes AtomGroups as input, not distances. It will iterate over the underlying trajectory. If you want the RDF between molecules then you need to select one specific atom in each molecule. More details in the user guide:

- RDF https://userguide.mdanalysis.org/stable/examples/analysis/structure/average_rdf.html

- atom-to-atom https://userguide.mdanalysis.org/stable/examples/analysis/structure/site_specific_rdf.html

earlier I was able to calculate the distance between the molecules using:

  l=20.1

  return MDAnalysis.lib.distances.distance_array(positions,positions,box=(l,l,l,90,90,90))

If you only want to do this for a single frame with your positions then you can use distance_array (or self_distance_array()) and histogram. See the actual code for InterRDF https://docs.mdanalysis.org/stable/_modules/MDAnalysis/analysis/rdf.html#InterRDF, specifically the _single_frame() method.

important criteria is that I have only positions of molecules in a  system

You mean you don’t have an AtomGroup?

If you don’t have a “topology” of your system, there’s not much that standard MDAnalysis will do for you. I would probably generate a minimal Universe (using Universe.empty()  and associate your coordinates with your new universe. Then you can use MDAnalysis as usual. If you need help with this approach, ask, and we’ll point you in the right direction.

Hope that helps,

Oliver

--

Oliver Beckstein (he/his/him)

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

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/

[Poojitha Ramachandra]

Hi Oliver,

Thank you for the quick response.

Could you let me know how I can create a dummy topology using positions just to calculate rdf easily.[I have positions of 500 molecules belonging only 1 kind of residue ]

Thanks.

[Oliver Beckstein]

A pretty extensive explanation is https://userguide.mdanalysis.org/1.0.0/examples/constructing_universe.html ; a shorter version is https://userguide.mdanalysis.org/1.0.0/universe.html#constructing-from-scratch

Something like

N = len(pos)

u = mda.Universe.empty(N, n_residues=N, trajectory=True)  # treat each coord as own atom and residue

u.load_new([pos])

might work. Note that if pos has shape (N, 3) then it represents a single frame and when loading with load_new() you probably need to make it a 1-frame trajectory so that’s why I put it into a list. If the MemoryReader complains, cast it to a np.array([pos]) or add the frame axis with pos[np.newaxis, :, 3].

Oliver

--
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/a3feef81-2139-423e-9e13-894fa4676d11n%40googlegroups.com.

--

Oliver Beckstein (he/his/him)

orbe...@gmail.com

[Poojitha Ramachandra]

HI,

I am using the api this way:

    from MDAnalysis.analysis.rdf import InterRDF

    N = pos.shape[0]

    print("size of array : ",N)

    u = MDAnalysis.Universe.empty(N, n_residues=N, trajectory=True)  # treat each coord as own atom and residue

    u.load_new(np.array([pos]))

    res=u.atoms

    rdf = InterRDF(res,res,range=(0,10.05),nbins=500,exclusion_block=(1,1))#

    rdf.run()

    plt.plot(rdf.bins, rdf.rdf)

But I get an odd graph. what do you think is the reason here?

You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/Dj7nzvDwfTc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/241B0848-A487-41EA-BB37-9589B4F27C88%40gmail.com.

--

Regards,

Poojitha Ramachandra

[Oliver Beckstein]

This looks reasonable at first glance. 

Compare pos and res.positions.  Are they the same?

Try 

D = self_distance_array(res.pos)

H, edges = np.histogram(D)

plt.plot(edges[:-1], H)

To see what should come out of rdf (ignoring normalization etc)

Oliver

Am 4/9/21 um 10:53 schrieb Poojitha Ramachandra <poojitha.r...@gmail.com>:



HI,

I am using the api this way:

    from MDAnalysis.analysis.rdf import InterRDF

    N = pos.shape[0]

    print("size of array : ",N)

    u = MDAnalysis.Universe.empty(N, n_residues=N, trajectory=True)  # treat each coord as own atom and residue

    u.load_new(np.array([pos]))

    res=u.atoms

    rdf = InterRDF(res,res,range=(0,10.05),nbins=500,exclusion_block=(1,1))#

    rdf.run()

    plt.plot(rdf.bins, rdf.rdf)

But I get an odd graph. what do you think is the reason here?

<image.png>

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAKER%2BOG5ygQg4k3Gh3gX_ciYYFDkyZ0fSx%2BD2eyqCoXa2vy-7Q%40mail.gmail.com.

[Poojitha Ramachandra]

u.dimensions is equal to [0.0,0.0,0.0,0.0,0.0,0.0] is this causing a problem? 

pos and res.positions are qual

    D = MDAnalysis.lib.distances.self_distance_array(res.positions)

    H, edges = np.histogram(D)

    plt.plot(edges[:-1], H)

gives this:

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/C717D309-CE11-4F2F-9C21-61C2F79F824E%40gmail.com.

--

Regards,

Poojitha Ramachandra

[Oliver Beckstein]

On Apr 9, 2021, at 1:31 PM, Poojitha Ramachandra <poojitha.r...@gmail.com> wrote:

u.dimensions is equal to [0.0,0.0,0.0,0.0,0.0,0.0] is this causing a problem?

That could be a problem; I think the RDF code always does PBC minimum image convention for distances. 

Try setting your dimensions in your custom universe and see if it helps. If so, you can raise an issue for adding a check in the RDF code if the box is sensible. 

pos and res.positions are qual

Good.

    D = MDAnalysis.lib.distances.self_distance_array(res.positions)

    H, edges = np.histogram(D)

    plt.plot(edges[:-1], H)

gives this:

That looks more reasonable. I’d check the box dimensions.

You can also try giving self_distance_array() the argument box=np.zeros(6) and see if that leads to the same effect. (RDF uses capped_distances() instead of self_distance_array() IIRC so you can try that, too)

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAKER%2BOG9xzZgZzTLP5zG%3DGL5or%2B_QbC25H3jOZZ3dsJ5fCSbhw%40mail.gmail.com.

[Poojitha Ramachandra]

hi,

the box dimensions was the issue.

I explicitly set the value u.dimensions=[20.1,20.1,20.1,90,90,90] and now I get a graph.

but how do I compare this, with the rdf generated from gromacs for the same data?

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/2CF48FA4-16A4-459C-AE31-A9C5D578E467%40mdanalysis.org.

--

Regards,

Poojitha Ramachandra

[Poojitha Ramachandra]

gentle reminder

[Oliver Beckstein]

Hi Poojitha Ramachandra,

The core developers are extremely busy at the moment with preparations for a new release, a workshop preparation, grant writing — and doing the work they actually get paid for in addition to all of that. Please be patient. You’re not being ignored, it’s just that at least core devs have absolutely no time to spare at the moment.

Maybe some of our users can help here.

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/b607d215-60f6-4400-9102-9b9b791f046fn%40googlegroups.com.