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)
