Radial Distribution Function

The radial distribution function measures the probability of finding a target atom at distance \(r\) from a reference atom relative to an ideal gas at the same effective target density [Hansen2013]. For histogram bin \(i\), PQAnalysis uses the standard simulation estimator [Allen2017],

\[g_i = \frac{H_i}{\rho_T N_R N_F \Delta V_i},\]

where \(H_i\) is the eligible pair count, \(\rho_T\) the target number density, \(N_R\) the number of reference atoms, \(N_F\) the number of frames and \(\Delta V_i\) the spherical-shell volume. Peaks mark preferred pair separations, minima separate coordination shells, and \(g(r) \approx 1\) is uncorrelated bulk-like pair density. The analytic schematic below (not simulation output) shades the first coordination shell up to the first minimum; the lower panel is the running coordination number \(N(r) = 4\pi\rho\int_0^r g(s)\,s^2\,\mathrm{d}s\) for \(\rho = 0.0334\) Å⁻³.

Input

traj_files = trajectory.xyz
reference_selection = O
target_selection = H
delta_r = 0.05
r_max = 8.0
out_file = rdf.dat

Saved as rdf.in, it runs with:

$ pqanalysis rdf rdf.in

The keys above are typical for a bulk trajectory. The bundled Example data fixture uses delta_r = 0.5 and r_max = 4.0 because it is one molecule in a 10 Å box. no_intra_molecular = True drops pairs inside the same molecule and needs restart_file and moldescriptor_file; naming both files enables it by itself. PQAnalysis infers the usual PQ companion filenames when they sit beside the trajectory. Selection strings are described in Atom selections; the full key table is on RDFInputFileReader.

Output

out_file has five columns: bin center, \(g(r)\), running coordination number, shell population and pair-count residual. Exact definitions are in RDF. The coordination number is read off at the row where \(g\) has its first minimum.

Python

from PQAnalysis.analysis import RDF
from PQAnalysis.io import TrajectoryReader

r, g, n, shell, residual = RDF(
    TrajectoryReader("examples/water/trajectory.xyz"),
    reference_species="O",
    target_species="H",
    delta_r=0.5,
    r_max=4.0,
).run()

rdf() runs an input file instead.

Before you trust it

• r_max above half the shortest box edge is clamped; on triclinic cells set it to half the smallest perpendicular width yourself.

• Bins near the first peak need about \(10^4\) pair counts for 1 % noise. Add frames rather than widening delta_r.

• The plateau at large \(r\) should sit at 1. If it does not, look at the volume distribution (NPT) or the selections.

• Intramolecular pairs are included unless no_intra_molecular is on; naming both topology files turns it on.

Each point is derived in RDF: Theory and Validity.