# Lesson 3 — How scientific workflows become software

## The calculation is a chain of scientific handoffs

We have chosen neutral, singlet nitrogen at a bond length of 1.85 Å. We also know that a molecular geometry, a set of orbitals, a many-electron state, and a Hamiltonian are different things. The software architecture makes those distinctions useful: each stage receives a particular scientific object and produces another one. A later stage should not need to know which internal loops produced the earlier object.

Imagine handing another colleague a wavefunction. A spreadsheet of unexplained numbers is insufficient. They need the orbital basis, electron occupations, coefficient conventions, and a way to determine which information is actually present. QDK Chemistry packages these together in data objects. An algorithm object is the separately configured processor that creates or transforms them.

A **class** describes a kind of object: `Structure` describes molecular structural data. An **object**, or instance, is one particular structure, such as our stretched nitrogen molecule. A **function** is a named operation with inputs and a result. A **method** is a function attached to an object or class: `wavefunction.get_orbitals()` retrieves the orbitals associated with that wavefunction. An **interface** is the agreement about which operations exist and what they accept and return. A **backend** is the implementation doing the work behind that agreement. These programming words matter because they tell us what can be exchanged without rewriting the science workflow.

The [architecture guide](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/docs/source/user/comprehensive/design/index.rst) expresses the intended separation as immutable data and configurable, stateless algorithms. Treat this as a design principle, not a claim that every operation is free of caches, randomness, external services, or implementation-specific behavior.

## Four decisions that must not be confused

Consider the first electronic calculation. Choosing **Hartree–Fock** selects a scientific approximation: represent the electrons by a single determinant and optimize its orbitals. Choosing **cc-pVDZ** selects the finite one-electron representation in which those orbitals are expanded. Choosing an SCF optimization strategy such as **DIIS** or the native **DIIS/GDM combination** selects a numerical route toward a stationary solution. Choosing implementation **`qdk`** rather than **`pyscf`** selects software that carries out the task.

Those choices interact, but they are not synonyms. A different backend can implement the same Hartree–Fock model; a different basis changes the model space even when the backend is unchanged. Two solvers with the same interface may use different convergence tests, integration grids, initial guesses, or orbital orderings. Interface compatibility means their objects fit the same handoff. It does not certify identical energies, spin behavior, numerical precision, or scientific conclusions.

Here is the data journey we will gradually make operational:

| Handoff | Scientific object being passed |
|---|---|
| Molecular input → SCF | `Structure`, charge, spin multiplicity, and a basis or orbital guess |
| SCF → classical analysis | A total energy and `Wavefunction`, with associated `Orbitals` |
| Orbital/model selection → integral construction | `Orbitals`, including the basis and relevant orbital-space information |
| Integral construction → many-electron or qubit methods | `Hamiltonian`, rather than just the SCF energy |
| Fermion-to-qubit translation | `QubitOperator`, carrying an encoding convention |
| Trial-state construction | `Circuit` for preparing a selected `Wavefunction` |
| Hamiltonian-based circuit construction | A `UnitaryRepresentation`, then circuit representations |
| Phase estimation and execution | A `QpeResult`; the circuit and workload also feed resource analysis |

This is not one compulsory straight-line script. Classical solvers and quantum solvers are alternative consumers of a model Hamiltonian. State preparation and Hamiltonian compilation are two branches that must agree when joined. Resource analysis can inspect circuits without pretending they were executed on a quantum device. Active-space selection, correlated wavefunctions, encodings, and execution deserve their own later lessons.

The concrete [Hamiltonian-constructor contract](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/include/qdk/chemistry/algorithms/hamiltonian.hpp) accepts `Orbitals` and returns `Hamiltonian`. The [unitary-builder interface](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/hamiltonian_unitary_builder/base.py) accepts `QubitOperator` and returns `UnitaryRepresentation`. The [phase-estimation interface](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/phase_estimation/base.py) accepts a preparation `Circuit` and `QubitOperator`, with optional noise information, and returns `QpeResult`. Notice that a common verb, `run`, does not imply common arguments.

## Asking the registry for the right tool

A **registry** is a catalog mapping names to implementations. A **factory** consults that catalog and constructs a new algorithm object. The Python entry point combines the two ideas:

```python
from qdk_chemistry.algorithms import (
    available, create, inspect_settings, show_default,
)

scf_names = available("scf_solver")
default_scf = show_default("scf_solver")
schema = inspect_settings("scf_solver", "qdk")
solver = create("scf_solver", "qdk", max_iterations=80)
```

These examples explain API use in an ordinary installed environment; no imports or calculations were executed while preparing this lesson. `available()` with no argument returns a dictionary of algorithm types and their available implementations. With `"scf_solver"`, it returns a list. `show_default` reports the configured default, not a benchmark-based recommendation. At this revision, `ScfSolverFactory.default_algorithm_name()` is `"qdk"`. The native factory also registers `"qdk_stabilized"`, a composite stability-checking workflow that we will study in Lesson 6. PySCF names appear only when that integration is successfully registered.

`create("scf_solver")`, `create("scf_solver", None)`, and an empty implementation name select the factory default. The default is explicitly declared: it is not whichever implementation happens to appear first in a list. `inspect_settings` returns records containing a setting name, expected Python type, default value, description, and constraints where present. `print_settings` provides a human-readable table. Discovery is safer than copying an old example table whose keys may no longer match the implementation.

The source of truth is the [Python registry](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/registry.py), backed for native algorithms by the [C++ factory implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/include/qdk/chemistry/algorithms/algorithm.hpp). The [registry tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_algorithms_registry.py) explicitly check default creation, named creation, keyword settings, and invalid-name errors. These are source-level checks of the intended behavior, not tests newly run for this course.

## The first real contract: molecule in, energy and state out

Using the `structure`, `charge`, and `spin_multiplicity` defined for our nitrogen example, the first handoff is:

```python
solver = create("scf_solver", "qdk", method="hf")
solver.settings().set("convergence_threshold", 1e-7)

energy, wavefunction = solver.run(
    structure,
    charge=charge,
    spin_multiplicity=spin_multiplicity,
    basis_or_guess="cc-pvdz",
)
orbitals = wavefunction.get_orbitals()
hamiltonian = create("hamiltonian_constructor", "qdk").run(orbitals)
```

The two names on the left of the SCF call **unpack a tuple**, an ordered bundle of results. The first item is a floating-point total energy in Hartree; the second is a `Wavefunction`. `orbitals` is obtained from that object, not mistaken for the second return value itself. The later Hamiltonian call returns one object, so it is not unpacked into the same pair.

The [SCF interface](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/include/qdk/chemistry/algorithms/scf.hpp) makes all four inputs explicit. Charge and multiplicity are run arguments, not universal settings. `basis_or_guess` accepts a basis-name string, an explicit `BasisSet`, or `Orbitals` containing a reusable orbital guess and its basis. An algorithm setting called `basis_set` should not be invented because an old demonstration table happens to contain it.

The distinction between arguments and settings is practical. Arguments identify the problem data for this run; settings configure how this algorithm tackles it. The same configured solver may be run again on compatible new inputs. Changing its numerical configuration after execution is a different matter.

## Settings are typed controls, not an unrestricted dictionary

`settings()` returns the algorithm-owned settings object. `get("max_iterations")` reads a value; `set` changes one existing key; `update` can apply several overrides; `to_dict()` gives a readable snapshot. Factory keyword arguments are forwarded to the same update mechanism. They are not miscellaneous constructor parameters forwarded to arbitrary third-party software.

The native base `ElectronicStructureSettings` defines `method="hf"`, `convergence_threshold=1e-7`, `max_iterations=50`, and `scf_type="auto"`. The native implementation extends that schema with settings such as `eri_method`, `enable_gdm`, and `nthreads`. These are real implementation-dependent controls. PySCF has its own extended schema, so an option accepted by one backend need not be accepted by another.

The C++ store distinguishes booleans, signed 64-bit integers, floating-point values, strings, several vector types, and nested `AlgorithmRef` values. A setting can have a numeric bound or an allowed-value list. For example, `max_iterations` must be at least one and native `eri_method` is one of `"direct"` and `"incore"`. Use `80` for an integer count and `1e-7` for a floating-point threshold. The Python binding deliberately rejects a boolean where an integer is expected, despite Python normally treating booleans as integer-like.

Importantly, the [actual `Settings::set` implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/data/settings.cpp) checks existing keys, type compatibility, and declared limits when setting values. The documentation phrase “validated at execution time” is not the whole story. Other checks require the molecule or a combination of settings and therefore happen only during execution. A valid string type cannot prove that a named basis file exists or that a requested method supports a particular spin case.

Both the [C++ algorithm base](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/include/qdk/chemistry/algorithms/algorithm.hpp) and [Python algorithm base](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/base.py) implement ordinary execution as “lock settings, then call `_run_impl`.” The leading underscore marks implementation machinery, not the method a course workflow should call directly.

Locking happens **before computation**, not only after success. Even a failing underlying run has already locked its settings. A later attempt to change them raises `SettingsAreLocked`; make a new solver to alter the configuration. Reusing a locked solver with the same settings and new appropriate arguments is not forbidden. The lock is a configuration safeguard, not a declaration that every algorithm object is thread-safe.

There are two precise architectural qualifications. The registry wraps algorithms to support optional cache and remote execution; a cache hit can return without entering the underlying `run`, so the lock behavior described above applies to actual underlying execution. Also, `Settings::lock()` sets that settings object’s lock flag; it is not a recursively enforced freeze of every separately reachable nested object. Configure the entire computation before running and do not mutate retained nested settings afterward. The intended usage discipline is stronger than a claim of universal deep immutability.

## Following the native call across the language boundary

Now trace the nitrogen call without getting lost in integral mathematics.

The Python name `ScfSolver` is re-exported from the compiled `_core._algorithms` extension. The [pybind11 SCF binding](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/pybind11/algorithms/scf.cpp) exposes `run(structure, charge, spin_multiplicity, basis_or_guess)` and converts the C++ pair of energy and wavefunction pointer into the Python tuple. **pybind11** is the bridging library that makes C++ classes callable from Python; it is not another chemical approximation.

The Python registry finds the C++ `ScfSolverFactory`. Its [registration source](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/scf.cpp) supplies a creator for `microsoft::ScfSolver`. That concrete implementation reports its name as `"qdk"`. After the base locks its settings, the [native adapter](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/scf.cpp) resolves the basis, converts structural data to the internal engine, selects the HF or KS solver, constructs an initial density, and runs the internal SCF calculation.

When that engine returns, the adapter checks convergence rather than silently labeling any final iterate a solution. It maps the internal shell ordering back to the QDK basis ordering. It creates `Orbitals` containing coefficients, orbital energies, overlap information, and the basis. It then constructs a canonical occupied `Configuration`, a `StateVectorContainer`, and finally a `Wavefunction`. The energy is taken from the engine’s total SCF energy. The KS path also packages its orbital determinant this way; packaging a KS reference does not make it an exact interacting many-electron state.

A C++ `shared_ptr` is a shared ownership handle: several objects can refer to the same data while its lifetime remains valid. It is not “copy all matrices now.” The binding uses lifetime-aware holders and returns the settings with a policy that keeps their owner alive. These engineering details let Python code hold a wavefunction or its orbitals without manually managing C++ memory. They do not excuse modifying somebody else’s scientific input.

The binding also has a **trampoline**, a bridge that lets a Python subclass implement a C++ virtual method such as `_run_impl`. That is how a Python-written SCF adapter can satisfy the native interface. The PySCF implementation inherits this bound `ScfSolver`, instead of requiring the entire chemistry pipeline to adopt PySCF-specific objects.

## Immutable scientific records, with carefully distinguished caches

The word immutable means that a data object’s defining scientific information is not normally rewritten after construction. To rotate orbitals or select a different space, an algorithm returns a new object. Keeping the old object allows comparisons and preserves the meaning of calculations that already used it.

The [Python `DataClass` base](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/data/base.py) marks an object initialized at the end of construction and rejects ordinary attribute assignment and deletion afterward. C++ data classes expose their own non-mutating interfaces. Do not generalize this into a promise that every nested Python dictionary, foreign-library object, or array in every subclass is recursively frozen.

Lazy derived information is different from defining information. A **lazy cache** computes a derived value when first requested and retains it to avoid repeating work. The native wavefunction/container API has caches for derived reduced-density information. A particularly visible example is [`Circuit.get_qir()`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/data/circuit.py): when a Q# factory is available, it can compile and internally cache QIR using an explicit internal attribute assignment, although user reassignment such as `circuit.qasm = ...` is rejected. The [circuit tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_circuit.py) check ordinary immutability.

“Self-contained” likewise needs interpretation. Some orbital properties are optional; not every wavefunction container can provide determinant coefficients or all reduced density matrices. A consumer must require the capabilities it actually uses. Common inheritance alone does not manufacture missing scientific data.

There is also a naming trap. A chemical `Configuration`, such as the occupation string `2du0`, records particle occupations in an orbital/mode convention. Algorithm configuration means settings such as a convergence threshold. They are unrelated uses of the same English word.

## Composing tools rather than hard-coding them

An `AlgorithmRef` records a nested algorithm type, implementation name, and optional settings. It describes what to construct later; it is not itself a running solver. For example:

```python
from qdk_chemistry.data import AlgorithmRef

stabilized = create("scf_solver", "qdk_stabilized", max_iterations=80)
stabilized.settings().set(
    "scf_solver",
    AlgorithmRef("scf_solver", "qdk", eri_method="incore"),
)
```

Here a composite SCF workflow will use a native SCF child with in-memory integrals. This illustrates software composition, not a recommendation that incore is best for nitrogen. `AlgorithmRef` can resolve default child settings through the registry and validate overrides. Its algorithm type is fixed; changing its implementation name re-resolves defaults and discards prior overrides. The [nested-settings tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_settings.py) cover construction, type protection, and multi-level references.

Read the composition policy too. In [`qdk_stabilized`](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/cpp/src/qdk/chemistry/algorithms/microsoft/stabilized_scf.cpp), common parent SCF settings are copied onto the newly created child after nested creation. Thus the parent’s `max_iterations` takes precedence over a competing child override of that common key. The child-only `eri_method` in this example is a cleaner illustration. This wrapper also creates a stability checker and can rerun SCF using modified orbitals; the physical meaning comes in Lesson 6.

Custom integrations register a **creator callable**, not one already-used shared solver: `register(lambda: MySolver())` means “construct a fresh instance when requested.” For native SCF integration, `MySolver` should implement the bound SCF interface, including the full four-argument contract and the energy/wavefunction pair. A new algorithm type needs its own factory. Existing names and aliases cannot be silently hijacked: registration raises `DuplicateRegistrationError` on conflicts.

## One genuinely Python-built quantum path

Classical C++ is not a mandatory intermediate language for every algorithm. Consider `create("state_prep", "dense_pure_state")`. This is a registered Python implementation of `StatePreparation`; its input is a coefficient-bearing `Wavefunction` and its output is `Circuit`.

The [common helper](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/state_preparation.py) reads coefficients and active determinants. It does not assume the coefficient list already sits at the correct bit-string indices. It derives each determinant’s index from occupation bits and scatters the coefficient into a dense array. Non-negligible imaginary coefficients are rejected by this real-amplitude path, and widths greater than 32 qubits are refused.

The [concrete implementation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/algorithms/state_preparation/dense_pure_state.py) constructs Q# parameters, including a reversed register-row map, and obtains a callable plus factory data. The [repository’s Q# operation](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/utils/qsharp/src/StatePreparation.qs) delegates dense amplitude loading to `Std.StatePreparation.PreparePureStateD`. Chemistry code therefore supplies scientific amplitudes and ordering; the external QDK/Q# machinery supplies language execution and compilation. The returned circuit declares `encoding="jordan-wigner"` and a register width. It need not already contain an OpenQASM string, despite a broad base-method docstring describing one.

This is a method with a real resource tradeoff, even before quantum execution: a dense array has \(2^q\) entries for \(q\) qubits. An illustrative four-spatial-orbital electron problem uses eight occupation qubits and 256 entries. At 32 qubits, just one float64 array occupies 32 GiB; the guard is not a promise that every accepted input is affordable. Our full 28-spatial-orbital nitrogen representation would need 56 occupation qubits and is outside this dense implementation’s limit. We are tracing an interface on a small example, not proposing this full-molecule call. Later state-preparation lessons will compare sparse and other approaches.

The [dense-preparation tests](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/tests/test_state_preparation_dense.py) compare prepared states against determinant-derived expectations, including bit ordering. That is why metadata matters: a circuit preparing Jordan–Wigner occupation bits is not automatically compatible with a Hamiltonian using parity or Bravyi–Kitaev encoding merely because both objects have the same qubit count.

## Plugins, identity, and useful failures

The [plugin API](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/plugins/base.py) exposes `QdkChemistryPlugin` with API version 1 and a `PluginRegistrar`. An external installed package can advertise an entry point in the `qdk_chemistry.plugins` group. Registration can contribute algorithms, new factories, data-class loaders, remote backends, and cache backends. The [loader](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/plugins/__init__.py) validates the plugin class/version and warns when an external plugin fails to load. Bundled optional integrations have a separate autoload path; missing optional imports are handled differently from duplicate registrations.

The [data registry](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/src/qdk_chemistry/data/registry.py) maps a stable serialized type name to a loader, rather than selecting a computational method. A plugin returning a new data class must supply that loading contract if another process is to reconstruct its results. Full serialization and remote execution come later.

Source presence is not installed availability. The [package metadata](https://github.com/microsoft/qdk-chemistry/blob/97a3ff4f/python/pyproject.toml) makes PySCF an optional integration and excludes that dependency from the `plugins` extra on Windows. An environment’s loaded registry is the relevant runtime catalog, not a promise that every file in this repository is usable on this machine. For this read-only preparation, we inspected source and metadata instead of importing the package, installing plugins, or triggering resource/stub initialization.

Algorithms also expose `hash(...)`, using the run’s algorithm identity, settings, and inputs to identify computations for caching. Native and Python hashing have implementation-specific encodings; do not assume interchangeable hashes across languages, versions, or different argument-passing forms. A short content hash is a bookkeeping identifier, not proof of chemical equivalence, absence of collisions, or bitwise reproducibility across machines. Preserve explicit geometry, units, charge, basis, method, implementation, settings, and dependency provenance too.

Failures should be read at the correct layer. `KeyError` from registry creation points to a missing type or implementation. `SettingNotFound` and `SettingTypeMismatch` identify configuration mistakes. An out-of-range value can raise `ValueError`. A native `RuntimeError` saying SCF did not converge means the computation did not satisfy its numerical stopping conditions; it is not a request to substitute a plausible energy. A missing circuit conversion dependency is a capability failure, not a bad molecule.

We can now explain why exchanging SCF engines need not rewrite downstream chemistry: each adapter returns the same scientific kind of object, with explicit conventions to check. We can also see what that architecture cannot decide for us. The next step is to make the molecule and its one-electron representation unambiguous before any solver begins work.