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 numpy as np
import freud
from bokeh.io import output_notebook
from bokeh.plotting import figure, show
import util
Loading BokehJS ...
def plot_hex_order_param(data_path, title):
    # Create hexatic object
    hex_order = freud.order.HexOrderParameter(rmax=1.2, k=6, n=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
    l_box = box_data[-1].tolist()
    l_pos = pos_data[-1]
    l_quat = quat_data[-1]
    l_ang = 2*np.arctan2(l_quat[:, 3], l_quat[:, 0])

    # Compute hexatic order for 6 nearest neighbors
    hex_order.compute(l_box, l_pos)
    psi_k = hex_order.psi
    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, l_pos[:, :2], l_ang)
    # 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],
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')