freud.pmft.PMFTXY¶
The PMFT returns the potential energy associated with finding a particle pair in a given spatial (positional and orientational) configuration. The PMFT is computed in the same manner as the RDF. The basic algorithm is described below:
for each particle i:
for each particle j:
v_ij = position[j] - position[i]
bin_x, bin_y = convert_to_bin(v_ij)
pcf_array[bin_y][bin_x]++
freud
uses spatial data structures and parallelism to optimize this algorithm.
The data sets used in this example are a system of hard hexagons, simulated in the NVT thermodynamic ensemble in HOOMD-blue, for a dense fluid of hexagons at packing fraction \(\phi = 0.65\) and solids at packing fractions \(\phi = 0.75, 0.85\).
[1]:
import freud
freud.set_num_threads(1)
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage.filters import gaussian_filter
# The following line imports this set of utility functions:
# https://github.com/glotzerlab/freud-examples/blob/master/util.py
import util
%matplotlib inline
matplotlib.rcParams.update(
{
"font.size": 20,
"axes.titlesize": 20,
"axes.labelsize": 20,
"xtick.labelsize": 16,
"ytick.labelsize": 16,
"savefig.pad_inches": 0.025,
"lines.linewidth": 2,
}
)
[2]:
def plot_pmft(data_path, phi):
# Create the pmft object
pmft = freud.pmft.PMFTXY(x_max=3.0, y_max=3.0, bins=300)
# Load the data
box_data = np.load(f"{data_path}/box_data.npy")
pos_data = np.load(f"{data_path}/pos_data.npy")
quat_data = np.load(f"{data_path}/quat_data.npy")
n_frames = pos_data.shape[0]
for i in range(1, n_frames):
# Read box, position data
box = box_data[i].tolist()
points = pos_data[i]
quats = quat_data[i]
angles = 2 * np.arctan2(quats[:, 3], quats[:, 0]) % (2 * np.pi)
pmft.compute(system=(box, points), query_orientations=angles, reset=False)
# Get the value of the PMFT histogram bins
pmft_arr = pmft.pmft.T
# Do some simple post-processing for plotting purposes
pmft_arr[np.isinf(pmft_arr)] = np.nan
dx = (2.0 * 3.0) / pmft.nbins[0]
dy = (2.0 * 3.0) / pmft.nbins[1]
nan_arr = np.where(np.isnan(pmft_arr))
for i in range(pmft.nbins[0]):
x = -3.0 + dx * i
for j in range(pmft.nbins[1]):
y = -3.0 + dy * j
if (x * x + y * y < 1.5) and (np.isnan(pmft_arr[j, i])):
pmft_arr[j, i] = 10.0
w = int(2.0 * pmft.nbins[0] / (2.0 * 3.0))
center = int(pmft.nbins[0] / 2)
# Get the center of the histogram bins
pmft_smooth = gaussian_filter(pmft_arr, 1)
pmft_image = np.copy(pmft_smooth)
pmft_image[nan_arr] = np.nan
pmft_smooth = pmft_smooth[center - w : center + w, center - w : center + w]
pmft_image = pmft_image[center - w : center + w, center - w : center + w]
x, y = pmft.bin_centers
reduced_x = x[center - w : center + w]
reduced_y = y[center - w : center + w]
# Plot figures
f = plt.figure(figsize=(12, 5), facecolor="white")
values = [-2, -1, 0, 2]
norm = matplotlib.colors.Normalize(vmin=-2.5, vmax=3.0)
n_values = [norm(i) for i in values]
colors = matplotlib.cm.viridis(n_values)
colors = colors[:, :3]
verts = util.make_polygon(sides=6, radius=0.6204)
lims = (-2, 2)
ax0 = f.add_subplot(1, 2, 1)
ax1 = f.add_subplot(1, 2, 2)
for ax in (ax0, ax1):
ax.contour(reduced_x, reduced_y, pmft_smooth, [9, 10], colors="black")
ax.contourf(
reduced_x, reduced_y, pmft_smooth, [9, 10], hatches="X", colors="none"
)
ax.plot(verts[:, 0], verts[:, 1], color="black", marker=",")
ax.fill(verts[:, 0], verts[:, 1], color="black")
ax.set_aspect("equal")
ax.set_xlim(lims)
ax.set_ylim(lims)
ax.xaxis.set_ticks([i for i in range(lims[0], lims[1] + 1)])
ax.yaxis.set_ticks([i for i in range(lims[0], lims[1] + 1)])
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax0.set_title(rf"PMFT Heat Map, $\phi = {phi}$")
im = ax0.imshow(
np.flipud(pmft_image),
extent=[lims[0], lims[1], lims[0], lims[1]],
interpolation="nearest",
cmap="viridis",
vmin=-2.5,
vmax=3.0,
)
ax1.set_title(rf"PMFT Contour Plot, $\phi = {phi}$")
ax1.contour(reduced_x, reduced_y, pmft_smooth, [-2, -1, 0, 2], colors=colors)
f.subplots_adjust(right=0.85)
cbar_ax = f.add_axes([0.88, 0.1, 0.02, 0.8])
f.colorbar(im, cax=cbar_ax)
plt.show()
65% density¶
The plot below shows the PMFT of hexagons at 65% density. The hexagons tend to be close to one another, in the darker regions (the lower values of the potential of mean force and torque).
The hatched region near the black hexagon in the center is a region where no data were collected: the hexagons are hard shapes and cannot overlap, so there is an excluded region of space close to the hexagon.
The ring around the hexagon where the PMFT rises and then falls corresponds to the minimum of the radial distribution function – particles tend to not occupy that region, preferring instead to be at close range (in the first neighbor shell) or further away (in the second neighbor shell).
[3]:
plot_pmft("data/phi065", 0.65)
75% density¶
As the system density is increased to 75%, the propensity for hexagons to occupy the six sites on the faces of their neighbors increases, as seen by the deeper (darker) wells of the PMFT. Conversely, the shapes strongly dislike occupying the yellow regions, and no particle pairs occupied the white region (so there is no data).
[4]:
plot_pmft("data/phi075", 0.75)
85% density¶
Finally, at 85% density, there is a large region where no neighbors can be found, and hexagons strictly occupy sites near those of the perfect hexagonal lattice, at the first- and second-neighbor shells. The wells are deeper and much more spatially confined that those of the systems at lower densities.
[5]:
plot_pmft("data/phi085", 0.85)