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