calc module

collRC is all about calculating and representing reaction coordinates from simulated colloidal ensembles. This module contains all of the framework to build up the various calculations needed to compute robust reaction coordinates:

The locality module processes particle positions to compute local environment metrics like nearest neighbors in 2D cartesian, 2D stretched, 2D projected, and 3D cartesian spaces.

The morphology module computes information like local density and gyration tensors of ensembles and ensemble subsets.

The orient_order module computes complex-valued local and global p-atic orientational order.

The bond_order module computes complex-valued bond orientational order in 2D cartesian, 2D stretched, 2D projected, and 3D cartesian spaces (via steinhardt calculations).

The clusters module computes identifies spatially correlated groups of particles and computes statistics on their sizes and shapes.

locality submodule

Contains many helper methods to compute geometric properties of particle configurations, including neighbor lists, handling periodic boundary conditions, and local tangent plane vectors on curved surfaces.

calc.locality.neighbors(pts: ndarray, neighbor_cutoff: float | None = None, num_closest: int | None = None) → ndarray

Determines neighbors in a configuration of particles based on a cutoff distance:

\[n_{jk} = \delta r_{jk} < r_{cut}\]

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions.

• neighbor_cutoff (scalar, optional) – specify the distance which defines neighbors. Defaults to halfway between the first coordination peaks for a perfect sixfold crystal.

• num_closest (int, optional) – specify the maximum number of neighbors (within the cutoff) per particle. a.k.a pick that many of the closest neighbors per particle.

Returns:

\([N,N]\) boolean array indicating which particles are neighbors

Return type:

ndarray

calc.locality.quat_to_angle(quat: ndarray) → ndarray

Finds the angle \(\theta\) about the z-axis from a quaternion representation in 2D: \(q = \cos(\theta/2) + \sin(\theta/2)\mathbf{k}\)

Parameters:

quat (ndarray) – an \([N,4]\) array of quaterions encoding particle orientation

Returns:

an \([N,]\) array of quadrant-corrected 2d angular orientations of the particles

Return type:

ndarray

calc.locality.stretched_neighbors(pts: ndarray, angles: ndarray, rx: float = 1.0, ry: float = 1.0, neighbor_cutoff: float = 2.6) → ndarray

Determines neighbors in a configuration of anisotropic particles based on a cutoff distance in the rotated/stretched frame of each particle using the equaitons defined in (Torrez-Diaz Soft Matter, 2022):

\[n_{jk} = \sqrt{\big(\mathbf{\delta r_jk} \cdot \hat{\mathbf{x_j}}/r_x\big)^2 + \big(\mathbf{\delta r_jk} \cdot \hat{\mathbf{y_j}}/r_y\big)^2} < n_{cut}\]

Where \(\hat{\mathbf{x_j}} = \cos(\theta_j)\hat{\mathbf{x}} + \sin(\theta_j)\hat{\mathbf{y}}\) and \(\hat{\mathbf{y_j}} = -\sin(\theta_j)\hat{\mathbf{x}} + \cos(\theta_j)\hat{\mathbf{y}}\) are the local unit vectors along the long (\(r_x\)) and short (\(r_y\)) axes of particle \(j\) rotated to it’s orientation \(\theta_j\), respectively.

Parameters:

• pts (ndarray) – an \([N,d]\) array of the positions each anisotropic particle in the configuration

• angles (ndarray) – an \([N,]\) array of the orientation of each anisotropic particle in the configuration

• rx (scalar, optional) – the radius of the long axis of the particle (insphere radius times aspect ratio), defaults to 1.0

• ry (scalar, optional) – the radius of the short axis of the partice (i.e. insphere radius), defaults to 1.0

• neighbor_cutoff (scalar, optional) – specify the dimemsionless stretched distance which defines neighbors. Defaults to 2.6.

Returns:

\([N,N]\) boolean array indicating which particles are neighbors

Return type:

ndarray

calc.locality.matrix_to_box(box: ndarray) → ndarray

\[\begin{split}\begin{bmatrix} L_x & xy*L_y & xz*L_z \\ 0 & L_y & yz*L_z \\ 0 & 0 & L_z \end{bmatrix} \rightarrow [L_x,L_y,L_z,xy,xz,yz]\end{split}\]

Parameters:

box (ndarray) – a matrix containing the basis vectors of a bounding box

Returns:

a length 6 list of box parameters

Return type:

ndarray

calc.locality.box_to_matrix(box: list) → ndarray

\[\begin{split}[L_x,L_y,L_z,xy,xz,yz] \rightarrow \begin{bmatrix} L_x & xy*L_y & xz*L_z \\ 0 & L_y & yz*L_z \\ 0 & 0 & L_z \end{bmatrix}\end{split}\]

Parameters:

box (array_like) – a length 6 list of box parameters

Returns:

a matrix containing the basis vectors for the equivalent hoomd box

Return type:

ndarray

calc.locality.expand_around_pbc(coords: ndarray, basis: ndarray, padfrac: float = 0.8, max_dist: float = None) → tuple[ndarray, ndarray]

Given a frame and a box basis matrix, returns a larger frame which includes surrounding particles from the nearest images, as well as the index relating padded particles back to their original image. This will enable methods like scipy.voronoi() to respect periodic boundary conditions.

Parameters:

• coords (ndarray) – a \([N,d]\) array of particle coordinates in d-dimensions

• basis (ndarray) – a \([d,d]\) matrix of basis vectors for the simulation box. A 2D box should have basis[2,2]=0, otherwise the 3D case is assumed and the function prepares 27 periodic images (rather than only 9 in 2D).

• padfrac (float, optional) – the number of extra particles, as a fraction of the total number, to include in the ‘pad’ of surrounding particles, defaults to 0.8

• max_dist (float, optional) – alternatively to padfrac, specify a maximum distance from original particles to include surrounding particles. defaults to None

Returns:

a \([N', d]\) array of particle coordinates in d-dimensions which respect periodic boundary conditions around the central N particles, as well as a \([N',]\) array of indices relating padded particles back to their original image

Return type:

ndarray, ndarray

calc.locality.padded_neighbors(pts: ndarray, basis: ndarray, neighbor_cutoff: float = np.float64(1.3660254037844386), padfrac: float = None) → ndarray

Determines neighbors in a configuration of particles based on a cutoff distance while respecting the periodic boundary condition using expand_around_pbc():

\[n_{jk} = \delta r_{jk} < r_{cut} \big|\big| \delta r_{jk'} < r_{cut}\]

For particles \(k'\) which are periodic images of particle \(k\).

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions.

• basis (ndarray) – a \([d,d]\) matrix of basis vectors for the simulation box

• neighbor_cutoff (scalar, optional) – specify the distance which defines neighbors. Defaults to halfway between the first coordination peaks for a perfect sixfold crystal.

• padfrac – the number of extra particles, as a fraction of the total number, to include in the ‘pad’ of surrounding particles, This may be computatinally faster than using expand_around_pbc() with max_dist. Defaults to None

Returns:

\([N,N]\) boolean array indicating which particles are neighbors

Return type:

ndarray

calc.locality.local_vectors(pts: ndarray, gradient: callable, ref: ndarray = array([0, -1, 0]))

Computes two orthogonal unit vectors tangent to the local surface at each point \(\mathbf{r}_j\):

\[\begin{split}\mathbf{e}_{1,j} = \frac{\nabla f(\mathbf{r}_j) \times \hat{\mathbf{\gamma}}}{|\nabla f(\mathbf{r}_j) \times \hat{\mathbf{\gamma}}|} \\ \mathbf{e}_{2,j} = \frac{\nabla f(\mathbf{r}_j) \times \mathbf{e}_1}{|\nabla f(\mathbf{r}_j) \times \mathbf{e}_1|}\end{split}\]

Given the implicit function \(f(x,y,z)=0\) defines the surface, its gradient \(\nabla f\) defines the normal vector, and \(\hat{\mathbf{\gamma}}\) is an arbitrary reference unit vector.

Parameters:

• pts (ndarray) – an \([N,3]\) array of positions at which to compute local tangent vectors

• gradient (callable) – a function which computes normal vector to a surface

• ref (ndarray, optional) – a reference vector to help define the local frame, defaults to \((0,-1,0)\) (so that \(e_1=\hat{x}\) usually)

Returns:

an \([N,2,3]\) array of two orthogonal unit vectors tangent to the local surface defined by the gradient at each x

Return type:

ndarray

calc.locality.tangent_connection(pts: ndarray, gradient: callable, ref: ndarray = array([0, -1, 0])) → ndarray

Computes the complex connection between local tangent planes at each pair of points. This factor takes the form

\[R_{jk} = e^{i \theta_{jk}} = \exp\bigg[i\tan^{-1}\bigg(\frac{e_{1,j} \cdot e_{2,k}}{e_{1,j} \cdot e_{1,k}}\bigg)\bigg]\]

where the local tangent vectors \(e_{1,j}\) and \(e_{2,j}\) are computed using local_vectors(). \(R_{jk}\) describes how to rotate complex numbers defined relative to \(e_{1,j}\) into those defined relative to \(e_{1,k}\). Due to parallel transport on curved surfaces, this factor is only really useful when \(j\) and \(k\) are nearby points on the surface, and is more and more approximate as the distance between points increases.

Parameters:

• pts (ndarray) – an \([N,3]\) array of positions at which to compute local tangent vectors

• gradient (callable) – a function which computes normal vector to a surface

• ref (ndarray, optional) – a reference vector to help define the local frame, defaults to \(-\hat{\mathbf{y}}\) so that \(e_1=\hat{\mathbf{x}}\) usually

Returns:

an \([N,N]\) complex representation of the connection between local tangent planes at each pair of points

Return type:

ndarray[complex]

morphology submodule

Contains methods to calculate morphological properties of particle ensembles and subsets of ensembles, including local area fraction, gyration tensors, and functions of their eigenvalues like gyration radius, acylindricity, asphericity and anisotropy.

calc.morphology.central_eta(pts: ndarray, box: list, ptcl_area: float = 0.7853981633974483, nbins: int = 3, bin_width: float = 4.0, jac: str = 'x')

Computes the average area fraction in the central region of a configuration of particles, accounting for the jacobian of the coordinate system:

\[\eta = \bigg\langle\frac{N\cdot A_p}{A_{bin}}\bigg\rangle_{{bins}}\]

Parameters:

• pts (ndarray) – an \([N,d]\) array of particle positions in any dimensions (though only the first two are used)

• box (array-like) – a list defining the simulation box in the gsd convention: [Lx, Ly, Lz, xy, xz, yz]

• ptcl_area (float, optional) – the area of a single particle, defaults to π/4 (for a circle of diameter 1.0)

• nbins (int, optional) – the number of histogram bins to average over in the center of the configuration, defaults to 3

• bin_width (float, optional) – the width of each histogram bin, defaults to 4

• jac (str, optional) – the type of jacobian to account for when computing area fraction. Options are ‘x’ (linear in x), ‘y’ (linear in y), and ‘r’ (radial). Defaults to ‘x’.

Returns:

the average area fraction in the central region of the configuration

Return type:

scalar

calc.morphology.gyration_radius(pts: ndarray) → float

Returns the radius of gyration according to the formula

\[R_g^2 \equiv \frac{1}{N}\sum_k|\mathbf{r}_k-\bar{\mathbf{r}}|^2 = \frac{1}{N^2}\sum_{k>j}|\mathbf{r}_j - \mathbf{r}_k|^2 \]

where the \(k>j\) in the summation index indicates that repeated pairs are not summed over

Parameters:

pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions

Returns:

the radius of gyration of the particles about their center of mass

Return type:

scalar

calc.morphology.gyration_tensor(pts: ndarray, ref: ndarray = None, weights: ndarray = None) → ndarray

Returns the gyration tensor (for principal moments analysis) of an ensemble of particles according to the formula

\[S_{mn} = \frac{1}{N}\sum_j r^{(j)}_m r^{(j)}_n\]

where the positions, \(\mathbf{r}^{(j)} = (r^{(j)}_1, r^{(j)}_2, \ldots, r^{(j)}_d)\), are defined in their center of mass reference frame. Alternatively users can pass in a weights array to compute a weighted gyration tensor:

\[S_{mn} = \frac{1}{\sum_j w^{(j)}}\sum_j w^{(j)}r^{(j)}_m r^{(j)}_n\]

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions,

• ref (ndarray , optional) – point in d-dimensional space from which to reference particle positions, defaults to the mean position of the points. Use this for constraining the center of mass to the surface of a manifold, for instance.

• weights (ndarray , optional) – \([N,]\) array of weights for each particle, defaults to None (equal weighting)

Returns:

the \([d,d]\) gyration tensor of the ensemble

Return type:

ndarray

calc.morphology.acylindricity(pts: ndarray = None, gyr: ndarray = None) → float

Returns the acylindricity of a set of points, defined as a function of the gyration tensor eigenvalues:

\[c \equiv \lambda_y^2 - \lambda_x^2\]

Where \(\lambda_x^2\leq\lambda_y^2\leq\lambda_z^2\) are the eigenvalues of the gyration tensor.

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions

• gyr (ndarray) – \([d,d]\) gyration tensor of the ensemble, optional

Returns:

the acylindricity of the particles

Return type:

scalar

calc.morphology.asphericity(pts: ndarray = None, gyr: ndarray = None) → float

Returns the asphericity of a set of points, defined as a function of the gyration tensor eigenvalues:

\[b \equiv \lambda_z^2 - \frac{{1}}{{2}}\big(\lambda_y^2 + \lambda_x^2\big)\]

Where \(\lambda_x^2\leq\lambda_y^2\leq\lambda_z^2\) are the eigenvalues of the gyration tensor.

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions

• gyr (ndarray) – \([d,d]\) gyration tensor of the ensemble, optional

Returns:

the asphericity of the particles

Return type:

scalar

calc.morphology.shape_anisotropy(pts: ndarray = None, gyr: ndarray = None) → float

Returns the anisotropy of a set of points, defined as a function of the gyration tensor eigenvalues:

\[\kappa^2 \equiv \frac{{3}}{{2}}\frac{{\lambda_z^4 + \lambda_y^4 + \lambda_x^4}}{{(\lambda_x^2 + \lambda_y^2 + \lambda_z^2)^2}} - \frac{{1}}{{2}} = \frac{{b^2 + (3/4)c^2}}{{R_g^4}}\]

Where \(\lambda_x^2\leq\lambda_y^2\leq\lambda_z^2\) are the eigenvalues of the gyration tensor.

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions

• gyr (ndarray) – \([d,d]\) gyration tensor of the ensemble, optional

Returns:

the anisotropy of the particles

Return type:

scalar

calc.morphology.circularity(pts: ndarray = None, gyr: ndarray = None) → float

Returns the ‘circularity’ of a 2D colloidal cluster as used in Zhang, Sci. Adv. 2020:

\[c = 1 - \frac{\lambda_x-\lambda_y}{\sqrt{\lambda_x^2 + \lambda_y^2}} \]

Where \(\lambda_x^2\geq\lambda_y^2\geq\lambda_z^2\simeq0\) are the eigenvalues of the gyration tensor (notably ordered opposite from the other methods defined above). When the cluster is circular the two principal moments are equal and so \(c=1\). When the cluster is a linear chain, the smaller principal moment approaches zero, and so \(c\to0\).

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions

• gyr (ndarray) – \([d,d]\) gyration tensor of the ensemble, optional

Returns:

the complex circularity of the particles

Return type:

scalar

calc.morphology.ellipticity(pts: ndarray = None, gyr: ndarray = None, ref: ndarray = array([1, 0, 0])) → complexfloating

Here we expand on the ‘circularity’ of a colloidal cluster as used in Zhang, Sci. Adv. 2020. This metric is calcuated using the principal moments of the gyration_tensor(). For 2D ensembles, after diagonalization:

\[\begin{split}S_{mn} \to \Lambda_{mn} = \begin{pmatrix} \lambda_x^2 & 0 & 0 \\ 0 & \lambda_y^2 & 0 \\ 0 & 0 & \lambda_z^2\simeq0 \end{pmatrix}\end{split}\]

With the prinicipal moments defined \(\lambda_x^2\geq\lambda_y^2\geq\lambda_z^2\). Using these principal moments and the major axis (largest eigenvector) \(\mathbf{e}_x\), we can define the circularity magnitude and orientation of the cluster relative to a reference vector, \(\hat{\mathbf{\gamma}}\):

\[1 - c = \frac{\lambda_x-\lambda_y}{\sqrt{\lambda_x^2 + \lambda_y^2}} \quad \phi = \arccos(\mathbf{e}_x \cdot \hat{\mathbf{\gamma}})\]

When the cluster is circular the two principal moments are equal and so \(1-c=0\), and the angle \(\phi\) is irrelevant. When the cluster is a linear chain, the smaller principal moment approaches zero, and so \(1-c\to1\) and the angle \(\phi\) points along the line. We can summarily express this as a complex number with 2-fold angular symmetry:

\[\varepsilon = (1-c) e^{2i\phi}\]

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions

• gyr (ndarray) – \([d,d]\) gyration tensor of the ensemble, optional

• ref (ndarray , optional) – reference direction (unit vector) in d-dimensional space to which to measure the major axis orientation, defaults to x-axis

Returns:

the complex ellipticity of the particles

Return type:

complex

orient_order submodule

Contains methods to calculate orientational order, meaning the order and structure inherent to to the orientations of the particles themselves. These methods generally cover nematic, tetract and hexatic order parameters, which can be generalized as ‘p-atic’ order.

\[S_2 \equiv \xi_2 \qquad T_4 \equiv \xi_4 \qquad H_6 \equiv \xi_6\]

calc.orient_order.global_patic(angles: ndarray, p: int = 2) → complexfloating

Computes the global p-atic order of a particle ensemble. Traditionally this is a scalar-valued order parameter which characterizes how aligned particles are with a director orientation:

\[|\xi_p| = \max_{\theta_p}\langle\cos(p(\theta_j-\theta_p))\rangle_j\]

But, it’s more convienient both to compute and represent this as a complex-valued quantity which also encodes the director angle:

\[\xi_p = \langle e^{i p \theta_j}\rangle = |\xi_p|e^{i p \theta_p}.\]

Parameters:

• angles (ndarray) – an \([N,]\) array of orientations of each particle in the frame

• p (int) – the p-atic order to compute (e.g. 2 for nematic, 4 for tetratic, 6 for hexatic)

Returns:

the complex value of the global p-atic order and the angle which defines it

Return type:

complex

calc.orient_order.local_patic(angles: ndarray, nei_bool: ndarray, p: int = 2) → ndarray

Computes the local p-atic order, per particle, of an ensemble. Quantifies the local orientational order of a system by calculating a director for particle \(j\) via its neighboring particle(s) \(k\), as defined by (Baron J. Chem. Phys. 2023).

\[|\xi_{p,j}| = \max_{\theta_{p,j}}\langle\cos(p(\theta_k-\theta_{p,j}))\rangle_{r_{jk}<6a_x}\]

Or, more succinctly:

\[\xi_{p,j} = \langle e^{ip\theta_k}\rangle_{j(r_{jk}<6a_x)} = |\xi_{p,j}|e^{i p \theta_{p,j}}\]

Parameters:

• angles (ndarray) – an \([N,]\) array of orientations of each particle in the frame

• nei_bool (ndarray) – a \([N,N]\) boolean array defining particle neighbors

• p (int, optional) – the p-atic order to compute (e.g. 2 for nematic, 4 for tetratic, 6 for hexatic), defaults to 2

Returns:

\([N,]\) array of the local p-atic parameter encoded with the associated angle of the director

Return type:

ndarray[complex]

bond_order submodule

Contains methods to calculate bond orientational order. First in 2D cartesian space, then in 2D stretched space for anisotropic particles, then in 2D spaces projected onto curved surfaces, and finally 3D cartesian spaces using Steinhardt calculations.

calc.bond_order.flat_bond_order(pts: ndarray, nei_bool: ndarray, order: int = 6, ret_global=False) → tuple[ndarray, complexfloating]

Calculates the local and global bond orientational order parameter of each particle in a 2D configuration with respect to the y axis. The local n-fold bond orientaitonal order for a particle \(j\) is:

\[\psi_{n,j} = \frac{1}{N_j}\sum_k\psi_{n,jk}=\frac{1}{N_j}\sum_ke^{in\theta_{jk}}\]

Where the sum is over all \(N_j\) neighboring particles and \(\theta_{jk}\) is the angle between particles \(j\) and \(k\). The global n-fold bond orientational order is:

\[\psi_{n,g} = \langle e^{in\theta_{jk}} \rangle_{jk} \simeq \langle \psi_{n,j} \rangle_{j}\]

Where the mean is over all unique bonds between particle pairs \(j\) and \(k\), or approximately the mean over all particles’ local bond order parameters.

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions, though the calculation only accesses the first two dimensions.

• nei_bool (ndarray) – a \([N,N]\) boolean array indicating neighboring particles.

• order (int, optional) – n-fold order defines the argument of the complex number used to calculate psi_n, defaults to 6

• ret_global (bool, optional) – whether to return the global bond order parameter, defaults to False

Returns:

\([N,]\) array of complex bond orientational order parameters, and (if ret_global) the global bond orientational order.

Return type:

ndarray[complex] (, complex)

calc.bond_order.stretched_bond_order(pts: ndarray, angles: ndarray, nei_bool: ndarray, rx: float = 1.0, ry: float = 1.0, order: int = 6, ret_global=False) → tuple[ndarray, complexfloating]

Computes the local and global stretched bond orientational order parameter. This calculation rotates coordinates into a frame of reference stretched according to the long and short axes of each particle according to equations given in (Torrez-Diaz Soft Matter, 2022):

\[\psi_{n,j} = \frac{1}{N_j}\sum_k\psi_{n,jk}=\frac{1}{N_j}\sum_ke^{in\theta^s_{jk}}e^{in\theta_j}\]

Where the sum is over all \(N_j\) neighboring particles, \(\theta_j\) is the orientaion of particle \(j\), and \(\theta^s_{jk}\) is the angle between particles \(j\) and \(k\) in the stretched coordiante system of particle \(j\):

\[\theta^s_{jk} = \tan^{-1}\bigg[\frac{\mathbf{\delta r_{jk}} \cdot \hat{\mathbf{y_j}}/r_y}{\mathbf{\delta r_{jk}} \cdot \hat{\mathbf{x_j}}/r_x}\bigg]\]

Where \(\hat{\mathbf{x_j}} = \cos(\theta_j)\hat{\mathbf{x}} + \sin(\theta_j)\hat{\mathbf{y}}\) and \(\hat{\mathbf{y_j}} = -\sin(\theta_j)\hat{\mathbf{x}} + \cos(\theta_j)\hat{\mathbf{y}}\) are the local unit vectors along the long (\(r_x\)) and short (\(r_y\)) axes of particle \(j\), respectively. The global n-fold stretched bond orientational order is:

\[\psi_{n,g} = \langle e^{in\theta^s_{jk}}e^{in\theta_j} \rangle_{jk} \simeq \langle \psi_{n,j} \rangle_j\]

Where the mean is over all unique bonds between particles \(j\) and \(k\), or approximately the mean over all particles’ local bond order parameters.

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions, though the calculation only access the first two dimensions.

• angles (ndarray) – the orientation of each anisotropic particle in the configuration

• nei_bool (ndarray) – a \([N,N]\) boolean array indicating neighboring particles.

• rx (scalar, optional) – the radius of the long axis of the particle (insphere radius times aspect ratio), defaults to 1.0

• ry (scalar, optional) – the radius of the short axis of the partice (i.e. insphere radius), defaults to 1.0

• order (int, optional) – n-fold order defines the argument of the complex number used to calculate psi_n, defaults to 6

• ret_global (bool, optional) – whether to return the global bond order parameter, defaults to False

Returns:

\([N,]\) array of complex bond orientational order parameters, and (if ret_global) the global bond orientational order.

Return type:

ndarray[complex] (, complex)

calc.bond_order.projected_bond_order(pts: ndarray, gradient: callable, nei_bool: ndarray, order: int = 6, ref: ndarray = array([0, -1, 0])) → tuple[ndarray, float]

Computes the local bond orientational order parameter of each particle in a 2D configuration projected onto the local tangent plane of a curved surface. The local n-fold bond orientaitonal order for a particle \(j\) is:

\[\psi_{n,j} = \frac{1}{N_j}\sum_k\psi_{n,jk}=\frac{1}{N_j}\sum_ke^{in\theta_{jk}}\]

Where the sum is over all \(N_j\) neighboring particles and \(\theta_{jk}\) is the angle between particles \(j\) and \(k\) projected onto the local tangent plane of particle \(j\), defined by a basis \(\{\hat{\mathbf{e}}_{1,j},\hat{\mathbf{e}}_{2,j}\}\) computed using the provided gradient function and the local_vectors() method.

\[\theta_{jk} = \tan^{-1}\big[ (\mathbf{r}_{jk}\cdot\hat{\mathbf{e}}_{2,j})\big/(\mathbf{r}_{jk}\cdot\hat{\mathbf{e}}_{1,j}) \big]\]

Due to parallel transport on curved surfaces, there is no well-defined global bond orientational order because of the lack of a common reference vector.

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions.

• gradient (callable) – a callable function that computes the local gradient at a given point

• nei_bool (ndarray) – a \([N,N]\) boolean array indicating neighboring particles.

• order (int, optional) – n-fold order defines the argument of the complex number used to calculate psi_n, defaults to 6

• ref (ndarray, optional) – a reference vector to help define the local frame, defaults to np.array([0,-1,0])

Returns:

\([N,]\) array of complex bond orientational order parameters, and the norm of their mean.

Return type:

ndarray[complex]

calc.bond_order.steinhardt_bond_order(pts: ndarray, nei_bool: ndarray, l: int = 6, ret_global: bool = False) → tuple[ndarray, complexfloating]

Calculates the local and global bond orientational Steinhardt order parameter of each particle in a 3D configuration using a superposition of complex spherical harmonics. The local l-fold bond orientaitonal order for a particle \(j\) is a list of \(2l+1\) complex-valued spherical harmonic coefficients:

\[\big\{q_{lm,j}\big\} = \big\{\frac{1}{N_j}\sum_kY_{lm}(\theta_{jk},\phi_{jk})\big\}\]

Where the sum is over all \(N_j\) neighboring particles and \(\theta_{jk}, \phi_{jk}\) are the angles (in spherical coordinates) between particles \(j\) and \(k\). Each particle has a rotationally invariant local bond order magnitude given by:

\[q_{l,j} = \sqrt{\frac{4\pi}{2l+1} \sum_m |q_{lm,j}|^2}\]

Which allows us to similarly define a rotationally-invariant global l-fold bond orientational order:

\[q_{l} = \sqrt{\frac{4\pi}{2l+1} \sum_m |\langle q_{lm,j}\rangle_j|^2}\]

Where we average each particle’s \(q_{lm}\) before summing their magnitudes over \(m\) to make it rotationally invariant.

Parameters:

• pts (ndarray) – \([N,d]\) array of particle positions in ‘d’ dimensions, though the calculation only access the first three dimensions.

• nei_bool (ndarray) – a \([N,N]\) boolean array indicating neighboring particles.

• l (int, optional) – The degree of the spherical harmonics used to calculate \(l\)-fold bond order parameter, defaults to 6

• ret_global (bool, optional) – whether to return the global bond order parameter, defaults to False

Returns:

\([N,2l+1]\) array of complex bond orientational order parameters, and (if ret_global) the global bond orientational order.

Return type:

ndarray[complex] (, complex)

calc.bond_order.crystal_connectivity(psis: ndarray, nei_bool: ndarray, crystallinity_threshold: float = 0.32, norm: float | None = 6, phase_rotate: ndarray | None = None, calc_3d=False) → ndarray

Computes the crystal connectivity of each particle in a 2D cartesian, 2D stretched, 2D projected, or 3D system. The crystal connectivity measures the similarity of bond-orientational order parameter over all pairs of neighboring particles in order to determine which particles are part of a definitive crystalline domain. In flat 2D the crystal connectivity of particle \(j\) is given by:

\[C_{n,j} = \frac{1}{N_C}\sum_k^{\text{nei}}\bigg[ \frac{\text{Re}\big[\psi_{n,j}\psi_{n,k}^*\big]}{|\psi_{n,j}\psi_{n,k}^*|} \geq \Theta_C \bigg]\]

Where the \(\psi\)’s are n-fold bond orientational order parameters for each particle, \(\Theta_{C}\) is a ‘crystallinity threshold’ used to determine whether two neighboring particles are part of the same crystalline domain, and \(N_C\) is a factor used to simply normalize \(C_{n,j}\) between zero and one.

On curved surfaces neighboring particles may have bond orientational order parameters defined in different local tangent planes. In this case the crystalline connectivity becomes:

\[C_{n,j} = \frac{1}{N_C}\sum_k^{\text{nei}}\bigg[ \frac{\text{Re}\big[\psi_{n,j}((R_{jk})^n\psi_{n,k})^*\big]}{|\psi_{n,j}((R_{jk})^n\psi_{n,k})^*|} \geq \Theta_C \bigg]\]

Where \(R_{jk}\) encodes the rotation between neighboring tangent planes, i.e. the output of tangent_connection(). \((R_{jk})^n\) is then used to rotate the n-fold bond orientational order of neighboring particles.

Optionally, users may pass in an \([N,2l+1]\) array of complex spherical harmonic Steinhardt order parameters (and pass calc_3d=True) for use in 3D systems. In this case, the inner product is computed over all \(2l+1\) components of the Steinhardt order parameter rather than a single complex bond-OP:

\[C_{l,j} = \frac{1}{N_C}\sum_k^{\text{nei}}\bigg[ \frac{\text{Re}\big[\sum_{m=-l}^{l}q_{lm,j}q_{lm,k}^*\big]}{\sqrt{\sum_m|q_{lm,j}|^2} \sqrt{\sum_m|q_{lm,k}|^2}} \geq \Theta_C \bigg]\]

Parameters:

• psis (ndarray) – \([N,]\) array of complex bond orientational order parameters. Optionally, an \([N,2l+1]\) array of complex spherical harmonic Steinhardt order parameters for use in 3D systems.

• nei_bool (ndarray) – a \([N,N]\) array indicating neighboring particles.

• crystallinity_threshold (float, optional) – the minimum inner product of adjacent complex bond-OPs needed in order to consider adjacent particles ‘connected’, defaults to 0.32

• norm (float | None, optional) – an optional factor to normalize the result, defaults to 6, but passing None will reference the connectivity value for a perfectly crystalline hexagon (i.e. equations 8 and 10 in the SI from (Juarez, Lab on a Chip 2012))

• phase_rotate (ndarray | None, optional) – an optional \([N,N]\) array of complex phase factors (\((R_{jk})^n\)) to include a rotation between neighboring particles’ \(\psi_{n,j}\), defaults to None

• calc_3d (bool, optional) – whether to compute the 3D version of the crystal connectivity using spherical harmonic Steinhardt order parameters, defaults to False

Returns:

\([N,]\) array of real crystal connectivities

Return type:

ndarray

calc.bond_order.weighted_chi_graph(psis: ndarray, nei_bool: ndarray, phase_correct: ndarray, normalize_weights: bool = True) → Graph

WIP represent colloidal systems as a graph

Parameters:

• psis (ndarray) – \([N,]\) array of complex bond orientational order parameters.

• nei_bool (ndarray) – a \([N,N]\) boolean array indicating neighboring particles.

• phase_correct (ndarray) – array of complex phase factors (\((R_{jk})^n\)) to include a rotation between neighboring particles’ \(\psi_{n,j}\).

• normalize_weights (bool, optional) – whether to normalize edge weights by the norm of the inner product, defaults to True

Returns:

a graph-tool Graph with weighted edges according to bond orientational similarity.

Return type:

Graph

calc.bond_order.global_bond_order(psis: ndarray, nei_bool: ndarray, phase_correct: ndarray, g: Graph | None = None, root: int | None = None, ret_graph: bool = False, normalize_weights: bool = True) → tuple[ndarray, ndarray, Graph | None, int | None]

WIP compute globally phase-corrected bond orientational order parameters using a graph representation of the system

Parameters:

• psis (ndarray) – \([N,]\) array of complex bond orientational order parameters.

• nei_bool (ndarray) – a \([N,N]\) boolean array indicating neighboring particles.

• phase_correct (ndarray) – array of complex phase factors (\((R_{jk})^n\)) to include a rotation between neighboring particles’ \(\psi_{n,j}\).

• g (Graph | None, optional) – an optional pre-computed graph-tool Graph with weighted edges according to bond orientational similarity, defaults to None

• root (int | None, optional) – an optional root particle index from which to compute the Dijkstra search, defaults to None

• ret_graph (bool, optional) – whether to return the graph used in the calculation, defaults to False

• normalize_weights (bool, optional) – whether to normalize edge weights by the norm of the inner product, defaults to True

Returns:

\([N,]\) array of globally phase-corrected complex bond orientational order parameters, \([N,]\) array of graph distances from the root particle, optionally the graph used in the calculation, and the root particle index.

Return type:

ndarray[complex], ndarray (, Graph, int)

clusters submodule

Contains methods to calculate cluster morphological properties, including cluster-averaged bond order parameters, centers of-mass, gyration tensors/radii, and anisotropies.

calc.clusters.c6_clusters(c6: ndarray, nei: ndarray[bool]) → ndarray[int]

Uses graph theory to identify clusters of particles based on their defect status (as defined by c6) and neighbor connectivity (as defined by nei).

Param:

c6: an \([N,]\) array of complex bond order parameters for each particle

Param:

nei: a \([N,N]\) boolean array defining particle neighbors

Returns:

an \([N,]\) array of cluster indices for each particle

Return type:

ndarray[int]

calc.clusters.cluster_averager(cidx: ndarray[int]) → tuple[ndarray[bool], ndarray[int]]

Computes cluster membership map and cluster sizes from cluster indices. These arrays can be used to compute cluster-averaged properties. cluster_averager is used internally by other cluster analysis functions to compute cluster-averaged quantities efficiently, any cluster-averaged quantity can be computed array-qise ysing c_map @ quantity / c_sizes.

Param:

cidx: an \([N,]\) array of cluster indices for each particle

Returns:

a tuple containing: - c_map: a \([N,C]\) boolean array mapping particles to clusters - c_sizes: a \([C,]\) array of cluster sizes

Return type:

tuple[np.ndarray[bool], np.ndarray[int]]

calc.clusters.cluster_com(pts: ndarray, cidx: ndarray[int]) → ndarray

Computes the center of mass for each cluster based on particle positions and their cluster indices.

Param:

pts: an \([N,d]\) array of particle positions in \(d\) dimensions

Param:

cidx: an \([N,]\) array of cluster indices for each particle

Returns:

an \([C,d]\) array of cluster centers of mass

Return type:

ndarray

calc.clusters.cluster_rg(pts: ndarray, cidx: ndarray[int]) → ndarray

Computes the radius of gyration using the gyration_radius() method for each cluster based on particle positions and their cluster indices.

Param:

pts: an \([N,d]\) array of particle positions in \(d\) dimensions

Param:

cidx: an \([N,]\) array of cluster indices for each particle

Returns:

an \([C,]\) array of cluster radii of gyration

Return type:

ndarray

calc.clusters.cluster_gyrate(pts: ndarray, cidx: ndarray[int]) → ndarray

Computes the gyration tensor, using the gyration_tensor() method, for each cluster based on particle positions and their cluster indices.

Param:

pts: an \([N,d]\) array of particle positions in \(d\) dimensions

Param:

cidx: an \([N,]\) array of cluster indices for each particle

Returns:

an \([C,d,d]\) array of cluster gyration tensors

Return type:

ndarray

calc.clusters.cluster_shape(pts: ndarray, cidx: ndarray[int]) → ndarray

Computes the shape anisotropy (\(\kappa^2\)) using the shape_anisotropy() method for each cluster based on particle positions and their cluster indices.

Param:

pts: an \([N,d]\) array of particle positions in \(d\) dimensions

Param:

cidx: an \([N,]\) array of cluster indices for each particle

Returns:

an \([C,]\) array of cluster shape anisotropies \(\kappa^2\)

Return type:

ndarray