04Specifying the molecule and representation

Before solving anything, specify what the numbers mean

Our next task is not to find an energy. It is to make the input so unambiguous that another colleague can calculate the same physical model. We need a nuclear framework, a total electron count and spin target, and a finite set of functions in which to represent the electrons. A molecule name alone supplies none of these precisely enough.

The running case remains neutral singlet N₂ at 1.85 Å, from the nitrogen tutorial source. The separate resource-estimation notebook uses 1.270025 Å. Those are different geometries, so its energies or resource results cannot be attached to this calculation.

All code below is instructional, not a record of calculations executed during course preparation. Numerical arithmetic is explicitly distinguished from electronic-structure results.

Read the XYZ as a contract, not a picture

The tutorial constructs its structure from this text:

from qdk_chemistry.data import Structure

structure = Structure.from_xyz("""2
Stretched N2 molecule for the ground-state QPE tutorial
N    0.000000    0.000000    0.000000
N    0.000000    0.000000    1.850000
""")
charge = 0
spin_multiplicity = 1

The first line is the number of atoms. The second is a comment line, which must be present even if empty. Each subsequent atom line gives an element symbol and three Cartesian coordinates. Cartesian here means distances along perpendicular x, y, and z axes; it does not refer to the later choice between Cartesian and spherical basis functions.

QDK Chemistry interprets these XYZ coordinates in Ångstrom and converts them to Bohr internally. A Bohr is a unit of length, not an atom or an orbital. The default CODATA 2022 constant is 0.529177210544 Å per Bohr. Thus the tutorial separation is, by unit conversion alone,

Here $R$ denotes the internuclear distance. This arithmetic does not require an SCF calculation. Passing a raw coordinate value of 1.85 to a Structure constructor instead would mean 1.85 Bohr, a different, much shorter bond.

The Structure interface uses Bohr for direct construction, coordinate getters, and JSON coordinates. For example, a direct-construction equivalent uses a two-row, three-column coordinate array with the second atom at approximately 3.495993333 Bohr. The exact converted value is preferable to rounding when preserving a reference calculation.

Structure.from_xyz reads a text string. Structure.from_xyz_file reads a file; the library’s filename convention is a name such as molecule.structure.xyz. This course does not create such a file. get_coordinates() returns all coordinates, while get_atom_coordinates(1) selects the second atom: indices start at zero.

Most importantly, the XYZ comment is not a program for setting charge, spin, basis, or method. The parser skips it. Writing “singlet, cc-pVDZ” in the comment does not configure SCF. XYZ also supplies no bond-order constraints, solvent model, thermodynamic state, active space, convergence tolerance, or error budget. Those belong elsewhere in the workflow.

Validation is layered, and chemical judgment remains necessary

Structure stores coordinates, element identities, masses, and nuclear charges. Its constructors check that the arrays have consistent sizes and that a nonempty coordinate matrix has three columns. Symbol lookup normalizes ordinary capitalization and rejects an unknown element. XYZ parsing reports missing atom counts, missing comment lines, or atom lines that cannot supply the required symbol and coordinates. The structure tests exercise malformed input and isotope lookup.

Those checks do not establish that the input is a chemically sensible molecule. The concrete _is_valid() checks dimensions, not every coordinate for finiteness, positive masses, plausible distances, or unwanted coincident atoms. The XYZ reader consumes the advertised atom records; it is not a full extended-XYZ metadata interpreter. Do not treat successful construction as a geometry quality certificate.

For ordinary calculations, we should independently establish that coordinates are finite, atoms are not accidentally superimposed, the atom count and element identities match the intended species, and the length units are known. Coincident nuclei are especially problematic because the nuclear repulsion contains division by their separation.

Charge and multiplicity are not properties stored on Structure in this interface. They are supplied to ScfSolver.run. This allows the same nuclear geometry to be studied in multiple charge or spin states without pretending that the atoms moved.

Let $Z_A$ be the physical nuclear charge of atom $A$, measured as a positive integer number of proton charges, and let $Q$ be the signed molecular charge. Positive $Q$ means electrons were removed. In an all-electron model,

For N₂, each nitrogen has $Z_A=7$. With $Q=0$, there are fourteen electrons. Making the molecule a cation changes the electron count through charge; it does not turn a nitrogen nucleus from charge seven into charge six.

The requested multiplicity is $M=2S+1$, where $S$ is the target total electronic spin. Conventional highest-spin-projection bookkeeping uses

Here $N$ is the number of electrons actually represented by the electronic solver. These counts must be nonnegative integers and fit the available orbitals. Thus an even-electron system requires odd multiplicity in this bookkeeping. Neutral singlet N₂ has seven alpha and seven beta electrons. Neutral N₂ with multiplicity two would imply half-integer population counts and is inconsistent.

The native SCF engine checks consistency of the resulting populations and the requested spin difference. Its high-level adapter also has a negative-multiplicity fallback to singlet for even electron count or doublet for odd count; that is a convenience policy, not a scientific inference about the ground state. We specify multiplicity explicitly rather than rely on backend-specific inference. Equal alpha and beta populations alone still do not guarantee a spin-pure correlated or unrestricted state.

Mass, isotope, and nuclear charge are different information

The default element masses come from the CIAAW-based element data. A standard atomic weight represents a reference isotopic composition; it is not necessarily the mass of one particular isotope. An isotope has the same number of protons but a different neutron count.

The string interface recognizes isotope forms such as "N15" and "C13", as well as "D" and "T" for hydrogen isotopes. For example, Structure.get_default_atomic_mass("N15") uses the stored nitrogen-15 mass, while its element identity remains nitrogen. One may also explicitly supply the masses array in atomic mass units. The stored isotope records and tests distinguish these values from standard atomic weights.

At the fixed-nucleus, nonrelativistic electronic level we are using, changing isotope mass while holding nuclear charges and positions fixed does not change the electronic Coulomb Hamiltonian. Mass matters when we later discuss nuclear vibration, mass-weighted Hessians, and dynamics. It should still be preserved as input provenance.

The data class permits custom floating-point nuclear charges, and data tests demonstrate that storage capability. The native molecular adapter is narrower: its charge conversion requires finite, nonnegative, integral values within a tolerance of $10^{-12}$. A data container accepting fractional charges does not imply the SCF backend implements a fractional-nucleus model. The PySCF structural conversion is based on element labels, so it should not be assumed to preserve arbitrary custom nuclear-charge fields either.

XYZ is a deliberately lossy interchange representation here. to_xyz() emits ordinary element symbols and six decimal places in Ångstrom. It does not export a separate custom-mass or custom-nuclear-charge array; an isotope stored as nitrogen plus a special mass will be printed as N. Consequently, an XYZ round trip should not be used as proof that all metadata or full coordinate precision survived. Structured serialization is the later provenance solution.

A basis set is a vocabulary for orbital shapes

Knowing where the nuclei are still leaves infinitely many possible functions of electron position. A basis set chooses a finite vocabulary of functions centered on those nuclei. The molecular orbital labeled $p$ is expanded as

$\mathbf r$ is an electron position, $\chi_\mu$ is basis function number $\mu$, $K$ is the basis-function count, and $C_{\mu p}$ is a molecular-orbital coefficient. These coefficients are optimized by SCF. They are not the determinant amplitudes from Lesson 2, nor the contraction coefficients used to define a basis function itself.

Most of this repository’s molecular machinery uses Gaussian-type orbitals. A Cartesian primitive centered at $\mathbf R_A=(X_A,Y_A,Z_A)$ has the schematic form

$a,b,c$ are nonnegative integers specifying the angular polynomial, $\alpha>0$ is an exponent, and $\mathcal N$ is a normalization factor. A large exponent makes a tight function near the nucleus; a small exponent makes a diffuse function reaching farther out. The total angular momentum label is $\ell=a+b+c$: s, p, d, f correspond to $\ell=0,1,2,3$.

Gaussians are computationally useful because products and integrals involving them have efficient algebraic recurrences. One Gaussian is usually a poor imitation of the radial shape of an atomic orbital, especially near the nuclear cusp and far from the nucleus. A contraction combines several primitives with fixed weights:

The $d_{k\mu}$ values define the basis, while the separate $C_{\mu p}$ values define the molecular orbitals in that basis. More primitives inside one contraction do not automatically mean more independently adjustable molecular-orbital coefficients.

The actual hydrogen STO-3G record in the distributed archive has exponents approximately 3.425250914, 0.6239137298, and 0.1688554040, with raw contraction coefficients 0.1543289673, 0.5353281423, and 0.4446345422. These are source data, not a fitted result from this session. The three primitives form one s basis function on hydrogen, not three independent orbitals. “STO-3G” refers to a three-Gaussian approximation to a Slater-type orbital shape.

Shells make the representation concrete

A Shell groups functions sharing an atom, angular momentum, and primitive/contraction data. A p shell has three angular orientations. A basis can contain several separate p shells with different radial flexibility; these are not three electrons placed in a shell.

The actual Shell definition stores atom_index, orbital_type, exponents, coefficients, and optional ECP radial powers. It checks matching exponent/coefficient lengths; ECP radial-power arrays must have compatible length too. The public OrbitalType labels extend through i, or $\ell=6$, with UL reserved for an ECP local-potential label. Such an enumeration does not promise that every compiled integral engine supports every calculation for every listed element and angular momentum.

The stored contraction coefficients are raw, unnormalized coefficients. Integral/backend conversion handles the required normalization conventions. Normalizing the list of coefficients as if it were a probability vector would be wrong: overlapping primitive functions must be normalized through their overlap integrals, not by dividing by the Euclidean length of the coefficient list. Individual contraction coefficients may legitimately be negative.

For example, an instructional explicit construction using those hydrogen parameters is:

import numpy as np
from qdk_chemistry.data import AOType, BasisSet, OrbitalType, Shell, Structure

hydrogen = Structure(np.array([[0.0, 0.0, 0.0]]), ["H"])
h_shell = Shell(
    0, OrbitalType.S,
    np.array([3.425250914, 0.6239137298, 0.1688554040]),
    np.array([0.1543289673, 0.5353281423, 0.4446345422]),
)
explicit_basis = BasisSet(
    BasisSet.custom_name, [h_shell], hydrogen, AOType.Spherical,
)

This demonstrates the representation, not a nitrogen model and not an executed hydrogen calculation. Shell is a construction record with writable fields in its binding; assembling it is different from modifying a completed scientific BasisSet. An actual custom basis requires physically appropriate positive exponents, supported angular momenta, matching geometry, and tested normalization/order conventions. The container’s structural validation is not an automatic basis-design service.

There are two angular-function choices. Spherical shells use $2\ell+1$ functions; Cartesian shells use $(\ell+1)(\ell+2)/2$. Both give one s and three p functions. A d shell has five spherical versus six Cartesian functions, and an f shell has seven versus ten. This changes the represented function space and counts, not merely a label printed on the result.

AOType.Spherical is the named-basis default. Although the data class can represent AOType.Cartesian, the current high-level molecular adapters should not be presented as supporting a complete Cartesian workflow. The native conversion serializes its internal basis as pure=true, and its reverse conversion explicitly rejects a non-pure basis. The PySCF conversion also retains Cartesian TODOs and constructs spherical output. We use spherical functions consistently instead of relying on unsupported conversion behavior.

What the familiar basis names buy us

A minimal basis supplies roughly the occupied atomic-orbital vocabulary. A split-valence basis gives valence orbitals additional radial flexibility. Double-zeta and triple-zeta indicate increasing radial flexibility in the family’s construction, not an exact universal multiplier of every function count. Polarization functions add angular flexibility, allowing the electron distribution to bend away from isolated-atom shapes. Diffuse functions add long-range flexibility, particularly relevant to anions, extended density tails, and some excited states. None of these additions changes the number of physical electrons.

Representative resources actually present include sto-3g, 6-31g, 6-31g*, cc-pvdz, cc-pvtz, aug-cc-pvdz, def2-svp, def2-tzvp, ano-rcc, and lanl2dz. The star in the Pople example denotes added polarization; aug- denotes diffuse augmentation in the correlation-consistent family. The cc-pVXZ family is designed for systematic valence-correlation convergence; selecting it does not itself turn on a correlated solver or freeze core electrons. The def2 family offers a different optimized hierarchy and includes ECP treatment for some heavier elements, not every element.

The basis/functionals catalog is helpful, but actual availability is determined by resources and supported element combinations. BasisSet.get_supported_basis_set_names() inspects the resource archives. BasisSet.get_supported_elements_for_basis_set("cc-pvdz") consults basis_summary.json. Listing an archive is not a proof that every element is covered, every high-angular-momentum integral is compiled, or an auxiliary fitting basis is suitable as the orbital basis. Archives with names such as -jkfit or -rifit serve auxiliary representations for integral approximations; they are not interchangeable substitutes for an orbital basis.

There is concrete external provenance here. The repository credits Basis Set Exchange as the source of the distributed basis and ECP definitions. BSE and the original basis authors provide exponents, contraction weights, angular-momentum records, and references. QDK Chemistry packages and parses those definitions; it does not invent them by fitting the molecule during SCF. The archived nitrogen cc-pVDZ and cc-pVTZ records retain the original dunning1989a reference key.

Ordinary named loading reads installed resources and extracts packaged archives; it is not a live scientific search of the web. For this read-only lesson preparation, the relevant archive members were inspected directly in memory, without calling the package loader or extracting files.

Count the nitrogen representation from its shells

The stored nitrogen cc-pVDZ record contains three s shells, two p shells, and one d shell. With spherical functions, one atom contributes

basis functions. Two atoms contribute 28. This explains the tutorial’s recorded 28-spatial-orbital result; it is not an unexplained magic number.

The archived nitrogen cc-pVTZ record contains four s, three p, two d, and one f shell, giving

functions per atom, or 60 for N₂. The 60 count is derived from the resource definition, not from a newly run calculation. The two archives are cc-pVDZ and cc-pVTZ.

Both models still have fourteen electrons and seven occupied spatial orbitals in a closed-shell RHF reference. The other orbitals are unoccupied, or virtual, in that reference. Virtual does not mean irrelevant: they permit orbital relaxation and provide destinations for excitations in correlated models. Before any active-space reduction, 28 spatial orbitals correspond to 56 electron spin modes, and 60 correspond to 120. That is a representation count, not a complete logical or physical resource estimate.

Overlap explains why more functions can become numerically difficult

Basis functions centered on different atoms generally overlap. Their overlap matrix is

The star denotes complex conjugation; for the real functions used in this native SCF path it has no effect. A diagonal entry measures a function’s norm, and an off-diagonal entry measures overlap between two different functions. An orthonormal set would have ones on the diagonal and zeros elsewhere. An atomic Gaussian basis usually does not.

If two functions are nearly duplicates, one combination almost cancels to zero. This is near linear dependence. Numerically, the overlap matrix then has a tiny eigenvalue, and dividing by its square root magnifies roundoff. Diffuse augmentation and geometries with closely overlapping functions can make this worse.

The native orthogonalization routine diagonalizes the overlap, retains sufficiently large eigenvalues, and builds an orthogonalizing transformation. The low-level default threshold is $10^{-6}$. If that low-level setting is negative, the routine instead uses the largest overlap eigenvalue divided by $10^9$. That low-level configuration is not a lindep_threshold option exposed by the high-level qdk SCF settings at this revision; do not invent such a factory keyword. The number of retained molecular orbitals can therefore be less than the original basis-function count.

This is why we inspect both basis.get_num_atomic_orbitals() and wavefunction.get_orbitals().get_num_molecular_orbitals(). Their common equality in small well-conditioned examples is not a universal invariant.

Choose a name, a constructed basis, or an orbital guess deliberately

The normal named-basis route is:

basis = BasisSet.from_basis_name("cc-pvdz", structure)
number_of_functions = basis.get_num_atomic_orbitals()
first_atom_shells = basis.get_shells_for_atom(0)
energy, reference = create("scf_solver", "qdk").run(
    structure, 0, 1, basis,
)

The same native calculation can instead receive "cc-pvdz" directly and let the solver load it. Using an explicit object makes the actual shell/ECP representation inspectable before solving. BasisSet.from_element_map assigns a basis for each element present; from_index_map assigns one for each atom index. Missing element or atom assignments are errors, not an invitation to guess a default. Mixed choices produce a basis named custom_basis_set; that label alone cannot reconstruct its contents.

For N₂, giving the two nitrogen atoms different basis quality would break the symmetric representation even when their physical identities are equivalent. Per-atom customization is a capability to use intentionally, not a shortcut to making one side of a bond cheaper without consequence. The basis tests demonstrate named, element-map, and index-map construction and their native SCF handoff.

An Orbitals input instead requests reuse of an initial orbital guess. The native adapter takes the basis from those orbitals, builds a starting density from occupied coefficient columns, and then optimizes again. It does not mean “return the old energy,” and it is not an automatic projection from cc-pVDZ to cc-pVTZ. Use a consistent structure, basis, orbital dimension, and spin convention. A geometry scan needs a correctly rebuilt or explicitly transferred guess, not a stale basis object whose centers still describe a different geometry.

Backend boundaries matter. At this revision PyscfScfSolver._build_scf accepts a name or orbital guess but explicitly raises NotImplementedError for any directly supplied BasisSet object, including one loaded from a standard basis name. Separate basis_to_pyscf_mol utilities can translate explicit shells; their existence does not remove that solver-level limitation. The PySCF solver also rebuilds its mean-field molecule using the basis name obtained from an orbital guess, so a custom guess basis must not be assumed to work as a named database basis. Spherical conventions, shell ordering, normalization, and the exact represented basis must be checked when transferring data.

Effective core potentials change the Hamiltonian, not the element

For heavier atoms, explicitly representing tightly bound core electrons can be expensive and may require relativistic modeling. An effective core potential, or ECP, replaces selected core electrons and their interaction with the remaining electrons by an effective operator. It is more than subtracting electrons: its angular channels and radial functions alter the one-electron Hamiltonian.

Let $n_A^{\mathrm{core}}$ be the number of electrons represented by an ECP on atom $A$. The implemented effective nuclear charge and explicit-electron count are

Structure.get_nuclear_charges() retains the physical stored charges. BasisSet.get_ecp_electrons() supplies per-atom replaced-core counts, and BasisSet.get_effective_nuclear_charges() performs the subtraction. Do not overwrite the structure with already-reduced charges and then ask the basis to subtract the ECP again.

A useful source-grounded example, separate from nitrogen, is neutral AgH with def2-SVP. Silver has 47 physical electrons when neutral and hydrogen one, so the physical count is 48. The packaged silver ECP replaces 28 core electrons; hydrogen has none. The model therefore explicitly represents $47-28+1=20$ electrons. Calling all nineteen retained silver electrons “the chemically active valence” would be misleading: an ECP partition and a later chemically chosen active space are different decisions.

The core–core repulsion convention changes with that effective Hamiltonian. For a separation $R$ in Bohr, the structure’s physical nuclear expression is $47/R$ Hartree for AgH, while basis.calculate_effective_nuclear_repulsion_energy() uses $19/R$. The native AgH test checks physical charges, effective charges, the Hamiltonian constant, and agreement between the mean-field energy and its Hamiltonian expectation. These are repository test assertions, not newly generated AgH results. Comparing all-electron and ECP absolute energies without reconciling their Hamiltonians and energy references is not meaningful.

The basis parser loads included ECP data automatically and currently accepts scalar_ecp records with one angular momentum per potential entry. It does not provide a general spin–orbit ECP treatment. ECP electron-count arrays must match the atom count and cannot exceed the corresponding nuclear charge. Native internal conversion preserves per-atom treatment: a test assigns an ECP to one Ag atom and an all-electron basis to another and verifies they remain different. One universal core count per element would lose that information.

In the native engine, Libint receives Gaussian basis/geometry data and evaluates ordinary overlap, kinetic, attraction, and electron-repulsion integrals. Libecpint receives ECP radial/angular data plus Gaussian shells and evaluates the extra ECP matrix elements. The one-electron engine uses both. Neither library decides which core electrons are safe to replace for the scientific question; that choice comes from the selected model and its validated ECP definition.

An ECP is also not the same as freezing occupied orbitals after an all-electron SCF calculation. Frozen-core correlation retains an all-electron reference and removes selected excitations from a later many-electron treatment. ECP selection changes the starting Hamiltonian and explicit-electron model itself.

Make a controlled comparison, then hand off the right model

Our nitrogen tutorial intentionally changes only the basis from cc-pVDZ to cc-pVTZ. Geometry, charge, spin target, HF approximation, and native solver remain fixed. Its guarded reference energies are approximately −108.418633697214 and −108.445215657498 Hartree. These are values recorded in the example source, not calculations performed here. The difference is a basis sensitivity within that HF comparison, not a measured error against exact chemistry.

More basis flexibility can lower the variational minimum when one space contains another and the same variational problem is solved to the corresponding minimum. Arbitrary named bases need not be nested, SCF can converge to different stationary solutions, and DFT has different functional error. Therefore “a larger basis must always give a lower and more correct energy” is not a sound general rule. Stretched N₂ can still need multiple determinants even after the one-electron representation is improved.

We now have a complete input handoff: the geometry and units, element/isotope information, physical and explicit electron counts, requested spin, actual basis shells and angular convention, any ECP definition, and the chosen algorithm implementation/settings. We have not yet solved the electronic problem. Lesson 5 takes exactly these inputs and explains how integrals, density matrices, and the self-consistent field loop produce the first electronic solution.