4.7.3. Radial Distribution Functions — MDAnalysis.analysis.rdf

This module contains two classes to calculate radial pair distribution functions (radial distribution functions or “RDF”). The RDF \(g_{ab}(r)\) between types of particles \(a\) and \(b\) is

\[g_{ab}(r) = (N_{a} N_{b})^{-1} \sum_{i=1}^{N_a} \sum_{j=1}^{N_b} \langle \delta(|\mathbf{r}_i - \mathbf{r}_j| - r) \rangle\]

which is normalized so that the RDF becomes 1 for large separations in a homogenous system. The RDF effectively counts the average number of \(b\) neighbours in a shell at distance \(r\) around a \(a\) particle and represents it as a density.

The radial cumulative distribution function is

\[G_{ab}(r) = \int_0^r \!\!dr' 4\pi r'^2 g_{ab}(r')\]

and the average number of \(b\) particles within radius \(r\)

\[N_{ab}(r) = \rho G_{ab}(r)\]

(with the appropriate density \(\rho\)). The latter function can be used to compute, for instance, coordination numbers such as the number of neighbors in the first solvation shell \(N(r_1)\) where \(r_1\) is the position of the first minimum in \(g(r)\).

4.7.3.1. Average radial distribution function

InterRDF is a tool to calculate average radial distribution functions between two groups of atoms. Suppose we have two AtomGroups A and B. A contains atom A1, A2, and B contains B1, B2. Give A and B to class:InterRDF, the output will be the average of RDFs between A1 and B1, A1 and B2, A2 and B1, A2 and B2. A typical application is to calculate the RDF of solvent with itself or with another solute.

class MDAnalysis.analysis.rdf.InterRDF(g1, g2, nbins=75, range=(0.0, 15.0), exclusion_block=None, **kwargs)[source]

Intermolecular pair distribution function

InterRDF(g1, g2, nbins=75, range=(0.0, 15.0))

Parameters

  • g1 (AtomGroup) – First AtomGroup

  • g2 (AtomGroup) – Second AtomGroup

  • nbins (int (optional)) – Number of bins in the histogram [75]

  • range (tuple or list (optional)) – The size of the RDF [0.0, 15.0]

  • exclusion_block (tuple (optional)) – A tuple representing the tile to exclude from the distance array. [None]

  • start (int (optional)) – The frame to start at (default is first)

  • stop (int (optional)) – The frame to end at (default is last)

  • step (int (optional)) – The step size through the trajectory in frames (default is every frame)

  • verbose (bool (optional)) – Show detailed progress of the calculation if set to True; the default is False.

Example

First create the InterRDF object, by supplying two AtomGroups then use the run() method

rdf = InterRDF(ag1, ag2)
rdf.run()

Results are available through the bins and rdf attributes:

plt.plot(rdf.bins, rdf.rdf)

The exclusion_block keyword allows the masking of pairs from within the same molecule. For example, if there are 7 of each atom in each molecule, the exclusion mask (7, 7) can be used.

New in version 0.13.0.

Parameters

  • trajectory (mda.Reader) – A trajectory Reader

  • verbose (bool, optional) – Turn on more logging and debugging, default False

run(start=None, stop=None, step=None, verbose=None)

Perform the calculation

Parameters

  • start (int, optional) – start frame of analysis

  • stop (int, optional) – stop frame of analysis

  • step (int, optional) – number of frames to skip between each analysed frame

  • verbose (bool, optional) – Turn on verbosity

4.7.3.2. Site-specific radial distribution function

InterRDF_s calculates site-specific radial distribution functions. Instead of two groups of atoms it takes as input a list of pairs of AtomGroup, [[A, B], [C, D], ...]. Give the same A and B to InterRDF_s, the output will be a list of RDFs between A1 and B1, A1 and B2, A2 and B1, A2 and B2 (and similarly for C and D). These site-specific radial distribution functions are typically calculated if one is interested in the solvation shells of a ligand in a binding site or the solvation of specific residues in a protein. A common use case is to choose A and C to be AtomGroups that only contain a single atom and W all solvent molecules: InterRDF_s(u, [[A, W], [B, W]]) will then produce the RDF of solvent around atom A[0] and around atom B[0].

class MDAnalysis.analysis.rdf.InterRDF_s(u, ags, nbins=75, range=(0.0, 15.0), density=True, **kwargs)[source]

Site-specific intermolecular pair distribution function

Parameters

  • u (Universe) – a Universe that contains atoms in ags

  • ags (list) – a list of pairs of AtomGroup instances

  • nbins (int (optional)) – Number of bins in the histogram [75]

  • range (tuple or list (optional)) – The size of the RDF [0.0, 15.0]

  • start (int (optional)) – The frame to start at (default is first)

  • stop (int (optional)) – The frame to end at (default is last)

  • step (int (optional)) – The step size through the trajectory in frames (default is every frame)

Example

First create the InterRDF_s object, by supplying one Universe and one list of pairs of AtomGroups, then use the run() method:

from MDAnalysisTests.datafiles import GRO_MEMPROT, XTC_MEMPROT
u = mda.Universe(GRO_MEMPROT, XTC_MEMPROT)

s1 = u.select_atoms('name ZND and resid 289')
s2 = u.select_atoms('(name OD1 or name OD2) and resid 51 and sphzone 5.0 (resid 289)')
s3 = u.select_atoms('name ZND and (resid 291 or resid 292)')
s4 = u.select_atoms('(name OD1 or name OD2) and sphzone 5.0 (resid 291)')
ags = [[s1, s2], [s3, s4]]

rdf = InterRDF_s(u, ags)
rdf.run()

Results are available through the bins and rdf attributes:

plt.plot(rdf.bins, rdf.rdf[0][0][0])

(Which plots the rdf between the first atom in s1 and the first atom in s2)

To generate the cumulative distribution function (cdf), use the get_cdf() method

cdf = rdf.get_cdf()

Results are available through the :attr:’cdf’ attribute:

plt.plot(rdf.bins, rdf.cdf[0][0][0])

(Which plots the cdf between the first atom in s1 and the first atom in s2)

New in version 0.19.0.

Parameters

  • trajectory (mda.Reader) – A trajectory Reader

  • verbose (bool, optional) – Turn on more logging and debugging, default False

get_cdf()[source]

Calculate the cumulative distribution functions (CDF) for all sites. Note that this is the actual count within a given radius, i.e., \(N(r)\). :returns: cdf – list of arrays with the same structure as rdf :rtype: list

run(start=None, stop=None, step=None, verbose=None)

Perform the calculation

Parameters

  • start (int, optional) – start frame of analysis

  • stop (int, optional) – stop frame of analysis

  • step (int, optional) – number of frames to skip between each analysed frame

  • verbose (bool, optional) – Turn on verbosity