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

Regarding computing CNA bond indices

I have loaded the dump file then applied the Create bonds modifier and then the common neighbor analysis modifier then i have applied the python script modifier and wrote the script as given in the manual as follows

from ovito.data import *

# 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("input/simulation.0.dump")

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

 

 

It shows the following error :

The Python script has exited with an error.

Traceback (most recent call last):

File "<string>", line 15, in <module>

File "C:\Program Files\Ovito Pro\plugins\python\ovito\io\import_file.py", line 120, in import_file

importer = FileImporter.autodetect_format(ovito.scene, first_location)

RuntimeError: File does not exist:

input/simulation.0.dump

 

 

What to do now how do i get a file with the CNA bond indices in it

 

Hi,

please note that there is a difference between stand-alone python scripts, that are to be executed from the terminal, and python script modifiers in the GUI. The first few pages in the scripting docs will give you an introduction into this topic: https://www.ovito.org/docs/current/python/introduction/running.php

The error message arises because you are trying to import a dump file called "simulation.0.dump" from within the python script modifier in the GUI, which is not possible. If you are working in the GUI, your pipeline is automatically built when you import your file. Simply remove that part from the script.

Since you said you already added the CreateBonds and Common Neighbor Analysis modifiers to your pipeline, it makes no sense to do that again in the Python script. You can shorten the script and copy and paste the following to the script editor into the GUI:

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)
    return (a[indices], counts)

def modify(frame, data):

    # Used below for enumerating the bonds of each particle:
    bond_enumerator = BondsEnumerator(data.particles.bonds)
    cna_indices = data.particles.bonds['CNA Indices']
    # 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")

-Constanze

 

 

 

Thanks a lot for the reply

but the code you suggested gives the following error

1 The Python script has exited with an error.

Traceback (most recent call last):

File "<string>", line 19, in modify

NameError: name 'cna_indices' is not defined

My mistake, I corrected that in line 14 in the script above.

Thankyou,I have generated the cna indices now i want to know the no .of pairs with the cna values 5 and 1 as their first and second indices and  the no .of pairs with the cna values 4 and 3 as their first and second indices in order to plot a graph

My notepad file has values shown like this

6 [4 4 4]:5 [4 5 5]:1 [6 6 6]:4 [6 7 7]:2 [7 9 9]:2
798 [4 4 4]:5 [4 5 5]:1 [6 6 6]:4 [6 7 7]:2 [7 9 9]:2
801 [4 4 4]:4 [4 5 5]:1 [5 6 6]:2 [6 6 6]:6 [7 9 9]:2
808 [4 4 4]:4 [4 5 5]:1 [5 6 6]:2 [6 6 6]:6 [7 9 9]:2
984 [4 4 4]:4 [5 6 6]:2 [6 6 6]:6 [6 7 7]:2
1939 [4 4 4]:3 [5 5 5]:6 [6 6 6]:7

how can i do this ???

Does Ovito has any option to do this ????

Please make sure you understand every single line of the script before you use it. In this example script, you loop over all particles and create a histogram of the CNA-indices of every particle and print that information along with the particle id. So the first line would be the triplet count of particle 6, the second one for particle 798 etc.

If you want the overall statistics of your system, you can modify your script such that you remove the part where you loop over all particles and directly create a histogram of the CNA indices.

Concerning your last question, yes this is doable with just a few lines of python code, the following code will only print (5,5,X) triplets and their frequencies.

# Print list of triplets with their respective counts.
for triplet, count in zip(unique_triplets, triplet_counts):
    if triplet[0] == 5 and triplet[1] == 5:
        sys.stdout.write("%s:%i \n" % (triplet, count))

So in summary your modifier script could look like this:

from ovito.data import *
import numpy
import sys

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)

def modify(frame, data):

    cna_indices = data.particles.bonds['CNA Indices']
    unique_triplets, triplet_counts = row_histogram(cna_indices)

    for triplet, count in zip(unique_triplets, triplet_counts):
        if triplet[0] == 5 and triplet[1] == 5:
            sys.stdout.write("%s:%i \n" % (triplet, count))
    
    for triplet, count in zip(unique_triplets, triplet_counts):
        if triplet[0] == 4 and triplet[1] == 3:
            sys.stdout.write("%s:%i \n" % (triplet, count))

-Constanze