diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 1288fd5..f49e1c3 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -35,7 +35,7 @@ jobs: with: path: freud-examples # Important! All steps that use the conda environment should define "shell: bash -l {0}". - - uses: conda-incubator/setup-miniconda@v2 + - uses: conda-incubator/setup-miniconda@v3 with: python-version: ${{ matrix.python-version }} mamba-version: "*" @@ -81,7 +81,7 @@ jobs: # notify developers if the scheduled check fails - name: Slack notification if: ${{ github.event_name == 'schedule' && (failure() || cancelled()) }} - uses: 8398a7/action-slack@v3.15.0 + uses: 8398a7/action-slack@v3.15.1 with: status: ${{ job.status }} fields: workflow,job,message,commit diff --git a/data/RAD.png b/data/RAD.png new file mode 100644 index 0000000..053a2dc Binary files /dev/null and b/data/RAD.png differ diff --git a/data/SANN.png b/data/SANN.png new file mode 100644 index 0000000..4eac2b2 Binary files /dev/null and b/data/SANN.png differ diff --git a/examples/Calculating Strain via Voxelization.ipynb b/examples/Calculating Strain via Voxelization.ipynb index c36a086..923e202 100644 --- a/examples/Calculating Strain via Voxelization.ipynb +++ b/examples/Calculating Strain via Voxelization.ipynb @@ -67,7 +67,9 @@ "\n", " match = re.match('^Lattice=\".*\"', header)\n", " boxstring = match.group(0)\n", - " boxes.append(np.array(str.split(boxstring[9:-1]), dtype=np.float).reshape((3, 3)).T)\n", + " boxes.append(\n", + " np.array(str.split(boxstring[9:-1]), dtype=np.float32).reshape((3, 3)).T\n", + " )\n", " ff.close()\n", "\n", "# find the average box\n", @@ -226,7 +228,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -240,7 +242,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.3" + "version": "3.10.10" }, "widgets": { "state": { diff --git a/examples/Visualizing 3D Voronoi and Voxelization.ipynb b/examples/Visualizing 3D Voronoi and Voxelization.ipynb index 79fa5d1..10de92f 100644 --- a/examples/Visualizing 3D Voronoi and Voxelization.ipynb +++ b/examples/Visualizing 3D Voronoi and Voxelization.ipynb @@ -210,7 +210,7 @@ "outputs": [], "source": [ "aq = freud.AABBQuery(box, positions)\n", - "voxel_neighbors = -np.ones(len(cubic_grid), dtype=np.int)\n", + "voxel_neighbors = -np.ones(len(cubic_grid), dtype=np.int32)\n", "for i, j, distance in aq.query(cubic_grid, {\"num_neighbors\": 1}):\n", " voxel_neighbors[i] = j" ] @@ -246,7 +246,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -260,7 +260,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.3" + "version": "3.10.10" } }, "nbformat": 4, diff --git a/module_intros/locality.Filter.ipynb b/module_intros/locality.Filter.ipynb new file mode 100644 index 0000000..e9095ea --- /dev/null +++ b/module_intros/locality.Filter.ipynb @@ -0,0 +1,412 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "80d44c4c", + "metadata": {}, + "source": [ + "# freud.locality.Filter\n", + "\n", + "In this notebook, we introduce the neighborlist filter concept, discuss the filtering methods implemented in freud, and demonstrate how to use them with good efficiency." + ] + }, + { + "cell_type": "markdown", + "id": "f78106b4", + "metadata": {}, + "source": [ + "## What is a NeighborList Filter?\n", + "\n", + "A neighborlist filter is a class which removes some bonds from a pre-defined neighborlist. A neighborlist filter can be thought of as a function for which an unfiltered `NeighborList` is given as the input and a filtered version of the input `NeighborList` is the output. Freud already has some methods in the `NeighborList` class like `filter` and `filter_r` which can achieve this concept, but often the full system definition is needed to decide which bonds to remove. This is often the case when typical neighbor finding algorithms which use `r_max` or `num_neighbors` label two particles as neighbors even though they are blocked by another particle in between them. Each of _freud_'s neighborlist filters defines the concept of \"blocking\" or being \"blocked\" differently and removes bonds between particles which are blocked by another particle." + ] + }, + { + "cell_type": "markdown", + "id": "5eb827ee", + "metadata": {}, + "source": [ + "The _freud_ library's neighborlist filter algorithms are the same as some well-known parameter-free neighbor finding methods like the Solid Angle Nearest Neighbor (SANN) method and the Relative Angular Distance (RAD) method. Each of these neighbor finding methods was created to find a set of neighbors in which no neighbor is blocked by any other neighbor. In that sense, _freud_'s neighborlist filters can be used to *find* neighbors as well as to *filter* neighbors. The difference between finding and filtering depends solely on the meaning attributed to the unfiltered neighborlist.\n", + "\n", + "The following sections give descriptions of the algorithms used by each neighborlist filter in _freud_ along with code examples which illustrate proper usage." + ] + }, + { + "cell_type": "markdown", + "id": "445531c4", + "metadata": {}, + "source": [ + "## The Solid Angle Nearest Neighbor (SANN) Method\n", + "\n", + "With the [Solid Angle Nearest Neighbor (SANN) method](https://doi.org/10.1063/1.4729313), we look for enough nearby neighbors to fully occupy the $4 \\pi$ solid angle distribution surrounding a particle and label all neighbors further away as blocked. Strictly speaking, we consider neighbors of a particle $i$ to consist of the nearest (i.e., closest) $m$ particles ${j}$ in neighborhood defined by shell radius $R_i^{(m)}$ such that the sum of their solid angles associated with $\\theta_{i,j}$ equals $4 \\pi$ :\n", + "$$ 4\\pi = \\sum_{j=1}^m 2 \\pi (1-\\cos{\\theta_{i,j}}) = \\sum_{j=1}^m 2 \\pi (1- r_{i,j}/R_i^{(m)})$$\n", + "\n", + "\n", + "\"Alt\n", + "The `freud.locality.FilterSANN` class implements this method, and an example usage is shown below:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "3164e5e7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import freud\n", + "\n", + "# make the system to operate on\n", + "N = 1000\n", + "L = 10\n", + "box, points = freud.data.make_random_system(L, N)\n", + "\n", + "# create the unfiltered neighborlist\n", + "nlist = (\n", + " freud.locality.AABBQuery(box, points)\n", + " .query(points, dict(r_max=4.9, exclude_ii=True))\n", + " .toNeighborList()\n", + ")\n", + "\n", + "# make the FilterSANN and call compute\n", + "sann = freud.locality.FilterSANN()\n", + "sann.compute((box, points), neighbors=nlist)\n", + "\n", + "# access the filtered neighborlist as a property at the end\n", + "sann.filtered_nlist" + ] + }, + { + "cell_type": "markdown", + "id": "f217ac3a", + "metadata": {}, + "source": [ + "Instead of having to do the neighbor query explicitly in your python script, freud will do it automatically if query arguents are passed to the `neighbors` argument of the filter's `compute` method. The following call to `compute` is equivalent to the call in the previous script:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "77e5a753", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# call compute\n", + "sann.compute((box, points), neighbors=dict(r_max=4.9, exclude_ii=True))\n", + "\n", + "# get the filtered neighborlist at the end\n", + "sann.filtered_nlist\n", + "\n", + "# get the unfiltered neighborlist automatically computed by freud\n", + "sann.unfiltered_nlist" + ] + }, + { + "cell_type": "markdown", + "id": "be265984", + "metadata": {}, + "source": [ + "Custom `NeighborList`'s are also supported as inputs to the `neighbors` argument, as shown below:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "f43a6be8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import numpy as np\n", + "\n", + "# custom neighborlist\n", + "M = 15\n", + "query_point_indices = np.arange(N).repeat(M)\n", + "point_indices = np.random.randint(low=0, high=N, size=M * N)\n", + "distances = np.linalg.norm(\n", + " box.wrap(points[point_indices, :] - points[query_point_indices, :]), axis=-1\n", + ")\n", + "nlist = freud.locality.NeighborList.from_arrays(\n", + " N, N, query_point_indices, point_indices, distances\n", + ")\n", + "\n", + "# call compute\n", + "sann.compute((box, points), neighbors=nlist)\n", + "\n", + "# get the filtered neighborlist at the end\n", + "sann.filtered_nlist" + ] + }, + { + "cell_type": "markdown", + "id": "64e7fc61", + "metadata": {}, + "source": [ + "If nothing is given to the `neighbors` argument, freud will automatically compute an all pairs neighborlist (excluding ii pairs) as the unfiltered neighborlist, equivalent to what is shown below:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "e4c01834", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# make all pairs neighborlist with ii pairs excluded\n", + "nlist = freud.locality.NeighborList.all_pairs((box, points))\n", + "\n", + "# call compute\n", + "sann.compute((box, points), neighbors=nlist)\n", + "\n", + "# get the filtered neighborlist at the end\n", + "sann.filtered_nlist" + ] + }, + { + "cell_type": "markdown", + "id": "085e3365", + "metadata": {}, + "source": [ + "## The Relative Angular Distance (RAD) Method\n", + "\n", + "With the [Relative Angular Distance (RAD) method](https://aip.scitation.org/doi/10.1063/1.4961439), we label a neighbor as blocked if the angle formed between it, the other neighbor, and a blocking particle is less than a given threshold. Strictly speaking, we begin by considering neighbors of particle $i$ starting with the closest neighbor and going radially outward. We label a potential neighbor $j$ blocked by a nearer neighbor particle $k$ if\n", + "\n", + "$$ \\frac{1}{r_{ij}^2} < \\frac{1}{r_{ik}^2} \\cos(\\theta_{jik})$$\n", + "\n", + "where $r_{\\alpha \\beta}$ is the distance between particles $\\alpha$ and $\\beta$ and $\\theta_{jik}$ is the angle centered at particle $i$ extending out to particles $j$ and $k$.\n", + "\n", + "\n", + "\"Alt" + ] + }, + { + "cell_type": "markdown", + "id": "e8b3246c", + "metadata": {}, + "source": [ + "There are two variants of the RAD method: RAD$_{open}$ and RAD$_{closed}$. In RAD$_{closed}$, we consider all neighbors further away than the first blocked neighbor $j$ to also be blocked. In other words, the RAD$_{closed}$ algorithm terminates the search for neighbors of particle $i$ and begins the search for neighbors of particle $i+1$ after it finds the first blocked neighbor of particle $i$. In RAD$_{open}$, we consider all neighbors regardless of whether or not closer neighbors are blocked. The flag determining which RAD algorithm we use is controlled by the `terminate_after_blocked` argument to the `FilterRAD` constructor, as shown below:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "68401679", + "metadata": {}, + "outputs": [], + "source": [ + "# default behavior is RAD closed\n", + "rad_closed = freud.locality.FilterRAD()\n", + "\n", + "# RAD open\n", + "rad_open = freud.locality.FilterRAD(terminate_after_blocked=False)" + ] + }, + { + "cell_type": "markdown", + "id": "f1e1cc26", + "metadata": {}, + "source": [ + "All patterns in the `compute` method established in the previous section for the `FilterSANN` class also apply to the `FilterRAD` class. " + ] + }, + { + "cell_type": "markdown", + "id": "53ef28bc", + "metadata": {}, + "source": [ + "## Incomplete Shells and Performance Considerations" + ] + }, + { + "cell_type": "markdown", + "id": "26c30fa6", + "metadata": {}, + "source": [ + "There are two more issues to discuss in this tutorial related to the input unfiltered neighborlist.\n", + "\n", + "In cases where the input neighborlist is sparse (i.e. relatively few neighbors), there may not be enough neighbors in the unfiltered neighborlist such that each particle has the proper number of neighbors in the filtered neighborlist according to the SANN (RAD) algorithm. In these cases it is impossible for `FilterSANN` (`FilterRAD`) to return a full SANN (RAD) neighborlist and we say those particles have \"incomplete shells\". The `FilterSANN` (`FilterRAD`) class will by default throw an error message detailing which particles do not have full neighbor shells. An example is shown in the code block below:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "136b5bbc", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Query point indices 0, 1, 2, 3, 4 do not have full neighbor shells.\n" + ] + } + ], + "source": [ + "# sparse system\n", + "N = 5\n", + "L = 100\n", + "system = freud.data.make_random_system(L, N)\n", + "\n", + "# try to compute SANN neighborlist\n", + "sann = freud.locality.FilterSANN()\n", + "try:\n", + " sann.compute(system)\n", + "except RuntimeError as e:\n", + " print(e)" + ] + }, + { + "cell_type": "markdown", + "id": "c9eed5eb", + "metadata": {}, + "source": [ + "This error can be downgraded to a warning through an optional argument `allow_incomplete_shell` to the filter class's constructor, as shown in the cell below. When downgraded, particles with incomplete shells will have the same number of neighbors in the filtered neighborlist as they had in the unfiltered neighborlist." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "24fc48ae", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "WARNING: Query point indices 0, 1, 2, 3, 4 do not have full neighbor shells.\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sann = freud.locality.FilterSANN(allow_incomplete_shell=True)\n", + "sann.compute(system)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "2ca97268", + "metadata": {}, + "source": [ + "In cases where the input neighborlist is dense (i.e. $\\approx N$ neighbors per particle), the filtering done by `FilterSANN` and `FilterRAD` can take a long time to complete and often most of the input neighbors will be blocked. To use these classes more efficiently, physical intuition about the system can be used to limit the number of neighbors in the unfiltered neighborlist. For example, plotting a $g(r)$ for the system and choosing an `r_max` slightly larger than fist peak distance can be used as a good starting point. Another option which does not rely on any physical intuition is to start with a small number of neighbors, and slowly increase the number until the filter class does not issue a `RuntimeError`. An example demonstrating this method is shown below:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "7fa6ac52", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "number of neighbors: 4\n", + "number of neighbors: 5\n", + "number of neighbors: 6\n", + "number of neighbors: 7\n", + "number of neighbors: 8\n", + "number of neighbors: 9\n", + "number of neighbors: 10\n", + "number of neighbors: 11\n", + "number of neighbors: 12\n", + "number of neighbors: 13\n", + "number of neighbors: 14\n" + ] + } + ], + "source": [ + "# dense system\n", + "N = 1000000\n", + "L = 10\n", + "system = freud.data.make_random_system(L, N)\n", + "\n", + "# start with a small number of neighbors\n", + "num_neighbors = 4\n", + "all_shells_full = False\n", + "rad = freud.locality.FilterRAD()\n", + "\n", + "# iterate and increase num_neighbors each time\n", + "while not all_shells_full:\n", + " print(\"number of neighbors:\", num_neighbors)\n", + " try:\n", + " rad.compute(system, dict(num_neighbors=num_neighbors, exclude_ii=True))\n", + " all_shells_full = True\n", + " except RuntimeError:\n", + " num_neighbors += 1" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}