Diffraction Module

Overview

freud.diffraction.DiffractionPattern

Computes a 2D diffraction pattern.

freud.diffraction.StaticStructureFactorDebye

Computes a 1D static structure factor using the Debye scattering equation.

freud.diffraction.StaticStructureFactorDirect

Computes a 1D static structure factor by operating on a \(k\) space grid.

Details

The freud.diffraction module provides functions for computing the diffraction pattern of particles in systems with long range order.

Stability

freud.diffraction is unstable. When upgrading from version 2.x to 2.y (y > x), existing freud scripts may need to be updated. The API will be finalized in a future release.

class freud.diffraction.DiffractionPattern(grid_size=512, output_size=None)

Bases: freud.util._Compute

Computes a 2D diffraction pattern.

The diffraction image represents the scattering of incident radiation, and is useful for identifying translational and/or rotational symmetry present in the system. This class computes the static structure factor \(S(\vec{k})\) for a plane of wavevectors \(\vec{k}\) orthogonal to a view axis. The view orientation \((1, 0, 0, 0)\) defaults to looking down the \(z\) axis (at the \(xy\) plane). The points in the system are converted to fractional coordinates, then binned into a grid whose resolution is given by grid_size. A higher grid_size will lead to a higher resolution. The points are convolved with a Gaussian of width \(\sigma\), given by peak_width. This convolution is performed as a multiplication in Fourier space. The computed diffraction pattern can be accessed as a square array of shape (output_size, output_size).

The \(\vec{k}=0\) peak is always located at index (output_size // 2, output_size // 2) and is normalized to have a value of \(S(\vec{k}=0) = 1\) (not \(N\), a common convention). The remaining \(\vec{k}\) vectors are computed such that each peak in the diffraction pattern satisfies the relationship \(\vec{k} \cdot \vec{R} = 2 \pi N\) for some integer \(N\) and lattice vector of the system \(\vec{R}\). See the reciprocal lattice Wikipedia page for more information.

This method is based on the implementations in the open-source GIXStapose application and its predecessor, diffractometer [JJ17].

Parameters
  • grid_size (unsigned int) – Resolution of the diffraction grid (Default value = 512).

  • output_size (unsigned int) – Resolution of the output diffraction image, uses grid_size if not provided or None (Default value = None).

compute(self, system, view_orientation=None, zoom=4, peak_width=1, reset=True)

Computes diffraction pattern.

Parameters
  • system – Any object that is a valid argument to freud.locality.NeighborQuery.from_system.

  • view_orientation ((\(4\)) numpy.ndarray, optional) – View orientation. Uses \((1, 0, 0, 0)\) if not provided or None (Default value = None).

  • zoom (float, optional) – Scaling factor for incident wavevectors (Default value = 4).

  • peak_width (float, optional) – Width of Gaussian convolved with points, in system length units (Default value = 1).

  • reset (bool, optional) – Whether to erase the previously computed values before adding the new computations; if False, will accumulate data (Default value = True).

property diffraction

(output_size, output_size) numpy.ndarray: Diffraction pattern.

grid_size

Resolution of the diffraction grid.

Type

int

property k_values

k-values.

Type

(output_size, ) numpy.ndarray

property k_vectors

(output_size, output_size, 3) numpy.ndarray: k-vectors.

output_size

Resolution of the output diffraction image.

Type

int

plot(self, ax=None, cmap='afmhot', vmin=4e-6, vmax=0.7)

Plot Diffraction Pattern.

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

  • cmap (str, optional) – Colormap name to use (Default value = 'afmhot').

  • vmin (float, optional) – Minimum of the color scale (Default value = 4e-6).

  • vmax (float, optional) – Maximum of the color scale (Default value = 0.7).

Returns

Axis with the plot.

Return type

(matplotlib.axes.Axes)

to_image(self, cmap='afmhot', vmin=4e-6, vmax=0.7)

Generates image of diffraction pattern.

Parameters
  • cmap (str, optional) – Colormap name to use (Default value = 'afmhot').

  • vmin (float, optional) – Minimum of the color scale (Default value = 4e-6).

  • vmax (float, optional) – Maximum of the color scale (Default value = 0.7).

Returns

RGBA array of pixels.

Return type

((output_size, output_size, 4) numpy.ndarray)

class freud.diffraction.StaticStructureFactorDebye

Bases: freud.diffraction._StaticStructureFactor

Computes a 1D static structure factor using the Debye scattering equation.

This computes the static structure factor \(S(k)\) at given \(k\) values by averaging over all \(\vec{k}\) vectors of the same magnitude. Note that freud employs the physics convention in which \(k\) is used, as opposed to the crystallographic one where \(q\) is used. The relation is \(k=2 \pi q\). The static structure factor calculation is implemented using the Debye scattering equation:

\[S(k) = \frac{1}{N} \sum_{i=0}^{N} \sum_{j=0}^{N} \text{sinc}(k r_{ij})\]

where \(N\) is the number of particles, \(\text{sinc}\) function is defined as \(\sin x / x\) (no factor of \(\pi\) as in some conventions). For more information see this Wikipedia article. For a full derivation see [FB09]. Note that the definition requires \(S(0) = N\).

The Debye implementation provides a much faster algorithm, but gives worse results than freud.diffraction.StaticStructureFactorDirect at low k values.

Note

This code assumes all particles have a form factor \(f\) of 1.

Partial structure factors can be computed by providing a set of query_points and the total number of points in the system N_total to the compute() method. The normalization criterion is based on the Faber-Ziman formalism. For particle types \(\alpha\) and \(\beta\), we compute the total scattering function as a sum of the partial scattering functions as:

\[S(k) - 1 = \sum_{\alpha}\sum_{\beta} \frac{N_{\alpha} N_{\beta}}{N_{total}^2} \left(S_{\alpha \beta}(k) - 1\right)\]
Parameters
  • bins (unsigned int) – Number of bins in \(k\) space.

  • k_max (float) – Maximum \(k\) value to include in the calculation.

  • k_min (float, optional) – Minimum \(k\) value included in the calculation. Note that there are practical restrictions on the validity of the calculation in the long wavelength regime, see min_valid_k (Default value = 0).

property S_k

Static structure factor \(S(k)\) values.

Type

(\(N_{bins}\),) numpy.ndarray

bin_centers

The centers of each bin of \(k\).

Type

numpy.ndarray

bin_edges

The edges of each bin of \(k\).

Type

numpy.ndarray

bounds

A tuple indicating upper and lower bounds of the histogram.

Type

tuple

compute(self, system, query_points=None, N_total=None, reset=True)

Computes static structure factor.

Example for a single component system:

>>> sf = freud.diffraction.StaticStructureFactorDebye(
...     bins=100, k_max=10, k_min=0
... )
>>> sf.compute((box, points))

Example for partial mixed structure factor for a multiple component system with types A and B:

>>> sf = freud.diffraction.StaticStructureFactorDebye(
...     bins=100, k_max=10, k_min=0
... )
>>> sf.compute(
...     system=(box, A_points),
...     query_points=B_points,
...     N_total=N_particles
... )
Parameters
  • system – Any object that is a valid argument to freud.locality.NeighborQuery.from_system. Note that box is allowed to change when accumulating average static structure factor.

  • query_points ((\(N_{query\_points}\), 3) numpy.ndarray, optional) – Query points used to calculate the partial structure factor. Uses the system’s points if None. See class documentation for information about the normalization of partial structure factors. If None, the full scattering is computed. (Default value = None).

  • N_total (int, optional) – Total number of points in the system. This is required if query_points are provided. See class documentation for information about the normalization of partial structure factors.

  • reset (bool, optional) – Whether to erase the previously computed values before adding the new computation; if False, will accumulate data (Default value: True).

k_max

Maximum value of k at which to calculate the structure factor.

Type

float

k_min

Minimum value of k at which to calculate the structure factor.

Type

float

property min_valid_k

Minimum valid value of k for the computed system box, equal to \(2\pi/(L/2)=4\pi/L\) where \(L\) is the minimum side length. For more information see [LP16].

Type

float

nbins

The number of bins in the histogram.

Type

int

plot(self, ax=None, **kwargs)

Plot static structure factor.

Note

This function plots \(S(k)\) for values above min_valid_k.

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)

class freud.diffraction.StaticStructureFactorDirect

Bases: freud.diffraction._StaticStructureFactor

Computes a 1D static structure factor by operating on a \(k\) space grid.

This computes the static structure factor \(S(k)\) at given \(k\) values by averaging over all \(\vec{k}\) vectors directions of the same magnitude. Note that freud employs the physics convention in which \(k\) is used, as opposed to the crystallographic one where \(q\) is used. The relation is \(k=2 \pi q\). This is implemented using the following formula:

\[S(\vec{k}) = \frac{1}{N} \sum_{i=0}^{N} \sum_{j=0}^N e^{i\vec{k} \cdot \vec{r}_{ij}}\]

where \(N\) is the number of particles. Note that the definition requires \(S(0) = N\).

This implementation provides a much slower algorithm, but gives better results than the freud.diffraction.StaticStructureFactorDebye method at low k values.

The \(\vec{k}\) vectors are sampled isotropically from a grid defined by the box’s reciprocal lattice vectors. This sampling of reciprocal space is based on the MIT licensed Dynasor library, modified to use parallelized C++ and to support larger ranges of \(k\) values. For more information see [FSEW21].

Note

This code assumes all particles have a form factor \(f\) of 1.

Partial structure factors can be computed by providing query_points and total number of points in the system N_total to the compute() method. The normalization criterion is based on the Faber-Ziman formalism. For particle types \(\alpha\) and \(\beta\), we compute the total scattering function as a sum of the partial scattering functions as:

\[S(k) - 1 = \sum_{\alpha}\sum_{\beta} \frac{N_{\alpha} N_{\beta}}{N_{total}^2} \left(S_{\alpha \beta}(k) - 1\right)\]
Parameters
  • bins (unsigned int) – Number of bins in \(k\) space.

  • k_max (float) – Maximum \(k\) value to include in the calculation.

  • k_min (float, optional) – Minimum \(k\) value included in the calculation. Note that there are practical restrictions on the validity of the calculation in the long wavelength regime, see min_valid_k (Default value = 0).

  • num_sampled_k_points (unsigned int, optional) – The desired number of \(\vec{k}\) vectors to sample from the reciprocal lattice grid. If set to 0, all \(\vec{k}\) vectors are used. If greater than 0, the \(\vec{k}\) vectors are sampled from the full grid with uniform radial density, resulting in a sample of num_sampled_k_points vectors on average (Default value = 0).

property S_k

Static structure factor \(S(k)\) values.

Type

(\(N_{bins}\),) numpy.ndarray

bin_centers

The centers of each bin of \(k\).

Type

numpy.ndarray

bin_edges

The edges of each bin of \(k\).

Type

numpy.ndarray

bounds

A tuple indicating upper and lower bounds of the histogram.

Type

tuple

compute(self, system, query_points=None, N_total=None, reset=True)

Computes static structure factor.

Example for a single component system:

>>> sf = freud.diffraction.StaticStructureFactorDirect(
...     bins=100, k_max=10, k_min=0
... )
>>> sf.compute((box, points))

Example for partial mixed structure factor for multiple component system AB:

>>> sf = freud.diffraction.StaticStructureFactorDirect(
...     bins=100, k_max=10, k_min=0
... )
>>> sf.compute(
...     (box, A_points),
...     query_points=B_points, N_total=N_particles
... )
Parameters
  • system – Any object that is a valid argument to freud.locality.NeighborQuery.from_system. Note that box is allowed to change when accumulating average static structure factor. For non-orthorhombic boxes the points are wrapped into a orthorhombic box.

  • query_points ((\(N_{query\_points}\), 3) numpy.ndarray, optional) – Query points used to calculate the partial structure factor. Uses the system’s points if None. See class documentation for information about the normalization of partial structure factors. If None, the full scattering is computed. (Default value = None).

  • N_total (int, optional) – Total number of points in the system. This is required if query_points are provided. See class documentation for information about the normalization of partial structure factors.

  • reset (bool, optional) – Whether to erase the previously computed values before adding the new computation; if False, will accumulate data (Default value = True).

k_max

Maximum value of k at which to calculate the structure factor.

Type

float

k_min

Minimum value of k at which to calculate the structure factor.

Type

float

property k_points

The \(\vec{k}\) points used in the calculation.

Type

numpy.ndarray

property min_valid_k

Minimum valid value of k for the computed system box, equal to \(2\pi/L\) where \(L\) is the minimum side length. For more information see [LP16].

Type

float

nbins

The number of bins in the histogram.

Type

int

num_sampled_k_points

The target number of \(\vec{k}\) points to use when constructing \(k\) space grid.

Type

int

plot(self, ax=None, **kwargs)

Plot static structure factor.

Note

This function plots \(S(k)\) for values above min_valid_k.

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)