Messages

1

Support Forum / Re: Average Quaternions for Each Grain

« on: October 20, 2018, 01:06:43 AM »

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.

2

Support Forum / Re: Average Quaternions for Each Grain

« 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


3

Support Forum / Re: Average Quaternions for Each Grain

« 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

4

Support Forum / Re: Average Quaternions for Each Grain

« 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