# freud¶

## Overview¶

The **freud** Python library provides a simple, flexible, powerful set of tools
for analyzing trajectories obtained from molecular dynamics or Monte Carlo
simulations. High performance, parallelized C++ is used to compute standard
tools such as radial distribution functions, correlation functions, order
parameters, and clusters, as well as original analysis methods including
potentials of mean force and torque (PMFTs) and local environment matching. The
**freud** library supports
many input formats
and outputs NumPy arrays, enabling integration
with the scientific Python ecosystem for many typical materials science
workflows.

## Resources¶

Reference Documentation: Examples, tutorials, topic guides, and package Python APIs.

Installation Guide: Instructions for installing and compiling

**freud**.freud-users Google Group: Ask questions to the

**freud**user community.GitHub repository: Download the

**freud**source code.Issue tracker: Report issues or request features.

## Citation¶

When using **freud** to process data for publication, please use this citation.

## Installation¶

The easiest ways to install **freud** are using pip:

```
pip install freud-analysis
```

or conda:

```
conda install -c conda-forge freud
```

**freud** is also available via containers for Docker or Singularity. If you need more detailed
information or wish to install **freud** from source, please refer to the
Installation Guide to compile
**freud** from source.

## Examples¶

The **freud** library is called using Python scripts. Many core features are
demonstrated in the freud documentation. The examples come in
the form of Jupyter notebooks, which can also be downloaded from the freud
examples repository or
launched interactively on Binder.
Below is a sample script that computes the radial distribution function for a
simulation run with HOOMD-blue and
saved into a GSD file.

```
import freud
import gsd.hoomd
# Create a freud compute object (RDF is the canonical example)
rdf = freud.density.RDF(bins=50, r_max=5)
# Load a GSD trajectory (see docs for other formats)
traj = gsd.hoomd.open('trajectory.gsd', 'rb')
for frame in traj:
rdf.compute(system=frame, reset=False)
# Get bin centers, RDF data from attributes
r = rdf.bin_centers
y = rdf.rdf
```

## Support and Contribution¶

Please visit our repository on GitHub for the library source code.
Any issues or bugs may be reported at our issue tracker, while questions and discussion can be directed to our user forum.
All contributions to **freud** are welcomed via pull requests!

## Table of Contents¶

### Introduction¶

The **freud** library is a Python package for analyzing particle simulations.
The package is designed to directly use numerical arrays of data, making it easy to use for a wide range of use-cases.
The most common use-case of **freud** is for computing quantities from molecular dynamics simulation trajectories, but it can be used for analyzing any type of particle simulation.
By operating directly on numerical arrays of data, **freud** allows users to parse custom simulation outputs into a suitable structure for input, rather than relying specific file types or data structures.

The core of **freud** is analysis of periodic systems, which are represented through the `freud.box.Box`

class.
The `freud.box.Box`

supports arbitrary triclinic systems for maximum flexibility, and is used throughout the package to ensure consistent treatment of these systems.
The package’s many methods are encapsulated in various *compute classes*, which perform computations and populate class attributes for access.
Of particular note are the various computations based on nearest neighbor finding in order to characterize particle environments.
Such methods are simplified and accelerated through a centralized neighbor finding interface defined in the `freud.locality.NeighborQuery`

family of classes in the `freud.locality`

module of freud.

### Installation¶

#### Installing freud¶

The **freud** library can be installed via conda or pip, or compiled from source.

##### Install via conda¶

The code below will install **freud** from conda-forge.

```
conda install -c conda-forge freud
```

##### Compile from source¶

The following are **required** for installing **freud**:

A C++14-compliant compiler

Python (>=3.6)

NumPy (>=1.14)

Cython (>=0.29.14)

scikit-build (>=0.10.0)

CMake (>=3.6.3)

Note

Depending on the generator you are using, you may require a newer version of CMake. In particular, on Windows Visual Studio 2017 requires at least CMake 3.7.1, while Visual Studio 2019 requires CMake 3.14. For more information on specific generators, see the CMake generator documentation.

The **freud** library uses scikit-build and CMake to handle the build process itself, while the other requirements are required for compiling code in **freud**.
These requirements can be met by installing the following packages from the conda-forge channel:

```
conda install -c conda-forge tbb tbb-devel numpy cython scikit-build cmake
```

All requirements other than TBB can also be installed via the Python Package Index

```
pip install numpy cython scikit-build cmake
```

Wheels for tbb and tbb-devel exist on PyPI, but only for certain operating systems, so your mileage may vary. For non-conda users, we recommend using OS-specific package managers (e.g. Homebrew for Mac) to install TBB. As in the snippets above, it may be necessary to install both both a TBB and a “devel” package in order to get both the headers and the shared libraries.

The code that follows builds **freud** and installs it for all users (append `--user`

if you wish to install it to your user site directory):

```
git clone --recurse-submodules https://github.com/glotzerlab/freud.git
cd freud
python setup.py install
```

You can also build **freud** in place so that you can run from within the folder:

```
# Run tests from the tests directory
python setup.py build_ext --inplace
```

Building **freud** in place has certain advantages, since it does not affect your Python behavior except within the **freud** directory itself (where **freud** can be imported after building).

###### CMake Options¶

The scikit-build tool allows setup.py to accept three different sets of options separated by `--`

, where each set is provided directly to scikit-build, to CMake, or to the code generator of choice, respectively.
For example, the command `python setup.py build_ext --inplace -- -DCOVERAGE=ON -G Ninja -- -j 4`

tell scikit-build to perform an in-place build, it tells CMake to turn on the `COVERAGE`

option and use Ninja for compilation, and it tells Ninja to compile with 4 parallel threads.
For more information on these options, see the scikit-build docs.

Note

The default CMake build configuration for freud is `ReleaseWithDocs`

(not a standard build configuration like `Release`

or `RelWithDebInfo`

).
On installation, `setup.py`

assumes `--build-type=ReleaseWithDocs`

by default if no build type is specified.
Using this build configuration is a workaround for this issue with scikit-build and Cython embedding docstrings.

In addition to standard CMake flags, the following CMake options are available for **freud**:

- --COVERAGE¶
Build the Cython files with coverage support to check unit test coverage.

The **freud** CMake configuration also respects the following environment variables (in addition to standards like `LD_LIBRARY_PATH`

).

- TBBROOT¶
The root directory where TBB is installed. Useful if TBB is installed in a non-standard location or cannot be located for some other reason. This variable is set by the

`tbbvars.sh`

script included with TBB when building from source.- TBB_INCLUDE_DIR¶
The directory where the TBB headers (e.g.

`tbb.h`

) are located. Useful if TBB is installed in a non-standard location or cannot be located for some other reason.

Note

**freud** makes use of git submodules. To manually update git submodules, execute:

```
git submodule update --init --recursive
```

#### Unit Tests¶

The unit tests for **freud** are included in the repository and are configured to be run using the Python `pytest`

library:

```
# Run tests from the tests directory
cd tests
python -m pytest .
```

Note that because **freud** is designed to require installation to run (i.e. it cannot be run directly out of the build directory), importing **freud** from the root of the repository will fail because it will try and import the package folder.
As a result, unit tests must be run from outside the root directory if you wish to test the installed version of **freud**.
If you want to run tests within the root directory, you can instead build **freud** in place:

```
# Run tests from the tests directory
python setup.py build_ext --inplace
```

This build will place the necessary files alongside the **freud** source files so that **freud** can be imported from the root of the repository.

#### Documentation¶

The documentation for **freud** is hosted online at ReadTheDocs.
You may also build the documentation yourself.

##### Building the documentation¶

The following are **required** for building **freud** documentation:

You can install these dependencies using conda:

```
conda install -c conda-forge sphinx sphinx_rtd_theme nbsphinx jupyter_sphinx sphinxcontrib-bibtex
```

or pip:

```
pip install sphinx sphinx-rtd-theme nbsphinx jupyter-sphinx sphinxcontrib-bibtex
```

To build the documentation, run the following commands in the source directory:

```
cd doc
make html
# Then open build/html/index.html
```

To build a PDF of the documentation (requires LaTeX and/or PDFLaTeX):

```
cd doc
make latexpdf
# Then open build/latex/freud.pdf
```

### Quickstart Guide¶

Once you have installed freud, you can start using **freud** with any simulation data that you have on hand.
As an example, we’ll assume that you have run a simulation using the HOOMD-blue and used the `hoomd.dump.gsd`

command to output the trajectory into a file `trajectory.gsd`

.
The GSD file format provides its own convenient Python file reader that offers access to data in the form of NumPy arrays, making it immediately suitable for calculation with **freud**. Many other file readers and data formats are supported, see Reading Simulation Data for freud for a full list and more examples.

We start by reading the data into a NumPy array:

```
import gsd.hoomd
traj = gsd.hoomd.open('trajectory.gsd', 'rb')
```

We can now immediately calculate important quantities.
Here, we will compute the radial distribution function \(g(r)\) using the `freud.density.RDF`

compute class.
Since the radial distribution function is in practice computed as a histogram, we must specify the histogram bin widths and the largest interparticle distance to include in our calculation.
To do so, we simply instantiate the class with the appropriate parameters and then perform a computation on the given data:

```
import freud
rdf = freud.density.RDF(bins=50, r_max=5)
rdf.compute(system=traj[-1])
```

We can now access the data through properties of the `rdf`

object.

```
r = rdf.bin_centers
y = rdf.rdf
```

Many classes in **freud** natively support plotting their data using Matplotlib <https://matplotlib.org/>:

```
import matplotlib as plt
fig, ax = plt.subplots()
rdf.plot(ax=ax)
```

You will note that in the above example, we computed \(g(r)\) only using the final frame of the simulation trajectory, `traj[-1]`

.
However, in many cases, radial distributions and other similar quantities may be noisy in simulations due to the natural fluctuations present.
In general, what we are interested in are *time-averaged* quantities once a system has equilibrated.
To perform such a calculation, we can easily modify our original calculation to take advantage of **freud**’s *accumulation* features.
To accumulate, just add the argument `reset=False`

with a supported compute object (such as histogram-like computations).
Assuming that you have some method for identifying the frames you wish to include in your sample, our original code snippet would be modified as follows:

```
import freud
rdf = freud.density.RDF(bins=50, r_max=5)
for frame in traj:
rdf.compute(frame, reset=False)
```

You can then access the data exactly as we previously did. And that’s it!

Now that you’ve seen a brief example of reading data and computing a radial distribution function, you’re ready to learn more.
If you’d like a complete walkthrough please see the Tutorial.
The tutorial walks through many of the core concepts in **freud** in greater detail, starting with the basics of the simulation systems we analyze and describing the details of the neighbor finding logic in **freud**.
To see specific features of **freud** in action, look through the Examples.
More detailed documentation on specific classes and functions can be found in the API documentation.

### Tutorial¶

This tutorial provides a complete introduction to **freud**.
Rather than attempting to touch on all features in **freud**, it focuses on common core concepts that will help understand how **freud** works with data and exposes computations to the user.
The tutorial begins by introducing the fundamental concepts of periodic systems as implemented in **freud** and the concept of `Compute classes`

, which consitute the primary API for performing calculations with **freud**.
The tutorial then discusses the most common calculation performed in **freud**, finding neighboring points in periodic systems.
The package’s neighbor finding tools are tuned for high performance neighbor finding, which is what enables most of other calculations in **freud**, which typically involve characterizing local environments of points in some way.
The next part of the tutorial discusses the role of histograms in **freud**, focusing on the common features and properties that all histograms share.
Finally, the tutorial includes a few more complete demonstrations of using **freud** that should provide reasonable templates for use with almost any other features in **freud**.

#### Periodic Boundary Conditions¶

The central goal of **freud** is the analysis of simulations performed in periodic boxes.
Periodic boundary conditions are ubiquitous in simulations because they permit the simulation of quasi-infinite systems with minimal computational effort.
As long as simulation systems are sufficiently large, i.e. assuming that points in the system experience correlations over length scales substantially smaller than the system length scale, periodic boundary conditions ensure that the system appears effectively infinite to all points.

In order to consistently define the geometry of a simulation system with periodic boundaries, **freud** defines the `freud.box.Box`

class.
The class encapsulates the concept of a triclinic simulation box in a right-handed coordinate system.
Triclinic boxes are defined as parallelepipeds: three-dimensional polyhedra where every face is a parallelogram.
In general, any such box can be represented by three distinct, linearly independent box vectors.
Enforcing the requirement of right-handedness guarantees that the box can be represented by a matrix of the form

where each column is one of the box vectors.

Note

All **freud** boxes are centered at the origin, so for a given box the
range of possible positions is \([-L/2, L/2)\).

As such, the box is characterized by six parameters: the box vector lengths \(L_x\), \(L_y\), and \(L_z\), and the tilt factors \(xy\), \(xz\), and \(yz\).
The tilt factors are directly related to the angles between the box vectors.
All computations in **freud** are built around this class, ensuring that they naturally handle data from simulations conducted in non-cubic systems.
There is also native support for two-dimensional (2D) systems when setting \(L_z = 0\).

Boxes can be constructed in a variety of ways.
For simple use-cases, one of the factory functions of the `freud.box.Box`

provides the easiest possible interface:

```
# Make a 10x10 square box (for 2-dimensional systems).
freud.box.Box.square(10)
# Make a 10x10x10 cubic box.
freud.box.Box.cube(10)
```

For more complex use-cases, the `freud.box.Box.from_box()`

method provides an interface to create boxes from any object that can easily be interpreted as a box.

```
# Create a 10x10 square box from a list of two items.
freud.box.Box.from_box([10, 10])
# Create a 10x10x10 cubic box from a list of three items.
freud.box.Box.from_box([10, 10, 10])
# Create a triclinic box from a list of six items (including tilt factors).
freud.box.Box.from_box([10, 5, 2, 0.1, 0.5, 0.7])
# Create a triclinic box from a dictionary.
freud.box.Box.from_box(dict(Lx=8, Ly=7, Lz=10, xy=0.5, xz=0.7, yz=0.2))
# Directly call the constructor.
freud.box.Box(Lx=8, Ly=7, Lz=10, xy=0.5, xz=0.7, yz=0.2, dimensions=3)
```

More examples on how boxes can be created may be found in the API documentation of the `Box`

class.

#### Compute Classes¶

Calculations in **freud** are built around the concept of `Compute classes`

, Python objects that encode a given method and expose it through a `compute`

method.
In general, these methods operate on a system composed of a triclinic box and a NumPy array of particle positions.
The box can be provided as any object that can be interpreted as a **freud** box (as demonstrated in the examples above).
We can look at the `freud.order.Hexatic`

order parameter calculator as an example:

```
import freud
positions = ... # Read positions from trajectory file.
op = freud.order.Hexatic(k=6)
op.compute(
system=({'Lx': 5, 'Ly': 5, 'dimensions': 2}, positions),
neighbors=dict(r_max=3)
)
# Plot the value of the order parameter.
from matplotlib import pyplot as plt
plt.hist(np.absolute(op.particle_order))
```

Here, we are calculating the hexatic order parameter, then using Matplotlib to plot.
The `freud.order.Hexatic`

class constructor accepts a single argument `k`

, which represents the periodicity of the calculation.
If you’re unfamiliar with this order parameter, the most important piece of information here is that many compute methods in **freud** require parameters that are provided when the `Compute class`

is constructed.

To calculate the order parameter we call `compute`

, which takes two arguments, a `tuple`

(box, points) and a `dict`

.
We first focus on the first argument.
The `box`

is any object that can be coerced into a `freud.box.Box`

as described in the previous section; in this case, we use a dictionary to specify a square (2-dimensional) box.
The `points`

must be anything that can be coerced into a 2-dimensional NumPy array of shape `(N, 3)`

In general, the points may be provided as anything that can be interpreted as an \(N\times 3\) list of positions; for more details on valid inputs here, see `numpy.asarray()`

.
Note that because the hexatic order parameter is designed for two-dimensional systems, the points must be provided of the form [x, y, 0] (i.e. the z-component must be 0).
We’ll go into more detail about the `(box, points)`

tuple soon, but for now, it’s sufficient to just think of it as specifying the system of points we want to work with.

Now let’s return to the `neighbors`

argument to `compute`

, which is a dictionary is used to determine which particle neighbors to use.
Many computations in **freud** (such as the hexatic order parameter) involve the bonds in the system (for example, the average length of bonds or the average number of bonds a given point has).
However, the concept of a bond is sufficiently variable between different calculations; for instance, should points be considered bonded if they are within a certain distance of each other?
Should every point be considered bonded to a fixed number of other points?

To accommodate this variability, **freud** offers a very general framework by which bonds can be specified, and we’ll go into more details in the next section.
In the example above, we’ve simply informed the `Hexatic`

class that we want it to define bonds as pairs of particles that are less than 3 distance units apart.
We then access the computed order parameter as `op.particle_order`

(we use `np.absolute()`

because the output is a complex number and we just want to see its magnitude).

##### Accessing Computed Properties¶

In general, `Compute classes`

expose their calculations using properties.
Any parameters to the `Compute`

object (e.g. `k`

in the above example) can typically be accessed as soon as the object is constructed:

```
op = freud.order.Hexatic(k=6)
op.k
```

Computed quantities can also be accessed in a similar manner, but only after the `compute`

method is called.
For example:

```
op = freud.order.Hexatic(k=6)
# This will raise an exception.
op.particle_order
op.compute(
system=({'Lx': 5, 'Ly': 5, 'dimensions': 2}, positions),
neighbors=dict(r_max=3)
)
# Now you can access this.
op.particle_order
```

Note

Most (but not all) of **freud**’s `Compute classes`

are Python wrappers
around high-performance implementations in C++. As a result, none of the
data or the computations is actually stored in the Python object. Instead,
the Python object just stores an instance of the C++ object that actually
owns all its data, performs calculations, and returns computed quantities
to the user. Python properties provide a nice way to hide this logic so
that the Python code involves just a few lines.

Compute objects is that they can be used many times to calculate quantities, and the most recently calculated output can be accessed through the property. If you need to perform a series of calculations and save all the data, you can also easily do that:

```
# Recall that lists of length 2 automatically convert to 2D freud boxes.
box = [5, 5]
op = freud.order.Hexatic(k=6)
# Assuming that we have a list of Nx3 NumPy arrays that represents a
# simulation trajectory, we can loop over it and calculate the order
# parameter values in sequence.
trajectory = ... # Read trajectory file into a list of positions by frame.
hexatic_values = []
for points in trajectory:
op.compute(system=(box, points), neighbors=dict(r_max=3))
hexatic_values.append(op.particle_order)
```

To make using **freud** as simple as possible, all `Compute classes`

are designed to return `self`

when compute is called.
This feature enables a very concise *method-chaining* idiom in **freud** where computed properties are accessed immediately:

```
particle_order = freud.order.Hexatic(k=6).compute(
system=(box, points)).particle_order
```

#### Finding Neighbors¶

Now that you’ve been introduced to the basics of interacting with **freud**, let’s dive into the central feature of **freud**: efficiently and flexibly finding neighbors in periodic systems.

##### Problem Statement¶

###### Neighbor-Based Calculations¶

As discussed in the previous section, a central task in many of the computations in **freud** is finding particles’ neighbors.
These calculations typically only involve a limited subset of a particle’s neighbors that are defined as characterizing its local environment.
This requirement is analogous to the force calculations typically performed in molecular dynamics simulations, where a cutoff radius is specified beyond which pair forces are assumed to be small enough to neglect.
Unlike in simulation, though, many analyses call for different specifications than simply selecting all points within a certain distance.

An important example is the calculation of order parameters, which can help characterize phase transitions.
Such parameters can be highly sensitive to the precise way in which neighbors are selected.
For instance, if a hard distance cutoff is imposed in finding neighbors for the hexatic order parameter, a particle may only be found to have five neighbors when it actually has six neighbors except the last particle is slightly outside the cutoff radius.
To accomodate such differences in a flexible manner, **freud** allows users to specify neighbors in a variety of ways.

###### Finding Periodic Neighbors¶

Finding neighbors in periodic systems is significantly more challenging than in aperiodic systems.
To illustrate the difference, consider the figure above, where the black dashed line indicates the boundaries of the system.
If this system were aperiodic, the three nearest neighbors for point 1 would be points 5, 6, and 7.
However, due to periodicity, point 2 is actually closer to point 1 than any of the others if you consider moving straight through the top (or equivalently, the bottom) boundary.
Although many tools provide efficient implementations of algorithms for finding neighbors in aperiodic systems, they seldom generalize to periodic systems.
Even more rare is the ability to work not just in *cubic* periodic systems, which are relatively tractable, but in arbitrary triclinic geometries as described in Periodic Boundary Conditions.
This is precisely the type of calculation **freud** is designed for.

##### Neighbor Querying¶

To understand how `Compute classes`

find neighbors in **freud**, it helps to start by learning about **freud**’s neighbor finding classes directly.
Note that much more detail on this topic is available in the Query API topic guide; in this section we will restrict ourselves to a higher-level overview.
For our demonstration, we will make use of the `freud.locality.AABBQuery`

class, which implements one fast method for periodic neighbor finding.
The primary mode of interfacing with this class (and other neighbor finding classes) is through the `query`

interface.

```
import numpy as np
import freud
# As an example, we randomly generate 100 points in a 10x10x10 cubic box.
L = 10
num_points = 100
# We shift all points into the expected range for freud.
points = np.random.rand(num_points)*L - L/2
box = freud.box.Box.cube(L)
aq = freud.locality.AABBQuery(box, points)
# Now we generate a smaller sample of points for which we want to find
# neighbors based on the original set.
query_points = np.random.rand(num_points/10)*L - L/2
distances = []
# Here, we ask for the 4 nearest neighbors of each point in query_points.
for bond in aq.query(query_points, dict(num_neighbors=4)):
# The returned bonds are tuples of the form
# (query_point_index, point_index, distance). For instance, a bond
# (1, 3, 0.2) would indicate that points[3] was one of the 4 nearest
# neighbors for query_points[1], and that they are separated by a
# distance of 0.2
# (i.e. np.linalg.norm(query_points[1] - points[3]) == 2).
distances.append(bond[2])
avg_distance = np.mean(distances)
```

Let’s dig into this script a little bit.
Our first step is creating a set of 100 points in a cubic box.
Note that the shifting done in the code above could also be accomplished using the `Box.wrap`

method like so: `box.wrap(np.random.rand(num_points)*L)`

.
The result would appear different, because if plotted without considering periodicity, the points would range from `-L/2`

to `L/2`

rather than from 0 to `L`

.
However, these two sets of points would be equivalent in a periodic system.

We then generate an additional set of `query_points`

and ask for neighbors using the `query`

method.
This function accepts two arguments: a set of points, and a `dict`

of **query arguments**.
Query arguments are a central concept in **freud** and represent a complete specification of the set of neighbors to be found.
In general, the most common forms of queries are those requesting either a fixed number of neighbors, as in the example above, or those requesting all neighbors within a specific distance.
For example, if we wanted to rerun the above example but instead find all bonds of length less than or equal to 2, we would simply replace the for loop above with:

```
for bond in aq.query(query_points, dict(r_max=2)):
distances.append(bond[2])
```

Query arguments constitute a powerful method for specifying a query request.
Many query arguments may be combined for more specific purposes.
A common use-case is finding all neighbors within a single set of points (i.e. setting `query_points = points`

in the above example).
In this situation, however, it is typically not useful for a point to find itself as a neighbor since it is trivially the closest point to itself and falls within any cutoff radius.
To avoid this, we can use the `exclude_ii`

query argument:

```
query_points = points
for bond in aq.query(query_points, dict(num_neighbors=4, exclude_ii=True)):
pass
```

The above example will find the 4 nearest neighbors to each point, excepting the point itself. A complete description of valid query arguments can be found in Query API.

##### Neighbor Lists¶

Query arguments provide a simple but powerful language with which to express neighbor finding logic.
Used in the manner shown above, `query`

can be used to express many calculations in a very natural, Pythonic way.
By itself, though, the API shown above is somewhat restrictive because the output of `query`

is a generator.
If you aren’t familiar with generators, the important thing to know is that they can be looped over, *but only once*.
Unlike objects like lists, which you can loop over as many times as you like, once you’ve looped over a generator once, you can’t start again from the beginning.

In the examples above, this wasn’t a problem because we simply iterated over the bonds once for a single calculation.
However, in many practical cases we may need to reuse the set of neighbors multiple times.
A simple solution would be to simply to store the bonds into a list as we loop over them.
However, because this is such a common use-case, **freud** provides its own containers for bonds: the `freud.locality.NeighborList`

.

Queries can easily be used to generate `NeighborList`

objects using their `toNeighborList`

method:

```
query_result = aq.query(query_points, dict(num_neighbors=4, exclude_ii))
nlist = query_result.toNeighborList()
```

The resulting object provides a persistent container for bond data.
Using `NeighborLists`

, our original example might instead look like this:

```
import numpy as np
import freud
L = 10
num_points = 100
points = np.random.rand(num_points)*L - L/2
box = freud.box.Box.cube(L)
aq = freud.locality.AABBQuery(box, points)
query_points = np.random.rand(num_points/10)*L - L/2
distances = []
# Here, we ask for the 4 nearest neighbors of each point in query_points.
query_result = aq.query(query_points, dict(num_neighbors=4)):
nlist = query_result.toNeighborList()
for (i, j) in nlist:
# Note that we have to wrap the bond vector before taking the norm;
# this is the simplest way to compute distances in a periodic system.
distances.append(np.linalg.norm(box.wrap(query_points[i] - points[j])))
avg_distance = np.mean(distances)
```

Note that in the above example we looped directly over the `nlist`

and recomputed distances.
However, the `query_result`

contained information about distances: here’s how we access that through the `nlist`

:

```
assert np.all(nlist.distances == distances)
```

The indices are also accessible through properties, or through a NumPy-like slicing interface:

```
assert np.all(nlist.query_point_indices == nlist[:, 0])
assert np.all(nlist.point_indices == nlist[:, 1])
```

Note that the `query_points`

are always in the first column, while the `points`

are in the second column.
`freud.locality.NeighborList`

objects also store other properties; for instance, they may assign different weights to different bonds.
This feature can be used by, for example, `freud.order.Steinhardt`

, which is typically used for calculating Steinhardt order parameters, a standard tool for characterizing crystalline order.
When provided appropriately weighted neighbors, however, the class instead computes Minkowski structure metrics, which are much more sensitive measures that can differentiate a wider array of crystal structures.

#### Pair Computations¶

Some computations in **freud** do not depend on bonds at all.
For example, `freud.density.GaussianDensity`

creates a “continuous equivalent” of a system of points by placing normal distributions at each point’s location to smear out its position, then summing the value of these distributions at a set of fixed grid points.
This calculation can be quite useful because it allows the application of various analysis tools like fast Fourier transforms, which require regular grids.
For the purposes of this tutorial, however, the importance of this class is that it is an example of a calculation where neighbors are unimportant: the calculation is performed on a per-point basis only.

The much more common pattern in **freud**, though, is that calculations involve the local neighborhoods of points.
To support efficient, flexible computations of such quantities, various `Compute classes`

essentially expose the same API as the query interface demonstrated in the previous section.
These `PairCompute classes`

are designed to mirror the querying functionality of **freud** as closely as possible.

As an example, let’s consider `freud.density.LocalDensity`

, which calculates the density of points in the local neighborhood of each point.
Adapting our code from the previous section, the simplest usage of this class would be as follows:

```
import numpy as np
import freud
L = 10
num_points = 100
points = np.random.rand(num_points)*L - L/2
box = freud.box.Box.cube(L)
# r_max specifies how far to search around each point for neighbors
r_max = 2
# For systems where the points represent, for instance, particles with a
# specific size, the diameter is used to add fractional volumes for
# neighbors that would be overlapping the sphere of radius r_max around
# each point.
diameter = 0.001
ld = freud.density.LocalDensity(r_max, diameter)
ld.compute(system=(box, points))
# Access the density.
ld.density
```

Using the same example system we’ve been working with so far, we’ve now calculated an estimate for the number of points in the neighborhood of each point.
Since we already told the computation how far to search for neighbors based on `r_max`

, all we had to do was pass a `tuple`

`(box, points)`

to compute indicate where the points were located.

##### Binary Systems¶

Imagine that instead of a single set of points, we actually had two different types of points and we were interested in finding the density of one set of points in the vicinity of the other. In that case, we could modify the above calculation as follows:

```
import numpy as np
import freud
L = 10
num_points = 100
points = np.random.rand(num_points)*L - L/2
query_points = np.random.rand(num_points/10)*L - L/2
r_max = 2
diameter = 0.001
ld = freud.density.LocalDensity(r_max, diameter)
ld.compute(system=(box, points), query_points=query_points)
# Access the density.
ld.density
```

The choice of names names here is suggestive of exactly what this calculation is now doing.
Internally, `freud.density.LocalDensity`

will search for all `points`

that are within the cutoff distance `r_max`

of every `query_point`

(essentially using the query interface we introduced previously) and use that to calculate `ld.density`

.
Note that this means that `ld.density`

now contains densities for every `query_point`

, i.e. it is of length 10, not 100.
Moreover, recall that one of the features of the querying API is the specification of whether or not to count particles as their own neighbors.
`PairCompute classes`

will attempt to make an intelligent determination of this for you; if you do not pass in a second set of `query_points`

, they will assume that you are computing with a single set of points and automatically exclude self-neighbors, but otherwise all neighbors will be included.

So far, we have included all points within a fixed radius; however, one might instead wish to consider the density in some shell, such as the density between 1 and 2 distance units away.
To address this need, you could simply adapt the call to `compute`

above as follows:

```
ld.compute(system=(box, points), query_points=query_points,
neighbors=dict(r_max=2, r_min=1))
```

The `neighbors`

argument to `PairCompute`

classes allows users to specify arbitary query arguments, making it possible to easily modify **freud** calculations on-the-fly.
The `neighbors`

argument is actually more general than query arguments you’ve seen so far: if query arguments are not precise enough to specify the exact set of neighbors you want to compute with, you can instead provide a `NeighborList`

directly

```
ld.compute(system=(box, points), query_points=query_points,
neighbors=nlist)
```

This feature allows users essentially arbitrary flexibility to specify the bonds that should be included in any bond-based computation.
A common use-case for this is constructing a `NeighborList`

using `freud.locality.Voronoi`

; Voronoi constructions provide a powerful alternative method of defining neighbor relationships that can improve the accuracy and robustness of certain calculations in **freud**.

You may have noticed in the last example that all the arguments are specified using keyword arguments.
As the previous examples have attempted to show, the `query_points`

argument defines a second set of points to be used when performing calculations on binary systems, while the `neighbors`

argument is how users can specify which neighbors to consider in the calculation.

The `system`

argument is what, to this point, we have been specifying as a `tuple`

`(box, points)`

.
However, we don’t have to use this tuple.
Instead, we can pass in any `freud.locality.NeighborQuery`

, the central class in **freud**’s querying infrastructure.
In fact, you’ve already seen examples of `freud.locality.NeighborQuery`

: the `freud.locality.AABBQuery`

object that we originally used to find neighbors.
There are also a number of other input types that can be converted via `freud.locality.NeighborQuery.from_system()`

, see also Reading Simulation Data for freud.
Since these objects all contain a `freud.box.Box`

and a set of points, they can be directly passed to computations:

```
aq = freud.locality.AABBQuery(box, points)
ld.compute(system=aq, query_points=query_points, neighbors=nlist)
```

For more information on why you might want to use `freud.locality.NeighborQuery`

objects instead of the tuples, see Using freud Efficiently.
For now, just consider this to be a way in which you can simplify your calls to many **freud** computes in one script by storing `(box, points)`

into another objects.

You’ve now covered the most important information needed to use **freud**!
To recap, we’ve discussed how **freud** handles periodic boundary conditions, the structure and usage of `Compute classes`

, and methods for finding and performing calculations with pairs of neighbors.
For more detailed information on specific methods in **freud**, see the Examples page or look at the API documentation for specific modules.

### Examples¶

Examples are provided as Jupyter notebooks in a separate freud-examples repository. These notebooks may be launched interactively on Binder or downloaded and run on your own system. Visualization of data is done via Matplotlib and Bokeh, unless otherwise noted.

#### Key concepts¶

There are a few critical concepts, algorithms, and data structures that are central to all of **freud**.
The `freud.box.Box`

class defines the concept of a periodic simulation box, and the `freud.locality`

module defines methods for finding nearest neighbors of particles.
Since both of these are used throughout **freud**, we recommend reading the Tutorial first, before delving into the workings of specific **freud** analysis modules.

##### freud.box.Box¶

In this notebook, we demonstrate the basic features of the `Box`

class, including wrapping particles back into the box under periodic boundary conditions. For more information, see the introduction to Periodic Boundary Conditions and the `freud.box`

documentation.

###### Creating a Box object¶

Boxes may be constructed explicitly using all arguments. Such construction is useful when performing *ad hoc* analyses involving custom boxes. In general, boxes are assumed to be 3D and orthorhombic unless otherwise specified.

```
[1]:
```

```
import freud.box
# All of the below examples are valid boxes.
box = freud.box.Box(Lx=5, Ly=6, Lz=7, xy=0.5, xz=0.6, yz=0.7, is2D=False)
box = freud.box.Box(1, 3, 2, 0.3, 0.9)
box = freud.box.Box(5, 6, 7)
box = freud.box.Box(5, 6, is2D=True)
box = freud.box.Box(5, 6, xy=0.5, is2D=True)
```

###### From another Box object¶

The simplest case is simply constructing one freud box from another.

**Note that all forms of creating boxes aside from the explicit method above use methods defined within the Box class rather than attempting to overload the constructor itself.**

```
[2]:
```

```
box = freud.box.Box(1, 2, 3)
box2 = freud.box.Box.from_box(box)
print("The original box: \n\t{}".format(box))
print("The copied box: \n\t{}\n".format(box2))
# Boxes are always copied by value, not by reference
box.Lx = 5
print("The original box is modified: \n\t{}".format(box))
print("The copied box is not: \n\t{}\n".format(box2))
# Note, however, that box assignment creates a new object that
# still points to the original box object, so modifications to
# one are visible on the other.
box3 = box2
print("The new copy: \n\t{}".format(box3))
box2.Lx = 2
print("The new copy after the original is modified: \n\t{}".format(box3))
print("The modified original box: \n\t{}".format(box2))
```

```
The original box:
freud.box.Box(Lx=1.0, Ly=2.0, Lz=3.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
The copied box:
freud.box.Box(Lx=1.0, Ly=2.0, Lz=3.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
The original box is modified:
freud.box.Box(Lx=5.0, Ly=2.0, Lz=3.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
The copied box is not:
freud.box.Box(Lx=1.0, Ly=2.0, Lz=3.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
The new copy:
freud.box.Box(Lx=1.0, Ly=2.0, Lz=3.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
The new copy after the original is modified:
freud.box.Box(Lx=2.0, Ly=2.0, Lz=3.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
The modified original box:
freud.box.Box(Lx=2.0, Ly=2.0, Lz=3.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
```

###### From a matrix¶

A box can be constructed directly from the box matrix representation described above using the `Box.from_matrix`

method.

```
[3]:
```

```
# Matrix representation. Note that the box vectors must represent
# a right-handed coordinate system! This translates to requiring
# that the matrix be upper triangular.
box = freud.box.Box.from_matrix([[1, 1, 0], [0, 1, 0.5], [0, 0, 0.5]])
print("This is a 3D box from a matrix: \n\t{}\n".format(box))
# 2D box
box = freud.box.Box.from_matrix([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
print("This is a 2D box from a matrix: \n\t{}\n".format(box))
# Automatic matrix detection using from_box
box = freud.box.Box.from_box([[1, 1, 0], [0, 1, 0.5], [0, 0, 0.5]])
print("The box matrix was automatically detected: \n\t{}\n".format(box))
# Boxes can be numpy arrays as well
import numpy as np
box = freud.box.Box.from_box(np.array([[1, 1, 0], [0, 1, 0.5], [0, 0, 0.5]]))
print("Using a 3x3 numpy array: \n\t{}".format(box))
```

```
This is a 3D box from a matrix:
freud.box.Box(Lx=1.0, Ly=1.0, Lz=0.5, xy=1.0, xz=0.0, yz=1.0, is2D=False)
This is a 2D box from a matrix:
freud.box.Box(Lx=1.0, Ly=1.0, Lz=0.0, xy=0.0, xz=0.0, yz=0.0, is2D=True)
The box matrix was automatically detected:
freud.box.Box(Lx=1.0, Ly=1.0, Lz=0.5, xy=1.0, xz=0.0, yz=1.0, is2D=False)
Using a 3x3 numpy array:
freud.box.Box(Lx=1.0, Ly=1.0, Lz=0.5, xy=1.0, xz=0.0, yz=1.0, is2D=False)
```

###### From a namedtuple or dict¶

A box can be also be constructed from any object that provides an attribute for `Lx`

, `Ly`

, `Lz`

, `xy`

, `xz`

, and `yz`

(or some subset), such as a `namedtuple`

. This method is suitable for passing in box objects constructed by some other program, for example.

```
[4]:
```

```
from collections import namedtuple
MyBox = namedtuple('mybox', ['Lx', 'Ly', 'Lz', 'xy', 'xz', 'yz'])
box = freud.box.Box.from_box(MyBox(Lx=5, Ly=3, Lz=2, xy=0, xz=0, yz=0))
print("Box from named tuple: \n\t{}\n".format(box))
box = freud.box.Box.from_box(MyBox(Lx=5, Ly=3, Lz=0, xy=0, xz=0, yz=0))
print("2D Box from named tuple: \n\t{}".format(box))
```

```
Box from named tuple:
freud.box.Box(Lx=5.0, Ly=3.0, Lz=2.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
2D Box from named tuple:
freud.box.Box(Lx=5.0, Ly=3.0, Lz=0.0, xy=0.0, xz=0.0, yz=0.0, is2D=True)
```

Similarly, construction is also possible using any object that supports key-value indexing, such as a dict.

```
[5]:
```

```
box = freud.box.Box.from_box(dict(Lx=5, Ly=3, Lz=2))
print("Box from dict: \n\t{}".format(box))
```

```
Box from dict:
freud.box.Box(Lx=5.0, Ly=3.0, Lz=2.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
```

###### From a list¶

Finally, boxes can be constructed from any simple iterable that provides the elements in the correct order.

```
[6]:
```

```
box = freud.box.Box.from_box((5, 6, 7, 0.5, 0, 0.5))
print("Box from tuple: \n\t{}\n".format(box))
box = freud.box.Box.from_box([5, 6])
print("2D Box from list: \n\t{}".format(box))
```

```
Box from tuple:
freud.box.Box(Lx=5.0, Ly=6.0, Lz=7.0, xy=0.5, xz=0.0, yz=0.5, is2D=False)
2D Box from list:
freud.box.Box(Lx=5.0, Ly=6.0, Lz=0.0, xy=0.0, xz=0.0, yz=0.0, is2D=True)
```

###### Convenience APIs¶

We also provide convenience constructors for common geometries, namely square (2D) and cubic (3D) boxes.

```
[7]:
```

```
cube_box = freud.box.Box.cube(L=5)
print("Cubic Box: \n\t{}\n".format(cube_box))
square_box = freud.box.Box.square(L=5)
print("Square Box: \n\t{}".format(square_box))
```

```
Cubic Box:
freud.box.Box(Lx=5.0, Ly=5.0, Lz=5.0, xy=0.0, xz=0.0, yz=0.0, is2D=False)
Square Box:
freud.box.Box(Lx=5.0, Ly=5.0, Lz=0.0, xy=0.0, xz=0.0, yz=0.0, is2D=True)
```

###### Export¶

If you want to export or display the box, you can export box objects into their matrix or dictionary representations, which provide completely specified descriptions of the box.

```
[8]:
```

```
cube_box = freud.box.Box.cube(L=5)
cube_box.to_matrix()
```

```
[8]:
```

```
array([[5., 0., 0.],
[0., 5., 0.],
[0., 0., 5.]])
```

```
[9]:
```

```
cube_box.to_dict()
```

```
[9]:
```

```
{'Lx': 5.0,
'Ly': 5.0,
'Lz': 5.0,
'xy': 0.0,
'xz': 0.0,
'yz': 0.0,
'dimensions': 3}
```

###### Using boxes¶

Given a freud box object, you can query it for all its attributes.

```
[10]:
```

```
box = freud.box.Box.from_matrix([[10, 0, 0], [0, 10, 0], [0, 0, 10]])
print("L_x = {}, L_y = {}, L_z = {}, xy = {}, xz = {}, yz = {}".format(
box.Lx, box.Ly, box.Lz, box.xy, box.xz, box.yz))
print("The length vector: {}".format(box.L))
print("The inverse length vector: ({:1.2f}, {:1.2f}, {:1.2f})".format(*[L for L in box.L_inv]))
```

```
L_x = 10.0, L_y = 10.0, L_z = 10.0, xy = 0.0, xz = 0.0, yz = 0.0
The length vector: [10. 10. 10.]
The inverse length vector: (0.10, 0.10, 0.10)
```

Boxes also support converting between fractional and absolute coordinates.

**Note that the origin in real coordinates is defined at the center of the box.** This means the fractional coordinate range \([0, 1]\) maps onto \([-L/2, L/2]\), not \([0, L]\).

```
[11]:
```

```
# Convert from fractional to absolute coordinates.
print(box.make_absolute([[0, 0, 0], [0.5, 0.5, 0.5], [0.8, 0.3, 1]]))
print()
# Convert from fractional to absolute coordinates and back.
print(box.make_fractional(box.make_absolute([[0, 0, 0], [0.5, 0.5, 0.5], [0.8, 0.3, 1]])))
```

```
[[-5. -5. -5.]
[ 0. 0. 0.]
[ 3. -2. 5.]]
[[0. 0. 0. ]
[0.5 0.5 0.5]
[0.8 0.3 1. ]]
```

Finally (and most critically for enforcing periodicity), boxes support wrapping vectors from outside the box into the box. The concept of periodicity and box wrapping is most easily demonstrated visually.

```
[12]:
```

```
# Construct the box and get points for plotting
Lx = Ly = 10
xy = 0.5
box = freud.box.Box.from_matrix([[Lx, xy*Ly, 0], [0, Ly, 0], [0, 0, 0]])
box.plot()
```

```
[12]:
```

```
<matplotlib.axes._subplots.AxesSubplot at 0x7f2729d20518>
```

With periodic boundary conditions, what this actually represents is an infinite set of these boxes tiling space. For example, you can locally picture this box as surrounding by a set of identical boxes.

```
[13]:
```

```
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(9, 6))
box.plot(ax=ax)
for image in [[-1, 0, 0], [1, 0, 0], [0, -1, 0], [0, 1, 0]]:
box.plot(ax=ax, image=image, linestyle='dashed', color='gray')
plt.show()
```

Any particles in the original box will also therefore be seen as existing in all the neighboring boxes.

```
[14]:
```

```
np.random.seed(0)
fractional_coords = np.zeros((5, 3))
fractional_coords[:, :2] = np.random.rand(5, 2)
particles = box.make_absolute(fractional_coords)
```

```
[15]:
```

```
fig, ax = plt.subplots(figsize=(9, 6))
# Plot the points in the original box.
box.plot(ax=ax)
ax.scatter(particles[:, 0], particles[:, 1])
# Plot particles in each of the periodic boxes.
for image in [[-1, 0, 0], [1, 0, 0], [0, -1, 0], [0, 1, 0]]:
box.plot(ax=ax, image=image, linestyle='dashed', color='gray')
particle_images = box.unwrap(particles, image)
ax.scatter(particle_images[:, 0], particle_images[:, 1])
plt.show()
```

Box wrapping takes points in the periodic images of a box, and brings them back into the original box. In this context, that means that if we apply wrap to each of the sets of particles plotted above, they should all overlap.

```
[16]:
```

```
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
images = [[-1, 0, 0], [1, 0, 0], [0, -1, 0], [0, 1, 0]]
# Plot particles in each of the periodic boxes.
for ax, image in zip(axes.flatten(), images):
box.plot(ax=ax)
box.plot(ax=ax, image=image, linestyle='dashed', color='gray')
particle_images = box.unwrap(particles, image)
ax.scatter(particle_images[:, 0],
particle_images[:, 1],
label='Images')
wrapped_particle_images = box.wrap(particle_images)
ax.scatter(wrapped_particle_images[:, 0],
wrapped_particle_images[:, 1],
label='Wrapped')
ax.tick_params(axis="both", which="both", labelsize=14)
ax.legend()
plt.show()
```

##### freud.locality.PeriodicBuffer: Unit Cell RDF¶

The `PeriodicBuffer`

class is meant to replicate points beyond a single image while respecting box periodicity. This example demonstrates how we can use this to compute the radial distribution function from a sample crystal’s unit cell.

```
[1]:
```

```
%matplotlib inline
import freud
import numpy as np
import matplotlib.pyplot as plt
```

Here, we create a box to represent the unit cell and put two points inside. We plot the box and points below.

```
[2]:
```

```
box = freud.box.Box(Lx=2, Ly=2, xy=np.sqrt(1/3), is2D=True)
points = np.array([[-0.5, -0.5, 0], [0.5, 0.5, 0]])
system = freud.AABBQuery(box, points)
system.plot(ax=plt.gca())
plt.show()
```

Next, we create a `PeriodicBuffer`

instance and have it compute the “buffer” points that lie outside the first periodicity. These positions are stored in the `buffer_points`

attribute. The corresponding `buffer_ids`

array gives a mapping from the index of the buffer particle to the index of the particle it was replicated from, in the original array of `points`

. Finally, the `buffer_box`

attribute returns a larger box, expanded from the original box to contain the replicated points.

```
[3]:
```

```
pbuff = freud.locality.PeriodicBuffer()
pbuff.compute(system=(box, points), buffer=6, images=True)
print(pbuff.buffer_points[:10], '...')
```

```
[[ 0.65470022 1.5 0. ]
[ 1.80940032 3.5 0. ]
[ 2.96410179 5.5 0. ]
[-3.96410131 -6.5 0. ]
[-2.80940104 -4.49999952 0. ]
[-1.65470016 -2.50000048 0. ]
[ 1.50000024 -0.5 0. ]
[ 2.65470076 1.5 0. ]
[ 3.80940032 3.5 0. ]
[ 4.96410179 5.5 0. ]] ...
```

Below, we plot the original unit cell and the replicated buffer points and buffer box.

```
[4]:
```

```
system.plot(ax=plt.gca())
plt.scatter(pbuff.buffer_points[:, 0], pbuff.buffer_points[:, 1])
pbuff.buffer_box.plot(ax=plt.gca(), linestyle='dashed', color='gray')
plt.show()
```

Finally, we can plot the radial distribution function (RDF) of this replicated system, using a value of `r_max`

that is larger than the size of the original box. This allows us to see the interaction of the original points with their replicated neighbors from the buffer.

```
[5]:
```

```
rdf = freud.density.RDF(bins=250, r_max=5)
rdf.compute(system=(pbuff.buffer_box, pbuff.buffer_points), query_points=points)
rdf.plot(ax=plt.gca())
plt.show()
```

##### freud.locality.Voronoi¶

The `freud.locality.Voronoi`

class uses voro++ to compute the Voronoi diagram of a set of points, **while respecting periodic boundary conditions** (which are not handled by `scipy.spatial.Voronoi`

, documentation).

These examples are two-dimensional (with \(z=0\) for all particles) for simplicity, but the `Voronoi`

class works for both 2D and 3D data.

```
[1]:
```

```
import numpy as np
import freud
import matplotlib
import matplotlib.pyplot as plt
```

First, we generate some sample points.

```
[2]:
```

```
points = np.array([
[-0.5, -0.5, 0],
[0.5, -0.5, 0],
[-0.5, 0.5, 0],
[0.5, 0.5, 0]])
plt.scatter(points[:, 0], points[:, 1])
plt.title('Points')
plt.xlim((-1, 1))
plt.ylim((-1, 1))
plt.gca().set_aspect('equal')
plt.show()
```

Now we create a box and a `Voronoi`

compute object.

```
[3]:
```

```
L = 2
box = freud.box.Box.square(L)
voro = freud.locality.Voronoi()
```

Next, we use the `compute`

method to determine the Voronoi polytopes (cells) and the `polytopes`

property to return their coordinates. Note that we use `freud`

’s *method chaining* here, where a compute method returns the compute object.

```
[4]:
```

```
cells = voro.compute((box, points)).polytopes
print(cells)
```

```
[array([[-1., -1., 0.],
[ 0., -1., 0.],
[ 0., 0., 0.],
[-1., 0., 0.]]), array([[ 0., -1., 0.],
[ 1., -1., 0.],
[ 1., 0., 0.],
[ 0., 0., 0.]]), array([[-1., 0., 0.],
[ 0., 0., 0.],
[ 0., 1., 0.],
[-1., 1., 0.]]), array([[0., 0., 0.],
[1., 0., 0.],
[1., 1., 0.],
[0., 1., 0.]])]
```

The `Voronoi`

class has built-in plotting methods for 2D systems.

```
[5]:
```

```
plt.figure()
ax = plt.gca()
voro.plot(ax=ax)
ax.scatter(points[:, 0], points[:, 1], s=10, c='k')
plt.show()
```

This also works for more complex cases, such as this hexagonal lattice.

```
[6]:
```

```
def hexagonal_lattice(rows=3, cols=3, noise=0, seed=None):
if seed is not None:
np.random.seed(seed)
# Assemble a hexagonal lattice
points = []
for row in range(rows*2):
for col in range(cols):
x = (col + (0.5 * (row % 2)))*np.sqrt(3)
y = row*0.5
points.append((x, y, 0))
points = np.asarray(points)
points += np.random.multivariate_normal(mean=np.zeros(3), cov=np.eye(3)*noise, size=points.shape[0])
# Set z=0 again for all points after adding Gaussian noise
points[:, 2] = 0
# Wrap the points into the box
box = freud.box.Box(Lx=cols*np.sqrt(3), Ly=rows, is2D=True)
points = box.wrap(points)
return box, points
```

```
[7]:
```

```
# Compute the Voronoi diagram and plot
box, points = hexagonal_lattice()
voro = freud.locality.Voronoi()
voro.compute((box, points))
voro
```

```
[7]:
```

For noisy data, we see that the Voronoi diagram can change substantially. We perturb the positions with 2D Gaussian noise. Coloring by the number of sides of each Voronoi cell, we can see patterns in the defects: 5-gons and 7-gons tend to pair up.

```
[8]:
```

```
# Compute the Voronoi diagram
box, points = hexagonal_lattice(rows=12, cols=8, noise=0.03, seed=2)
voro = freud.locality.Voronoi()
voro.compute((box, points))
# Plot Voronoi with points and a custom cmap
plt.figure()
ax = plt.gca()
voro.plot(ax=ax, cmap='RdBu')
ax.scatter(points[:, 0], points[:, 1], s=2, c='k')
plt.show()
```

We can also compute the volumes of the Voronoi cells. Here, we plot them as a histogram:

```
[9]:
```

```
plt.hist(voro.volumes)
plt.title('Voronoi cell volumes')
plt.show()
```

The `Voronoi`

class also computes a `freud.locality.NeighborList`

, where particles are neighbors if they share an edge in the Voronoi diagram. The `NeighborList`

effectively represents the bonds in the Delaunay triangulation. The neighbors are **weighted** by the length (in 2D) or area (in 3D) between them. The neighbor weights are stored in `voro.nlist.weights`

.

```
[10]:
```

```
nlist = voro.nlist
line_data = np.asarray([[points[i],
points[i] + box.wrap(points[j] - points[i])]
for i, j in nlist])[:, :, :2]
line_collection = matplotlib.collections.LineCollection(line_data, alpha=0.2)
plt.figure()
ax = plt.gca()
voro.plot(ax=ax, cmap='RdBu')
ax.add_collection(line_collection)
plt.show()
```

#### Analysis Modules¶

These introductory examples showcase the functionality of specific modules in **freud**, showing how they can be used to perform specific types of analyses of simulations.

##### freud.cluster.Cluster and freud.cluster.ClusterProperties¶

The `freud.cluster`

module determines clusters of points and computes cluster quantities like centers of mass, gyration tensors, and radii of gyration. The example below generates random points, and shows that they form clusters. This case is two-dimensional (with \(z=0\) for all particles) for simplicity, but the cluster module works for both 2D and 3D simulations.

```
[1]:
```

```
import numpy as np
import freud
import matplotlib.pyplot as plt
```

First, we generate a box and random points to cluster.

```
[2]:
```

```
box = freud.Box.square(L=6)
points = np.empty(shape=(0, 2))
for center_point in [(-1.8, 0), (1.5, 1.5), (-0.8, -2.8), (1.5, 0.5)]:
points = np.concatenate(
(points, np.random.multivariate_normal(mean=center_point, cov=0.08*np.eye(2), size=(100,))))
points = np.hstack((points, np.zeros((points.shape[0], 1))))
points = box.wrap(points)
system = freud.AABBQuery(box, points)
system.plot(ax=plt.gca(), s=10)
plt.title('Raw points before clustering', fontsize=20)
plt.gca().tick_params(axis='both', which='both', labelsize=14, size=8)
plt.show()
```

Now we create a box and a cluster compute object.

```
[3]:
```

```
cl = freud.cluster.Cluster()
```

Next, we use the `computeClusters`

method to determine clusters and the `clusterIdx`

property to return their identities. Note that we use `freud`

’s *method chaining* here, where a compute method returns the compute object.

```
[4]:
```

```
cl.compute(system, neighbors={'r_max': 1.0})
print(cl.cluster_idx)
```

```
[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
```

```
[5]:
```

```
fig, ax = plt.subplots(1, 1, figsize=(9, 6))
for cluster_id in range(cl.num_clusters):
cluster_system = freud.AABBQuery(system.box, system.points[cl.cluster_keys[cluster_id]])
cluster_system.plot(ax=ax, s=10, label="Cluster {}".format(cluster_id))
print("There are {} points in cluster {}.".format(len(cl.cluster_keys[cluster_id]), cluster_id))
ax.set_title('Clusters identified', fontsize=20)
ax.legend(loc='best', fontsize=14)
ax.tick_params(axis='both', which='both', labelsize=14, size=8)
plt.show()
```

```
There are 200 points in cluster 0.
There are 100 points in cluster 1.
There are 100 points in cluster 2.
```

We may also compute the clusters’ centers of mass and gyration tensor using the `ClusterProperties`

class.

```
[6]:
```

```
clp = freud.cluster.ClusterProperties()
clp.compute(system, cl.cluster_idx);
```

Plotting these clusters with their centers of mass, with size proportional to the number of clustered points:

```
[7]:
```

```
fig, ax = plt.subplots(1, 1, figsize=(9, 6))
for i in range(cl.num_clusters):
cluster_system = freud.AABBQuery(system.box, system.points[cl.cluster_keys[i]])
cluster_system.plot(ax=ax, s=10, label="Cluster {}".format(i))
for i, c in enumerate(clp.centers):
ax.scatter(c[0], c[1], s=len(cl.cluster_keys[i]),
label="Cluster {} Center".format(i))
plt.title('Center of mass for each cluster', fontsize=20)
plt.legend(loc='best', fontsize=14)
plt.gca().tick_params(axis='both', which='both', labelsize=14, size=8)
plt.gca().set_aspect('equal')
plt.show()
```

The 3x3 gyration tensors \(G\) can also be computed for each cluster. For this two-dimensional case, the \(z\) components of the gyration tensor are zero. The gyration tensor can be used to determine the principal axes of the cluster and radius of gyration along each principal axis. Here, we plot the gyration tensor’s eigenvectors with length corresponding to the square root of the eigenvalues (the singular values).

```
[8]:
```

```
fig, ax = plt.subplots(1, 1, figsize=(9, 6))
for i in range(cl.num_clusters):
cluster_system = freud.AABBQuery(system.box, system.points[cl.cluster_keys[i]])
cluster_system.plot(ax=ax, s=10, label="Cluster {}".format(i))
for i, c in enumerate(clp.centers):
ax.scatter(c[0], c[1], s=len(cl.cluster_keys[i]),
label="Cluster {} Center".format(i))
for cluster_id in range(cl.num_clusters):
com = clp.centers[cluster_id]
G = clp.gyrations[cluster_id]
evals, evecs = np.linalg.eig(G[:2, :2])
arrows = np.sqrt(evals) * evecs
for arrow in arrows.T:
plt.arrow(com[0], com[1], arrow[0], arrow[1], width=0.05, color='k')
plt.title('Eigenvectors of the gyration tensor for each cluster', fontsize=20)
plt.legend(loc='best', fontsize=14)
ax.tick_params(axis='both', which='both', labelsize=14, size=8)
ax.set_aspect('equal')
plt.show()
```

##### freud.diffraction.DiffractionPattern¶

The `freud.diffraction.DiffractionPattern`

class computes a diffraction pattern, which is a 2D image of the static structure factor \(S(\vec{k})\) of a set of points.

```
[1]:
```

```
import freud
import matplotlib.pyplot as plt
import numpy as np
import rowan
```

First, we generate a sample system, a face-centered cubic crystal with some noise.

```
[2]:
```

```
box, points = freud.data.UnitCell.fcc().generate_system(num_replicas=10, sigma_noise=0.02)
```

Now we create a `DiffractionPattern`

compute object.

```
[3]:
```

```
dp = freud.diffraction.DiffractionPattern(grid_size=512, output_size=512)
```

Next, we use the `compute`

method and plot the result. We use a view orientation with the identity quaternion `[1, 0, 0, 0]`

so the view is aligned down the z-axis.

```
[4]:
```

```
fig, ax = plt.subplots(figsize=(4, 4), dpi=150)
dp.compute((box, points), view_orientation=[1, 0, 0, 0])
dp.plot(ax)
plt.show()
```

We can also use a random quaternion for the view orientation to see what the diffraction looks like from another axis.

```
[5]:
```

```
fig, ax = plt.subplots(figsize=(4, 4), dpi=150)
np.random.seed(0)
view_orientation = rowan.random.rand()
dp.compute((box, points), view_orientation=view_orientation)
print('Looking down the axis:', rowan.rotate(view_orientation, [0, 0, 1]))
dp.plot(ax)
plt.show()
```

```
Looking down the axis: [0.75707404 0.33639217 0.56007071]
```

The `DiffractionPattern`

object also provides \(\vec{k}\) vectors in the original 3D space and the magnitudes of \(k_x\) and \(k_y\) in the 2D projection along the view axis.

```
[6]:
```

```
print('Magnitudes of k_x and k_y along the plot axes:')
print(dp.k_values[:5], '...', dp.k_values[-5:])
```

```
Magnitudes of k_x and k_y along the plot axes:
[-25.6 -25.5 -25.4 -25.3 -25.2] ... [25.1 25.2 25.3 25.4 25.5]
```

```
[7]:
```

```
print('3D k-vectors corresponding to each pixel of the diffraction image:')
print('Array shape:', dp.k_vectors.shape)
print('Center value: k =', dp.k_vectors[dp.output_size//2, dp.output_size//2, :])
print('Top-left value: k =', dp.k_vectors[0, 0, :])
```

```
3D k-vectors corresponding to each pixel of the diffraction image:
Array shape: (512, 512, 3)
Center value: k = [0. 0. 0.]
Top-left value: k = [ 19.03669015 -29.77826902 -7.84723661]
```

We can also measure the diffraction of a random system (note: this is an ideal gas, not a liquid-like system, because the particles have no volume exclusion or repulsion). Note that the peak at \(\vec{k} = 0\) persists. The diffraction pattern returned by this class is normalized by dividing by \(N^2\), so \(S(\vec{k}=0) = 1\) after normalization.

```
[8]:
```

```
box, points = freud.data.make_random_system(box_size=10, num_points=10000)
fig, ax = plt.subplots(figsize=(4, 4), dpi=150)
dp.compute((box, points))
dp.plot(ax)
plt.show()
```

##### freud.density.CorrelationFunction¶

###### Orientational Ordering in 2D¶

The `freud.density`

module is intended to compute a variety of quantities that relate spatial distributions of particles with other particles. This example shows how correlation functions can be used to measure orientational order in 2D.

```
[1]:
```

```
import numpy as np
import freud
import matplotlib.pyplot as plt
import matplotlib.cm
from matplotlib.colors import Normalize
```

This helper function will make plots of the data we generate in this example.

```
[2]:
```

```
def plot_data(title, points, angles, values, box, cf, s=200):
cmap = matplotlib.cm.viridis
norm = Normalize(vmin=-np.pi/4, vmax=np.pi/4)
plt.figure(figsize=(16, 6))
plt.subplot(121)
for point, angle, value in zip(points, angles, values):
plt.scatter(point[0], point[1], marker=(4, 0, np.rad2deg(angle)+45),
edgecolor='k', c=[cmap(norm(angle))], s=s)
plt.title(title)
plt.gca().set_xlim([-box.Lx/2, box.Lx/2])
plt.gca().set_ylim([-box.Ly/2, box.Ly/2])
plt.gca().set_aspect('equal')
sm = plt.cm.ScalarMappable(cmap='viridis', norm=norm)
sm.set_array(angles)
plt.colorbar(sm)
plt.subplot(122)
plt.title('Orientation Spatial Autocorrelation Function')
cf.plot(ax=plt.gca())
plt.xlabel(r'$r$')
plt.ylabel(r'$C(r)$')
plt.show()
```

First, let’s generate a 2D structure with perfect orientational order and slight positional disorder (the particles are not perfectly on a grid, but they are perfectly aligned). The color of the particles corresponds to their angle of rotation, so all the particles will be the same color to begin with.

We create a `freud.density.CorrelationFunction`

object to compute the correlation functions. Given a particle orientation \(\theta\), we compute its complex orientation value (the quantity we are correlating) as \(s = e^{4i\theta}\), to account for the fourfold symmetry of the particles. We will compute the correlation function \(C(r) = \left\langle s^*_1(0) \cdot s_2(r) \right\rangle\) by taking the average over all particle pairs and binning the results into a histogram by the
distance \(r\) between the particles.

When we compute the correlations between particles, the complex conjugate of the `values`

array is used internally for the query points. This way, if \(\theta_1\) is close to \(\theta_2\), then we get \(\left(e^{4i\theta_1}\right)^* \cdot \left(e^{4i\theta_2}\right) = e^{4i(\theta_2-\theta_1)} \approx e^0 = 1\).

This system has perfect spatial correlation of the particle orientations, so we see \(C(r) = 1\) for all values of \(r\).

```
[3]:
```

```
def make_particles(L, repeats):
uc = freud.data.UnitCell.square()
return uc.generate_system(num_replicas=repeats, scale=L/repeats, sigma_noise=5e-3*L)
# Make a small system
box, points = make_particles(L=5, repeats=20)
# All the particles begin with their orientation at 0
angles = np.zeros(len(points))
values = np.array(np.exp(angles * 4j))
# Create the CorrelationFunction compute object and compute the correlation function
cf = freud.density.CorrelationFunction(bins=25, r_max=box.Lx/2.01)
cf.compute(system=(box, points), values=values,
query_points=points, query_values=values)
plot_data('Particles before introducing Orientational Disorder',
points, angles, values, box, cf)
```

Now we will generate random angles from \(-\frac{\pi}{4}\) to \(\frac{\pi}{4}\), which orients our squares randomly. The four-fold symmetry of the squares means that the space of unique angles is restricted to a range of \(\frac{\pi}{2}\). Again, we compute a complex value for each particle, \(s = e^{4i\theta}\).

Because we have purely random orientations, we expect no spatial correlations in the plot above. As we see, \(C(r) \approx 0\) for all \(r\).

```
[4]:
```

```
# Change the angles to values randomly drawn from a uniform distribution
angles = np.random.uniform(-np.pi/4, np.pi/4, size=len(points))
values = np.exp(angles * 4j)
# Recompute the correlation functions
cf.compute(system=(box, points), values=values,
query_points=points, query_values=values)
plot_data('Particles with Orientational Disorder',
points, angles, values, box, cf)
```

The plot below shows what happens when we intentionally introduce a correlation length by adding a spatial pattern to the particle orientations. At short distances, the correlation is very high. As \(r\) increases, the oppositely-aligned part of the pattern some distance away causes the correlation to drop.

```
[5]:
```

```
# Use angles that vary spatially in a pattern
angles = np.pi/4 * np.cos(2*np.pi*points[:, 0]/box.Lx)
values = np.exp(angles * 4j)
# Recompute the correlation functions
cf.compute(system=(box, points), values=values,
query_points=points, query_values=values)
plot_data('Particles with Spatially Correlated Orientations',
points, angles, values, box, cf)
```

In the larger system shown below, we see the spatial autocorrelation rise and fall with damping oscillations.

```
[6]:
```

```
# Make a large system
box, points = make_particles(L=10, repeats=40)
# Use angles that vary spatially in a pattern
angles = np.pi/4 * np.cos(8*np.pi*points[:, 0]/box.Lx)
values = np.exp(angles * 4j)
# Create a CorrelationFunction compute object
cf = freud.density.CorrelationFunction(bins=25, r_max=box.Lx/2.01)
cf.compute(system=(box, points), values=values,
query_points=points, query_values=values)
plot_data('Larger System with Spatially Correlated Orientations',
points, angles, values, box, cf, s=80)
```

##### freud.density.GaussianDensity¶

The `freud.density`

module is intended to compute a variety of quantities that relate spatial distributions of particles with other particles. In this notebook, we demonstrate `freud`

’s Gaussian density calculation, which provides a way to interpolate particle configurations onto a regular grid in a meaningful way that can then be processed by other algorithms that require regularity, such as a Fast Fourier Transform.

```
[1]:
```

```
import numpy as np
from scipy import stats
import freud
import matplotlib.pyplot as plt
```

To illustrate the basic concept, consider a toy example: a simple set of point particles with unit mass on a line. For analytical purposes, the standard way to accomplish this would be using Dirac delta functions.

```
[2]:
```

```
n_p = 10000
np.random.seed(129)
x = np.linspace(0, 1, n_p)
y = np.zeros(n_p)
points = np.random.rand(10)
y[(points*n_p).astype('int')] = 1
plt.plot(x, y);
plt.show()
```

However, delta functions can be cumbersome to work with, so we might instead want to smooth out these particles. One option is to instead represent particles as Gaussians centered at the location of the points. In that case, the total particle density at any point in the interval \([0, 1]\) represented above would be based on the sum of the densities of those Gaussians at those points.

```
[3]:
```

```
# Note that we use a Gaussian with a small standard deviation
# to emphasize the differences on this small scale
dists = [stats.norm(loc=i, scale=0.1) for i in points]
y_gaussian = 0
for dist in dists:
y_gaussian += dist.pdf(x)
plt.plot(x, y_gaussian);
plt.show()
```

The goal of the GaussianDensity class is to perform the same interpolation for points on a 2D or 3D grid, accounting for Box periodicity.

```
[4]:
```

```
N = 1000 # Number of points
L = 10 # Box length
box, points = freud.data.make_random_system(L, N, is2D=True, seed=0)
aq = freud.AABBQuery(box, points)
gd = freud.density.GaussianDensity(L*L, L/3, 1)
gd.compute(aq)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
aq.plot(ax=axes[0])
gd.plot(ax=axes[1])
plt.show()
```

The effects are much more striking if we explicitly construct our points to be centered at certain regions.

```
[5]:
```

```
N = 1000 # Number of points
L = 10 # Box length
box = freud.box.Box.square(L)
centers = np.array([[L/4, L/4, 0],
[-L/4, L/4, 0],
[L/4, -L/4, 0],
[-L/4, -L/4, 0]])
points = []
for center in centers:
points.append(np.random.multivariate_normal(center, cov=np.diag([1, 1, 0]), size=(int(N/4),)))
points = box.wrap(np.concatenate(points))
aq = freud.AABBQuery(box, points)
gd = freud.density.GaussianDensity(L*L, L/3, 1)
gd.compute(aq)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
aq.plot(ax=axes[0])
gd.plot(ax=axes[1])
plt.show()
```

##### freud.density.LocalDensity¶

The `freud.density`

module is intended to compute a variety of quantities that relate spatial distributions of particles with other particles. In this notebook, we demonstrate `freud`

’s local density calculation, which can be used to characterize the particle distributions in some systems. In this example, we consider a toy example of calculating the particle density in the vicinity of a set of other points. This can be visualized as, for example, billiard balls on a table with certain
regions of the table being stickier than others. In practice, this method could be used for analyzing, *e.g*, binary systems to determine how densely one species packs close to the surface of the other.

```
[1]:
```

```
import numpy as np
import freud
import matplotlib.pyplot as plt
from matplotlib import patches
```

```
[2]:
```

```
# Define some helper plotting functions.
def add_patches(ax, points, radius=1, fill=False, color="#1f77b4", ls="solid", lw=None):
"""Add set of points as patches with radius to the provided axis"""
for pt in points:
p = patches.Circle(pt, fill=fill, linestyle=ls, radius=radius,
facecolor=color, edgecolor=color, lw=lw)
ax.add_patch(p)
def plot_lattice(box, points, radius=1, ls="solid", lw=None):
"""Helper function for plotting points on a lattice."""
fig, ax = plt.subplots(1, 1, figsize=(9, 9))
box.plot(ax=ax)
add_patches(ax, points, radius, ls=ls, lw=lw)
return fig, ax
```

Let us consider a set of regions on a square lattice.

```
[3]:
```

```
area = 2
radius = np.sqrt(area/np.pi)
spot_area = area*100
spot_radius = np.sqrt(spot_area/np.pi)
num = 6
scale = num*4
uc = freud.data.UnitCell(freud.Box.square(1), [[0.5, 0.5, 0]])
box, spot_centers = uc.generate_system(num, scale=scale)
fig, ax = plot_lattice(box, spot_centers, spot_radius, ls="dashed", lw=2.5)
plt.tick_params(axis="both", which="both", labelsize=14)
plt.show()
```

Now let’s add a set of points to this box. Points are added by drawing from a normal distribution centered at each of the regions above. For demonstration, we will assume that each region has some relative “attractiveness,” which is represented by the covariance in the normal distributions used to draw points. Specifically, as we go up and to the right, the covariance increases proportional to the distance from the lower right corner of the box.

```
[4]:
```

```
points = []
fractional_distances_to_corner = np.linalg.norm(box.make_fractional(spot_centers), axis=-1)
cov_basis = 20 * fractional_distances_to_corner
for i, p in enumerate(spot_centers):
np.random.seed(i)
cov = cov_basis[i]*np.diag([1, 1, 0])
points.append(
np.random.multivariate_normal(p, cov, size=(50,)))
points = box.wrap(np.concatenate(points))
```

```
[5]:
```

```
fig, ax = plot_lattice(box, spot_centers, spot_radius, ls="dashed", lw=2.5)
plt.tick_params(axis="both", which="both", labelsize=14)
add_patches(ax, points, radius, True, 'k', lw=None)
plt.show()
```

We see that the density decreases as we move up and to the right. In order to compute the actual densities, we can leverage the `LocalDensity`

class. The class allows you to specify a set of query points around which the number of other points is computed. These other points can, but need not be, distinct from the query points. In our case, we want to use the blue regions as our query points with the small black dots as our data points.

When we construct the `LocalDensity`

class, there are two arguments. The first is the radius from the query points within which particles should be included in the query point’s counter. The second is the circumsphere diameter of the **data points**, not the query points. This distinction is critical for getting appropriate density values, since these values are used to actually check cutoffs and calculate the density.

```
[6]:
```

```
density = freud.density.LocalDensity(spot_radius, radius)
density.compute(system=(box, points), query_points=spot_centers);
```

```
[7]:
```

```
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
for i, data in enumerate([density.num_neighbors, density.density]):
poly = np.poly1d(np.polyfit(cov_basis, data, 1))
axes[i].tick_params(axis="both", which="both", labelsize=14)
axes[i].scatter(cov_basis, data)
x = np.linspace(*axes[i].get_xlim(), 30)
axes[i].plot(x, poly(x), label="Best fit")
axes[i].set_xlabel("Covariance", fontsize=16)
axes[0].set_ylabel("Number of neighbors", fontsize=16);
axes[1].set_ylabel("Density", fontsize=16);
plt.show()
```

As expected, we see that increasing the variance in the number of points centered at a particular query point decreases the total density at that point. The trend is noisy since we are randomly sampling possible positions, but the general behavior is clear.

##### freud.density.RDF: Accumulating g(r) for a Fluid¶

The `freud.density`

module is intended to compute a variety of quantities that relate spatial distributions of particles with other particles. This example demonstrates the calculation of the radial distribution function \(g(r)\) for a fluid, averaged over multiple frames.

```
[1]:
```

```
import numpy as np
import freud
import matplotlib.pyplot as plt
data_path = "data/phi065"
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
def plot_rdf(box_arr, points_arr, prop, r_max=10, bins=100, label=None, ax=None):
"""Helper function for plotting RDFs."""
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(12, 8))
ax.set_title(prop, fontsize=16)
rdf = freud.density.RDF(bins, r_max)
for box, points in zip(box_arr, points_arr):
rdf.compute(system=(box, points), reset=False)
if label is not None:
ax.plot(rdf.bin_centers, getattr(rdf, prop), label=label)
ax.legend()
else:
ax.plot(rdf.bin_centers, getattr(rdf, prop))
return ax
```

Here, we show the difference between the RDF of one frame and an accumulated (averaged) RDF from several frames. Including more frames makes the plot smoother.

```
[2]:
```

```
# Compute the RDF for the last frame
box_arr = [box_data[-1].tolist()]
pos_arr = [pos_data[-1]]
ax = plot_rdf(box_arr, pos_arr, 'rdf', label='One frame')
# Compute the RDF for the last 20 frames
box_arr = [box.tolist() for box in box_data[-20:]]
pos_arr = pos_data[-20:]
ax = plot_rdf(box_arr, pos_arr, 'rdf', label='Last 20 frames', ax=ax)
plt.show()
```

The difference between `accumulate`

(which should be called on a series of frames) and `compute`

(meant for a single frame) is more striking for smaller bin sizes, which are statistically noisier.

```
[3]:
```

```
# Compute the RDF for the last frame
box_arr = [box_data[-1].tolist()]
pos_arr = [pos_data[-1]]
ax = plot_rdf(box_arr, pos_arr, 'rdf', bins=1000, label='One frame')
# Compute the RDF for the last 20 frames
box_arr = [box.tolist() for box in box_data[-20:]]
pos_arr = pos_data[-20:]
ax = plot_rdf(box_arr, pos_arr, 'rdf', bins=1000, label='Last 20 frames', ax=ax)
plt.show()
```

##### freud.density.RDF: Choosing Bin Widths¶

The `freud.density`

module is intended to compute a variety of quantities that relate spatial distributions of particles with other particles. This example demonstrates the calculation of the radial distribution function \(g(r)\) using different bin sizes.

```
[1]:
```

```
import numpy as np
import freud
import matplotlib.pyplot as plt
```

```
[2]:
```

```
# Define some helper plotting functions.
def plot_rdf(box, points, prop, r_max=3.5, bins_array=[20, 75, 3000]):
"""Helper function for plotting RDFs."""
fig, axes = plt.subplots(1, len(bins_array), figsize=(16, 3))
for i, bins in enumerate(bins_array):
rdf = freud.density.RDF(bins, r_max)
rdf.compute(system=(box, points))
axes[i].plot(rdf.bin_centers, getattr(rdf, prop))
axes[i].set_title("Bin width: {:.3f}".format(r_max/bins), fontsize=16)
plt.show()
```

To start, we construct and visualize a set of points sitting on a simple square lattice.

```
[3]:
```

```
box, points = freud.data.UnitCell.square().generate_system(5, scale=2)
aq = freud.AABBQuery(box, points)
aq.plot(ax=plt.gca())
plt.show()
```

If we try to compute the RDF directly from this, we will get something rather uninteresting since we have a perfect crystal. Indeed, we will observe that as we bin more and more finely, we approach the true behavior of the RDF for perfect crystals, which is a simple delta function.

```
[4]:
```

```
plot_rdf(box, points, 'rdf')
```

In these RDFs, we see two sharply defined peaks, with the first corresponding to the nearest neighbors on the lattice (which are all at a distance 2 from each other), and the second, smaller peak caused by the particles on the diagonal (which sit at distance \(\sqrt{2^2+2^2} \approx 2.83\).

However, in more realistic systems, we expect that the lattice will not be perfectly formed. In this case, the RDF will exhibit more features. To demonstrate this fact, we reconstruct the square lattice of points from above, but we now introduce some noise into the system.

```
[5]:
```

```
box, clean_points = freud.data.UnitCell.square().generate_system(10, scale=2, sigma_noise=0)
box, noisy_points = freud.data.UnitCell.square().generate_system(10, scale=2, sigma_noise=0.1)
aq_clean = freud.AABBQuery(box, clean_points)
aq_clean.plot(ax=plt.gca(), c='k', s=3)
aq_noisy = freud.AABBQuery(box, noisy_points)
deviations = np.linalg.norm(box.wrap(noisy_points-clean_points), axis=-1)
_, sc = aq_noisy.plot(ax=plt.gca(), c=deviations)
cbar = plt.colorbar(sc)
cbar.set_label('Distance from lattice site')
plt.show()
```

```
[6]:
```

```
plot_rdf(box, noisy_points, 'rdf')
```

In this RDF, we see the same rough features as we saw with the perfect lattice. However, the signal is much noisier, and in fact we see that increasing the number of bins essentially leads to overfitting of the data. As a result, we have to be careful with how we choose to bin our data when constructing the RDF object.

An alternative route for avoiding this problem can be using the cumulative RDF instead. The relationship between the cumulative RDF and the RDF is akin to that between a cumulative density and a probability density function, providing a measure of the total density of particles experienced up to some distance rather than the value at that distance. Just as a CDF can help avoid certain mistakes common to plotting a PDF, plotting the cumulative RDF may be helpful in some cases. Here, we see that
decreasing the bin size slightly alters the features of the plot, but only in very minor way (*i.e.* decreasing the smoothness of the line due to small jitters).

```
[7]:
```

```
plot_rdf(box, noisy_points, 'n_r')
```

##### freud.environment.AngularSeparation¶

The `freud.environment`

module analyzes the local environments of particles. The `freud.environment.AngularSeparation`

class enables direct measurement of the relative orientations of particles.

```
[1]:
```

```
import freud
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['axes.titlepad'] = 20
from mpl_toolkits.mplot3d import Axes3D
import rowan # for quaternion math, see rowan.readthedocs.io for more information.
```

In order to work with orientations in freud, we need to do some math with quaternions. If you are unfamiliar with quaternions, you can read more about their definition and how they can be used to represent rotations. For the purpose of this tutorial, just consider them as 4D vectors, and know that the set of normalized (*i.e.* unit norm) 4D vectors can be used to represent
rotations in 3D. In fact, there is a 1-1 mapping between normalized quaternions and 3x3 rotation matrices. Quaternions are more computationally convenient, however, because they only require storing 4 numbers rather than 9, and they can be much more easily chained together. The rowan library (rowan.readthedocs.io) defines many useful operations using quaternions, such as the rotations of vectors using quaternions instead of matrices.

###### Neighbor Angles¶

One usage of the AngularSeparation class is to compute angles between neighboring particles. To show how this works, we generate a simple configuration of particles with random orientations associated with each point.

```
[2]:
```

```
uc = freud.data.UnitCell.sc()
box, positions = uc.generate_system(5)
N = len(positions)
# Generate random, correlated particle orientations by taking identity
# quaternions and slightly rotating them in a random direction
np.random.seed(0)
interpolate_amount = 0.2
identity_quats = np.array([[1, 0, 0, 0]] * N)
ref_orientations = rowan.interpolate.slerp(
identity_quats, rowan.random.rand(N), interpolate_amount)
orientations = rowan.interpolate.slerp(
identity_quats, rowan.random.rand(N), interpolate_amount)
```

```
[3]:
```

```
# To show orientations, we use arrows rotated by the quaternions.
ref_arrowheads = rowan.rotate(ref_orientations, np.array([1, 0, 0]))
arrowheads = rowan.rotate(orientations, np.array([1, 0, 0]))
fig = plt.figure(figsize=(12, 6))
ref_ax = fig.add_subplot(121, projection='3d')
ax = fig.add_subplot(122, projection='3d')
ref_ax.quiver3D(positions[:, 0], positions[:, 1], positions[:, 2],
ref_arrowheads[:, 0], ref_arrowheads[:, 1], ref_arrowheads[:, 2])
ax.quiver3D(positions[:, 0], positions[:, 1], positions[:, 2],
arrowheads[:, 0], arrowheads[:, 1], arrowheads[:, 2])
ref_ax.set_title("Reference orientations", fontsize=16)
ax.set_title("Orientations", fontsize=16)
plt.show()
```

We can now use the AngularSeparation class to compare the orientations in these two systems.

```
[4]:
```

```
# For simplicity, we'll assume that our "particles" are completely
# asymmetric, i.e. there are no rotations that map the particle
# back onto itself. If we had a regular polyhedron, then we would
# want to specify all the quaternions that rotate that polyhedron
# onto itself.
equiv_orientations = np.array([[1, 0, 0, 0]])
ang_sep = freud.environment.AngularSeparationNeighbor()
ang_sep.compute(system=(box, positions),
orientations=orientations,
query_points=positions,
query_orientations=ref_orientations,
equiv_orientations=equiv_orientations,
neighbors={'num_neighbors': 12})
# Convert angles from radians to degrees and plot histogram
neighbor_angles = np.rad2deg(ang_sep.angles)
plt.hist(neighbor_angles)
plt.title('Histogram of angular separations between neighbors')
plt.xlabel('Angular separation (degrees)')
plt.ylabel('Frequency')
plt.show()
```

###### Global Angles¶

Alternatively, the AngularSeparationGlobal class can also be used to compute the orientation of all points in the system relative to some global set of orientations. In this case, we simply provide a set of global quaternions that we want to consider. For simplicity, let’s consider \(180^\circ\) rotations about each of the coordinate axes, which have very simple quaternion representations.

```
[5]:
```

```
global_orientations = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
ang_sep = freud.environment.AngularSeparationGlobal()
ang_sep.compute(global_orientations, ref_orientations, equiv_orientations)
global_angles = np.rad2deg(ang_sep.angles)
```

```
[6]:
```

```
plt.hist(global_angles[:, 0])
plt.title('Histogram of angular separation relative to identity quaternion')
plt.xlabel('Angular separation (degrees)')
plt.ylabel('Frequency')
plt.show()
```

As a simple check, we can ensure that for the identity quaternion \((1, 0, 0, 0)\), which performs a \(0^\circ\) rotation, the angles between the reference orientations and that quaternion are equal to the original angles of rotation of those quaternions (*i.e.* how much those orientations were already rotated relative to the identity).

```
[7]:
```

```
ref_axes, ref_angles = rowan.to_axis_angle(ref_orientations)
np.allclose(global_angles[:, 0], np.rad2deg(ref_angles), rtol=1e-4)
```

```
[7]:
```

```
True
```

##### freud.environment.BondOrder¶

###### Computing Bond Order Diagrams¶

The `freud.environment`

module analyzes the local environments of particles. In this example, the `freud.environment.BondOrder`

class is used to plot the bond order diagram (BOD) of a system of particles.

```
[1]:
```

```
import numpy as np
import freud
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
```

Our sample data will be taken from an face-centered cubic (FCC) structure. The array of points is rather large, so that the plots are smooth. Smaller systems may need to gather data from multiple frames in order to smooth the resulting array’s statistics, by computing multiple times with `reset=False`

.

```
[2]:
```

```
uc = freud.data.UnitCell.fcc()
box, points = uc.generate_system(40, sigma_noise=0.05)
```

Now we create a `BondOrder`

compute object and create some arrays useful for plotting.

```
[3]:
```

```
n_bins_theta = 100
n_bins_phi = 100
bod = freud.environment.BondOrder((n_bins_theta, n_bins_phi))
phi = np.linspace(0, np.pi, n_bins_phi)
theta = np.linspace(0, 2*np.pi, n_bins_theta)
phi, theta = np.meshgrid(phi, theta)
x = np.sin(phi) * np.cos(theta)
y = np.sin(phi) * np.sin(theta)
z = np.cos(phi)
```

Next, we use the `compute`

method and the `bond_order`

property to return the array.

```
[4]:
```

```
bod_array = bod.compute(system=(box, points), neighbors={'num_neighbors': 12}).bond_order
# Clean up polar bins for plotting
bod_array = np.clip(bod_array, 0, np.percentile(bod_array, 99))
plt.imshow(bod_array.T)
plt.show()
```

This code shows the bond order diagram on a sphere as the sphere is rotated. The code takes a few seconds to run, so be patient.

```
[5]:
```

```
fig = plt.figure(figsize=(12, 8))
for plot_num in range(6):
ax = fig.add_subplot(231 + plot_num, projection='3d')
ax.plot_surface(x, y, z, rstride=1, cstride=1, shade=False,
facecolors=matplotlib.cm.viridis(bod_array / np.max(bod_array)))
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)
ax.set_axis_off()
# View angles in degrees
view_angle = 0, plot_num*15
ax.view_init(*view_angle)
plt.show()
```

We can also use a custom neighbor query to determine bonds. For example, we can filter for a range of bond lengths. Below, we only consider neighbors between \(r_{min} = 2.5\) and \(r_{max} = 3\) and plot the resulting bond order diagram.

```
[6]:
```

```
bod_array = bod.compute(system=(box, points), neighbors={'r_max': 3.0, 'r_min': 2.5}).bond_order
# Clean up polar bins for plotting
bod_array = np.clip(bod_array, 0, np.percentile(bod_array, 99))
plt.imshow(bod_array.T)
plt.show()
```

##### freud.environment.EnvironmentCluster¶

The `freud.environment.EnvironmentCluster`

class finds and clusters local environments, as determined by the vectors pointing to neighbor particles. Neighbors can be defined by a cutoff distance or a number of nearest-neighbors, and the resulting `freud.locality.NeighborList`

is used to enumerate a set of vectors, defining an “environment.” These environments are compared with the environments of neighboring particles to form spatial clusters, which usually correspond to grains, droplets, or
crystalline domains of a system. `EnvironmentCluster`

has several parameters that alter its behavior, please see the documentation or helper functions below for descriptions of these parameters.

In this example, we cluster the local environments of hexagons. Clusters with 5 or fewer particles are colored dark gray.

```
[1]:
```

```
import numpy as np
import freud
from collections import Counter
import matplotlib.pyplot as plt
def get_cluster_arr(system, num_neighbors, threshold,
registration=False, global_search=False):
"""Computes clusters of particles' local environments.
Args:
system:
Any object that is a valid argument to
:class:`freud.locality.NeighborQuery.from_system`.
num_neighbors (int):
Number of neighbors to consider in every particle's local environment.
threshold (float):
Maximum magnitude of the vector difference between two vectors,
below which we call them matching.
global_search (bool):
If True, do an exhaustive search wherein the environments of
every single pair of particles in the simulation are compared.
If False, only compare the environments of neighboring particles.
registration (bool):
Controls whether we first use brute force registration to
orient the second set of vectors such that it minimizes the
RMSD between the two sets.
Returns:
tuple(np.ndarray, dict): array of cluster indices for every particle
and a dictionary mapping from cluster_index keys to vector_array)
pairs giving all vectors associated with each environment.
"""
# Perform the env-matching calcuation
neighbors = {'num_neighbors': num_neighbors}
match = freud.environment.EnvironmentCluster()
match.compute(system, threshold, neighbors=neighbors,
registration=registration, global_search=global_search)
return match.cluster_idx, match.cluster_environments
def color_by_clust(cluster_index_arr, no_color_thresh=1,
no_color='#333333', cmap=plt.get_cmap('viridis')):
"""Takes a cluster_index_array for every particle and returns a
dictionary of (cluster index, hexcolor) color pairs.
Args:
cluster_index_arr (numpy.ndarray):
The array of cluster indices, one per particle.
no_color_thresh (int):
Clusters with this number of particles or fewer will be
colored with no_color.
no_color (color):
What we color particles whose cluster size is below no_color_thresh.
cmap (color map):
The color map we use to color all particles whose
cluster size is above no_color_thresh.
"""
# Count to find most common clusters
cluster_counts = Counter(cluster_index_arr)
# Re-label the cluster indices by size
color_count = 0
color_dict = {cluster[0]: counter for cluster, counter in
zip(cluster_counts.most_common(),
range(len(cluster_counts)))}
# Don't show colors for clusters below the threshold
for cluster_id in cluster_counts:
if cluster_counts[cluster_id] <= no_color_thresh:
color_dict[cluster_id] = -1
OP_arr = np.linspace(0.0, 1.0, max(color_dict.values())+1)
# Get hex colors for all clusters of size greater than no_color_thresh
for old_cluster_index, new_cluster_index in color_dict.items():
if new_cluster_index == -1:
color_dict[old_cluster_index] = no_color
else:
color_dict[old_cluster_index] = cmap(OP_arr[new_cluster_index])
return color_dict
```

We load the simulation data and call the analysis functions defined above. Notice that we use 6 nearest neighbors, since our system is made of hexagons that tend to cluster with 6 neighbors.

```
[2]:
```

```
ex_data = np.load('data/MatchEnv_Hexagons.npz')
box = ex_data['box']
positions = ex_data['positions']
orientations = ex_data['orientations']
aq = freud.AABBQuery(box, positions)
cluster_index_arr, cluster_envs = get_cluster_arr(
aq, num_neighbors=6, threshold=0.2,
registration=False, global_search=False)
color_dict = color_by_clust(cluster_index_arr, no_color_thresh=5)
colors = [color_dict[i] for i in cluster_index_arr]
```

Below, we plot the resulting clusters. The colors correspond to the cluster size.

```
[3]:
```

```
plt.figure(figsize=(12, 12), facecolor='white')
aq.plot(ax=plt.gca(), c=colors, s=20)
plt.title('Clusters Colored by Particle Local Environment')
plt.show()
```

##### freud.environment.LocalDescriptors: Steinhardt Order Parameters from Scratch¶

The `freud.environment`

module analyzes the local environments of particles. The `freud.environment.LocalDescriptors`

class is a useful tool for analyzing identifying crystal structures in a rotationally invariant manner using local particle environments. The primary purpose of this class is to compute spherical harmonics between neighboring particles in a way that orients particles correctly relative to their local environment, ensuring that global orientational shifts do not change the
output.

```
[1]:
```

```
import freud
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
```

###### Computing Spherical Harmonics¶

To demonstrate the basic application of the class, let’s compute the spherical harmonics between neighboring particles. For simplicity, we consider points on a simple cubic lattice.

```
[2]:
```

```
uc = freud.data.UnitCell.sc()
box, points = uc.generate_system(5)
system = freud.AABBQuery(box, points)
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
system.plot(ax=ax)
ax.set_title("Simple cubic crystal", fontsize=16)
plt.show()
```

Now, let’s use the class to compute an array of spherical harmonics for the system. The harmonics are computed for each bond, where a bond is defined by a pair of particles that are determined to lie within each others’ nearest neighbor shells based on a standard neighbor list search. The number of bonds and spherical harmonics to calculate is configurable.

```
[3]:
```

```
num_neighbors = 6
l_max = 12
nlist = system.query(points, {'num_neighbors': num_neighbors, 'exclude_ii': True}).toNeighborList()
ld = freud.environment.LocalDescriptors(l_max, mode='global')
ld.compute(system, neighbors=nlist);
```

###### Accessing the Data¶

The resulting spherical harmonic array has a shape corresponding to the number of neighbors. We can now extract the spherical harmonics corresponding to a particular \((l, m)\) pair using the ordering used by the `LocalDescriptors`

class: increasing values of \(l\), and for each \(l\), the nonnegative \(m\) values followed by the negative values.

```
[4]:
```

```
sph_raw = np.mean(ld.sph, axis=0)
count = 0
sph = np.zeros((l_max+1, l_max+1), dtype=np.complex128)
for l in range(l_max+1):
for m in range(l+1):
sph[l, m] = sph_raw[count]
count += 1
for m in range(-l, 0):
sph[l, m] = sph_raw[count]
count += 1
```

###### Using Spherical Harmonics to Compute Steinhardt Order Parameters¶

The raw per bond spherical harmonics are not typically useful quantities on their own. However, they can be used to perform sophisticated crystal structure analyses with different methods; for example, the pythia library uses machine learning to find patterns in the spherical harmonics computed by this class. In this notebook, we’ll use the quantities for a more classical application: the computation of Steinhardt order parameters. The order
parameters \(q_l\) provide a rotationally invariant measure of the system that can for some structures, provide a unique identifying fingerprint. They are a particularly useful measure for various simple cubic structures such as structures with underlying simple cubic, BCC, or FCC lattices. The `freud`

library actually provides additional classes to efficiently calculate these order parameters directly, but they also provide a reasonable demonstration here.

For more information on Steinhardt order parameters, see the original paper or the `freud.order.Steinhardt`

documentation.

```
[5]:
```

```
def get_ql(num_particles, descriptors, nlist, weighted=False):
"""Given a set of points and a LocalDescriptors object (and the
underlying NeighborList), compute the per-particle Steinhardt ql
order parameter for all :math:`l` values up to the maximum quantum
number used in the computation of the descriptors."""
qbar_lm = np.zeros((num_particles, descriptors.sph.shape[1]),
dtype=np.complex128)
for i in range(num_particles):
indices = nlist.query_point_indices == i
Ylms = descriptors.sph[indices, :]
if weighted:
weights = nlist.weights[indices, np.newaxis]
weights /= np.sum(weights)
num_neighbors = 1
else:
weights = np.ones_like(Ylms)
num_neighbors = descriptors.sph.shape[0]/num_particles
qbar_lm[i, :] = np.sum(Ylms * weights, axis=0)/num_neighbors
ql = np.zeros((qbar_lm.shape[0], descriptors.l_max+1))
for i in range(ql.shape[0]):
for l in range(ql.shape[1]):
for k in range(l**2, (l+1)**2):
ql[i, l] += np.absolute(qbar_lm[i, k])**2
ql[i, l] = np.sqrt(4*np.pi/(2*l + 1) * ql[i, l])
return ql
ld_ql = get_ql(len(points), ld, nlist)
```

Since `freud`

provides the ability to calculate these parameter as well, we can directly check that our answers are correct. *Note: More information on the ``Steinhardt`` class can be found in the documentation or in the ``Steinhardt`` example.*

```
[6]:
```

```
L = 6
steinhardt = freud.order.Steinhardt(l=L)
steinhardt.compute(system, neighbors=nlist)
if np.allclose(steinhardt.ql, ld_ql[:, L]):
print("Our manual calculation matches the Steinhardt class!")
```

```
Our manual calculation matches the Steinhardt class!
```

For a brief demonstration of why the Steinhardt order parameters can be useful, let’s look at the result of thermalizing our points and recomputing this measure.

```
[7]:
```

```
sigmas = [0.03, 0.05, 0.1]
systems = []
nlists = []
for sigma in sigmas:
box, points = uc.generate_system(5, sigma_noise=sigma)
system = freud.AABBQuery(box, points)
systems.append(system)
nlists.append(
system.query(
points, {'num_neighbors': num_neighbors, 'exclude_ii': True}
).toNeighborList()
)
```

```
[8]:
```

```
fig = plt.figure(figsize=(14, 6))
axes = []
for i, v in enumerate(sigmas):
ax = fig.add_subplot("1{}{}".format(len(sigmas), i+1), projection='3d')
systems[i].plot(ax=ax)
ax.set_title("$\sigma$ = {}".format(v), fontsize=16);
plt.show()
```

If we recompute the Steinhardt OP for each of these data sets, we see that adding noise has the effect of smoothing the order parameter such that the peak we observed for the perfect crystal is no longer observable.

```
[9]:
```

```
ld_qls = []
for i, sigma in enumerate(sigmas):
ld = freud.environment.LocalDescriptors(l_max, mode='global')
ld.compute(systems[i], neighbors=nlists[i])
ld_qls.append(get_ql(len(systems[i].points), ld, nlists[i]))
```

```
[10]:
```

```
fig, ax = plt.subplots()
for i, ld_ql in enumerate(ld_qls):
lim_out = ax.hist(ld_ql[:, L], label="$\sigma$ = {}".format(sigmas[i]), density=True)
if i == 0:
# Can choose any element, all are identical in the reference case
ax.vlines(ld_ql[:, L][0], 0, np.max(lim_out[0]), label='Reference')
ax.set_title("Histogram of $q_{L}$ values".format(L=L), fontsize=16)
ax.set_ylabel("Frequency", fontsize=14)
ax.set_xlabel("$q_{L}$".format(L=L), fontsize=14)
ax.legend(fontsize=14)
plt.show()
```

This type of identification process is what the LocalDescriptors data outputs may be used for. In the case of Steinhardt OPs, it provides a simple fingerprint for comparing thermalized systems to a known ideal structure to measure their similarity.

For reference, we can also check these values against the `Steinhardt`

class again.

```
[11]:
```

```
for i, (system, nlist) in enumerate(zip(systems, nlists)):
steinhardt = freud.order.Steinhardt(l=L)
steinhardt.compute(system, nlist)
if np.allclose(steinhardt.particle_order, ld_qls[i][:, L]):
print("Our manual calculation matches the Steinhardt class!")
```

```
Our manual calculation matches the Steinhardt class!
Our manual calculation matches the Steinhardt class!
Our manual calculation matches the Steinhardt class!
```

##### freud.interface.Interface¶

###### Locating Particles on Interfacial Boundaries¶

The `freud.interface`

module compares the distances between two sets of points to determine the interfacial particles.

```
[1]:
```

```
import freud
import numpy as np
import matplotlib.pyplot as plt
```

To make a pretend data set, we create a large number of **blue (-1)** particles on a square grid. Then we place grain centers on a larger grid and draw grain radii from a normal distribution. We color the particles **yellow (+1)** if their distance from a grain center is less than the grain radius.

```
[2]:
```

```
np.random.seed(0)
system_size = 100
num_grains = 4
uc = freud.data.UnitCell.square()
box, points = uc.generate_system(num_replicas=system_size, scale=1)
_, centroids = uc.generate_system(num_replicas=num_grains, scale=system_size/num_grains)
system = freud.AABBQuery(box, points)
values = np.array([-1 for p in points])
grain_radii = np.abs(np.random.normal(size=num_grains**2, loc=5, scale=2))
for center, radius in zip(centroids, grain_radii):
for i, j, dist in system.query(center, {'r_max': radius}):
values[j] = 1
plt.figure(figsize=(9, 9))
system.plot(ax=plt.gca(), c=values, cmap='cividis', s=12)
plt.title('System of two particle types')
plt.show()
```

This system is **phase-separated** because the yellow particles are generally near one another, and so are the blue particles.

We can use `freud.interface.InterfaceMeasure`

to label the particles on either side of the yellow-blue boundary. The class can tell us how many points are on either side of the interface:

```
[3]:
```

```
iface = freud.interface.Interface()
iface.compute((box, points[values > 0]), points[values < 0], neighbors={'r_max': 1.5})
print('There are {} query points (blue) on the interface.'.format(iface.query_point_count))
print('There are {} points (yellow) on the interface.'.format(iface.point_count))
```

```
There are 856 query points (blue) on the interface.
There are 724 points (yellow) on the interface.
```

Now we can plot the particles on the interface. We color the outside of the interface dark blue and the inside of the interface yellow.

```
[4]:
```

```
plt.figure(figsize=(9, 9))
interface_values = np.zeros(len(points))
interface_values[np.where(values < 0)[0][iface.query_point_ids]] = -1
interface_values[np.where(values > 0)[0][iface.point_ids]] = 1
system.plot(ax=plt.gca(), c=interface_values, cmap='cividis', s=12)
plt.title('Particles on the interface between types')
plt.show()
```

##### freud.order.Hexatic: Hard Hexagons¶

###### Hexatic Order Parameter¶

The hexatic order parameter measures how closely the local environment around a particle resembles perfect \(k\)-atic symmetry, *e.g.* how closely the environment resembles hexagonal/hexatic symmetry for \(k=6\). The order parameter is given by:

where \(\theta_{ij}\) is the angle between the vector \(\vec{r}_{ij}\) and \(\left(1, 0\right)\).

The pseudocode is given below:

```
for each particle i:
neighbors = nearestNeighbors(i, n):
for each particle j in neighbors:
r_ij = position[j] - position[i]
theta_ij = arctan2(r_ij.y, r_ij.x)
psi_array[i] += exp(complex(0,k*theta_ij))
```

The data sets used in this example are a system of hard hexagons, simulated in the NVT thermodynamic ensemble in HOOMD-blue, for a dense fluid of hexagons at packing fraction \(\phi = 0.65\) and solids at packing fractions \(\phi = 0.75, 0.85\).

```
[1]:
```

```
import numpy as np
import freud
from bokeh.io import output_notebook
from bokeh.plotting import figure, show
# The following line imports this set of utility functions:
# https://github.com/glotzerlab/freud-examples/blob/master/util.py
import util
output_notebook()
```

```
[2]:
```

```
def plot_hex_order_param(data_path, title):
# Create hexatic object
hex_order = freud.order.Hexatic(k=6)
# Load the data
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
quat_data = np.load("{}/quat_data.npy".format(data_path))
# Grab data from last frame
box = box_data[-1].tolist()
points = pos_data[-1]
quats = quat_data[-1]
angles = 2*np.arctan2(quats[:, 3], quats[:, 0])
# Compute hexatic order for 6 nearest neighbors
hex_order.compute(system=(box, points), neighbors={'num_neighbors': 6})
psi_k = hex_order.particle_order
avg_psi_k = np.mean(psi_k)
# Create hexagon vertices
verts = util.make_polygon(sides=6, radius=0.6204)
# Create array of transformed positions
patches = util.local_to_global(verts, points[:, :2], angles)
# Create an array of angles relative to the average
relative_angles = np.angle(psi_k) - np.angle(avg_psi_k)
# Plot in bokeh
p = figure(title=title)
p.patches(xs=patches[:, :, 0].tolist(), ys=patches[:, :, 1].tolist(),
fill_color=[util.cubeellipse(x) for x in relative_angles],
line_color="black")
util.default_bokeh(p)
show(p)
```

```
[3]:
```

```
plot_hex_order_param('data/phi065', 'Hexatic Order Parameter, 0.65 density')
```