freud.density module is intended to compute a variety of quantities that relate spatial distributions of particles with other particles. In this notebook, we demonstrate
freud’s Gaussian density calculation, which provides a way to interpolate particle configurations onto a regular grid in a meaningful way that can then be processed by other algorithms that require regularity, such as a Fast Fourier Transform.
import freud import matplotlib.pyplot as plt import numpy as np from scipy import stats
To illustrate the basic concept, consider a toy example: a simple set of point particles with unit mass on a line. For analytical purposes, the standard way to accomplish this would be using Dirac delta functions.
n_p = 10000 np.random.seed(129) x = np.linspace(0, 1, n_p) y = np.zeros(n_p) points = np.random.rand(10) y[(points * n_p).astype("int")] = 1 plt.plot(x, y) plt.show()
However, delta functions can be cumbersome to work with, so we might instead want to smooth out these particles. One option is to instead represent particles as Gaussians centered at the location of the points. In that case, the total particle density at any point in the interval \([0, 1]\) represented above would be based on the sum of the densities of those Gaussians at those points.
# Note that we use a Gaussian with a small standard deviation # to emphasize the differences on this small scale dists = [stats.norm(loc=i, scale=0.1) for i in points] y_gaussian = 0 for dist in dists: y_gaussian += dist.pdf(x) plt.plot(x, y_gaussian) plt.show()
The goal of the GaussianDensity class is to perform the same interpolation for points on a 2D or 3D grid, accounting for Box periodicity.
N = 1000 # Number of points L = 10 # Box length box, points = freud.data.make_random_system(L, N, is2D=True, seed=0) aq = freud.AABBQuery(box, points) gd = freud.density.GaussianDensity(L * L, L / 3, 1) gd.compute(aq) fig, axes = plt.subplots(1, 2, figsize=(14, 6)) aq.plot(ax=axes) gd.plot(ax=axes) plt.show()
The effects are much more striking if we explicitly construct our points to be centered at certain regions.
N = 1000 # Number of points L = 10 # Box length box = freud.box.Box.square(L) centers = np.array( [[L / 4, L / 4, 0], [-L / 4, L / 4, 0], [L / 4, -L / 4, 0], [-L / 4, -L / 4, 0]] ) points =  for center in centers: points.append( np.random.multivariate_normal( center, cov=np.diag([1, 1, 0]), size=(int(N / 4),) ) ) points = box.wrap(np.concatenate(points)) aq = freud.AABBQuery(box, points) gd = freud.density.GaussianDensity(L * L, L / 3, 1) gd.compute(aq) fig, axes = plt.subplots(1, 2, figsize=(14, 6)) aq.plot(ax=axes) gd.plot(ax=axes) plt.show()