Forum Navigation
You need to log in to create posts and topics.

Frequent Voronoi polyhedra at each frame

I am trying to calculate the most frequent polyhedra in my system using the OVITO script:

https://www.ovito.org/docs/current/python/introduction/examples/batch_scripts/voronoi_indices.php

 

I can easily do it for a single frame, but when I try to calculate this at each frame, I get some error message.

Could you please help.

 

Thanks.

 

 

Hi,

can you show me your code and the error message?

-Constanze

I just got the code running that calculates Voronoi polyhedra and indices for all the 1000 frames. I am trying to write the most frequent polyhedra at each frame to a file with this code:

 

def render(args):
    # # This demo code prints the current animation frame into the upper left corner of the viewport.
    text1 = "Frame {}".format(args.frame)
    args.painter.drawText(10, 10 + args.painter.fontMetrics().ascent(), text1)

    # Also print the current number of particles into the lower left corner of the viewport.
    pipeline = args.scene.selected_pipeline
    if pipeline:
        data = pipeline.compute(args.frame)
        num_particles = data.particles.count if data.particles else 0
        text2 = "{} particles".format(num_particles)
        args.painter.drawText(10, args.painter.window().height() - 10, text2)

# 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("snap.1.lmp")

# Set atomic radii (required for polydisperse Voronoi tessellation).

atom_types = pipeline.source.data.particles['Particle Type'].types
atom_types[0].radius = 1.70 # Ta atomic radius (atom type 1 in input file)
atom_types[1].radius = 0.66 # O atomic radius (atom type 2 in input file)
atom_types[2].radius = 1.75 # Zr atomic radius (atom type 3 in input file)

# Set up the Voronoi analysis modifier.
voro = VoronoiAnalysisModifier(
    compute_indices = True,
    use_radii = True,
    edge_threshold = 1.0 )

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

out1 = open("a.dat",'a') #Output file for holding the 1000 sets of most frequent VP

# 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)))
    #Write the output also to a file
    print("%s\t%i\t(%.1f %%)" % (tuple(unique_indices[i]),counts[i],100.0*float(counts[i])/len(voro_indices)),file=out1)

out1.close()
    

 

This however only writes a single set of the most frequent polyhedra.

I don't understand why this issue arises.

 

Thanks,

 

Thanks for uploading the script, that makes it easier to help. Your script only computes the Voronoi statistics for the first frame of your trajectory.

Just to make sure, "snap.1.lmp" is a complete trajectory, right? A file sequence named "snap.1.lmp", "snap.2.lmp" etc., should be imported like this:

pipeline = import_file("snap.*.lmp")

Please note that when you call pipeline.compute() and don’t specify a frame number, the frame that will be evaluated will always be frame 0.
https://www.ovito.org/docs/current/python/modules/ovito_pipeline.php#ovito.pipeline.Pipeline.compute

If you would like to evaluate the complete trajectory, a for-loop of the following form is used to iterate over all frames:

for frame in range(pipeline.source.num_frames): 
    data = pipeline.compute(frame)

https://www.ovito.org/docs/current/python/introduction/pipelines.php

Could you explain why you are using lines 1 to 13? This is a python script overlay, that you would use to include additional information into your rendered images. Since your script does not render any images or animations, you can remove this part.

If you would like to write all histograms to a single file you could modify the last part of your script like this:

out1 = open("a.dat",'a') 

for frame in range(pipeline.source.num_frames):
    data = pipeline.compute(frame)
    voro_indices = data.particles["Voronoi Index"]
    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)))
        #Write the output also to a file
        print("%s\t%i\t(%.1f %%)" % (tuple(unique_indices[i]),counts[i],100.0*float(counts[i])/len(voro_indices)),file=out1)

out1.close()

 

-Constanze

One thing, I'd like to add here. If you prefer to work in the GUI, you can do that. First, set the atom radii and then add the Voronoi analysis modifier to calculate the Voronoi indices.

Then, in a next step, you could use a Python script modifier to calculate the Voronoi statistics and then make the histogram available in the Data Inspector as a "Data Table". Copy and paste the following script into the script editor:

from ovito.data import *
import numpy 

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


def modify(frame, data):
    # Compute frequency histogram.
    voro_indices = data.particles["Voronoi Index"]
    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)))    
    
    #Store information as Data Table (see Data inspector)
    table = DataTable(identifier='voronoi-statistics',title='Voronoi Statistics', plot_mode=DataTable.PlotMode.BarChart)
    table.x = table.create_property('Voronoi polyhedra', data=numpy.arange(10))
    for i in range(10):
        #For the plot I only used <n3 n4 n5 n6 n7>, adapt to your needs
        vp = "<{}>".format( str(unique_indices[i][2:7]).replace(']', '').replace('[', '') )
        table.x.types.append(ElementType(id=i, name=vp))
  
    table.y = table.create_property('Counts', data=counts[:10])
    data.objects.append(table)

That way, you can easily export the time evolution of that histogram directly through the Export function in the Data inspector (as marked in blue in the screenshot).

Uploaded files:
  • You need to login to have access to uploads.

And by the way, you can use the "Data Table" object you have generated to plot in the Viewport. Here is in example script of how to use the Python Script viewport layer. https://www.ovito.org/docs/current/viewport_layers.python_script.php

The script I have used is attached below.

Sorry for the flood of information, I just wanted to give you a quick overview of what OVITO Pro is capable of. 😉

 

Uploaded files:
  • You need to login to have access to uploads.

Thank you for your help. I will try all of these and get back to you if I get into any issue.