mod model_hessian

module model_hessian

A model Hessian, and the quench depth it predicts for free.

The delayed-acceptance first stage needs the energy a relaxation will recover from an unrelaxed structure. That quantity has a closed form to leading order. Expanding about the minimum the descent will reach,

E(x) - E(Q(x)) ~ 1/2 g^T H^-1 g

with g the gradient at x and H the Hessian. Regressing a scalar proxy for that, which is what the surrogate did with |g|^2 / 2 lambda, collapses an anisotropic quantity onto one curvature and was measured to leave the second stage rejecting 65 per cent of what the first passed.

The true H costs second derivatives. A model Hessian costs nothing: the empirical pairwise force constants that quantum chemistry codes use to start a quasi-Newton optimiser, where the constant between two points falls off exponentially with their separation. Lindh and co-workers set them for exactly this purpose, and the same construction appears as the exponential preconditioner for atomistic optimisation. It is a function of the geometry alone, so evaluating it charges nothing on the ledger.

The operator

Stretch terms only, which is what a preconditioner needs and what a cluster of one species has:

H = sum_{i<j} k(r_ij) (u_ij u_ij^T) (x) [[1, -1], [-1, 1]]

with u_ij the unit vector along the pair and k the force constant. The form is a weighted graph Laplacian with a directional projector on each edge, so it is positive semidefinite and its nullspace is the three translations, which the solver projects out.

Applied matrix free, so a product costs O(N^2) arithmetic and no storage beyond the structure.

Why this is not a preconditioner here

It is the same object used differently. A preconditioner uses H^-1 g as a direction; this uses g^T H^-1 g as a number, the depth, and hands it to a decision about whether the real relaxation is worth paying for.

Variables

const ALPHA: f64

Decay of the force constant with separation, in units of the spacing.

A pair at twice the nearest-neighbour distance contributes about e^-3 of what a contacting pair does, which is the range over which a displacement is actually resisted.

const FLOOR: f64

Floor added to the diagonal, as a fraction of the contact force constant.

The stretch-only operator does not resist bending, so its nullspace holds every transverse motion as well as the rigid ones, and a gradient lying in it drives the solve to infinity: without this the depth of a chain pushed sideways came back as 5e19. The true Hessian has no such nullspace, because bending costs energy; the shift stands in for the angle terms the model omits, and is what makes the operator invertible rather than merely semidefinite.

const K0: f64

Force-constant scale, in the units the structure’s own spacing sets.

The magnitude does not matter for the decision the depth feeds, since the surrogate regresses on it and absorbs any constant, but the shape does: what carries information is that near pairs are stiff and far pairs are soft.

Functions

fn apply(x: ArrayView1<f64>, n: usize, v: ArrayView1<f64>, scale: f64) -> Array1<f64>

Applies the model Hessian to v.

Matrix free: each pair contributes a rank-one term along its own direction, so nothing is stored and a product is one pass over the pairs.

fn depth(x: ArrayView1<f64>, n: usize, gradient: ArrayView1<f64>, iters: usize) -> f64

Predicted energy a relaxation from x will recover, from the model Hessian.

Solves H d = g by conjugate gradients in the space orthogonal to the translations and returns 1/2 g^T d. Charges nothing: every operation is on the geometry and the gradient the caller already holds.

iters bounds the solve. A handful is enough, because the number feeds a decision rather than a step, and the surrogate that consumes it regresses away any systematic bias a truncated solve leaves.

fn spacing(x: ArrayView1<f64>, n: usize) -> f64

Mean nearest-neighbour distance, which sets the scale everything is in.