# 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 Scripting usage, 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.data.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,
edge_threshold = 0.1
)
pipeline.modifiers.append(voro)

# Let OVITO compute the results.
data = pipeline.compute()

# Access computed Voronoi indices.
# This is an (N) x (M) array, where M is the maximum face order.
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.particles.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.particles.bonds)

# 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 the creation of particles, a simulation cell, and bonds on the fly without loading them from an external simulation file. This approach can be used to implement custom data importers or dynamically generate atomic structures within OVITO, which can then be further processed or exported to a file.

The script creates different data objects and adds them to a new `DataCollection`. Finally, a `Pipeline` is created and a `StaticSource` object is used to make the `DataCollection` its data source.

```from ovito.data import *
from ovito.pipeline import *

# Create the data collection containing a Particles object:
particles = Particles()
data = DataCollection()
data.objects.append(particles)

# 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 = particles.create_property('Position', data=pos)

# Create the particle type property and insert two atom types:
type_prop = 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 = 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]]
particles.bonds = Bonds()
particles.bonds.create_property('Topology', data=bond_topology)

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

## 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¶

This example presents a user-defined modifier function for calculating the mean square displacement (MSD) for a system of moving particles. OVITO provides the built-in Displacement Vectors modifier, which calculates the individual displacement of each particle. It stores its results in the `"Displacement Magnitude"` particle property. So all our user-defined modifier function needs to do is to sum up the squared displacement magnitudes and divide by the number of particles:

```def modify(frame, data):

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

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

# Output MSD value as a global attribute:
data.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 function from within a non-interactive batch script, which is run with the `ovitos` interpreter. Then we have to insert the `CalculateDisplacementsModifier` programmatically:

```from ovito.io import import_file, export_file
from ovito.modifiers import 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, data):

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

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

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

# Insert user-defined modifier function into the data pipeline.
pipeline.modifiers.append(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/attr",
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, 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)

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.

### 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.

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))
```

### 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, data):

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

# Create 'Selection' output particle property
selection = data.particles_.create_property('Selection')

# Prepare neighbor finder
finder = CutoffNeighborFinder(overlap_distance, data)

with selection:

# Iterate over all particles
for index in range(data.particles.count):

# Update progress display in the status bar.
yield (index / data.particles.count)

# 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.

```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¶

Ths user-defined viewport overlay function demonstrates how to use the Matplotlib Python module to render the radial distribution function, which is dynamically computed by a `CoordinationAnalysisModifier` in the data pipeline, on top the three-dimensional visualization.

```import matplotlib
matplotlib.use('Agg') # Activate 'agg' backend for off-screen plotting.
import matplotlib.pyplot as plt
import PyQt5.QtGui

def render(args):
# Request the output data collection from the current pipeline:
data = args.scene.selected_pipeline.compute()
# Look up the DataSeries object generated by the CoordinationAnalysisModifier:
if 'coordination-rdf' not in data.series:
raise RuntimeError('No RDF data found')
rdf_data = data.series['coordination-rdf'].as_table()

#  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 RDF histogram data
ax.bar(rdf_data[:,0], rdf_data[:,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¶

```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.positions
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 the particle.
if radius <= 0 and data.particles.particle_types is not None:
particle_type = data.particles.particle_types[pindex]

# Calculate screen-space size of the particle in pixels.

# Draw a dashed circle around the particle.
pen = QPen(Qt.DashLine)
pen.setWidth(3)
pen.setColor(QColor(0,0,255))
args.painter.setPen(pen)

# 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)