freud.order.Hexatic: Hard Hexagons

Hexatic Order Parameter

The hexatic order parameter measures how closely the local environment around a particle resembles perfect \(k\)-atic symmetry, e.g. how closely the environment resembles hexagonal/hexatic symmetry for \(k=6\). The order parameter is given by:

\[\psi_k \left( i \right) = \frac{1}{n} \sum \limits_j^n e^{k i \theta_{ij}}\]

where \(\theta_{ij}\) is the angle between the vector \(\vec{r}_{ij}\) and \(\left(1, 0\right)\).

The pseudocode is given below:

for each particle i:
    neighbors = nearestNeighbors(i, n):
    for each particle j in neighbors:
        r_ij = position[j] - position[i]
        theta_ij = arctan2(r_ij.y, r_ij.x)
        psi_array[i] += exp(complex(0,k*theta_ij))

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\).

import freud
import numpy as np
from import output_notebook
from bokeh.plotting import figure, show

# The following line imports this set of utility functions:
import util

Loading BokehJS ...
def plot_hex_order_param(data_path, title):
    # Create hexatic object
    hex_order = freud.order.Hexatic(k=6)

    # 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")

    # Grab data from last frame
    box = box_data[-1].tolist()
    points = pos_data[-1]
    quats = quat_data[-1]
    angles = 2 * np.arctan2(quats[:, 3], quats[:, 0])

    # Compute hexatic order for 6 nearest neighbors
    hex_order.compute(system=(box, points), neighbors={"num_neighbors": 6})
    psi_k = hex_order.particle_order
    avg_psi_k = np.mean(psi_k)

    # Create hexagon vertices
    verts = util.make_polygon(sides=6, radius=0.6204)
    # Create array of transformed positions
    patches = util.local_to_global(verts, points[:, :2], angles)
    # Create an array of angles relative to the average
    relative_angles = np.angle(psi_k) - np.angle(avg_psi_k)
    # Plot in bokeh
    p = figure(title=title)
        xs=patches[:, :, 0].tolist(),
        ys=patches[:, :, 1].tolist(),
        fill_color=[util.cubeellipse(x) for x in relative_angles],
plot_hex_order_param("data/phi065", "Hexatic Order Parameter, 0.65 density")

As the density increases to \(\phi=0.75\), the shapes are forced to align more closely so that they may tile space effectively.

plot_hex_order_param("data/phi075", "Hexatic Order Parameter, 0.75 density")