Hessians

Analytic nuclear Hessians of an integral are its second derivative with respect to the Cartesian coordinates of two atoms $A$, $B$ (which may be the same atom), e.g.

\[\frac{\partial^2 S_{\mu\nu}}{\partial \mathbf{R}_A \partial \mathbf{R}_B}, \qquad \frac{\partial^2 (\mu\nu|\lambda\sigma)}{\partial \mathbf{R}_A \partial \mathbf{R}_B}\]

returned as two extra trailing length-3 axes ($x,y,z$ for $A$, then $B$). These feed analytic vibrational frequency calculations and second-order response (CPHF/CPKS) machinery.

Units: `R_A`/`R_B` are in bohr, not Angstrom

As with Gradients, a BasisSet stores and accepts nuclear coordinates in Angstrom, but every Hessian here is with respect to bohr displacement of each nucleus. So ∇2overlap(bset, iA, iB) gives ∂²S/∂R_iA∂R_iB per bohr², not per Angstrom². To convert to a per-Angstrom² derivative, divide by Molecules.bohr_to_angstrom^2 (≈0.529177 Å/bohr, squared).

Throughout this page the examples use

julia> using GaussianBasis, Molecules

julia> water = Molecules.parse_string("""
       O        0.000000000000     -0.143225816552      0.000000000000
       H        1.638036840407      1.136548822547     -0.000000000000
       H       -1.638036840407      1.136548822547     -0.000000000000
       """);

julia> bset = BasisSet("sto-3g", water);   # 5 shells, 7 basis functions

Overlap

\[\frac{\partial^2 S_{\mu\nu}}{\partial \mathbf{R}_A \partial \mathbf{R}_B}\]

A full dense Hessian is computed directly from a basis set and two atom indices:

julia> hS = ∇2overlap(bset, 1, 2);   # ∂²S/∂R_O ∂R_H(first)

julia> size(hS)
(7, 7, 3, 3)

An in-place (mutating) version writes into a preallocated array:

julia> out = zeros(bset.nbas, bset.nbas, 3, 3);

julia> ∇2overlap!(out, bset, 1, 2);

julia> out == hS
true

A per-shell-pair option is also available. This is usually what integral-direct and CPHF code wants, since it never materializes the full array:

julia> b = ∇2overlap(bset, 1, 2, 4, 3)   # 4 and 3 are shell indexes
1×3×3×3 Array{Float64, 4}:
[:, :, 1, 1] =
 -0.00951483  -0.0358454  0.0
...

For maximum efficiency, use the mutating shell-pair form with a preallocated output and a caller-owned scratch buffer, which makes the call allocation-free:

julia> Nmax = maximum(num_basis, bset.shells);

julia> scratch = Vector{Float64}(undef, 9*Nmax^2);

julia> out = zeros(num_basis(bset[4]), num_basis(bset[3]), 3, 3);

julia> ∇2overlap!(out, bset, 1, 2, 4, 3; scratch=scratch);

julia> out == b
true

GaussianBasis.∇2overlap — Method

∇2overlap(BS::BasisSet, iA, iB) -> Array{Float64,4}

Second derivative (Hessian) of the AO overlap matrix S w.r.t. atoms iA,iB's three Cartesian coordinates each, ∂²S/∂R_iA∂R_iB, with R_iA/R_iB in bohr (see Hessians for units). Returns a dense nbas × nbas × 3 × 3 array. For repeated calls, see ∇2overlap!.

level	example	notes
dense, allocating	`∇2overlap(bs, A, B)`	fresh `nbas × nbas × 3 × 3` array; not offered for the 4-center ERI
dense, preallocated	`∇2overlap!(out, bs, A, B)`	reuses `out`; zeroes it for you
shell pair/triple/quartet	`∇2overlap!(out, bs, A, B, i, j)`	validates sizes, resolves shell membership, returns the free zero when the block vanishes
bare primitive	`∇2overlap_μμ!(out, bs, i, j)`	one libcint call, no checks at all

routine	keyword	build it with
`∇2overlap!`, `∇2kinetic!`, `∇2ERI_2e2c!` (shell pair)	`scratch`	`Vector{Float64}(undef, 9*Nmax^2)`
`∇2ERI_2e3c!` (shell triple)	`scratch`, `Bmerged`	`merge_basis`
`∇2ERI_2e4c!` (shell quartet)	`buf` (positional)	`Vector{Float64}(undef, 9*Nmax^4)`

primitive	libcint	both/one derivative on
`_μμ!`	`ipip1`	both on the FIRST shell argument
`_PP!`	`ipip2`	both on the auxiliary shell
`_μν!`	`ipvip1`	one on each of the two bra shells
`_μP!`	`ip1ip2`	one on the first bra shell, one on the auxiliary

primitive	libcint	placement
`_μμ!`	`ipip1`	both derivatives on shell-argument position 1
`_μν!`	`ipvip1`	one on position 1, one on position 2 (same pair)
`_μλ!`	`ip1ip2`	one on position 1, one on position 3 (across pairs)