# Lesson 5 — Building the first electronic solution

## Turn the molecular model into something we can solve

We now have a specific problem: neutral singlet N₂, nuclei fixed 1.85 Å apart, initially represented in spherical cc-pVDZ functions. There are fourteen electrons and 28 one-electron basis functions. Our first calculation will optimize one determinant, not search all possible many-electron states.

That approximation is **Hartree–Fock**, or HF. Its orbitals adapt to the average field created by all occupied orbitals, including the exchange effect required by antisymmetry. Because that field depends on the orbitals we are trying to find, the numerical problem is circular. **Self-consistent field**, or SCF, is the process of resolving that circular dependence.

A future normal-environment call for the already-defined nitrogen `structure` is:

```python
from qdk_chemistry.algorithms import create

solver = create(
    "scf_solver", "qdk",
    method="hf",
    scf_type="auto",
    eri_method="direct",
    convergence_threshold=1e-7,
    max_iterations=50,
)
energy, reference = solver.run(structure, 0, 1, "cc-pvdz")
orbitals = reference.get_orbitals()
```

This is an annotated use of the actual API, not an executed calculation. The `qdk` implementation is a classical native electronic-structure solver. No quantum processor is involved in producing this reference. The same SCF interface also offers Kohn–Sham density functional theory, or KSDFT, which we will derive after understanding HF.

## The physical operator becomes integrals

In atomic units, the fixed-nucleus all-electron Hamiltonian contains electron kinetic energy, attraction to nuclei, electron–electron repulsion, and a constant nuclear repulsion:

\[
\hat H = \sum_i\left[-\frac12\nabla_i^2
 -\sum_A\frac{Z_A}{|\mathbf r_i-\mathbf R_A|}\right]
 +\sum_{i<j}\frac{1}{|\mathbf r_i-\mathbf r_j|}+E_{\mathrm{nuc}}.
\]

\(i,j\) label electrons; \(A\) labels a nucleus. \(\mathbf r_i\) is an electron position and \(\mathbf R_A\) a fixed nuclear position. \(Z_A\) is nuclear charge. The operator \(\nabla_i^2\), called the Laplacian, adds second spatial derivatives; the minus one-half factor gives electronic kinetic energy in atomic units. The nuclear term is attractive, so it has a minus sign. Electron repulsion is positive. The restriction \(i<j\) counts each electron pair once.

An **integral** turns one of these continuous-space expressions into a number for a particular pair or quartet of basis functions. We use real Gaussian basis functions here, matching the native SCF matrices. With \(\chi_\mu\) denoting AO basis function \(\mu\), the one-electron matrices are

\[
T_{\mu\nu}=\int\chi_\mu(\mathbf r)
 \left(-\frac12\nabla^2\right)\chi_\nu(\mathbf r)\,d\mathbf r,
\]

\[
V_{\mu\nu}=\int\chi_\mu(\mathbf r)
 \left[-\sum_A\frac{Z_A}{|\mathbf r-\mathbf R_A|}\right]
 \chi_\nu(\mathbf r)\,d\mathbf r,
\qquad h=T+V.
\]

\(h\) is the **one-electron core Hamiltonian matrix**. “Core” in this phrase means the kinetic-plus-external-potential operator, not “only the chemically inert core electrons.” All explicitly modeled electrons interact with it. The many-electron \(\hat H\) is not this small matrix.

The overlap matrix from Lesson 4 is \(S_{\mu\nu}=\int\chi_\mu\chi_\nu\). It is dimensionless; the energy integrals are in Hartree. The four-index electron-repulsion integral, using the **chemist ordering** defined by this equation, is

\[
(\mu\nu|\lambda\kappa)=
\iint \chi_\mu(\mathbf r_1)\chi_\nu(\mathbf r_1)
 \frac{1}{|\mathbf r_1-\mathbf r_2|}
 \chi_\lambda(\mathbf r_2)\chi_\kappa(\mathbf r_2)
 \,d\mathbf r_1\,d\mathbf r_2.
\]

The first two functions share one electron coordinate; the last two share the other. “Two-electron integral” describes this operator structure, not a simulation containing only two electrons. An entire molecule’s electron interaction is assembled from these entries.

**Libint** evaluates these Gaussian integrals through efficient recurrences. We do not put a dense three-dimensional mesh around nitrogen to evaluate ordinary Coulomb integrals by brute force. QDK Chemistry supplies normalized shell data, centers, operator choices, and screening policies; Libint supplies numerical integral blocks. The repository wraps it in [`OneBodyIntegral` engines](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/util/int1e.cpp) and electron-repulsion engines.

With an ECP, this same stage uses effective nuclear attraction plus an extra ECP operator: \(h=T+V_{\mathrm{eff}}+V_{\mathrm{ECP}}\). **Libecpint** evaluates the extra ECP matrix elements from the Gaussian shells and radial/angular potential data. The actual [native setup](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp) constructs `H_ = T + V + ECP`, replicated across spin-density blocks as needed. N₂ in the tutorial has no ECP; this extra term is zero there.

## Occupied orbitals produce a density matrix

For a closed-shell restricted HF reference, alpha and beta electrons occupy the same spatial orbitals in pairs. If \(C_{\mu i}\) is the coefficient of basis function \(\mu\) in occupied spatial orbital \(i\), the spin-summed density matrix is

\[
P_{\mu\nu}=2\sum_{i\in\mathrm{occupied}} C_{\mu i}C_{\nu i}.
\]

The factor two is essential: both spin projections occupy each spatial orbital. For our initial N₂ model, the coefficient matrix has 28 AO rows and, if no dependent directions are removed, 28 MO columns. Seven columns are occupied. Multiplying the 28-by-7 occupied block by its transpose gives the 28-by-28 density matrix.

The density matrix is not the many-electron wavefunction written as a 28-by-28 table. It compactly describes one-electron information. It generates the spatial electron density,

\[
\rho(\mathbf r)=\sum_{\mu\nu}P_{\mu\nu}
\chi_\mu(\mathbf r)\chi_\nu(\mathbf r).
\]

This density integrates to the electron count. In a nonorthogonal AO basis the corresponding matrix check is \(N=\operatorname{Tr}(PS)\), not generally \(\operatorname{Tr}(P)\). A **trace** adds the diagonal entries of a matrix. The overlap matrix is needed because our basis functions are not independent orthonormal bins.

For unrestricted calculations we keep separate \(P^\alpha\) and \(P^\beta\), each built without the factor two; their sum gives the total density. The [native density-update code](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp) implements precisely these occupancy factors. This is one of the conventions an adapter must get right: a factor-of-two error is physical, not just a formatting discrepancy.

## Coulomb and exchange build the mean field

The current density and electron-repulsion integrals define two important matrices. The **Coulomb matrix** is

\[
J[P]_{\mu\nu}=\sum_{\lambda\kappa}P_{\lambda\kappa}
(\mu\nu|\lambda\kappa).
\]

It describes the electrostatic field generated by the electron distribution. The **exchange matrix** uses a different pairing of basis indices:

\[
K[P]_{\mu\nu}=\sum_{\lambda\kappa}P_{\lambda\kappa}
(\mu\lambda|\nu\kappa).
\]

Exchange comes from the antisymmetric determinant, not a new classical force or literal swapping of electrons along trajectories. It contributes only between equal-spin electrons. Together with Coulomb it cancels the unphysical interaction of a spin orbital with itself in HF, but does not account for all correlated avoidance between different electrons.

For a restricted closed-shell calculation using the spin-summed \(P\), the **Fock matrix** is

\[
F[P]=h+J[P]-\frac12K[P].
\]

The one-half exchange factor reflects our spin-summed density convention. It is not present in the same form when spin densities are stored separately:

\[
F^\alpha=h+J[P^\alpha+P^\beta]-K[P^\alpha],
\qquad
F^\beta=h+J[P^\alpha+P^\beta]-K[P^\beta].
\]

Thus each spin sees Coulomb repulsion from both spins, but exchange from its own spin sector. The native [Fock update](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/scf_impl.cpp) follows these two cases: `J - 0.5*K` for restricted closed-shell density and the sum of spin Coulomb blocks minus the corresponding exchange block for spin-split density.

A Fock matrix is a density-dependent **effective one-electron operator**. Diagonalizing a 28-by-28 Fock matrix is not diagonalizing the enormously larger fourteen-electron Hamiltonian. That difference is the source of HF’s tractability and of its limitation.

## Orbital optimization is a constrained problem

We want a low-energy determinant, but cannot obtain it by multiplying every orbital by zero or allowing occupied orbitals to collapse onto one identical function. Orbitals must remain orthonormal. In AO coefficients, the constraint is

\[
C^{\mathsf T}SC=I,
\]

where the superscript \(\mathsf T\) means transpose and \(I\) is the identity matrix. For complex orbitals one would use a conjugate transpose; the traced native SCF path uses real coefficients.

Stationarity of the HF energy under this constraint leads to the generalized eigenvalue problem

\[
FC=SC\varepsilon.
\]

\(\varepsilon\) is a diagonal matrix of orbital energies. The equation is generalized because \(S\) is not the identity. It finds orbitals that are eigenvectors of the current mean field in the correct overlap metric. Since \(F\) itself depends on occupied columns of \(C\), solving it once is generally insufficient.

The actual numerical method transforms to an orthonormal representation. Write the overlap eigendecomposition as \(S=V s V^{\mathsf T}\), where \(s\) contains overlap eigenvalues. After removing directions below the linear-dependence threshold, construct \(X=V_{\mathrm{kept}}s_{\mathrm{kept}}^{-1/2}\). Then

\[
X^{\mathsf T}SX=I,\qquad F^{\prime}=X^{\mathsf T}FX,
\qquad F^{\prime}Y=Y\varepsilon,\qquad C=XY.
\]

Each column of \(Y\) is an ordinary eigenvector in the transformed basis. The native [eigenproblem implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp) uses this transformation and LAPACK’s symmetric eigensolver on the CPU path. **Eigen** supplies matrix representations and operations; **BLAS/LAPACK** supply optimized linear algebra. They do not choose the electronic approximation.

## The SCF loop as it actually happens

First the native adapter resolves the basis and spin case, constructs the internal solver, and builds overlap and one-electron integrals. Named or explicit-basis input selects its `Atom` initial-density path. This [initial-guess implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/guess.cpp) assembles atomic density blocks, using available atomic-guess files or generating them when absent. Generated guesses use [spherically averaged atomic HF](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/asahf.cpp): atomic angular degeneracies are treated without prematurely choosing a preferred bonding direction. This is an initial guess, not a statement that the molecule is a collection of noninteracting atoms.

An `Orbitals` input instead selects `UserProvided` density. The adapter takes occupied coefficient columns according to the requested electron populations. A restricted reference gives \(2C_{\mathrm{occ}}C_{\mathrm{occ}}^{\mathsf T}\); unrestricted input gives separate spin densities. This makes the restart mechanics concrete: we reuse orbital information to construct the starting field, then solve again.

At each outer iteration, the engine builds Coulomb and exchange contributions for the current density, updates the Fock matrix, evaluates total energy, and checks density and stationarity residuals. The numerical SCF algorithm then chooses improved orbitals and occupations. In the DIIS path it extrapolates a Fock matrix, solves the generalized eigenproblem, fills the prescribed lowest orbital columns, and constructs a new density. **Aufbau occupation** means this prescribed filling of the lowest orbitals in each spin channel; it is not a guarantee that we have discovered the globally correct electronic state.

The native default has `enable_gdm=True`, so later iterations may update orbitals through direct minimization instead of a fresh diagonalization at every step. It is therefore more accurate to say “update the occupied orbital subspace” than to insist all SCF implementations are only repeated diagonalization. After GDM, the native code rebuilds the field and updates canonical orbital eigenvalues.

The high-level native adapter exposes the choice between this atom-based start and an orbital-provided start through `basis_or_guess`; lower-level density-start options such as `Core` and `SOAD` exist, but are not additional arbitrary keys accepted by the high-level `qdk` settings. A source file containing an internal capability does not automatically make it a public workflow option.

The default iteration limit of fifty is the molecular outer-loop limit. Atomic-guess work and optimization substeps have their own mechanics; it is not fifty matrix multiplications or a wall-clock budget. When native SCF cannot satisfy its stopping criteria within that limit, it throws rather than returning an unconverged energy as a successful result.

## RHF, UHF, ROHF, RKS, and UKS are different restrictions

For a closed-shell molecule, **restricted HF (RHF)** uses one common set of spatial orbitals with paired occupations. **Unrestricted HF (UHF)** allows different alpha and beta orbital shapes. **Restricted open-shell HF (ROHF)** keeps one common spatial-orbital set but permits doubly occupied and singly occupied orbitals. These are different spaces over which the determinant is optimized.

The implementation name remains `"qdk"` or `"pyscf"`; RHF is not a required separate registry name. The method string and `scf_type` select the case:

| Request | Native HF / KSDFT | PySCF adapter HF / KSDFT |
|---|---|---|
| `auto`, singlet | RHF / RKS | RHF / RKS |
| `auto`, open shell | UHF / UKS | UHF / UKS |
| `restricted`, singlet | RHF / RKS | RHF / RKS |
| `restricted`, open shell | ROHF / rejected for DFT | ROHF / ROKS |
| `unrestricted` | UHF / UKS | UHF / UKS |

Here RKS and UKS are restricted and unrestricted Kohn–Sham orbital treatments. ROKS denotes the restricted open-shell Kohn–Sham route selected by this PySCF adapter; it is not a promise that all low-spin open-shell chemistry can be represented by one determinant.

The native [dispatch](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp) supports ROHF even though some overview documentation lists only RHF and UHF. The [native tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_scf.py) cover ROHF with DIIS, incore integrals, and GDM, and explicitly expect rejection of restricted open-shell DFT. The [PySCF tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_pyscf_plugin.py) include a lithium ROKS example. These capabilities were verified by reading implementation and tests, not by executing them here.

ROHF is not simply “average alpha and beta Fock matrices everywhere.” The native code partitions common orbitals into doubly occupied, singly occupied, and virtual blocks. It constructs an effective Fock matrix from block-specific combinations of the spin Focks: for example, doubly-to-singly occupied coupling uses the beta Fock, and singly-occupied-to-virtual coupling uses the alpha Fock. It then solves in the shared orbital space. At this revision that [ROHF construction](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp) requires equal AO and retained MO dimensions, so removal of linearly dependent AO directions can make this route fail.

UHF can lower the energy by letting opposite spins localize differently, which can help a stretched bond. But the determinant may no longer be an eigenstate of total spin. **Spin contamination** means its \(\langle\hat S^2\rangle\) differs from the intended \(S(S+1)\), in units where \(\hbar=1\). Equal alpha/beta counts fix a spin projection, not necessarily a pure singlet. Nor does requesting unrestricted orbitals automatically break symmetry: both adapters warn when a closed-shell unrestricted calculation is requested without a spin-broken orbital guess. Our default singlet nitrogen calculation is RHF, not secretly a broken-symmetry UHF calculation.

## Total energy is not a sum of orbital energies

Each occupied orbital energy includes interaction with the mean field of the other occupied electrons. Adding all orbital energies therefore counts pair interactions incorrectly. For the restricted HF density convention above, the correct expression is

\[
E_{\mathrm{HF}}=E_{\mathrm{nuc}}+
\frac12\sum_{\mu\nu}P_{\mu\nu}(h_{\mu\nu}+F_{\mu\nu}).
\]

Equivalently, the one-electron contribution is \(\sum P_{\mu\nu}h_{\mu\nu}\), and the interaction contribution is \(\tfrac12\sum P_{\mu\nu}(F_{\mu\nu}-h_{\mu\nu})\). The native engine records these separately before adding the nuclear constant. UHF and ROHF use the corresponding sum over spin blocks with their spin-specific densities and Focks.

A deliberately **toy algebra example**, not a molecule calculation, makes double counting visible. Take one normalized spatial basis function occupied by two opposite-spin electrons. Assign \(h=-1.0\) Hartree, \((11|11)=0.4\) Hartree, and a fixed nuclear constant of 0.1 Hartree. Then \(P=2\), \(J=K=0.8\), and \(F=-1+0.8-0.4=-0.6\) Hartree. The electronic energy is \(\tfrac12(2)(-1-0.6)=-1.6\) Hartree and the total is −1.5 Hartree. The occupied orbital-energy sum is only \(2(-0.6)=-1.2\) Hartree; adding the nuclear constant to that incorrect sum does not fix it.

The real API already returns total SCF energy. For the all-electron N₂ example, this includes the nuclear repulsion computed with physical charges. With an ECP it uses the effective-charge core–core expression and ECP one-electron terms described in Lesson 4. Do not add nuclear repulsion again, and do not add an imagined energy for removed core electrons. An ECP total energy belongs to its effective Hamiltonian’s reference convention.

## Kohn–Sham DFT changes the energy model

HF restricts the interacting wavefunction to one determinant. Kohn–Sham DFT takes another route: use auxiliary noninteracting orbitals to represent the electron density, and express the missing physics through an **exchange-correlation functional**. A functional takes a function, such as a spatial density, as input and returns a number.

For a simple spin-unpolarized, non-hybrid discussion, the energy decomposition is

\[
E_{\mathrm{KS}}[\rho]=T_s[\rho]
+\int v_{\mathrm{ext}}(\mathbf r)\rho(\mathbf r)\,d\mathbf r
+E_H[\rho]+E_{\mathrm{xc}}[\rho]+E_{\mathrm{nuc}},
\]

\[
E_H[\rho]=\frac12\iint
\frac{\rho(\mathbf r)\rho(\mathbf r^{\prime})}
{|\mathbf r-\mathbf r^{\prime}|}\,d\mathbf r\,d\mathbf r^{\prime}.
\]

\(T_s\) is the kinetic energy of the auxiliary orbital system, \(v_{\mathrm{ext}}\) is nuclear attraction, and \(E_H\) is the classical Coulomb energy. \(E_{\mathrm{xc}}\) includes exchange, electron correlation, and the remaining kinetic-energy difference not already described by this decomposition. Exact theory motivates the separation, but useful calculations require an approximation for \(E_{\mathrm{xc}}\).

For a local or semilocal non-hybrid functional, the effective matrix has the form \(F_{\mathrm{KS}}=h+J+V_{\mathrm{xc}}\). The XC potential comes from differentiating the chosen energy functional with respect to the density. Spin-polarized variants depend separately on alpha and beta densities. The SCF loop still updates density and field until consistent; changing from HF to DFT changes what is being made self-consistent, not the need for iteration.

A **local-density approximation** uses the density at each point. A **GGA**, or generalized-gradient approximation, also uses spatial density gradients; PBE is an example. A **meta-GGA** can use quantities such as kinetic-energy density built from orbital gradients. A **hybrid** mixes some HF-style nonlocal exchange with density-functional contributions; B3LYP and PBE0 are examples of that category. A **range-separated hybrid** changes the exchange mixing with electron separation; CAM-B3LYP is an implemented and tested example here. These labels distinguish approximations, not automatically increasing levels of guaranteed accuracy.

In both adapters, `method="hf"` selects HF. Any other method string is treated as an XC functional name. Use `method="pbe"` or `method="b3lyp"`, not a generic `method="dft"` plus an invented `functional` key. Native processing normalizes names and consults the supported functional map; unsupported names fail.

The native [KS energy code](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf/ks_impl.cpp) deliberately subtracts the XC potential contribution from its half-weighted Fock contraction and adds the separately integrated XC energy. A potential matrix is an energy derivative, not the XC energy itself. For hybrids, the exchange engine incorporates the functional’s exchange fractions; range-separated exchange can require both ordinary Coulomb and error-function-screened integrals. The [native range-separated tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/tests/test_scf.cpp) compare direct and incore CAM-B3LYP calculations.

A KS determinant is an auxiliary/reference object, not the exact correlated interacting wavefunction. In particular, evaluating the ordinary electronic Hamiltonian in that determinant does not generally reproduce the reported DFT energy. The native [Hamiltonian constructor](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/hamiltonian.cpp) constructs physical one- and two-electron integrals in the supplied orbitals; it does not automatically insert the chosen XC functional as a new many-electron interaction. This distinction matters when KS orbitals become a starting representation for later quantum calculations.

## Numerical quadrature and the actual XC library boundary

DFT introduces numerical integration over space. On an atom-centered grid, a continuous integral is approximated by a weighted sum,

\[
E_{\mathrm{xc}}\approx\sum_g w_g\,
 f_{\mathrm{xc}}\bigl(\rho(\mathbf r_g),\nabla\rho(\mathbf r_g),\ldots\bigr).
\]

\(g\) labels a grid point, \(\mathbf r_g\) its position, and \(w_g\) its integration weight. Here \(f_{\mathrm{xc}}\) denotes an energy density, so the formula does not require an extra density factor hidden in the notation. The additional arguments depend on the functional family. More accurate grids generally cost more basis evaluations and functional evaluations; grid convergence is a numerical issue distinct from choosing the physical XC approximation.

**GauXC** manages molecular grids, basis evaluation, integration work, and XC matrix/energy assembly in the native path. The current [input defaults](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/include/qdk/chemistry/scf/util/gauxc_util.h) are `UltraFineGrid`, Mura–Knowles radial quadrature, and `Unpruned` grids, with a basis-screening tolerance of `1e-10`. These are low-level defaults. The high-level native SCF settings do not expose PySCF’s `xc_grid` option or arbitrary GauXC configuration keywords.

A crucial repository-specific correction concerns **Libxc**. Libxc is a community library for evaluating XC functional energy densities and derivatives; it is not a complete SCF solver or a grid generator. Documentation credits it in the DFT ecosystem. However, this revision’s concrete [GauXC adapter](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/exc/GauXC/gauxc_impl.cpp) constructs `GauXC::functional_type` using **`ExchCXX::Backend::builtin`**, and the [dependency configuration](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/cmake/third_party.cmake) explicitly sets `EXCHCXX_ENABLE_LIBXC OFF`. We must not teach that every native DFT call directly runs Libxc. For this path, ExchCXX’s built-in functional implementation sits inside the GauXC integration machinery. A Libxc-style functional catalog in a resource JSON file is not proof that every entry is active in this compiled backend.

Density conventions are adapted at this boundary. Native RKS stores the total density, but the GauXC call used here expects the alpha half-density, so the adapter passes half of it. For UKS, the adapter forms total density and alpha-minus-beta density for GauXC, then converts the resulting potentials back to alpha/beta form. Those conversions illustrate why “both libraries accept a matrix” is an inadequate interoperability specification.

## Direct versus incore: choose where the cost lives

With \(K\) AO basis functions, a full four-index integral tensor has \(K^4\) entries before symmetry/storage reductions. The native setting `eri_method` selects `"direct"` or `"incore"`. Both describe integral evaluation/storage strategies, not different electron-correlation models.

**Direct**, the default, evaluates needed shell quartets while forming J and K and avoids retaining the entire four-index tensor. It repeats integral work as the density changes but saves large tensor storage. It still needs matrices, shell-pair data, engines, and sometimes per-thread accumulation buffers; “direct” does not mean memory-free.

**Incore** generates and stores integrals for later contractions. In this implementation the [stored arrays](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/INCORE/incore_impl.cpp) have four-index size, partitionable by the first index in the MPI design. They are not universally compressed to the smallest symmetry-packed form. Range-separated exchange may require an additional tensor. CPU builds retain them on the host; GPU-enabled branches can transfer storage to the device.

Illustrative single-tensor float64 storage, calculated only by multiplying counts by eight bytes, is about 4.69 MiB for 28 AOs and 98.88 MiB for 60 AOs. At 300 AOs it is 64.8 decimal GB. These figures are not total peak process memory, timings, or measured nitrogen resource use. They omit other arrays, scratch buffers, possible extra range-separated tensors, and parallel replication.

The direct [Libint2 implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/eri/LIBINT2_DIRECT/libint2_direct.cpp) avoids negligible work through multiple screening layers. It first rejects insignificant shell pairs by overlap. It uses Schwarz bounds and density-block magnitudes to reject quartets whose field contribution is too small. At an individual-integral level the relevant bound has the form

\[
|(\mu\nu|\lambda\kappa)|\le
\sqrt{(\mu\nu|\mu\nu)(\lambda\kappa|\lambda\kappa)}.
\]

The native high-level `eri_threshold` default is `-1.0`. Its actual dispatch rule uses a positive supplied value, otherwise chooses `convergence_threshold * 1e-5`, giving `1e-12` at the default SCF tolerance. `shell_pair_threshold` defaults to `1e-12`. These are screening controls, not certified final energy-error bounds. Making screening too loose can change the field or prevent agreement with a tighter reference.

The incore factory does not receive these direct-screening settings in the same way; it constructs its ERI object with a zero base threshold and generates its integral arrays through its own Libint route. Do not assume tweaking `eri_threshold` has identical semantics under both choices.

`eri_use_atomics=False` selects thread-local accumulation buffers in the direct engine; `True` uses shared-memory atomic updates. “Atomic” here is a parallel-programming term for an indivisible update, not an atomic-orbital approximation. Thread-local buffers trade memory for reduced contention and controlled accumulation. Floating-point reduction order can still affect last digits across thread counts or implementations.

The ordinary Fock build has quartic worst-case integral work, while dense eigensolution has cubic cost in matrix dimension. Screening changes practical cost substantially. The high-level route does not turn on density fitting merely because internal auxiliary-basis and `do_dfj` facilities exist; its default is `do_dfj=false`. Later factorization lessons will explain other representations rather than conflate them with this setting.

## Which stopping criterion did we actually request?

The scalar `convergence_threshold` is not a promised absolute energy error in Hartree. Native dispatch sets the orbital-gradient threshold to that value and the density threshold to one hundred times it. At the default `1e-7`, the latter is `1e-5`.

For ordinary RHF/UHF, the [convergence code](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/scf_algorithm.cpp) checks a density-change norm and the antisymmetric residual corresponding to

\[
R=FPS-SPF.
\]

A zero residual expresses that the occupied density and its field are mutually stationary in the nonorthogonal basis. The implementation normalizes the density norm and maximum-entry residual measure by the AO dimension and requires both to pass. ROHF uses its effective Fock/total-density construction. Energy changes are logged, but the native boolean stopping test is not simply “energy changed by less than the tolerance.”

The first few native Fock updates use the full density; subsequently the code can accumulate contributions from density differences. `fock_reset_steps`, default `1073741824`, controls periodic full rebuilds, not the SCF iteration budget. A default of \(2^{30}\) makes periodic resets effectively absent over an ordinary short run, although initial full builds still occur.

The default DIIS/GDM combination is selected by `enable_gdm=True`; setting it false selects DIIS alone. **DIIS** uses recent residuals and Fock matrices to extrapolate a better field. **GDM**, geometric direct minimization, searches over orbital rotations while preserving orthonormality. Its detailed geometry, line search, and stability analysis belong to Lesson 6.

One preview requires exact semantics: the [actual hybrid switch](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/scf_algorithm/diis_gdm.cpp) occurs when the DIIS step limit is exceeded, or after more than one step when the absolute energy change is **below** `energy_thresh_diis_switch`, default `1e-3`. Some prose/comments say “exceeds”; the conditional uses `<`. `gdm_max_diis_iteration` and `gdm_bfgs_history_size_limit` default to fifty and control the DIIS-to-GDM transition limit and retained optimization history respectively. `level_shift` is applied by DIIS only for a positive value; the default `-1.0` does not activate it. A documented-looking key such as `max_scf_steps` is not in the actual high-level native settings schema here.

Convergence remains distinct from stability. A stationary determinant can be a local minimum, a saddle point, or the wrong spin-constrained solution for the question. Numerical convergence also says nothing by itself about finite-basis or electron-correlation error.

## The PySCF adapter implements the same handoff, not identical numerics

The optional implementation is created as `create("scf_solver", "pyscf")`. The [adapter code](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/plugins/pyscf/scf_solver.py) constructs a PySCF `Mole`, passes the run’s charge, converts multiplicity to PySCF’s spin difference using `multiplicity - 1`, and sets `unit="Bohr"`. The atom-label helper prints coordinates read from `Structure`; a stale helper docstring saying Ångstrom does not describe this path.

The adapter chooses PySCF RHF, UHF, ROHF, RKS, UKS, or ROKS as shown earlier. For DFT it assigns the functional string to `mf.xc` and the integer `xc_grid` setting to `mf.grids.level`. This setting defaults to three, accepts zero through nine, and is a PySCF grid-density control; it is not numerically equivalent to GauXC’s `UltraFineGrid` label.

It maps `max_iterations` to `mf.max_cycle`. More subtly, it maps `convergence_threshold` to **`mf.conv_tol = 0.1 * convergence_threshold`**, an energy-change tolerance. It does not explicitly set `mf.conv_tol_grad` to the same number. Therefore identical QDK setting names across the two adapters do not establish identical stopping tests.

Without an orbital guess, `_run_scf` calls `mf.kernel()` and uses PySCF’s initial-guess behavior. With compatible QDK orbitals it reconstructs an initial density through the conversion utilities and calls `mf.kernel(dm0=dm0)`. The resulting orbital coefficients, energies, overlap, basis and ECP records are converted back into QDK objects.

The energy returned by `kernel` is passed through unchanged as the **SCF total energy**. The adapter neither subtracts nuclear repulsion nor replaces the passed charge/spin with settings values. Several docstrings claim otherwise; tracing `_build_scf`, `_run_scf`, and `_to_wavefunction` resolves the discrepancy. Do not copy those descriptions into a scientific analysis.

There is another important contract limit: this regular adapter does not explicitly inspect `mf.converged` before returning the tuple. PySCF can provide a final iterate without the same fail-on-nonconvergence behavior as the native wrapper. A returned object is therefore not uniformly a convergence certificate across backends. Logs, diagnostic behavior, and independent checks belong in the trustworthiness step, rather than assuming every adapter enforces identical success conditions.

Basis interoperability also remains bounded. Direct custom `BasisSet` input is explicitly unimplemented in `PyscfScfSolver`, even though separate conversion utilities support explicit shell data. The initial-guess route must match the basis, geometry and spin treatment. PySCF loads named bases and ECP mappings through its own ecosystem; equal basis-name strings alone are weaker evidence than matching the actual converted definitions and conventions.

PySCF is an external electronic-structure engine. Its internal integral, XC, and linear-algebra dependencies should be recorded from the installed PySCF environment rather than attributed to QDK’s native Libint/GauXC path. In particular, Libxc’s role in a PySCF installation is determined by that installation; this adapter does not select or verify a universal XC-library backend. The [package metadata](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/pyproject.toml) also makes PySCF optional and excludes the PySCF dependency from the Windows `plugins` extra. No such optional package was imported or installed for this lesson.

## Parallel software is another boundary, not another physical theory

Native Gaussian integrals, dense matrix algebra, and XC quadrature can all use classical parallel libraries. OpenMP distributes work across CPU threads where enabled; MPI supports distributed processes in compiled paths; BLAS may have its own threads. These are performance choices, not quantum gates.

The high-level native `nthreads` setting affects OpenMP only when compiled with OpenMP, and only a positive value requests a change. The default `-1` leaves the runtime’s current setting in place; it is not proof that every physical CPU core will be used. The [current build defaults](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/CMakeLists.txt) disable GPU and MPI by default and do not enable OpenMP in the Windows default configuration. Nothing in a source walkthrough determines which options an already-installed binary has.

Nested parallelism can oversubscribe a machine: many grid threads can each invoke a threaded matrix library. The [GauXC boundary](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/exc/GauXC/gauxc_impl.cpp) uses a [scoped BLAS-thread guard](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf/src/util/blas_threads.cpp) to request one BLAS thread during that work when the linked library’s control API is supported, then restore it. Performance and reproducibility checks must therefore include build and library context, not just the molecular input.

## What we have obtained, and what still needs to be trusted

The native adapter maps internal shells back into the public AO order, creates restricted or unrestricted `Orbitals`, and attaches a single canonical occupied `Configuration` through `StateVectorContainer` to a `Wavefunction`. In the HF case, the determinant has one configuration coefficient, while many MO coefficients were optimized. Those are entirely different quantities. The public `run` result is the energy/wavefunction pair, not the whole internal `SCFResult` containing every diagnostic field.

The output supplies a classical reference energy, orbitals for later transformation or integral construction, an initial many-electron reference, and a basis for classical diagnostics. It has not supplied correlation energies, a reliable dissociation curve, an active-space choice, a quantum circuit, or a hardware resource estimate. HF cannot represent all of the multiconfigurational structure that becomes important in stretched nitrogen. DFT incorporates approximated correlation through its functional but does not thereby provide a universally accurate ground-state wavefunction or a variational upper-bound guarantee.

The functional handoff is now clear: chemical input and BSE definitions become integrals through Libint/Libecpint; classical SCF and linear algebra optimize orbitals; KSDFT additionally uses functional kernels and GauXC or PySCF integration machinery; the resulting QDK objects feed later chemistry and quantum algorithms. Resource-estimation and error-correction specialists will eventually need the logical workload and error targets derived from that model. Better error correction cannot repair an inappropriate basis, a wrong spin state, or an unconverged classical reference.

Our next lesson therefore examines **trustworthiness of the first electronic solution**: how DIIS and GDM navigate the orbital problem, how failures are diagnosed, and why an SCF stability check asks a different question from an SCF convergence check.