Obtaining Radial Distribution Functions using MDAnalysis

Question

I am running a simple benzene simulation in GROMOS54a7. I want to calculate the RDF of the center of masses of each benzene molecule, using MDAnalysis 1.0.0.

Is this possible? I have create the rdf for the C molecules g_cc(r) using the following code in a Jupyter Notebook:

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

u = MDAnalysis.Universe("739-c6h6-MolDynamics.tpr","739-c6h6-MolDynamics_good-pbc.xtc")
s1 = u.select_atoms("resid 0 and type CAro")
s2 = u.select_atoms("not (resid 0) and type CAro")
rdf = mda.InterRDF(s1, s2)
rdf.run()

I want to take each benzene molecule (each benzene molecule is a residue in my simulation), calculate its COM and run a script like the one above on it. Is it possible to do something like that?

A general question about RDFs: does the method I have used above construct an RDF using every frame of my trajectory? I don't know if this was made clear in the documentation, so I apologize if this is an obvious question.

Thank you for any advice you have!

mateuszb · Accepted Answer · 2021-02-24T10:47:05.247

It would be useful to make it possible to use CG groups as native atoms in order to reuse the analysis tools in MDAnalysis.

Here is a quick fix that mimics the MDAnalysis group and presents a new positions property. The new positions provides the centre of geometry instead of the actual positions. I also overwrite the len to convey that only one bead is being used for the CG element.

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

u = mda.Universe('sys_solv.pdb','prod.dcd')

class CG:
    def __init__(self, ag):
        self.ag = ag
        self.universe = self.ag.universe
        self.trajectory = self.ag.universe.trajectory

    @property
    def positions(self):
        return np.array([self.ag.center_of_geometry()])

    def __len__(self):
        return 1

cg_selection = u.select_atoms('resid 1')
cg_atom = CG(cg_selection.atoms)
waters = u.select_atoms('name O')

rdf = MDAnalysis.analysis.rdf.InterRDF(cg_atom, waters)
rdf.run() 
plt.plot(rdf.bins, rdf.rdf)

Verification: I selected a single atom for the CG bead and it reproduces the original RDF.

MDAnalysis uses the entire trajectory. You might in the documentation the start/stop/step parameters for the .run() function which allows you to narrow down which frames to use specifically.

Obtaining Radial Distribution Functions using MDAnalysis

1 Answers1