Examples are given here in abbreviated form. We invite you to look over our other examples which explain these calculations in greater detail, as well as read the freud documentation to better understand these examples
# import bokeh for plotting; use the inline (local) files
# so that an internet connection is not needed
from bokeh.resources import INLINE
from bokeh.io import output_notebook
output_notebook(resources=INLINE)
from bokeh.plotting import figure, output_file, show
from bokeh.layouts import gridplot
# import numpy
import numpy as np
# import freud and set number of parallel threads
from freud import parallel
parallel.setNumThreads(4)
# import ipywidgets for display purposes
from ipywidgets import IntProgress
from IPython.display import display
# import time to display computation time
import time
from aps_helpers import *
While Freud does not provide any methods for visualization, Freud allows you to enhance your visualizations in ways which make them easier to understand. Here we demonstrate how to visualize a trajectory for a simple demonstration system.
Freud does not contain any file readers or writers. For this demonstration I have saved the relevant information as numpy binary files.
# load the data
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))
quat_data = np.load("{}/quat_data.npy".format(data_path))
n_frames = pos_data.shape[0]
Now we visualize our system, using bokeh to draw our shapes
# import freud box class
from freud import box
# grab data from last frame
l_box = box_data[-1]
l_pos = pos_data[-1]
l_quat = quat_data[-1]
l_ang = 2*np.arctan2(np.copy(l_quat[:,3]), np.copy(l_quat[:,0]))
# create box
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
side_length = max(fbox.Lx, fbox.Ly)
l_min = -side_length / 2.0
l_min *= 1.1
l_max = -l_min
# take local vertices and rotate, translate into system coordinates
patches = local_to_global(verts, l_pos[:,0:2], l_ang)
# plot
p = figure(title="System Visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),
fill_color=(42,126,187), line_color="black", line_width=1.5, legend="hexagons")
# box display
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
# create a grid for this
Just viewing particle positions and orientations doesn't give us much information. Let's consider two different ways of looking at our system in more detail: the number of nearest neighbors a particle has, and the angle of the hexatic order parameter for a particle.
First, let's view the number of nearest neighbors each particle has. We'll start by just showing a single particle with its 6 nearest neighbors:
# import freud locality object
from freud import locality
# create freud nearest neighbor object
# set number of neighbors
n_neigh = 6
# create freud nearest neighbors object
nn = locality.NearestNeighbors(rmax=1.5,n_neigh=n_neigh,strict_cut=False)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = nn.getNeighborList()
# get the number of particles
num_particles = nn.getNRef()
# get the neighbors for particle 1000
pidx = 1000
n_idxs = n_list[pidx]
# get position, orientation for the central particle
center_pos = np.zeros(shape=(1, 3),dtype=np.float32)
center_ang = np.zeros(shape=(1),dtype=np.float32)
center_pos[:] = l_pos[pidx]
center_ang[:] = l_ang[pidx]
# get the positions, orientations for the neighbor particles
neigh_pos = np.zeros(shape=(n_neigh, 3),dtype=np.float32)
neigh_ang = np.zeros(shape=(n_neigh),dtype=np.float32)
neigh_pos[:] = l_pos[n_idxs]
neigh_ang[:] = l_ang[n_idxs]
# render in bokeh
# create array of transformed positions
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
n_patches = local_to_global(verts, neigh_pos[:,0:2], neigh_ang)
# turn into list of colors
# bokeh (as of this version) requires hex colors, so convert rgb to hex
center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
# plot
p = figure(title="Nearest Neighbors visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
p.patches(xs=n_patches[:,:,0].tolist(), ys=n_patches[:,:,1].tolist(),
fill_color=neigh_color.tolist(), line_color="black", legend="neighbors")
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="centers")
# box display
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.0209500789642334
Note that there appears to be an outlier particle, one that is much farther away from the central particle than the others. We can use freud to detect and not count that particle:
# set number of neighbors
n_neigh = 6
# create freud nearest neighbors object
nn = locality.NearestNeighbors(rmax=1.5,n_neigh=n_neigh,strict_cut=True)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = nn.getNeighborList()
# get the number of particles
num_particles = nn.getNRef()
# get the neighbors for particle 1000
pidx = 1000
n_idxs = n_list[pidx]
# clip padded values
n_idxs = n_idxs[np.where(n_idxs < num_particles)]
n_neigh = len(n_idxs)
# get position, orientation for the central particle
center_pos = np.zeros(shape=(1, 3),dtype=np.float32)
center_ang = np.zeros(shape=(1),dtype=np.float32)
center_pos[:] = l_pos[pidx]
center_ang[:] = l_ang[pidx]
# get the positions, orientations for the neighbor particles
neigh_pos = np.zeros(shape=(n_neigh, 3),dtype=np.float32)
neigh_ang = np.zeros(shape=(n_neigh),dtype=np.float32)
neigh_pos[:] = l_pos[n_idxs]
neigh_ang[:] = l_ang[n_idxs]
# render in bokeh
# create array of transformed positions
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
n_patches = local_to_global(verts, neigh_pos[:,0:2], neigh_ang)
# turn into list of colors
# bokeh (as of this version) requires hex colors, so convert rgb to hex
center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
# plot
p = figure(title="Nearest Neighbors visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
p.patches(xs=n_patches[:,:,0].tolist(), ys=n_patches[:,:,1].tolist(),
fill_color=neigh_color.tolist(), line_color="black", legend="neighbors")
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="centers")
# box display
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.00730586051940918
# create freud nearest neighbor object
# set number of neighbors
n_neigh = 8
# create freud nearest neighbors object
# set rmax to some value
rmax=1.7
nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = np.copy(nn.getNeighborList())
# get the number of particles
num_particles = nn.getNRef()
# now for array manipulation magic
# create an integer array of the same shape as the neighbor list array
int_arr = np.ones(shape=n_list.shape, dtype=np.int32)
# "search" for non-indexed particles (missing neighbors)
# while it would be most accurate to use the UINTMAX value
# provided by nn.getUINTMAX(), but this works just as well
int_arr[n_list > (num_particles-1)] = 0
# sum along particle index axis to determine the number of neighbors per particle
n_neighbors = np.sum(int_arr, axis=1)
# find the complement (if desired) to find number of missing neighbors per particle
n_deficits = n_neigh - n_neighbors
p = figure(title="Nearest Neighbors visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
for k in np.unique(n_neighbors):
# find particles with k neighbors
c_idxs = np.copy(np.where(n_neighbors==k)[0])
center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)
center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)
center_pos = l_pos[c_idxs]
center_ang = l_ang[c_idxs]
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.006757020950317383
In addition to viewing the number of neighbors to find defects, freud can also compute a variety of order parameters. In this example we color each hexagon by its hexatic order parameter angle so as to better see the overall orientation of neighbors in the system.
# import freud order module
from freud import order
# create hexatic order parameter object
hex_order = order.HexOrderParameter(rmax=1.2, k=6, n=6);
# compute hexatic order for 6 nearest neighbors
start_time = time.time()
hex_order.compute(fbox, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get values from freud object
psi_k = hex_order.getPsi()
avg_psi_k = np.mean(psi_k)
# render in bokeh
# create array of transformed positions
patches = local_to_global(verts, l_pos[:,0:2], l_ang)
# create an array of angles relative to the average
a = np.angle(psi_k) - np.angle(avg_psi_k)
# turn into an rgb array of tuples
color = [tuple(cubeellipse(x)) for x in a]
# bokeh (as of this version) requires hex colors, so convert rgb to hex
hex_color = ["#{0:02x}{1:02x}{2:02x}".format(clamp(r), clamp(g), clamp(b)) for (r,g,b) in color]
# plot
p = figure(title="Hexatic Order Parameter visualization")
p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),
fill_color=hex_color, line_color="black")
default_bokeh(p)
show(p)
time to calc 1 frame = 0.046263933181762695
# load the data
data_path = "ex_data/phi075"
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
quat_data = np.load("{}/quat_data.npy".format(data_path))
n_frames = pos_data.shape[0]
# grab data from last frame
l_box = box_data[-1]
l_pos = pos_data[-1]
l_quat = quat_data[-1]
l_ang = 2*np.arctan2(np.copy(l_quat[:,3]), np.copy(l_quat[:,0]))
# create box
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
side_length = max(fbox.Lx, fbox.Ly)
l_min = -side_length / 2.0
l_min *= 1.1
l_max = -l_min
# compute hexatic order for 6 nearest neighbors
start_time = time.time()
hex_order.compute(fbox, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get values from freud object
psi_k = hex_order.getPsi()
avg_psi_k = np.mean(psi_k)
# render in bokeh
# create array of transformed positions
patches = local_to_global(verts, l_pos[:,0:2], l_ang)
# create an array of angles relative to the average
a = np.angle(psi_k) - np.angle(avg_psi_k)
# turn into an rgb array of tuples
color = [tuple(cubeellipse(x)) for x in a]
# bokeh (as of this version) requires hex colors, so convert rgb to hex
hex_color = ["#{0:02x}{1:02x}{2:02x}".format(clamp(r), clamp(g), clamp(b)) for (r,g,b) in color]
# plot
p = figure(title="Hexatic Order Parameter visualization")
p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),
fill_color=hex_color, line_color="black")
default_bokeh(p)
show(p)
time to calc 1 frame = 0.030714988708496094
In addtion to computing per-particle quanities and order parameters useful for visualization of simulations, freud can also compute global quantities for simulations, useful for determining thermodynamic quantities. Here we demonstrate the calculation of the Radial Distribution Function (RDF or $ g(r) $).
First, we compute the RDF for a single trajectory frame:
from freud import box, density
# reload the data for the lower-density system
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))
quat_data = np.load("{}/quat_data.npy".format(data_path))
n_frames = pos_data.shape[0]
# create the rdf object
rdf = density.RDF(rmax=5.0, dr=0.01)
# 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
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)
We often want to average our data over a number of simulation frames. Rather than use numpy or other libraries to do this, freud provides for averaging your RDF histograms within its computation loop:
rdf.resetRDF()
# compute the rdf for for all frames except the first (your syntax will vary based on your reader)
myProgressBar = IntProgress(min=1,max=n_frames)
display(myProgressBar)
start_time = time.time()
for i in range(1, n_frames):
# read box, position data
myProgressBar.value = i
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)
stop_time = time.time()
print("time to calc {} frames = {}".format(n_frames-1, stop_time-start_time))
print("speed of calc: {} (frames/sec)".format((n_frames-1)/(stop_time-start_time)))
# get the center of the histogram bins
r = rdf.getR()
# get the value of the histogram bins
y = rdf.getRDF()
# 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)
time to calc 399 frames = 6.444867134094238 speed of calc: 61.909732458134755 (frames/sec)
rdf.resetRDF()
# load the data
data_path = "ex_data/phi075"
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)
myProgressBar = IntProgress(min=1,max=n_frames)
display(myProgressBar)
start_time = time.time()
for i in range(1, n_frames):
# read box, position data
myProgressBar.value = i
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)
stop_time = time.time()
print("time to calc {} frames = {}".format(n_frames-1, stop_time-start_time))
print("speed of calc: {} (frames/sec)".format((n_frames-1)/(stop_time-start_time)))
# get the center of the histogram bins
r = rdf.getR()
# get the value of the histogram bins
y = rdf.getRDF()
# 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)
time to calc 399 frames = 7.012582063674927 speed of calc: 56.89772987710393 (frames/sec)