2
« 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
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