freud documentation¶
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, and clusters, as well as original analysis methods including potentials of mean force and torque (PMFTs) and local environment matching. The freud library uses single-precision NumPy arrays for input and output, enabling integration with the scientific Python ecosystem for many typical materials science workflows.
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:
Python (2.7+ required, 3.5+ recommended)
The following are optional for installing freud:
Cython (0.28+ required): The freud repository contains Cython-generated
*.cpp
files in thefreud/
directory that can be used directly. However, Cython is necessary if you wish to recompile these files.
For conda users, 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
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). Additionally, due to limitations inherent to the distutils/setuptools infrastructure, building extension modules can only be parallelized using the build_ext subcommand of setup.py, not with install. As a result, it will be faster to manually run build_ext and then install (which normally calls build_ext under the hood anyway) the built packages. In general, the following options are available for setup.py in addition to the standard setuptools options (notes are included to indicate which options are only available for specific subcommands such as build_ext):
- --PRINT-WARNINGS
Specify whether or not to print compilation warnings resulting from the build even if the build succeeds with no errors.
- --ENABLE-CYTHON
Rebuild the Cython-generated C++ files. If there are any unexpected issues with compiling the C++ shipped with the build, using this flag may help. It is also necessary any time modifications are made to the Cython files.
- -j
Compile in parallel. This affects both the generation of C++ files from Cython files and the subsequent compilation of the source files. In the latter case, this option controls the number of Python modules that will be compiled in parallel.
- --TBB-ROOT
The root directory where TBB is installed. Useful if TBB is installed in a non-standard location or cannot be located by Python for some other reason. Note that this information can also be provided using the environment variable TBB_ROOT. The options –TBB-INCLUDE and –TBB-LINK will take precedence over –TBB-ROOT if both are specified.
- --TBB-INCLUDE
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 by Python for some other reason. Note that this information can also be provided using the environment variable TBB_ROOT. The options –TBB-INCLUDE and –TBB-LINK will take precedence over –TBB-ROOT if both are specified.
- --TBB-LINK
The directory where the TBB shared library (e.g. libtbb.so or libtbb.dylib) is located. Useful if TBB is installed in a non-standard location or cannot be located by Python for some other reason. Note that this information can also be provided using the environment variable TBB_ROOT. The options –TBB-INCLUDE and –TBB-LINK will take precedence over –TBB-ROOT if both are specified.
The following additional arguments are primarily useful for developers:
- --COVERAGE
Build the Cython files with coveragerc support to check unit test coverage.
- --NTHREAD
Specify the number of threads to allocate to compiling each module. This option is primarily useful for rapid development, particularly when all changes are in one module. While the -j option will not help parallelize this case, this option allows compilation of multiple source files belonging to the same module in parallel.
Note
freud makes use of submodules. If you ever wish to manually update these, you can execute:
git submodule update --init
Unit Tests¶
The unit tests for freud are included in the repository and are configured to be run using the Python unittest
library:
# Run tests from the tests directory
cd tests
python -m unittest discover .
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, but you may also build the documentation yourself:
Building the documentation¶
The following are required for building freud documentation:
You can install sphinx using conda
conda install sphinx
or from PyPi
pip install sphinx
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
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 [Matplotlib] and Bokeh [Bokeh], unless otherwise noted.
Key concepts¶
There are a few critical concepts, algorithms, and data structures that are central to all of freud. The box module defines the concept of a periodic simulation box, and the locality module defines methods for finding nearest neighbors for particles. Since both of these are used throughout freud, we recommend familiarizing yourself with these first, before delving into the workings of specific freud analysis modules.
Box¶
The goal of freud
is to perform generic analyses of particle simulations. Such simulations are always conducted within some region representing physical space; in freud
, these regions are known as simulation boxes, or simply boxes. An important characteristic of many simulations is that the simulation box is periodic, i.e. particles can travel and interact across system boundaries (for more information, see the Wikipedia
page). Simulations frequently use periodic boundary conditions to effectively simulate infinite systems without actually having to include an infinite number of particles. In such systems, a box in N dimensions can be represented by N linearly independent vectors.
The Box
class provides the standard API for such simulation boxes throughout freud
. The class represents some 2- or 3-dimensional region of space, and it provides utility functions for interacting with this space, including the ability to wrap vectors outside this box into the box according to periodic boundary conditions. Boxes are represented according to the HOOMD-blue convention for boxes. According to this convention, a 3D
(2D) simulation box is fully defined by 3 (2) linearly independent vectors, which are represented by 3 (2) characteristic lengths and 3 (1) tilt factors indicating how these vectors are angled with respect to one another. With this convention, a generic box is represented by the following \(3\times3\) matrix:
where \(xy\), \(xz\), and \(yz\) are the tilt factors. Note that this convention imposes the requirement that the box vectors form a right-handed coordinate system, which manifests itself in the form of an upper (rather than lower) triangular box matrix.
In this notebook, we demonstrate the basic features of the Box
class, particularly the facility for wrapping particles back into the box under periodic boundary conditions. For more information, see the freud.box
documentation.
Box Creation¶
There are many ways to construct a box. We demonstrate all of these below, with some discussion of when they might be useful.
Default (full) API¶
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 a 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 a namedtuple with the appropriate entries. Any other object that provides a similar API for attribute-based access of \(L_x\), \(L_y\), \(L_z\), \(xy\), \(xz\), and \(yz\) (or some subset) will work equally well. 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', 'dimensions'])
box = freud.box.Box.from_box(MyBox(Lx=5, Ly=3, Lz=2, xy=0, xz=0, yz=0, dimensions=3))
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, dimensions=2))
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 namedtuple representations, which provide completely specified descriptions of the box. Note that the namedtuple type used by freud boxes, the BoxTuple, is simply an internal representation.
[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_tuple()
[9]:
BoxTuple(Lx=5.0, Ly=5.0, Lz=5.0, xy=0.0, xz=0.0, yz=0.0)
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.Linv]))
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 to and from fractional 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]:
# Conversion to coordinate representation from fractions.
print(box.makeCoordinates([0, 0, 0]))
print(box.makeCoordinates([0.5, 0.5, 0.5]))
print(box.makeCoordinates([0.8, 0.3, 1]))
print()
# Conversion to and from coordinate representation, resulting
# in the input fractions.
print(box.makeFraction(box.makeCoordinates([0, 0, 0])))
print(box.makeFraction(box.makeCoordinates([0.5, 0.5, 0.5])))
print("[{:1.1f}, {:1.1f}, {:1.1f}]".format(*box.makeFraction(box.makeCoordinates([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.0]
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]:
# We define box plot generation separately
from util import box_2d_to_points
# 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]])
points = box_2d_to_points(box)
[13]:
from matplotlib import pyplot as plt
fig, ax = plt.subplots(figsize=(9, 6))
ax.plot(points[:, 0], points[:, 1], color='k')
plt.show()
<Figure size 900x600 with 1 Axes>
[14]:
plt.figure()
plt.plot(points[:, 0], points[:, 1])
plt.show()

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.
[15]:
fig, ax = plt.subplots(figsize=(9, 6))
ax.plot(points[:, 0], points[:, 1], color='k')
ax.plot(points[:, 0] + Lx, points[:, 1], linestyle='dashed', color='k')
ax.plot(points[:, 0] - Lx, points[:, 1], linestyle='dashed', color='k')
ax.plot(points[:, 0] + xy*Ly, points[:, 1] + Ly, linestyle='dashed', color='k')
ax.plot(points[:, 0] - xy*Ly, points[:, 1] - Ly, linestyle='dashed', color='k')
plt.show()

Any particles in the original box will also therefore be seen as existing in all the neighboring boxes.
[16]:
np.random.seed(0)
tmp = np.random.rand(5, 2)
origin = np.array(box.makeCoordinates([0, 0, 0]))
u = np.array(box.makeCoordinates([1, 0, 0])) - origin
v = np.array(box.makeCoordinates([0, 1, 0])) - origin
particles = u*tmp[:, [0]] + v*tmp[:, [1]]
[17]:
fig, ax = plt.subplots(figsize=(9, 6))
# Plot the boxes.
ax.plot(points[:, 0], points[:, 1], color='k')
ax.plot(points[:, 0] + Lx, points[:, 1], linestyle='dashed', color='k')
ax.plot(points[:, 0] - Lx, points[:, 1], linestyle='dashed', color='k')
ax.plot(points[:, 0] + xy*Ly, points[:, 1] + Ly, linestyle='dashed', color='k')
ax.plot(points[:, 0] - xy*Ly, points[:, 1] - Ly, linestyle='dashed', color='k')
# Plot the points in the original box.
ax.plot(particles[:, 0] + origin[0], particles[:, 1] + origin[1],
linestyle='None', marker='.', color='#1f77b4')
# Define the different origins.
origins = []
origins.append(np.array(box.makeCoordinates([-1, 0, 0])))
origins.append(np.array(box.makeCoordinates([1, 0, 0])))
origins.append(np.array(box.makeCoordinates([0, -1, 0])))
origins.append(np.array(box.makeCoordinates([0, 1, 0])))
# Plot particles in each of the periodic boxes.
for o in origins:
ax.plot(particles[:, 0] + o[0], particles[:, 1] + o[1],
linestyle='None', marker='.', color='#1f77b4')
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.
[18]:
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# Plot the boxes.
for i, ax in enumerate(axes.flatten()):
ax.plot(points[:, 0], points[:, 1], color='k')
ax.plot(points[:, 0] + Lx, points[:, 1], linestyle='dashed', color='k')
ax.plot(points[:, 0] - Lx, points[:, 1], linestyle='dashed', color='k')
ax.plot(points[:, 0] + xy*Ly, points[:, 1] + Ly, linestyle='dashed', color='k')
ax.plot(points[:, 0] - xy*Ly, points[:, 1] - Ly, linestyle='dashed', color='k')
# Plot the points relative to origin i.
o = origins[i]
ax.plot(particles[:, 0] + o[0], particles[:, 1] + o[1],
linestyle='None', marker='.', label='Original')
# Now wrap these points and plot them.
wrapped_particles = box.wrap(particles + o)
ax.plot(wrapped_particles[:, 0], wrapped_particles[:, 1],
linestyle='None', marker='.', label='Wrapped')
ax.tick_params(axis="both", which="both", labelsize=14)
ax.legend(fontsize=14)
plt.show()

ParticleBuffer - Unit Cell RDF¶
The ParticleBuffer
class is meant to replicate particles 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]:
import freud
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from util import box_2d_to_points
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.asarray([[-0.5, -0.5, -0.5], [0.5, 0.5, 0.5]])
corners = box_2d_to_points(box)
ax = plt.gca()
box_patch = plt.Polygon(corners[:, :2])
patch_collection = matplotlib.collections.PatchCollection([box_patch], edgecolors='black', alpha=0.4)
ax.add_collection(patch_collection)
plt.scatter(points[:, 1], points[:, 2])
plt.show()

Next, we create a ParticleBuffer
instance and have it compute the “buffer” particles that lie outside the first periodicity. These positions are stored in the buffer_positions
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.box.ParticleBuffer(box)
pbuff.compute(points, 6, images=True)
print(pbuff.buffer_particles[:10], '...')
[[ 0.6547002 1.5 0. ]
[ 1.8094003 3.5 0. ]
[ 2.9641018 5.5 0. ]
[-3.9641013 -6.5 0. ]
[-2.809401 -4.4999995 0. ]
[-1.6547002 -2.5000005 0. ]
[ 1.5000002 -0.5 0. ]
[ 2.6547008 1.5 0. ]
[ 3.8094003 3.5 0. ]
[ 4.964102 5.5 0. ]] ...
Below, we plot the original unit cell and the replicated buffer points and buffer box.
[4]:
plt.scatter(points[:, 0], points[:, 1])
plt.scatter(pbuff.buffer_particles[:, 0], pbuff.buffer_particles[:, 1])
box_patch = plt.Polygon(corners[:, :2])
buff_corners = box_2d_to_points(pbuff.buffer_box)
buff_box_patch = plt.Polygon(buff_corners[:, :2])
patch_collection = matplotlib.collections.PatchCollection(
[box_patch, buff_box_patch], facecolors=['blue', 'orange'],
edgecolors='black', alpha=0.2)
plt.gca().add_collection(patch_collection)
plt.show()

Finally, we can plot the radial distribution function (RDF) of this replicated system, using a value of rmax
that is larger than the size of the original box. This allows us to see the interaction of the original particles in ref_points
with their replicated neighbors from the buffer in points
.
[5]:
rdf = freud.density.RDF(rmax=5, dr=0.02)
rdf.compute(pbuff.buffer_box, ref_points=points, points=pbuff.buffer_particles)
plt.plot(rdf.R, rdf.RDF)
plt.show()

LinkCell¶
Many of the most powerful analyses of particle simulations involve some characterization of the local environments of particles. Whether the analyses involve finding clusters, identifying interfaces, computing order parameters, or something else entirely, they always require finding particles in proximity to others so that properties of the local environment can be computed. The freud.locality.NeighborList
and freud.locality.LinkCell
classes are the fundamental building blocks for this
type of calculation. The NeighborList
class is essentially a container for particle pairs that are determined to be adjacent to one another. The LinkCell
class implements the standard linked-list cell algorithm, in which a cell list is computed using linked lists to store the particles in each cell. In this notebook, we provide a brief demonstration of how this data structure works and how it
is used throughout freud.
We begin by demonstrating how a cell list works, which is essentially by dividing space into fixed width cells.
[1]:
from __future__ import division
import freud
import numpy as np
from matplotlib import pyplot as plt
import timeit
# place particles 0 and 1 in cell 0
# place particle 2 in cell 1
# place particles 3,4,5 in cell 3
# and no particles in cells 4,5,6,7
particles = np.array([[-0.5, -0.5, 0],
[-0.6, -0.6, 0],
[0.5, -0.5, 0],
[-0.5, 0.5, 0],
[-0.6, 0.6, 0],
[-0.7, 0.7, 0]], dtype='float32')
L = 2 # The box size
r_max = 1 # The cell width, and the nearest neighbor distance
box = freud.box.Box.square(L)
lc = freud.locality.LinkCell(box, r_max)
lc.compute(box, particles)
for c in range(0, lc.num_cells):
print("The following particles are in cell {}: {}".format(c, ', '.join([str(x) for x in lc.itercell(c)])))
The following particles are in cell 0: 0, 1
The following particles are in cell 1: 2
The following particles are in cell 2: 3, 4, 5
The following particles are in cell 3:
[2]:
from matplotlib import patches
from matplotlib import cm
cmap = cm.get_cmap('plasma')
colors = [cmap(i/lc.num_cells) for i in range(lc.num_cells)]
fig, ax = plt.subplots(figsize=(9, 6))
ax.scatter(particles[:, 0], particles[:, 1])
ax.set_xlim([-1, 1])
ax.set_ylim([-1, 1])
corners = [(-1, -1), (0, -1), (-1, 0), (0, 0)]
handles = []
labels = []
for i, corner in enumerate(corners):
p = patches.Rectangle(corner, 1, 1, color=colors[i], alpha=0.3)
ax.add_patch(p)
handles.append(p)
labels.append("Cell {}".format(i))
ax.tick_params(axis='both', which='both', labelsize=14)
fig.legend(handles, labels, fontsize=16)
fig.subplots_adjust(right=0.8)
for i, p in enumerate(particles):
ax.text(p[0]+0.05, p[1]-0.03, "Particle {}".format(i), fontsize=14)

The principle behind a cell list is that depending on how close particles have to be to be considered neighbors, we can construct a cell list of an appropriate width such that a given particle’s neighbors can always be found by only looking in the neighboring cells, saving us the work of checking all the other particles in the system. We can now extract the NeighborList object computed using this cell list for finding particle neighbors.
[3]:
nlist = lc.nlist
for i in set(nlist.index_i):
js = nlist.index_j[nlist.index_i == i]
print("The particles within a distance 1 of particle {} are: {}".format(
i, ', '.join([str(j) for j in js])))
The particles within a distance 1 of particle 0 are: 1, 4, 5
The particles within a distance 1 of particle 1 are: 0, 2, 3, 4, 5
The particles within a distance 1 of particle 2 are: 1
The particles within a distance 1 of particle 3 are: 1, 4, 5
The particles within a distance 1 of particle 4 are: 0, 1, 3, 5
The particles within a distance 1 of particle 5 are: 0, 1, 3, 4
Finally, we can easily check this computation manually by just computing particle distances. Note that we need to be careful to make sure that we properly respect the box periodicity, which means that interparticle distances should be calculated according to the minimum image convention. In essence, this means that since the box is treated as being infinitely replicated in all directions, we have to ensure that each particle is only interacting with the closest copy of another particle. We can easily enforce this here by making sure that particle distances are never large than half the box length in any given dimension.
[4]:
def compute_distances(box, positions):
"""Compute pairwise particle distances, taking into account PBCs.
Args:
box (:class:`freud.box.Box`): The simulation box the particles live in.
positions (:class:`np.ndarray`): The particle positions.
"""
# First we shift all the particles so that the coordinates lie from
# [0, L] rather than [-L/2, L/2].
positions[:, 0] = np.mod(positions[:, 0]+box.Lx/2, box.Lx)
positions[:, 1] = np.mod(positions[:, 1]+box.Ly/2, box.Ly)
positions[:, 0] = np.mod(positions[:, 0]+box.Lx/2, box.Lx)
positions[:, 1] = np.mod(positions[:, 1]+box.Ly/2, box.Ly)
# To apply minimum image convention, we check if the distance is
# greater than half the box length in either direction, and if it
# is, we replace it with L-distance instead. We use broadcasting
# to get all pairwise positions, then modify the pos2 array where
# the distance is found to be too large for a specific pair.
pos1, pos2 = np.broadcast_arrays(positions[np.newaxis, :, :], positions[:, np.newaxis, :])
vectors = pos1 - pos2
pos2[:, :, 0] = np.where(np.abs(vectors[:, :, 0]) > box.Lx/2,
box.Lx - np.abs(pos2[:, :, 0]),
pos2[:, :, 0])
pos2[:, :, 1] = np.where(np.abs(vectors[:, :, 1]) > box.Ly/2,
box.Ly - np.abs(pos2[:, :, 1]),
pos2[:, :, 1])
distances = np.linalg.norm(pos1 - pos2, axis=-1)
return distances
[5]:
pairwise_distances = compute_distances(box, particles)
for i in range(pairwise_distances.shape[0]):
js = np.where(pairwise_distances[i, :] < r_max)
print("The particles within a distance 1 of particle {} are: {}".format(
i, ', '.join([str(j) for j in js[0] if not j==i])))
The particles within a distance 1 of particle 0 are: 1, 4, 5
The particles within a distance 1 of particle 1 are: 0, 2, 3, 4, 5
The particles within a distance 1 of particle 2 are: 1
The particles within a distance 1 of particle 3 are: 1, 4, 5
The particles within a distance 1 of particle 4 are: 0, 1, 3, 5
The particles within a distance 1 of particle 5 are: 0, 1, 3, 4
For larger systems, however, such pairwise calculations would quickly become prohibitively expensive. The primary benefit of the LinkCell object is that it can dramatically improve this cost.
[6]:
log_Ns = np.arange(5, 12)
lc_times = []
naive_times = []
for log_N in log_Ns:
print("Running for log_N = {}".format(log_N))
particles = np.random.rand(int(2**log_N), 3)*L-L/2
particles[:, 0] = 0
lc_times.append(timeit.timeit("lc.compute(box, particles)", number=10, globals=globals()))
naive_times.append(timeit.timeit("compute_distances(box, particles)", number=10, globals=globals()))
Running for log_N = 5
Running for log_N = 6
Running for log_N = 7
Running for log_N = 8
Running for log_N = 9
Running for log_N = 10
Running for log_N = 11
[7]:
fig, ax = plt.subplots()
ax.plot(2**log_Ns, lc_times, label="LinkCell")
ax.plot(2**log_Ns, naive_times, label="Naive Calculation")
ax.legend(fontsize=14)
ax.tick_params(axis='both', which='both', labelsize=14)
ax.set_xlabel("Number of particles", fontsize=14)
ax.set_ylabel("Time to compute (s)", fontsize=14);

Nearest Neighbors¶
One of the basic computations required for higher-level computations (such as the hexatic order parameter) is finding the nearest neighbors of a particle. This tutorial will show you how to compute the nearest neighbors and visualize that data.
The algorithm is straightforward:
for each particle i:
for each particle j in neighbor_cells(i):
r_ij = position[j] - position[i]
r = sqrt(dot(r_ij, r_ij))
l_r_array.append(r)
l_n_array.append(j)
# sort by distance
sort(n_array, r_array)
neighbor_array[i] = n_array[:k]
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 a solid at packing fractions \(\phi = 0.75\).
[1]:
from bokeh.io import output_notebook
output_notebook()
from bokeh.models import Legend
from bokeh.plotting import figure, output_file, show
from bokeh.layouts import gridplot
import numpy as np
import time
from freud import parallel, box, locality
parallel.setNumThreads(4)
import util
# Create hexagon vertices
verts = util.make_polygon(sides=6, radius=0.6204)
# Define colors for our system
c_list = ["#30A2DA", "#FC4F30", "#E5AE38", "#6D904F", "#9757DB",
"#188487", "#FF7F00", "#9A2C66", "#626DDA", "#8B8B8B"]
c_dict = dict()
c_dict[6] = c_list[0]
c_dict[5] = c_list[1]
c_dict[4] = c_list[2]
c_dict[3] = c_list[7]
c_dict[7] = c_list[4]
def render_plot(p, fbox):
# Display box
corners = util.box_2d_to_points(fbox)
p.patches(xs=[corners[:-1, 0]], ys=[corners[:-1, 1]],
fill_color=(0, 0, 0, 0), line_color="black", line_width=2)
p.legend.location = 'bottom_center'
p.legend.orientation = 'horizontal'
util.default_bokeh(p)
show(p)
[2]:
# Load the data
data_path = "data/phi065"
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))
n_frames = pos_data.shape[0]
Viewing our system¶
Before proceeding, we should probably view our system first. freud
does not make any assumptions about your data and is not specifically designed for any one visualization package. Here we use bokeh to render our system. Bokeh is not appropriate for real-time interaction with your simulation data, nor is it appropriate for 3D data, but is perfectly fine for rendering individual simulation frames, so we will use it here
[3]:
# Grab data from last frame
l_box = box_data[-1].tolist()
l_pos = pos_data[-1]
l_quat = quat_data[-1]
l_ang = 2*np.arctan2(np.copy(l_quat[:, 3]), np.copy(l_quat[:, 0]))
# Create box
fbox = box.Box.from_box(l_box)
side_length = max(fbox.Lx, fbox.Ly)
l_max = side_length / 2.0
l_max *= 1.1
# Take local vertices and rotate, translate into system coordinates
patches = util.local_to_global(verts, l_pos[:, :2], l_ang)
# Plot
p = figure(title="System Visualization",
x_range=(-l_max, l_max),
y_range=(-l_max, l_max))
p.patches(xs=patches[:, :, 0].tolist(), ys=patches[:, :, 1].tolist(),
fill_color=(42, 126, 187), line_color="black", line_width=1.5, legend="Hexagons")
render_plot(p, fbox)
By eye, we can see regions where the hexagons appear close-packed, as well as regions where there are vacancies. We will be using the Nearest Neighbor object to investigate this in our system.
The Nearest Neighbor object¶
This module will give the indices of the \(k\) particles which are nearest to another particle. freud
provides two different modes by which to compute the nearest neighbors, selected by the strict_cut
variable:
strict_cut=False
(default): The value forrmax
is expanded until every particle has \(k\) nearest neighborsstrict_cut=True
: thermax
value is not expanded, so that any “vacancies” in the number of neighbors found are filled withUINTMAX
strict_cut=False
¶
First we show how to use the strict_cut=False
mode to find the neighbors of a specific particle
[4]:
# Create freud nearest neighbor object
n_neigh = 6
nn = locality.NearestNeighbors(rmax=1.5, n_neigh=n_neigh, strict_cut=False)
# Compute nearest neighbors
nn.compute(fbox, l_pos, l_pos)
# Get the NeighborList
n_list = nn.nlist
# Get the neighbors for particle 0
pidx = 0
n_idxs = n_list.index_j[np.where(n_list.index_i == pidx)[0]]
# Get position, orientation for the central particle
center_pos = l_pos[np.newaxis, pidx]
center_ang = l_ang[np.newaxis, pidx]
# Get the positions, orientations for the neighbor particles
neigh_pos = np.zeros(shape=(n_neigh, 3), dtype=np.float32)
neigh_ang = np.zeros(shape=(n_neigh), dtype=np.float32)
neigh_pos[:] = l_pos[n_idxs]
neigh_ang[:] = l_ang[n_idxs]
# Create array of transformed positions
c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
n_patches = util.local_to_global(verts, neigh_pos[:, 0:2], neigh_ang)
# Create array of colors
center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
# Plot
p = figure(title="Nearest Neighbors Visualization",
x_range=(-l_max, l_max), y_range=(-l_max, l_max))
p.patches(xs=n_patches[:, :, 0].tolist(), ys=n_patches[:, :, 1].tolist(),
fill_color=neigh_color.tolist(), line_color="black", legend="Neighbors")
p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="Center")
render_plot(p, fbox)
Notice that nearest neighbors properly handles periodic boundary conditions.
We do the same thing below, but for a particle/neighbors not spanning the box.
[5]:
# Get the neighbors for particle 1000
pidx = 1000
segment_start = n_list.segments[pidx]
segment_end = segment_start + n_list.neighbor_counts[pidx]
n_idxs = n_list.index_j[segment_start:segment_end]
# Get position, orientation for the central particle
center_pos = l_pos[np.newaxis, pidx]
center_ang = l_ang[np.newaxis, pidx]
# Get positions, orientations for the neighbors and one non-neighbor
neigh_pos = l_pos[n_idxs]
neigh_ang = l_ang[n_idxs]
non_neigh_pos = neigh_pos[np.newaxis, -1]
non_neigh_ang = neigh_ang[np.newaxis, -1]
# Create array of transformed positions
c_patches = util.local_to_global(verts, center_pos[:, :2], center_ang)
n_patches = util.local_to_global(verts, neigh_pos[:, :2], neigh_ang)
non_n_patches = util.local_to_global(verts, non_neigh_pos[:, :2], non_neigh_ang)
# Create array of colors
center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
# Color the last particle differently
non_neigh_color = np.array([c_list[-2] for _ in range(non_neigh_pos.shape[0])])
# Plot
p = figure(title="Nearest Neighbors Visualization",
x_range=(-l_max, l_max), y_range=(-l_max, l_max))
p.patches(xs=n_patches[:, :, 0].tolist(), ys=n_patches[:, :, 1].tolist(),
fill_color=neigh_color.tolist(), line_color="black", legend="Neighbors")
p.patches(xs=non_n_patches[:, :, 0].tolist(), ys=non_n_patches[:, :, 1].tolist(),
fill_color=non_neigh_color.tolist(), line_color="black", legend="Non-neighbor")
p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="Center")
render_plot(p, fbox)
Notice that while freud
found the 6 nearest neighbors, one of the particles isn’t really in a neighbor position (which we have colored purple). How do we go about finding particles with a deficit or surplus of neighbors?
strict_cut=True
¶
Now for strict_cut=True
. This mode allow you to find particles which have fewer than the specified number of particles. For this system, we’ll search for 8 neighbors, so that we can display particles with both a deficit and a surplus of neighbors.
[6]:
# Create freud nearest neighbors object
n_neigh = 8
rmax = 1.7
nn = locality.NearestNeighbors(rmax=rmax, n_neigh=n_neigh, strict_cut=True)
# Compute nearest neighbors
nn.compute(fbox, l_pos, l_pos)
# Get the neighborlist
n_list = nn.nlist
# Get the number of particles
num_particles = nn.n_ref
p = figure(title="Nearest Neighbors visualization",
x_range=(-l_max, l_max), y_range=(-l_max, l_max))
for k in np.unique(n_list.neighbor_counts):
# Find particles with k neighbors
c_idxs = np.where(n_list.neighbor_counts == k)[0]
center_pos = l_pos[c_idxs]
center_ang = l_ang[c_idxs]
c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
render_plot(p, fbox)
Visualize each set of k values independently¶
[7]:
for k in np.unique(n_list.neighbor_counts):
p = figure(title="Nearest Neighbors: k={}".format(k),
x_range=(-l_max, l_max), y_range=(-l_max, l_max))
# Find particles with k neighbors
c_idxs = np.where(n_list.neighbor_counts == k)[0]
center_pos = l_pos[c_idxs]
center_ang = l_ang[c_idxs]
c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
patches = util.local_to_global(verts, l_pos[:, 0:2], l_ang)
center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
p.patches(xs=patches[:, :, 0].tolist(), ys=patches[:, :, 1].tolist(),
fill_color=(128, 128, 128, 0.1), line_color=(0, 0, 0, 0.1), legend="other")
p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
render_plot(p, fbox)
[8]:
for k in np.unique(n_list.neighbor_counts):
p = figure(title="Nearest Neighbors: k={}".format(k),
x_range=(-l_max, l_max), y_range=(-l_max, l_max))
# Find particles with k neighbors
c_idxs = np.copy(np.where(n_list.neighbor_counts == k)[0])
center_pos = l_pos[c_idxs]
center_ang = l_ang[c_idxs]
neigh_pos = np.zeros(shape=(k*len(c_idxs), 3), dtype=np.float32)
neigh_ang = np.zeros(shape=(k*len(c_idxs)), dtype=np.float32)
for i, pidx in enumerate(c_idxs):
# Create a list of positions, angles to draw
segment_start = n_list.segments[pidx]
segment_end = segment_start + n_list.neighbor_counts[pidx]
n_idxs = n_list.index_j[segment_start:segment_end]
for j, nidx in enumerate(n_idxs):
neigh_pos[k*i+j] = l_pos[nidx]
neigh_ang[k*i+j] = l_ang[nidx]
c_patches = util.local_to_global(verts, center_pos[:, 0:2], center_ang)
n_patches = util.local_to_global(verts, neigh_pos[:, 0:2], neigh_ang)
patches = util.local_to_global(verts, l_pos[:, 0:2], l_ang)
center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
p.patches(xs=patches[:, :, 0].tolist(), ys=patches[:, :, 1].tolist(),
fill_color=(128, 128, 128, 0.1), line_color=(0, 0, 0, 0.1), legend="other")
p.patches(xs=n_patches[:, :, 0].tolist(), ys=n_patches[:, :, 1].tolist(),
fill_color=neigh_color.tolist(), line_color="black", legend="neighbors")
p.patches(xs=c_patches[:, :, 0].tolist(), ys=c_patches[:, :, 1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
render_plot(p, fbox)