03How 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 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:

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 accepts Orbitals and returns Hamiltonian. The unitary-builder interface accepts QubitOperator and returns UnitaryRepresentation. The phase-estimation interface 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:

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, backed for native algorithms by the C++ factory implementation. The registry tests 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:

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 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 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 and Python algorithm base 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 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 supplies a creator for microsoft::ScfSolver. That concrete implementation reports its name as "qdk". After the base locks its settings, the native adapter 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 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(): 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 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:

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 cover construction, type protection, and multi-level references.

Read the composition policy too. In qdk_stabilized, 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 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 constructs Q# parameters, including a reversed register-row map, and obtains a callable plus factory data. The repository’s Q# operation 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 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 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 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 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 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.