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

`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