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[0].radius = 1.35        # Cu atomic radius (atom type 1 in input file)
atypes[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
)
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[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 %)

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[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 = 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")

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.

from ovito import *
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[0] = (1.0, 1.5, 0.3)
pos_prop.marray[1] = (7.0, 4.2, 6.0)
pos_prop.marray[2] = (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[0] = 1  # First atom is Cu
type_prop.marray[1] = 2  # Second atom is Ni
type_prop.marray[2] = 2  # Third atom is Ni

# Create a user-defined particle property.
my_prop = ParticleProperty.create_user('My property', 'float', num_particles)
my_prop.marray[0] = 3.141
my_prop.marray[1] = 0.0
my_prop.marray[2] = 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()
data.add(pos_prop)
data.add(type_prop)
data.add(my_prop)
data.add(cell)
data.add(bonds)

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

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