07Moving nuclei and interpreting electron distribution

An energy is only one kind of answer

Our fixed-geometry nitrogen calculation answers an energy question at 1.85 Å. It does not yet tell us how a chosen analysis partitions the electrons between atoms, which way the nuclei would be pushed, or where a particular electronic model predicts an equilibrium bond length. Those questions require different quantities and, sometimes, a different calculation.

We will keep three operations separate. Population analysis interprets an existing electronic state. Nuclear differentiation asks how its fully specified energy model changes with nuclear coordinates. Geometry optimization repeatedly uses energies and derivatives to choose new coordinates. None of these operations automatically adds electron correlation or validates the electronic branch it receives.

All examples are unexecuted API illustrations. Source tests are evidence of intended implementation behavior, not calculations run for this lesson. If we relax the N₂ geometry, the result would be a new scientific case, not a replacement for the original fixed 1.85-Å benchmark.

Atomic populations are a convention for assigning shared density

The electronic density extends across the molecule; it does not arrive with a unique label saying which atom owns each electron. An atomic population assigns part of that density to an atom using a specified rule. A partial charge is the corresponding nuclear charge minus assigned electron population. Such charges are useful summaries, but are not unique atom-by-atom observables independent of basis and analysis method.

The actual public algorithm is population_analyzer. Both the native "qdk" implementation and optional "pyscf" implementation expose Mulliken analysis, selected by method="mulliken". The native schema allows only that method. A method name is not an implementation name: create("population_analyzer", "mulliken") is not the native factory contract. Do not assume that every population scheme offered by the external PySCF project is exposed through this adapter.

Mulliken analysis uses the AO density matrix $P$ and overlap matrix $S$. Let $A$ denote an atom and $\mu\in A$ denote basis functions assigned to that atom. Its electron population is

This allocates overlap-associated density between the atoms carrying the basis functions. Summing all atomic populations recovers $\operatorname{Tr}(PS)$, the total represented electron count. Summing only diagonal entries of $P$ would ignore the overlap convention and generally be wrong.

As a purely illustrative two-function example, assign one basis function to each of two centers and let

Each diagonal entry of $PS$ is one, so the two Mulliken populations sum to two. The trace of $P$ alone is only five-thirds. These are invented matrices showing the mathematics, not computed H₂ or N₂ data.

For an all-electron calculation the partial charge is $q_A=Z_A-N_A$, in units of the positive elementary charge. A positive value means fewer assigned electrons than nuclear charge. The molecular sum should agree with the specified total charge when all represented electrons are included. Individual Mulliken values need not behave like nonnegative probabilities, and they can be sensitive to diffuse functions or a changed basis. A large change after basis augmentation can be a partitioning effect rather than a large physical transfer of charge.

The implementation uses the actual one-electron density, including correlation

A correlated wavefunction need not have an AO density of the simple HF form $2C_{\mathrm{occ}}C_{\mathrm{occ}}^{\mathsf T}$. Its one-particle reduced density matrix, abbreviated 1-RDM, is a compact representation of one-electron information. In an MO basis we denote it by $\gamma$: diagonal entries are orbital populations, and off-diagonal entries retain coherence between orbitals. Its full construction from many-electron amplitudes belongs to Lesson 9.

For a real restricted orbital representation, the transformation to the AO density is $P=C\gamma C^{\mathsf T}$. Unrestricted calculations transform alpha and beta blocks separately and add their AO densities. Using only the diagonal occupations in a non-natural MO basis would discard information needed for the population analysis.

The native analyzer reads the active-space 1-RDM, embeds it at the actual active orbital indices, and restores inactive occupied orbitals. For ordinary restricted electron orbitals those inactive diagonal entries are two; for each unrestricted spin block they are one. It then forms the AO density, multiplies by the overlap matrix, and uses the basis’s AO-to-atom mapping to assign the diagonal contributions.

This supports more than the initial HF determinant, provided the required density information exists. It requires real-valued 1-RDM data in the native implementation, appropriate spin-resolved blocks for unrestricted orbitals, and molecular orbitals carrying both a basis/structure and overlap matrix. Incompatible active-space dimensions or missing data are explicit errors. For ModelOrbitals, the native implementation instead reports populations from the site-basis density diagonal; that is a lattice-model interpretation, not atom-centered Gaussian analysis.

The PySCF adapter similarly builds the density from the supplied wavefunction’s 1-RDM, calls PySCF’s mulliken_pop, and sums AO populations using PySCF’s atom slices. It returns electron populations, not the partial-charge array returned by a differently named external routine. The native tests specifically check overlap contraction, heteronuclear atom assignment, unrestricted density, correlated off-diagonal density, and inactive-core restoration.

For our existing reference, an ordinary API use would be:

from qdk_chemistry.algorithms import create

populations = create("population_analyzer", "qdk", method="mulliken").run(reference)
basis = reference.get_orbitals().get_basis_set()
effective_charges = basis.get_effective_nuclear_charges()
partial_charges = [float(z) - n for z, n in zip(effective_charges, populations)]

For the all-electron N₂ case, effective and physical nuclear charges coincide. A density respecting interchange of the two nitrogen atoms should give equal populations; with fourteen represented electrons, that symmetry would imply seven on each. This is a symmetry expectation, not a measured result here.

The effective-charge expression also handles ECP bookkeeping correctly: compare the explicitly represented electron population with $Z_A^{\mathrm{eff}}$, or equivalently restore removed core electrons before subtracting from physical $Z_A$. Subtracting only explicit electrons from a heavy atom’s unreduced physical charge would manufacture a spurious large positive charge. Frozen inactive orbitals and ECP-replaced core electrons are separate accounting issues.

The energy surface has nuclear coordinates as its arguments

Now regard the electronic energy, including the appropriate nuclear/core constant, as a function of all nuclear positions. With $N_{\mathrm{at}}$ atoms, collect the Cartesian coordinates into

A component of the nuclear gradient is $g_a=\partial E/\partial R_a$. It tells us how the energy changes for a small displacement in coordinate $a$. The corresponding mechanical force is $F_a=-g_a$: nuclei accelerate downhill in energy, not uphill. Coordinates are in Bohr, energies in Hartree, gradients in Hartree/Bohr.

The nuclear Hessian is

It is a $d$-by-$d$ matrix, in Hartree/Bohr², describing how one coordinate’s gradient changes when another coordinate moves. It is not the electronic orbital Hessian used in Lesson 6. The latter varied orbitals at fixed nuclei; this one varies nuclei with the electronic calculation defined at each geometry.

For our two-atom case the gradient has six entries and the Hessian is six by six. A diatomic bond scan can be described by one chosen bond-length coordinate, but the Cartesian derivative API still returns all six nuclear components. In an isolated, well-converged calculation translational invariance implies the net force should be nearly zero. This provides a diagnostic, not an excuse to discard a nonzero net force before investigating it.

The derivative objects preserve the geometry they refer to

NuclearGradients stores atom-major components: x, y, z for atom zero, then x, y, z for atom one. get_values() gives the flat vector, as_matrix() gives an atom-by-three matrix, and get_atom_gradient(i) retrieves one atom. It is a gradient object, so a forces array requires explicit negation.

NuclearHessian stores the matrix in the same coordinate order. get_matrix() returns it, and get_atom_pair_block(i, j) extracts the three-by-three block coupling atoms $i$ and $j$. Both objects retain a Structure identifying the geometry. Constructors reject a null structure or incorrect array shape; a correctly shaped matrix is not thereby a physically correct, converged Hessian.

Both support JSON and HDF5 with a serialization version, the associated structure, and explicit units. The gradient serializer writes hartree/bohr; the Hessian serializer writes hartree/bohr^2. File names follow the nuclear_gradients and nuclear_hessian type suffixes. No files are created in this preparation; this is the persistence contract for later use.

The computational interface is

energy, gradients, hessian, wavefunction = calculator.run(
    structure, charge=0, spin_multiplicity=1,
    seed_or_basis="cc-pvdz", n_inactive_orbitals=0,
)

hessian and wavefunction are optional outputs and may be None. Returning a gradient does not imply a Hessian was computed. These names and the keyword seed_or_basis come from the actual Python binding, rather than an assumed common name shared by every algorithm.

Analytic derivatives account for moving basis functions

An analytic gradient differentiates the mathematical energy expression and its integral ingredients. It does not mean an exact solution of the many-electron problem; an analytic HF gradient remains the derivative of an HF model.

At an accurately stationary variational solution, many first-order orbital-response terms can be eliminated through the stationarity conditions. But the Gaussian basis is attached to moving atoms. Moving a nucleus changes overlap and basis-dependent integrals as well as the explicit nuclear potential. The additional moving-basis contribution is commonly called the Pulay term. Keeping just electrostatic attraction derivatives and forgetting basis motion gives the wrong finite-basis force.

For a closed-shell SCF discussion, the gradient has the schematic structure

where $W$ is an energy-weighted occupied density matrix. The formula is a conceptual separation, not a complete replacement for the spin- and functional-dependent implementation.

The native gradient assembly includes nuclear repulsion, kinetic and attraction derivatives, overlap derivatives contracted with an energy-weighted quantity, ECP derivatives when present, electron-repulsion derivatives, and XC gradients for DFT. Libint provides ordinary Gaussian-integral derivatives; Libecpint provides ECP derivatives; GauXC supplies the native XC-gradient integration boundary. As established earlier, this native XC route uses ExchCXX built-in kernels rather than an assumed direct Libxc call.

The registered analytic implementation is create("nuclear_derivative_calculator", "qdk"), with alias "analytical_gradient". Its current wrapper requires its energy_calculator reference to be exactly scf_solver/qdk, requires n_inactive_orbitals=0, and rejects compute_hessian=True. There is no automatic numeric Hessian fallback inside that analytic implementation. It cannot be redirected to pyscf, qdk_stabilized, or a correlated calculator merely by changing the nested reference.

The wrapper uses native SCF’s analytic-gradient entry point and forwards its method settings. Representative native tests compare analytic and numeric GGA, hybrid-GGA, open-shell hybrid-GGA, and range-separated DFT gradients. Those establish concrete tested paths, not a blanket claim that every reference restriction or every external functional is validated. A returned derivative remains subject to the underlying SCF support and convergence requirements. Analytic second derivatives would require additional electronic-response machinery; this public analytic calculator does not supply them.

Finite differences trade extra electronic calculations for generality

The default derivative implementation is qdk_finite_difference, with alias "numeric". It approximates derivatives by evaluating the total energy at displaced structures. The step $h$ is a nuclear displacement in Bohr, not an SCF mixing factor.

For coordinate direction $a$, central differences give

where $e_a$ is a unit coordinate vector. The diagonal Hessian element is

For different coordinates $a,b$, the mixed element uses four doubly displaced geometries:

with, for example, $E_{+-}=E(R+h e_a-h e_b)$. The implementation computes these expressions directly from energies, reuses the single-coordinate energies for diagonal elements, and writes symmetric off-diagonal pairs. It is not a Hessian obtained by finite-differencing analytic gradients.

finite_difference_step defaults to 1e-3 Bohr with a declared allowed range of 1e-8 to 1.0. For a smooth energy surface the central formulas have truncation error of order $h^2$. But energy noise is amplified by approximately $1/h$ for gradients and $1/h^2$ for Hessians. Reducing the step can therefore make the result worse unless the electronic calculations are sufficiently accurate and stay on a consistent branch.

For illustration, an energy uncertainty scale of $10^{-8}$ Hartree with $h=10^{-3}$ Bohr can generate a gradient-noise scale around $10^{-5}$ Hartree/Bohr and a Hessian-noise scale around $10^{-2}$ Hartree/Bohr², with coefficients depending on the differences combined. These are order-of-magnitude arithmetic estimates, not error bars produced by the library.

The cost is explicit. With $d=3N_{\mathrm{at}}$, a gradient uses $1+2d$ energy evaluations, including the central energy. With a full Hessian this implementation uses $1+2d^2$ in total. For N₂, $d=6$, so those counts are 13 and 73. They count complete energy evaluations, not SCF iterations. Each displaced calculation can have its own convergence cost. The water derivative test checks shapes, finite values, near-zero total gradient, and symmetric numeric Hessian data without implying that every step choice is reliable.

Which energy methods can finite differences actually call?

The shared derivative settings contain an energy_calculator AlgorithmRef. Unlike many fixed-type nested settings, this reference is intentionally allowed to change algorithm type. Nevertheless, the energy-dispatch code handles specific families, not every possible registered algorithm.

For scf_solver, it directly calls the selected solver and returns its total energy and wavefunction. The default is qdk; an optional registered PySCF SCF implementation can be selected for the numeric route, subject to its own conventions. orbital_solver is ignored for this direct SCF path.

For multi_configuration_scf, the numeric route supplies reference orbitals and active alpha/beta counts to an MCSCF calculator, such as the available PySCF implementation. For multi_configuration_calculator, it first builds a Hamiltonian using the nested hamiltonian_constructor, then calls a calculator such as macis_cas or macis_asci. Those correlated methods are taught later. A coupled-cluster calculator is not silently supported merely because it also reports an energy: unhandled algorithm types raise an explicit error.

The default orbital_solver and hamiltonian_constructor are native qdk. Optional localize_reference_orbitals=True acts only on multireference paths and uses a nested localizer, default qdk_pipek_mezey, on the current active orbital indices. This is a representation-control step, not a replacement for choosing a scientifically appropriate active space.

n_inactive_orbitals subtracts that many doubly occupied orbitals from each spin population to determine active counts. The helper validates positive multiplicity, nonnegative explicit electron count after ECP accounting, compatible charge/spin parity, and that inactive count does not exceed either spin population. The multireference paths currently require positive active counts for both spins. Even when the selected SCF path does not use an active space, invalid inactive-count input is still rejected; for our first examples it is zero.

Seed behavior is especially important for finite differences. The central energy may use an appropriate Orbitals or Wavefunction seed. Displaced evaluations intentionally regenerate orbitals from a named basis rather than reuse coefficient matrices whose basis centers belong to the undisplaced geometry. reuse_seed_active_space=True preserves active/inactive indices for multireference evaluations; it does not preserve displaced coefficients or guarantee chemical orbital continuity. This will motivate Lesson 8.

A seed whose basis has the explicit custom_basis_set name is rejected for displaced numeric evaluations. Other nonstandard names must actually be resolvable when reloaded. A basis object carrying custom shells but masquerading as a standard name should not be assumed to retain those custom shells after displacement. With an optional SCF backend, numeric derivative validation also needs a basis representation the QDK counting path can understand, not only a name recognized by the external engine.

The default suppress_child_algorithm_logging=True temporarily raises the process-global native logger threshold during child energy evaluations. For convergence diagnosis, setting it false preserves the current child logging. That makes a large numeric calculation easier to diagnose without pretending a logging setting changes its physics.

First obtain derivatives without changing the benchmark geometry

A useful analytic-gradient request for the nitrogen model is:

from qdk_chemistry.algorithms import create
from qdk_chemistry.data import AlgorithmRef

analytic = create(
    "nuclear_derivative_calculator", "qdk",
    energy_calculator=AlgorithmRef(
        "scf_solver", "qdk", method="hf", convergence_threshold=1e-9,
    ),
)
energy_at_stretch, gradients, no_hessian, reference_at_stretch = analytic.run(
    structure, 0, 1, "cc-pvdz", 0,
)
forces = -gradients.as_matrix()

The tighter electronic tolerance is a deliberate numerical choice for derivative work, not a new HF approximation. forces would have two rows and three columns. No values or force directions are being reported here.

If a Hessian is required, choose a separate numeric calculator:

numeric = create(
    "nuclear_derivative_calculator", "qdk_finite_difference",
    compute_hessian=True, finite_difference_step=1e-3,
    energy_calculator=AlgorithmRef(
        "scf_solver", "qdk", method="hf", convergence_threshold=1e-9,
    ),
)
energy_check, numeric_gradients, nuclear_hessian, numeric_reference = numeric.run(
    structure, 0, 1, "cc-pvdz", 0,
)

These are complete requested operations, not guarantees that this step and tolerance suffice near every electronic instability. Compare numeric gradients with analytic ones where supported and repeat a step-size study before interpreting delicate curvatures. The release example shows the same separation between analytic gradients and numeric Hessians on its own small H₂ example.

Geometry optimization asks a different scientific question

A geometry optimization repeatedly chooses nuclear displacements to reduce the energy or locate a specified saddle point. This differs from a scan, where we prescribe a sequence of geometries or bond lengths and evaluate them. It also differs from dynamics, where forces, masses, initial conditions, and time integration describe motion. An optimization step is not a physical time step.

Near a candidate geometry, a local quadratic approximation combines the gradient with a curvature model. A minimum-search method moves in a downhill direction and limits the step when the local model is not trustworthy. A transition-state search instead seeks a first-order saddle: one meaningful downhill curvature direction and uphill curvature in the remaining internal directions. It therefore cannot be understood as merely “minimize harder.”

At a stationary geometry the nuclear gradient is small, but classification needs curvature. For an isolated nonlinear molecule, three translations and three rotations produce zero modes ideally; a linear molecule has three translations and two independent rotations. Thus linear N₂ has one internal vibrational mode after removing five rigid-motion directions. Small numerical deviations in those nominal zero modes must not automatically be interpreted as reaction-coordinate instabilities.

For harmonic vibrational interpretation, one forms the mass-weighted matrix $M^{-1/2}H^{\mathrm{nuc}}M^{-1/2}$, where $M$ contains each atomic mass repeated for its three coordinates. Its internal eigenvectors describe coordinated motions, and eigenvalues relate to squared harmonic frequencies after unit conversion. A negative internal eigenvalue is associated with an imaginary harmonic frequency. This is the conceptual next analysis of a Hessian, not a claim that the data object automatically computes normal modes, vibrational thermochemistry, tunneling, or reaction rates.

Geometry convergence must also be conditional on a well-defined electronic energy branch. If one displaced SCF point is a different spin-broken solution or an unconverged iterate, its difference from its neighbor need not be a derivative of one smooth surface. A stabilized electronic solver can help find locally stable branches, but switching branches can introduce its own discontinuities. Keep electronic stability and nuclear stationary-point classification separate.

geomeTRIC supplies the nuclear optimizer; QDK supplies the engine

The available implementation here is the optional geometry_optimizer/geometric plugin. The C++ factory declares "geometric" as its default name but has no built-in concrete optimizer registration; the plugin registers it when the external geometric package is available. A default name therefore does not prove availability on this environment.

The plugin’s optimizer setting defaults to "tric". TRIC means translation–rotation–internal coordinates: molecular or fragment translation and rotation are represented explicitly alongside internal bond/angle-like motion, improving the coordinate description for many optimizations. Internal coordinates describe molecular distortions more naturally than independently stepping every x/y/z entry. They change the optimization coordinates, not the underlying energy model.

The adapter accepts the forwarded coordinate-system names tric, tric-p, dlc, hdlc, prim, and cart. Delocalized internal coordinates combine primitive internal-coordinate directions; hybrid schemes combine coordinate descriptions; cart uses Cartesian coordinates. The important API distinction is that these are optimizer choices passed to geomeTRIC, not separately registered algorithms named geometric_tric or geometric_cart. Their detailed coordinate construction is external geomeTRIC functionality, not reimplemented in this repository.

The adapter constructs a geomeTRIC molecule whose stored geometry is in Ångstrom using geomeTRIC’s own bohr2ang constant. The custom engine receives optimizer coordinates in Bohr, constructs a new QDK Structure preserving element/mass/nuclear-charge information, and calls the configured derivative calculator. It returns a Hartree energy and atom-major Hartree/Bohr gradient to geomeTRIC. It does not negate the gradient into forces because the optimizer asks for an energy gradient.

Final geomeTRIC geometries are converted back from Ångstrom to Bohr using that same external constant; the last engine geometry can be used as a fallback. Using one consistent conversion constant at both boundaries avoids tiny artificial mismatches with cached evaluations. The plugin tests explicitly check the Bohr engine contract and final-coordinate conversion.

The engine retains the original optimization input on repeated calls; it is not a general orbital/basis transport service. A basis-name string, as in our example, avoids assuming an arbitrary original orbital seed or custom basis has automatically followed every new geometry. The selected derivative calculator must be able to rebuild the scientific model at those coordinates.

The geometry controls and their exact units

The shared settings use a nested derivative_calculator, default nuclear_derivative_calculator/qdk_finite_difference, max_iterations=300, and compute_hessian=False. The iteration cap has a declared range from one to one million. This cap counts geometry steps, not the electronic evaluations inside each derivative call.

The geomeTRIC extension supplies these default convergence thresholds:

RMS means root mean square, a combined measure of many components; the maximum criterion prevents one large component from being hidden by the average. By default geomeTRIC requires the five criteria together. Each exposed threshold has a declared range of zero to one in its stated units. transition_state=False selects minimization; setting it true forwards geomeTRIC’s transition-state search flag. print_level=0 controls external verbosity.

compute_hessian has a precise separate meaning. During geometry iterations the adapter sets the derivative calculator’s Hessian flag false. After optimization, if a final Hessian was requested it creates a fresh derivative calculator and evaluates it at the optimized structure. This final reporting option is not a request to use that Hessian at every optimizer step, nor an explicit initial transition-state Hessian input.

Consequently, selecting analytic qdk derivatives and also compute_hessian=True on the optimizer fails at the final derivative request: the analytic implementation rejects Hessians. Use numeric derivatives for that combined workflow, or optimize with analytic gradients and explicitly request a separate numeric Hessian afterward.

The return contract is (energy, optimized_structure, hessian, wavefunction), with the latter two optional. If a matching final engine evaluation is cached and no Hessian is requested, the adapter reuses its energy/wavefunction; otherwise it evaluates the final geometry again. It does not return a new detailed optimization-history object or an independent convergence certificate. Underlying derivative/optimizer errors remain important; a plausible final geometry alone does not establish a transition state. A missing coordinate payload can fall back to the engine’s last coordinates, so meaningful success interpretation remains tied to the upstream run and requested criteria.

A controlled relaxation fork, not a silent change of benchmark

Here is a source-backed configuration pattern adapted to our nitrogen variables. It requests analytic gradients during optimization and leaves the original structure untouched:

optimizer = create(
    "geometry_optimizer", "geometric",
    optimizer="tric", max_iterations=50, compute_hessian=False,
    derivative_calculator=AlgorithmRef(
        "nuclear_derivative_calculator", "qdk",
        energy_calculator=AlgorithmRef(
            "scf_solver", "qdk", method="hf", convergence_threshold=1e-9,
        ),
    ),
)
relaxed_energy, relaxed_structure, no_hessian, relaxed_reference = optimizer.run(
    structure, 0, 1, "cc-pvdz", 0,
)

This is not an executed result, a recommendation that RHF is adequate for stretched N₂, or a promise of the experimental equilibrium distance. It defines an optimization on a particular model and orbital branch. To classify the resulting stationary point, apply a supported nuclear Hessian calculation and electronic checks separately.

The repository’s geomeTRIC release example instead starts H₂ at 1.4 Bohr and uses numeric derivatives with a 1e-2-Bohr step. Its smoke test compares to an HF/STO-3G reference bond length. Those are a different molecule, basis, and starting geometry, not evidence for our nitrogen equilibrium. The ordinary external optimizer uses working files and logging infrastructure; none of it was invoked in this read-only course preparation.

What these outputs can and cannot justify

A Mulliken population can explain how one chosen convention partitions density and can help compare related calculations with consistent bases. A nuclear gradient identifies first-order forces on a chosen electronic energy surface. A Hessian describes local nuclear curvature. A converged geometry optimizer can locate a stationary structure within its model and search setup. None alone establishes a full reaction mechanism, a thermodynamic free energy, or a reaction rate.

For the original fixed-geometry QPE problem, we may deliberately keep 1.85 Å because the scientific objective is a difficult stretched-bond electronic energy, not equilibrium optimization. Population and force diagnostics are still useful, but they do not authorize relabeling a relaxed calculation as the same benchmark. The extra classical cost also belongs to any honest end-to-end workflow accounting: 73 energy evaluations for a numeric diatomic Hessian are not one electronic calculation simply because one high-level method returned the result.

We have also encountered a representation problem that will grow along a bond scan. Orbital indices can remain numerically identical while their chemical character changes. The next lesson develops canonical, natural, and localized orbital descriptions so those handoffs become chemically intelligible, without assuming that a rotation automatically preserves every truncated many-electron approximation.