O1NumHess: Constant-Gradient Hessian Reconstruction

• O1NumHess is a seminumerical algorithm that reconstructs full molecular nuclear Hessians using a constant (O(1)) number of gradient evaluations.
• It exploits a local-plus-low-rank (ODLR) decomposition to compress the Hessian information, significantly reducing computational cost compared to traditional methods.
• The method is applicable to vibrational analysis, reaction-path following, and nonadiabatic dynamics, offering efficient scalability for large molecular systems.

O1NumHess is a seminumerical algorithm for reconstructing the full nuclear Hessian of a molecular system from a number of gradient evaluations that does not grow with system size, asymptotically $O(1)$ in the number of gradients. It was introduced as a fast and accurate alternative to conventional seminumerical Hessians, which require gradients at $O(N_{\mathrm{atom}})$ displaced geometries, and it relies on the off-diagonal low-rank (ODLR) structure of molecular Hessians: blocks associated with two distant groups of atoms have low rank rather than sparsity (Wang et al., 11 Aug 2025). In this formulation, O1NumHess is a specific Hessian-reconstruction method for molecular quantum chemistry, not merely a generic label for first-order numerical Hessians.

1. Problem setting and motivation

For nuclear coordinates $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$, with $n=3N_{\text{atom}}$, the target object is the nuclear Hessian

$$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$

This matrix is required for harmonic vibrational frequencies and normal modes, IR/Raman spectra, thermochemical quantities, reaction-path following and transition-state theory, and vibronic spectra and nonadiabatic dynamics via multiple-state Hessians (Wang et al., 11 Aug 2025).

The standard computational routes are analytic Hessians, numerical Hessians from energy finite differences, and seminumerical Hessians from finite differences of gradients. In the standard seminumerical scheme, one displaces one Cartesian coordinate at a time and evaluates

$$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$

which requires $6N_{\text{atom}}$ gradient evaluations in the double-sided variant or $3N_{\text{atom}}+1$ in the single-sided variant (Wang et al., 11 Aug 2025). O1NumHess was proposed to remove this linear growth in gradient calls while still recovering the full Hessian.

The method is motivated by two practical bottlenecks. First, gradients can be expensive for TDDFT, post-HF methods, and large systems. Second, analytic Hessians, although typically fastest and most accurate when available, are difficult to implement, memory-hungry because of CP-SCF with $3N_{\text{atom}}$ right-hand sides, and often unavailable for excited states and higher-level correlated methods (Wang et al., 11 Aug 2025). O1NumHess targets precisely this regime: full Hessian recovery when gradients are available but conventional seminumerical or analytic Hessians are unattractive.

2. Structural basis: locality, low rank, and information compression

The mathematical starting point is the linearized relation between displaced gradients and the Hessian. For displaced geometries

$$\boldsymbol{\xi}^{(j)} = \boldsymbol{\xi}^{(0)} + \Delta\boldsymbol{\xi}^{(j)},$$

the gradients satisfy, to linear order,

$O(N_{\mathrm{atom}})$0

After normalization by the Euclidean norm of each displacement, this becomes

$O(N_{\mathrm{atom}})$1

If $O(N_{\mathrm{atom}})$2 has full rank $O(N_{\mathrm{atom}})$3, then $O(N_{\mathrm{atom}})$4, which is the conventional picture and requires at least $O(N_{\mathrm{atom}})$5 displacements (Wang et al., 11 Aug 2025).

O1NumHess departs from this dense-matrix viewpoint by using a structural decomposition. The paper argues that the molecular Hessian can be approximated as

$O(N_{\mathrm{atom}})$6

The local part arises mainly from second derivatives of nuclear–nuclear interactions and local electronic contributions; the long-range part is dominated by a small set of near-degenerate electronic states and is therefore low rank. In block form, if atoms are partitioned into distant groups, the corresponding off-diagonal Hessian blocks are not assumed to be small, but they are assumed to have rapidly decaying singular values (Wang et al., 11 Aug 2025).

This shifts the effective number of degrees of freedom from $O(N_{\mathrm{atom}})$7 for a generic dense symmetric matrix to $O(N_{\mathrm{atom}})$8: a local block structure contributes $O(N_{\mathrm{atom}})$9 parameters, and a rank-$\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$0 symmetric correction $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$1 contributes $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$2 when $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$3 (Wang et al., 11 Aug 2025). Since each gradient is an $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$4-vector, a constant number of gradients then provides $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$5 scalar information, which is sufficient, in principle, to determine the independent parameters of the Hessian.

A common misconception is that O1NumHess assumes Hessian sparsity. The data support a narrower statement: the method assumes a local plus low-rank decomposition and, more specifically, ODLR structure. This distinction matters in systems such as conjugated polyenes or excited states, where the Hessian is not sparse but off-diagonal blocks remain numerically low rank (Wang et al., 11 Aug 2025).

3. Reconstruction algorithm

The algorithm has two principal reconstruction stages: a penalized local fit and an iterative low-rank correction.

The first stage estimates a local Hessian $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$6 by minimizing

$\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$7

subject to $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$8 (Wang et al., 11 Aug 2025). The Hadamard-weight matrix $\boldsymbol{\xi}=(\xi_1,\dots,\xi_n)$9 penalizes distant couplings according to

$n=3N_{\text{atom}}$0

with defaults $n=3N_{\text{atom}}$1 Bohr, $n=3N_{\text{atom}}$2, and $n=3N_{\text{atom}}$3 a.u. (Wang et al., 11 Aug 2025). Very distant pairs are set exactly to zero when

$n=3N_{\text{atom}}$4

This produces a sparse, symmetric positive definite linear system in the unknown non-far Hessian entries. Its dimension is $n=3N_{\text{atom}}$5, with $n=3N_{\text{atom}}$6, and it is solved iteratively using only matrix–vector products (Wang et al., 11 Aug 2025).

The second stage restores the missing long-range information by iterative low-rank updates. Gradients and displacements are first rescaled,

$n=3N_{\text{atom}}$7

with $n=3N_{\text{atom}}$8, რათა large forces are deemphasized and low-frequency information is emphasized (Wang et al., 11 Aug 2025). Starting from $n=3N_{\text{atom}}$9, the update is

$$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$0

followed by symmetrization,

$$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$1

Each correction has rank at most $$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$2, so the iteration explicitly enforces a local-plus-low-rank form (Wang et al., 11 Aug 2025).

The paper states that convergence is checked on the residual $$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$3, with stopping when the norm is below $$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$4 or its relative change per iteration is below $$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$5 (Wang et al., 11 Aug 2025). Because the true Hessian is only approximately ODLR and finite-difference or SCF noise is present, the residual is interpreted as a measure of numerical error rather than expected to vanish identically.

4. Displacement design and the meaning of “$$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$6 gradients”

The constant-gradient claim depends on a specific deterministic construction of collective displacement directions. O1NumHess begins from a model Hessian, $$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$7, built cheaply in internal coordinates using a modified Swart model and then transformed to Cartesian form (Wang et al., 11 Aug 2025). This model is used only to construct informative displacements.

For non-linear molecules, the initial displacement set contains seven directions: three translations, three rotations about principal axes of inertia, and one symmetric breathing mode (Wang et al., 11 Aug 2025). The translations require no gradients because the forces are zero in the absence of an external field. Rotational gradients at a non-equilibrium geometry are obtained analytically from the gradient at the reference geometry, so they do not require additional electronic-structure calls. The breathing mode is treated with double-sided finite differences because bond-stretch anharmonicity is strong; the remaining directions use single-sided differentiation (Wang et al., 11 Aug 2025).

Additional directions are generated through local completeness. For each atom $$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$8, the method defines a neighborhood $$H_{ij} \equiv \frac{\partial^2 E}{\partial \xi_i \partial \xi_j}.$$9 using the near-range cutoff. It projects out the subspace already spanned by current displacements restricted to that neighborhood, diagonalizes the residual local model Hessian, and takes the leading eigenvector as a raw local mode. Local modes from different neighborhoods are then combined into a global mode by choosing signs $$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$0 greedily to maximize the norm of the global displacement; the resulting vector is orthonormalized and scaled so that its maximum component equals the finite-difference step $$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$1 Bohr (Wang et al., 11 Aug 2025).

The critical asymptotic observation is that each neighborhood has bounded size for fixed cutoffs, so each local neighborhood requires only a bounded number of independent displacement directions. Once all neighborhoods are locally complete, adding more atoms does not require new global directions. The paper therefore states that, for sufficiently large systems,

$$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$2

and the total number of gradient evaluations, including equilibrium, breathing-mode doubling, and any post-processing directions, converges to a constant of order 100–120 (Wang et al., 11 Aug 2025).

This clarifies the phrase “$$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$3 gradients.” It does not mean a single gradient evaluation. It means that the number of gradients saturates at a system-size-independent constant for sufficiently large systems. The paper also states that O1NumHess “requires only about 100 gradients for sufficiently large systems” (Wang et al., 11 Aug 2025).

A further robustness step addresses spurious imaginary modes. After a first Hessian is built, the method diagonalizes it, selects a small number of negative-eigenvalue eigenvectors, orthogonalizes them against the displacement set, adds them as extra displacements, and recomputes the Hessian. The number of such additional directions is described as small, so the overall gradient count remains $$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$4 (Wang et al., 11 Aug 2025).

5. Accuracy, scaling, and implementations

The reported numerical behavior is benchmarked on long alkanes and polyenes, WCCR10 transition-metal reactions, and S30L-CI non-covalent complexes (Wang et al., 11 Aug 2025). The paper states that frequency, zero-point energy, enthalpy, and Gibbs free energy errors are only about two times those of conventional double-sided seminumerical Hessians. For large systems, O1NumHess is described as always faster than the conventional numerical Hessian algorithm, frequently even faster than the analytic Hessian, and it requires only about 100 gradients for sufficiently large systems (Wang et al., 11 Aug 2025).

For model systems, the reported vibrational-frequency mean absolute deviations are a few $$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$5 for well-behaved ground-state chains, with $$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$6-C$$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$7H$$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$8 giving $$H_{ij} \approx \frac{g_i(\boldsymbol{\xi}+\eta \mathbf{e}_j) - g_i(\boldsymbol{\xi}-\eta \mathbf{e}_j)}{2\eta},$$9–$6N_{\text{atom}}$0 depending on $6N_{\text{atom}}$1, and tens of $6N_{\text{atom}}$2 for challenging conjugated excited states such as C$6N_{\text{atom}}$3H$6N_{\text{atom}}$4 S$6N_{\text{atom}}$5 and T$6N_{\text{atom}}$6 (Wang et al., 11 Aug 2025). For WCCR10, the MADs of reaction enthalpies and Gibbs free energies are reported as about twice those of conventional double-sided seminumerical Hessians, with reaction $6N_{\text{atom}}$7 MAD $6N_{\text{atom}}$8. For S30L-CI, the complexation $6N_{\text{atom}}$9 MAD is reported as $3N_{\text{atom}}+1$0 (Wang et al., 11 Aug 2025).

In computational complexity terms, conventional double-sided seminumerical Hessians require $3N_{\text{atom}}+1$1 gradients, whereas O1NumHess saturates at roughly $3N_{\text{atom}}+1$2–$3N_{\text{atom}}+1$3 gradients for large systems (Wang et al., 11 Aug 2025). The reconstruction overhead is described as small relative to the gradient computations. The paper also reports that memory use is typically $3N_{\text{atom}}+1$4–$3N_{\text{atom}}+1$5 lower than for analytic Hessians on the tested systems, with smaller savings of $3N_{\text{atom}}+1$6–$3N_{\text{atom}}+1$7 in PCM cases (Wang et al., 11 Aug 2025).

The method is implemented in an open-source Python 3 core library, O1NumHess, together with a quantum-chemistry wrapper, O1NumHess_QC. The wrapper provides nuclear distance matrices, model Hessians, and interfaces to electronic-structure programs; the method is fully interfaced with BDF through the keyword O1NumHess in $$bdfopt</code>, and also interfaced with <a href="https://www.emergentmind.com/topics/orca-fda986d8-f1de-4492-b755-5a7c44800551" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">ORCA</a> through Python (<a href="/papers/2508.07544" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">Wang et al., 11 Aug 2025</a>). The core is described as method-agnostic: it assumes only a gradient function for arbitrary coordinate vectors and a distance matrix between coordinates, and it does not assume that coordinates are grouped in threes or that the target function is an electronic energy (<a href="/papers/2508.07544" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">Wang et al., 11 Aug 2025</a>).</p> <h2 class='paper-heading' id='interpretation-limitations-and-relation-to-other-hessian-methodologies'>6. Interpretation, limitations, and relation to other Hessian methodologies</h2> <p>The defining assumptions are spatial locality, low-rank off-diagonal blocks, and the availability of reliable gradients (<a href="/papers/2508.07544" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">Wang et al., 11 Aug 2025</a>). The paper explicitly notes several problematic regimes: very small systems, where constant overhead can dominate; highly delocalized systems such as metallic clusters or extended$$3N_{\text{atom}}+1$8-systems with dense low-lying excited-state manifolds, where ODLR ranks may be higher; pathological regions of the potential-energy surface such as conical intersections or very flat surfaces; and noisy gradients from poor SCF convergence or inadequate grids, whose errors may be amplified by the reconstruction (Wang et al., 11 Aug 2025).

The article’s terminology can be confused with another usage of “O1 numerical Hessian.” In matrix-based derivative-free optimization, an “order-1 numerical Hessian” denotes a Hessian approximation with $3N_{\text{atom}}+1$9 error in the sampling radius. The generalized simplex Hessian in “A matrix algebra approach to approximate Hessians” is such an order-1 numerical Hessian, constructed from function values alone, while its centered variant is order-2 (Hare et al., 2023). This is conceptually distinct from O1NumHess, which is a seminumerical molecular Hessian algorithm using gradient evaluations and ODLR reconstruction rather than a generic order-of-accuracy label.

This distinction suggests a broader classification of numerical Hessian methods. One family, exemplified by generalized simplex Hessians, uses only function evaluations and pseudoinverse-based matrix formulas, yielding order-$3N_{\text{atom}}$0 or order-$3N_{\text{atom}}$1 approximation guarantees for full or partial Hessians under rank and smoothness assumptions (Hare et al., 2023). O1NumHess belongs to a different family: it uses gradients at collective displacements, exploits problem-specific Hessian structure, and aims to reconstruct the full molecular Hessian with a constant number of gradient calls asymptotically (Wang et al., 11 Aug 2025). A plausible implication is that the name “O1NumHess” should be read historically as a project-specific algorithm name rather than as a statement about first-order truncation accuracy.

Within computational chemistry and beyond, the method is presented as generic for twice-differentiable multivariate functions whenever a distance between coordinates can be defined, far-coordinate Hessian blocks are numerically low rank, and gradients are available at arbitrary points (Wang et al., 11 Aug 2025). The paper specifically mentions classical atomistic potentials, PDE discretizations with integral operators of Calderón–Zygmund type, and machine-learning models with localized basis functions plus global pooling as plausible settings where analogous structure may arise (Wang et al., 11 Aug 2025). This suggests that O1NumHess is best understood not only as a molecular Hessian algorithm, but also as a structured inverse problem: recovering a full dense second-derivative operator from few directional gradient measurements under local-plus-low-rank priors.