RadialDistributionFunction¶

class mdhelper.analysis.structure.RadialDistributionFunction(ag1: AtomGroup, ag2: AtomGroup = None, n_bins: int = 201, range: tuple[float] = (0.0, 15.0), *, drop_axis: int | str = None, norm: str = 'rdf', exclusion: tuple[int] = None, groupings: str | tuple[str] = 'atoms', reduced: bool = False, n_batches: int = None, parallel: bool = False, verbose: bool = True, **kwargs)[source]¶

Bases: DynamicAnalysisBase

Serial and parallel implementations to calculate the radial distribution function (RDF) \(g_{ij}(r)\) between types \(i\) and \(j\) and its related properties for two- and three-dimensional systems.

The RDF is given by

\[\begin{split}g_{ij}^\mathrm{3D}(r)=\frac{V}{4\pi r^2N_iN_j} \sum_{\alpha=1}^{N_i}\sum_{\beta=1}^{N_j}\left\langle \delta\left(|\mathbf{r}_\alpha-\mathbf{r}_\beta|-r\right) \right\rangle\\ g_{ij}^\mathrm{2D}(r)=\frac{A}{2\pi rN_iN_j} \sum_{\alpha=1}^{N_i}\sum_{\beta=1}^{N_j}\left\langle \delta\left(|\mathbf{r}_\alpha-\mathbf{r}_\beta|-r\right) \right\rangle\end{split}\]

where \(V\) and \(A\) are the system volume and area, \(N_i\) and \(N_j\) are the number of particles, and \(\mathbf{r}_\alpha\) and \(\mathbf{r}_\beta\) are the positions of particles \(\alpha\) and \(\beta\) belonging to species \(i\) and \(j\), respectively. The RDF is normalized such that \(\lim_{r\rightarrow\infty}g_{ij}(r)=1\) in a homogeneous system.

(A closely related quantity is the single particle density \(n_{ij}(r)=\rho_jg_{ij}(r)\), where \(\rho_j\) is the number density of species \(j\).)

The cumulative RDF is

\[\begin{split}G_{ij}^\mathrm{3D}(r)=4\pi\int_0^rR^2g_{ij}(R)\,dR\\ G_{ij}^\mathrm{2D}(r)=2\pi\int_0^rRg_{ij}(R)\,dR\end{split}\]

and the average number of \(j\) particles found within radius \(r\) is

\[N_{ij}(r)=\rho_jG_{ij}(r)\]

The expression above can be used to compute the coordination numbers (number of neighbors in each solvation shell) by setting \(r\) to the \(r\)-values where \(g_{ij}(r)\) is locally minimized, which signify the solvation shell boundaries.

The RDF can also be used to obtain other relevant structural properties, such as

the potential of mean force

\[w_{ij}(r)=-k_\mathrm{B}T\ln{g_{ij}(r)}\]

where \(k_\mathrm{B}\) is the Boltzmann constant and \(T\) is the system temperature, and
the (partial) static structure factor (see calculate_structure_factor() for the possible definitions).

Parameters:

ag1MDAnalysis.AtomGroup

First atom group \(i\).

ag2MDAnalysis.AtomGroup

Second atom group \(j\).

n_binsint, default: 201

Number of histogram bins \(N_\mathrm{bins}\).

rangearray-like, default: (0.0, 15.0)

Range of radii values. The upper bound should be less than half the largest system dimension.

Shape: \((2,)\).

Reference unit: \(\mathrm{Å}\).

drop_axisint or str, keyword-only, default: 2

Axis in three-dimensional space to ignore in the two-dimensional analysis.

Valid values: 0 or x for the \(x\)-axis, 1 or y for the \(y\)-axis, and 2 or z for the \(z\)-axis.

normstr, keyword-only, default: "rdf"

Determines how the radial histograms are normalized.

Valid values:

norm="rdf": The radial distribution function \(g_{ij}(r)\) is computed.
norm="density": The single particle density \(n_{ij}(r)\) is computed.
norm=None: The raw particle pair count in the radial histogram bins is returned.

exclusionarray-like, keyword-only, optional

Tiles to exclude from the interparticle distances. The groupings parameter dictates what a tile represents.

Shape: \((2,)\).

Example: (1, 1) to exclude self-interactions.

groupingsstr or array-like, keyword-only, default: "atoms"

Determines whether the centers of mass are used in lieu of individual atom positions. If groupings is a str, the same value is used for all groups.

Note

If the desired grouping is not "atoms",

the trajectory file should have segments (or chains) containing residues (or molecules) and residues containing atoms, and
residues and segments should be locally unwrapped at the simulation box edges, if not already, using MDAnalysis.transformations.wrap.unwrap, MDAnalysis.core.groups.AtomGroup.unwrap(), or MDAnalysis.lib.mdamath.make_whole().

Valid values:

"atoms": Atom positions (generally or for coarse-grained simulations).
"residues": Residues’ centers of mass (for atomistic simulations).
"segments": Segments’ centers of mass (for atomistic polymer simulations).

reducedbool, keyword-only, default: False

Specifies whether the data is in reduced units.

n_batchesint, keyword-only, optional

Number of batches to divide the histogram calculation into. This is useful for large systems that cannot be processed in a single pass.

Note

If you use too few bins and too many batches, the histogram counts may be off by a few due to the floating-point nature of the cutoffs. However, when the RDF is averaged over a long trajectory with many particles, the difference should be negligible.

parallelbool, keyword-only, default: False

Determines whether the analysis is performed in parallel.

verbosebool, keyword-only, default: True

Determines whether detailed progress is shown.

**kwargs

Additional keyword arguments to pass to MDAnalysis.analysis.base.AnalysisBase.

Attributes:

universeMDAnalysis.Universe

MDAnalysis.core.universe.Universe object containing all information describing the system.

results.unitsdict

Reference units for the results. For example, to get the reference units for results.bins, call results.units["results.bins"].

results.edgesnumpy.ndarray

Edges of the histogram bins.

Shape: \((N_\mathrm{bins}+1,)\).

Reference unit: \(\textrm{Å}\).

results.binsnumpy.ndarray

Centers of the histogram bins.

Shape: \((N_\mathrm{bins},)\).

Reference unit: \(\textrm{Å}\).

results.countsnumpy.ndarray

Raw particle pair counts in the radial histogram bins.

Shape: \((N_\mathrm{bins},)\).

results.rdfnumpy.ndarray

One of

norm="rdf": the radial distribution function \(g_{ij}(r)\),
norm="density": the single particle density \(n_{ij}(r)\), or
norm=None: the raw particle pair count in the radial histogram bins.

Shape: \((N_\mathrm{bins},)\).

results.coordination_numbersnumpy.ndarray

Coordination numbers \(n_k\). Only available after running calculate_coordination_numbers().

results.pmfnumpy.ndarray

Potential of mean force \(w(r)\). Only available after running calculate_pmf().

Shape: \((N_\mathrm{bins},)\).

Reference unit: \(\mathrm{kJ/mol}\).

results.wavenumbersnumpy.ndarray

Wavenumbers \(q\). Only available after running calculate_structure_factor().

Reference unit: \(\textrm{Å}^{-1}\).

results.ssfnumpy.ndarray

(Partial) static structure factor. Only available after running calculate_structure_factor().

Shape: Same as results.wavenumbers.

Methods

`calculate_coordination_numbers`	Calculates the coordination numbers \(n_k\).
`calculate_pmf`	Calculates the potential of mean force \(w_{ij}(r)\).
`calculate_structure_factor`	Computes the (partial) static structure factor \(S_{ij}(q)\) using the radial histogram bins \(r\) and the radial distribution function \(g_{ij}(r)\) for an isotropic fluid.
`run`	Performs the calculation.
`save`	Saves results to a binary or archive file in NumPy format.

calculate_coordination_numbers(rho: float, *, n_coord_nums: int = 2, threshold: float = 0.1) → None[source]¶

Calculates the coordination numbers \(n_k\).

If the radial distribution function \(g_{ij}(r)\) does not contain \(k\) local minima, this method will return numpy.nan for the coordination numbers that could not be calculated.

Parameters:

rhofloat

Number density \(\rho_j\) of species \(j\).

Reference unit: \(\mathrm{nm}^{-3}\).

n_coord_numsint, keyword-only, default: 2

Number of coordination numbers to calculate.

thresholdfloat, keyword-only, default: 0.1

Minimum \(g_{ij}(r)\) value for a local minimum to be considered the boundary of a radial shell.

calculate_pmf(temperature: float | Quantity | Quantity) → None[source]¶

Calculates the potential of mean force \(w_{ij}(r)\).

Parameters:

temperaturefloat or openmm.unit.Quantity

System temperature \(T\).

Note

If reduced=True was set in the RDF constructor, temperature should be equal to the energy scale. When the Lennard-Jones potential is used, it generally means that \(T^*=1\), or temperature=1.

Reference unit: \(\mathrm{K}\).

calculate_structure_factor(rho: float, x_i: float = None, x_j: float = None, q: ndarray[float] = None, *, q_lower: float = None, q_upper: float = None, n_q: int = 1000, formalism: str = 'FZ') → None[source]¶

Computes the (partial) static structure factor \(S_{ij}(q)\) using the radial histogram bins \(r\) and the radial distribution function \(g_{ij}(r)\) for an isotropic fluid.

Parameters:

rhofloat

Bulk number density \(\rho\).

Reference unit: \(\mathrm{Å}^{-3}\).

x_ifloat, default: 1

Number concentration of species \(i\). Required if the two atom groups are not identical.

x_jfloat, optional

Number concentration of species \(j\). Required if the two atom groups are not identical.

qnumpy.ndarray, optional

Wavenumbers \(q\).

Reference unit: \(\mathrm{Å}^{-1}\).

q_lowerfloat, keyword-only, optional

Lower bound for the wavenumbers \(q\). Has no effect if q is specified.

Reference unit: \(\mathrm{Å}^{-1}\).

q_upperfloat, keyword-only, optional

Upper bound for the wavenumbers \(q\). Has no effect if q is specified.

Reference unit: \(\mathrm{Å}^{-1}\).

n_qint, keyword-only, default: 1_000

Number of wavenumbers \(q\) to generate. Has no effect if q is specified.

formalism :`str`, keyword-only, default: :code:`”FZ”`

Formalism to use for the partial structure factor. Has no effect if the two atom groups are the same.

Valid values:

"general": A general formalism similar to that of the static structure factor, except the second term is multiplied by \(x_ix_j\).
"FZ": Faber–Ziman formalism.
"AL": Ashcroft–Langreth formalism.