OVITO => Support Forum => Topic started by: mward2 on July 20, 2018, 06:33:31 PM

Title: Average Quaternions for Each Grain
Post by: mward2 on July 20, 2018, 06:33:31 PM
Hi all,

I want to get a characteristic quaternion for each grain. I am unsure how to go about this seeing as the orientations (W,X,Y,Z) that can be exported from Polyhedral Template Matching are not all similar within a grain. Any direction would be great.


Title: Re: Average Quaternions for Each Grain
Post by: Alexander Stukowski on July 23, 2018, 02:54:47 PM

Yes, the PTM function calculates the local lattice orientations, but averaging the corresponding quaternion values, which all get mapped to the fundamental zone, is non-trivial (the correct way depends on the lattice symmetry). And if you have a polycrystalline system, another challenge is to correctly identify the individual grains and assign the atoms to them first. Currently, Ovito doesn't have a built-in function for that, but this is an area of active research and, in fact, Peter Larsen, the author of the PTM function, is working on a powerful algorithm to perform the grain segmentation automatically. But we haven't made this function publicly available so far, because it's still very experimental. If you have an interesting application or test case for this new type of grain analysis function, please get in touch with me or Peter via email.

If you don't care about the grain segmentation problem, and simply want to average all the quaternions calculated by PTM for a known set of atoms, then it is certainly possible to write a small Python routine for Ovito that calculates the mean orientation value. Probably, Peter has already written some code for that in the past and I should be able to dig it out. Let me know if you want to take this path instead of waiting for a proper grain segmentation function to become available in Ovito in the future.

Title: Re: Average Quaternions for Each Grain
Post by: aadrych on October 04, 2018, 10:00:59 AM
Dear Dr Stukowski,

I was reading the forum for an answer to my question and I came across this thread.  I am also looking at trying to get an average orientation for a group of hcp atoms because I find that the orientation varies across the atoms.  I have tried plotting the orientations of all the atoms in the group to see if they cluster around a typical orientation but I find that even when looking at a group of atoms with the [0001] lattice orientation pointing parallel with the z axis of the reference frame I do not get the [0001] crystal orientation output via the quaternions.
So my questions are how do I get the quaternions to be output with reference to the hcp crystal structure? and how can I take an average of the grains without the requirement of the grain segmentation function?

Many thanks,
Title: Re: Average Quaternions for Each Grain
Post by: pmla on October 04, 2018, 04:58:45 PM
Hi Alex,

The HCP structure poses an additional challenge: each layer in stacking sequence has a different orientation relative to the hexagonal template used in PTM.
The issue is described very well in this article:
See page 6, where it is described as an 'undocumented bonus-feature' of PTM.

Here is a python script which will align the orientations output by PTM to that expected of a HCP crystal:

Code: [Select]
from ovito.data import *
import math
import numpy as np

hcp_gcs = np.array([
[             1,               0,               0,               0],
[           0.5,             0.5,             0.5,             0.5],
[           0.5,            -0.5,            -0.5,            -0.5],
[             0,   np.sqrt(2./3),  -np.sqrt(1./6),  -np.sqrt(1./6)],
[             0,   np.sqrt(1./6),  -np.sqrt(2./3),   np.sqrt(1./6)],
[             0,   np.sqrt(1./6),   np.sqrt(1./6),  -np.sqrt(2./3)],
[             0,  -np.sqrt(1/3.),  -np.sqrt(1/3.),  -np.sqrt(1/3.)],
[ np.sqrt(3./4),  np.sqrt(1./12),  np.sqrt(1./12),  np.sqrt(1./12)],
[ np.sqrt(3./4), -np.sqrt(1./12), -np.sqrt(1./12), -np.sqrt(1./12)],
[             0,  1 / np.sqrt(2), -1 / np.sqrt(2),               0],
[             0,               0,  1 / np.sqrt(2), -1 / np.sqrt(2)],
[             0,  1 / np.sqrt(2),              -0, -1 / np.sqrt(2)] ])

inv_hcp_gcs = np.copy(hcp_gcs)
inv_hcp_gcs[:,0] = -inv_hcp_gcs[:,0]

def quat_rot(r, a):

    b0 = r[0] * a[0] - r[1] * a[1] - r[2] * a[2] - r[3] * a[3]
    b1 = r[0] * a[1] + r[1] * a[0] + r[2] * a[3] - r[3] * a[2]
    b2 = r[0] * a[2] - r[1] * a[3] + r[2] * a[0] + r[3] * a[1]
    b3 = r[0] * a[3] + r[1] * a[2] - r[2] * a[1] + r[3] * a[0]
    return np.array([b0, b1, b2, b3])
def hcp_crystalline(q):
    index = np.argmax(np.abs(np.dot(inv_hcp_gcs, q)))
    return quat_rot(q, hcp_gcs[index])
def modify(frame, input, output):

    qs = input.particle_properties.orientation.array
    fundamental_property = output.create_particle_property(ParticleProperty.Type.Orientation)
    mapped = np.zeros_like(qs)

    for i, q in enumerate(qs):
        x, y, z, w = q
        qptm = [w, x, y, z]
        qf = hcp_crystalline(qptm)
        w, x, y, z = qf
        qf = np.array([x, y, z, w])
        mapped[i,:] = qf
    fundamental_property.marray[:] = mapped

The script stores the modified orientations in the 'fundamental_property' particle property.

This should result in a single cluster of orientations, at least for a perfect HCP crystal.  You may have further challenges if your orientations are spread out.

Title: Re: Average Quaternions for Each Grain
Post by: Constanze Kalcher on October 04, 2018, 05:00:57 PM
Thank you Peter!  :)
Title: How to get hcp orientation out of quartenions?
Post by: aadrych on October 05, 2018, 02:41:53 PM
Hi all,

I have been trying to get the lattice orientation of a group of atoms by collecting the quaternions output from the polyhedral template matching option in Ovito.  However when I loaded the quaternions for a crystal with the [0001] vector point parallel with the z reference frame and plotted it on a pole figure the orientation was not near the [0001] crystal direction.
I was wondering why this may be the case? Perhaps the quaternions are calculated from a different reference frame ? And if there was a way to find the average orientation for a group of atoms rather than an orientation per atom?

Any help would be much appreciated and thank you in advance.

Title: Re: How to get hcp orientation out of quartenions?
Post by: Alexander Stukowski on October 05, 2018, 04:49:08 PM
Hi Alex,

Please take a look at this updated version of the PTM doc page:


Peter Larsen (user ‘pmla’ in this forum) has included a table that precisely defines the reference orientation for each lattice type. I’m not sure, but I was convinced that for HCP the reference orientation is indeed defined such that [0001] is parallel to the z-axis.

A possible source of error (please check yourself) may be the representation of quaternions. OVITO employs a (x,y,z,w) representation, but other codes use a (w,x,y,z) representation.

If you want to calculate an average orientation for a crystal, you can do that by computing the arithmetic mean of the atomic quaternions as long as their orientations are not close to a boundary of the fundamental zone, where symmetry flips occur. Currently, OVITO doesn’t provide a function for computing an average orientation from the per-atom orientations which takes into account the lattice symmetries. I hope we can add something like this in the future.

Let me mention that the existing DXA function of OVITO also calculates mean orientations of crystallites as a side product of the dislocation identification. These grain orientations are currently not accessible via the public Python API though. The only way to obtain them is to export the DXA results to a .ca file. Here you will find the mean orientation matrices for each crystallite identified by the DXA.

Title: Re: Average Quaternions for Each Grain
Post by: aadrych on October 10, 2018, 09:44:04 AM
Hi Alex and Peter,

Thank you for posting replies so quickly.

I tried your python code Peter and I can now see that the colouring and the output quaternions are the same for all the atoms in bulk rather than varying dependent on the z-position of the atoms, which was more like I was expecting.  However the output quaternion for an atom in a bulk with the [0001] crystal orientation pointing parallel with the z-dimension comes out as [0.279848,-0.364705,-0.115917,0.880476] ( written as qx,qy,qz,qw ) which as a euler angle with a bunge notation (ZXZ rotation) should be 0,0,0,0 but instead the calculated euler angle is [-125, 45, -90 ].  I have checked that my conversion calcualtion from quartenion to euler uses the same input of qx, qy, qz and qw so I am unsure why the output quaternions do not come out as expected.

Is there no other rotation that has been applied to the quaternions ?

Many thanks in advance,

Title: Re: Average Quaternions for Each Grain
Post by: pmla on October 10, 2018, 09:55:51 PM
Hi Alex,

The output quaternions describe rotations of the reference templates.  Unfortunately, when I wrote PTM I did not select an intuitive alignment for the HCP template.

Nonetheless, you can still map the orientation into your preferred frame of reference using the reference templates.
The template coordinates are described here:

and are also in code form here:

Title: Re: Average Quaternions for Each Grain
Post by: aadrych on October 16, 2018, 08:31:13 AM
Hi Peter,

I'm sorry but I'm struggling to understand this. 

So the reference template isn't the hcp lattice but the fundamental zone?  I tried plotting the reference template coordinates from the online help page but this does not look like a hcp lattice.

I am unsure of how to align the reference template because I am not sure of what orientation the reference template is in ?

I have instead tried realigning the rotated atomic coordinates using the output quaternions, back to the atomic positions that I know lie with the [0001] orientation pointing parallel with the z-axis but as of yet I am still unsuccessful. 

Any guidance would be much appreciated to get the quaternions to show the correct orientation of the atoms.

Many thanks in advance.