""" ovito_localization_factor.py OVITO Pro "Python script" modifier that computes the strain Localization Factor (LF) of Duan, Sansoz & Pan, Int. J. Mech. Sci. 317 (2026) 111534, Eq. (1): LF = (1/V) * sum_{ithr, eps_j>thr} [ (eps_i/eps_avg - 1) * (eps_j/eps_avg - 1) / r_ij ] Only particle pairs within a user-set cutoff distance are summed, using OVITO's CutoffNeighborFinder (a cell-list / spatial-binning search), giving ~O(N) scaling for large atomistic configurations, as recommended in Section 3.4 of the paper. HOW TO USE IN OVITO PRO ------------------------ 1. Load your simulation/trajectory and make sure a per-particle strain property exists (e.g. computed via a "Compute property" modifier, the "Atomic strain" modifier, or imported from file). By default this script looks for a property called "Strain"; change STRAIN_PROPERTY_NAME below (or set it as a global variable, see "Parameters" section) to match whatever your strain property is named (e.g. "Shear Strain", "von Mises Strain", a custom property, etc.). If the property has multiple components (e.g. a strain tensor), set STRAIN_COMPONENT accordingly, or precompute a scalar (e.g. von Mises) beforehand with a "Compute property" modifier. 2. In OVITO Pro: Add modification > Python script, then either: (a) paste the contents of this file into the modifier's code editor, or (b) point the modifier to this file on disk ("Load script file"). 3. The modifier adds: - A new particle property "LFCharge" (the normalized "point charge" q_i = eps_i/eps_avg - 1 for points above threshold, 0 otherwise) — useful for visualizing which particles are contributing. - Global attributes: "LocalizationFactor" (LF, Eq. 1) "ModifiedLocalizationFactor" (MLF = LF * average strain) "AverageStrain" (eps_avg used for normalization) "NumActivePoints" (number of points above threshold) These attributes appear in the OVITO data inspector ("Attributes" tab) and can be plotted vs. time/frame using "Add modification > Time series" in OVITO Pro, or exported as a table via File > Export File > "Table of values" / "Calculation results text file". PARAMETERS (edit these to suit your data) ------------------------------------------ """ import numpy as np from ovito.data import DataCollection, CutoffNeighborFinder from ovito.pipeline import ModifierInterface # ---------------------------------------------------------------------- # Parameters - edit these for your dataset # ---------------------------------------------------------------------- # Name of the per-particle strain property already present in the pipeline. STRAIN_PROPERTY_NAME = "Shear Strain" # If the strain property has multiple components (e.g. a tensor or vector), # set the component index to use (0-based). Leave as None for a scalar # property. STRAIN_COMPONENT = None # Cutoff distance for the neighbor search, in the same length units as the # particle coordinates (e.g. Angstrom). Per the paper (Sec. 3.4), choose a # cutoff large enough to span the separation between localized high-strain # regions in your configuration. CUTOFF_DISTANCE = 15.0 # Strain threshold above which points are considered "high strain" and # contribute to the LF sum. Set to None to use the average strain # (as done throughout the paper). STRAIN_THRESHOLD = None # Sample volume V used for normalization. Set to None to use the # simulation cell volume reported by OVITO (data.cell.volume). SAMPLE_VOLUME = None def _get_strain_array(data: DataCollection) -> np.ndarray: """Extract a 1D numpy array of per-particle strain values.""" if STRAIN_PROPERTY_NAME not in data.particles: raise KeyError( f"Particle property '{STRAIN_PROPERTY_NAME}' not found. " f"Available properties: {list(data.particles.keys())}. " "Set STRAIN_PROPERTY_NAME at the top of this script to match " "your data (e.g. add an 'Atomic strain' or 'Compute property' " "modifier upstream to create it)." ) prop = data.particles[STRAIN_PROPERTY_NAME] arr = np.asarray(prop) if arr.ndim == 1: return arr.astype(float) if STRAIN_COMPONENT is None: raise ValueError( f"Property '{STRAIN_PROPERTY_NAME}' has {arr.shape[1]} " "components; set STRAIN_COMPONENT to select which one to use, " "or precompute a scalar (e.g. von Mises strain) beforehand." ) return arr[:, STRAIN_COMPONENT].astype(float) def compute_localization_factor( data: DataCollection, epsilon: np.ndarray, cutoff: float, volume: float = None, threshold: float = None, ): """ Core LF calculation using OVITO's CutoffNeighborFinder for an efficient, cell-list-based neighbor search (periodic boundary conditions handled automatically by OVITO from the DataCollection). Parameters ---------- data : OVITO DataCollection (positions, cell, and PBC are read from here by CutoffNeighborFinder automatically) epsilon : (N,) ndarray of per-particle strain values cutoff : neighbor search cutoff distance (same units as coordinates) volume : sample volume for normalization; defaults to cell.volume threshold : strain threshold for "active" points; defaults to eps_avg Returns ------- lf : float mlf : float eps_avg : float n_active : int charge : (N,) ndarray per-particle "point charge" q_i (0 if inactive) """ n_points = len(epsilon) eps_avg = float(np.mean(epsilon)) if eps_avg == 0: raise ValueError("Average strain is zero; cannot normalize Eq. (1).") thr = eps_avg if threshold is None else threshold active_mask = epsilon > thr n_active = int(np.sum(active_mask)) charge = np.zeros(n_points, dtype=float) charge[active_mask] = epsilon[active_mask] / eps_avg - 1.0 if volume is None: volume = data.cell.volume if volume <= 0: raise ValueError( "Sample volume is zero or undefined; set SAMPLE_VOLUME " "explicitly or ensure the simulation cell is defined." ) lf_sum = 0.0 if n_active >= 2: # CutoffNeighborFinder(cutoff, data) — OVITO reads positions, cell, # and PBC directly from the DataCollection; no extra arguments. finder = CutoffNeighborFinder(cutoff, data) active_indices = np.nonzero(active_mask)[0] # To avoid double counting (i