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:
The first peak has an intensity of ~ 0.0004.
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.