Code examples

This page presents various example scripts in the following categories. Make sure to read the intro section of each category first:

Batch scripts

As described in Running scripts, batch scripts carry out sequences of program actions in a non-interactive manner and are typically executed from the terminal using the ovitos script interpreter. The following examples demonstrate how this can be used to automate various tasks and accomplish new things that cannot even be done with the graphical program version.

Example B1: 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.
pipeline = import_file("../data/CuZr_metallic_glass.dump.gz")

# Set atomic radii (required for polydisperse Voronoi tessellation).
atom_types = pipeline.source.compute().particles['Particle Type'].types
atom_types[0].radius = 1.35   # Cu atomic radius (atom type 1 in input file)
atom_types[1].radius = 1.55   # Zr atomic radius (atom type 2 in input file)

# Set up the Voronoi analysis modifier.
voro = VoronoiAnalysisModifier(
    compute_indices = True,
    use_radii = True,
    edge_count = 6, # Length after which Voronoi index vectors are truncated
    edge_threshold = 0.1
)
pipeline.modifiers.append(voro)
                      
# Let OVITO compute the results.
data = pipeline.compute()

# Make sure we did not lose information due to truncated Voronoi index vectors.
if data.attributes['Voronoi.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(data.attributes['Voronoi.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 = data.attributes['Voronoi.max_face_order']
    #   data = pipeline.compute()

# Access computed Voronoi indices.
# This is an (N) x (edge_count) array.
voro_indices = data.particles['Voronoi Index']

# 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[1])
    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 %)

Example B2: 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 BondsEnumerator 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.
pipeline = import_file("../data/NanocrystallinePd.dump.gz")

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

# Compute CNA indices on the basis of the created bonds.
pipeline.modifiers.append(CommonNeighborAnalysisModifier(
    mode = CommonNeighborAnalysisModifier.Mode.BondBased))
                      
# Let OVITO's data pipeline do the heavy work.
data = pipeline.compute()

# The 'CNA Indices' bond property is a a two-dimensional array 
# containing the three CNA indices computed for each bond in the system.
cna_indices = data.bonds['CNA Indices']

# 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[1])
    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 = BondsEnumerator(data)

# Loop over particles and print their CNA indices.
for particle_index in range(data.particles.count):
    
    # 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")

Example B3: 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 StaticSource instance. Finally, a Pipeline is created and the StaticSource is set as its data source.

import ovito
from ovito.data import DataCollection, SimulationCell, ParticleType
from ovito.pipeline import StaticSource, Pipeline

# Create the data collection (initially empty):
data = DataCollection()

# XYZ coordinates of the three atoms to create:
pos = [(1.0, 1.5, 0.3),
       (7.0, 4.2, 6.0),
       (5.0, 9.2, 8.0)]

# Create the particle position property:
pos_prop = data.particles.create_property('Position', data=pos)

# Create the particle type property and insert two atom types:
type_prop = data.particles.create_property('Particle Type')
type_prop.types.append(ParticleType(id = 1, name = 'Cu', color = (0.0,1.0,0.0)))
type_prop.types.append(ParticleType(id = 2, name = 'Ni', color = (0.0,0.5,1.0)))
with type_prop:
    type_prop[0] = 1  # First atom is Cu
    type_prop[1] = 2  # Second atom is Ni
    type_prop[2] = 2  # Third atom is Ni

# Create a user-defined particle property with some data:
my_data = [3.141, -1.2, 0.23]
my_prop = data.particles.create_property('My property', data=my_data)

# Create the simulation box:
cell = SimulationCell(pbc = (False, False, False))
with cell:
    cell[...] = [[10,0,0,0],
                 [0,10,0,0],
                 [0,0,10,0]]
cell.vis.line_width = 0.1
data.objects.append(cell)

# Create 3 bonds between particles:
bond_topology = [[0,1], [1,2], [2,0]]
data.bonds.create_property('Topology', data=bond_topology)

# Create a pipeline, set source and insert it into the scene:
pipeline = Pipeline(source = StaticSource())
pipeline.source.assign(data)
pipeline.add_to_scene()

User-defined modifier functions

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.

Example M1: Calculating mean square displacement

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 the 
    # CalculateDisplacementsModifier that precedes this user-defined modifier in the 
    # data pipeline:
    displacement_magnitudes = input.particles['Displacement']

    # 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 Table of Values 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 a data pipeline.
pipeline = import_file("input/simulation.dump")

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

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

    # Access the per-particle displacement magnitudes computed by the 
    # CalculateDisplacementsModifier that precedes this user-defined modifier in the 
    # data pipeline:
    displacement_magnitudes = input.particles['Displacement']

    # 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.
pipeline.modifiers.append(PythonScriptModifier(function = modify))

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

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 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 status text in the status bar:
    yield 'Calculating order parameters'

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

        # Loop over all particles.
        for i in range(len(order_params)):
            
            # Update progress indicator in the status bar
            yield (i/len(order_params))
            
            # 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.

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 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[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, input, output):
    """ The user-defined modifier function """
    
    # Input:
    orientations = input.particles['Orientation']
    
    # Output:
    output.particles.create_property('Color', data=quaternions_to_colors(orientations))

Example M4: Finding overlapping particles

This example shows how to write a user-defined modifier function that searches for pairs of particles whose distance of separation is within the specified cutoff distance. Then one of the two particles in the pair is selected by the modifier. Subsequently, the user may apply the DeleteSelectedModifier to remove these selected particles from the system and eliminate any potential overlaps among particles.

The modifier function below makes use of the CutoffNeighborFinder utility class, which allows finding neighboring particles that are within a certain range of a central particles. The modifier produces the standard output particle property Selection.

from ovito.data import CutoffNeighborFinder

# Control parameter:
overlap_distance = 2.5

# The user-defined modifier function:
def modify(frame, input, output):

    # Show this text in the status bar while the modifier function executes:
    yield "Selecting overlapping particles"

    # Create 'Selection' output particle property
    selection = output.particles.create_property('Selection')
    
    # Prepare neighbor finder
    finder = CutoffNeighborFinder(overlap_distance, input)
    
    # Request write access to the output property array
    with selection:
    
        # Iterate over all particles
        for index in range(len(selection)):
        
            # Progress bar update
            yield (index / len(selection))
            
            # Iterate over all closeby particles around the current center particle
            for neigh in finder.find(index):
                
                # Once we find a neighbor which hasn't been marked yet, 
                # mark the current center particle. This test is to ensure that we
                # always select only one of the particles in a close pair.
                if selection[neigh.index] == 0:
                    selection[index] = 1
                    break

User-defined overlay functions

OVITO allows you to implement your type of viewport overlay by writing a Python function that gets called every time a viewport image is being rendered.

Example O1: Scale bar

The following script renders a scale bar into the viewport (with a fixed length of 4 nm, as shown in the example picture). You can copy/paste the source code into the script input field and adjust the parameters in the code as needed.

../_images/python_script_scale_bar_example.png
from PyQt5.QtCore import *
from PyQt5.QtGui import *

# Parameters:
bar_length = 40   # Simulation units (e.g. Angstroms)
bar_color = QColor(0,0,0)
label_text = "{} nm".format(bar_length/10)
label_color = QColor(255,255,255)

# This function is called by OVITO on every viewport update.
def render(args):
    if args.is_perspective: 
        raise Exception("This overlay only works with non-perspective viewports.")
        
    # Compute length of bar in screen space
    screen_length = args.project_size((0,0,0), bar_length)

    # Define geometry of bar in screen space
    height = 0.07 * args.painter.window().height()
    margin = 0.02 * args.painter.window().height()
    rect = QRectF(margin, margin, screen_length, height)

    # Render bar rectangle
    args.painter.fillRect(rect, bar_color)

    # Render text label
    font = args.painter.font()
    font.setPixelSize(height)
    args.painter.setFont(font)
    args.painter.setPen(QPen(label_color))
    args.painter.drawText(rect, Qt.AlignCenter, label_text)

Example O2: Including data plots in rendered images

The following script demonstrates how to use the Matplotlib Python module to render a histogram on top the three-dimensional visualization. The histogram data is dynamically computed by a HistogramModifier in the data pipeline in this example.

../_images/python_script_scale_bar_example.png
from ovito.modifiers import HistogramModifier
import matplotlib
matplotlib.use('Agg') # Activate 'agg' backend for off-screen plotting.
import matplotlib.pyplot as plt
import PyQt5.QtGui

def render(args):
    # Look up the HistogramModifier in the pipeline whose data we will plot.
    hist_modifier = None
    for mod in args.scene.pipelines[0].modifiers:
        if isinstance(mod, HistogramModifier):
            hist_modifier = mod
            break
    if not hist_modifier: raise RuntimeError('Histogram modifier is not present')

    #  Compute plot size in inches (DPI determines label size)
    dpi = 80
    plot_width = 0.5 * args.size[0] / dpi
    plot_height = 0.5 * args.size[1] / dpi
    
    # Create matplotlib figure:
    fig, ax = plt.subplots(figsize=(plot_width,plot_height), dpi=dpi)
    fig.patch.set_alpha(0.5)
    plt.title('Coordination')
    
    # Plot histogram data
    ax.bar(hist_modifier.histogram[:,0], hist_modifier.histogram[:,1])
    plt.tight_layout()
    
    # Render figure to an in-memory buffer.
    buf = fig.canvas.print_to_buffer()
    
    # Create a QImage from the memory buffer
    res_x, res_y = buf[1]
    img = PyQt5.QtGui.QImage(buf[0], res_x, res_y, PyQt5.QtGui.QImage.Format_RGBA8888)
    
    # Paint QImage onto viewport canvas 
    args.painter.drawImage(0, 0, img)	

Example O3: Highlight a particle

../_images/python_script_highlight_example.png
from ovito.vis import ParticlesVis
from PyQt5.QtCore import *
from PyQt5.QtGui import *

def render(args):
    
    # Get output data collection of first scene pipeline.
    data = args.scene.pipelines[0].compute()
    positions = data.particles['Position']
    pindex = 0 # The index of the particle to be highlighted
    
    # Project center point of particle.
    xy = args.project_point(positions[pindex])
    if xy is None: return
    
    # Determine display radius of particle.
    if 'Radius' in data.particles:
        radius = data.particles['Radius'][pindex]
    elif 'Particle Type' in data.particles:
        particle_type = data.particles['Particle Type'][pindex]
        radius = data.particles['Particle Type'].type_by_id(particle_type).radius
    else:
        radius = node.get_vis(ParticlesVis).radius

    # Calculate screen-space size of particle in pixels.
    screen_radius = args.project_size(positions[pindex], radius)

    # Draw a dashed circle around the particle.
    pen = QPen(Qt.DashLine)
    pen.setWidth(3)
    pen.setColor(QColor(0,0,255))
    args.painter.setPen(pen)	
    args.painter.drawEllipse(QPointF(xy[0], xy[1]), screen_radius, screen_radius)
    
    # Draw an arrow pointing at the particle.
    arrow_shape = QPolygonF()
    arrow_shape.append(QPointF(0,0))
    arrow_shape.append(QPointF(10,10))
    arrow_shape.append(QPointF(10,5))
    arrow_shape.append(QPointF(40,5))
    arrow_shape.append(QPointF(40,-5))
    arrow_shape.append(QPointF(10,-5))
    arrow_shape.append(QPointF(10,-10))
    args.painter.setPen(QPen())
    args.painter.setBrush(QBrush(QColor(255,0,0)))
    args.painter.translate(QPointF(xy[0], xy[1]))
    args.painter.rotate(-45.0)
    args.painter.translate(QPointF(screen_radius,0))
    args.painter.scale(2,2)
    args.painter.drawPolygon(arrow_shape)