{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction\n", "\n", "In this example we will be computing the *radial distribution function* TODO add link to documentation...when that documentation is online...get that added to the website!\n", "\n", "To calculate the radial distribution function TODO (add in math to actual section?), we will need to follow the following workflow:\n", "\n", "1. Load data\n", "2. Create Freud RDF object\n", "3. Compute the RDF\n", "4. Extract data from the Freud RDF object\n", "5. (optional) Plot data\n", "\n", "# Required Packages\n", "\n", "We will use the following python packages, please install with your favorite package manager:\n", "\n", "* `numpy` - general scientific goodness\n", "* `bokeh` - plotting\n", "* `matplotlib` - plotting" ] }, { "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 = \"1\";\n", "\n", " if (typeof (window._bokeh_onload_callbacks) === \"undefined\" || force !== \"\") {\n", " window._bokeh_onload_callbacks = [];\n", " window._bokeh_is_loading = undefined;\n", " }\n", "\n", "\n", " \n", " if (typeof (window._bokeh_timeout) === \"undefined\" || force !== \"\") {\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", " Bokeh.$(\"#3a89b414-a2e0-4f5c-8398-a440a3ea11fc\").text(\"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(\"3a89b414-a2e0-4f5c-8398-a440a3ea11fc\");\n", " if (element == null) {\n", " console.log(\"Bokeh: ERROR: autoload.js configured with elementid '3a89b414-a2e0-4f5c-8398-a440a3ea11fc' 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.3.min.js', 'https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.3.min.js'];\n", "\n", " var inline_js = [\n", " function(Bokeh) {\n", " Bokeh.set_log_level(\"info\");\n", " },\n", " \n", " function(Bokeh) {\n", " \n", " Bokeh.$(\"#3a89b414-a2e0-4f5c-8398-a440a3ea11fc\").text(\"BokehJS is loading...\");\n", " },\n", " function(Bokeh) {\n", " console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-0.12.3.min.css\");\n", " Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-0.12.3.min.css\");\n", " console.log(\"Bokeh: injecting CSS: https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.3.min.css\");\n", " Bokeh.embed.inject_css(\"https://cdn.pydata.org/bokeh/release/bokeh-widgets-0.12.3.min.css\");\n", " }\n", " ];\n", "\n", " function run_inline_js() {\n", " \n", " if ((window.Bokeh !== undefined) || (force === \"1\")) {\n", " for (var i = 0; i < inline_js.length; i++) {\n", " inline_js[i](window.Bokeh);\n", " }if (force === \"1\") {\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) {\n", " var cell = $(\"#3a89b414-a2e0-4f5c-8398-a440a3ea11fc\").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.plotting import figure, output_file, show\n", "from bokeh.layouts import gridplot\n", "import numpy as np\n", "from freud import parallel, box, density\n", "parallel.setNumThreads(4)\n", "\n", "def default_bokeh(p):\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\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Set Up Freud\n", "\n", "Now we need to import freud to our environment. In this case, we will import specific modules rather than import the whole of freud. Feel free to just use\n", "\n", " import freud\n", "\n", "as you see fit." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from freud import parallel, box, density\n", "parallel.setNumThreads(4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We imported the following modules:\n", "\n", "* parallel - allows you to set the number of threads you will use during calculations. Yes, Freud is parallelized!\n", " - By default freud will use the maximum number of threads available on your machine. Here, I set it to 4, please set as needed on your machine. If working on a shared machine, consider using a moderate number of threads.\n", " - Don't forget to benchmark your calculations, or refer to the freud benchmarks to determine the optimal number of cores to use. More is not always better!\n", "* box - module used to create the simulation boxes used in freud calculations\n", "* density - module containing density-related calculations, including the radial distribution function" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Loading Data\n", "\n", "Freud makes no assumptions about your data, and doesn't provide a method to load specified formats, everything passed into freud must be a `np.ndarray` of type required by freud (see documentation). In this example numpy binary files containing data are used.\n", "\n", "# Create the RDF Object\n", "\n", "Each module in Freud contains a number of available computations. In order to perform any computations, you need to create an instance of the object. Please refer to the documentation for what arguments need to be supplied for the compute module you wish to use.\n", "\n", "Let's create the RDF object. For your convenience, we use the help function to show what needs to be supplied:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on class RDF in module freud._freud:\n", "\n", "class RDF(builtins.object)\n", " | Computes RDF for supplied data\n", " | \n", " | The RDF (:math:`g \\left( r \\right)`) is computed and averaged for a given set of reference points in a sea of data points. Providing the same points calculates them against themselves. Computing the RDF results in an rdf array listing the value of the RDF at each given :math:`r`, listed in the r array.\n", " | \n", " | The values of :math:`r` to compute the rdf are set by the values of rmax, dr in the constructor. rmax sets the maximum\n", " | distance at which to calculate the :math:`g \\left( r \\right)` while dr determines the step size for each bin.\n", " | \n", " | .. moduleauthor:: Eric Harper \n", " | \n", " | .. note::\n", " | 2D: RDF properly handles 2D boxes. Requires the points to be passed in [x, y, 0]. Failing to z=0 will lead to undefined behavior.\n", " | \n", " | :param rmax: maximum distance to calculate\n", " | :param dr: distance between histogram bins\n", " | :type rmax: float\n", " | :type dr: float\n", " | \n", " | Methods defined here:\n", " | \n", " | __new__(*args, **kwargs) from builtins.type\n", " | Create and return a new object. See help(type) for accurate signature.\n", " | \n", " | accumulate(...)\n", " | RDF.accumulate(self, box, ndarray ref_points, ndarray points)\n", " | \n", " | Calculates the rdf and adds to the current rdf histogram.\n", " | \n", " | :param box: simulation box\n", " | :param ref_points: reference points to calculate the local density\n", " | :param points: points to calculate the local density\n", " | :type box: :py:class:`freud.box.Box`\n", " | :type ref_points: :class:`numpy.ndarray`, shape=(:math:`N_{particles}`, 3), dtype= :class:`numpy.float32`\n", " | :type points: :class:`numpy.ndarray`, shape=(:math:`N_{particles}`, 3), dtype= :class:`numpy.float32`\n", " | \n", " | compute(...)\n", " | RDF.compute(self, box, ndarray ref_points, ndarray points)\n", " | \n", " | Calculates the rdf for the specified points. Will overwrite the current histogram.\n", " | \n", " | :param box: simulation box\n", " | :param ref_points: reference points to calculate the local density\n", " | :param points: points to calculate the local density\n", " | :type box: :py:meth:`freud.box.Box`\n", " | :type ref_points: :class:`numpy.ndarray`, shape=(:math:`N_{particles}`, 3), dtype= :class:`numpy.float32`\n", " | :type points: :class:`numpy.ndarray`, shape=(:math:`N_{particles}`, 3), dtype= :class:`numpy.float32`\n", " | \n", " | getBox(...)\n", " | RDF.getBox(self)\n", " | \n", " | :return: Freud Box\n", " | :rtype: :py:class:`freud.box.Box`\n", " | \n", " | getNr(...)\n", " | RDF.getNr(self)\n", " | \n", " | :return: histogram of cumulative rdf values\n", " | :rtype: :class:`numpy.ndarray`, shape=(:math:`N_{bins}`, 3), dtype= :class:`numpy.float32`\n", " | \n", " | getR(...)\n", " | RDF.getR(self)\n", " | \n", " | :return: values of the histogram bin centers\n", " | :rtype: :class:`numpy.ndarray`, shape=(:math:`N_{bins}`, 3), dtype= :class:`numpy.float32`\n", " | \n", " | getRDF(...)\n", " | RDF.getRDF(self)\n", " | \n", " | :return: histogram of rdf values\n", " | :rtype: :class:`numpy.ndarray`, shape=(:math:`N_{bins}`, 3), dtype= :class:`numpy.float32`\n", " | \n", " | reduceRDF(...)\n", " | RDF.reduceRDF(self)\n", " | \n", " | Reduces the histogram in the values over N processors to a single histogram. This is called automatically by\n", " | :py:meth:`freud.density.RDF.getRDF()`, :py:meth:`freud.density.RDF.getNr()`.\n", " | \n", " | resetRDF(...)\n", " | RDF.resetRDF(self)\n", " | \n", " | resets the values of RDF in memory\n", "\n" ] } ], "source": [ "help(density.RDF)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like the constructor takes in `rmax` and `dr` as parameters, or the maximum distance to calculate the rdf, and the size of the rdf histogram bin e.g. `r_max = 10` and `dr = 0.1`:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "rdf = density.RDF(rmax=10.0, dr=0.1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Compute the RDF\n", "\n", "It's time to actually compute the rdf, so here we go, what is happening is briefly explained in the comments" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "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", "n_frames = pos_data.shape[0]\n", "\n", "# compute the rdf for the last frame\n", "# read box, position data\n", "l_box = box_data[-1]\n", "l_pos = pos_data[-1]\n", "# create the freud box object\n", "fbox = box.Box(Lx=l_box[\"Lx\"], Ly=l_box[\"Ly\"], is2D=True)\n", "# compute\n", "rdf.compute(fbox, l_pos, l_pos)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# I computed the RDF, now what?\n", "\n", "You might be thinking to yourself \"ok, now what?\" First, you need to get your data out of freud. Second, you need to visualize that data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Getting your data out of Freud\n", "\n", "That sweet, sweet RDF data is currently in C++, but we want it in python, so let's get it out of there!" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# get the center of the histogram bins\n", "r = rdf.getR()\n", "# get the value of the histogram bins\n", "y = rdf.getRDF()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualizing your data\n", "\n", "Remember, Freud makes no assumptions about your data or how you want to visualize it. Want to use matplotlib like we do here? Go right ahead. Want to be a heathen and use Excel? Go right ahead, no one is stopping you (except your conscience...don't use Excel)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# create bokeh plot\n", "p = figure(title=\"RDF\", x_axis_label='r', y_axis_label='g(r)')\n", "p.line(r, y, legend=\"g(r)\", line_width=2)\n", "default_bokeh(p)\n", "show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Accumulate, or how I stopped worrying and let Freud do trivial tasks for me\n", "\n", "A lot of times we want to average our data over many simualtion frames. Freud does this for you with the `accumulate` method: " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "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", "n_frames = pos_data.shape[0]\n", "\n", "# compute the rdf for for all frames except the first (your syntax will vary based on your reader)\n", "for i in range(1, n_frames):\n", " # read box, position data\n", " l_box = box_data[i]\n", " l_pos = pos_data[i]\n", " # create the freud box object\n", " fbox = box.Box(Lx=l_box[\"Lx\"], Ly=l_box[\"Ly\"], is2D=True)\n", " # accumulate\n", " rdf.accumulate(fbox, l_pos, l_pos)\n", "\n", "# get the center of the histogram bins\n", "r = rdf.getR()\n", "# get the value of the histogram bins\n", "y = rdf.getRDF()\n", "# create bokeh figure\n", "p = figure(title=\"RDF\", x_axis_label='r', y_axis_label='g(r)')\n", "p.circle(r, y, legend=\"g(r)\", line_width=2)\n", "p.line(r, y, legend=\"g(r)\", line_width=2)\n", "default_bokeh(p)\n", "show(p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# What's the difference?\n", "\n", "Let's plot together and find out:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# reset the rdf; required if not using compute\n", "rdf.resetRDF()\n", "# compute the rdf for for all frames except the first (your syntax will vary based on your reader)\n", "for i in range(1, n_frames):\n", " # read box, position data\n", " l_box = box_data[i]\n", " l_pos = pos_data[i]\n", " # create the freud box object\n", " fbox = box.Box(Lx=l_box[\"Lx\"], Ly=l_box[\"Ly\"], is2D=True)\n", " # accumulate\n", " rdf.accumulate(fbox, l_pos, l_pos)\n", "\n", "# get the center of the histogram bins\n", "r_avg = rdf.getR()\n", "# get the value of the histogram bins\n", "y_avg = rdf.getRDF()\n", "\n", "# compute the rdf for the last frame\n", "# read box, position data\n", "l_box = box_data[-1]\n", "l_pos = pos_data[-1]\n", "# create the freud box object\n", "fbox = box.Box(Lx=l_box[\"Lx\"], Ly=l_box[\"Ly\"], is2D=True)\n", "# compute\n", "rdf.compute(fbox, l_pos, l_pos)\n", "# get the center of the histogram bins\n", "r = rdf.getR()\n", "# get the value of the histogram bins\n", "y = rdf.getRDF()\n", "# create bokeh plot\n", "p0 = figure(title=\"RDF\", x_axis_label='r', y_axis_label='g(r)')\n", "p0.circle(r, y, legend=\"Compute\")\n", "p0.line(r, y, legend=\"Compute\", line_width=2)\n", "p0.square(r_avg, y_avg, legend=\"Accumulate\", fill_color=None, line_color=\"red\")\n", "p0.line(r_avg, y_avg, legend=\"Accumulate\", line_dash=[4,4], line_width=2, line_color=\"red\")\n", "\n", "default_bokeh(p0)\n", "\n", "p1 = figure(title=\"RDF\", x_axis_label='r', y_axis_label='g(r)')\n", "p1.line(r, y, legend=\"Compute\", line_width=2)\n", "p1.line(r_avg, y_avg, legend=\"Accumulate\", line_width=2, line_color=\"red\")\n", "\n", "default_bokeh(p1)\n", "\n", "grid = gridplot([p0, p1], ncols=2, plot_width=400, plot_height=400)\n", "\n", "show(grid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Wait a second, there's no difference?!\n", "\n", "Right you are, but there should be..." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# reset the rdf; required if not using compute\n", "rdf.resetRDF()\n", "# compute the rdf for for all frames except the first (your syntax will vary based on your reader)\n", "for i in range(1, n_frames):\n", " # read box, position data\n", " l_box = box_data[i]\n", " l_pos = pos_data[i]\n", " # create the freud box object\n", " fbox = box.Box(Lx=l_box[\"Lx\"], Ly=l_box[\"Ly\"], is2D=True)\n", " # accumulate\n", " rdf.accumulate(fbox, l_pos, l_pos)\n", "\n", "# USE NP.COPY!\n", "# get the center of the histogram bins\n", "r_avg = np.copy(rdf.getR())\n", "# get the value of the histogram bins\n", "y_avg = np.copy(rdf.getRDF())\n", "\n", "# compute the rdf for the last frame\n", "# read box, position data\n", "l_box = box_data[-1]\n", "l_pos = pos_data[-1]\n", "# create the freud box object\n", "fbox = box.Box(Lx=l_box[\"Lx\"], Ly=l_box[\"Ly\"], is2D=True)\n", "# compute\n", "rdf.compute(fbox, l_pos, l_pos)\n", "# get the center of the histogram bins\n", "r = rdf.getR()\n", "# get the value of the histogram bins\n", "y = rdf.getRDF()\n", "# create bokeh plot\n", "p0 = figure(title=\"RDF\", x_axis_label='r', y_axis_label='g(r)')\n", "p0.circle(r, y, legend=\"Compute\")\n", "p0.line(r, y, legend=\"Compute\", line_width=2)\n", "p0.square(r_avg, y_avg, legend=\"Accumulate\", fill_color=None, line_color=\"red\")\n", "p0.line(r_avg, y_avg, legend=\"Accumulate\", line_dash=[4,4], line_width=2, line_color=\"red\")\n", "\n", "default_bokeh(p0)\n", "\n", "p1 = figure(title=\"RDF\", x_axis_label='r', y_axis_label='g(r)')\n", "p1.line(r, y, legend=\"Compute\", line_width=2)\n", "p1.line(r_avg, y_avg, legend=\"Accumulate\", line_width=2, line_color=\"red\")\n", "\n", "default_bokeh(p1)\n", "\n", "grid = gridplot([p0, p1], ncols=2, plot_width=400, plot_height=400)\n", "\n", "show(grid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Pointers vs. Copies\n", "\n", "By default Freud returns a NumPy array *as a pointer*. This is done for speed, but can result in the above problem. Please be sure to copy your data out as needed. A more in-depth discussion can be found in the [NumPy arrays from pointers example notebook](CopyVsPointers.ipynb)." ] } ], "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.2" } }, "nbformat": 4, "nbformat_minor": 0 }