Examples¶

This page provides a collection of example scripts:

Computing Voronoi indices¶

This script demonstrates the use of the Voronoi analysis modifier. The script calculates the distribution of Voronoi coordination polyhedra in an amorphous structure.

A Voronoi polyhedron is expressed in terms of the Schlaefli notation, which is a vector of indices (n1, n2, n3, n4, n5, n6, ...), where ni is the number of polyhedron faces with i edges/vertices.

The script computes the distribution of these Voronoi index vectors and lists the 10 most frequent polyhedron types in the dataset. In the case of a Cu64%-Zr36% bulk metallic glass, the most frequent polyhedron type is the icosahedron. It has 12 faces with five edges each. Thus, the corresponding Voronoi index vector is:

(0, 0, 0, 0, 12, 0, ...)

Python script:

# Import OVITO modules.
from ovito.io import *
from ovito.modifiers import *

# Import NumPy module.
import numpy

# Load a simulation snapshot of a Cu-Zr metallic glass.
node = import_file("../data/CuZr_metallic_glass.dump.gz")

# Set atomic radii (required for polydisperse Voronoi tessellation).
atypes = node.source.particle_properties.particle_type.type_list
atypes.radius = 1.35        # Cu atomic radius (atom type 1 in input file)
atypes.radius = 1.55        # Zr atomic radius (atom type 2 in input file)

# Set up the Voronoi analysis modifier.
voro = VoronoiAnalysisModifier(
compute_indices = True,
edge_count = 6, # Length after which Voronoi index vectors are truncated
edge_threshold = 0.1
)
node.modifiers.append(voro)

# Let OVITO compute the results.
node.compute()

# Make sure we did not lose information due to truncated Voronoi index vectors.
if voro.max_face_order > voro.edge_count:
print("Warning: Maximum face order in Voronoi tessellation is {0}, "
"but computed Voronoi indices are truncated after {1} entries. "
"You should consider increasing the 'edge_count' parameter to {0}."
.format(voro.max_face_order, voro.edge_count))
# Note that it would be possible to automatically increase the 'edge_count'
# parameter to 'max_face_order' here and recompute the Voronoi tessellation:
#   voro.edge_count = voro.max_face_order
#   node.compute()

# Access computed Voronoi indices as NumPy array.
# This is an (N)x(edge_count) array.
voro_indices = node.output.particle_properties['Voronoi Index'].array

# This helper function takes a two-dimensional array and computes a frequency
# histogram of the data rows using some NumPy magic.
# It returns two arrays (of equal length):
#    1. The list of unique data rows from the input array
#    2. The number of occurences of each unique row
# Both arrays are sorted in descending order such that the most frequent rows
# are listed first.
def row_histogram(a):
ca = numpy.ascontiguousarray(a).view([('', a.dtype)] * a.shape)
unique, indices, inverse = numpy.unique(ca, return_index=True, return_inverse=True)
counts = numpy.bincount(inverse)
sort_indices = numpy.argsort(counts)[::-1]
return (a[indices[sort_indices]], counts[sort_indices])

# Compute frequency histogram.
unique_indices, counts = row_histogram(voro_indices)

# Print the ten most frequent histogram entries.
for i in range(10):
print("%s\t%i\t(%.1f %%)" % (tuple(unique_indices[i]),
counts[i],
100.0*float(counts[i])/len(voro_indices)))

Program output:

(0, 0, 0, 0, 12, 0)     12274   (11.4 %)
(0, 0, 0, 2, 8, 2)      7485    (6.9 %)
(0, 0, 0, 3, 6, 4)      5637    (5.2 %)
(0, 0, 0, 1, 10, 2)     4857    (4.5 %)
(0, 0, 0, 3, 6, 3)      3415    (3.2 %)
(0, 0, 0, 2, 8, 1)      2927    (2.7 %)
(0, 0, 0, 1, 10, 5)     2900    (2.7 %)
(0, 0, 0, 1, 10, 4)     2068    (1.9 %)
(0, 0, 0, 2, 8, 6)      2063    (1.9 %)
(0, 0, 0, 2, 8, 5)      1662    (1.5 %)

Computing CNA bond indices¶

The following script demonstrates how to use the CreateBondsModifier to create bonds between particles. The structural environment of each created bond is then characterized with the help of the CommonNeighborAnalysisModifier, which computes a triplet of indices for each bond from the topology of the surrounding bond network. The script accesses the computed CNA bond indices in the output DataCollection of the modification pipeline and exports them to a text file. The script enumerates the bonds of each particle using the Bonds.Enumerator helper class.

The generated text file has the following format:

Atom    CNA_pair_type:Number_of_such_pairs ...

1       [4 2 1]:2  [4 2 2]:1 [5 4 3]:1
2       ...
...

Python script:

# Import OVITO modules.
from ovito.io import *
from ovito.modifiers import *
from ovito.data import *

# Import standard Python and NumPy modules.
import sys
import numpy

# Load the simulation dataset to be analyzed.
node = import_file("../data/NanocrystallinePd.dump.gz")

# Create bonds.
node.modifiers.append(CreateBondsModifier(cutoff = 3.5))

# Compute CNA indices on the basis of the created bonds.
node.modifiers.append(
CommonNeighborAnalysisModifier(mode = CommonNeighborAnalysisModifier.Mode.BondBased))

# Let OVITO's data pipeline do the heavy work.
node.compute()

# A two-dimensional array containing the three CNA indices
# computed for each bond in the system.
cna_indices = node.output.bond_properties['CNA Indices'].array

# This helper function takes a two-dimensional array and computes the frequency
# histogram of the data rows using some NumPy magic.
# It returns two arrays (of same length):
#    1. The list of unique data rows from the input array
#    2. The number of occurences of each unique row
def row_histogram(a):
ca = numpy.ascontiguousarray(a).view([('', a.dtype)] * a.shape)
unique, indices, inverse = numpy.unique(ca, return_index=True, return_inverse=True)
counts = numpy.bincount(inverse)
return (a[indices], counts)

# Used below for enumerating the bonds of each particle:
bond_enumerator = Bonds.Enumerator(node.output.bonds)

# Loop over particles and print their CNA indices.
for particle_index in range(node.output.number_of_particles):

# Print particle index (1-based).
sys.stdout.write("%i " % (particle_index+1))

# Create local list with CNA indices of the bonds of the current particle.
bond_index_list = list(bond_enumerator.bonds_of_particle(particle_index))
local_cna_indices = cna_indices[bond_index_list]

# Count how often each type of CNA triplet occurred.
unique_triplets, triplet_counts = row_histogram(local_cna_indices)

# Print list of triplets with their respective counts.
for triplet, count in zip(unique_triplets, triplet_counts):
sys.stdout.write("%s:%i " % (triplet, count))

# End of particle line
sys.stdout.write("\n")

Writing a custom modifier for calculating the mean square displacement¶

OVITO allows you to implement your own type of analysis modifier by writing a Python function that gets called every time the data pipeline is evaluated. This user-defined function has access to the positions and other properties of particles and can output information and results as new properties or global attributes.

As a first simple example, we look at the calculation of the mean square displacement (MSD) in a system of moving particles. OVITO already provides the built-in Displacement Vectors modifier, which calculates the displacement of every particle. It stores its results in the "Displacement Magnitude" particle property. So all our custom analysis modifier needs to do is to sum up the squared displacement magnitudes and divide by the number of particles:

def modify(frame, input, output):

# Access the per-particle displacement magnitudes computed by an existing
# Displacement Vectors modifier that precedes this custom modifier in the
# data pipeline:
displacement_magnitudes = input.particle_properties.displacement_magnitude.array

# Compute MSD:
msd = numpy.sum(displacement_magnitudes ** 2) / len(displacement_magnitudes)

# Output MSD value as a global attribute:
output.attributes["MSD"] = msd

When used within the graphical program, the MSD value computed by this custom modifier may be exported to a text file as a function of simulation time using OVITO’s standard file export feature (Select Calculation Results Text File as output format).

Alternatively, we can make use of the custom modifier from within a non-interactive batch script, which is executed by the ovitos interpreter. Then we have to insert the CalculateDisplacementsModifier programmatically:

from ovito.io import import_file, export_file
from ovito.modifiers import PythonScriptModifier, CalculateDisplacementsModifier
import numpy

# Load input data and create an ObjectNode with a data pipeline.
node = import_file("simulation.dump", multiple_frames = True)

# Calculate per-particle displacements with respect to initial simulation frame:
dmod = CalculateDisplacementsModifier()
node.modifiers.append(dmod)

# Define the custom modifier function:
def modify(frame, input, output):

# Access the per-particle displacement magnitudes computed by an existing
# Displacement Vectors modifier that precedes this custom modifier in the
# data pipeline:
displacement_magnitudes = input.particle_properties.displacement_magnitude.array

# Compute MSD:
msd = numpy.sum(displacement_magnitudes ** 2) / len(displacement_magnitudes)

# Output MSD value as a global attribute:
output.attributes["MSD"] = msd

# Insert custom modifier into the data pipeline.
node.modifiers.append(PythonScriptModifier(function = modify))

# Export calculated MSD value to a text file and let OVITO's data pipeline do the rest:
export_file(node, "msd_data.txt",
format = "txt",
columns = ["Timestep", "MSD"],
multiple_frames = True)

Implementing an advanced analysis modifier¶

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 rfcc). The sum of the distances over the 12 neighbors, phi = 1/12*sum(| ri - rfcc |), 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, input, output):

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

# Create output property.
order_param = output.create_user_particle_property(
"Order Parameter", "float").marray

# Prepare neighbor lists.
neigh_finder = NearestNeighborFinder(num_neighbors, input)

# Loop over all input particles
nparticles = input.number_of_particles
for i in range(nparticles):

# Update progress percentage indicator
yield (i/nparticles)

oparam = 0.0	# The order parameter of the current atom

# 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 particle property.
order_param[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.

Creating particles and bonds programmatically¶

The following script demonstrates how to create particles, a simulation cell, and bonds on the fly without loading an external simulation file. This approach can be used to implement custom data importers or dynamically generate atomic structures, which can then be further processed with OVITO or exported to a file.

The script creates different data objects and adds them to a new DataCollection instance. Finally, an ObjectNode is created and the DataCollection is set as its data source.

import ovito
from ovito.data import *

# The number of particles we are going to create.
num_particles = 3

# Create the particle position property.
pos_prop = ParticleProperty.create(ParticleProperty.Type.Position, num_particles)
pos_prop.marray = (1.0, 1.5, 0.3)
pos_prop.marray = (7.0, 4.2, 6.0)
pos_prop.marray = (5.0, 9.2, 8.0)

# Create the particle type property and insert two atom types.
type_prop = ParticleProperty.create(ParticleProperty.Type.ParticleType, num_particles)
type_prop.type_list.append(ParticleType(id = 1, name = 'Cu', color = (0.0,1.0,0.0)))
type_prop.type_list.append(ParticleType(id = 2, name = 'Ni', color = (0.0,0.5,1.0)))
type_prop.marray = 1  # First atom is Cu
type_prop.marray = 2  # Second atom is Ni
type_prop.marray = 2  # Third atom is Ni

# Create a user-defined particle property.
my_prop = ParticleProperty.create_user('My property', 'float', num_particles)
my_prop.marray = 3.141
my_prop.marray = 0.0
my_prop.marray = 0.0

# Create the simulation box.
cell = SimulationCell()
cell.matrix = [[10,0,0,0],
[0,10,0,0],
[0,0,10,0]]
cell.pbc = (True, True, True)
cell.display.line_width = 0.1

# Create bonds between particles.
bonds = Bonds()
bonds.add_full(0, 1)    # Creates two half bonds 0->1 and 1->0.
bonds.add_full(1, 2)    # Creates two half bonds 1->2 and 2->1.
bonds.add_full(2, 0)    # Creates two half bonds 2->0 and 0->2.

# Create a data collection to hold the particle properties, bonds, and simulation cell.
data = DataCollection()

# Create a node and insert it into the scene.
node = ovito.ObjectNode()
node.source = data

# Select the new node and adjust cameras of all viewports to show it.
ovito.dataset.selected_node = node
for vp in ovito.dataset.viewports:
vp.zoom_all()

Visualizing local lattice orientations using particle coloring¶

The Polyhedredral Template Matching (PTM) function of OVITO allows computing the local lattice orientation for each atom in a (poly)crystal. The computed local orientations are stored by the modifier as quaternions, i.e. as 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 != 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.