Grain Size Determination

The freud.density module is intended to compute a variety of quantities that relate spatial distributions of particles with other particles. This example shows how correlation functions can be used to approximate the grain size within a system.

import freud
import numpy as np
import matplotlib.pyplot as plt

To show the correlation function, we need a pretend data set. We start with a phase separated Ising lattice and assign type values, either blue (-1) or red (+1), to generate “grains.”

To make the pretend data set, we create a large number of blue (-1) particles on a square grid. Then we place grain centers on a larger grid and draw grain radii from a normal distribution. We color the particles red (+1) if their distance from a grain center is less than the grain radius.

# Set up the system
box = freud.box.Box.square(L=10)
dx = 0.15
num_grains = 4
dg = box.Lx/num_grains
points = np.array([[i, j, 0]
                   for j in np.arange(-box.Ly/2, box.Ly/2, dx)
                   for i in np.arange(-box.Lx/2, box.Lx/2, dx)])
values = np.array([-1 for p in points])
centroids = [[i*dg + 0.5*dg, j*dg + 0.5*dg, 0]
             for i in range(num_grains) for j in range(num_grains)]
grain_radii = np.abs(np.random.normal(size=num_grains**2, loc=0.25*dg, scale=0.05*dg))
for center, radius in zip(centroids, grain_radii):
    lc = freud.locality.LinkCell(box, radius).compute(box, points, [center])
    for i in lc.nlist.index_i:
        values[i] = 1

plt.figure(figsize=(8, 8))
plt.scatter(points[values > 0, 0],
            points[values > 0, 1],
            marker='o', color='red', s=25)
plt.scatter(points[values < 0, 0],
            points[values < 0, 1],
            marker='o', color='blue', s=25)

This system is phase-separated because the red particles are generally near one another, and so are the blue particles. However, there is some length scale at which the phase separation is no longer visible. Imagine looking at this image from a far distance away: the red and blue regions would be indistinguishable, and the picture would appear purple.

We can use the correlation lengths between the grain centers and the surrounding red and blue particles to find where the grains end. First, we need to locate the centers of the grains using freud.cluster.

cl = freud.cluster.Cluster(box=box, rcut=2*dx)
cluster_idx = cl.computeClusters(points[values > 0, :]).cluster_idx

Now we’ll find the cluster centroids (the grain centers).

clp = freud.cluster.ClusterProperties(box=box)
clp.computeProperties(points[values > 0, :], cl.cluster_idx)
freud.cluster.ClusterProperties(box=freud.box.Box(Lx=10.0, Ly=10.0, Lz=0.0, xy=0.0, xz=0.0, yz=0.0, is2D=True))

Now we create a freud.density.FloatCF compute object, to compute the correlation function.

fcf = freud.density.FloatCF(rmax=0.5*box.Lx/num_grains, dr=dx)

Now we compute the product of type values. Knowing that two sites with the same type will have a product of 1, we can estimate the radius of a grain as the smallest value of \(r\) such that \(\langle p \ast q \rangle (r) < 0\), where the set \(p\) represents our red and blue particles and \(q\) represents the grain centers.

plt.figure(figsize=(10, 5))
weights = []
sizes = []
measured_grain_radii = []
for cluster_id in range(cl.num_clusters):
                ref_points=points,  # all points in the system
                ref_values=values,   # all types in the system
    # get the center of the histogram bins
    r = fcf.R
    # get the value of the histogram bins
    y = fcf.RDF
    grainsize = y[y > 0][-1] * (r[y < 0][0]-r[y > 0][-1]) / \
        (y[y > 0][-1]-y[y < 0][0]) + r[y > 0][-1]

    plt.plot(r, y)
    plt.scatter(x=grainsize, y=0, label='Cluster {}, r={}'.format(
        cluster_id, round(grainsize, 3)))

plt.title("Grain Center Spatial Correlation", fontsize=16)
plt.xlabel(r"$r$", fontsize=14)
plt.ylabel(r"$\langle p \ast q \rangle (r)$", fontsize=14)
plt.tick_params(which="both", axis="both", labelsize=12)

This table shows the grain radii we expected compared to the ones we computed.

print('Actual Radius\tMeasured Radius')
for actual, measured in zip(sorted(grain_radii), sorted(measured_grain_radii)):
    print('{:.4f}\t\t{:.4f}'.format(actual, measured))
Actual Radius   Measured Radius
0.4628          0.4821
0.5246          0.5573
0.5345          0.5684
0.5457          0.5684
0.6608          0.6471
0.6705          0.6471
0.6825          0.6831
0.6879          0.6831
0.6939          0.6831
0.6953          0.7226
0.7101          0.7330
0.7298          0.7383
0.7510          0.7525
0.7731          0.7640
0.7952          0.7712
0.8106          0.8115