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
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
and the average number of \(b\) particles within radius \(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)\).
We provide options for calculating the density of particle \(b\) in a shell at distance \(r\) around a \(a\) particle, which is
- class MDAnalysis.analysis.rdf.InterRDF(g1, g2, nbins=75, range=(0.0, 15.0), norm='rdf', exclusion_block=None, **kwargs)[source]
Radial distribution function
InterRDF
is a tool to calculate average radial distribution functions between two groups of atoms. Suppose we have two AtomGroupsA
andB
.A
contains atomA1
,A2
, andB
containsB1
,B2
. GivenA
andB
toInterRDF
, the output will be the average of RDFs betweenA1
andB1
,A1
andB2
,A2
andB1
,A2
andB2
. A typical application is to calculate the RDF of solvent with itself or with another solute.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
results.rdf
andresults.count
.- Parameters:
g1 (AtomGroup) – First AtomGroup
g2 (AtomGroup) – Second AtomGroup
nbins (int) – Number of bins in the histogram
norm (str, {'rdf', 'density', 'none'}) –
For ‘rdf’ calculate \(g_{ab}(r)\). For ‘density’ the single particle density \(n_{ab}(r)\) is computed. ‘none’ computes the number of particles occurences in each spherical shell.
New in version 2.3.0.
exclusion_block (tuple) – A tuple representing the tile to exclude from the distance array.
verbose (bool) – Show detailed progress of the calculation if set to True
- results.bins
numpy.ndarray
of the centers of the nbins histogram bins.New in version 2.0.0.
- Type:
- 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.- Type:
- results.edges
numpy.ndarray
of the nbins + 1 edges of the histogram bins.New in version 2.0.0.
- Type:
- 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.- Type:
- results.rdf
numpy.ndarray
of the radial distribution function values for theresults.bins
.New in version 2.0.0.
- Type:
- 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.- Type:
- results.count
numpy.ndarray
representing the radial histogram, i.e., the raw counts, for allresults.bins
.New in version 2.0.0.
- Type:
- 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.- Type:
Example
First create the
InterRDF
object, by supplying two AtomGroups then use therun()
methodrdf = InterRDF(ag1, ag2) rdf.run()
Results are available through the
results.bins
andresults.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
.
- class MDAnalysis.analysis.rdf.InterRDF_s(u, ags, nbins=75, range=(0.0, 15.0), norm='rdf', density=False, **kwargs)[source]
Site-specific radial distribution function
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 sameA
andB
toInterRDF_s
, the output will be a list of individual RDFs betweenA1
andB1
,A1
andB2
,A2
andB1
,A2
andB2
(and similarly forC
andD
). 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.- Parameters:
u (Universe) –
a Universe that contains atoms in ags
Deprecated since version 2.3.0: This parameter is superflous and will be removed in MDAnalysis 3.0.0.
nbins (int) – Number of bins in the histogram
norm (str, {'rdf', 'density', 'none'}) –
For ‘rdf’ calculate \(g_{ab}(r)\). For ‘density’ the single particle density \(n_{ab}(r)\) is computed. ‘none’ computes the number of particles occurences in each spherical shell.
New in version 2.3.0.
density (bool) –
False: calculate \(g_{ab}(r)\); True: calculate the true single particle density \(n_{ab}(r)\). density overwrites the norm parameter.
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.
Deprecated since version 2.3.0: Instead of density=True use norm=’density’
- 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.
- Type:
- 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.- Type:
- results.edges
array of the
nbins + 1
edges of the histogram bins; all individual site-specific RDFs have the same bins.New in version 2.0.0.
- Type:
- 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.- Type:
- results.rdf
list
of the site-specific radial distribution functions if norm=’rdf’ or density functions for thebins
if norm=’density’. The list containslen(ags)
entries. Each entry for thei
-th pair [A, B] = ags[i] in ags is anumpy.ndarray
with shape(len(A), len(B))
, i.e., a stack of RDFs. For example,results.rdf[i][0, 2]
is the RDF between atomsA[0]
andB[2]
.New in version 2.0.0.
- Type:
- 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.- Type:
- results.count
list
of the site-specific radial histograms, i.e., the raw counts, for allresults.bins
. The data have the same structure asresults.rdf
except that the arrays contain the raw counts.New in version 2.0.0.
- Type:
- 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.- Type:
- results.cdf
list
of the site-specific cumulative counts, for allresults.bins
. The data have the same structure asresults.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.
- Type:
- 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.- Type:
Example
First create the
InterRDF_s
object, by supplying one Universe and one list of pairs of AtomGroups, then use therun()
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
andresults.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 ins2
)To generate the cumulative distribution function (cdf) in the sense of “particles within radius \(r\)”, i.e., \(N_{ab}(r)\), use the
get_cdf()
methodcdf = 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 ins2
)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
.Changed in version 2.3.0: Introduce norm and exclusion_blocks attributes.
Deprecated since version 2.3.0: Instead of density=True use norm=’density’
Deprecated since version 2.3.0: The universe parameter is superflous.
- 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: