Hessian¶
The Hessian job is part of eOn Goswami [HS_Gos25]. The older description of that client is Chill et al. Chill et al. [HS_CWT+14].
The Hessian matrix (second derivative of the potential energy with respect to atomic coordinates) is used for:
Vibrational frequency analysis: Eigenvalues of the mass-weighted Hessian give squared vibrational frequencies. Positive eigenvalues correspond to stable modes; negative eigenvalues indicate saddle point character.
hTST prefactors: The harmonic transition state theory rate constant requires the product of frequencies at the minimum and saddle point (see prefactor).
Saddle verification: A first-order saddle point has exactly one negative Hessian eigenvalue.
How it works¶
eOn computes the Hessian numerically with a selectable finite-difference scheme
on the mobile (displaced) atoms: non-fixed atoms in pos.con, optionally
restricted by [Hessian] phva_atoms (comma-separated indices, or All). That
list is the hybrid/PHVA-class active set: atoms moved in FD. The remaining
atoms stay frozen (Li & Jensen, Theor. Chem. Acc. 107, 211, 2002).
Step size is Main.finite_difference (historically also called finite-difference
displacement; default (0.01,\text{Å})).
FD schemes ([Hessian] fd_scheme)¶
one_sided (default) — forward difference, (\sim M) force evaluations for (M = 3 N_\text{mobile}) active DOF (plus one base gradient):
$$H_{ij} \approx -\frac{F_j(x + h e_i) - F_j(x)}{h}$$
central — central difference, (\sim 2M) force evaluations:
$$H_{ij} \approx -\frac{F_j(x + h e_i) - F_j(x - h e_i)}{2h}$$
fourth — fourth-order central difference, (\sim 4M) force evaluations
(aliases fourth_order and central4):
$$H_{ij} \approx -\frac{-F_j(x+2h e_i)+8F_j(x+h e_i)-8F_j(x-h e_i)+F_j(x-2h e_i)}{12h}$$
One-sided is preferred for classical EAM AKMC cost; central reduces (O(h)) bias
when validating prefactors. Fourth-order cancels the (O(h^2)) term as well.
The same fd_scheme = fourth value selects that stencil for Lanczos and
Davidson Hessian-vector products. one_sided and central leave those
products on the forward difference. The rgsaddle min-mode
(min_mode_method = xtsci) takes forward differences for one_sided and
central differences for central and fourth; see Optimizer. Dimer
rotation, conjugate gradients, and LBFGS keep their own real steps. Complex-step is not a scheme: cutoffs and
external potentials are not holomorphic. The assembled matrix is symmetrized
((H+H^T)/2) before diagonalization. Mass-weighting uses
(\tilde H_{ij} = H_{ij}/\sqrt{m_i m_j}).
Hessian-range coloring¶
Coloring is exact only when every Hessian block between two atoms farther
apart than some distance is zero. A potential reports that distance as
Potential::hessianRange. The value is in the same length unit as the
positions. 0, or any non-finite value, means the range is unknown and
coloring stays off. One column is then built per coordinate.
An interaction cutoff is not this distance. A pair potential is zero past
its cutoff, so the range equals the cutoff. An atom-centred many-body
potential (EAM, Gupta, Stillinger-Weber, Tersoff) is zero only past twice
the cutoff, because the density or the bond order at an atom depends on
every neighbour inside the cutoff. A flexible water model does not have a
hard intramolecular extent, so a cutoff on the oxygen-oxygen interaction
does not make the blocks zero past any fixed distance. A message-passing
model is zero past n_layers times the edge cutoff, and only when that
layer count is part of the kernel.
Coloring is opt-in. RgpotAdapter reads HessianRangeOf<Kernel> and does
not look at config().cutoff. A kernel with a cutoff field and no
specialization stays at 0. The specialization returns the cutoff for a
pair kernel, twice the cutoff for an atom-centred many-body kernel, the
cutoff plus twice a hard intramolecular extent for a site model, or
n_layers times the edge cutoff for message passing. Fe-He is the
exception: its range is the larger of its pair cutoffs and twice its
largest density cutoff. Anything that is not exactly derivable returns 0.
A potential that is not an rgpot kernel opts in the same way, by
overriding hessianRange. An override of the deprecated
finiteCutoff() still supplies that distance. Everything else stays at 0.
The table of kernels eOn can wrap, at the current rgpot pin and at rgpot
v3.5.0, is in the developer page
Hessian range.
Mobile atoms are greedy-colored so that two atoms share a color only when
their closed neighborhoods of that range are disjoint. One displacement of
a whole color along x, y, or z fills every column of that color. Central
differences use two force evaluations per Cartesian direction per color,
and fourth-order differences use four, independent of the number of atoms.
Catch2 Colored FD Hessian matches serial central difference uses four
atoms on a line, spaced 1.2 Å, with a Lennard-Jones cutoff of 1.5 Å.
Coloring still reduces the force calls: three colors, so the central
scheme takes 1 + 2 × 3 × 3 = 19 force calls. The per-column loop takes
1 + 2 × 12 = 25. The colored matrix matches the serial central matrix
to 1e-8. Catch2 Embedding-chain coloring matches serial only at twice the cutoff is a planar zig-zag of eight atoms, bond 0.75 Å, with the
embedding cutoff at 1 Å. The central scheme takes 49 force calls total
and 43 once colored at 2 r_c. The fourth-order scheme takes 97
and 85. Coloring at r_c takes 19 and 37. The colored matrix matches
the per-column matrix to 1e-8. The short-range matrix misses by
0.1875 against max |H| 1.041667, and an off-axis entry reaches
0.541. The same central numbers hold with remove_net_force on and
one atom fixed, which is the path that still colors.
Catch2 Opted-in kernels match the serial matrix at the declared range
runs each opted-in kernel on its own zig-zag. The colored matrix matches
the per-column matrix exactly in these runs. A block past the interaction
cutoff is nonzero for every many-body kernel, and every block past the
declared range is 0. ZBL is a pair kernel, so that past-cutoff window is
empty. Its nearest-neighbour block at 1.5 Å is nonzero and every block
past cut_global is 0. The per-kernel call counts and the short-range
errors are in the developer page Hessian range.
Catch2 Coloring falls back when every mobile atom needs its own color
is four atoms inside a 20 Å Lennard-Jones cutoff, and the same cluster
with a batching spring of range 30 Å. Both need four colors. Both take
the column path, 25 force calls, which is one base gradient and two
central samples on 12 coordinates. The spring records a batch call, so
the colored loop did not run.
Catch2 Fourth-order stencil matches a complex-step oracle evaluates z^4 + 0.3 z^2 at 0.8. The fourth-order value stays
within 1e-10 of the complex-step derivative. The central value is more
than 1e-3 away, and farther from the oracle than the fourth-order
value. The same matrix is built one column at a time when the range is
unknown, when [Main] remove_net_force couples every atom, or when a
column checkpoint is in use.
When checkpoint_path is empty and the potential evaluates batches, the
displaced structures go out together, at most 32 columns per call.
RgpotPot takes that path once [RgpotPot] ranks_per_image has more than
one calculator group. Structure j runs on group j modulo the group
count, and every rank receives every result. A prefactor Hessian uses the
same call.
Column resume¶
Long partial Hessians can checkpoint FD columns to checkpoint_path (e.g.
hessian.ckpt) with resume = true. Interrupted jobs continue from the next
column; a successful run removes the checkpoint.
Usage¶
[Main]
job = hessian
finite_difference = 0.01
[Hessian]
phva_atoms = All
fd_scheme = one_sided
# fd_scheme = central
# fd_scheme = fourth
# resume = true
# checkpoint_path = hessian.ckpt
zero_freq_value = 1e-6
The output results.dat records force-call counts; hessian.dat contains the
mass-weighted matrix. [Main] quiet = true suppresses the write log line only.
Eigenvalues (squared frequencies) are obtained by diagonalizing the symmetrized
matrix (ColMajor eigen solve in the client).
The job drops each eigenvalue at or below zero_freq_value (default
1e-6). A free cluster has 6 of those eigenvalues, a linear cluster has
5, a periodic cell has 3, and a structure with any fully fixed atom has
none. The job still drops a different count, and it logs a warning. It
does not stop for that warning.
With [Hessian] write_modes = true (the default), modes.con holds one frame
per eigenvalue, lowest first. Each frame is the input structure with the mode
in its displacements section: the Cartesian displacement (x = q/\sqrt{m})
of the mass-weighted eigenvector (q), scaled to unit norm, and zero on atoms
outside the mobile set. The frame metadata carries mode_eigenvalue in
eV/(Ų amu), hbar_omega in eV and wavenumber in cm⁻¹. The last two are
negative for an imaginary mode. Any readcon reader gets the mode as an
N x 3 array beside the positions:
import readcon
for frame in readcon.read_con("modes.con"):
print(frame.metadata["wavenumber"], frame.disp.shape)
Free/fixed versus active (PHVA)¶
free/fixed on pos.con / Matter is the optimizer mask: which atoms may
move during geometry steps. Under an active-volume setup that set can be large
(the entire movable shell).
phva_atoms is the PHVA active set (Li & Jensen, Theor. Chem. Acc.
107, 211, 2002): atoms that are displaced when building a partial
Hessian or a matrix-free Hessian-vector product. It is always intersected with
free flags so frozen atoms are never moved. Default All means
active = free (historical behavior).
The same polarity applies to matrix-free min-mode. Free/fixed is left unchanged; an explicit list sets the Krylov / Ritz dimension to (3 N_\mathrm{active}) rather than (3 N_\mathrm{free}):
[Lanczos]
phva_atoms = 12, 15, 18, 21
tolerance = 0.01
max_iterations = 20
[Davidson]
phva_atoms = 12, 15, 18, 21
In pyeonclient:
import pyeonclient as pyec
mobile = pyec.resolve_mobile_atoms(matter, "12,15,18,21")
lanczos.compute(matter, direction, atoms=mobile)
# or: params.lanczos_phva_atoms = "12,15,18,21"
# lanczos.compute(matter, direction)
References¶
Samuel T Chill, Matthew Welborn, Rye Terrell, Liang Zhang, Jean-Claude Berthet, Andreas Pedersen, Hannes Jónsson, and Graeme Henkelman. EON: software for long time simulations of atomic scale systems. Modelling and Simulation in Materials Science and Engineering, 22(5):055002, July 2014. doi:10.1088/0965-0393/22/5/055002.
Rohit Goswami. Efficient exploration of chemical kinetics. October 2025. arXiv:2510.21368, doi:10.48550/arXiv.2510.21368.