LAMMPS quarterions, Rowan and Freud

[ygk2zhang]

Hi,

I would like to make sure that there is no mistake. I have a system which looks like this.

I use Lammps to dump centre of mass of each hexagon and corresponding quarterions:
centre_mass x, y and z and quatw,quati,quatj,quatk = quaternion components for body
I convert the centre of mass xyz to points and quartw,I,j, k to rowanist which will be used to compute orientations.

However, the results look weird, so I want to know if I need to correct something.

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:
orientations = rowan.rotate(rowanist,ax)
nop = freud.order.Nematic(ax)
nop.compute(orientations)
print(f"For axis {ax}, the value of the order parameter is {nop.order:0.3f}. \n Director: {nematic.director}")

Results:
For axis [1, 0, 0], the value of the order parameter is 0.996.
Director: [ 0.02584443 -0.00645558 -0.9996453 ]
For axis [0, 1, 0], the value of the order parameter is 0.263.
Director: [ 0.02584443 -0.00645558 -0.9996453 ]
For axis [0, 0, 1], the value of the order parameter is 0.263.
Director: [ 0.02584443 -0.00645558 -0.9996453 ]
For axis [1, 1, 0], the value of the order parameter is 0.574.
Director: [ 0.02584443 -0.00645558 -0.9996453 ]
For axis [1, 0, 1], the value of the order parameter is 0.579.
Director: [ 0.02584443 -0.00645558 -0.9996453 ]
For axis [0, 1, 1], the value of the order parameter is 0.263.
Director: [ 0.02584443 -0.00645558 -0.9996453 ]
For axis [1, 1, 1], the value of the order parameter is 0.450.
Director: [ 0.02584443 -0.00645558 -0.9996453 ]

I don't see a strong preference of the system toward [1, 0,0] so makes me wonder I did something wrong.

Best Regards,

[DomFijan]

If you are unsure if the conventions for quaternions are the same between rowan and Lammps, you can check it in rowans docs here . In particular it says:

Quaternions are encoded as arrays of shape (..., 4), with the convention that the final dimension of an array ``(a, b, c, d) represents the quaternion a + bi + cj + dk.

It seems you are using the old version of Freud. For old API the compute expects QUATERNIONS not orientations.

If you use the new freud nematic API (v3.0+), the Nematic constructor no longer requires a "director". Instead you can just provide the orientations to the compute.

[ygk2zhang]

Hi DomFijan

Thank you for your response. I am using the new Freud hence I modified the code a little.
In the old one:
the same code will look like this

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(rowanist)
print(f"For axis {ax}, the value of the order parameter is {nop.order:0.3f}.")

[ygk2zhang]

Hi. I know what is the problem is. I am not sure how to fix it.

Here is the code:

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

%matplotlib inline

orientations = rowan.rotate(rowanist, [1, 0, 0])
nematic = freud.order.Nematic()
nematic.compute(orientations)
# make a 3D plot of the system
fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(111, projection="3d")
ax.set_title(f"Nematic order parameter: {nematic.order:.3f} \n Nem Director: {nematic.director}")
ax.quiver3D(
    pointss[:, 0],
    pointss[:, 1],
    pointss[:, 2],
    orientations[:, 0],
    orientations[:, 1],
    orientations[:, 2],
    length=60,
    normalize=False,
    color="black",
)
ax.quiver3D(
    -L / 2 * nematic.order,
    -L / 2 * 1.2,
    -L / 2,
    nematic.director[0],
    nematic.director[1],
    nematic.director[2],
    length=60,
    normalize=False,
    color="blue",
)

I am not sure why [1 0 0] vector is along the z-axis.

[DomFijan]

The [1,0,0] vector should be the starting, reference orientation of your shapes you simulated in lammps. The quaternion defines the rotation of this starting orientation into a new orientation (presumably at a given simulation frame/snapshot) for each particle. The nematic director has nothing to do with this starting orientation of your particles.

The computed nematic director is [0,0,-1] - blue arrow, which is along the z axis. Based on the figure, your orientations seem to be pointing along the z axis, so the result makes sense. Two things that might have gone wrong here is the starting orientation of your particles ([1,0,0]), and computation of the orientations. You should go back to your simulation script and double check what your starting orientation was, and to confirm that the orientations that you are getting with rowan.rotate are correct. It is possible that you have to use a different quaternion convention.

Another issue here might be that nematic order might not be the best order parameter to use here. There are many equivalent rotations in hexagons, so these should be somehow taken into consideration, and are not taken into consideration by nematic order parameter. Perhaps something in the environment module would work better?

[ygk2zhang]

Thank you very much for your help. Now, the results make sense to me. I will look into the environment module.