Savi-Bhransha: Graph-Theoretic Dislocation-Loop Characterization in Crystals

Given the atomic coordinates of a simulation snapshot, the first step is to identify the defect content and assign each atom to its nearest ideal lattice site. Each atom is mapped to its nearest ideal lattice site via modular arithmetic: given axis periods $\mathbf{p}$ (lattice parameters) and fractional origin $\mathbf{o}$ which is similar to the offset or shift used in building the simulation box,

where $\{\mathbf{s}_{k}\}$ are the $M$ sub-lattice offsets within one unit cell ($M{=}2$ BCC, $M{=}4$ FCC/HCP) and $d_{\mathrm{PBC}}$ is the minimum-image distance. The fractional anchor is $a_{i}=\lfloor(r_{i}-o_{i}p_{i})/p_{i}\rfloor+s_{k^{*},i}/p_{i}+o_{i}$. The resulting fractional anchors are PBC-folded into $[\mathbf{o},\,\mathbf{o}+N_{\mathrm{cell}})$ and sorted lexicographically.

Lattice-site occupancy is then recovered by co-traversing the sorted anchor list against an on-the-fly enumeration of the ideal sites in the same order, so the full reference lattice is never held in memory. A site absent from the atom list is a vacancy; a site claimed by more than one atom is an interstitial. This recovers the same lattice-site occupancy as the Wigner-Seitz construction [29] without an explicit reference frame or spatial index; total cost is $O(N\log N)$ [4]. For each over-occupied site we identify the triad: the vacancy at the lattice position and two displaced atoms forming the dumbbell. The axis $\mathbf{d}$ connecting the two interstitial atoms defines a line primitive $\mathcal{L}(\mathbf{d})=\{\mathbf{r}_{\text{vac}}+\lambda\mathbf{d}\mid\lambda\in\mathbb{R}\}$. In addition, atoms whose displacement from their nearest lattice site exceeds a threshold (typically $0.3$–$0.4\times$ the nearest-neighbor spacing [9]) are paired with that lattice site and treated as displaced-atom–vacancy lines. These threshold-based pairs carry zero net defect content, but they capture the diffuse displacement field around extended defects and help in identifying morphological details in later steps. They typically appear in one of three configurations: coincident with a triad-bearing dumbbell line forming a crowdion, coincident with a true vacancy line in vacancy loops, or as non-aligned interface disturbances surrounding a well-defined dislocation or other defect structure.

Each dumbbell axis $\mathbf{d}$ is assigned to a low-index crystallographic direction family. This orientation governs both graph construction (Section 2.3) and, once loops are identified, determines the Burgers vector. The identification rests on the Mura–Willis representation of a closed loop of Burgers vector $\mathbf{b}$ bounding a habit surface $S$,

with $\Omega(\mathbf{r})$ the solid angle subtended by $S$ at $\mathbf{r}$ and $\mathbf{u}_{\mathrm{LI}}(\mathbf{r})$ the smooth line-integral remainder governed by the elastic Green’s tensor [28, 2]. The first term is everywhere parallel to $\mathbf{b}$ and carries the entire $\mathbf{b}$-jump across $S$; $\mathbf{u}_{\mathrm{LI}}$ contributes the smaller Poisson-type bulging. The inserted material at the core therefore aligns with $\mathbf{b}$ and is observed atomistically as a dumbbell or crowdion triad parallel to $\mathbf{b}$, so the per-primitive family assignment reads the local Burgers direction directly.

The relevant crystallographic direction families depend on crystal structure:

BCC: $\langle 111\rangle$ (close-packed, produces highly mobile loops when clustered), $\langle 100\rangle$ (non-close-packed, sessile loops), and $\langle 110\rangle$ (appears at interfaces between domains or in C15 ring structures); the higher-index $\langle 112\rangle$ and $\langle 221\rangle$ families are also optionally retained in the candidate set so that the less common dumbbell orientations encountered in dense cascade debris can still be assigned rather than being forced into a nearby low-index family.

FCC: $\langle 110\rangle$, $\langle 112\rangle$, $\langle 221\rangle$, $\langle 100\rangle$, and $\langle 111\rangle$ families, enabling separation of perfect and partial/faulted loop populations within the same crystallographic framework.

HCP: Three families expressed in Miller-Bravais notation $[UVTW]$:

• •

  $\langle a\rangle$-type: $\frac{a}{3}[11\bar{2}0]$ and symmetry-equivalent directions;

• •

  $\langle c\rangle$-type: $[0001]$ along the $c$-axis;

• •

  $\langle c+a\rangle$-type: $\frac{1}{3}[11\bar{2}3]$ (Type II) and $\frac{1}{3}[10\bar{1}1]$ (Type I).

Each dumbbell is assigned to a family using angular matching against pre-computed standard directions, with an integer family code enabling efficient computation in vectorized operations. The families can be added to or removed from the candidate set as needed.

With orientation families assigned, loop identification proceeds in two stages: collinear defect lines are first merged into extended line primitives, and a graph over these primitives then resolves loops and other morphological units. Both stages operate within a single defect cluster at a time; clusters are obtained beforehand by grouping point defects whose pairwise separation lies within the second-nearest-neighbor (2NN) distance, which keeps every subsequent pairwise computation local.

The two stages share a common pair geometry, parameterized by an inter-line distance and an inter-axis angle. The natural distance measure depends on the relative orientation of the two lines. When axes $\mathbf{d}_{i},\mathbf{d}_{j}$ are nominally parallel (the merging step and the parallel-edge criterion below), we use the perpendicular distance between the axes treated as infinite lines,

which reduces to the standard point-to-line distance in the strictly parallel limit. To suppress thermal scatter and local deviations from the ideal orientation, both axes are first projected onto the assigned family direction $\hat{\mathbf{d}}$ before $d_{\perp,ij}$ is evaluated. For non-parallel pairs (the C15 ring criterion), the perpendicular distance between infinite axes is poorly conditioned, and we instead use the underlying point distance between the two primitives—lattice-site centers for SIA-bearing dumbbell lines, vacancy positions for vacancy-centered lines—denoted $\|\mathbf{r}_{i}-\mathbf{r}_{j}\|$. The relative orientation in both cases is the unsigned angular difference $\theta_{ij}=\arccos(|\hat{\mathbf{d}}_{i}\cdot\hat{\mathbf{d}}_{j}|)$.

The habit plane of a dislocation loop is obtained by a weighted principal-component fit to its constituent line centroids $\{\mathbf{r}_{i}\}$, with weights $w_{i}$ proportional to the constituent-defect count of each line. The weighted covariance

is diagonalized to give eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}$ and corresponding eigenvectors $\mathbf{e}_{\alpha}$; the habit normal $\mathbf{n}=\mathbf{e}_{3}$ minimizes the weighted sum of squared point-to-plane distances. Two scalar diagnostics accompany the fit: the in-plane scatter

and the planarity ratio $R_{\text{planar}}=1-\lambda_{3}/\lambda_{1}$, which approaches unity for a well-defined planar loop and falls toward zero for three-dimensional or disordered components.

The fitted normal is then matched to a low-index plane family. For each candidate $(hkl)$ the Cartesian normal $\mathbf{n}_{hkl}=B(h,k,l)^{\top}$ is computed from the reciprocal basis and the deviation $\theta_{\text{dev}}=\arccos(|\mathbf{n}\cdot\mathbf{n}_{hkl}|)$ is taken as the match metric. Candidate families are $\{110\},\{111\},\{100\}$ for BCC and FCC, and $(0001)$ basal, $\{10\bar{1}0\}$ prismatic, $\{10\bar{1}1\}$ pyramidal-I, and $\{11\bar{2}2\}$ pyramidal-II for HCP. Each assignment carries a confidence label derived from the deviation: high ($<5^{\circ}$), medium ($5^{\circ}$–$15^{\circ}$), low ($15^{\circ}$–$25^{\circ}$), very low ($>25^{\circ}$).

With the Burgers family of each loop component already fixed by the connected-component analysis (Section 2.3.2), we now refine to a specific Burgers state from the admissible crystallographic set of the parent lattice and report the corresponding magnitude $|\mathbf{b}|$. Material dependence enters only through lattice geometry and the allowed Burgers set: for BCC and HCP we use the conventional canonical set for the structure, while for FCC the admissible set additionally includes fault-related partial dislocation loops and stacking faults [21, 24].

Within a single direction family the specific Burgers state is disambiguated using the resolved habit plane, the discrete defect content carried by the merged line primitives, and, where the local signal is sufficiently clear, the displaced-atom geometry associated with the component. In FCC, for example, both $\frac{1}{2}\langle 110\rangle$ and $\frac{1}{6}\langle 110\rangle$ states belong to the $\langle 110\rangle$ family; their crystallographic distinction is reflected in the resolved loop habit, with perfect loops associated with glide on $\{111\}$ and stair-rod character associated with $\{100\}$ junction geometry.

When the extra atom count per constituent line exceeds unity, it signifies a translation step of more than shortest lattice translation in that direction. The Burgers magnitude follows from the assigned crystallographic type as

where $f$ is the fractional coefficient and $\mathbf{d}$ is the primitive direction vector of the corresponding family.

Separating defects at the loop perimeter from those in the interior is essential for comparison with thermal-transport measurements, where phonon scattering depends predominantly on the boundary population [34, 7]. For each line $i$ in a parallel component we count its same-family neighbors within the graph-edge distance,

where $\mathbb{1}[\cdot]$ is the indicator function (unity when its condition holds, zero otherwise), so $n_{i}$ counts the line primitives that are both within the graph-edge distance $\delta_{\perp}$ of line $i$ and share its direction family. We then classify $i$ as boundary when $n_{i}<n_{\text{surf}}$. The threshold reflects the in-plane coordination of a fully embedded line on the relevant habit:

so that close-packed and basal-$\langle a\rangle$ loops carry their respective characteristic interior coordinations. The resulting boundary count $N_{\text{boundary}}$ and bulk count $N_{\text{bulk}}=N_{\text{total}}-N_{\text{boundary}}$ are object-level descriptors that can be compared directly with the boundary-sensitive observables of TGS and related experiments.

With the Burgers vector $\mathbf{b}$ and habit normal $\mathbf{n}$ in hand, the edge/screw character is resolved along the loop perimeter. The boundary defect coordinates are first projected onto the habit plane,

and expressed in an in-plane orthonormal basis $(\mathbf{u}_{1},\mathbf{u}_{2})$ as $(x_{i},y_{i})=(\mathbf{r}^{\prime}_{i}\cdot\mathbf{u}_{1},\mathbf{r}^{\prime}_{i}\cdot\mathbf{u}_{2})$. An $\alpha$-shape (concave hull) is then fitted to these 2D points, with $\alpha$ chosen tight enough to track the interatomic spacing of cascade defects, yielding an ordered perimeter sequence $\{\mathbf{r}^{\prime}_{k}\}$. The perimeter and enclosed area follow as

and the equivalent circular diameter $d_{\text{circ}}=2\sqrt{A/\pi}$ summarizes the loop size.

For each perimeter segment, the local tangent $\bm{\xi}_{k}=(\mathbf{r}^{\prime}_{k+1}-\mathbf{r}^{\prime}_{k})/|\mathbf{r}^{\prime}_{k+1}-\mathbf{r}^{\prime}_{k}|$ defines the unsigned angle $\phi_{k}\in[0,90^{\circ}]$ with the Burgers vector via $\cos\phi_{k}=|\bm{\xi}_{k}\cdot\hat{\mathbf{b}}|$ (the absolute value renders the decomposition invariant to traversal direction). The classical dislocation character decomposition then gives

with $f_{\text{screw}}+f_{\text{edge}}=1$. Segments are labeled screw ($f_{\text{screw}}>0.8$), edge ($f_{\text{edge}}>0.8$), or mixed, and a length-weighted mean

characterizes the overall loop.

In FCC metals with low stacking-fault energy, a perfect $\frac{1}{2}\langle 110\rangle\{111\}$ dislocation spontaneously dissociates into two Shockley partials of Burgers vector $\frac{1}{6}\langle 112\rangle$ bounded by a ribbon of intrinsic stacking fault [19]. The object-level loop descriptors do not by themselves separate the dissociated state from the undissociated one: both carry the same net $\frac{1}{2}\langle 110\rangle$ Burgers sum, and the dominant line-direction family inferred from the triad primitives can remain $\langle 110\rangle$ when the two partials sit close together. An optional analysis stage resolves the partial structure where it is present.

The two partial-dislocation cores that bound the stacking fault sit on two distinct $\{111\}$ layers, populated by dumbbell or crowdion triads with surviving displaced atoms; the stacking-fault ribbon between them carries no inserted material and is signalled instead by non-surviving displaced-atom and vacancy pairs. Operationally, the dumbbells of each $\langle 110\rangle$-family SIA cluster are partitioned into bands across the candidate $\{111\}$ habit normal, and the partial structure is read from the populations of those bands: a single populated band keeps the perfect $\frac{1}{2}\langle 110\rangle$ label; two well-populated outer bands (with an optional sparse middle band from transient SF-region atoms) are reported as a Shockley pair emitting two $\frac{1}{6}\langle 112\rangle$ partials on $\{111\}$; and three distinct $\{111\}$ planes with largely disjoint band populations identify a stair-rod junction emitting two Shockley pairs together with a $\frac{1}{6}\langle 110\rangle$ stair-rod along the intersection line.

A residual subset of clusters has all vacancy lattice sites coplanar on the primary $\{111\}$ but is nonetheless dissociated: the two SIA atoms of each dumbbell straddle their vacancy by $\pm\tfrac{1}{2}L\,\hat{b}$ along the modal $\langle 110\rangle$, so the parallel half-vectors (vacancy $\to$ SIA) sit on one $\{111\}$ layer and the anti-parallel half-vectors on an adjacent layer even when the vacancies themselves are coplanar. A secondary partition that projects the SIA atom positions onto the primary $\{111\}$ normal and bins them into adjacent integer layer indices recovers these Shockley dissociations whose two partial cores have glided onto adjacent $\{111\}$ layers while the dumbbell–vacancy line remains in the original plane.

The cost of the Savi-Bhransha pipeline decouples from the simulation-cell size after the initial defect-identification step. Defect identification itself is the only stage that touches every atom; line merging, graph construction, component classification, and morphological characterisation all operate solely on the $N_{\text{defect}}$ line primitives within an individual cluster. Because cascade simulation boxes are typically chosen large enough to suppress finite-size artefacts, the clusters they contain remain modest: a $100{\times}100{\times}100$ unit-cell box with $\sim$2 M atoms produces clusters of at most $10$–$10^{2}$ line primitives even at high PKA energies, so per-cluster analysis is essentially free relative to the defect-identification pass. Material dependence enters only through lattice geometry—nearest-neighbour distances and symmetry families—and every threshold in the pipeline is scaled from these lattice quantities rather than chosen as a bare cut-off.

The dominant cost within each cluster is pairwise line-distance evaluation. We expose a single computational kernel with three back-ends that share the same vectorized formulation: a CPU back-end using array broadcasting for the pairwise computation, a parallel-CPU back-end that distributes the per-line neighbor search across cores via just-in-time compilation, and a CUDA back-end that runs the same kernel on the GPU for the largest clusters. The resulting adjacency graph is sparse—each line has at most four to eight parallel neighbors set by the lattice—so it contains $O(N_{\text{defect}})$ edges rather than $O(N_{\text{defect}}^{2})$. Orientation families are carried as integer codes throughout, which keeps every inner loop free of string comparisons and amenable to vectorization.

Memory follows the same decoupling: only the per-cluster $O(N_{\text{defect}})$ line primitives and their sparse adjacency are held, in contrast to the $O(N_{\text{atoms}})$ connectivity data that mesh-based extraction must maintain over the entire cell. At the largest HCP system tested (172,M atoms), Savi-Bhransha completes the analysis in approximately 5,min using 24,GB, compared with 185,GB for DXA; the full benchmark is reported in Section 3.

The cascade trajectories analysed in Section 3 were generated with LAMMPS [33, 43] for two distinct loading regimes: single-PKA cascades for surveying primary damage as a function of energy and crystal structure, and successive collision cascades (SCC) for tracking microstructural evolution under accumulated dose. Both regimes share the same per-cascade procedure—NPT equilibration at 300 K and zero pressure, an NVE cascade phase with an adaptive timestep that keeps the fastest atom’s displacement below $0.1$ Å per step, and Lindhard–Scharff electronic stopping treated as a frictional drag [23]—and use ZBL-stiffened pair interactions at short range [55]. Material dependence enters through the choice of interatomic potential and threshold displacement energy $E_{d}$.

For BCC W we use two complementary cascade datasets: (i) a 100-cascade single-PKA energy sweep (10–150 keV, two potentials, 10 replicas per energy–potential combination) that fixes the morphology of primary damage as a function of PKA energy, and (ii) a successive-cascade trajectory at 50 keV up to 0.2 dpa that follows the same morphology under accumulated dose. We discuss the two in turn.

We apply Savi-Bhransha to 25 single-cascade simulations in hexagonal close-packed (HCP) zirconium spanning 10–125 keV using the Starikov–Smirnova ADP interatomic potential [38], with five replicas per energy. Figure 4 summarizes the extracted loop morphology.

The FCC FeNiCr alloy has a low stacking-fault energy, favouring Shockley-partial pairs that bound an intrinsic stacking fault [19]. Using Savi-Bhransha we report each surviving SIA cluster in two complementary views: a perfect-Burgers view and a partial-resolved view that decomposes $\langle 110\rangle$-family clusters into Heidenreich–Shockley partials wherever the local dumbbell geometry supports it (Section 2.8). The dataset is 250 successive 20 keV cascades in a $120^{3}$-unit-cell supercell ($E_{d}=40$ eV), accumulated to $7.205\times 10^{-3}$ dpa under the standard NRT expression $\nu_{\mathrm{NRT}}=0.8E/(2E_{d})$; vacancy loops and non-$\langle 110\rangle$ clusters pass through both views unchanged. Figure 5 summarises the loop morphology.

In this section we present qualitative and quantitative comparisons between Savi-Bhransha and the widely used Dislocation Extraction Algorithm (DXA) [40] as implemented in OVITO [42], drawing on the single-cascade BCC W and HCP Zr datasets and on a representative full-cell snapshot from the FCC FeNiCr trajectory. DXA, which has become the de-facto standard for recovering dislocation networks from atomistic configurations, is taken as the reference; the discussion focuses on where the two views agree, where they differ, and what the differences imply for the choice of observables in cascade-microstructure analysis. Both algorithms operate directly on molecular-dynamics snapshots in which the atoms vibrate continuously about their nominal positions, so the comparison is inherently sensitive to snapshot-to-snapshot fluctuations and to the parameter choices of either method—a point we return to repeatedly below.