jaxdem.utils.dynamicalMatrix#

Functions for calculating the dynamical matrix (hessian of the potential energy w.r.t. the system’s generalized coordinates).

  • non_bonded_hessian() Works for spheres and for deformable-particle (DP) particles with using the positions of the spheres in both cases as the coordinates.

  • bonded_hessian() Calculates the contribution from system.bonded_force_model (DP internal elastic / plastic energies: em, ec, eb, el, gamma). Adds to the non-bonded hessian by linearity: H_total = H_non_bonded + H_bonded.

  • clump_non_bonded_hessian() Works for rigid clumps using the center of mass coordinates and infinitesimal rotations \((\delta r, \omega)\). Can also define it in terms of a scaling of the rotations \((\delta r, R \omega)\).

  • Works unchanged for any existing (or future) force model that correctly implements energy.

  • The bonded and non-bonded contributions stay algorithmically independent; the total dynamical matrix is the simple sum.

Functions

bonded_hessian(state, system)

Dense (N*dim, N*dim) hessian of the bonded potential energy \(\partial^2 U_{bonded} / \partial r \, \partial r\).

clump_non_bonded_hessian(state, system[, ...])

Dense (n_clumps*group_dim, n_clumps*group_dim) clump hessian.

non_bonded_hessian(state, system[, cutoff, ...])

Assemble the dense (N*dim, N*dim) non-bonded hessian of \(\partial^2 U / \partial r \, \partial r\) where r concatenates all sphere positions.

pair_non_bonded_hessian(state, system[, ...])

Per-pair non-bonded hessian blocks for every neighbor-list pair.

zero_mode_mask(eigenvalues[, rel_gap])

Boolean mask identifying numerically-zero eigenvalues via gap detection.
jaxdem.utils.dynamicalMatrix.pair_non_bonded_hessian(state: State, system: System, cutoff: float | None = None, max_neighbors: int | None = None) → tuple['State', 'System', jax.Array, jax.Array][source]#

Per-pair non-bonded hessian blocks for every neighbor-list pair.

Returns:

  • state (State) – Potentially updated state (after neighbor-list rebuild).

  • system (System) – Potentially updated system.

  • pair_ids (jax.Array) – (M, 2) int array of (i, j) sphere index pairs, where M = state.N * max_neighbors. Padding pairs have j == -1 and corresponding blocks entries are zero.

  • blocks (jax.Array) – (M, 2*dim, 2*dim) hessian blocks. blocks[k, :dim, :dim] is \(\partial^2 \phi / \partial r_i^2\), the upper-right sub-block is \(\partial^2 \phi / (\partial r_i \partial r_j)\), and so on.

jaxdem.utils.dynamicalMatrix.non_bonded_hessian(state: State, system: System, cutoff: float | None = None, max_neighbors: int | None = None) → tuple['State', 'System', jax.Array][source]#

Assemble the dense (N*dim, N*dim) non-bonded hessian of \(\partial^2 U / \partial r \, \partial r\) where r concatenates all sphere positions.

Scatters the per-pair blocks from pair_non_bonded_hessian() into a dense matrix. Symmetric by construction; translational null modes (sum of rows in each block-row == 0) follow from each pair’s 4-block structure.

Each pair appears twice in the neighbor list ((i, j) and (j, i)) and each lane scatters the full 4-quadrant hessian into (i, i), (i, j), (j, i), (j, j), so every slot accumulates two identical contributions. We divide by 2 at the end to undo this, matching the 0.5 * sum convention used in NeighborList.compute_potential_energy().

See pair_non_bonded_hessian() for how padding entries are handled — padding blocks are zero and contribute nothing.

jaxdem.utils.dynamicalMatrix.bonded_hessian(state: State, system: System) → tuple['State', 'System', jax.Array][source]#

Dense (N*dim, N*dim) hessian of the bonded potential energy \(\partial^2 U_{bonded} / \partial r \, \partial r\).

Returns zeros when system.bonded_force_model is None. The total dynamical matrix is the sum bonded_hessian + non_bonded_hessian (by linearity).

Note: this function does not apply a 0.5 double-count correction. The bonded potential energy is a plain sum(over bonds/elements)—each bond appears exactly once, so the hessian of the total energy is the direct sum of per-bond hessians with no correction. The 0.5 factor in non_bonded_hessian() undoes neighbor-list double-counting, which does not apply here.

jaxdem.utils.dynamicalMatrix.clump_non_bonded_hessian(state: State, system: System, cutoff: float | None = None, max_neighbors: int | None = None, rotation_scale: jax.Array | None = None) → tuple['State', 'System', jax.Array][source]#

Dense (n_clumps*group_dim, n_clumps*group_dim) clump hessian.

Generalized coordinates per clump are (δr_c, ω) – translation + small-rotation tangent (scalar in 2D, 3-vector in 3D). The hessian is accumulated by summing per-sphere-pair contributions via the rigid-body chain rule r_i = r_cI + p_i^lab; see _pair_clump_non_bonded_hessian_block().

Parameters:

rotation_scale (jax.Array, optional) – (n_clumps,) array of length scales (e.g., bounding-sphere radii) for the R·θ convention. When provided, the rotation rows/columns of clump I are divided by rotation_scale[I] on output. Default (None) returns the angle-based hessian (ω in radians).

jaxdem.utils.dynamicalMatrix.zero_mode_mask(eigenvalues: Array, rel_gap: float = 10000.0) → Array[source]#

Boolean mask identifying numerically-zero eigenvalues via gap detection.

Given a 1-D array of eigenvalues (in any order), we sort by \(|\lambda|\) ascending and look for the largest relative gap \(|\lambda_{k+1}| / |\lambda_k|\). Entries below the first gap that exceeds rel_gap are flagged as numerically zero. The mask is returned aligned with the original eigenvalue ordering, so it can be used directly to slice eigenvectors too (e.g. evecs[:, ~mask] for finite modes).

This is robust to problem scale: a hessian with |λ_max| ~ 1 and zero modes at ~ 1e-16 gives the same mask as one with |λ_max| ~ 1e6 and zero modes at ~ 1e-10, because the criterion is the ratio between successive magnitudes, not an absolute threshold.

Parameters:

  • eigenvalues (jax.Array) – 1-D array of eigenvalues (e.g. from jax.numpy.linalg.eigvalsh() or jax.numpy.linalg.eigh()). Ordering is not required.

  • rel_gap (float, optional) – Minimum ratio |λ_{k+1}| / |λ_k| that counts as the zero / finite boundary. Default 1e4 (four orders of magnitude). True zero modes at machine precision vs. real modes of order \(k \cdot \mathrm{overlap}\) in a jammed spring packing typically sit 10-14 orders of magnitude apart, so the threshold is not sensitive.

Returns:

Boolean array of the same shape as eigenvalues; True where the eigenvalue is below the first large relative gap. If no gap larger than rel_gap is found (all eigenvalues are comparable in magnitude), returns all-False.

Return type:

jax.Array