#Import necessary OVITO libraries
from ovito.io import import_file
from ovito.modifiers import VoronoiAnalysisModifier, ExpressionSelectionModifier
import numpy as np

#This is the list of temperatures to be used for each frame.
temperatures = [
300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000,
1050, 1100, 1150, 1200, 1250, 1300, 1350, 1400, 1450, 1500, 1550, 1600, 1650,
1700, 1750, 1800, 1850, 1900, 1950, 2000, 2050, 2100, 2150, 2200, 2250, 2300,
2350, 2400, 2450, 2500
]

#Load your simulation trajectory file.
pipeline = import_file(“…/cooling.lammpstrj”)

#1. Perform Voronoi tessellation analysis.
pipeline.modifiers.append(VoronoiAnalysisModifier(compute_indices=True))

#2. Define the selection expression based on your Voronoi index criteria.
selection_expression = “(VoronoiIndex.3 == 0 && VoronoiIndex.4 == 4 && VoronoiIndex.5 == 4 && VoronoiIndex.6 == 7) ||
(VoronoiIndex.3 == 0 && VoronoiIndex.4 == 4 && VoronoiIndex.5 == 4 && VoronoiIndex.6 == 6) ||
(VoronoiIndex.3 == 0 && VoronoiIndex.4 == 4 && VoronoiIndex.5 == 4 && VoronoiIndex.6 == 5)”

#3. Apply the selection modifier to tag the matching atoms.
pipeline.modifiers.append(ExpressionSelectionModifier(expression=selection_expression))

#Print the header for the output data
print(“Temperature(K) Fraction_Voronoi(%)”)

#Loop through each frame (timestep) of the simulation
#The loop will stop if it runs out of frames or temperatures.
num_steps_to_process = min(pipeline.source.num_frames, len(temperatures))
for frame in range(num_steps_to_process):
# Process the data pipeline for the current frame
data = pipeline.compute(frame)

# Get the temperature from the predefined list using the frame index.
temperature = temperatures[frame]

# Get the total number of atoms
total_atoms = data.particles.count

if total_atoms > 0:
    # --- MODIFIED SECTION ---
    # Instead of reading a summary attribute, we will count the selected particles manually.
    # The 'ExpressionSelectionModifier' creates a 'Selection' property for each particle (1=selected, 0=not).
    
    # Get the 'Selection' property array from all particles.
    selection_property = data.particles['Selection']
    
    # Count how many particles are non-zero (i.e., selected).
    num_selected = np.count_nonzero(selection_property)
    # --- END OF MODIFIED SECTION ---

    # Calculate the fraction as a percentage
    fraction_percentage = (num_selected / total_atoms) * 100