Potential issue with view_orientation in freud.diffraction.DiffractionPattern.compute

[jennyfothergill]

Describe the bug
While testing out some diffraction utilities I noticed that I wasn't getting the symmetry I expect for certain view orientations. (Please correct me if I am completely misunderstanding how the view_orientation kwarg should be used!)

To Reproduce
(The gsd file I am using can be found here.)
I want to calculate the diffraction pattern of a simple cubic structure (see gsd file) from the face, edge, and corner views. I am using the following code to calculate the rotation quaternion from one vector to another:

def q_from_vectors(b, a=np.array([0,0,1])):                              
    """Calculate the quaternion representing the rotation from a to b."""
    q = np.empty(4)                                                      
    q[:3] = np.cross(a,b)                                                
    q[3] = np.dot(a,b) + np.sqrt(np.linalg.norm(a)**2 * np.linalg.norm(b)**2)                
    q /= np.linalg.norm(q)                                                         
    return q     

and to plot the diffraction pattern I am using this code

gsdfile = "sc10.gsd"

with gsd.hoomd.open(gsdfile) as f:
    snap = f[-1]
points = snap.particles.position
box = freud.Box.from_box(snap.configuration.box)
dp = freud.diffraction.DiffractionPattern(grid_size=1024, output_size=1024)

q = q_from_vectors(np.array([1,1,1]))
dp.compute((box, points), view_orientation=q)
dp.plot();

(q=[-0.32505758 0.32505758 0. 0.88807383])
I thought the rotation quaternion from [0,0,1] to [1,1,1] should result in a corner-on view which should have six-fold symmetry, but I don't see that, instead I see:

However, if I change the quaternion calculation to ignore the vector lengths:

def q_from_vectors(b, a=np.array([0,0,1])):                              
    """Calculate the quaternion representing the rotation from a to b."""
    q = np.empty(4)                                                      
    q[:3] = np.cross(a,b)                                                
    q[3] = np.dot(a,b)                
    q /= np.linalg.norm(q)                                                         
    return q     

(q=[-0.57735027 0.57735027 0. 0.57735027])
I do see six-fold symmetry:

To my knowledge this is the quaternion for double the required rotation? Please forgive me if I am just making a silly mistake...

System configuration (please complete the following information):

• OS: macOS
• Version of Python: 3.7.10
• Version of freud: 2.5.1

[bdice]

@jennyfothergill This is a long reply with many separate ideas. Apologies in advance for the brain dump.

Possible cause: quaternion math is tricky

From a quick glance, this might just be an issue with the conventions for quaternions. Some quaternion conventions use the format [x, y, z, w] and some use [w, x, y, z].

I recommend using rowan for quaternion calculations, as shown below. I think rowan's conventions are the most common conventions in our field, and all the software from our group relies on rowan for quaternion math.

rowan docs:

• https://rowan.readthedocs.io/en/latest/package-rowan.html#rowan.rotate
• https://rowan.readthedocs.io/en/latest/package-rowan.html#rowan.vector_vector_rotation

I edited the API for your snippet a little but I left the implementation untouched.

import numpy as np
import rowan


def q_from_vectors(a, b):
    """Calculate the quaternion representing the rotation from a to b."""
    q = np.empty(4)
    q[:3] = np.cross(a, b)
    q[3] = np.dot(a, b) + np.sqrt(np.linalg.norm(a) ** 2 * np.linalg.norm(b) ** 2)
    q /= np.linalg.norm(q)
    return q


target_axis = np.array([1, 1, 1])
unit_z = np.array([0, 0, 1])

q_jenny = q_from_vectors(target_axis, unit_z)

q_rowan = rowan.vector_vector_rotation(target_axis, unit_z)

np.set_printoptions(precision=5, suppress=True)
print("Jenny's quaternion:", q_jenny)
print("rowan's quaternion:", q_rowan)

print("Jenny's view axis:", rowan.rotate(q_jenny, target_axis))
print("rowan's view axis:", rowan.rotate(q_rowan, target_axis))  # should be aligned with z

Output:

Jenny's quaternion: [ 0.32506 -0.32506  0.       0.88807]
rowan's quaternion: [0.      0.32506 0.32506 0.88807]
Jenny's view axis: [-1.73205  0.       0.     ]
rowan's view axis: [0.      0.      1.73205]

(Fun fact: my spell checker wants to correct "quaternions" to "consternation." Seems accurate.)

freud DiffractionPattern

As for the diffraction pattern itself, we know there are some issue with freud's diffraction code that we're actively working on with @AKerr9509 this summer (#777, #723, #789, #750). It's possible that the DiffractionPattern docs could be improved by including a snippet demonstrating how to generate the quaternion you're seeking with rowan.

I built the latest master version of freud, which includes #777 (merged today) and is newer than what is present in v2.5.1. Here is a code snippet and the results I get, which are closer to your expectation:

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

gsdfile = "sc10.gsd"

with gsd.hoomd.open(gsdfile) as f:
    snap = f[-1]
points = snap.particles.position
box = freud.Box.from_box(snap.configuration.box)
dp = freud.diffraction.DiffractionPattern(grid_size=1024, output_size=1024)

target_axis = np.array([1, 1, 1])
unit_z = np.array([0, 0, 1])
q = rowan.vector_vector_rotation(target_axis, unit_z)
dp.compute((box, points), view_orientation=q)

fig, ax = plt.subplots(figsize=(5, 5), dpi=150)
dp.plot(ax=ax)
plt.show()

The rotation of the diffraction pattern is governed by rowan's choice of quaternion for the vector-vector rotation, which is arbitrary (vector-vector rotation has a degree of freedom that is undefined).

[jennyfothergill]

@bdice Thank you SO MUCH for the detailed reply. :) I wasn't aware of the differences in quaternion convention--using rowan will make this much simpler going forward. I agree that an example code snippet showing how to calculate the quaternion for the view_orientation kwarg would be very helpful to users (me, haha).
Again, thanks!

[bdice]

@jennyfothergill Glad I could help! 😄