13
$\begingroup$

I have been running MD simulations on water droplets with an ion in them. I am trying to implement a program that could calculate the radial distribution function from ion to oxygen (or ion to hydrogen)

I use the MDAnalysis Python program to read the trajectory file and provide coordinates of atoms, frame by frame. Following the example given here, I have written a Python code (comments are added to explain each line):

import MDAnalysis as mda
import matplotlib.pyplot as plt
import numpy as np

trj300 = mda.Universe("file.psf","trajectory.dcd") #load files
oxy = trj300.select_atoms('name OH2')
cal = trj300.select_atoms('name CAL') # ion is calcium here

mybinsize = 80 # set up the bins for the histogram
counts300 = np.zeros(mybinsize,dtype=np.int_) #numpy histogram generates same no. of counts as binsize argument
for frame1 in trj300.trajectory[5000:]: # discard initial 100 ns
    dist_vec = oxy.positions - cal.position # from each oxygen coordinate, subtract ion coodinate
    dist_scal = np.linalg.norm(dist_vec,axis=1) # get list of ion-oxygen distances
    tmp_hist, lengths300 = np.histogram(dist_scal,bins=mybinsize,range=(0,18)) # calculate histogram
    counts300 += tmp_hist # collect the histogram counts of each frame into one array

# normalize the count by dividing by number of frames
counts300_mod = counts300/len(trj300.trajectory[5000:])
# calculate the volume of each shell (determined by bin width)
shell_volumes300 = (4/3) * np.pi * (lengths300[1:]**3 - lengths300[:-1]**3)
# normalize the count by dividing each by shell volume (because the shell volume changes with radius)
counts300_mod = counts300_mod/shell_volumes300
# normalize by dividing by the number of selected oxygen atoms
counts300_mod = counts300_mod/len(oxy) # -> should this be done??

Then I plot with pyplot:

plt.plot((lengths300[:-1]+lengths300[1:])/2,counts300_mod,color='red')
plt.show()

and get (click for larger image):

calcium-oxygen radial distribution function

The first peak has an intensity of ~ 0.0004.

However, when I compare it to the one I generated from VMD, the shape matches, but the y-axis does not match (click for larger images):

VMD g(r) plugin

calculated calcium-oxygen rdf from VMD

Here, the first peak appears at the same location on x-axis (~2.5 Angstrom) but it's intensity on y-axis is 0.015, it does not match.

This makes me believe that I have gotten the formulas for normalization of RDF wrong. Clearly the same pattern is visible, which means histogramming and shell-volume normalization was done correctly. Which formulas should I use to calculate the RDF? Have I missed a nomralization procedure in my algorithm?

Note: I have tried removing the normalization for frame number, or the normalization by the number of oxygen atoms, but neither gives the same value as VMD.

If someone knows how to accomplish this with or without Python, I'd be interested in their answers. The programming language or scripting language is not the most important aspect of this question. I want to know the correct algorithm (i.e. formula) for RDF calculation when there is no PBC.

Improve this question
$\endgroup$
17
  • 1
    $\begingroup$ Seems I'm the only one that gave a +1! Welcome to the club of users that have such high rep that no one wants to upvote you anymore 😂😂😂. We've missed out on a lot of HNQ opportunities because of that!. Also related: electronics.meta.stackexchange.com/q/1198/192433 and electronics.meta.stackexchange.com/a/3852/192433. $\endgroup$
    – Nike Dattani
    Aug 2, 2021 at 5:37
  • 2
    $\begingroup$ @SRMaiti Maybe you can try different analysis tools to confirm these results. I would recommend using TRAVIS (travis-analyzer.de). I don't think you should rely too much on VMD for the absolute numbers. $\endgroup$
    – mykd
    Aug 25, 2021 at 10:41
  • 1
    $\begingroup$ @SRMaiti Okay. I've changed the title to say "RDF" instead of "Radial Distribution Function" so that anyone who works with RDFs will see their favorite acronym jump out at them. Based on your first comment, and your most recent comment, it looks like the Python aspect isn't so important, so I've removed that too in order to make the question appeal to a wider audience (non-Python users too). If there's any edit you can make to help the community answer this, it would be appreciated (you seem to have fortified your understanding of things but haven't edited the question since the last comment)! $\endgroup$
    – Nike Dattani
    Dec 31, 2021 at 1:39
  • 2
    $\begingroup$ If the cumulative distributions match, then the issue could be in the bin widths: g(r) as plotted is in units of per-box-width. Debugging-wise, what happens if you set the bin width in the python computation to be the same as in VMD, or alternatively standardize both plots of g(r) to be in units of per-angstrom or something like that by dividing out the dr in the g(r) plots? $\endgroup$
    – Andrey Poletayev
    Jan 14, 2022 at 19:03
  • 2
    $\begingroup$ Since the comment chain is getting long, I'd like to remind everyone not to create a new chat room (i.e. don't click the button when the system tells you about moving to chat). The existing room is perfectly appropriate. $\endgroup$
    – Nike Dattani
    Dec 29, 2022 at 16:42

1 Answer 1

Reset to default
0
$\begingroup$

An RDF is the ratio of local density to bulk density of your system. You can look for the equation in number of good books (e.g. Allen and Tidsley).

For a simulated system with a defined periodic boundary, you need to consider the effect of PBCs otherwise you might end up getting incorrect distances for atoms present at the surface with elongated bonds.

To solve this you can again look for a very basic code available in the Allen and Tidsely book.

The results that you are getting might differ due to the size of bin one chooses for binning the atoms at a particular distance.

Improve this answer
$\endgroup$
3
  • 4
    $\begingroup$ As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center. $\endgroup$
    – Community Bot
    Dec 30, 2022 at 13:37
  • 1
    $\begingroup$ Ideally answers here will provide more info than just pointing to an external reference (e.g. include the relevant equation, provide some commentary clarifying the descriptions in the resources). If it is not possible to reproduce (with citation) or paraphrase a resource, one should try to provide a precise reference (i.e full title of book, chapter and/or page number) $\endgroup$
    – Tyberius ♦
    Dec 30, 2022 at 15:42
  • $\begingroup$ The questions says: "This makes me believe that I have gotten the formulas for normalization of RDF wrong. [...] Which formulas should I use to calculate the RDF?" Can you write the formula for the normalization/calculating-the-RDF in your answer, rather than referring the user to a book which they might not have available? "To solve this you can again look for a very basic code available in the Allen and Tidsely book." Can you provide that "basic" code? Your 2nd paragraph talks about PBC but "I want to know the correct algorithm (i.e. formula) for RDF calculation when there is no PBC" $\endgroup$
    – Nike Dattani
    Dec 30, 2022 at 15:47

Not the answer you're looking for? Browse other questions tagged .