# How to calculate radial distribution functions for MD trajectories without PBC, directly from coordinates with python?

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

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()


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:

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.

• I guess the python part is not particularly important, as I am trying to understand if I got the formula for RDF correct. Jul 30 at 14:31
• 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. Aug 2 at 5:37
• @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. Aug 2 at 13:47
• @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.
– mykd
Aug 25 at 10:41