Messages

1

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

2

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