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

Question

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

Answer 1

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
u = mda.Universe("15_ws.psf","100H2O_prelim.dcd") #load files

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