Example M3: Color mapping to visualize local lattice orientation

The Polyhedredral Template Matching (PTM) function of OVITO allows computing the local lattice orientation for each atom in a (poly)crystal. The modifier stores the computed local orientations as quaternions, i.e., rotations within the fundamental zone, in the particle property named Orientation. Each per-particle quaternion can be translated into an RGB color to visualize the local lattice orientation. This can be achieved by inserting a custom Python modifier into the pipeline, which translates the output of the PTM modifier into RGB values and stores them in the Color particle property.

In the graphical OVITO version, simply insert a new Python modifier and copy/paste the following script into the source code window:

from ovito.data import *
import math
import numpy as np

def quaternions_to_colors(qs):
    """ Takes a list of quaternions (Nx4 array) and returns a list of
        corresponding RGB colors (Nx3 array) """
 
    if len(qs.shape) != 2:
        raise RuntimeError("qs must be a 2-dimensional array")
 
    if qs.shape[1] != 4:
        raise RuntimeError("qs must be a n x 4 dimensional array")
 
    # Project quaternions into Rodrigues space: rs = (qs.X/qs.W, qs.Y/qs.W, qs.Z/qs.W)
    # Note that the qs.W may be zero for particles for which no lattice orientation
    # could be computed by the PTM modifier.
    rs = np.zeros_like(qs[:,:3])
    np.divide(qs[:,0], qs[:,3], out=rs[:,0], where=qs[:,3] != 0)
    np.divide(qs[:,1], qs[:,3], out=rs[:,1], where=qs[:,3] != 0)
    np.divide(qs[:,2], qs[:,3], out=rs[:,2], where=qs[:,3] != 0)

    # Compute vector lengths rr = norm(rs)
    rr = np.linalg.norm(rs, axis=1)
    rr = np.maximum(rr, 1e-9) # hack
 
    # Normalize Rodrigues vectors.
    rs[:,0] /= rr
    rs[:,1] /= rr
    rs[:,2] /= rr
 
    theta = 2 * np.arctan(rr)
    rs[:,0] *= theta
    rs[:,1] *= theta
    rs[:,2] *= theta
    
    # Normalize values.
    rs += math.radians(62.8)
    rs[:,0] /= 2*math.radians(62.8)
    rs[:,1] /= 2*math.radians(62.8)
    rs[:,2] /= 2*math.radians(62.8)
    
    return rs
    
def modify(frame, data):
    """ The user-defined modifier function """
    
    # Input:
    orientations = data.particles['Orientation']
    
    # Output:
    data.particles_.create_property('Color', data=quaternions_to_colors(orientations))