Diffraction Module¶
Overview
Computes a 2D diffraction pattern. 

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

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:
_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 highergrid_size
will lead to a higher resolution. The points are convolved with a Gaussian of width \(\sigma\), given bypeak_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) = N\), per 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 opensource GIXStapose application and its predecessor, diffractometer [JJ17].
Note
freud only supports diffraction patterns for cubic boxes.
 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 orNone
(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 orNone
(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.
 property k_values¶
kvalues.
 Type:
(
output_size
,)numpy.ndarray
 property k_vectors¶
kvectors.
 Type:
(
output_size
,output_size
, 3)numpy.ndarray
 plot(self, ax=None, cmap='afmhot', vmin=None, vmax=None)¶
Plot Diffraction Pattern.
 Parameters:
ax (
matplotlib.axes.Axes
, optional) – Axis to plot on. IfNone
, make a new figure and axis. (Default value =None
)cmap (str, optional) – Colormap name to use (Default value =
'afmhot'
).vmin (float) – Minimum of the color scale. Uses
4e6 * N_points
if not provided orNone
(Default value =None
).vmax (float) – Maximum of the color scale. Uses
0.7 * N_points
if not provided orNone
(Default value =None
).
 Returns:
Axis with the plot.
 Return type:
 to_image(self, cmap='afmhot', vmin=None, vmax=None)¶
Generates image of diffraction pattern.
 Parameters:
 Returns:
RGBA array of pixels.
 Return type:
((output_size, output_size, 4)
numpy.ndarray
)
 class freud.diffraction.StaticStructureFactorDebye¶
Bases:
_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\).
Note
For 2D systems freud uses the Bessel function \(J_0\) instead of the \(\text{sinc}\) function in the equation above. See [Wie12] for more information. For users wishing to calculate the structure factor of quasi 2D systems (i.e. a 2D simulation is used to model a real system such as particles on a 2D interface or similar) the 3D formula should be used. In these cases users should use a 3D box with its longest dimension being in the zdirection and particle positions of the form \((x, y, 0)\).
This implementation uses an evenly spaced number of \(k\) points between k_min` and
k_max
. Ifk_min
is set to 0 (the default behavior), the computed structure factor will show \(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 systemN_total
to thecompute()
method. The normalization criterion is based on the FaberZiman 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:
num_k_values (unsigned int) – Number of values to use 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
 compute(self, system, query_points=None, N_total=None, reset=True)¶
Computes static structure factor.
Example for a single component system:
>>> box, points = freud.data.make_random_system(10, 100, seed=0) >>> sf = freud.diffraction.StaticStructureFactorDebye( ... num_k_values=100, k_max=10, k_min=0 ... ) >>> sf.compute((box, points)) freud.diffraction.StaticStructureFactorDebye(...)
Example for partial mixed structure factor for a multiple component system with types A and B:
>>> N_particles = 100 >>> box, points = freud.data.make_random_system( ... 10, N_particles, seed=0 ... ) >>> A_points = points[:N_particles//2] >>> B_points = points[N_particles//2:] >>> sf = freud.diffraction.StaticStructureFactorDebye( ... num_k_values=100, k_max=10, k_min=0 ... ) >>> sf.compute( ... system=(box, A_points), ... query_points=B_points, ... N_total=N_particles ... ) freud.diffraction.StaticStructureFactorDebye(...)
 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 ifNone
. See class documentation for information about the normalization of partial structure factors. IfNone
, 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_values¶
The \(k\) values for the calculation.
 Type:
 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:
 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. IfNone
, make a new figure and axis. (Default value =None
) Returns:
Axis with the plot.
 Return type:
 class freud.diffraction.StaticStructureFactorDirect¶
Bases:
_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
Currently 2D boxes are not supported for this method. Use Debye instead.
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 systemN_total
to thecompute()
method. The normalization criterion is based on the FaberZiman 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:
 bin_edges¶
The edges of each bin of \(k\).
 Type:
 compute(self, system, query_points=None, N_total=None, reset=True)¶
Computes static structure factor.
Example for a single component system:
>>> box, points = freud.data.make_random_system(10, 100, seed=0) >>> sf = freud.diffraction.StaticStructureFactorDirect( ... bins=100, k_max=10, k_min=0 ... ) >>> sf.compute((box, points)) freud.diffraction.StaticStructureFactorDirect(...)
Example for partial mixed structure factor for multiple component system with types A and B:
>>> N_particles = 100 >>> box, points = freud.data.make_random_system( ... 10, N_particles, seed=0 ... ) >>> A_points = points[:N_particles//2] >>> B_points = points[N_particles//2:] >>> 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 ... ) freud.diffraction.StaticStructureFactorDirect(...)
 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 nonorthorhombic 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 ifNone
. See class documentation for information about the normalization of partial structure factors. IfNone
, 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).
 property k_points¶
The \(\vec{k}\) points used in the calculation.
 Type:
 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:
 num_sampled_k_points¶
The target number of \(\vec{k}\) points to use when constructing \(k\) space grid.
 Type:
 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. IfNone
, make a new figure and axis. (Default value =None
) Returns:
Axis with the plot.
 Return type: