Author Topic: Average Quaternions for Each Grain  (Read 1440 times)

mward2

  • Newbie
  • *
  • Posts: 1
Average Quaternions for Each Grain
« 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.

Thanks,

Madeline

Alexander Stukowski

  • Administrator
  • Hero Member
  • *****
  • Posts: 501
Re: Average Quaternions for Each Grain
« Reply #1 on: July 23, 2018, 02:54:47 PM »
Madeline,

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.

-Alex

aadrych

  • Newbie
  • *
  • Posts: 6
Re: Average Quaternions for Each Grain
« Reply #2 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,
Alex

pmla

  • Newbie
  • *
  • Posts: 4
Re: Average Quaternions for Each Grain
« Reply #3 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:
https://journals.aps.org/prmaterials/abstract/10.1103/PhysRevMaterials.2.083606
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.

Peter
« Last Edit: October 04, 2018, 05:13:26 PM by pmla »

Constanze Kalcher

  • Administrator
  • Full Member
  • *****
  • Posts: 113
Re: Average Quaternions for Each Grain
« Reply #4 on: October 04, 2018, 05:00:57 PM »
Thank you Peter!  :)

aadrych

  • Newbie
  • *
  • Posts: 6
How to get hcp orientation out of quartenions?
« Reply #5 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.

Alex

Alexander Stukowski

  • Administrator
  • Hero Member
  • *****
  • Posts: 501
Re: How to get hcp orientation out of quartenions?
« Reply #6 on: October 05, 2018, 04:49:08 PM »
Hi Alex,

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

http://www.ovito.org/manual_testing/particles.modifiers.polyhedral_template_matching.html

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.

-Alex
« Last Edit: October 05, 2018, 05:00:47 PM by Alexander Stukowski »

aadrych

  • Newbie
  • *
  • Posts: 6
Re: Average Quaternions for Each Grain
« Reply #7 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,

Alex

pmla

  • Newbie
  • *
  • Posts: 4
Re: Average Quaternions for Each Grain
« Reply #8 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:
http://www.ovito.org/manual_testing/particles.modifiers.polyhedral_template_matching.html

and are also in code form here:
https://gitlab.com/stuko/ovito/blob/master/src/3rdparty/ptm/ptm_constants.h

Peter

aadrych

  • Newbie
  • *
  • Posts: 6
Re: Average Quaternions for Each Grain
« Reply #9 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.

Alex

pmla

  • Newbie
  • *
  • Posts: 4
Re: Average Quaternions for Each Grain
« Reply #10 on: October 18, 2018, 07:35:55 PM »
Hi Alex,

Sorry for the slow reply.
For the HCP structure, the reference template contains the 12 nearest neighbours of a central atom.  It is in fact a small piece of the HCP lattice.
Anyway, I have made a change to the previous script.  The identity quaternion should now represent the orientation with [0001] parallel to the z-axis.
Hopefully this will solve your problem.

Peter

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)

    qfix = np.array([-0.880476239217, 0.279848142333, -0.364705199631, -0.115916895959])
    qfix /= np.linalg.norm(qfix)

    for i, q in enumerate(qs):
        x, y, z, w = q
        qptm = [w, x, y, z]
        qf = hcp_crystalline(qptm)
        qf = quat_rot(qf, qfix)

        w, x, y, z = qf
        qf = np.array([x, y, z, w])
        mapped[i,:] = qf
    fundamental_property.marray[:] = mapped

aadrych

  • Newbie
  • *
  • Posts: 6
Re: Average Quaternions for Each Grain
« Reply #11 on: October 19, 2018, 12:09:24 PM »
Hi Peter,

Thank you for looking in to this!
The quaternion for the [0001] crystal orientation with 0 misorientation around the z-axis now makes sense but unfortunately the quaternion output for any misorientation of the lattice does not come out as I would expect.

I have several lattices with the same crystal orientation i.e. the [0001] crystal direction is pointing parallel with the z-axis of the reference but I rotate the lattice around some angle theta around the z-axis.  If the quaternion mapping was working I would expect the quaternion to remain the same because the orientation of the crystal hasn't changed? I'm not sure how the quaternions map this to the reference structure though? Or whether it's because of the position of the neighbouring atoms?

Could you explain a bit more about how you chose the position of the neighbouring atoms?

I guess that if the lattice is rotated around the z axis the position relative to the neighbouring atoms would change and perhaps this checks as an orientation change?

Any help would be much appreciated.

Thanks in advance,
Alex   
 
« Last Edit: October 19, 2018, 03:00:49 PM by aadrych »

pmla

  • Newbie
  • *
  • Posts: 4
Re: Average Quaternions for Each Grain
« Reply #12 on: October 20, 2018, 01:06:43 AM »
Quote
I have several lattices with the same crystal orientation i.e. the [0001] crystal direction is pointing parallel with the z-axis of the reference but I rotate the lattice around some angle theta around the z-axis.  If the quaternion mapping was working I would expect the quaternion to remain the same because the orientation of the crystal hasn't changed?
If you have rotated the lattice, the orientation has changed.  You can think of the orientation as a rotation.
A quaternion is just a more convenient way of encoding a rotation/orientation.  You can read more about it here:
https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation

The section titled "Quaternion-derived rotation matrix" shows the relationship between quaternions and rotation matrices.

aadrych

  • Newbie
  • *
  • Posts: 6
Re: Average Quaternions for Each Grain
« Reply #13 on: October 23, 2018, 04:15:36 PM »
Hi Peter,

Ahh I see now.  I was using the wrong rotation matrix which didn't help me trying to visualise the misorientation of the matrix.

Thank you for helping!

Many thanks,
Alex