freud: A Software Suite for High Throughput Analysis of Particle Simulation Data

TL;DR

freud is a versatile Python library that efficiently analyzes particle simulation data using parallelized C++ routines, supporting various methods and flexible data input for high-throughput and adaptive analysis.

Contribution

It introduces a flexible, high-performance Python package that directly analyzes numerical data arrays, supporting a wide range of particle analysis methods with a uniform API.

Findings

Supports standard and advanced particle analysis methods.

Compatible with on-the-fly analysis in simulation workflows.

Efficiently handles large datasets on various computing platforms.

Abstract

The freud Python package is a powerful library for analyzing simulation data. Written with modern simulation and data analysis workflows in mind, freud provides a Python interface to fast, parallelized C++ routines that run efficiently on laptops, workstations, and supercomputing clusters. The package provides the core tools for finding particle neighbors in periodic systems, and offers a uniform API to a wide variety of methods implemented using these tools. As such, freud users can access standard methods such as the radial distribution function as well as newer, more specialized methods such as the potential of mean force and torque and local crystal environment analysis with equal ease. While many comparable tools place a heavy emphasis on reading and operating on trajectory file formats, freud instead accepts numerical arrays of data directly as inputs. By remaining agnostic to its data source, freud is suitable for analyzing any coarse-grained particle simulation, regardless of the original data representation or simulation method. When used for on-the-fly analysis in conjunction with scriptable simulation software such as HOOMD-blue, freud enables smart simulations that adapt to the current state of the system, allowing users to study phenomena such as nucleation and growth.

Equations

$$MSD(i,m)=\frac{1}{N_f-m}\sum _{k=0}^{N_f-m-1}\Vert ({\vec{r}}_i(k+m)-{\vec{r}}_i(k)){\Vert }^2$$

$$MSD(m)=\frac{1}{N_p}\sum _{i=1}^{N_p}MSD(i,m)$$

$$Z$$

$$=\int d\Delta {\xi }_{12}J(\Delta {\xi }_{12})e^{-\beta U(\Delta {\xi }_{12})}{e}^{-\beta {\tilde{F}}_{12}(\Delta {\xi }_{12})}$$

$$Z\equiv \int d\Delta {\xi }_{12}{e}^{-\beta F_{12}(\Delta {\xi }_{12})}$$

$$\beta F_{12}(\Delta {\xi }_{12})=\beta U(\Delta {\xi }_{12})-\log J(\Delta {\xi }_{12})+\beta {\tilde{F}}_{12}(\Delta {\xi }_{12})$$

$$F_{12}(\Delta {\xi }_{12})=-k_{B}T\log (H(d(\Delta {\xi }_{12}))J(\Delta {\xi }_{12}))+{\tilde{F}}_{12}(\Delta {\xi }_{12})$$

Full text

PROGRAM SUMMARY

Program Title: freud

Licensing provisions: BSD 3-Clause

Programming language: Python, C++

Nature of problem: Simulations of coarse-grained, nano-scale, and colloidal particle systems typically require analyses specialized to a particular system. Certain more standardized techniques – including correlation functions, order parameters, and clustering – are computationally intensive tasks that must be carefully implemented to scale to the larger systems common in modern simulations.

Solution method: freud performs a wide variety of particle system analyses, offering a Python API that interfaces with many other tools in computational molecular sciences via NumPy array inputs and outputs. The algorithms in freud leverage parallelized C++ to scale to large systems and enable real-time analysis. The library’s broad set of features encode few assumptions compared to other analysis packages, enabling analysis of a broader class of data ranging from biomolecular simulations to colloidal experiments.

Unusual features:

freud provides very fast, parallel implementations of standard analysis methods like RDFs and correlation functions. 2. 2.

freud includes the reference implementation for the potential of mean force and torque (PMFT). 3. 3.

freud provides various novel methods for characterizing particle environments, including the calculation of descriptors useful for machine learning.

Additional comments: The source code is hosted on GitHub (https://github.com/glotzerlab/freud), and documentation is available online (https://freud.readthedocs.io/). The package may be installed via pip install freud-analysis or conda install -c conda-forge freud.

1 Introduction

Molecular simulation is a crucial pillar in the investigation of scientific phenomena. Increased computational resources, better algorithms, and new hardware architectures have made it possible to simulate complex systems over longer timescales than ever before [1, 2, 3, 4, 5]. The sheer volume of data necessitates computationally efficient analysis tools, while the diversity of data requires flexible tools that can be adapted for specific systems. Additionally, to support scientists with limited prior computing experience, tools must be usable without extensive knowledge of the underlying code.

Numerous software packages that satisfy these requirements have been developed in recent years. Tools such as MDTraj [6], MDAnalysis [7], LOOS [8], MMTK [9], and VMD [10] provide efficient implementations of various standard analysis methods. Although powerful, such tools are generally limited in scope to all-atom simulations, particularly biomolecular simulations. This focus is manifested not only through the features these tools provide, but also in their general design philosophies.

Perhaps the most pronounced characteristic of such tools is a strong emphasis on trajectory management, which includes parsing trajectory files and supporting extensive topology selection features to enable, for instance, selecting all residues or atoms in a protein backbone. Although such tools are crucial for working with topologies in atomistic simulations, they are frequently cumbersome for working with coarse-grained simulation data where the trivial selection (all particles in the system) is the most common selection for various analyses. Moreover, such topology selection tools make assumptions that are inappropriate for non-atomistic systems: “bonding” in colloidal systems, for instance, is typically based on whether two particles are found to be in the same neighborhood by some distance-based metric, not by the presence of a true chemical bond. Since such determination of nearest neighbors is highly dynamic and parameter-dependent, it must be calculated on-the-fly and cannot be stored in a trajectory.

Another inconvenient but almost universal implementation choice is to directly tie analysis methods to trajectories by writing code that acts directly on some in-memory representation of a trajectory. This direct linkage is generally inflexible because it inhibits pre-processing of the data before running the analysis, which is often crucial to analyzing more specialized systems. More importantly, existing tools emphasize implementations of highly specific analyses involving, for instance, hydrogen bonding and protein secondary structure (using, e.g., DSSP [11]), which are far less useful for analyzing non-biomolecular systems. The predominant analyses of coarse-grained, colloidal-scale, or nanoparticle simulations usually involve measurements like numbers of nearest neighbors, diffraction patterns, or bond-orientational order parameters. These analyses bear little relation to the analyses performed for atomistic systems. These considerations suggest a need for a different type of analysis package that offers different methods than most existing tools.

In this paper we introduce freud, an open-source simulation analysis toolkit that addresses these needs. All inputs to and outputs from freud are numerical arrays of data, and the package makes no reference to predefined notions of atoms or molecules. As a result, freud can analyze particle-based data from both experiments and simulations regardless of the specific tools, methods, or software that were used to generate it. The package provides a Python Application Programming Interface (API) for accessing fast methods implemented in C++, and it implements numerous specific methods such as radial distribution functions and correlation functions that are common in the field of soft-matter physics (see fig. 1). Prior works have used freud for: determining spatial correlation functions and potentials of mean force and torque (PMFTs) in two dimensions [1]; calculating Steinhardt order parameters for identifying solid-like particles [12, 13]; computing spherical harmonics for machine learning on crystal structures [14]; optimizing pair potentials for designing complex crystals [15]; calculating strain fields by finding neighbors of particles against a uniform grid [16]; finding PMFTs in depletion-mediated self-assembly of hard cuboctahedra [17]; measuring rotational degrees of freedom in entropically ordered systems [18]; umbrella sampling of solid-solid phase transitions using Steinhardt order parameters [19]; evaluating PMFTs in analysis of two-dimensional shape allophiles [20]; and more. The freud library is designed to work well with coarse-grained particle models, such as those used in simulations of anisotropic nanoparticles, colloidal crystals, and polymers, and its methods are particularly useful for studies of phase transitions and critical phenomena in such systems. The package is likely to be of greatest interest to scientific communities in materials science, chemical engineering, and physics, though many of its analysis methods would be useful in generic particle systems. The freud library also integrates well into the scientific Python ecosystem, especially in data pipelines for machine learning and visualization [21].

The paper is organized as follows. We first address the core design principles that went into building freud in section 2. Section 3 focuses more specifically on the details of the code, including information on class structures. Section 4 describes the various analysis methods in freud and details their uses. Finally, in section 5 we provide some example code demonstrating the usage of freud111The code for these examples and many others is available at https://github.com/glotzerlab/freud-examples and in our online documentation.. The figures in this paper are rendered using Matplotlib [22] unless otherwise noted.

2 Design

Many of the best known tools for analyzing molecular simulations are built into either simulation toolkits (such as LAMMPS [23], GROMACS [24], or the cpptraj [25] plugin to Amber [26]) or visualization toolkits (such as VMD [10], PyMOL [27], or OVITO [28]). Although most of these have introduced varying degrees of scripting support over the years, the analyses built into simulation toolkits are primarily focused on performing one-shot analyses on trajectory files directly from the command line. The visualization toolkits tend to have more full-featured scripting interfaces, but they are frequently difficult (if not impossible) to use outside their own sandboxed environments, complicating or even prohibiting integration with other tools. More recently, many newer tools such as MDTraj [6], MDAnalysis [7], LOOS [8], and Pteros [29] have aimed to decouple analysis from simulation and visualization, making scriptability a primary focus to increase flexibility. Among such tools, Python is the most common language of choice due to its ease of use and the fact that it can be naturally extended with high performance languages like C, C++, and FORTRAN.

freud follows in the footsteps of these tools, providing a full-featured Python API to access all of its routines. However, while most other tools focus on calculating properties involving molecular topologies, freud is fundamentally designed for analyzing the local neighborhoods of particles, particularly where such local analyses provide global insight about the system. Such analyses are typically far more varied and system-dependent than the standard analyses of molecular topologies and therefore require more flexible tools. To meet this need, freud eschews any form of trajectory object encoding system topology and is instead designed such that each analysis method is an independent Python class that performs computations directly on NumPy arrays [30] of data.

This design makes it possible to use a far wider range of data with freud than is possible with tools that are tied to simulation trajectories. For instance, calculating Voronoi diagrams and computing spatial correlation functions with freud is possible for essentially arbitrary spatial data, not just the result of a molecular simulation. Another major benefit is that NumPy arrays (the de facto standard for numerical data in Python) can be easily passed between multiple tools, making freud equally easy to use for one-off analyses or as part of a larger analysis pipeline involving many steps and using various software packages. As a result, freud is a much more flexible choice both for analyzing disparate sources of data and for incorporating into Python workflows. For example, most of freud’s analyses can be used within the OVITO visualization environment for real-time visualization with almost no noticeable performance cost.

Producing such array data from simulation trajectories for input to freud is straightforward because high quality file parsers with Python APIs already exist for all common trajectory file formats. Through integration with tools like MDAnalysis, GSD [31], and garnett [32], freud can be used with data from over 25 different file formats, including common formats like DCD, XTC, and TRJ files. freud integrates with the trajectory objects produced by many of these tools, but if necessary, users can read trajectories into arrays and modify them before passing the data to freud for analysis. By using data read by other tools, freud’s analyses can be made aware of molecular topology if needed, but only when the analysis method requires such information. Similarly, since the outputs of freud’s analyses are also NumPy arrays, they can be passed to almost any tool in the scientific Python software stack. For example, constructing a Pandas [33] DataFrame from the outputs of any freud analysis requires just one line of code, and it immediately enables writing the output to text, CSV, or HDF5, or saving into an SQL database.

Beyond differences in trajectory and data handling, the most significant design choice in freud stems from the most common pattern followed by its many analysis methods. Since the first task in characterizing local neighborhoods is often the identification of neighboring particles, freud provides efficient methods for finding neighbors in arbitrary system geometries. The nearest-neighbor finding routines are designed to be as fast and flexible as possible, supporting various algorithms optimized for different system configurations and offering different criteria for neighbor selection. In general, queries can be based on either a cutoff distance or a desired number of nearest neighbors. These tools are optimized to provide cheap access to neighbors even in highly performance-critical loops in C++ analysis routines, but the package’s system box representation and neighbor finding tools also have Python APIs, so users can implement custom analyses directly in Python (for an example, see section 5.4).

The analysis methods in freud are essentially independent tools that make use of these objects to efficiently perform various calculations. These features are all presented with a common API, easing the transition between the different types of analyses needed for different simulations. All methods in freud are accelerated through extensive parallelization.

3 Implementation

The freud package is entirely object-oriented, with two core C++ classes: the Box class, which encapsulates all logic associated with periodicity in arbitrary triclinic boxes (boxes with 3 linearly independent basis vectors); and the NeighborQuery class, which facilitates efficiently finding, storing, and iterating over nearest neighbors. In keeping with the Python ethos, box objects in freud may be constructed from a variety of inputs. Any method in freud that accepts a box object also accepts a number of objects that can be interpreted as a box, such as a $3\times 3$ NumPy array of box vectors or a list of three numbers representing the edge lengths of an orthorhombic box. There are two subclasses of NeighborQuery in freud that each implement different neighbor search algorithms: one implements a bounding volume hierarchy (BVH) [37], while the other implements a cell list [34]. The NeighborList class is a lightweight storage mechanism for NeighborQuery results that accelerates performing multiple analyses on the same set of neighbor pairs.

The analysis methods in freud are encapsulated by Compute classes, which are loosely defined as classes providing a compute method that populates class attributes after performing some computation. Compute classes, such as the density module’s RDF class, are usually configured with constructor arguments, after which they can be used multiple times to perform distinct calculations. Some classes in freud (e.g. the radial distribution function (RDF), PMFT, or bond-orientational order diagram) represent histogram-like quantities, and therefore allow the user to specify reset=False as an argument to compute in order to accumulate and average data over many calls.

Compute classes can be divided into two groups, those that depend on finding neighbors and those that do not. A majority of calculations in freud require neighbors, and the compute methods of such classes all share two arguments, system and neighbors (in addition to analysis-specific arguments like particle orientations for PMFTs; such arguments are also typically NumPy arrays). The system parameter accepts a NeighborQuery or any object that can be interpreted as a tuple (box, points), where the box is any valid box-like object (as described above) and the points argument is anything that can be interpreted as an $N\times 3$ NumPy array of positions. Valid systems include simulation frame objects from tools such as MDAnalysis, GSD, garnett, OVITO, or the particle simulation engine HOOMD-blue [35, 36, 37].

When performance is critical, providing a NeighborQuery object is advantageous because many compute methods can reuse these neighbor search data structures. For all other system inputs, freud internally constructs a NeighborQuery if the compute method requires neighbor pairs. The neighbors argument is a dictionary of query arguments, such as dict(num_neighbors=12) or dict(

r_max=3.0) (the complete specification for freud’s Query API is provided in the documentation). Alternatively, users may precompute a NeighborList and provide it as the neighbors. In this case, whether system is a NeighborQuery or not has no impact on performance because the calculation will be carried out directly on the provided set of neighbor pairs and no additional spatial searches are required. Figure 2 shows a flowchart demonstrating how these classes and data structures are used.

Some methods in freud do not operate on neighboring pairs of particles. For those that still depend on particle positions (such as GaussianDensity), the first argument is still any valid system, but no neighbors are provided. Some methods do not depend on positions at all; for instance, the Nematic order parameter only requires particle orientations. In such cases, the user can simply pass that quantity alone to the calculation. This mode of operation is particularly useful when performing high-throughput analysis of large files; using file formats like GSD that permit reading only certain properties of the trajectory, users can minimize I/O operations by only reading the required arrays from memory.

All Compute classes use efficient, thread-parallel C++ implementations for performance-critical components. The Python bindings for these C++ classes are generated using Cython [38], and the C++ methods are mirrored in Python using thin Cython classes that dispatch calls to the underlying C++ class instances. The Cython classes have limited responsibilities: managing the memory of the underlying C++ instances, sanitizing inputs when necessary, and providing transparent access via memory views on C++ arrays.

The main exception to this design is the msd module, which is implemented in pure Python in freud. The mean squared displacement (MSD) is a measure of, on average, how far particles move in a given window of time. In a simulation trajectory of $N_{f}$ frames, the MSD of particle $i$ over a window of length $m$ frames is given by:

[TABLE]

Therefore, the total MSD is given by:

[TABLE]

Direct computation of the MSD is an $\mathcal{O}(N_{p}N_{f}^{2})$ operation, but by using a fast Fourier transform (FFT) this cost can be reduced to $\mathcal{O}(N_{p}N_{f}\log(N_{f}))$ [39]. When using this approach, the FFT is responsible for most of the computation time, and since packages like NumPy [30] and SciPy [40] already expose fast C and FORTRAN FFT routines to Python, freud simply leverages them directly and implements the rest of the MSD in pure Python.

Calculations in freud are generally parallelized over particles (e.g. the Nematic order parameter class) or over pairs of particles (e.g. computing inter-particle distances for an RDF with the RDF class). Both the number of particles and the number of particle-particle pairs increase with system size, ensuring that the work is load-balanced well among threads because the number of threads is much less than the number of particles or pairs. Parallelism in freud is accomplished using Intel Threading Building Blocks (TBB) [41]. Analysis routines are written as lambda functions operating on a particle or a pair of particles; freud provides wrappers that then automatically parallelize these functions appropriately using TBB. Modern compilers aggressively inline such lambda functions, thereby optimizing away any additional cost that could arise from the extra function calls. freud uses thread-local storage extensively to avoid any parallel writes to data containers. For histograms that accumulate over many frames of simulation data, freud performs reduction over thread-local containers lazily.

Currently, freud is at version 2.1.0 and supports Python versions 3.5.0 or greater. The package is distributed through the Python Package Index (PyPI) and the conda-forge channel of the Anaconda package manager [42], making it easy to install on any Unix-based operating system (e.g. Linux or macOS). Builds for the Windows operating system are also available on conda-forge. freud depends on NumPy and TBB libraries, which are automatically installed with freud. The library can also be compiled from source using a C++11 compliant compiler. Compilation requires NumPy and TBB headers as well as a Cython installation. Code documentation is written using Google-style docstrings rendered using Sphinx and hosted on ReadTheDocs. The freud library is released open source under the BSD 3-Clause License, and the source code is available in a GitHub repository [43]. Continuous integration testing is performed using CircleCI.

4 Features

4.1 General Utilities

The general utilities in freud are contained in two modules: box and locality. The box module contains the core Box class. The locality module contains the NeighborQuery abstract class, which defines the standardized query API. NeighborQuery results (neighboring particle pairs) can be obtained dynamically or stored in the NeighborList class provided by the locality module.

Box periodicity is built in at the lowest level of the NeighborQuery subclasses, which are highly optimized for this use case. The AABBQuery subclass implements a tree data structure of Axis-Aligned Bounding Boxes (AABBs), a type of BVH which greatly accelerates the process of finding particles’ neighbors [37, 13]. A second approach is implemented in the LinkCell subclass [34], which uses linked cell lists to find particle neighbors. Both of these classes can find neighboring particle pairs based on a distance cutoff or a desired number of neighbors, and both were adapted from HOOMD-blue.

As a performance benchmark, we compare freud’s AABBQuery class with the cKDTree implementation in SciPy [40]. As part of SciPy, this implementation is the most readily available alternative to AABBQuery. Figure 3 shows that freud’s AABBQuery routines clearly outperform the cKDTree as system sizes increase to thousands of points. Moreover, we note that while freud supports general triclinic boxes, SciPy’s cKDTree only supports periodic orthorhombic boxes (i.e. a cuboid, a rectangular prism where all angles are right angles).

In addition to these performance gains, the NeighborQuery objects in freud are designed to interoperate seamlessly with analysis routines. Since analyses in freud are written in C++, using a Python API to find neighbors and then pass them into other C++ routines would waste time in unnecessary memory transfers. Furthermore, while a Python API should make certain promises, such as sorting the order of the resulting neighbors, the analyses using neighbors simply loop over all pairs and therefore do not require any such extra work. To avoid this cost, the features of the NeighborQuery classes are directly accessible in C++ in the form of iterators that lazily produce neighbors. In practice, using the NeighborQuery classes in this manner speeds up computations by a factor of two or more depending on the system size. To make use of these iterators, developers implement analysis methods as lambda functions that are passed as arguments to freud’s internal TBB wrappers that apply these functions to neighbor pairs in parallel.

The final feature of the locality module is the Voronoi class, which uses the voro++ library [44] to generate Voronoi diagrams for systems of particles. Voronoi diagrams are a standard method for characterizing the local geometric arrangements in the system, and they also provide a parameter-free method for defining nearest-neighbor relationships [45]. The Voronoi class produces a NeighborList object that can then be used as the neighbors argument for other compute classes.

4.2 Analysis Modules

The remaining modules in freud are independent of one another and contain groups of classes that implement related features. While some of freud’s features are unique, many others are standard techniques. However, implementations of these methods commonly lack support for periodicity. For example, the SciPy library [40] has functions for computing Voronoi diagrams and correlation functions, but these are restricted to aperiodic systems.

The cluster module of freud can be used to find clusters of particles—where cluster membership is defined by neighbor bonds—and then compute properties of these clusters such as gyration tensors. The density module contains features for calculating radial distribution functions as well as spatial correlation functions of arbitrary quantities. Additionally, the density module can estimate local particle density and interpolate particle density onto a regular grid suitable for, e.g., computing discrete Fourier transforms. The interface module provides a quick tool for identifying interfaces between two mutually exclusive sets of points (e.g. a solid and a liquid phase). The msd module enables the calculation of mean squared displacements of particles over the course of a trajectory.

The order module is the most extensive one in freud, containing a large number of different order parameters commonly used to measure ordering and identify phase transitions in crystalline systems. Of particular note are the bond-orientational order parameters $Q_{l}$ and $W_{l}$ [47] and the cubatic order parameter [48] (see fig. 4). The module also contains the nematic order parameter for identifying orientationally ordered, translationally disordered phases, as well as a solid-liquid order parameter for identifying generic ordered phases [49].

The other features of freud are analysis methods developed by researchers in our group and not yet implemented anywhere else. In particular, the pmft and environment modules implement features unique to freud that we now discuss in greater detail.

4.3 Potentials of Mean Force and Torque

The potential of mean force and torque (PMFT) is a generalization of the classical potential of mean force (PMF) that was recently developed to quantify directional entropic forces that emerge in crowded systems [50, 51]. Given the canonical partition function as a function of particle positions $\{q\}$ and particle orientations $\{Q\}$, ref. [51] derives the PMFT by separating out a component corresponding to the relative coordinates of a pair of particles $\Delta\xi_{12}$:

[TABLE]

where $J$ is the Jacobian transforming to the local coordinate system and $\tilde{F}_{12}$ is the free energy of the other particles, which have been integrated over in eq. 4. The PMFT $F_{12}$ is defined by the relation

[TABLE]

Combining eqs. 4 and 5 gives an expression for the PMFT

[TABLE]

In hard particle systems governed exclusively by excluded volume interactions, the potential energy term becomes an infinite Heaviside function $H$ and the PMFT can be simplified to

[TABLE]

To contextualize the PMFT, we note that if in eq. 3 we redefine $\Delta\xi_{12}$ to only include the center-to-center distance of the pair of particles and otherwise follow the same steps, the resulting potential $F_{12}$ reduces to the classical PMF with the usual RDF relation $g(r)=e^{-\beta F_{12}(r)}$. This suggests that although the PMFT is a function of all degrees of freedom required to characterize the relative configuration of a pair of particles, examining a more limited coordinate system can still be informative. Figure 5 shows two examples of PMFTs that contain more information than a PMF without containing all available degrees of freedom. In the 2D PMFT in the left panel, the orientation of the second particle is ignored, but its angular position relative to the reference particle is sufficient to illustrate the clear preference for facet-to-facet alignment. Similarly, the right panel ignores the orientation of the second particle (which encodes three degrees of freedom in 3D, as represented by e.g. Euler angles), but once again the preference for facet-to-facet alignment is clear. For an example of a case where analyzing the full, high-dimensional PMFT is necessary, see ref. [53].

freud calculates the PMFT by accumulating a histogram of the configurations of all other particles and then taking the negative logarithm of the counts. PMFTs may be accumulated over many frames to generate smoother energy surfaces. This method of computing the PMFT is very similar to that of computing an RDF, so we compare their scaling behavior in a two-dimensional system in fig. 6. The calculations scale almost identically to many threads, with a constant scaling factor between them. There are two components contributing to the absolute difference in their performance: 1) the extra operations required to compute the orientation of the local coordinate system in the PMFT, and 2) the extra cost of binning in multiple dimensions.

We also tested performance as a function of the parameters of these two methods, namely the number of bins and the maximum interparticle distance. In the latter case, both methods show the expected quadratic growth, since the number of particles included in the calculation increases as the square of this cutoff distance. The behavior with respect to the number of bins is more interesting: this parameter has no effect on performance until it becomes sufficiently large, at which point performance begins to degrade. This degradation can be understood as the result of two things: 1) poor cache performance as the histograms become too large to fit in memory, and 2) increasing costs of reduction, which can eventually affect performance. Since the PMFT shown is a two-dimensional histogram, the number of bins that can be used along each dimension before experiencing this performance drop is commensurately smaller than can be used with the RDF; this effect would be even more pronounced for a three-dimensional PMFT.

4.4 Local Environments

The environment module provides methods for characterizing the local environments of particles that we now illustrate in greater detail.

4.4.1 Bond-Orientational Order Diagrams

The BondOrder class enables the calculation of bond-orientational order diagrams (BOODs) [54, 55, 56, 57, 58]. Inspired by the bond-orientational order parameters defined by Steinhardt et al. [47], BOODs characterize the local ordering of systems by calculating the vectors between all neighboring particles in a system and then projecting these vectors onto a sphere. One example of how BOODs can be used is to identify $n$-fold ordering in a system; in simple crystal structures with $n$-fold coordination, the BOOD will show $n$ peaks corresponding to the average location of nearest neighbors.

In addition to the standard BOOD calculation, the class offers some additional modes of operation that can be useful in specific cases. One mode involves finding the positions of nearest neighbors in the local coordinate system of a given particle rather than the global coordinate system, which can prevent misidentifying systems with multiple grains [54]. Another mode modifies the BOOD to help identify plastic crystals, which appear crystalline due to having translational order but lack orientational ordering. In this mode, the positions of the nearest neighbors of each particle are modified by the relative orientations of the neighboring particles, creating a BOOD in which positional ordering will no longer appear except when orientational ordering is also present.

4.4.2 Spherical Harmonic Descriptors

The BOOD is closely related to the Steinhardt order parameters $Q_{l}$ and $W_{l}$, which measure rotational order in a system using spherical harmonics [47]. While the BOOD is essentially a histogram of nearest-neighbor bonds, the Steinhardt order parameters take this one step further, measuring $l$-fold order by constructing scalar quantities from rotationally invariant combinations of spherical harmonics of degree $l$ calculated from the locations of nearest-neighbor bonds. However, spherical harmonic representations can also be used in a variety of different ways. For example, distinguishing different grains of the same crystal structure could be done using descriptors that are not rotationally invariant. Alternatively, we can often obtain rotationally-invariant descriptions of local environments for crystal structure identification via the principal axes of the moment of inertia tensor of the environments, or by using particle orientations of anisotropic particles [14, 48, 59]. To support such spherical harmonic analyses, the LocalDescriptors class in freud computes spherical harmonics characterizing particle neighborhoods. These harmonics can then be combined in arbitrary ways to generate custom descriptors of local particle environments. Such descriptors have proven useful in identifying multiple complex crystals (see fig. 7). One method for identifying these structures is to use the information contained in this array of spherical harmonics as a set of per-particle features in an artificial neural network (ANN) for structure classification [14].

4.4.3 Environment Matching

Methods like the spherical harmonic descriptors and the BOOD characterize ordering in systems by calculating system-averaged quantities from neighbor bonds. The EnvironmentCluster class takes a different approach by defining environments according to the nearest neighbors of each particle and performing point set registration to identify and cluster similar environments [60]. This type of analysis is particularly useful because it emphasizes local information for each particle. As a result, it can be used for tasks such as identifying different Wyckoff positions in a crystal. The complementary EnvironmentMotifMatch class can be used to match clusters to specific motifs, allowing deeper analysis of a given structural motif.

Since this method performs a direct pairwise comparison of all motifs, it is substantially more expensive than common methods used for structure identification. For instance, the performance is at least an order of magnitude slower than the implementation of Polyhedral Template Matching in OVITO and Common Neighbor Analysis. Unlike these methods, however, the environment matching algorithm does not depend on previous knowledge of possible structures, and instead infers possible structures entirely from the local motifs present in a system. Moreover, it can be tuned much more finely than the other methods, allowing not only the identification of crystal structures present, but also the precise identification of stacking faults like those found in fig. 8. As a result, it is a good complement to existing methods for identifying crystal structures in various systems.

4.4.4 Angular Separation

The AngularSeparation classes provide a way to characterize typical particle orientations in a system. The AngularSeparationGlobal class allows comparison of particle orientations to a set of reference orientations, which can be used to characterize orientational order relative to the reference input. This metric can be used as an order parameter for measuring orientational disorder in plastic crystals, which exhibit translational order and orientational disorder [61]. Alternatively, the AngularSeparationNeighbor computes minimum separation angles between neighboring particles, allowing a more fine-grained analysis of the orientational ordering in local motifs. Both of these methods account for symmetry by accepting an array of equivalent quaternions corresponding to all symmetry-preserving transformations of the particle (i.e. the particle’s point group).

4.5 Data Generation and Plotting

For the purposes of teaching via code examples and testing freud’s analyses, freud includes the freud.data module. It includes the UnitCell class for representing arbitrary unit cells with user-provided box vectors and basis positions. The UnitCell class includes class methods that generate common crystal structures like face-centered cubic, body-centered cubic, and simple cubic. The data module also includes a method for generating a random system with uniformly distributed points in a periodic box.

Many analysis modules in freud implement a plot() method, which can be used for rapid data visualization. The Compute classes (e.g. instances of freud.density.RDF) also define a repr_png() method that allows their data to be automatically plotted in IPython environments (such as Jupyter notebooks) using Matplotlib [22], when the last line in a code cell returns that analysis object.

5 Examples

In this section, we demonstrate the use of freud in conjunction with the broader scientific software ecosystem. The code for these examples and many others is available at https://github.com/glotzerlab/freud-examples.

5.1 RDF and MSD from LAMMPS simulation

Here, we consider the problem of calculating the RDF and the MSD of a system simulated using LAMMPS (version 5 Jun 2019) [23]. LAMMPS is a standard tool for particle simulation used in many fields, and it supports multiple output formats, including those used by other simulation codes (e.g., the DCD format from CHARMM [62] and the XTC format from GROMACS [24]). In this case, we demonstrate the case of using the output of a custom dump format in LAMMPS, which allows users to dump selected quantities into a text file. Although the default XYZ file format lacks sufficient information to calculate an MSD, the necessary particle image information can be included as shown.

⬇

import numpy as np

import freud

\par# For the MSD we also need images, which can be dumped

using the LAMMPS dump custom command as follows:

dump 2 all custom 100 output_custom.xyz x y z ix iy iz

\par# We read the number of particles, the system box, and the

particle positions into 3 separate arrays.

N = int(np.genfromtxt(

’output_custom.xyz’, skip_header=3, max_rows=1))

box_data = np.genfromtxt(

’output_custom.xyz’, skip_header=5, max_rows=3)

data = np.genfromtxt(

’output_custom.xyz’, skip_header=9,

invalid_raise=False)

\par# Remove the unwanted text rows

data = data[~np.isnan(data).all(axis=1)].reshape(-1, N, 6)

\parbox = freud.box.Box.from_box(

box_data[:, 1] - box_data[:, 0])

\par# We shift the system by half the box lengths to match the

freud coordinate system, which is centered at the origin.

Since all methods support periodicity, this shift is

simply for consistency but does not affect any analyses.

data[…, :3] -= box.L/2

rdf = freud.density.RDF(bins=100, r_max=4, r_min=1)

for frame in data:

rdf.compute(system=(box, frame[:, :3]), reset=False)

\parmsd = freud.msd.MSD(box)

msd.compute(

positions=data[:, :, :3], images=data[:, :, 3:])

\par# The object contains all the data we need to plot the RDF.

from matplotlib import pyplot as plt

plt.plot(rdf.bin_centers, rdf.rdf)

plt.title(’Radial Distribution Function’)

plt.xlabel(’$r$’)

plt.ylabel(’$g(r)$’)

plt.show()

If our trajectory was stored in a DCD file, we could modify our code above to read the input data using MDAnalysis (version 0.20.1):

⬇

reader = MDAnalysis.coordinates.DCD.DCDReader(

’output.dcd’)

rdf = freud.density.RDF(bins=100, r_max=4, r_min=1)

for frame in reader:

rdf.compute(system=frame, reset=False)

…

5.2 On-the-fly analysis with HOOMD-blue

A major strength of freud is that it can also be used for on-the-fly analysis. For example, freud can be used to terminate a simulation based on some additional condition, or to log a quantity at a higher frequency than we want to save the full system trajectory. In our previous example, we demonstrated the calculation of an RDF using freud. An RDF can be noisy when calculated with limited data, so we would like to average it over a large number of simulation frames; however, storing many frames can lead to unreasonably large simulation trajectory files. Using the simulation engine HOOMD-blue (v2.9.0), we can accumulate RDF data during a simulation without storing the entire output. Additionally, we show that we can log an order parameter over the course of the simulation:

⬇

import hoomd

from hoomd import hpmc

import freud

import numpy as np

\parhoomd.context.initialize()

system = hoomd.init.create_lattice(

hoomd.lattice.sc(a=1), n=10)

mc = hpmc.integrate.sphere(seed=42, d=0.1, a=0.1)

mc.shape_param.set(’A’, diameter=0.5)

\parrdf = freud.density.RDF(bins=50, r_max=4)

q6 = freud.order.Steinhardt(l=6)

\pardef calc_rdf(timestep):

snap = system.take_snapshot()

rdf.compute(system=snap, reset=False)

\pardef calc_Q6(timestep):

snap = system.take_snapshot()

q6.compute(system=snap,

neighbors={’num_neighbors’: 12})

return np.mean(q6.particle_order)

\par# Equilibrate the system before accumulating the RDF.

hoomd.run(1e4)

hoomd.analyze.callback(calc_rdf, period=100)

\parlogger = hoomd.analyze.log(filename=’output.log’,

quantities=[’q6’],

period=100,

header_prefix=’#’,

overwrite=True)

\parlogger.register_callback(’q6’, calc_Q6)

\parhoomd.run(1e4)

\par# Store the computed RDF in a file.

np.savetxt(’rdf.csv’,

np.vstack((rdf.bin_centers, rdf.rdf)).T,

delimiter=’,’, header=’r, g(r)’)

5.3 Analyzing Atomistic Trajectories from GROMACS

As discussed in sections 1 and 2, freud’s design focus differs from that of many similar tools in the lack of focus on trajectory management. The example below is based on a simulation trajectory of water molecules in a box generated using GROMACS (version 2020) [24]. We use MDTraj (version 1.9.3) [6] to read in an XTC trajectory file and then compute an RDF of the oxygen atoms in the water molecules using freud. In the process, we demonstrate how the sophisticated subsetting functionality offered by tools like MDTraj can be replicated with Python code, which is very useful when such subsets must be computed from coarse-grained trajectories with highly customized topology definitions that standard trajectory management tools cannot handle.

⬇

import mdtraj

import freud

import numpy as np

\partraj = mdtraj.load_xtc(

’output/prd.xtc’, top=’output/prd.gro’)

bins = 300

r_max = 1

r_min = 0.01

\par# Expression selection, a common feature of analysis tools

for atomistic systems, can be used to identify all

oxygen atoms.

oxygen_pairs = traj.top.select_pairs(’name O’, ’name O’)

\par# We can directly use the above selection in freud.

oxygen_indices = traj.top.select(’name O’)

\par# Alternatively, we can subset directly using Python logic.

Such selectors require the user to define the nature of

the selection, but can be more precisely tailored to a

specific system.

oxygen_indices = [atom.index for atom in traj.top.atoms

if atom.name == ’O’]

\parfreud_rdf = freud.density.RDF(

bins=bins, r_min=r_min, r_max=r_max)

for system in zip(np.asarray(traj.unitcell_vectors),

traj.xyz[:, oxygen_indices, :]):

freud_rdf.compute(system, reset=False)

\par# We can plot these RDFs to verify that they are equivalent.

from matplotlib import pyplot as plt

fig, ax = plt.subplots()

ax.plot(freud_rdf.bin_centers, freud_rdf.rdf, ’o’,

label=’freud’, alpha=0.5)

ax.plot(*mdtraj_rdf, ’x’, label=’mdtraj’, alpha=0.5)

ax.set_xlabel(’$r$’)

ax.set_ylabel(’$g(r)$’)

ax.set_title(’Radial Distribution Function’)

ax.legend()

5.4 Common Neighbor Analysis

Common Neighbor Analysis (CNA) [63] is a standard technique for analyzing the local neighborhoods of particles in a crystal. The method involves a classification of local neighborhoods based on a number of features. Using freud’s NeighborList, however, the method is straightforward to implement in Python.

We first consider the simpler problem of identifying all common neighbors between any pair of points. This is equivalent to searching for the second-nearest neighbor pairs, which can be done using freud as follows (note that this code is primarily written for clarity and could easily be optimized):

⬇

import freud

import numpy as np

from collections import defaultdict

\par# Use a face-centered cubic (fcc) system.

box, points = freud.data.UnitCell.fcc().generate_system(4)

aq = freud.AABBQuery(box, points)

query = aq.query(

points, {’num_neighbors’: 12, ’exclude_ii’: True})

nl = query.toNeighborList()

\par# Get all sets of common neighbors.

common_neighbors = defaultdict(list)

for i, p in enumerate(points):

selection1 = nl.query_point_indices == i

for j in nl.point_indices[selection1]:

selection2 = nl.query_point_indices == j

for k in nl.point_indices[selection2]:

if i != k:

common_neighbors[(i, k)].append(j)

Our dictionary common_neighbors now contains lists of common neighbors j for every pair of points (i, k). This information could itself be useful for performing some analysis on the system. If we are interested in actually implementing CNA, then we need to use this information to build local graphs, which we can do with the networkx (version 2.4) [64] Python package. Combined with the code above, the CNA algorithm can be implemented as follows:

⬇

import networkx as nx

from collections import Counter

\pardiagrams = defaultdict(list)

particle_counts = defaultdict(Counter)

\parfor (a, b), neighbors in common_neighbors.items():

Build up the graph of connections between the

common neighbors of a and b.

g = nx.Graph()

for i in neighbors:

for j in set(nl.point_indices[

nl.query_point_indices == i]

).intersection(neighbors):

g.add_edge(i, j)

\par# Define the four identifiers for a CNA diagram:

1. 1 if particles are bonded, 0 if not.

2. Number of shared neighbors.

3. Number of bonds between shared neighbors.

4. Index guaranteeing diagram uniqueness.

diagram_type = 2 - int(

b in nl.point_indices[nl.query_point_indices == a])

key = (diagram_type, len(neighbors),

g.number_of_edges())

If we’ve seen any neighborhood graphs with this

signature, we explicitly check if the two graphs are

identical to determine whether to save this one.

Otherwise, we always add the new graph.

if key in diagrams:

isomorphs = [

nx.is_isomorphic(g, h) for h in diagrams[key]]

if any(isomorphs):

idx = isomorphs.index(True)

else:

diagrams[key].append(g)

idx = diagrams[key].index(g)

else:

diagrams[key].append(g)

idx = diagrams[key].index(g)

cna_signature = key + (idx,)

particle_counts[a].update([cna_signature])

In this code, we are looping over all pairs of previously identified second neighbor shells, and finding bonds between the common neighbors of these pairs. The graph of these bonds then uniquely identifies a new environment.

6 Conclusion

freud is a high-performance Python library for analyzing particle simulations. Among simulation analysis packages, freud is unique due to its emphasis on coarse-grained simulations and its flexibility. Its high-performance C++ back-end makes freud a suitable solution for large-scale, high-throughput simulation analysis, while its simple, compact API is highly amenable to integration with other tools for, e.g., machine learning applications. The package’s API also promotes the prototyping of new analyses directly in Python, and the intuitive design of freud’s internals ensures that translating such analyses into C++ is a relatively painless process.

The package’s design is general enough to work with any particle-based system. However, freud is primarily targeted at communities of materials scientists, chemical engineers, and physicists analyzing molecular dynamics and Monte Carlo for which existing tools are too specialized to be convenient. Since it makes no assumptions about the types of its input data or the system topology, freud can be used with arbitrary simulation outputs based on topologies defined by the user. As a result, freud can find wide use in these areas to simplify workflows that require consideration of periodic systems without the complexity associated with specific atomistic features. Contributions to this open-source toolkit are highly encouraged as new methods are developed in future research applications.

Acknowledgments

Support for the design and development of freud has evolved over time and with programmatic research directions. Conceptualization and early implementations were supported in part by the DOD/ASD(R&E) under Award No. N00244-09-1-0062 and also by the National Science Foundation, Integrative Graduate Education and Research Traineeship, Award # DGE 0903629 (E.S.H. and M.P.S.). A majority of the code development including all public code releases was supported by the National Science Foundation, Division of Materials Research under a Computational and Data-Enabled Science & Engineering Award # DMR 1409620 (2014-2018) and the Office of Advanced Cyberinfrastructure Award # OAC 1835612 (2018-2021). V.R. holds the 2019-2020 J. Robert Beyster Computational Innovation Graduate Fellowship at the Uiniversity of Michigan. B.D. acknowledges fellowship support from the National Science Foundation under ACI-1547580, S212: Impl: The Molecular Sciences Software Institute and an earlier National Science Foundation Graduate Research Fellowship Grant # DGE 1256260 (2016-2019) [65, 66]. M.P.S. also acknowledges support from the University of Michigan Rackham Predoctoral Fellowship program. Computational resources and services supported in part by Advanced Research Computing at the University of Michigan, Ann Arbor.

Declarations of interest: none.