HIP Hessians: Direct Learning of Hessian Matrices

• HIP Hessians are a framework that directly predicts full molecular Hessian matrices using SE(3)-equivariant graph neural networks without relying on derivatives.
• It constructs symmetric 3N×3N Hessians by assembling irreducible representation features, ensuring correct rotational transformation and physical symmetry.
• The method accelerates Hessian-dependent tasks, such as transition-state search and geometry optimization, with significant improvements in accuracy and computation speed.

HIP Hessians, in the sense of “Hessian Interatomic Potentials,” are a direct-learning framework for predicting molecular Hessian matrices from an SE(3)-equivariant graph neural network without using automatic differentiation, finite differences, or coupled-perturbed solvers (Burger et al., 25 Sep 2025). The method targets the full molecular Hessian needed in transition-state search, vibrational analysis, geometry optimization, and zero-point energy corrections, and does so by constructing equivariant, symmetric $3N\times 3N$ outputs from irreducible-representation features already present in message-passing architectures (Burger et al., 25 Sep 2025). Within the broader landscape of learned Hessians, HIP differs from approaches that recover Hessians by differentiating a learned energy surface, such as NewtonNet-based transition-state optimization, where the Hessian is obtained analytically from a $C^2$-continuous neural potential trained only on energies and forces (Yuan et al., 2024).

1. Concept, scope, and nomenclature

HIP is the acronym used for “Hessian Interatomic Potentials” in “Shoot from the HIP: Hessian Interatomic Potentials without derivatives” (Burger et al., 25 Sep 2025). Its central claim is that Hessians can be predicted directly from a deep learning model in one forward pass, rather than derived from an energy model by second-order differentiation. The paper motivates this by the observation that molecular Hessians are computationally expensive and scale poorly with system size under both quantum-mechanical methods and neural-network differentiation pipelines (Burger et al., 25 Sep 2025).

The target application domain is molecular simulation and electronic-structure workflows in which second derivatives of the potential energy are operationally central. The paper explicitly validates the approach on transition-state search, accelerated geometry optimization, zero-point energy corrections, and vibrational analysis benchmarks (Burger et al., 25 Sep 2025). This places HIP in the same problem family as recent work on learned Hessians for saddle-point optimization, where full Hessian information improves RS-PRFO-based transition-state searches relative to quasi-Newton updates (Yuan et al., 2024).

The terminology is potentially ambiguous. In another paper, “HIP Hessians” refers to local layer-wise Hessians of neural networks, defined blockwise with respect to each layer’s parameters and used as diagnostics of optimization geometry, overfitting, and generalization (Bolshim et al., 20 Oct 2025). In molecular ML, however, HIP denotes Hessian Interatomic Potentials (Burger et al., 25 Sep 2025). The two usages are unrelated beyond the shared emphasis on second-order structure.

2. Hessian object and symmetry constraints

HIP predicts the molecular Hessian

$$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$

where $I,J=1,\dots,N$ index atoms and $\alpha,\beta\in\{x,y,z\}$ index Cartesian components (Burger et al., 25 Sep 2025). The output is therefore a real symmetric $3N\times 3N$ matrix.

Two structural constraints organize the construction. First, the Hessian must satisfy symmetry,

$\mathbf{H}=\mathbf{H}^\top.$

Second, under a global rotation $\mathbf{Q}$, each atom-pair block must transform as a Cartesian rank-2 tensor,

$$\mathbf{H}_{I,J} \xrightarrow[]{\mathbf{Q}} \mathbf{Q}\,\mathbf{H}_{I,J}\,\mathbf{Q}^\top,$$

and at the full-matrix level

$$\mathbf{H} \xrightarrow[]{\mathbf{Q}} (\mathbf{I}_{N}\otimes \mathbf{Q})\,\mathbf{H}\,(\mathbf{I}_{N}\otimes \mathbf{Q})^\top.$$

The paper emphasizes that HIP enforces these properties by construction rather than through penalty terms (Burger et al., 25 Sep 2025).

A notable conceptual point is that Hessians are treated differently from forces. The paper states that Hessians do not need to be conservative in the same sense as forces; they only need to be symmetric and transform correctly. This is the rationale for predicting Hessians directly rather than deriving them from a scalar energy model (Burger et al., 25 Sep 2025). By contrast, differentiable-energy approaches infer Hessians from a smooth learned potential, as in NewtonNet, where smooth activations and a polynomial cutoff render the energy surface $C^2$0-continuous so that forces and Hessians can be obtained analytically by automatic differentiation (Yuan et al., 2024).

3. Equivariant construction from irreducible representations

HIP starts from an equivariant message-passing backbone; in the reported implementation the backbone is EquiformerV2 with four transformer layers (Burger et al., 25 Sep 2025). The key observation is representation-theoretic: if the backbone provides equivariant features with angular-momentum channels up to $C^2$1, then those features can be assembled into Cartesian $C^2$2 Hessian blocks while preserving SE(3) equivariance and matrix symmetry (Burger et al., 25 Sep 2025).

The readout begins by forming atom-pair features through message passing,

$C^2$3

These pair features are then projected to the compact irrep set $C^2$4,

$C^2$5

A Clebsch–Gordan tensor-product expansion reconstructs an intermediate $C^2$6 block,

$C^2$7

The final prediction is explicitly symmetrized,

$C^2$8

The paper identifies this final symmetrization as the mechanism enforcing $C^2$9 (Burger et al., 25 Sep 2025).

The requirement of $$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$0 features is not incidental. Because $$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$1 decomposes into $$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$2, irreps through $$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$3 are needed to span all components of a Cartesian $$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$4 tensor block. The paper states explicitly that irreps only up to $$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$5 would be insufficient to represent the $$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$6 content needed for a full Hessian readout (Burger et al., 25 Sep 2025). This is the mathematical basis of the direct blockwise construction.

4. Supervision, datasets, and loss design

HIP is trained on explicit Hessian labels rather than on energies and forces alone (Burger et al., 25 Sep 2025). The standard elementwise losses are

$$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$7

but the main experiments augment this with a subspace loss over the lowest eigenmodes,

$$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$8

leading to the combined objective

$$\mathbf{H}_{IJ}^{\alpha,\beta}=\frac{\partial^2 E}{\partial \mathbf{R}_I^{\alpha}\,\partial \mathbf{R}_J^{\beta}},$$9

The reported choice is $I,J=1,\dots,N$0, justified by the fact that the Hessian has 5 or 6 rigid-body zero modes and that the low-frequency subspace is crucial for ZPE, transition-state search, and frequency analysis (Burger et al., 25 Sep 2025).

The main training resource is the HORM dataset, described as the largest Hessian database used in the paper. It contains reactive molecular geometries with 5–30 atoms, sampled from Transition1x and RGD1, with Hessian, energy, and force labels recomputed at $I,J=1,\dots,N$1B97X/6-31G$I,J=1,\dots,N$2 (Burger et al., 25 Sep 2025). The training set includes 1,725,362 molecules from Transition1x and 60,000 from RGD1, with 50,844 Transition1x samples used for validation (Burger et al., 25 Sep 2025). To isolate Hessian prediction from energy/force-model differences, the backbone is frozen and only the Hessian head is trained, so HIP-EquiformerV2 and the AD-based baseline share the same underlying energy/force model (Burger et al., 25 Sep 2025).

Within the broader ecosystem of Hessian-supervised molecular learning, HIP is complemented by dataset-building efforts such as Hessian QM9, which provides 41,645 molecules with equilibrium geometries, vibrational modes, vibrational frequencies, and numerical Hessian matrices in vacuum, water, THF, and toluene at the $I,J=1,\dots,N$3B97X/6-31G$I,J=1,\dots,N$4 level (Williams et al., 2024). That work showed that including Hessians in the loss of an E(3)-equivariant message-passing GNN reduced vibrational-frequency MAE from $I,J=1,\dots,N$5 to $I,J=1,\dots,N$6 for motions above $I,J=1,\dots,N$7 (Williams et al., 2024). HIP takes a different step: it uses Hessian labels not merely to refine an energy model, but to learn Hessians directly (Burger et al., 25 Sep 2025).

5. Reported empirical performance

On the HORM-Transition1x validation set, HIP-EquiformerV2 is reported to outperform the AD-based EquiformerV2 baseline on all listed Hessian metrics while reducing per-molecule latency from $I,J=1,\dots,N$8 ms to $I,J=1,\dots,N$9 ms (Burger et al., 25 Sep 2025).

Metric	HIP-EquiformerV2	EquiformerV2 (AD)
Hessian MAE $\alpha,\beta\in\{x,y,z\}$0	0.030	0.074
Eigenvalue MAE $\alpha,\beta\in\{x,y,z\}$1	0.063	0.242
First-eigenvector cosine similarity	0.870	0.541
First-eigenvalue MAE $\alpha,\beta\in\{x,y,z\}$2	0.130	0.324
Time per molecule	38.5 ms	633.0 ms

The text summarizes these validation gains as “2x lower Hessian MAE, 3.5x lower eigenvalue MAE, and 1.5x higher eigenvector cosine similarity” (Burger et al., 25 Sep 2025). For speed and memory, the paper reports that direct Hessian prediction is $\alpha,\beta\in\{x,y,z\}$3–$\alpha,\beta\in\{x,y,z\}$4 faster than AD depending on molecule size, at least $\alpha,\beta\in\{x,y,z\}$5 faster in batched prediction, and $\alpha,\beta\in\{x,y,z\}$6–$\alpha,\beta\in\{x,y,z\}$7 lower in memory usage on small molecules with 5–30 atoms (Burger et al., 25 Sep 2025). The batched-AD disadvantage is attributed to the fact that AD effectively treats a batch as one larger problem, so Hessian memory scales quadratically with batch size (Burger et al., 25 Sep 2025).

Downstream results are similarly reported in Hessian-dependent tasks. For geometry optimization, using HIP Hessians in exact second-order RFO gives convergence in the fewest steps in the median case, and HIP-initialized BFGS performs best among the hybrid methods; finite-difference and AD Hessians are said to frequently fail to converge in this setting (Burger et al., 25 Sep 2025). For zero-point energy corrections, HIP-EquiformerV2 attains ZPE MAE $\alpha,\beta\in\{x,y,z\}$8 eV and $\alpha,\beta\in\{x,y,z\}$9ZPE MAE $3N\times 3N$0 eV, with standard deviations $3N\times 3N$1 and $3N\times 3N$2, respectively (Burger et al., 25 Sep 2025). For transition-state search in the ReactBench workflow, GSM success is identical between HIP and AD versions of the same backbone because GSM uses only energies and forces, but improvements appear in RS-P-RFO refinement and Hessian-based frequency analysis (Burger et al., 25 Sep 2025). In the main table, HIP-EquiformerV2 achieves 92% frequency-analysis accuracy for classifying extrema, compared with 71%–82% for the AD baselines (Burger et al., 25 Sep 2025).

The ablation on the loss function further emphasizes the role of low-mode supervision. The MAE+subspace loss yields 750 TS successes, 705 RFO convergences, and 698 cases passing both, versus 740/693/683 for plain MAE and 718/669/659 for MSE (Burger et al., 25 Sep 2025). For frequency analysis, MAE+subspace achieves 93% true positive rate, 6% false positive rate, 7% false negative rate, 94% true negative rate, and 92% overall accuracy (Burger et al., 25 Sep 2025). These results are aligned with the paper’s claim that the low-frequency subspace, rather than uniform elementwise accuracy alone, is decisive for many downstream chemical tasks.

6. Vibrational analysis, transition states, and relation to other Hessian-learning strategies

The vibrational-analysis pipeline reported for HIP removes the 5 or 6 rigid-body modes by mass weighting and Eckart projection before computing frequencies, ZPE, and the number of negative eigenvalues (Burger et al., 25 Sep 2025). The appendix details the construction: center-of-mass coordinates, mass weighting

$3N\times 3N$3

followed by translation and rotation vectors, orthonormalization, and projection onto the vibrational subspace (Burger et al., 25 Sep 2025). This makes the subspace loss operationally consistent with the evaluation procedures used for spectroscopy and transition-state diagnostics.

HIP is best understood relative to two neighboring strategies. The first is differentiable-energy learning, exemplified by NewtonNet for transition-state optimization, where the model is trained on energies and forces, uses smooth activations and a polynomial cutoff to remain $3N\times 3N$4-continuous, and supplies analytic Hessians by automatic differentiation at every optimization step (Yuan et al., 2024). That approach reduced optimization steps by about $3N\times 3N$5 relative to quasi-Newton optimization and reduced total gradient calls by nearly $3N\times 3N$6, while making Hessian evaluation at least three orders of magnitude faster than DFT analytical Hessians (Yuan et al., 2024). The second strategy is explicit Hessian supervision of an energy model, as in Hessian QM9, where second-derivative labels improve local PES curvature and vibrational frequencies across solvent environments (Williams et al., 2024). HIP departs from both by removing the requirement that Hessians be obtained from a learned scalar potential (Burger et al., 25 Sep 2025).

A plausible implication is that direct Hessian predictors and differentiable-energy Hessian generators occupy complementary niches. Direct prediction avoids the $3N\times 3N$7 Hessian-vector-product burden of AD Hessians and the $3N\times 3N$8 cost of analytic or finite-difference quantum-chemistry Hessians, while differentiable-energy models preserve exact derivative consistency with a learned scalar energy surface (Burger et al., 25 Sep 2025, Yuan et al., 2024). The available data do not claim equivalence between these design philosophies, but they do show that both can materially improve second-order chemistry workflows.

7. Implementation details, limitations, and open questions

The reported HIP implementation uses EquiformerV2 with four layers, 128 sphere channels, 4 attention heads, 64-dimensional attention channels, 128-dimensional FFN hidden size, SiLU activations, a Gaussian distance basis of size 512, cutoff radius 12 Å, cutoff Hessian 12 Å, and $3N\times 3N$9 spherical harmonics with $\mathbf{H}=\mathbf{H}^\top.$0; the Hessian head has three Hessian-specific layers (Burger et al., 25 Sep 2025). Training is performed on a single H100 GPU for roughly six days using AdamW with learning rate $\mathbf{H}=\mathbf{H}^\top.$1, batch size 128 for training and 256 for validation, gradient clipping at 0.1, StepLR decay with step size 10 and gamma 0.85, AMSGrad, zero weight decay, and 300 epochs with a limit of 1600 training batches per epoch, described as corresponding to about 35 full epochs (Burger et al., 25 Sep 2025). The transition-state workflows use pysisyphus for RS-RFO, RS-P-RFO, and IRC, and ReactBench for the benchmark setup (Burger et al., 25 Sep 2025).

The paper also states several limitations. HIP is trained only on relatively small organic molecules because DFT Hessian labels are expensive to obtain, and scaling to larger systems, materials, biomolecules, or proteins is left for future work (Burger et al., 25 Sep 2025). The reported models train only the Hessian head with a fixed backbone rather than in an end-to-end fashion, and the authors explicitly expect improvements from joint training of energies, forces, and Hessians, longer training, and additional hyperparameter tuning (Burger et al., 25 Sep 2025). Locality is built into the head through an interaction cutoff, so the method relies on the empirically and theoretically supported sparsity of molecular Hessians arising from localized electronic interactions (Burger et al., 25 Sep 2025). Within the cutoff the pairwise construction has $\mathbf{H}=\mathbf{H}^\top.$2 memory and compute, but the paper states that for larger systems the effective scaling becomes closer to $\mathbf{H}=\mathbf{H}^\top.$3 because of sparsity (Burger et al., 25 Sep 2025).

The codebase and model weights are reported as open sourced at https://github.com/BurgerAndreas/hip (Burger et al., 25 Sep 2025). This operational detail matters because HIP is presented not only as a modeling idea but as a concrete replacement for derivative-based Hessian generation in downstream molecular workflows. Within the evidence provided, its defining proposition is that Hessians need not remain a secondary quantity derived from energies: they can be the primary prediction target of an equivariant molecular model (Burger et al., 25 Sep 2025).