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

[1]:
import numpy as np
import freud
from bokeh.io import output_notebook
from bokeh.plotting import figure, show
# The following line imports this set of utility functions:
# https://github.com/glotzerlab/freud-examples/blob/master/util.py
import util
output_notebook()
Loading BokehJS ...
[2]:
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("{}/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))

    # 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)
    p.patches(xs=patches[:, :, 0].tolist(), ys=patches[:, :, 1].tolist(),
              fill_color=[util.cubeellipse(x) for x in relative_angles],
              line_color="black")
    util.default_bokeh(p)
    show(p)
[3]:
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.

[4]:
plot_hex_order_param('data/phi075', 'Hexatic Order Parameter, 0.75 density')

As the density increases to \(\phi=0.85\), the alignment becomes even stronger and defects are no longer visible.

[5]:
plot_hex_order_param('data/phi085', 'Hexatic Order Parameter, 0.85 density')