Dear kalcher,

Thank you for your detailed reply. Regarding the calculation of local entropy, I would like to help me some questions. I am trying to use the script below to perform a multi-dimensional (1D, 2D) calculation of the local entropy of a certain atom type in a certain area. In the code below, I don’t know how to create the property ‘Entropy’ on the SpatialBinningModifier to achieve this goal. Any help is greatly appreciated. Also I know I’m calling ‘Entropy’ without defining it in my script, which is an obvious mistake. Do you have any suggestions for this?

# -*- coding: utf-8 -*-

“”"

Created on Wed Nov 16 10:46:40 2022

@author: ms

“”"

import numpy as np

from ovito.data import CutoffNeighborFinder, DataCollection

from ovito.io import import_file

from ovito.modifiers import SpatialBinningModifier,ExpressionSelectionModifier,ComputePropertyModifier

import pandas as pd

import time

def compute_local_entropy(frame: int, data: DataCollection, cutoff=5.0, sigma=0.2, use_local_density=True, compute_average=True, average_cutoff=5.0):

assert(cutoff > 0.0)

assert(sigma > 0.0 and sigma < cutoff)

assert(average_cutoff > 0)

```
global_rho = data.particles.count / data.cell.volume
finder = CutoffNeighborFinder(cutoff, data)
local_entropy = np.empty(data.particles.count)
nbins = int(cutoff / sigma) + 1
r = np.linspace(0.0, cutoff, num=nbins)
rsq = r**2
prefactor = rsq * (4 * np.pi * global_rho * np.sqrt(2 * np.pi * sigma**2))
prefactor[0] = prefactor[1]
for particle_index in range(data.particles.count):
r_ij = finder.neighbor_distances(particle_index)
r_diff = np.expand_dims(r, 0) - np.expand_dims(r_ij, 1)
g_m = np.sum(np.exp(-r_diff**2 / (2.0*sigma**2)), axis=0) / prefactor
if use_local_density:
local_volume = 4/3 * np.pi * cutoff**3
rho = len(r_ij) / local_volume
g_m *= global_rho / rho
else:
rho = global_rho
integrand = np.where(g_m >= 1e-10, (g_m * np.log(g_m) - g_m + 1.0) * rsq, rsq)
local_entropy[particle_index] = -2.0 * np.pi * rho * np.trapz(integrand, r)
data.particles_.create_property('Entropy', data=local_entropy)
if compute_average:
data.apply(ComputePropertyModifier(
output_property = 'Entropy',
operate_on = 'particles',
cutoff_radius = average_cutoff,
expressions = ['Entropy / (NumNeighbors + 1)'],
neighbor_expressions = ['Entropy / (NumNeighbors + 1)']))
```

start = time.perf_counter()

binx = 81

binz = 1

pipeline = import_file(“./mydump-nve.lammpstrj”, multiple_frames=True)

print(“Number of MD frames:”, pipeline.source.num_frames)

pipeline.modifiers.append(ExpressionSelectionModifier(expression=‘Position.Y >0 && Position.Y <=3 && ParticleType !=1 && ParticleType !=2 && ParticleType !=3’))

# Prepare a spatial binning modifier but do not append it to the pipeline yet

binning_modifier = SpatialBinningModifier(

property=‘Entropy’,

direction=SpatialBinningModifier.Direction.XZ,

bin_count=(binx, binz),

reduction_operation=SpatialBinningModifier.Operation.Mean,

only_selected=True

)

mean = np.zeros((binx, binz))

num_frames = pipeline.source.num_frames

half_frames = int(num_frames * 0.5)

for frame in range(half_frames, num_frames):

data = pipeline.compute(frame)

compute_local_entropy(frame, data)

```
# Apply the binning modifier locally
binning_modifier.modify(data)
grid = np.array(data.grids["binning"]["Entropy"]).reshape((binx, binz))
mean += grid
```

mean /= (num_frames - half_frames)

end = time.perf_counter()

print('Running time of CPU: s Seconds' (end - start))

mydump-nve.lammpstrj (888.9 KB)