Example M2: Custom order parameter calculation

In the paper [Phys. Rev. Lett. 86, 5530] an order parameter is specified as a means of labeling an atom in the simulation as belonging to either the liquid or solid fcc crystal phase. In the following we will develop a custom analysis modifier for OVITO, which calculates this per-atom order parameter.

The order parameter is defined as follows (see the paper for details): For any of the 12 nearest neighbors of a given atom one can compute the distance the neighbor makes from the ideal fcc positions of the crystal in the given orientation (denoted by vector r_fcc). The sum of the distances over the 12 neighbors, phi = 1/12*sum(| r_i - r_fcc |), acts as an “order parameter” for the central atom.

Calculating this parameter involves finding the 12 nearest neighbors of each atom and, for each of these neighbors, determining the closest ideal lattice vector. To find the neighbors, OVITO provides the NearestNeighborFinder utility class. It directly provides the vectors from the central atom to its nearest neighbors.

Let us start by defining some inputs for the order parameter calculation at the global scope:

from ovito.data import NearestNeighborFinder
import numpy as np

# The lattice constant of the FCC crystal:
lattice_parameter = 3.6 

# The list of <110> ideal neighbor vectors of the reference lattice (FCC):
reference_vectors = np.asarray([
    (0.5, 0.5, 0.0),
    (-0.5, 0.5, 0.0),
    (0.5, -0.5, 0.0),
    (-0.5, -0.5, 0.0),
    (0.0, 0.5, 0.5),
    (0.0, -0.5, 0.5),
    (0.0, 0.5, -0.5),
    (0.0, -0.5, -0.5),
    (0.5, 0.0, 0.5),
    (-0.5, 0.0, 0.5),
    (0.5, 0.0, -0.5),
    (-0.5, 0.0, -0.5)
])
# Rescale ideal lattice vectors with lattice constant.
reference_vectors *= lattice_parameter

# The number of neighbors to take into account per atom:
num_neighbors = len(reference_vectors)

The actual modifier function needs to create an output particle property, which will store the calculated order parameter of each atom. Two nested loops run over all input atoms and their 12 nearest neighbors respectively.

def modify(frame, data):

    # Show a status text in the status bar:
    yield 'Calculating order parameters'

    # Create output particle property.
    order_params = data.particles_.create_property(
        'Order Parameter', dtype=float, components=1)
    
    # Prepare neighbor lists.
    neigh_finder = NearestNeighborFinder(num_neighbors, data)
    
    # Request write access to the output property array.
    with order_params:

        # Loop over all particles.
        for i in range(data.particles.count):
            
            # Update progress indicator in the status bar
            yield (i/data.particles.count)
            
            # Stores the order parameter of the current atom
            oparam = 0.0	
            
            # Loop over neighbors of current atom.
            for neigh in neigh_finder.find(i):
                
                # Compute squared deviation of neighbor vector from every 
                # reference vector.
                squared_deviations = np.linalg.norm(
                    reference_vectors - neigh.delta, axis=1) ** 2
                
                # Sum up the contribution from the best-matching vector.
                oparam += np.min(squared_deviations)

            # Store result in output array.
            order_params[i] = oparam / num_neighbors		

Note that the yield statements in the modifier function above are only needed to support progress feedback in the graphical version of OVITO and to give the pipeline system the possibility to interrupt the long-running calculation when needed.