10
$\begingroup$

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?

RDF of O-O

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

$\endgroup$
4
  • $\begingroup$ P.S. I cannot seem to be able to install MDTraj so I cannot use that for the analysis. $\endgroup$
    – S R Maiti
    Jul 1 at 8:50
  • $\begingroup$ 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 $\endgroup$
    – Vasista
    Jul 1 at 10:05
  • $\begingroup$ 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. $\endgroup$
    – Vasista
    Jul 1 at 10:14
  • $\begingroup$ @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. $\endgroup$
    – S R Maiti
    Jul 1 at 11:16
9
$\begingroup$

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)
$\endgroup$
3
  • $\begingroup$ +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 $\endgroup$
    – S R Maiti
    Jul 2 at 3:39
  • $\begingroup$ @ShoubhikRMaiti I've made an edit, does it answer your question? $\endgroup$
    – Hebo
    Jul 2 at 12:55
  • $\begingroup$ Yes that solves it, thanks! $\endgroup$
    – S R Maiti
    Jul 3 at 4:12

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .