Order Module#

Overview

freud.order.ContinuousCoordination

Computes the continuous local coordination number.

freud.order.Cubatic

Compute the cubatic order parameter [HAG15] for a system of particles using simulated annealing instead of Newton-Raphson root finding.

freud.order.Nematic

Compute the nematic order parameter for a system of particles.

freud.order.Hexatic

Calculates the \(k\)-atic order parameter for 2D systems.

freud.order.Steinhardt

Compute one or more of the rotationally invariant Steinhardt order parameter \(q_l\) or \(w_l\) for a set of points [Ste83].

freud.order.SolidLiquid

Identifies solid-like clusters using dot products of \(q_{lm}\).

freud.order.RotationalAutocorrelation

Calculates a measure of total rotational autocorrelation.

Details

The freud.order module contains functions which compute order parameters for the whole system or individual particles. Order parameters take bond order data and interpret it in some way to quantify the degree of order in a system using a scalar value. This is often done through computing spherical harmonics of the bond order diagram, which are the spherical analogue of Fourier Transforms.

class freud.order.ContinuousCoordination#

Bases: _PairCompute

Computes the continuous local coordination number.

The ContinuousCoordination class implements extensions of the Voronoi discrete coordination number to the real numbers. The formulas for the various implementations are:

Power:

\[CN_p = N^{2.0 - m} \sum_{i=1}^{k} {\left[\left(\frac{V_i}{V}\right)^{m}\right]}^{-1}\]

Log:

\[CN_{log} = \frac{-1}{\log{N}} \sum_{i=1}^{k}\log{\left(\frac{V_i}{V}\right)}\]

Exponential:

\[CN_{exp} = \sum_{i=1}^{k}\exp{\left(\frac{V_i}{V} - \frac{1}{N} \right)}\]

where \(k\) is the number of neighbors a particle has, \(V_i\) is the volume of the pyramid (or area of the triangle in 2D) whose base is the Voronoi polytope facet between the central particle and neighbor \(i\) and whose height is half the distance vector, and \(V\) is the volume/area of the Voronoi polytope.

Note

When using multiple powers, space them out to avoid high correlations. A minimum spacing of 2.0 is recommended with even larger values leading to less correlation.

Parameters:

  • powers (list[float], optional) – The powers to compute the continuous coordination number for. The default value indicates only compute for power 2.0. (Default value: None)

  • compute_log (bool, optional) – Whether to compute the log continuous coordination number. (Default value: True)

  • compute_exp (bool, optional) – Whether to compute the exp continuous coordination number. (Default value: True)

compute(self, system=None, voronoi=None)#

Calculates the local coordination number for the specified points.

Example:

>>> import freud
>>> box, points = freud.data.make_random_system(10, 100, seed=0)
>>> # Compute ContinuousCoordination
>>> coord = freud.order.ContinuousCoordination([2, 4], True)
>>> coord.compute(system=(box, points))
freud.order.ContinuousCoordination(...)
Parameters:

  • system (optional) – Any object that is a valid argument to freud.locality.NeighborQuery.from_system. (Default value: None).

  • voronoi (freud.locality.Voronoi, optional) – A precomputed Voronoi compute object. If provided, the object is assumed to have been computed already, and system is ignored. (Default value: None).

compute_exp#

Whether to compute the exponential coordination number.

Type:

bool

compute_log#

Whether to compute the log continuous coordination number.

Type:

bool

property coordination#

coordination of points per query point.

Coordination numbers are in order of selected powers, log, and exp.

Type:

(\((N_{points}, N_{coord}\)) numpy.ndarray)

default_query_args#

No default query arguments.

number_of_coordinations#

The number of coordination numbers computed.

Type:

int

property particle_order#

coordination of points per query point.

Coordination numbers are in order of selected powers, log, and exp.

Type:

(\((N_{points}, N_{coord}\)) numpy.ndarray

powers#

The powers to compute the continuous coordination number.

Changes to this property are not reflected when computing coordination numbers.

Type:

list[float]

class freud.order.Cubatic#

Bases: _Compute

Compute the cubatic order parameter [HAG15] for a system of particles using simulated annealing instead of Newton-Raphson root finding.

Parameters:

  • t_initial (float) – Starting temperature.

  • t_final (float) – Final temperature.

  • scale (float) – Scaling factor to reduce temperature.

  • n_replicates (unsigned int, optional) – Number of replicate simulated annealing runs. (Default value = 1).

  • seed (unsigned int, optional) – Random seed to use in calculations. If None, system time is used. (Default value = None).

compute(self, orientations)#

Calculates the per-particle and global order parameter.

Parameters:

orientations ((\(N_{particles}\), 4) numpy.ndarray) – Orientations as angles to use in computation.

property cubatic_tensor#

Rank 4 homogeneous tensor representing the optimal system-wide coordinates.

Type:

\(\left(3, 3, 3, 3 \right)\) numpy.ndarray

property global_tensor#

Rank 4 tensor corresponding to the global orientation. Computed from all orientations.

Type:

\(\left(3, 3, 3, 3 \right)\) numpy.ndarray

n_replicates#

Number of replicate simulated annealing runs.

Type:

unsigned int

property order#

Cubatic order parameter of the system.

Type:

float

property orientation#

The quaternion of global orientation.

Type:

\(\left(4 \right)\) numpy.ndarray

property particle_order#

Order parameter.

Type:

\(\left(N_{particles} \right)\) numpy.ndarray

scale#

The scale.

Type:

float

seed#

Random seed to use in calculations.

Type:

unsigned int

t_final#

The value of the final temperature.

Type:

float

t_initial#

The value of the initial temperature.

Type:

float

class freud.order.Hexatic#

Bases: _PairCompute

Calculates the \(k\)-atic order parameter for 2D systems.

The \(k\)-atic order parameter (called the hexatic order parameter for \(k = 6\)) is analogous to Steinhardt order parameters, and is used to measure order in the bonds of 2D systems.

The \(k\)-atic order parameter for a particle \(i\) and its \(N_b\) neighbors \(j\) is given by:

\(\psi_k \left( i \right) = \frac{1}{N_b} \sum \limits_{j=1}^{N_b} e^{i k \phi_{ij}}\)

The parameter \(k\) governs the symmetry of the order parameter and typically matches the number of neighbors to be found for each particle. The quantity \(\phi_{ij}\) is the angle between the vector \(r_{ij}\) and \(\left(1, 0\right)\).

If the weighted mode is enabled, contributions of each neighbor are weighted. Neighbor weights \(w_{ij}\) default to 1 but are defined for a freud.locality.NeighborList from freud.locality.Voronoi or one with user-provided weights. The formula is modified as follows:

\(\psi'_k \left( i \right) = \frac{1}{\sum_{j=1}^{N_b} w_{ij}} \sum \limits_{j=1}^{N_b} w_{ij} e^{i k \phi_{ij}}\)

The hexatic order parameter as written above is complex-valued. The magnitude of the complex value, np.abs(hex_order.particle_order), is frequently what is desired when determining the \(k\)-atic order for each particle. The complex phase angle np.angle(hex_order.particle_order) indicates the orientation of the bonds as an angle measured counterclockwise from the vector \(\left(1, 0\right)\). The complex valued order parameter is not rotationally invariant because of this phase angle, but the magnitude is rotationally invariant.

Note

2D: freud.order.Hexatic is only defined for 2D systems. The points must be passed in as [x, y, 0].

Parameters:

  • k (unsigned int, optional) – Symmetry of order parameter (Default value = 6).

  • weighted (bool, optional) – Determines whether to use neighbor weights in the computation of spherical harmonics over neighbors. If enabled and used with a Voronoi neighbor list, this results in the 2D Minkowski Structure Metrics \(\psi'_k\) [MKSTM13] (Default value = False).

compute(self, system, neighbors=None)#

Calculates the hexatic order parameter.

Example:

>>> box, points = freud.data.make_random_system(
...     box_size=10, num_points=100, is2D=True, seed=0)
>>> # Compute the hexatic (6-fold) order for the 2D system
>>> hex_order = freud.order.Hexatic(k=6)
>>> hex_order.compute(system=(box, points))
freud.order.Hexatic(...)
>>> print(hex_order.particle_order)
[...]
Parameters:
default_query_args#

The default query arguments are {'mode': 'nearest', 'num_neighbors': self.k}.

k#

Symmetry of the order parameter.

Type:

unsigned int

property particle_order#

Order parameter.

Type:

\(\left(N_{particles} \right)\) numpy.ndarray

plot(self, ax=None)#

Plot order parameter distribution.

Parameters:

ax (matplotlib.axes.Axes, optional) – Axis to plot on. If None, make a new figure and axis (Default value = None).

Returns:

Axis with the plot.

Return type:

(matplotlib.axes.Axes)

weighted#

Whether neighbor weights were used in the computation.

Type:

bool

class freud.order.Nematic#

Bases: _Compute

Compute the nematic order parameter for a system of particles.

Note

In some cases, such as HOOMD-blue simulations, orientations are represented as quaternions which rotate a particle from an initial reference orientation. To convert from the quaternion representation to the orientation vectors taken by freud, one should rotate the reference orientation by the particle’s orientation quaternion. An example of this process using the rowan library is shown below.

Example:

>>> import rowan
>>> import numpy as np
>>> import freud
>>> quats = np.array([[1, 0, 0, 0]] * 100)
>>> starting_orientation = np.array([1, 1, 0])
>>> orientations = rowan.rotate(quats,starting_orientation)
>>> nematic = freud.order.Nematic()
>>> nematic.compute(orientations)
freud.order.Nematic()
compute(self, orientations)#

Calculates the per-particle and global order parameter.

Example:

>>> orientations = np.array([[1, 0, 0]] * 100)
>>> nematic = freud.order.Nematic()
>>> nematic.compute(orientations)
freud.order.Nematic()
>>> print(nematic.order)
1.0
Parameters:

orientations (\(\left(N_{particles}, 3 \right)\) numpy.ndarray) – Orientation vectors for which to calculate the order parameter.

property director#

The average nematic director.

Type:

\(\left(3 \right)\) numpy.ndarray

property nematic_tensor#

3x3 matrix corresponding to the average particle orientation.

Type:

\(\left(3, 3 \right)\) numpy.ndarray

property order#

Nematic order parameter of the system.

Type:

float

property particle_tensor#

\(\left(N_{particles}, 3, 3 \right)\) numpy.ndarray: One 3x3 matrix per-particle corresponding to each individual particle orientation.

class freud.order.RotationalAutocorrelation#

Bases: _Compute

Calculates a measure of total rotational autocorrelation.

For any calculation of rotational correlations of extended (i.e. non-point) particles, encoding the symmetries of these particles is crucial to appropriately determining correlations. For systems of anisotropic particles in three dimensions, representing such equivalence can be quite mathematically challenging. This calculation is based on the hyperspherical harmonics as laid out in [KDvAG19]. Generalizations of spherical harmonics to four dimensions, hyperspherical harmonics provide a natural basis for periodic functions in 4 dimensions just as harmonic functions (sines and cosines) or spherical harmonics do in lower dimensions. The idea behind this calculation is to embed orientation quaternions into a 4-dimensional space and then use hyperspherical harmonics to find correlations in a symmetry-aware fashion.

The choice of the hyperspherical harmonic parameter \(l\) determines the symmetry of the functions. The output is not a correlation function, but rather a scalar value that measures total system orientational correlation with an initial state. As such, the output can be treated as an order parameter measuring degrees of rotational (de)correlation. For analysis of a trajectory, the compute call needs to be done at each trajectory frame.

Parameters:

l (int) – Order of the hyperspherical harmonic. Must be a positive, even integer.

compute(self, ref_orientations, orientations)#

Calculates the rotational autocorrelation function for a single frame.

Parameters:

  • ref_orientations ((\(N_{orientations}\), 4) numpy.ndarray) – Orientations for the initial frame.

  • orientations ((\(N_{orientations}\), 4) numpy.ndarray) – Orientations for the frame of interest.

l#

The azimuthal quantum number, which defines the order of the hyperspherical harmonic.

Type:

int

property order#

Autocorrelation of the system.

Type:

float

property particle_order#

Rotational autocorrelation values calculated for each orientation.

Type:

(\(N_{orientations}\)) numpy.ndarray

class freud.order.SolidLiquid#

Bases: _PairCompute

Identifies solid-like clusters using dot products of \(q_{lm}\).

The solid-liquid order parameter [FHND10, tW95] uses a Steinhardt-like approach to identify solid-like particles. First, a bond parameter \(q_l(i, j)\) is computed for each neighbor bond.

If normalize_q is true (default), the bond parameter is given by \(q_l(i, j) = \frac{\sum \limits_{m=-l}^{l} \text{Re}~q_{lm}(i) q_{lm}^*(j)} {\sqrt{\sum \limits_{m=-l}^{l} \lvert q_{lm}(i) \rvert^2} \sqrt{\sum \limits_{m=-l}^{l} \lvert q_{lm}(j) \rvert^2}}\)

If normalize_q is false, then the denominator of the above expression is left out.

Next, the bonds are filtered to keep only “solid-like” bonds with \(q_l(i, j)\) above a cutoff value \(q_{threshold}\).

If a particle has more than \(S_{threshold}\) solid-like bonds, then the particle is considered solid-like. Finally, solid-like particles are clustered.

Parameters:

  • l (unsigned int) – Spherical harmonic quantum number l.

  • q_threshold (float) – Value of dot product threshold when evaluating \(q_l(i, j)\) to determine if a bond is solid-like. For \(l=6\), 0.7 is generally good for FCC or BCC structures [FHND10].

  • solid_threshold (unsigned int) – Minimum required number of adjacent solid-like bonds for a particle to be considered solid-like for clustering. For \(l=6\), 6-8 is generally good for FCC or BCC structures.

  • normalize_q (bool) – Whether to normalize the dot product (Default value = True).

property cluster_idx#

\(\left(N_{particles}\right)\) numpy.ndarray: Solid-like cluster indices for each particle.

property cluster_sizes#

The sizes of all clusters.

Type:

\((N_{clusters}, )\) np.ndarray

compute(self, system, neighbors=None)#

Compute the order parameter.

Parameters:
default_query_args#

No default query arguments.

l#

Spherical harmonic quantum number l.

Type:

unsigned int

property largest_cluster_size#

The largest cluster size.

Type:

unsigned int

property nlist#

Neighbor list of solid-like bonds.

Type:

freud.locality.NeighborList

normalize_q#

Whether the dot product is normalized.

Type:

bool

property num_connections#

The number of solid-like bonds for each particle.

Type:

\(\left(N_{particles}\right)\) numpy.ndarray

property particle_harmonics#

\(\left(N_{particles}, 2*l+1\right)\) numpy.ndarray: The raw array of overline{q}_{lm}(i). The array is provided in the order given by fsph: \(m = 0, 1, ..., l, -1, ..., -l\).

plot(self, ax=None)#

Plot solid-like cluster distribution.

Parameters:

ax (matplotlib.axes.Axes, optional) – Axis to plot on. If None, make a new figure and axis (Default value = None).

Returns:

Axis with the plot.

Return type:

(matplotlib.axes.Axes)

q_threshold#

Value of dot product threshold.

Type:

float

property ql_ij#

Bond dot products \(q_l(i, j)\). Indexed by the elements of self.nlist.

Type:

\(\left(N_{bonds}\right)\) numpy.ndarray

solid_threshold#

Value of number-of-bonds threshold.

Type:

float

class freud.order.Steinhardt#

Bases: _PairCompute

Compute one or more of the rotationally invariant Steinhardt order parameter \(q_l\) or \(w_l\) for a set of points [Ste83].

Implements the local rotationally invariant \(q_l\) or \(w_l\) order parameter described by Steinhardt.

First, we describe the computation of \(q_l(i)\). For a particle \(i\), we calculate the quantity \(q_{lm}\) by summing the spherical harmonics between particle \(i\) and its neighbors \(j\) in a local region:

\[q_{lm}(i) = \frac{1}{N_b} \sum \limits_{j=1}^{N_b} Y_{lm}(\theta(\vec{r}_{ij}), \phi(\vec{r}_{ij}))\]

Then the \(q_l\) order parameter is computed by combining the \(q_{lm}\) in a rotationally invariant fashion to remove local orientational order:

\[q_l(i) = \sqrt{\frac{4\pi}{2l+1} \sum \limits_{m=-l}^{l} |q_{lm}(i)|^2 }\]

If the wl parameter is True, this class computes the quantity \(w_l\), defined as a weighted average over the \(q_{lm}(i)\) values using Wigner 3-j symbols (related to Clebsch-Gordan coefficients). The resulting combination is rotationally invariant:

\[\begin{split}w_l(i) = \sum \limits_{m_1 + m_2 + m_3 = 0} \begin{pmatrix} l & l & l \\ m_1 & m_2 & m_3 \end{pmatrix} q_{lm_1}(i) q_{lm_2}(i) q_{lm_3}(i)\end{split}\]

If wl is True, then setting the wl_normalize parameter to True will normalize the \(w_l\) order parameter as follows (if wl=False, wl_normalize has no effect):

\[\begin{split}w_l(i) = \frac{ \sum \limits_{m_1 + m_2 + m_3 = 0} \begin{pmatrix} l & l & l \\ m_1 & m_2 & m_3 \end{pmatrix} q_{lm_1}(i) q_{lm_2}(i) q_{lm_3}(i)} {\left(\sum \limits_{m=-l}^{l} |q_{lm}(i)|^2 \right)^{3/2}}\end{split}\]

If average is True, the class computes a variant of this order parameter that performs an average over the first and second shell combined [LD08]. To compute this parameter, we perform a second averaging over the first neighbor shell of the particle to implicitly include information about the second neighbor shell. This averaging is performed by replacing the value \(q_{lm}(i)\) in the original definition by \(\overline{q}_{lm}(i)\), the average value of \(q_{lm}(k)\) over all the \(N_b\) neighbors \(k\) of particle \(i\), including particle \(i\) itself:

\[\overline{q}_{lm}(i) = \frac{1}{N_b} \sum \limits_{k=0}^{N_b} q_{lm}(k)\]

If weighted is True, the contributions of each neighbor are weighted. Neighbor weights \(w_{ij}\) are defined for a freud.locality.NeighborList obtained from freud.locality.Voronoi or one with user-provided weights, and default to 1 if not otherwise provided. The formulas are modified as follows, replacing \(q_{lm}(i)\) with the weighted value \(q'_{lm}(i)\):

\[q'_{lm}(i) = \frac{1}{\sum_{j=1}^{N_b} w_{ij}} \sum \limits_{j=1}^{N_b} w_{ij} Y_{lm}(\theta(\vec{r}_{ij}), \phi(\vec{r}_{ij}))\]

Note

The value of per-particle order parameter will be set to NaN for particles with no neighbors. We choose this value rather than setting the order parameter to 0 because in more complex order parameter calculations (such as when computing the \(w_l\)), it is possible to observe a value of 0 for the per-particle order parameter even with a finite number of neighbors. If you would like to ignore this distinction, you can mask the output order parameter values using NumPy: numpy.nan_to_num(particle_order).

Parameters:

  • l (unsigned int or sequence of unsigned int) – One or more spherical harmonic quantum number l’s used to compute the Steinhardt order parameter.

  • average (bool, optional) – Determines whether to calculate the averaged Steinhardt order parameter (Default value = False).

  • wl (bool, optional) – Determines whether to use the \(w_l\) version of the Steinhardt order parameter (Default value = False).

  • weighted (bool, optional) – Determines whether to use neighbor weights in the computation of spherical harmonics over neighbors. If enabled and used with a Voronoi neighbor list, this results in the 3D Minkowski Structure Metrics \(q'_l\) [MKSTM13] (Default value = False).

  • wl_normalize (bool, optional) – Determines whether to normalize the \(w_l\) version of the Steinhardt order parameter (Default value = False).

average#

Whether the averaged Steinhardt order parameter was calculated.

Type:

bool

compute(self, system, neighbors=None)#

Compute the order parameter.

Example:

>>> box, points = freud.data.make_random_system(10, 100, seed=0)
>>> ql = freud.order.Steinhardt(l=6)
>>> ql.compute((box, points), {'r_max':3})
freud.order.Steinhardt(l=6, ...)
Parameters:
default_query_args#

No default query arguments.

l#

Spherical harmonic quantum number l.

Type:

unsigned int

property order#

The system wide normalization of the order parameter, computed by averaging the \(q_{lm}\) values (or \(\overline{q}_{lm}\) values if average is enabled) over all particles before computing the rotationally-invariant order parameter.

Type:

float

property particle_harmonics#

\(\left(N_{particles}, 2l+1\right)\) numpy.ndarray: The raw array of \(q_{lm}(i)\). The array is provided in the order given by fsph: \(m = 0, 1, ..., l, -1, ..., -l\).

property particle_order#

\(\left(N_{particles}, N_l \right)\) numpy.ndarray: Variant of the Steinhardt order parameter for each particle (filled with nan for particles with no neighbors).

plot(self, ax=None)#

Plot order parameter distribution.

Parameters:

ax (matplotlib.axes.Axes, optional) – Axis to plot on. If None, make a new figure and axis (Default value = None).

Returns:

Axis with the plot.

Return type:

(matplotlib.axes.Axes)

property ql#

\(\left(N_{particles}, N_l\right)\) numpy.ndarray: \(q_l\) Steinhardt order parameter for each particle (filled with nan for particles with no neighbors). This is always available, no matter which other options are selected. It obeys the weighted argument but otherwise returns the “plain” \(q_l\) regardless of average, wl, wl_normalize.

weighted#

Whether neighbor weights were used in the computation.

Type:

bool

wl#

Whether the \(w_l\) version of the Steinhardt order parameter was used.

Type:

bool