Recover a sparse Hessian

A cutoff potential has a zero Hessian block between atoms that do not interact. sparse_hessian recovers those blocks directly. Build an atom-level SparsityPattern, colour it at distance 2, and read SparseHessian. The same value implements rgmin::ApplyHessian. to_dense is the dense export.

The colouring is greedy in smallest-last order. Two atoms share a colour only when they are not neighbours and share no neighbour, so one gradient difference recovers every block touched by that colour. Colours are independent. A block of them is one call to eval_displacements, and the rows of that call run together. After the batch returns, each colour is written back in colour index order. The block holds as many colours as SparseHessianOpts::memory_budget can store as geometry plus gradient buffers (256 MiB by default), and every colour when the whole problem fits. The budget is not the thread count. Central differences use both sides. Forward differences use one side and a single centre evaluation. The reported evaluation count is n_colours * 3 * sides, plus one for a forward centre. The displacement h is SparseHessianOpts::step. A central difference truncates as O(h²). A gradient accurate to eps_gadds roundoff O(eps_g/ h). Those terms meet near (3 eps_g/ |g'''|)^(1/3). On noise-free f64 gradients the Lennard-Jones and Gupta sweeps are smallest at 1e-6, and at 1e-4 they are 9.05e-7 and 8.13e-8. The default 1e-4 matches Index1Config’s fd_dr. An ECnoise estimate plugs in through the same field. Each geometry is its own chunk.

After recovery the blocks are replaced by (H + H^T) / 2. Powell and Toint (1979) cover both this symmetric direct read and the substitution. Star colouring (Coleman and Moré 1984) is the direct method that uses symmetry and stays a later option.

The pattern is symmetric and off-diagonal. The diagonal block is implied. from_positions bins atoms with linkcell’s cell primitives and keeps every pair that has an image inside the cutoff. The pinned linkcell crate exposes k-nearest queries, not a cutoff pair list, so the cutoff walk lives in this crate. An absent cell is an open cluster and the distance is Euclidean. Pass r_c for a pair potential, 2 r_c for an atom-centred model, and r_c times the layer count for a message-passing model, plus SparseHessianOpts::search_skin, which is max(2 * step, skin). That pad keeps an unmoved neighbour from entering the interaction range during a probe. A cell smaller than the cutoff does not store a self loop: a periodic image of an atom is that atom.

lowest_eigenpairs runs a thick-restart Lanczos on matvec. It does not call the surface. A Ritz pair is kept when ||H y - theta y|| is at most 1e-8 times an estimate of ||H||, and the residuals are returned with the pairs. The iteration cap is an error, not an unconverged vector. Inner products are a fixed-order dot: serial below 65,536 elements, and above that a parallel sum of 4,096-element chunks added back in chunk index order. Axpy updates go through rgmin::vecops. At dimension 96 or below the basis is the whole space. project_excluded removes the rows returned by PointSurface::excluded_modes. The recovered matrix is frozen at that geometry. apply_hessian and matvec panic if the vector length does not match. SparseHessianOpts and LowestEigenpairs are non-exhaustive. Build the opts with Default and with_difference, with_step, with_memory_budget, and with_skin.

one_atom_hessian is the per-atom baseline that eOn uses when a potential has no finite cutoff. eOn colours a finite cutoff itself. prefactor_hessian is the compressed seed. Index1Config is non-exhaustive. with_sparse_pattern passes that seed to Index1Session on the first step when no Hessian was supplied. The field defaults to none, so an index-1 search that does not set it still builds a dense central-difference Hessian. Index1Config::default() is unchanged.

The partition is Curtis, Powell, and Reid, IMA J. Appl. Math. 13, 117 (1974), https://doi.org/10.1093/imamat/13.1.117. The distance-2 rule is Gebremedhin, Manne, and Pothen, SIAM Rev. 47, 629 (2005), https://doi.org/10.1137/S0036144504444711. Powell and Toint, SIAM J. Numer. Anal. 16, 1060 (1979), https://doi.org/10.1137/0716078, cover the symmetric direct read and the substitution. Star colouring is Coleman and Moré, Math. Program. 28, 243 (1984), https://doi.org/10.1007/BF02612334.

let opts = SparseHessianOpts::default();
let pattern = SparsityPattern::from_positions(x.view(), Some(&cell), cutoff, opts.search_skin())?;
let hessian = sparse_hessian(&surface, x.view(), &pattern, opts)?;
let dense = hessian.to_dense();
let pairs = hessian.lowest_eigenpairs(1)?;
let _ = (dense, pairs.values, pairs.vectors, pairs.residuals);