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!
To calculate the radial distribution function TODO (add in math to actual section?), we will need to follow the following workflow:
We will use the following python packages, please install with your favorite package manager:
numpy
- general scientific goodnessbokeh
- plottingmatplotlib
- plottingfrom bokeh.io import output_notebook
output_notebook()
from bokeh.plotting import figure, output_file, show
from bokeh.layouts import gridplot
import numpy as np
from freud import parallel, box, density
parallel.setNumThreads(4)
def default_bokeh(p):
p.title.text_font_size = "18pt"
p.title.align = "center"
p.xaxis.axis_label_text_font_size = "14pt"
p.yaxis.axis_label_text_font_size = "14pt"
p.xaxis.major_tick_in = 10
p.xaxis.major_tick_out = 0
p.xaxis.minor_tick_in = 5
p.xaxis.minor_tick_out = 0
p.yaxis.major_tick_in = 10
p.yaxis.major_tick_out = 0
p.yaxis.minor_tick_in = 5
p.yaxis.minor_tick_out = 0
p.xaxis.major_label_text_font_size = "12pt"
p.yaxis.major_label_text_font_size = "12pt"
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
import freud
as you see fit.
from freud import parallel, box, density
parallel.setNumThreads(4)
We imported the following modules:
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.
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.
Let's create the RDF object. For your convenience, we use the help function to show what needs to be supplied:
help(density.RDF)
Help on class RDF in module freud._freud: class RDF(builtins.object) | Computes RDF for supplied data | | 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. | | The values of :math:`r` to compute the rdf are set by the values of rmax, dr in the constructor. rmax sets the maximum | distance at which to calculate the :math:`g \left( r \right)` while dr determines the step size for each bin. | | .. moduleauthor:: Eric Harper <harperic@umich.edu> | | .. note:: | 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. | | :param rmax: maximum distance to calculate | :param dr: distance between histogram bins | :type rmax: float | :type dr: float | | Methods defined here: | | __new__(*args, **kwargs) from builtins.type | Create and return a new object. See help(type) for accurate signature. | | accumulate(...) | RDF.accumulate(self, box, ndarray ref_points, ndarray points) | | Calculates the rdf and adds to the current rdf histogram. | | :param box: simulation box | :param ref_points: reference points to calculate the local density | :param points: points to calculate the local density | :type box: :py:class:`freud.box.Box` | :type ref_points: :class:`numpy.ndarray`, shape=(:math:`N_{particles}`, 3), dtype= :class:`numpy.float32` | :type points: :class:`numpy.ndarray`, shape=(:math:`N_{particles}`, 3), dtype= :class:`numpy.float32` | | compute(...) | RDF.compute(self, box, ndarray ref_points, ndarray points) | | Calculates the rdf for the specified points. Will overwrite the current histogram. | | :param box: simulation box | :param ref_points: reference points to calculate the local density | :param points: points to calculate the local density | :type box: :py:meth:`freud.box.Box` | :type ref_points: :class:`numpy.ndarray`, shape=(:math:`N_{particles}`, 3), dtype= :class:`numpy.float32` | :type points: :class:`numpy.ndarray`, shape=(:math:`N_{particles}`, 3), dtype= :class:`numpy.float32` | | getBox(...) | RDF.getBox(self) | | :return: Freud Box | :rtype: :py:class:`freud.box.Box` | | getNr(...) | RDF.getNr(self) | | :return: histogram of cumulative rdf values | :rtype: :class:`numpy.ndarray`, shape=(:math:`N_{bins}`, 3), dtype= :class:`numpy.float32` | | getR(...) | RDF.getR(self) | | :return: values of the histogram bin centers | :rtype: :class:`numpy.ndarray`, shape=(:math:`N_{bins}`, 3), dtype= :class:`numpy.float32` | | getRDF(...) | RDF.getRDF(self) | | :return: histogram of rdf values | :rtype: :class:`numpy.ndarray`, shape=(:math:`N_{bins}`, 3), dtype= :class:`numpy.float32` | | reduceRDF(...) | RDF.reduceRDF(self) | | Reduces the histogram in the values over N processors to a single histogram. This is called automatically by | :py:meth:`freud.density.RDF.getRDF()`, :py:meth:`freud.density.RDF.getNr()`. | | resetRDF(...) | RDF.resetRDF(self) | | resets the values of RDF in memory
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
:
rdf = density.RDF(rmax=10.0, dr=0.1)
It's time to actually compute the rdf, so here we go, what is happening is briefly explained in the comments
data_path = "ex_data/phi065"
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
n_frames = pos_data.shape[0]
# compute the rdf for the last frame
# read box, position data
l_box = box_data[-1]
l_pos = pos_data[-1]
# create the freud box object
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
# compute
rdf.compute(fbox, l_pos, l_pos)
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.
That sweet, sweet RDF data is currently in C++, but we want it in python, so let's get it out of there!
# get the center of the histogram bins
r = rdf.getR()
# get the value of the histogram bins
y = rdf.getRDF()
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)
# create bokeh plot
p = figure(title="RDF", x_axis_label='r', y_axis_label='g(r)')
p.line(r, y, legend="g(r)", line_width=2)
default_bokeh(p)
show(p)
A lot of times we want to average our data over many simualtion frames. Freud does this for you with the accumulate
method:
data_path = "ex_data/phi065"
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
n_frames = pos_data.shape[0]
# compute the rdf for for all frames except the first (your syntax will vary based on your reader)
for i in range(1, n_frames):
# read box, position data
l_box = box_data[i]
l_pos = pos_data[i]
# create the freud box object
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
# accumulate
rdf.accumulate(fbox, l_pos, l_pos)
# get the center of the histogram bins
r = rdf.getR()
# get the value of the histogram bins
y = rdf.getRDF()
# create bokeh figure
p = figure(title="RDF", x_axis_label='r', y_axis_label='g(r)')
p.circle(r, y, legend="g(r)", line_width=2)
p.line(r, y, legend="g(r)", line_width=2)
default_bokeh(p)
show(p)
Let's plot together and find out:
# reset the rdf; required if not using compute
rdf.resetRDF()
# compute the rdf for for all frames except the first (your syntax will vary based on your reader)
for i in range(1, n_frames):
# read box, position data
l_box = box_data[i]
l_pos = pos_data[i]
# create the freud box object
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
# accumulate
rdf.accumulate(fbox, l_pos, l_pos)
# get the center of the histogram bins
r_avg = rdf.getR()
# get the value of the histogram bins
y_avg = rdf.getRDF()
# compute the rdf for the last frame
# read box, position data
l_box = box_data[-1]
l_pos = pos_data[-1]
# create the freud box object
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
# compute
rdf.compute(fbox, l_pos, l_pos)
# get the center of the histogram bins
r = rdf.getR()
# get the value of the histogram bins
y = rdf.getRDF()
# create bokeh plot
p0 = figure(title="RDF", x_axis_label='r', y_axis_label='g(r)')
p0.circle(r, y, legend="Compute")
p0.line(r, y, legend="Compute", line_width=2)
p0.square(r_avg, y_avg, legend="Accumulate", fill_color=None, line_color="red")
p0.line(r_avg, y_avg, legend="Accumulate", line_dash=[4,4], line_width=2, line_color="red")
default_bokeh(p0)
p1 = figure(title="RDF", x_axis_label='r', y_axis_label='g(r)')
p1.line(r, y, legend="Compute", line_width=2)
p1.line(r_avg, y_avg, legend="Accumulate", line_width=2, line_color="red")
default_bokeh(p1)
grid = gridplot([p0, p1], ncols=2, plot_width=400, plot_height=400)
show(grid)
Right you are, but there should be...
# reset the rdf; required if not using compute
rdf.resetRDF()
# compute the rdf for for all frames except the first (your syntax will vary based on your reader)
for i in range(1, n_frames):
# read box, position data
l_box = box_data[i]
l_pos = pos_data[i]
# create the freud box object
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
# accumulate
rdf.accumulate(fbox, l_pos, l_pos)
# USE NP.COPY!
# get the center of the histogram bins
r_avg = np.copy(rdf.getR())
# get the value of the histogram bins
y_avg = np.copy(rdf.getRDF())
# compute the rdf for the last frame
# read box, position data
l_box = box_data[-1]
l_pos = pos_data[-1]
# create the freud box object
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
# compute
rdf.compute(fbox, l_pos, l_pos)
# get the center of the histogram bins
r = rdf.getR()
# get the value of the histogram bins
y = rdf.getRDF()
# create bokeh plot
p0 = figure(title="RDF", x_axis_label='r', y_axis_label='g(r)')
p0.circle(r, y, legend="Compute")
p0.line(r, y, legend="Compute", line_width=2)
p0.square(r_avg, y_avg, legend="Accumulate", fill_color=None, line_color="red")
p0.line(r_avg, y_avg, legend="Accumulate", line_dash=[4,4], line_width=2, line_color="red")
default_bokeh(p0)
p1 = figure(title="RDF", x_axis_label='r', y_axis_label='g(r)')
p1.line(r, y, legend="Compute", line_width=2)
p1.line(r_avg, y_avg, legend="Accumulate", line_width=2, line_color="red")
default_bokeh(p1)
grid = gridplot([p0, p1], ncols=2, plot_width=400, plot_height=400)
show(grid)
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.