freud.order.Nematic#

Nematic Order Parameter#

The freud.order module provids the tools to calculate various order parameters that can be used to identify phase transitions. This notebook demonstrates the nematic order parameter, which can be used to identify systems with strong orientational ordering but no translational ordering. For this example, we’ll start with a set of random positions in a 3D system, each with a fixed, assigned orientation. Then, we will show how deviations from these orientations are exhibited in the order parameter.

[1]:
import freud
import matplotlib.pyplot as plt
import numpy as np
import rowan  # for quaternion math, see rowan.readthedocs.io for more information.
from mpl_toolkits.mplot3d import Axes3D

In order to work with orientations in freud, we need to do some math with quaternions. If you are unfamiliar with quaternions, you can read more about their definition and how they can be used to represent rotations. For the purpose of this tutorial, just consider them as 4D vectors, and know that the set of normalized (i.e. unit norm) 4D vectors can be used to represent rotations in 3D. In fact, there is a 1-1 mapping between normalize quaternions and 3x3 rotation matrices. Quaternions are more computationally convenient, however, because they only require storing 4 numbers rather than 9, and they can be much more easily chained together. For our purposes, you can largely ignore the contents of the next cell, other than to note that this is how we perform rotations of vectors using quaternions instead of matrices.

[2]:
# Random positions are fine for this. Order is measured
# in terms of similarity of orientations, not positions.
L = 10
N = 100
box, points = freud.data.make_random_system(L, N, seed=0)
orientations = np.array([[1, 0, 0, 0]] * N)
[3]:
# To show orientations, we use arrows rotated by the quaternions.
arrowheads = rowan.rotate(orientations, [1, 0, 0])

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.quiver3D(
    points[:, 0],
    points[:, 1],
    points[:, 2],
    arrowheads[:, 0],
    arrowheads[:, 1],
    arrowheads[:, 2],
)
ax.set_title("Orientations", fontsize=16);
../../../_images/gettingstarted_examples_module_intros_order.Nematic_4_0.png

The nematic order parameter provides a measure of how much of the system is aligned with respect to some provided reference vector. As a result, we can now compute the order parameter for a few simple cases. Since our original system is oriented along the x-axis, we can immediately test for that, as well as orientation along any of the other coordinate axes.

[4]:
nop = freud.order.Nematic([1, 0, 0])
nop.compute(orientations)
print(f"The value of the order parameter is {nop.order}.")
The value of the order parameter is 1.0.

In general, the nematic order parameter is defined as the eigenvalue corresponding to the largest eigenvector of the nematic tensor, which is also computed by this class and provides an average over the orientations of all particles in the system. As a result, we can also look at the intermediate results of our calculation and see how they are related. To do so, let’s consider a more interesting system with random orientations.

[5]:
# We rotate identity quaternions slightly, in a random direction
np.random.seed(0)
interpolate_amount = 0.3
identity_quats = np.array([[1, 0, 0, 0]] * N)
orientations = rowan.interpolate.slerp(
    identity_quats, rowan.random.rand(N), interpolate_amount
)
[6]:
# To show orientations, we use arrows rotated by the quaternions.
arrowheads = rowan.rotate(orientations, [1, 0, 0])

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.quiver3D(
    points[:, 0],
    points[:, 1],
    points[:, 2],
    arrowheads[:, 0],
    arrowheads[:, 1],
    arrowheads[:, 2],
)
ax.set_title("Orientations", fontsize=16);
../../../_images/gettingstarted_examples_module_intros_order.Nematic_9_0.png

First, we see that for this nontrivial system the order parameter now depends on the choice of director.

[7]:
axes = [[1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0], [1, 0, 1], [0, 1, 1], [1, 1, 1]]
for ax in axes:
    nop = freud.order.Nematic(ax)
    nop.compute(orientations)
    print(f"For axis {ax}, the value of the order parameter is {nop.order:0.3f}.")
For axis [1, 0, 0], the value of the order parameter is 0.600.
For axis [0, 1, 0], the value of the order parameter is 0.586.
For axis [0, 0, 1], the value of the order parameter is 0.587.
For axis [1, 1, 0], the value of the order parameter is 0.591.
For axis [1, 0, 1], the value of the order parameter is 0.589.
For axis [0, 1, 1], the value of the order parameter is 0.573.
For axis [1, 1, 1], the value of the order parameter is 0.578.

Furthermore, increasing the amount of variance in the orientations depresses the value of the order parameter even further.

[8]:
interpolate_amount = 0.4
orientations = rowan.interpolate.slerp(
    identity_quats, rowan.random.rand(N), interpolate_amount
)

arrowheads = rowan.rotate(orientations, [1, 0, 0])
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.quiver3D(
    points[:, 0],
    points[:, 1],
    points[:, 2],
    arrowheads[:, 0],
    arrowheads[:, 1],
    arrowheads[:, 2],
)
ax.set_title("Orientations", fontsize=16)

axes = [[1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0], [1, 0, 1], [0, 1, 1], [1, 1, 1]]
for ax in axes:
    nop = freud.order.Nematic(ax)
    nop.compute(orientations)
    print(f"For axis {ax}, the value of the order parameter is {nop.order:0.3f}.")
For axis [1, 0, 0], the value of the order parameter is 0.451.
For axis [0, 1, 0], the value of the order parameter is 0.351.
For axis [0, 0, 1], the value of the order parameter is 0.342.
For axis [1, 1, 0], the value of the order parameter is 0.374.
For axis [1, 0, 1], the value of the order parameter is 0.391.
For axis [0, 1, 1], the value of the order parameter is 0.316.
For axis [1, 1, 1], the value of the order parameter is 0.344.
../../../_images/gettingstarted_examples_module_intros_order.Nematic_13_1.png

Finally, we can look at the per-particle quantities and build them up to get the actual value of the order parameter.

[9]:
# The per-particle values averaged give the nematic tensor
print(np.allclose(np.mean(nop.particle_tensor, axis=0), nop.nematic_tensor))
print("The nematic tensor:")
print(nop.nematic_tensor)

eig = np.linalg.eig(nop.nematic_tensor)
print("The eigenvalues of the nematic tensor:")
print(eig[0])
print("The eigenvectors of the nematic tensor:")
print(eig[1])

# The largest eigenvalue
print(
    "The largest eigenvalue, {:0.3f}, is equal to the order parameter {:0.3f}.".format(
        np.max(eig[0]), nop.order
    )
)
True
The nematic tensor:
[[ 0.0115407   0.21569438  0.14729623]
 [ 0.21569438  0.02040018  0.14309749]
 [ 0.14729623  0.14309748 -0.03194092]]
The eigenvalues of the nematic tensor:
[ 0.34387365 -0.20013455 -0.14373913]
The eigenvectors of the nematic tensor:
[[ 0.6173224   0.73592573 -0.27807635]
 [ 0.6237324  -0.6732561  -0.3970945 ]
 [ 0.47944868 -0.07169023  0.87463677]]
The largest eigenvalue, 0.344, is equal to the order parameter 0.344.