I have been running MD simulations on water droplets with an ion in them. I am trying to implement a python 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:

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:


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: 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 the python code?

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.

  • $\begingroup$ I guess the python part is not particularly important, as I am trying to understand if I got the formula for RDF correct. $\endgroup$
    – S R Maiti
    Jul 30 at 14:31
  • 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$ Aug 2 at 5:37
  • $\begingroup$ @NikeDattani Yeah I was surprised that it did not get more than two upvotes after all this time :-| Those pages you linked were very interesting (love the vote compensation idea lol). And I see this type of thing more often in Chem.SE where getting more than 1 upvote is difficult nowadays. $\endgroup$
    – S R Maiti
    Aug 2 at 13:47
  • 1
    $\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 at 10:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.