Gradients

Overview

Analytic nuclear gradients of an integral are its first derivative with respect to the Cartesian coordinates of a single atom A, e.g.

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

returned as an extra length-3 trailing axis ($x,y,z$). For example:

julia> bset = BasisSet("sto-3g", """
              H 0.0 0.0 0.0
              H 0.7 0.0 0.0   
              """)  # Note that the `x` is the bonding axis
julia> ∇overlap(bset, 1) # Derivative w.r.t the position of the first atom
2×2×3 Array{Float64, 3}:
[:, :, 1] =
 0.0       0.347146
 0.347146  0.0

[:, :, 2] =
 0.0  0.0
 0.0  0.0

[:, :, 3] =
 0.0  0.0
 0.0  0.0

The first sub-array ([:, :, 1]) is the derivative with respect to the x position of the requested atom. In this simple example, displacements along y and z are rotations and do not affect overlap. Note also that the x sub-array is zero for shells on the same atom (such as diagonals).

Units: `R_A` is in bohr, not Angstrom

A BasisSet stores and accepts nuclear coordinates in Angstrom (see Quickstart), but every gradient/Hessian in this package is computed and returned in atomic units (bohr). To convert to a per-Angstrom derivative, divide by Molecules.bohr_to_angstrom (≈0.529177 Å/bohr), i.e. multiply by ≈1.8897 bohr/Å.

Overlap

The derivative of the overlap integral is defined as follows:

\[\frac{\partial S_{\mu\nu}}{\partial \mathbf{R}_A} = \frac{\partial\langle \chi_\mu | \chi_\nu \rangle}{\partial \mathbf{R}_A} \]

If the shells are on the same atoms ($\mu, \nu \in A$) or if neither shell is on the atom ($\mu, \nu \notin A$), the derivative is exacly zero.

Functions follow the same pattern found in One-Electron Integrals, except you must include a mandatory argument A which indicates the atom for each derivatives are evaluated.

julia> dr = ∇overlap(bset, 2, 1, 2) # Second atom, first shell, second shell
julia> dr[1,1,:] # Shows derivatives along x,y, and z
3-element Vector{Float64}:
 -0.3471463589514361
  0.0
  0.0
GaussianBasis.∇overlap — Function
∇overlap(BS::BasisSet, A) -> Array{Float64,3}
∇overlap(BS::BasisSet, A::Int, i::Int, j::Int) -> Array{Float64,3}

Gradient of the AO overlap matrix S w.r.t. atom A's three Cartesian coordinates, ∂S/∂R_A, with R_A in bohr (see Gradients for units).

Methods

  • ∇overlap(BS, A): full, dense nbas × nbas × 3 array.
  • ∇overlap(BS, A, i, j): just the (Ni,Nj,3) block for shells i,j of BS (shell indices, not AO indices).

For repeated calls, see ∇overlap!, which writes into a preallocated array instead of allocating.

source
GaussianBasis.∇overlap! — Function
∇overlap!(out, BS::BasisSet{LCint}, A::Int, i::Int, j::Int; scratch=nothing)
∇overlap!(out, BS::BasisSet, A)

Mutating counterpart of ∇overlap: writes into the caller-supplied out instead of allocating. The shell-pair form is the primitive the full-tensor form builds on – a direct libcint call, exactly zero (returned without any libcint call) whenever i,j are BOTH on atom A or BOTH off it: S_ij doesn't depend on R_A at all if neither shell sits there, and translating i,j together with atom A (both belong to it) leaves their relative separation – and so S_ij – unchanged.

Only implemented for the LCint backend – there is no ACSint fallback for gradients.

Methods

  • ∇overlap!(out, BS, A, i, j): out must be (Ni,Nj,3), the block for shells i,j of BS (shell indices, not AO indices).
  • ∇overlap!(out, BS, A): out must be a dense nbas × nbas × 3 array.

The shell-pair form needs a 3*Ni*Nj scratch vector whenever the derivative falls on shell j (see below). It allocates one per call unless you pass scratch; in a loop over shell pairs, hand it a buffer sized from 3*Nmax^2 to make the call allocation-free.

source

Kinetic

The derivative of the kinetic energy integral is defined as follows:

\[\frac{\partial T_{\mu\nu}}{\partial \mathbf{R}_A} = \langle \frac{\partial \chi_\mu}{\partial \mathbf{R}_A} | -\frac{1}{2}\nabla^2 | \chi_\nu \rangle + \langle \chi_\mu | -\frac{1}{2}\nabla^2 | \frac{\partial \chi_\nu}{\partial \mathbf{R}_A} \rangle\]

As with overlap, if the shells are on the same atom ($\mu, \nu \in A$) or if neither shell is on the atom ($\mu, \nu \notin A$), the derivative is exactly zero.

Functions follow the same pattern found in One-Electron Integrals, except you must include a mandatory argument A which indicates the atom for which derivatives are evaluated.

julia> dr = ∇kinetic(bset, 2, 1, 2) # Second atom, first shell, second shell
julia> dropdims(dr, dims=(1,2)) # Shows derivatives along x, y, and z
3-element Vector{Float64}:
 -0.33446292736778443
 -0.0
 -0.0
GaussianBasis.∇kinetic — Function
∇kinetic(BS::BasisSet, A) -> Array{Float64,3}
∇kinetic(BS::BasisSet, A::Int, i::Int, j::Int) -> Array{Float64,3}

Gradient of the AO kinetic energy matrix T w.r.t. atom A's three Cartesian coordinates, ∂T/∂R_A, with R_A in bohr (see Gradients for units).

Methods

  • ∇kinetic(BS, A): full, dense nbas × nbas × 3 array.
  • ∇kinetic(BS, A, i, j): just the (Ni,Nj,3) block for shells i,j of BS (shell indices, not AO indices).

For repeated calls, see ∇kinetic!, which writes into a preallocated array instead of allocating.

source
GaussianBasis.∇kinetic! — Function
∇kinetic!(out, BS::BasisSet{LCint}, A::Int, i::Int, j::Int; scratch=nothing)
∇kinetic!(out, BS::BasisSet, A)

Mutating counterpart of ∇kinetic: writes into the caller-supplied out instead of allocating. The shell-pair form is the primitive the full-tensor form builds on – see ∇overlap! for the identical structure and zero-block cases, T in place of S.

Only implemented for the LCint backend – there is no ACSint fallback for gradients.

Methods

  • ∇kinetic!(out, BS, A, i, j): out must be (Ni,Nj,3), the block for shells i,j of BS (shell indices, not AO indices).
  • ∇kinetic!(out, BS, A): out must be a dense nbas × nbas × 3 array.

As with ∇overlap!, the shell-pair form takes an optional scratch vector (>= 3*Ni*Nj elements) to stay allocation-free in a loop.

source

Nuclear

The nuclear attraction integral sums the potential of every nucleus $C$ in the molecule:

\[V_{\mu\nu} = -\sum_C Z_C \langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \chi_\nu \rangle\]

Its derivative w.r.t. atom $C$ has three contributing terms instead of two: the usual pair of shell-center ("bra"/"ket") derivatives from moving $\chi_\mu$/$\chi_nu$ summed over every nucleus C plus an extra term from the operator itself, since the C = A piece of the potential, $Z_A/|r-R_A|$, depends on $R_A$ too:

\[\frac{\partial V_{\mu\nu}}{\partial \mathbf{R}_A} = -\sum_C Z_C \left[ \langle \frac{\partial \chi_\mu}{\partial \mathbf{R}_A} | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \chi_\nu \rangle + \langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \frac{\partial \chi_\nu}{\partial \mathbf{R}_A} \rangle \right] - Z_A \langle \chi_\mu | \frac{\partial}{\partial \mathbf{R}_A}\frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle\]

${\mu,\nu} \in A$

For the case in which both shells are on the atom, we have

\[\frac{\partial V_{\mu\nu}}{\partial \mathbf{R}_A} = -\sum_{C\notin A} Z_C \left[ \langle \frac{\partial \chi_\mu}{\partial \mathbf{R}_A} | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \chi_\nu \rangle + \langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \frac{\partial \chi_\nu}{\partial \mathbf{R}_A} \rangle \right] - Z_A \frac{\partial}{\partial \mathbf{R}_A}\langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle\]

but translation invariance requires

\[\frac{\partial}{\partial \mathbf{R}_A}\langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle = 0\]

Hence,

\[\frac{\partial V_{\mu\nu}}{\partial \mathbf{R}_A} = -\sum_{C\notin A} Z_C \left[ \langle \frac{\partial \chi_\mu}{\partial \mathbf{R}_A} | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \chi_\nu \rangle + \langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \frac{\partial \chi_\nu}{\partial \mathbf{R}_A} \rangle \right]\]

${\mu,\nu} \notin A$

When no shell is at the atom, the formula simplifies to

\[\frac{\partial V_{\mu\nu}}{\partial \mathbf{R}_A} = - Z_A \langle \chi_\mu | \frac{\partial}{\partial \mathbf{R}_A}\frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle\]

Translational invariance still applies, but with a caveat:

\[\frac{\partial}{\partial \mathbf{R}_A}\langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle = 0 \\ \langle \frac{\partial \chi'_\mu}{\partial \mathbf{R}_A} | \frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle + \langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \frac{\partial \chi'_\nu}{\partial \mathbf{R}_A} \rangle + \langle \chi_\mu | \frac{\partial}{\partial \mathbf{R}_A}\frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle = 0\]

The trick here is that $\chi_\mu'$ and $\chi_\nu'$ are moving with the atom $A$. This is artificial, but computable. Hence:

\[\frac{\partial V_{\mu\nu}}{\partial \mathbf{R}_A} = Z_A \left[ \langle \frac{\partial \chi'_\mu}{\partial \mathbf{R}_A} | \frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle + \langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \frac{\partial \chi'_\nu}{\partial \mathbf{R}_A} \rangle\right]\]

$\mu \in A$ and $\nu \; \notin A$

\[\frac{\partial V_{\mu\nu}}{\partial \mathbf{R}_A} = -\sum_C Z_C \langle \frac{\partial \chi_\mu}{\partial \mathbf{R}_A} | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \chi_\nu \rangle - Z_A \langle \chi_\mu | \frac{\partial}{\partial \mathbf{R}_A}\frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle\]

Unlike overlap and kinetic, this is never exactly zero. The first two (shell-derivative) terms still vanish unless μ/ν sits on atom A – same translational-invariance argument as before – but the third (operator-derivative) term survives regardless of which atoms μ/ν sit on, since nucleus A's potential reaches every AO pair. For example, with shells 1 and 1 both sitting on atom 1 (not atom 2):

julia> dropdims(∇overlap(bset, 2, 1, 1), dims=(1,2))  # shell-derivative terms only -- exactly zero
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> dropdims(∇nuclear(bset, 2, 1, 1), dims=(1,2))  # operator-derivative term survives
3-element Vector{Float64}:
 0.36265343354910295
 0.0
 0.0

Reducing it to two libcint calls

The operator-derivative term above isn't something libcint can compute directly – cint1e_ipnuc_sph! only differentiates a shell, never the potential's own center. So instead the code rewrites the whole three-term expression as a sum of two shell-derivative pieces, split by whether each nucleus's charge is A's or not.

Split the sum at C = A:

\[V_{\mu\nu} = V^A_{\mu\nu} + V^{\neg A}_{\mu\nu}, \qquad V^A_{\mu\nu} \equiv -Z_A \langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_A|} | \chi_\nu \rangle, \qquad V^{\neg A}_{\mu\nu} \equiv -\sum_{C \notin A} Z_C \langle \chi_\mu | \frac{1}{|\mathbf{r}-\mathbf{R}_C|} | \chi_\nu \rangle\]

Define the two libcint-computable pieces. For a shell pair (p,q) and a charge set (either "just A" or "every nucleus except A"), one cint1e_ipnuc_sph! call (using a charge array with every other nucleus's charge zeroed out) gives the bra-shell electron-coordinate derivative of the very same quantity cint1e_nuc_sph! itself would return for that charge set – V^A_{pq} or V^{¬A}_{pq} from above, nothing separately signed:

\[K_A(p,q) \equiv \frac{\partial V^A_{pq}}{\partial \mathbf{r}_p}, \qquad K_{\neg A}(p,q) \equiv \frac{\partial V^{\neg A}_{pq}}{\partial \mathbf{r}_p}\]

This is the one rule every ip-prefixed libcint kernel follows, uniformly across ovlp, kin, and nuc: cint1e_ipXXX_sph! is the plain electron-coordinate derivative of whatever cint1e_XXX_sph! itself returns, no hidden extra sign. The translation identity ∂χ/∂R_p = -∂χ/∂r then converts this to the nuclear-coordinate derivative, exactly the same single flip used for overlap and kinetic:

\[\frac{\partial V^A_{pq}}{\partial \mathbf{R}_p} = -K_A(p,q), \qquad \frac{\partial V^{\neg A}_{pq}}{\partial \mathbf{R}_p} = -K_{\neg A}(p,q)\]

∇overlap_μ!/∇kinetic_μ! apply that same flip by hand (.*= -1.0) because S/T carry no leading sign of their own to fold it into; ∇nuclear_μ! doesn't need to, only because V^A/V^{¬A} already carry a leading minus (from Z being stored positive but the potential being attractive) for the flip to land inside of – not because cint1e_ipnuc_sph! is doing anything libcint's other kernels don't.

Piece 1 – V^{¬A}. R_A never appears in this operator, only possibly in the shell centers, so plain chain rule plus the relation above gives

\[\frac{\partial V^{\neg A}_{\mu\nu}}{\partial \mathbf{R}_A} = \frac{\partial V^{\neg A}_{\mu\nu}}{\partial \mathbf{R}_\mu}\Big|_{\mu \in A} + \frac{\partial V^{\neg A}_{\mu\nu}}{\partial \mathbf{R}_\nu}\Big|_{\nu \in A} = -K_{\neg A}(\mu,\nu)\Big|_{\mu \in A} \;-\; K_{\neg A}(\nu,\mu)^\top\Big|_{\nu \in A}\]

(each term present only if that shell sits on A.)

Piece 2 – V^A. Here R_A can appear in up to three places: the operator itself, and (if the shells sit there) R_μ, R_ν. Treat V^A(R_μ,R_ν,R_A) = -Z_A⟨χ_μ|1/|r-R_A||χ_ν⟩ as a function of three independent positions – shifting all three together leaves the integral unchanged (same argument as overlap's two-position version, one position more), so

\[\frac{\partial V^A_{\mu\nu}}{\partial \mathbf{R}_\mu} + \frac{\partial V^A_{\mu\nu}}{\partial \mathbf{R}_\nu} + \frac{\partial V^A_{\mu\nu}}{\partial \mathbf{R}_A} = 0\]

The full derivative with respect to atom A's position – dragging along whichever shells are attached to it, plus the operator's own dependence, substituted from the invariance above:

\[\frac{dV^A_{\mu\nu}}{d\mathbf{R}_A} = \underbrace{\frac{\partial V^A_{\mu\nu}}{\partial \mathbf{R}_\mu}\Big|_{\mu \in A}}_{\text{shell-drag, if any}} + \underbrace{\frac{\partial V^A_{\mu\nu}}{\partial \mathbf{R}_\nu}\Big|_{\nu \in A}}_{\text{shell-drag, if any}} + \underbrace{\left(-\frac{\partial V^A_{\mu\nu}}{\partial \mathbf{R}_\mu}-\frac{\partial V^A_{\mu\nu}}{\partial \mathbf{R}_\nu}\right)}_{\partial V^A_{\mu\nu}/\partial \mathbf{R}_A\text{, substituted}}\]

The shell-drag term and the substituted term cancel exactly whenever that shell is on A, and survive (sign-flipped) whenever it isn't:

\[\frac{dV^A_{\mu\nu}}{d\mathbf{R}_A} = K_A(\mu,\nu)\Big|_{\mu \notin A} \;+\; K_A(\nu,\mu)^\top\Big|_{\nu \notin A}\]

Add the two pieces. Writing them side by side gives the closed form ∇nuclear! actually evaluates:

\[\frac{\partial V_{\mu\nu}}{\partial \mathbf{R}_A} = f(\mu) + f(\nu)^\top, \qquad f(p) = \begin{cases} +K_A(p,q) & p \notin A \\ -K_{\neg A}(p,q) & p \in A \end{cases}\]

which needs exactly two libcint calls per shell pair (K_A and K_{¬A}, each also computed with p,q swapped and transposed for the f(ν) piece) – never three, and never a dedicated operator-derivative kernel.

Functions follow the same pattern found in One-Electron Integrals, except you must include a mandatory argument A which indicates the atom for which derivatives are evaluated.

julia> dr = ∇nuclear(bset, 2, 1, 2) # Second atom, first shell, second shell
julia> dropdims(dr, dims=(1,2)) # Shows derivatives along x, y, and z
3-element Vector{Float64}:
 0.9707070249006159
 0.0
 0.0
GaussianBasis.∇nuclear — Function
∇nuclear(BS::BasisSet, A) -> Array{Float64,3}
∇nuclear(BS::BasisSet, A::Int, i::Int, j::Int) -> Array{Float64,3}

Gradient of the AO nuclear attraction matrix V w.r.t. atom A's three Cartesian coordinates, ∂V/∂R_A (derivative of both the shell centers and the potential itself, since moving atom A also moves its nuclear charge), with R_A in bohr (see Gradients for units).

Methods

  • ∇nuclear(BS, A): full, dense nbas × nbas × 3 array.
  • ∇nuclear(BS, A, i, j): just the (Ni,Nj,3) block for shells i,j of BS (shell indices, not AO indices).

For repeated calls, see ∇nuclear!, which writes into a preallocated array instead of allocating.

source
GaussianBasis.∇nuclear! — Function
∇nuclear!(out, BS::BasisSet{LCint}, A::Int, i::Int, j::Int; scratch=nothing, charges=nothing)
∇nuclear!(out, BS::BasisSet, A)

Mutating counterpart of ∇nuclear: writes into the caller-supplied out instead of allocating. The shell-pair form is the primitive the full-tensor form builds on. Unlike ∇overlap!/∇kinetic!, no case is ever a free/skippable zero: V_ij sums the potential over every nucleus, so even a shell pair with neither i nor j on atom A still has a nonzero derivative through the Z_A/|r-R_A| operator term itself moving.

Every case reduces to the same two-term rule. Writing K_X(p,q) for ∇nuclear_μ! – the kernel differentiated w.r.t. its FIRST shell argument, over the nuclei that X leaves un-zeroed – and assigning each shell a (sign, charge set) of (+, no_A) when it sits on A and (-, only_A) when it does not:

out = s_i * K_{X_i}(i,j)  +  s_j * K_{X_j}(j,i)ᵀ

(transpose over the two AO axes, per Cartesian direction). Expanding the four on/off combinations recovers the familiar cases, e.g. i,j both on A gives K_notA(i,j) + K_notA(j,i)ᵀ, and both off gives -K_A(i,j) - K_A(j,i)ᵀ. Both terms are the same primitive, differing only in argument order and which charge set is passed.

Only implemented for the LCint backend – there is no ACSint fallback for gradients.

Methods

  • ∇nuclear!(out, BS, A, i, j): out must be (Ni,Nj,3), the block for shells i,j of BS (shell indices, not AO indices).
  • ∇nuclear!(out, BS, A): out must be a dense nbas × nbas × 3 array.

In a loop over shell pairs, pass charges (from nuclear_charge_sets, which depends only on BS/A) and a scratch vector of at least 3*Ni*Nj elements; the call is then allocation-free. Without them it rebuilds both charge arrays and a scratch buffer every time.

source

Two-electron four centers

GaussianBasis.∇ERI_2e4c — Method
∇ERI_2e4c(BS::BasisSet, iA) -> Array{Float64,5}

Gradient of the full two-electron four-center integral tensor (ij|kl) w.r.t. atom iA's three Cartesian coordinates, with R_iA in bohr (see Gradients for units). Returns a dense nbas × nbas × nbas × nbas × 3 array respecting the same 8-fold permutational symmetry as ERI_2e4c. This is the full, uncompressed tensor – for large basis sets prefer ∇sparseERI_2e4c, which screens and stores only the unique elements. For a single shell quartet, see ∇ERI_2e4c(BS,iA,i,j,k,l). For repeated calls, see ∇ERI_2e4c!.

source
GaussianBasis.∇sparseERI_2e4c — Function
∇sparseERI_2e4c(BS::BasisSet, iA, cutoff=1e-12; ij_vals=nothing, σvals=nothing)

Derivative (w.r.t. atom iA's three Cartesian directions, in bohr – see Gradients for units) of the unique (permutation-compressed) two-electron four-center integrals, Schwarz-screened the same way sparseERI_2e4c screens the energy integrals (see schwarz_bounds). ij_vals/σvals default to a fresh schwarz_bounds(BS) call if not supplied – pass a precomputed pair in when calling this repeatedly across atoms to avoid recomputing it every time.

source

Two-electron three centers

GaussianBasis.∇ERI_2e3c — Method
∇ERI_2e3c(BS1::BasisSet, BS2::BasisSet, iA) -> Array{Float64,4}

Gradient of the full two-electron three-center integral tensor (μν|P) (BS1=regular basis, BS2=auxiliary/fitting basis) w.r.t. atom iA's three Cartesian coordinates, with R_iA in bohr (see Gradients for units). iA indexes into BS1.atoms. Returns a dense BS1.nbas × BS1.nbas × BS2.nbas × 3 array, symmetric under μ↔ν swap. For repeated calls, see ∇ERI_2e3c!.

source
GaussianBasis.∇ERI_2e3c! — Function
∇ERI_2e3c!(out, BS1::BasisSet, BS2::BasisSet, iA; Bmerged=nothing)

Mutating counterpart of ∇ERI_2e3c: writes the dense BS1.nbas × BS1.nbas × BS2.nbas × 3 gradient into out instead of allocating it. out is zeroed first, so a reused buffer is safe.

Bmerged depends only on BS1/BS2, never on iA. Callers looping over atoms should build it once and pass it in – it is a small fraction of the runtime but over half of this routine's allocation.

source

Two-electron two centers

GaussianBasis.∇ERI_2e2c — Method
∇ERI_2e2c(BS::BasisSet, iA) -> Array{Float64,3}

Gradient of the full two-electron two-center integral matrix (P|Q) (the density-fitting Coulomb metric J_PQ) w.r.t. atom iA's three Cartesian coordinates, with R_iA in bohr (see Gradients for units). Returns a dense nbas × nbas × 3 array, symmetric under P↔Q swap. For repeated calls, see ∇ERI_2e2c!.

source
GaussianBasis.∇ERI_2e2c! — Method
∇ERI_2e2c!(out, BS::BasisSet, iA)

Mutating counterpart of ∇ERI_2e2c: writes the dense nbas × nbas × 3 gradient of the DF Coulomb metric (P|Q) into out instead of allocating it. out is zeroed first, so a reused buffer is safe.

source

Choosing a level

Every gradient here comes in the same four levels, from most convenient to fastest. All of them ultimately go through the same libcint primitive, so they agree to the last bit – they differ only in how much bookkeeping they do for you.

levelexamplenotes
dense, allocating∇overlap(bs, A)returns a fresh nbas × nbas × 3 array
dense, preallocated∇overlap!(out, bs, A)reuses out; zeroes it for you
shell pair/quartet∇overlap!(out, bs, A, i, j)validates sizes, resolves shell membership, returns the free zero when the block vanishes
bare primitive∇overlap_μ!(out, bs, i, j)one libcint call, no checks at all

Two things are worth knowing before dropping a level.

The dense drivers exploit symmetry that a per-pair loop must reproduce itself. S, T, V, (P|Q) and (μν|P) are all symmetric under the bra-shell swap, so the dense builders visit only i <= j and write the mirror from the same block. A loop that calls the shell-pair form for every ordered (i,j) gets the right answer at twice the libcint cost, with nothing to warn it. Loop i <= j and mirror, and the shell-pair form costs essentially nothing over the dense driver (measured ~1.0x).

Hoistable state should be hoisted. Several routines rebuild per-call state that depends only on the basis and the atom. Pass it in and the call becomes allocation-free:

routinekeywordbuild it with
∇overlap!, ∇kinetic!, ∇nuclear! (shell pair)scratchVector{Cdouble}(undef, 3*Nmax^2)
∇nuclear! (shell pair)chargesnuclear_charge_sets
∇ERI_2e3c!Bmergedmerged BasisSet of the two bases
∇sparseERI_2e4cij_vals, σvalsschwarz_bounds

Shell-pair/shell-quartet-level forms

For callers building an integral-direct gradient or CPHF loop that never wants the full dense array materialized:

GaussianBasis.∇ERI_2e4c — Method
∇ERI_2e4c(BS::BasisSet, iA::Int, i::Int, j::Int, k::Int, l::Int)
∇ERI_2e4c(BS::BasisSet, on_A::NTuple{4,Bool}, i::Int, j::Int, k::Int, l::Int)

Two forms. The iA::Int form is the convenient, standalone-safe one: computes on_A (which of shells i,j,k,l sit on atom iA, via === – see on_atom_flags) and returns the free zero (no libcint call) when i,j,k,l are ALL on atom iA or ALL off it. The on_A::NTuple{4,Bool} form is the unchecked core: it does whatever on_A says unconditionally, no membership computation and no zero-skip, for callers (like Fermi.jl's gradient loop) that have already screened at a higher level and precomputed on_A once outside a hot loop – calling this form on an all-on/all-off on_A does NOT raise an error, it just wastefully computes an answer of exactly zero the long way, so getting that screening right is entirely the caller's responsibility for this form.

source
GaussianBasis.∇ERI_2e4c! — Method
∇ERI_2e4c!(out, BS::BasisSet, on_A::NTuple{4,Bool}, i, j, k, l, buf)

Scratch-buffer-accepting core: identical math to the 7-argument form above, but takes a caller-owned buf (sized >= 3*Nmax^4, Nmax = the largest num_basis over any shell the caller will ever pass) instead of allocating one fresh every call – zero-allocation, for callers in a hot per-quartet loop. Not thread-safe to share: each concurrent caller (e.g. each worker task) needs its own buf.

source

Bare libcint primitives

The lowest level. Each does exactly one libcint call and applies the sign flip to the nuclear-coordinate convention – no bounds checking, no zero-block skipping, no output-size validation. Every one of them differentiates its first shell argument, so any other center is reached by passing that shell first and permuting the result yourself.

GaussianBasis.∇overlap_μ! — Function
∇overlap_μ!(out, BS::BasisSet{LCint}, i::Int, j::Int)

Libcint call for the overlap gradient block differentiated with respect to shell i (the "μ" AO), for shells i,j of BS, already sign-flipped to the nuclear-coordinate convention: libcint's cint1e_ipovlp_sph! computes ∂χ/∂r (w.r.t. the electron coordinate), but a GTO shell only depends on r and its center R through r-R, so ∂χ/∂R = -∂χ/∂r – this writes -cint1e_ipovlp_sph!(...), i.e. ∂S/∂R_i directly, matching every other ∇-prefixed function in this module.

No bounds checking, no zero-block skipping, no output-size validation. i,j must be valid shell indices (1:BS.nshells) and out must be exactly (Ni,Nj,3) – violating either is undefined behavior at the C level, not a catchable Julia error: an out-of-range shell index can segfault the whole process, and an undersized out can silently write past its end into unrelated memory. Prefer ∇overlap!/∇overlap unless you're writing your own integral-direct loop and have already validated i,j and out yourself.

source
GaussianBasis.∇nuclear_μ! — Function
∇nuclear_μ!(out, BS::BasisSet{LCint}, charge_atm, i::Int, j::Int)

Same as ∇overlap_μ! (same sign-flip reasoning, same lack of bounds checking), for a nuclear-attraction-type potential differentiated with respect to shell i (the "μ" AO) instead – except here the operator itself is a potential that can move with a nucleus, not a fixed multiplicative operator like 1 (overlap) or -∇²/2 (kinetic). charge_atm fixes WHICH potential: it's a libcint atom-description array (same layout as BS.lib.atm, 6 Cint slots per atom) with a chosen subset of nuclear charges zeroed out, so the raw libcint call sums Zc/|r-Rc| only over the nuclei left un-zeroed. Like ∇overlap_μ!/∇kinetic_μ!, libcint's cint1e_ipnuc_sph! computes ∂/∂r (w.r.t. the electron coordinate), so this writes -cint1e_ipnuc_sph!(...), i.e. ∂V^{charge}/∂R_i directly for the potential charge_atm encodes (see Nuclear for the derivation).

No bounds checking, no zero-block skipping, no output-size validation, and charge_atm isn't validated either – same segfault/heap-corruption risks as ∇overlap_μ! apply here too.

source
GaussianBasis.∇ERI_2e2c_μ! — Function
∇ERI_2e2c_μ!(out, BS::BasisSet{LCint}, i::Int, j::Int)

Libcint call for the 2-center (P|Q) gradient block differentiated with respect to shell i, sign-flipped to the nuclear-coordinate convention – the analogue of ∇overlap_μ! for the DF Coulomb metric.

Differentiates its FIRST shell argument, so Q is reached by swapping the two arguments, which returns a (Nj,Ni,3) block to be transposed in. out may be a contiguous view.

No bounds checking or output-size validation – same risks as ∇overlap_μ!.

source
GaussianBasis.∇ERI_2e3c_μ! — Function
∇ERI_2e3c_μ!(out, BS::BasisSet{LCint}, i::Int, j::Int, k::Int)
∇ERI_2e3c_P!(out, BS::BasisSet{LCint}, i::Int, j::Int, k::Int)

Libcint calls for the 3-center (μν|P) gradient block, sign-flipped to the nuclear-coordinate convention – the analogues of ∇overlap_μ! for this integral.

BS is the merged basis (regular shells followed by auxiliary ones, as ERI_2e3c!(out, BS, i, j, k) also expects), so an auxiliary shell k is addressed as k + BS1.nshells.

Two primitives are needed rather than one, because the three centers do not live in interchangeable slots: μ and ν are the two shells of the bra, while P sits alone in the ket. ∇ERI_2e3c_μ! (libcint's ip1) differentiates the FIRST shell argument, so ν is reached by swapping the first two arguments – the same trick used throughout this package. ∇ERI_2e3c_P! (libcint's ip2) differentiates the third center; no argument permutation can express it in terms of ip1.

out comes back in libcint's raw layout for the shell order as passed, and may be a contiguous view.

No bounds checking or output-size validation – same risks as ∇overlap_μ!.

source
GaussianBasis.∇ERI_2e4c_μ! — Function
∇ERI_2e4c_μ!(out, BS::BasisSet{LCint}, i::Int, j::Int, k::Int, l::Int)

Libcint call for the 4-center ERI gradient block differentiated with respect to shell i (the "μ" shell), for shells i,j,k,l of BS, already sign-flipped to the nuclear-coordinate convention – the direct analogue of ∇overlap_μ!, and the only place this file touches libcint.

Always differentiates its FIRST shell argument. To differentiate any other center, pass that shell first and permute the result, exactly as the 1e gradients do with ∇overlap_μ!(out, BS, j, i). out comes back in libcint's raw (Ni,Nj,Nk,Nl,3) layout for the shell order as passed, so a permuted call returns a permuted block – see ∇ERI_2e4c! for the four index maps that fold those permutations into the scatter.

Handles the 1-based to 0-based shell index conversion and passes the indices as an SVector, so no caller-owned shls buffer is needed and the call is allocation-free. out may be a contiguous view (e.g. view(buf, 1:3*Nijkl)).

No bounds checking, no zero-block skipping, no output-size validation – same segfault/heap-corruption risks as ∇overlap_μ!.

source
GaussianBasis.nuclear_charge_sets — Function
nuclear_charge_sets(BS::BasisSet{LCint}, A::Int) -> (only_A, no_A)

The two charge-fudged copies of BS.lib.atm that the nuclear gradient needs: only_A keeps atom A's nuclear charge and zeroes every other, no_A does the reverse. Passing either to ∇nuclear_μ! restricts the potential Zc/|r-Rc| it differentiates to that subset of nuclei.

Depends only on BS and A, never on the shell pair, so callers looping over shell pairs for one atom should build it once and pass it through ∇nuclear!'s charges keyword instead of paying for it per call.

source
GaussianBasis.schwarz_bounds — Function
schwarz_bounds(BS::BasisSet)

Per-shell-pair Cauchy-Schwarz screening bound σ_ij := sqrt(max|(ij|ij)|), used by both sparseERI_2e4c and ∇sparseERI_2e4c to skip shell quartets whose bound falls below a cutoff. Independent of which atom is being differentiated – callers making several ∇sparseERI_2e4c calls (e.g. looping over atoms) should compute it once and pass it through via the ij_vals/σvals kwargs rather than recomputing it each time.

source

Example

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);

julia> gS = ∇overlap(bset, 2);  # ∂S/∂R for the first H atom

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

julia> gS[2, 6, :]   # ⟨2s_O| shell | 1s_H(2)⟩ block, (x, y, z) components
3-element Vector{Float64}:
 -0.05625663685821103
 -0.0439525017729466
  0.0

julia> gERI = ∇ERI_2e4c(bset, 2);

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

For repeated calls, the mutating ∇overlap!, ∇kinetic!, ∇nuclear!, ∇ERI_2e4c!, ∇ERI_2e3c!, and ∇ERI_2e2c! forms write into a preallocated output array instead of reallocating.