# Example B1: Computing Voronoi indices

This script demonstrates the use of the `VoronoiAnalysisModifier`. 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("input/simulation.0.dump")

    # The LAMMPS dump file imported above contains only numeric atom type IDs but
    # no chemical element names or atom radius information. That's why we explicitly set the
    # atomic radii of Cu & Zr atoms now (required for polydisperse Voronoi tessellation).
    def assign_particle_radii(frame, data):
        atom_types = data.particles_.particle_types_
        atom_types.type_by_id_(1).radius = 1.35   # Cu atomic radius assigned to atom type 1
        atom_types.type_by_id_(2).radius = 1.55   # Zr atomic radius assigned to atom type 2
    pipeline.modifiers.append(assign_particle_radii)

    # Set up the Voronoi analysis modifier.
    voro = VoronoiAnalysisModifier(
        compute_indices = True,
        use_radii = 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 %)