# How should I generate a radial distribution function (RDF) for an MD trajectory with python?

I am attempting to analyse the radial distribution function (oxygen to oxygen) for a ~400ps simulation on a droplet of water (MD done with NAMD). I am using the python library MDAnalysis to do this. However, the problem is that the default RDF function (MDAnalysis.analysis.rdf.InterRDF) from MDAnalysis uses periodic boundary conditions. The water droplet simulation I am doing is in vacuum i.e. with no boundary conditions.

The calculation of the RDF seems to have a step where the total number is divided by the box volume to get the density. So, when I am using the function on my system, python is giving me a divide by zero error, as the box volume is zero. The RDF I got is zero throughout because of this. How should I fix this?

I have considered writing my python code to get the RDF, because MDAnalysis provides the coordinates of each atoms for each frame. I was sort of able to do this when there was only one reference atom (i.e. RDF from one atom) (in the jupyter notebook), but I don't know how to get the RDF when there are multiple reference atoms. There are a lot of formulas for RDF that I could find, but I am not sure which one to use, or how to implement them in python. If anyone could help me with this, that would also be greatly appreciated.

I cannot share my raw data as it is an ongoing research project, but I can share the Jupyter Notebook I used for the analysis: https://github.com/ShoubhikRaj/molecular-modelling/blob/main/matter-modelling/rdf-mdanalysis/RDF_analysis-Copy1.ipynb

• P.S. I cannot seem to be able to install MDTraj so I cannot use that for the analysis. Jul 1 at 8:50
• I have no expertise in MDAnalysis, but have a good exposure in tcl scripting and VMD for analysis. Did you try using the VMD plugin for g(r). If not here is the link: ks.uiuc.edu/Research/vmd/plugins/gofrgui Jul 1 at 10:05
• I think there is an option for use PBC. If you uncheck it, maybe it will not take PBC I am not sure whether it works or not, but I would give it a try. Jul 1 at 10:14
• @AdupaVasista I know that the VMD plugin can do this, however, the problem is that I also have to (additionally) calculate the RDF from the centre of mass of the system as well. The g(r) plugin seems to be unable to do anything with the centre of mass, so I had to use MDAnalysis. For consistency, I am now calculating every RDF with MDAnalysis. Jul 1 at 11:16

You could fool the MDAnalysis RDF function InterRDF by setting an artificial cell of a size large enough so that the periodic boundary conditions don't cause any atom to be counted twice.

This can be done by adding the following lines before computing the rdf. It basically creates a cell twice as large as the largest distance in the selected frame of the universe.

d = np.max(u.atoms.positions) - np.min(u.atoms.positions)
u.dimensions = [2*d, 2*d, 2*d, 90, 90, 90]


If your system size fluctuates a lot you can be on the safe size by taking the largest possible distance accross the entire trajectory.

d = np.max([np.max(u.atoms.positions) - np.min(u.atoms.positions) for ts in u.trajectory])


You should only pay attention to the fact that this rdf isn't normalised as you take an arbitrary density when defining the cell. However it shouldn't be too hard to use a prefactor to normalise it if you have a particular definition of the density for your system. The density of the MDAnalysis arbitrary cell can be obtained with

    density = u.atoms.n_atoms / MDAnalysis.lib.mdamath.box_volume(u.dimensions)


The quantity rdf * density is invariant from the size of the arbitrary cell (if it's large enough).

A minimal working example inspired by your notebook could be:

import MDAnalysis as mda
from MDAnalysis.analysis import rdf
import matplotlib.pyplot as plt
import numpy as np

# New two lines to create a large cell
d = np.max(u.atoms.positions) - np.min(u.atoms.positions)
u.dimensions = [2*d, 2*d, 2*d, 90, 90, 90]

oxy = u.select_atoms('name OH2') # get oxygens
hyd = u.select_atoms('type HT') # get hydrogens
O_O = rdf.InterRDF(oxy,oxy,verbose=True)
O_O.run()
plt.plot(O_O.bins, O_O.rdf)
# or if you're interested in the density invariant quantity
density = u.atoms.n_atoms / MDAnalysis.lib.mdamath.box_volume(u.dimensions)
plt.plot(O_O.bins, density * O_O.rdf)

• +1 Thanks! Is there any way to remove the normalisation by density that is done? Because if I put a dummy cell in, then the density is arbitray as you say. I don't really have a predefined density for the system Jul 2 at 3:39