# How to calculate RDFs for MD trajectories without PBC, directly from coordinates? [closed]

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

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): 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):  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.

• 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, 2021 at 5:37
• @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, 2021 at 10:41
• @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)! Dec 31, 2021 at 1:39
• 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? Jan 14, 2022 at 19:03
• 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. Dec 29, 2022 at 16:42

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.