spatialstats.paircount

Pair-count correlation functions

Routines to calculate 2-point correlation functions \(G(\mathbf{r}), \ S(\mathbf{q})\) for point and rod-like particle pairs. Can reduce to the usual isotropic \(G(r)\) with the corresponding fourier representation \(S(q)\) and even compute multipoles \(G_{\ell}(r)\) and \(S_{\ell}(q)\) in angular dependencies \(G(r, \hat{\mathbf{r}}\cdot\hat{\mathbf{z}})\) and \(S(q, \hat{\mathbf{q}}\cdot\hat{\mathbf{z}})\).

Adapted from https://github.com/wenyan4work/point_cloud.

corr(positions, boxsize, weights=None, z=1, orientations=None, rmin=None, rmax=None, nr=100, nphi=None, ntheta=None, cos=True, int=<class 'numpy.int32'>, float=<class 'numpy.float64'>, bench=False, **kwargs)

Compute the 2-point correlaton function \(G(\mathbf{r})\) for a set of \(N\) point-like or rod-like particles \(\mathbf{r}_i\) in a 2D or 3D periodic box, where \(\mathbf{r} = (r, \phi, \theta)\) are the spherical coordinates for displacement vectors between particle pairs \(\mathbf{r} = \mathbf{r}_i - \mathbf{r}_j\). \(\phi\) is the azimuthal angle and \(\theta\) is the inclination angle.

If weight = None, \(G(\mathbf{r}) = g(\mathbf{r})\), i.e. the spatial distribution function. This is computed as \(g(\mathbf{r}) = \langle \delta(\mathbf{r}_j - \mathbf{r}_i) \rangle\), where \(\langle\cdot\rangle\) is an average over particle pair displacements \(\mathbf{r}_j - \mathbf{r}_i\) in a periodic box summed over each choice of origin \(\mathbf{r}_i\). The normalization of \(G\) is chosen so that \(g(r)\) decays to 1.

Generally, if the weights argument is a vector defined for all particles \(\mathbf{w}_i\), the pair correlation function is computed as \(G(\mathbf{r}) = \langle (\mathbf{w}_i \cdot \mathbf{w}_j)^z \rangle\), where \(z\) is some exponent.

If particles orientations \(\mathbf{p}_i\) are included, define \((r, \phi, \theta)\) as the rotated coordinate system with \(\mathbf{p}_i\) pointed in the \(+z\) direction.

Note

Reduces to the 1D radial distribution function \(g(r)\) when nphi = None, ntheta = None, and weights = None.

Parameters

positions (np.ndarray, shape (N, ndim)) – Particle positions \(\mathbf{r}_i\) in 2D or 3D for \(N\) particles. Passed to scipy.spatial.cKDTree.
boxsize (list of float) – The rectangular domain over which to apply periodic boundary conditions. Passed to scipy.spatial.cKDTree.
weights (np.ndarray, shape (N, ndim), optional) – Particle vectors \(\mathbf{w}_i\) over which to calculate pair correlation function.
z (int) – Exponent in averaging \(\langle (\mathbf{w}_i \cdot \mathbf{w}_j)^z \rangle\).
orientations (np.ndarray, shape (N, ndim), optional) – Particle orientation vectors \(\mathbf{p}_i\). Vectors should be unitary, but they will be normalized automatically.
rmin (float, optional) – Minimum \(r\) value in \(g(r, \phi, \theta)\).
rmax (float, optional) – Cutoff radius for KDTree search and maximum \(r\) value in \(g(r, \phi, \theta)\). Default is half the maximum dimension of boxsize.
nr (int, optional) – Number of points to bin in \(r\).
nphi (int, optional) – Number of points to bin in \(\phi\).
ntheta (int, optional) – Number of points to bin in \(\cos\theta\) or \(\theta\).
cos (bool, optional) – Choose whether to bin in slices of \(\cos\theta\) or \(\theta\). Default is \(\cos\theta\).
int (np.dtype, optional) – Integer type for pair counting array. Lets the user relax memory requirements.
float (np.dtype, optional) – Floating-point type for displacement buffers. Lets the user relax memory requirements.
bench (bool, optional) – Print message for time of calculation.

Returns

g (np.ndarray, shape (nr, nphi, ntheta)) – Radial distribution function \(G(r, \phi, \cos\theta)\). If the user does not bin for a certain coordinate, \(G\) will not be 3 dimensional (e.g. if nphi = None, \(G\) will be shape (nr, ntheta)).
r (np.ndarray, shape (nr,)) – Left edges of radial bins \(r\).
phi (np.ndarray, shape (nphi,)) – Left edges of angular bins \(\phi \in [-\pi, \pi)\).
theta (np.ndarray, shape (ntheta,)) – Left edges of angular bins \(\theta \in [0, \pi)\) or \(\cos\theta \in [-1, 1)\). Not returned for 2D datasets.
vol (np.ndarray, shape g.shape) – Bin volume used to normalize \(G\).

fourier_multipoles(g_ell, poles, r, N, boxsize, dq=None, nq=None)

Compute fourier space multipoles \(S_{\ell}(q)\) from real space multipoles \(G_{\ell}(r)\) using the transform

\[S_{\ell}(q) = 4\pi\rho (-i)^{\ell} \int dr \ r^2 j_{\ell}(qr) G_{\ell}(r).\]

in 3D and

\[S_{\ell}(q) = 2\pi\rho (-i)^{\ell} \int dr \ r J_{\ell}(qr) G_{\ell}(r)\]

in 2D. \(\rho = N / V\) is the density. Note that \(G_{\ell}\) must decay to zero for this to be well-defined.

Normalization is chosen based on the isotropic structure factor \(s(q)\) defined by

\[s(q) = 1 + 4\pi\rho \int dr \ r^2 j_0(qr) [g(r) - 1]\]

, where \(g(r)\) is the radial distribution function.

Parameters

g_ell (np.ndarray, shape (nr,)) – Real space multipoles \(G_{\ell}(r)\) returned by spatialstats.particles.multipoles. Can be a list if poles is too.
poles (int) – Multipole indices of g_ell.
r (np.ndarray, shape (nr,)) – Radial bins \(r\) to integrate over.
N (int) – The number of particles \(N\).
boxsize (list of float) – The rectangular domain over which to apply periodic boundary conditions. Simply used to calculate \(\rho = N/V\).
dq (float) – Spacing of \(q\) bins. Default is the fundamental mode given by dq = np.pi / r.max().
nq (int) – Number of \(q\) bins. Default is r.size.

Returns

s_ells (np.ndarray, shape (nq,)) – Multipoles \(S_{\ell}(q)\). Returns a list if poles is a list.
q (np.ndarray, shape (nq,)) – Wavenumber bins \(q\).

multipoles(g, costheta, poles=0)

Compute multipoles of a correlation function of the form \(G(r, \cos\theta)\) defined by

\[G_{\ell}(r) = \frac{2\ell+1}{2} \int_{-1}^1 d\cos\theta \ G(r, \cos\theta)\mathcal{L}_{\ell}(\cos\theta).\]

Parameters

g (np.ndarray, shape (nr, ntheta)) – 2D correlation function \(G(r, \cos\theta)\) returned by spatialstats.particles.corr.
costheta (np.ndarray, shape (ntheta,)) – Array of \(\cos\theta\) bins used as integration domain in scipy.integrate.trapezoid.
ell (int or list) – Degree of multipole to compute. For example, poles = [0, 1, 2] computes the monopole, dipole, and quadrapole.

Returns

g_ells – Multipoles \(G_{\ell}(r)\). If type(poles) = list, then g_ells is a list.

Return type

np.ndarray, shape (nr,)