biotite.structure.rdf

biotite.structure.rdf(center, atoms, selection=None, interval=(0, 10), bins=100, box=None, periodic=False)[source]

Compute the radial distribution function g(r) (RDF) for one or multiple given central positions based on a given system of particles.

Parameters
centerAtom or AtomArray or AtomArrayStack or ndarray, dtype=float

Coordinates or atoms(s) to use as origin(s) for RDF calculation.

  • If a single Atom or a ndarray with shape (3,) is given, the RDF is only calculated for this position.

  • If an AtomArray or a ndarray with shape (n,3) is given, the calculated RDF histogram is an average over n postions.

  • If an AtomArrayStack or a ndarray with shape (m,n,3) is given, different centers are used for each model m. The calculated RDF histogram is an average over m models and n positions. This requires atoms to be an AtomArrayStack.

atomsAtomArray or AtomArrayStack

The distribution is calculated based on these atoms. When an an AtomArrayStack is provided, the RDF histogram is averaged over all models. Please not that atoms must have an associated box, unless box is set.

selectionndarray, dtype=bool, shape=(n,), optional

Boolean mask for atoms to limit the RDF calculation to specific atoms.

intervaltuple, optional

The range in which the RDF is calculated.

binsint or sequence of scalars or str, optional

Bins for the RDF.

  • If bins is an int, it defines the number of bins for the given interval.

  • If bins is a sequence, it defines the bin edges, ignoring the interval parameter. The output bins has the length of this input parameter reduced by one.

  • If bins is a string, it defines the function used to calculate the bins.

See numpy.histogram() for further details.

boxndarray, shape=(3,3) or shape=(m,3,3), optional

If this parameter is set, the given box is used instead of the box attribute of atoms. Must have shape (3,3) if atoms is an AtomArray or (m,3,3) if atoms is an AtomArrayStack, respectively.

periodicbool, optional

Defines if periodic boundary conditions are taken into account.

Returns
binsndarray, dtype=float, shape=n

The centers of the histogram bins. The length of the array is given by the bins input parameter.

rdfndarry, dtype=float, shape=n

RDF values for every bin.

Notes

Since the RDF depends on the average particle density of the system, this function strictly requires an box.

Examples

Calculate the oxygen-oxygen radial distribution function of water. The range of the histogram starts at 0.2 Å, in order to ignore the counts for the density for each oxygen to itself.

>>> from os.path import join
>>> waterbox = load_structure(join(path_to_structures, "waterbox.gro"))
>>> oxygens = waterbox[:, waterbox.atom_name == 'OW']
>>> bins, g_r = rdf(oxygens, oxygens, bins=49, interval=(0.2, 10), periodic=True)

Print the RDF depending on the radius. Bins are in Å.

>>> for x, y in zip(bins, g_r):
...     print(f"{x:.2f} {y:.2f}")
0.30 0.00
0.50 0.00
0.70 0.04
0.90 0.02
1.10 0.03
1.30 0.06
1.50 0.03
1.70 0.04
1.90 0.04
2.10 0.04
2.30 0.04
2.50 0.16
2.70 1.99
2.90 2.22
3.10 1.34
3.30 1.04
3.50 0.97
3.70 0.94
3.90 0.97
4.10 0.94
4.30 0.98
4.50 0.97
4.70 0.96
4.90 0.99
5.10 0.99
5.30 1.02
5.50 1.02
5.70 0.99
5.90 0.98
6.10 0.98
6.30 0.99
6.50 1.02
6.70 1.02
6.90 1.00
7.10 1.01
7.30 1.01
7.50 1.00
7.70 1.01
7.90 0.99
8.10 0.99
8.30 0.99
8.50 0.99
8.70 0.99
8.90 1.00
9.10 1.01
9.30 1.01
9.50 1.00
9.70 1.00
9.90 0.99

Find the radius for the first solvation shell. In this simple case, the density peak is identified by finding the maximum of the function.

>>> peak_position = np.argmax(g_r)
>>> print(f"{bins[peak_position]/10:.2f} nm")
0.29 nm