{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "*Datasets used in this example are a system of hard hexagons, simulated in the NVT thermodynamic ensemble in HOOMD-Blue, for a dense fluid (phi065) and a solid (phi075)*" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "
\n", " \n", " Loading BokehJS ...\n", "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "\n", "(function(global) {\n", " function now() {\n", " return new Date();\n", " }\n", "\n", " var force = true;\n", "\n", " if (typeof (window._bokeh_onload_callbacks) === \"undefined\" || force === true) {\n", " window._bokeh_onload_callbacks = [];\n", " window._bokeh_is_loading = undefined;\n", " }\n", "\n", "\n", " \n", " if (typeof (window._bokeh_timeout) === \"undefined\" || force === true) {\n", " window._bokeh_timeout = Date.now() + 5000;\n", " window._bokeh_failed_load = false;\n", " }\n", "\n", " var NB_LOAD_WARNING = {'data': {'text/html':\n", " \"
\\n\"+\n", " \"

\\n\"+\n", " \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n", " \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n", " \"

\\n\"+\n", " \"\\n\"+\n", " \"\\n\"+\n", " \"from bokeh.resources import INLINE\\n\"+\n", " \"output_notebook(resources=INLINE)\\n\"+\n", " \"\\n\"+\n", " \"
\"}};\n", "\n", " function display_loaded() {\n", " if (window.Bokeh !== undefined) {\n", " document.getElementById(\"2efe779a-d382-41d7-8033-02826955a85f\").textContent = \"BokehJS successfully loaded.\";\n", " } else if (Date.now() < window._bokeh_timeout) {\n", " setTimeout(display_loaded, 100)\n", " }\n", " }\n", "\n", " function run_callbacks() {\n", " window._bokeh_onload_callbacks.forEach(function(callback) { callback() });\n", " delete window._bokeh_onload_callbacks\n", " console.info(\"Bokeh: all callbacks have finished\");\n", " }\n", "\n", " function load_libs(js_urls, callback) {\n", " window._bokeh_onload_callbacks.push(callback);\n", " if (window._bokeh_is_loading > 0) {\n", " console.log(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n", " return null;\n", " }\n", " if (js_urls == null || js_urls.length === 0) {\n", " run_callbacks();\n", " return null;\n", " }\n", " console.log(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n", " window._bokeh_is_loading = js_urls.length;\n", " for (var i = 0; i < js_urls.length; i++) {\n", " var url = js_urls[i];\n", " var s = document.createElement('script');\n", " s.src = url;\n", " s.async = false;\n", " s.onreadystatechange = s.onload = function() {\n", " window._bokeh_is_loading--;\n", " if (window._bokeh_is_loading === 0) {\n", " console.log(\"Bokeh: all BokehJS libraries loaded\");\n", " run_callbacks()\n", " }\n", " };\n", " s.onerror = function() {\n", " console.warn(\"failed to load library \" + url);\n", " };\n", " console.log(\"Bokeh: injecting script tag for BokehJS library: \", url);\n", " document.getElementsByTagName(\"head\")[0].appendChild(s);\n", " }\n", " };var element = document.getElementById(\"2efe779a-d382-41d7-8033-02826955a85f\");\n", " if (element == null) {\n", " console.log(\"Bokeh: ERROR: autoload.js configured with elementid '2efe779a-d382-41d7-8033-02826955a85f' but no matching script tag was found. \")\n", " return false;\n", " }\n", "\n", " var js_urls = [\"https://cdn.pydata.org/bokeh/release/bokeh-0.12.4.min.js\", \"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.4.min.js\"];\n", "\n", " var inline_js = [\n", " function(Bokeh) {\n", " Bokeh.set_log_level(\"info\");\n", " },\n", " \n", " function(Bokeh) {\n", " \n", " document.getElementById(\"2efe779a-d382-41d7-8033-02826955a85f\").textContent = \"BokehJS is loading...\";\n", " },\n", " function(Bokeh) {\n", " console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-0.12.4.min.css\");\n", " Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-0.12.4.min.css\");\n", " console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.4.min.css\");\n", " Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.4.min.css\");\n", " }\n", " ];\n", "\n", " function run_inline_js() {\n", " \n", " if ((window.Bokeh !== undefined) || (force === true)) {\n", " for (var i = 0; i < inline_js.length; i++) {\n", " inline_js[i](window.Bokeh);\n", " }if (force === true) {\n", " display_loaded();\n", " }} else if (Date.now() < window._bokeh_timeout) {\n", " setTimeout(run_inline_js, 100);\n", " } else if (!window._bokeh_failed_load) {\n", " console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n", " window._bokeh_failed_load = true;\n", " } else if (force !== true) {\n", " var cell = $(document.getElementById(\"2efe779a-d382-41d7-8033-02826955a85f\")).parents('.cell').data().cell;\n", " cell.output_area.append_execute_result(NB_LOAD_WARNING)\n", " }\n", "\n", " }\n", "\n", " if (window._bokeh_is_loading === 0) {\n", " console.log(\"Bokeh: BokehJS loaded, going straight to plotting\");\n", " run_inline_js();\n", " } else {\n", " load_libs(js_urls, function() {\n", " console.log(\"Bokeh: BokehJS plotting callback run at\", now());\n", " run_inline_js();\n", " });\n", " }\n", "}(this));" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from bokeh.io import output_notebook\n", "output_notebook()\n", "from bokeh.models import Legend\n", "from bokeh.plotting import figure, output_file, show\n", "from bokeh.layouts import gridplot\n", "import numpy as np\n", "import time\n", "from freud import parallel\n", "parallel.setNumThreads(4)\n", "\n", "# define vertices for hexagons\n", "verts = [[0.537284965911771, 0.31020161970069976],\n", " [3.7988742065678664e-17, 0.6204032394013997],\n", " [-0.5372849659117709, 0.31020161970070004],\n", " [-0.5372849659117711, -0.31020161970069976],\n", " [-1.1396622619703597e-16, -0.6204032394013997],\n", " [0.5372849659117711, -0.3102016197006997]]\n", "verts = np.array(verts)\n", "\n", "# define colors for our system\n", "c_list = [\"#30A2DA\", \"#FC4F30\", \"#E5AE38\", \"#6D904F\", \"#9757DB\",\n", " \"#188487\", \"#FF7F00\", \"#9A2C66\", \"#626DDA\", \"#8B8B8B\"]\n", "c_dict = dict()\n", "c_dict[6] = c_list[0]\n", "c_dict[5] = c_list[1]\n", "c_dict[4] = c_list[2]\n", "c_dict[3] = c_list[7]\n", "c_dict[7] = c_list[4]\n", "\n", "def default_bokeh(p):\n", " \"\"\"\n", " wrapper which takes the default bokeh outputs and changes them to more sensible values\n", " \"\"\"\n", " p.title.text_font_size = \"18pt\"\n", " p.title.align = \"center\"\n", "\n", " p.xaxis.axis_label_text_font_size = \"14pt\"\n", " p.yaxis.axis_label_text_font_size = \"14pt\"\n", "\n", " p.xaxis.major_tick_in = 10\n", " p.xaxis.major_tick_out = 0\n", " p.xaxis.minor_tick_in = 5\n", " p.xaxis.minor_tick_out = 0\n", "\n", " p.yaxis.major_tick_in = 10\n", " p.yaxis.major_tick_out = 0\n", " p.yaxis.minor_tick_in = 5\n", " p.yaxis.minor_tick_out = 0\n", "\n", " p.xaxis.major_label_text_font_size = \"12pt\"\n", " p.yaxis.major_label_text_font_size = \"12pt\"\n", "\n", "def cubeellipse(theta, lam=0.5, gamma=1., s=4.0, r=1., h=1.):\n", " \"\"\"Create an RGB colormap from an input angle theta. Takes lam (a list of\n", " intensity values, from 0 to 1), gamma (a nonlinear weighting power),\n", " s (starting angle), r (number of revolutions around the circle), and\n", " h (a hue factor).\"\"\"\n", " import numpy\n", " lam = lam**gamma\n", "\n", " a = h*lam*(1 - lam)*.5\n", " v = numpy.array([[-.14861, 1.78277], [-.29227, -.90649], [1.97294, 0.]], dtype=numpy.float32)\n", " ctarray = numpy.array([numpy.cos(theta*r + s), numpy.sin(theta*r + s)], dtype=numpy.float32)\n", " # convert to 255 rgb\n", " ctarray = (lam + a*v.dot(ctarray)).T\n", " ctarray *= 255\n", " ctarray = ctarray.astype(dtype=np.int32)\n", " return ctarray\n", "\n", "def local_to_global(verts, positions, orientations):\n", " \"\"\"\n", " Take a list of vertices, positions, and orientations and create\n", " a list of vertices in the \"global coordinate system\" for plotting\n", " in bokeh\n", " \"\"\"\n", " num_particles = len(positions)\n", " num_verts = len(verts)\n", " # create list of vertices in the \"local reference frame\" i.e.\n", " # centered at (0,0)\n", " l_verts = np.zeros(shape=(num_particles, num_verts, 2), dtype=np.float32)\n", " l_verts[:] = verts\n", " # create array of rotation matrices\n", " rot_mat = np.zeros(shape=(num_particles, 2, 2), dtype=np.float32)\n", " for i, theta in enumerate(orientations):\n", " rot_mat[i] = [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]\n", " # rotate; uses einsum for speed; please see numpy documentation\n", " # for more information\n", " r_verts = np.einsum(\"lij,lkj->lki\", rot_mat, l_verts)\n", " # now translate to global coordinates\n", " # need to create a position array with same shape as vertex array\n", " l_pos = np.zeros(shape=(num_particles, num_verts, 2), dtype=np.float32)\n", " for i in range(num_particles):\n", " for j in range(len(verts)):\n", " l_pos[i,j] = positions[i]\n", " # translate\n", " output_array = np.add(r_verts, l_pos)\n", " return output_array\n", "\n", "def clamp(x):\n", " \"\"\"\n", " limit values between 0 and 255\n", " http://stackoverflow.com/questions/3380726/converting-a-rgb-color-tuple-to-a-six-digit-code-in-python\n", " \"\"\"\n", " return max(0, min(x, 255))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Nearest Neighbors\n", "\n", "One of the basic computations required for higher-level computations (such as the [hexatic order parameter](AIChE Hexatic.ipynb)) is finding the nearest neighbors of a particle. This tutorial will show you how to compute the nearest neighbors and visualize that data.\n", "\n", "The algorithm is straightforward:\n", "\n", "~~~\n", "for each particle i:\n", " for each particle j in neighbor_cells(i):\n", " r_ij = position[j] - position[i]\n", " r = sqrt(dot(r_ij, r_ij))\n", " l_r_array.append(r)\n", " l_n_array.append(j)\n", " # sort by distance\n", " sort(n_array, r_array)\n", " neighbor_array[i] = n_array[:k]\n", "~~~" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# import the box and locality modules;\n", "# Nearest Neighbors calculation is provided by the locality module\n", "from freud import box, locality\n", "# load the data\n", "data_path = \"ex_data/phi065\"\n", "box_data = np.load(\"{}/box_data.npy\".format(data_path))\n", "pos_data = np.load(\"{}/pos_data.npy\".format(data_path))\n", "quat_data = np.load(\"{}/quat_data.npy\".format(data_path))\n", "n_frames = pos_data.shape[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Viewing our system\n", "\n", "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" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# grab data from last frame\n", "l_box = box_data[-1]\n", "l_pos = pos_data[-1]\n", "l_quat = quat_data[-1]\n", "l_ang = 2*np.arctan2(np.copy(l_quat[:,3]), np.copy(l_quat[:,0]))\n", "\n", "# create box\n", "fbox = box.Box(Lx=l_box[\"Lx\"], Ly=l_box[\"Ly\"], is2D=True)\n", "side_length = max(fbox.Lx, fbox.Ly)\n", "l_min = -side_length / 2.0\n", "l_min *= 1.1\n", "l_max = -l_min\n", "\n", "# take local vertices and rotate, translate into system coordinates\n", "patches = local_to_global(verts, l_pos[:,0:2], l_ang)\n", "\n", "# plot\n", "p = figure(title=\"System Visualization\",x_range=(l_min, l_max), y_range=(l_min, l_max))\n", "p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),\n", " fill_color=(42,126,187), line_color=\"black\", line_width=1.5, legend=\"hexagons\")\n", "# box display\n", "p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],\n", " ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],\n", " fill_color=(0,0,0,0), line_color=\"black\", line_width=2)\n", "p.legend.location='bottom_center'\n", "p.legend.orientation='horizontal'\n", "default_bokeh(p)\n", "show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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\n", "\n", "## The Nearest Neighbor object\n", "\n", "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:\n", "\n", "* `strict_cut=False` (*default*): The value for `rmax` is expanded until every particle has \\\\(k\\\\) nearest neighbors\n", "* `strict_cut=True`: the `rmax` value is not expanded, so that any \"vacancies\" in the number of neighbors found are filled with `UINTMAX`\n", "\n", "### `strict_cut=False`\n", "\n", "First we show how to use the `strict_cut=False` mode to find the neighbors of a specific particle" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time to calc 1 frame = 0.023049116134643555\n" ] }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# create freud nearest neighbor object\n", "# set number of neighbors\n", "n_neigh = 6\n", "# create freud nearest neighbors object\n", "nn = locality.NearestNeighbors(rmax=1.5,n_neigh=n_neigh,strict_cut=False)\n", "\n", "# compute nearest neighbors for 6 nearest neighbors\n", "start_time = time.time()\n", "nn.compute(fbox, l_pos, l_pos)\n", "stop_time = time.time()\n", "print(\"time to calc 1 frame = {}\".format(stop_time-start_time))\n", "# get the neighborlist\n", "n_list = nn.getNeighborList()\n", "# get the number of particles\n", "num_particles = nn.getNRef()\n", "# get the neighbors for particle 0\n", "pidx = 0\n", "n_idxs = n_list[pidx]\n", "\n", "# get position, orientation for the central particle\n", "center_pos = np.zeros(shape=(1, 3),dtype=np.float32)\n", "center_ang = np.zeros(shape=(1),dtype=np.float32)\n", "center_pos[:] = l_pos[pidx]\n", "center_ang[:] = l_ang[pidx]\n", "\n", "# get the positions, orientations for the neighbor particles\n", "neigh_pos = np.zeros(shape=(n_neigh, 3),dtype=np.float32)\n", "neigh_ang = np.zeros(shape=(n_neigh),dtype=np.float32)\n", "neigh_pos[:] = l_pos[n_idxs]\n", "neigh_ang[:] = l_ang[n_idxs]\n", "\n", "# render in bokeh\n", "# create array of transformed positions\n", "c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)\n", "n_patches = local_to_global(verts, neigh_pos[:,0:2], neigh_ang)\n", "# turn into list of colors\n", "# bokeh (as of this version) requires hex colors, so convert rgb to hex\n", "center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])\n", "neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])\n", "\n", "# plot\n", "p = figure(title=\"Nearest Neighbors visualization\",x_range=(l_min, l_max), y_range=(l_min, l_max))\n", "p.patches(xs=n_patches[:,:,0].tolist(), ys=n_patches[:,:,1].tolist(),\n", " fill_color=neigh_color.tolist(), line_color=\"black\", legend=\"neighbors\")\n", "p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),\n", " fill_color=center_color.tolist(), line_color=\"black\", legend=\"centers\")\n", "# box display\n", "p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],\n", " ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],\n", " fill_color=(0,0,0,0), line_color=\"black\", line_width=2)\n", "p.legend.location='bottom_center'\n", "p.legend.orientation='horizontal'\n", "default_bokeh(p)\n", "show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that nearest neighbors properly handles periodic boundary conditions\n", "\n", "We do the same thing below, but for a particle/neighbors not spanning the box" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time to calc 1 frame = 0.02012920379638672\n" ] }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# create freud nearest neighbor object\n", "# set number of neighbors\n", "n_neigh = 6\n", "# create freud nearest neighbors object\n", "nn = locality.NearestNeighbors(rmax=1.5,n_neigh=n_neigh,strict_cut=False)\n", "\n", "# compute nearest neighbors for 6 nearest neighbors\n", "start_time = time.time()\n", "nn.compute(fbox, l_pos, l_pos)\n", "stop_time = time.time()\n", "print(\"time to calc 1 frame = {}\".format(stop_time-start_time))\n", "# get the neighborlist\n", "n_list = nn.getNeighborList()\n", "# get the number of particles\n", "num_particles = nn.getNRef()\n", "# get the neighbors for particle 1000\n", "pidx = 1000\n", "n_idxs = n_list[pidx]\n", "\n", "# get position, orientation for the central particle\n", "center_pos = np.zeros(shape=(1, 3),dtype=np.float32)\n", "center_ang = np.zeros(shape=(1),dtype=np.float32)\n", "center_pos[:] = l_pos[pidx]\n", "center_ang[:] = l_ang[pidx]\n", "\n", "# get positions, orientations for the neighbor particles\n", "neigh_pos = np.zeros(shape=(n_neigh, 3),dtype=np.float32)\n", "neigh_ang = np.zeros(shape=(n_neigh),dtype=np.float32)\n", "neigh_pos[:] = l_pos[n_idxs]\n", "neigh_ang[:] = l_ang[n_idxs]\n", "\n", "non_neigh_pos = np.zeros(shape=(1, 3),dtype=np.float32)\n", "non_neigh_ang = np.zeros(shape=(1),dtype=np.float32)\n", "non_neigh_pos[:] = neigh_pos[-1]\n", "non_neigh_ang[:] = neigh_ang[-1]\n", "\n", "# render in bokeh\n", "# create array of transformed positions\n", "c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)\n", "n_patches = local_to_global(verts, neigh_pos[:,0:2], neigh_ang)\n", "non_n_patches = local_to_global(verts, non_neigh_pos[:,0:2], non_neigh_ang)\n", "# turn into list of colors\n", "# bokeh (as of this version) requires hex colors, so convert rgb to hex\n", "center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])\n", "neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])\n", "non_neigh_color = np.array([c_list[-2] for _ in range(non_neigh_pos.shape[0])])\n", "# color the last particle differently\n", "\n", "# plot\n", "p = figure(title=\"Nearest Neighbors visualization\",x_range=(l_min, l_max), y_range=(l_min, l_max))\n", "p.patches(xs=n_patches[:,:,0].tolist(), ys=n_patches[:,:,1].tolist(),\n", " fill_color=neigh_color.tolist(), line_color=\"black\", legend=\"neighbors\")\n", "p.patches(xs=non_n_patches[:,:,0].tolist(), ys=non_n_patches[:,:,1].tolist(),\n", " fill_color=non_neigh_color.tolist(), line_color=\"black\", legend=\"non-neighbor\")\n", "p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),\n", " fill_color=center_color.tolist(), line_color=\"black\", legend=\"centers\")\n", "# box display\n", "p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],\n", " ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],\n", " fill_color=(0,0,0,0), line_color=\"black\", line_width=2)\n", "p.legend.location='bottom_center'\n", "p.legend.orientation='horizontal'\n", "default_bokeh(p)\n", "show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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?\n", "\n", "### `strict_cut=True`\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time to calc 1 frame = 0.004876852035522461\n" ] }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# create freud nearest neighbor object\n", "# set number of neighbors\n", "n_neigh = 8\n", "# create freud nearest neighbors object\n", "# set rmax to some value\n", "rmax=1.7\n", "nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)\n", "\n", "# compute nearest neighbors for 6 nearest neighbors\n", "start_time = time.time()\n", "nn.compute(fbox, l_pos, l_pos)\n", "stop_time = time.time()\n", "print(\"time to calc 1 frame = {}\".format(stop_time-start_time))\n", "# get the neighborlist\n", "n_list = np.copy(nn.getNeighborList())\n", "# get the number of particles\n", "num_particles = nn.getNRef()\n", "# now for array manipulation magic\n", "# create an integer array of the same shape as the neighbor list array\n", "int_arr = np.ones(shape=n_list.shape, dtype=np.int32)\n", "# \"search\" for non-indexed particles (missing neighbors)\n", "# while it would be most accurate to use the UINTMAX value\n", "# provided by nn.getUINTMAX(), but this works just as well\n", "int_arr[n_list > (num_particles-1)] = 0\n", "# sum along particle index axis to determine the number of neighbors per particle\n", "n_neighbors = np.sum(int_arr, axis=1)\n", "# find the complement (if desired) to find number of missing neighbors per particle\n", "n_deficits = n_neigh - n_neighbors\n", "\n", "p = figure(title=\"Nearest Neighbors visualization\",x_range=(l_min, l_max), y_range=(l_min, l_max))\n", "for k in np.unique(n_neighbors):\n", " # find particles with k neighbors\n", " c_idxs = np.copy(np.where(n_neighbors==k)[0])\n", " center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)\n", " center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)\n", " center_pos = l_pos[c_idxs]\n", " center_ang = l_ang[c_idxs]\n", " c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)\n", " center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])\n", " p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),\n", " fill_color=center_color.tolist(), line_color=\"black\", legend=\"k={}\".format(k))\n", "p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],\n", " ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],\n", " fill_color=(0,0,0,0), line_color=\"black\", line_width=2)\n", "p.legend.location='bottom_center'\n", "p.legend.orientation='horizontal'\n", "default_bokeh(p)\n", "show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualize each set of k values independently" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time to calc 1 frame = 0.0025110244750976562\n" ] }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# create freud nearest neighbor object\n", "# set number of neighbors\n", "n_neigh = 8\n", "# create freud nearest neighbors object\n", "# set rmax to some value\n", "rmax=1.7\n", "nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)\n", "\n", "# compute nearest neighbors for 6 nearest neighbors\n", "start_time = time.time()\n", "nn.compute(fbox, l_pos, l_pos)\n", "stop_time = time.time()\n", "print(\"time to calc 1 frame = {}\".format(stop_time-start_time))\n", "# get the neighborlist\n", "n_list = np.copy(nn.getNeighborList())\n", "# get the number of particles\n", "num_particles = nn.getNRef()\n", "# now for array manipulation magic\n", "# create an integer array of the same shape as the neighbor list array\n", "int_arr = np.ones(shape=n_list.shape, dtype=np.int32)\n", "# \"search\" for non-indexed particles (missing neighbors)\n", "# while it would be most accurate to use the UINTMAX value\n", "# provided by nn.getUINTMAX(), but this works just as well\n", "int_arr[n_list > (num_particles-1)] = 0\n", "# sum along particle index axis to determine the number of neighbors per particle\n", "n_neighbors = np.sum(int_arr, axis=1)\n", "# find the complement (if desired) to find number of missing neighbors per particle\n", "n_deficits = n_neigh - n_neighbors\n", "\n", "for k in np.unique(n_neighbors):\n", " p = figure(title=\"Nearest Neighbors: k={}\".format(k),\n", " x_range=(l_min, l_max), y_range=(l_min, l_max))\n", " # find particles with k neighbors\n", " c_idxs = np.copy(np.where(n_neighbors==k)[0])\n", " center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)\n", " center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)\n", " center_pos = l_pos[c_idxs]\n", " center_ang = l_ang[c_idxs]\n", " c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)\n", " # take local vertices and rotate, translate into system coordinates\n", " patches = local_to_global(verts, l_pos[:,0:2], l_ang)\n", " center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])\n", " p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),\n", " fill_color=(128,128,128,0.1), line_color=(0,0,0,0.1), legend=\"other\")\n", " p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),\n", " fill_color=center_color.tolist(), line_color=\"black\", legend=\"k={}\".format(k))\n", " p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],\n", " ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],\n", " fill_color=(0,0,0,0), line_color=\"black\", line_width=2)\n", " p.legend.location='bottom_center'\n", " p.legend.orientation='horizontal'\n", " default_bokeh(p)\n", " show(p)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time to calc 1 frame = 0.002343893051147461\n" ] }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# create freud nearest neighbor object\n", "# set number of neighbors\n", "n_neigh = 8\n", "# create freud nearest neighbors object\n", "# set rmax to some value\n", "rmax=1.7\n", "nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)\n", "\n", "# compute nearest neighbors for 6 nearest neighbors\n", "start_time = time.time()\n", "nn.compute(fbox, l_pos, l_pos)\n", "stop_time = time.time()\n", "print(\"time to calc 1 frame = {}\".format(stop_time-start_time))\n", "# get the neighborlist\n", "n_list = np.copy(nn.getNeighborList())\n", "# get the number of particles\n", "num_particles = nn.getNRef()\n", "# now for array manipulation magic\n", "# create an integer array of the same shape as the neighbor list array\n", "int_arr = np.ones(shape=n_list.shape, dtype=np.int32)\n", "# \"search\" for non-indexed particles (missing neighbors)\n", "# while it would be most accurate to use the UINTMAX value\n", "# provided by nn.getUINTMAX(), but this works just as well\n", "int_arr[n_list > (num_particles-1)] = 0\n", "# sum along particle index axis to determine the number of neighbors per particle\n", "n_neighbors = np.sum(int_arr, axis=1)\n", "# find the complement (if desired) to find number of missing neighbors per particle\n", "n_deficits = n_neigh - n_neighbors\n", "\n", "for k in np.unique(n_neighbors):\n", " p = figure(title=\"Nearest Neighbors: k={}\".format(k),\n", " x_range=(l_min, l_max), y_range=(l_min, l_max))\n", " # find particles with k neighbors\n", " c_idxs = np.copy(np.where(n_neighbors==k)[0])\n", " center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)\n", " center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)\n", " center_pos = l_pos[c_idxs]\n", " center_ang = l_ang[c_idxs]\n", "\n", " neigh_pos = np.zeros(shape=(k*len(c_idxs), 3),dtype=np.float32)\n", " neigh_ang = np.zeros(shape=(k*len(c_idxs)),dtype=np.float32)\n", " for i, pidx in enumerate(c_idxs):\n", " # create a list of positions, angles to draw\n", " n_idxs = n_list[pidx,:k]\n", " for j, nidx in enumerate(n_idxs):\n", " neigh_pos[k*i+j] = l_pos[nidx]\n", " neigh_ang[k*i+j] = l_ang[nidx]\n", " c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)\n", " n_patches = local_to_global(verts, neigh_pos[:,0:2], neigh_ang)\n", " patches = local_to_global(verts, l_pos[:,0:2], l_ang)\n", " center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])\n", " neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])\n", " p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),\n", " fill_color=(128,128,128,0.1), line_color=(0,0,0,0.1), legend=\"other\")\n", " p.patches(xs=n_patches[:,:,0].tolist(), ys=n_patches[:,:,1].tolist(),\n", " fill_color=neigh_color.tolist(), line_color=\"black\", legend=\"neighbors\")\n", " p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),\n", " fill_color=center_color.tolist(), line_color=\"black\", legend=\"k={}\".format(k))\n", " p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],\n", " ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],\n", " fill_color=(0,0,0,0), line_color=\"black\", line_width=2)\n", " p.legend.location='bottom_center'\n", " p.legend.orientation='horizontal'\n", " default_bokeh(p)\n", " show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " ## Compare against a solid/crystalline system\n", " \n", " Visualize in the same way, but with a solid system" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time to calc 1 frame = 0.0023381710052490234\n" ] }, { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data_path = \"ex_data/phi075\"\n", "box_data = np.load(\"{}/box_data.npy\".format(data_path))\n", "pos_data = np.load(\"{}/pos_data.npy\".format(data_path))\n", "quat_data = np.load(\"{}/quat_data.npy\".format(data_path))\n", "n_frames = pos_data.shape[0]\n", "\n", "# grab data from last frame\n", "l_box = box_data[-1]\n", "l_pos = pos_data[-1]\n", "l_quat = quat_data[-1]\n", "l_ang = 2*np.arctan2(np.copy(l_quat[:,3]), np.copy(l_quat[:,0]))\n", "\n", "# create box\n", "fbox = box.Box(Lx=l_box[\"Lx\"], Ly=l_box[\"Ly\"], is2D=True)\n", "side_length = max(fbox.Lx, fbox.Ly)\n", "l_min = -side_length / 2.0\n", "l_min *= 1.1\n", "l_max = -l_min\n", "# create freud nearest neighbor object\n", "# set number of neighbors\n", "n_neigh = 8\n", "# create freud nearest neighbors object\n", "# set rmax to some value\n", "rmax=1.65\n", "nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)\n", "\n", "# compute nearest neighbors for 6 nearest neighbors\n", "start_time = time.time()\n", "nn.compute(fbox, l_pos, l_pos)\n", "stop_time = time.time()\n", "print(\"time to calc 1 frame = {}\".format(stop_time-start_time))\n", "# get the neighborlist\n", "n_list = np.copy(nn.getNeighborList())\n", "# get the number of particles\n", "num_particles = nn.getNRef()\n", "# now for array manipulation magic\n", "# create an integer array of the same shape as the neighbor list array\n", "int_arr = np.ones(shape=n_list.shape, dtype=np.int32)\n", "# \"search\" for non-indexed particles (missing neighbors)\n", "# while it would be most accurate to use the UINTMAX value\n", "# provided by nn.getUINTMAX(), but this works just as well\n", "int_arr[n_list > (num_particles-1)] = 0\n", "# sum along particle index axis to determine the number of neighbors per particle\n", "n_neighbors = np.sum(int_arr, axis=1)\n", "# find the complement (if desired) to find number of missing neighbors per particle\n", "n_deficits = n_neigh - n_neighbors\n", "\n", "p = figure(title=\"Nearest Neighbors visualization\",x_range=(l_min, l_max), y_range=(l_min, l_max))\n", "for k in np.unique(n_neighbors):\n", " # find particles with k neighbors\n", " c_idxs = np.copy(np.where(n_neighbors==k)[0])\n", " center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)\n", " center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)\n", " center_pos = l_pos[c_idxs]\n", " center_ang = l_ang[c_idxs]\n", " c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)\n", " center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])\n", " p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),\n", " fill_color=center_color.tolist(), line_color=\"black\", legend=\"k={}\".format(k))\n", "p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],\n", " ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],\n", " fill_color=(0,0,0,0), line_color=\"black\", line_width=2)\n", "p.legend.location='bottom_center'\n", "p.legend.orientation='horizontal'\n", "default_bokeh(p)\n", "show(p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.5.3" } }, "nbformat": 4, "nbformat_minor": 0 }