4.7.2.1. 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)\).

In InterRDF_s, we provide an option density. When density is False, it will return the RDF \(g_{ab}(r)\); when density is True, it will return the density of particle \(b\) in a shell at distance \(r\) around a \(a\) particle, which is

\[n_{ab}(r) = \rho g_{ab}(r)\]

4.7.2.1.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. Given A and B to 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

The radial distribution function is calculated by histogramming distances between all particles in g1 and g2 while taking periodic boundary conditions into account via the minimum image convention.

The exclusion_block keyword may be used to exclude a set of distances from the calculations.

Results are available in the attributes rdf and count.

Parameters

  • g1 (AtomGroup) – First AtomGroup

  • g2 (AtomGroup) – Second AtomGroup

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

  • range (tuple or list (optional)) – The size of the RDF

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

  • verbose (bool (optional)) – Show detailed progress of the calculation if set to True

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 results.bins and results.rdf attributes:

plt.plot(rdf.results.bins, rdf.results.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.

Changed in version 1.0.0: Support for the start, stop, and step keywords has been removed. These should instead be passed to InterRDF.run().

Changed in version 2.0.0: Store results as attributes bins, edges, rdf and count of the results attribute of AnalysisBase.

results.bins

numpy.ndarray of the centers of the nbins histogram bins.

New in version 2.0.0.

bins

Alias to the results.bins attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.bins instead.

results.edges

numpy.ndarray of the nbins + 1 edges of the histogram bins.

New in version 2.0.0.

edges

Alias to the results.edges attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.edges instead.

results.rdf

numpy.ndarray of the radial distribution function values for the results.bins.

New in version 2.0.0.

rdf

Alias to the results.rdf attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.rdf instead.

results.count

numpy.ndarray representing the radial histogram, i.e., the raw counts, for all results.bins.

New in version 2.0.0.

count

Alias to the results.count attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.count instead.

run(start=None, stop=None, step=None, frames=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

  • frames (array_like, optional) – array of integers or booleans to slice trajectory

  • verbose (bool, optional) – Turn on verbosity

Changed in version 2.1.0: Added ability to iterate through trajectory by passing a list of frame indices

4.7.2.1.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], ...]. Given the same A and B to InterRDF_s, the output will be a list of individual 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.

class MDAnalysis.analysis.rdf.InterRDF_s(u, ags, nbins=75, range=(0.0, 15.0), density=False, **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

  • range (tuple or list (optional)) – The size of the RDF

  • density (bool (optional)) –

    False: calculate \(g_{ab}(r)\); True: calculate the true single particle density \(n_{ab}(r)\).

    New in version 1.0.1: This keyword was available since 0.19.0 but was not documented. Furthermore, it had the opposite meaning. Since 1.0.1 it is officially supported as documented.

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 results.bins and results.rdf attributes:

plt.plot(rdf.results.bins, rdf.results.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) in the sense of “particles within radius \(r\)”, i.e., \(N_{ab}(r)\), use the get_cdf() method

cdf = rdf.get_cdf()

Results are available through the results.cdf attribute:

plt.plot(rdf.results.bins, rdf.results.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.

Changed in version 1.0.0: Support for the start, stop, and step keywords has been removed. These should instead be passed to InterRDF_s.run().

Changed in version 2.0.0: Store results as attributes bins, edges, rdf, count and cdf of the results attribute of AnalysisBase.

results.bins

numpy.ndarray of the centers of the nbins histogram bins; all individual site-specific RDFs have the same bins.

New in version 2.0.0.

bins

Alias to the results.bins attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.bins instead.

results.edges

numpy.ndarray of the nbins + 1 edges of the histogram bins; all individual site-specific RDFs have the same bins.

New in version 2.0.0.

edges

Alias to the results.edges attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.edges instead.

results.rdf

list of the site-specific radial distribution functions or density functions for the bins. The list contains len(ags) entries. Each entry for the i-th pair [A, B] = ags[i] in ags is a numpy.ndarray with shape (len(A), len(B)), i.e., a stack of RDFs. For example, results.rdf[i][0, 2] is the RDF between atoms A[0] and B[2].

New in version 2.0.0.

rdf

Alias to the results.rdf attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.rdf instead.

results.count

list of the site-specific radial histograms, i.e., the raw counts, for all results.bins. The data have the same structure as results.rdf except that the arrays contain the raw counts.

New in version 2.0.0.

count

Alias to the results.count attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.count instead.

results.cdf

list of the site-specific cumulative counts, for all results.bins. The data have the same structure as results.rdf except that the arrays contain the cumulative counts.

This attribute only exists after get_cdf() has been run.

New in version 2.0.0.

cdf

Alias to the results.cdf attribute.

Deprecated since version 2.0.0: This attribute will be removed in 3.0.0. Use results.cdf instead.

get_cdf()[source]

Calculate the cumulative counts for all sites.

This is the cumulative count within a given radius, i.e., \(N_{ab}(r)\).

The result is returned and also stored in the attribute results.cdf.

Returns

cdf – list of arrays with the same structure as results.rdf

Return type

list

run(start=None, stop=None, step=None, frames=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

  • frames (array_like, optional) – array of integers or booleans to slice trajectory

  • verbose (bool, optional) – Turn on verbosity

Changed in version 2.1.0: Added ability to iterate through trajectory by passing a list of frame indices