Author Topic: How to find out the distribution of Voronoi indices? I got some error in script  (Read 251 times)

RMCtestFYP

  • Newbie
  • *
  • Posts: 18
Hello, I am now trying to find the 10 most frequency Voronoi indices.
(The sample .cfg file have been uploaded for reference.)

I encounter some problems.
I copy the script from this example. https://ovito.org/manual/python/introduction/examples.html
And make some change, I highlight them in red colour below.
But some error occur.

Here is the possible problems I guess.
1. The .cfg file have some error. The cell size is wrong, or the coordinate of particles is not in the range [0 1].

Sorry for my silly question, since I have few experience on coding.
I search some key word in google, but I still don't know the reason.
I try the script several time and encounter several error too.
I don't know how to fix the problems :'( Are there any problems in the 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("C:/Users/wai/Desktop/example.cfg.cfg")

# Set atomic radii (required for polydisperse Voronoi tessellation).
atypes = node.source.particle_properties.particle_type.type_list
atypes[1].radius = 1.35        # Pd
atypes[2].radius = 1.55        # Cu
atypes[3].radius = 1.55        # P


# 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),
                                 counts,
                                 100.0*float(counts)/len(voro_indices)))

« Last Edit: July 31, 2018, 05:41:30 PM by RMCtestFYP »

Constanze Kalcher

  • Administrator
  • Jr. Member
  • *****
  • Posts: 81
Hi,

there are some changes in the python programming interface of OVITO 3.0 as compared to OVITO 2.X and you're using the example script from the OVITO 2.X manual.
Let me refer you to the updated script in the OVITO 3.0 documentation.

However, from your screen-shot I can see that you're using the graphical user interface, where you can make your life a lot easier by calculating the Voronoi indices by inserting
the Voronoi analysis modifier in the modification pipeline (see attached screen-shot). Afterwards, you can still use a python script modifier function to print the 10 most frequent Voronoi motifs, e.g. like this:

Code: [Select]
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, input, output):
   
    voro_indices = input.particles['Voronoi Index']
    # 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)))
   

-Constanze


RMCtestFYP

  • Newbie
  • *
  • Posts: 18
Thank you very much again! ;)
This method is much convenient, which can give the result very fast.
For this function, I still have some questions.

(1) These two method gives the same result?

(2) As picture 'Result.jpg', is this result reasonable for metallic glass? The value is very low.
Am I do anything wrong? :'(

(2) I have tried the new code from OVITO 3.0.0-dev211 documentation.
But I still got error when using the new code.
I type wrong code on atom type? ???

The thing I have changed is:
pipeline = import_file("C:/Users/wai/Desktop/example.cfg.cfg")

atom_types[0].radius = 1.35   # Pd atomic radius (atom type 1 in input file)
atom_types[1].radius = 1.55   # Cu atomic radius (atom type 2 in input file)
atom_types[2].radius = 1.55   # P atomic radius (atom type 3 in input file)



Code:
Code: [Select]
# 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("C:/Users/wai/Desktop/example.cfg.cfg")

# Set atomic radii (required for polydisperse Voronoi tessellation).
atom_types = pipeline.source.compute().particles['Particle Type'].types
atom_types[0].radius = 1.35   # Pd atomic radius (atom type 1 in input file)
atom_types[1].radius = 1.55   # Cu atomic radius (atom type 2 in input file)
atom_types[2].radius = 1.55   # P atomic radius (atom type 3 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)))

Constanze Kalcher

  • Administrator
  • Jr. Member
  • *****
  • Posts: 81
Dear RMCtestFYP,

you can't run whole batch scripts from within a python script modifier in the graphical user interface (GUI) of OVITO.
Have a look at the documentation part about running scripts.

So either you execute the voronoi example script from within the terminal, e.g.
Code: [Select]
ovitos voronoi.py

or you follow the GUI procedure which I explained in my previous post. Note that if you want to reproduce what's happening in the example batch script, you need to activate
Use particle radii in the Voronoi analysis modifier.

The particle radii can be edited when you click on -> Particle Types in the Data source panel (see screenshot).

Concerning your results I can only suggest you to consult the literature. My expertise is mostly on CuZr metallic glasses.

-Constanze

RMCtestFYP

  • Newbie
  • *
  • Posts: 18
OH... it seems that it is too complicated for me... I don't understand :'(
But it' okay, I still can use the code you provide, it works, thank you! :)

1. For the script you provide, it is to find the distribution of Voronoi indices for all types of atoms?

For example, my input data contains 3 types of atoms
if I want to find the distribution of Voronoi indices of type-1 atom only (e.g. The Voronoi polyhedron around type-atom),
Can I find it by OVITO?


2 To find the Voronoi indice, do I need to click 'Use particle radii' in 'Voronoi analysis'?

Code: [Select]
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, input, output):
   
    voro_indices = input.particles['Voronoi Index']
    # 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)))



3. how to know that whether the script is batch scripts?

4. DO OVITO able to find the distribution of polyhedrons?

Thank you for you reply! Sorry if I ask some silly questions... since I am a super rookie in coding :'(
« Last Edit: August 02, 2018, 03:39:32 AM by RMCtestFYP »

Constanze Kalcher

  • Administrator
  • Jr. Member
  • *****
  • Posts: 81
Hi,
let me try to answer your questions in the order you posted them.  :)

1. Yes, but it's not the script that controls that, it's the Voronoi analysis modifier. You can simply use the Select type modifier in conjunction with the option "Use only selected particles" in the Voronoi analysis modifier.

2. No you don't need to. It's up to you. If you don't use that option then all particles will be treated as if they have the same size during the Voronoi tesselation. Otherwise, bigger atoms will take more space so to say.

3. In a non-interactive batch script (=outside the graphical user interface) you import your data from an external file (as shown below) to start the pipeline and then you add all your modifiers similarly to what you would do in the graphical user interface.

The python script modifier in the graphical user interface is only meant to define modifier functions
Code: [Select]
def modify(frame, input, output):
    ....
and is part of the pipeline. It gets the your imported data as input. That's why it failed when you tried to import your data (twice).
Code: [Select]
node = import_file("C:/Users/wai/Desktop/example.cfg.cfg")

4. Well, that's what the function called "modify" your using here does, right?
Code: [Select]
unique_indices, counts = row_histogram(voro_indices)

-Constanze

RMCtestFYP

  • Newbie
  • *
  • Posts: 18
Thank you for you patient reply! :)
It seems that learning more concept of coding is a must.
Be frankly, I still don't know something you said. It looks a bit complicated, haha :)