Self-consistent field (SCF) solver

The ScfSolver algorithm in QDK/Chemistry performs Self-Consistent Field (SCF) calculations to optimize molecular orbitals for a given molecular structure. Following QDK/Chemistry’s algorithm design principles, it takes a Structure instance, total molecular charge and multiplicity, and a desired basis representation as input ,and produces an Orbitals instance and the associated energy as output. Its primary purpose is to find the best single-particle orbitals within a mean-field approximation. For Hartree-Fock (HF) theory, it yields the mean field energy, which misses electron correlation and typically requires post-HF methods for accurate energetics. For Density Functional Theory (DFT), some correlation effects are included through the exchange-correlation functional.

Overview

SCF theory encompasses both HF and DFT methods in quantum chemistry. Both methods rely on a single Slater determinant representation of the many-electron wavefunction, using molecular orbitals that are optimized to minimize the electronic energy. This single-determinant approach is a key simplification that makes these methods computationally efficient but limits their ability to capture certain correlation effects. The SCF procedure iteratively refines these orbitals until self-consistency is achieved.

The orbitals from SCF calculations typically serve as input for these post-SCF methods which capture correlation effects. SCF methods thus serve as the foundation for more advanced quantum and classical electronic structure calculations and provide essential insights into molecular properties, reactivity, and spectroscopic characteristics.

Running an SCF calculation

This section demonstrates how to setup, configure, and run a SCF calculation. The run method returns two values: a scalar representing the converged SCF energy and a Wavefunction object containing the optimized molecular orbitals.

Input requirements

The ScfSolver requires several inputs to perform a calculation:

Structure

A Structure instance defining the molecular geometry (atomic positions and element types).

Charge

The total molecular charge (integer). A neutral molecule has charge 0.

Spin multiplicity

The spin multiplicity of the system, defined as \(2S + 1\) where \(S\) is the total spin. Common values are 1 (singlet), 2 (doublet), 3 (triplet), etc.

Basis set or initial guess

This required input specifies the atomic orbital basis for the calculation and can be provided in several forms:

String

A standard basis set name (e.g., "sto-3g", "def2-svp", "cc-pvdz"). See the basis set documentation for available options.

BasisSet object

A BasisSet instance for custom basis sets. See the BasisSet documentation for details.

Orbitals object

A Orbitals instance which provides both the basis set and an initial orbital guess for the SCF optimization.

Creating a SCF solver

C++ API

#include <iostream>
#include <qdk/chemistry.hpp>
#include <string>
using namespace qdk::chemistry::algorithms;
using namespace qdk::chemistry::data;

// Create default ScfSolver instance
auto scf_solver = ScfSolverFactory::create();

Python API

import numpy as np
from pathlib import Path
from qdk_chemistry.algorithms import create
from qdk_chemistry.data import Structure

# Create the default ScfSolver instance
scf_solver = create("scf_solver")

Configuring settings

Settings can be modified using the settings() object. See Available settings below for a complete list of options.

C++ API

  // Set the method
  // Note the following line is optional, since hf is the default method
  scf_solver->settings().set("method", "hf");

Python API

# Configure the SCF solver using the settings interface
# Note that the following line is optional, since hf is the default method
scf_solver.settings().set("method", "hf")

Running the calculation

C++ API

  // Load structure from XYZ file
  auto structure = Structure::from_xyz_file("../data/h2.structure.xyz");

  // Run the SCF calculation
  auto [E_scf, wfn] = scf_solver->run(structure, 0, 1, "def2-tzvpp");
  auto scf_orbitals = wfn->get_orbitals();
  std::cout << "SCF Energy: " << E_scf << " Hartree" << std::endl;

Python API

# Load structure from XYZ file
structure = Structure.from_xyz_file(Path(__file__).parent / "../data/h2.structure.xyz")

# Run scf
E_scf, wfn = scf_solver.run(
    structure, charge=0, spin_multiplicity=1, basis_or_guess="def2-tzvpp"
)
scf_orbitals = wfn.get_orbitals()

print(f"SCF Energy: {E_scf:.10f} Hartree")

Available settings

The ScfSolver accepts a range of settings to control its behavior. All implementations share a common base set of settings from ElectronicStructureSettings:

Setting	Type	Default	Description
method	string	"hf"	The method to use: "hf" for Hartree-Fock, or a DFT functional name (e.g., "b3lyp", "pbe")
convergence_threshold	float	1e-7	Convergence tolerance for orbital gradient norm
max_iterations	int	50	Maximum number of SCF iterations (must be ≥ 1)

See Settings for a more general treatment of settings in QDK/Chemistry.

Available implementations

QDK/Chemistry’s ScfSolver provides a unified interface to SCF calculations across various quantum chemistry packages. You can discover available implementations programmatically:

C++ API

  auto names = ScfSolverFactory::available();
  for (const auto& name : names) {
    std::cout << name << std::endl;
  }

Python API

from qdk_chemistry.algorithms import registry  # noqa: E402

print(registry.available("scf_solver"))  # ['pyscf', 'qdk']

QDK (Native)

Factory name: "qdk" (default)

The native QDK/Chemistry implementation provides high-performance SCF calculations using the built-in quantum chemistry engine.

Capabilities

• Restricted Hartree-Fock (RHF) and Unrestricted Hartree-Fock (UHF)

• Restricted Kohn-Sham (RKS) and Unrestricted Kohn-Sham (UKS) DFT

• Extensive library of basis sets including Pople, Dunning, and Karlsruhe families

• Full range of exchange-correlation functionals for DFT - Optimization algorithms including the direct inversion in the iterative subspace (DIIS) method [Pul82], and the geometric direct minimization (GDM) method [VanVoorhisHG02]

SCF Convergence Algorithms in QDK

Achieving stable SCF convergence is a non-trivial problem in computational chemistry. QDK/Chemistry implements two complementary algorithms that can be used independently or in combination.

Direct Inversion in the Iterative Subspace (DIIS)

DIIS is an extrapolation technique that accelerates SCF convergence by constructing an optimal linear combination of previous Fock matrices [Pul82]. DIIS is highly effective for well-behaved systems, often achieving convergence in low number of iterations. However, it can fail for challenging cases such as open-shell systems or molecules with near-degenerate orbitals, where the error surface is highly nonlinear.

Geometric Direct Minimization (GDM)

When DIIS encounters difficulties, the GDM algorithm provides a robust alternative [VanVoorhisHG02]. Rather than extrapolating Fock matrices, GDM directly minimizes the energy with respect to orbital rotation parameters using a quasi-Newton optimization approach.

The key insight of GDM is to parameterize orbital changes through unitary rotations, which converts the constrained optimization problem of determining the energy-minimizing set of orthonormal orbitals into an unconstrained optimization over exponentials [Hig05] of anti-Hermitian matrices. This allows the use of standard nonlinear optimization techniques while preserving orbital orthonormality.

The GDM algorithm then proceeds via a slightly modified [VanVoorhisHG02] BFGS optimization [LN89] which smoothly converges to a nearby energy minimum. If provided a guess close to the true minimum, GDM can converge in a similar number of iterations as DIIS, but it is more robust for difficult cases. However, if initialized further from the minimum, GDM may converge to local minima, which may require additional strategies (e.g. Stability analysis) to ensure the global minimum is found. This may be overcome in many cases by combining GDM with DIIS in a hybrid approach.

Hybrid DIIS-GDM Strategy

By default, the native QDK implementation uses DIIS alone (enable_gdm=False). When enabled, the hybrid strategy (enable_gdm=True) provides enhanced robustness:

1. Start with DIIS for rapid initial convergence

2. Monitor energy changes; if the energy change exceeds energy_thresh_diis_switch (default: \(10^{-3}\) Ha), switch to GDM

3. Once switched, continue with GDM until convergence

This hybrid approach combines the speed of DIIS for typical systems with the robustness of GDM for challenging cases.

Settings

Setting	Type	Default	Description
method	string	"hf"	Method: "hf" for Hartree-Fock, or a DFT functional name
convergence_threshold	float	1e-7	Convergence tolerance for orbital gradient norm
max_iterations	int	50	Maximum number of SCF iterations
max_scf_steps	int	100	Maximum number of overall SCF steps
enable_gdm	bool	False	Enable geometric direct minimization (GDM) algorithm
gdm_max_diis_iteration	int	50	Maximum DIIS iterations in GDM
gdm_bfgs_history_size_limit	int	50	BFGS history size limit for GDM
energy_thresh_diis_switch	float	0.001	Energy threshold for DIIS switch
level_shift	float	-1.0	Level shift parameter (negative = auto)
eri_threshold	float	-1.0	Electron repulsion integral threshold (negative = auto)
eri_use_atomics	bool	False	Use atomic operations for ERI computation
fock_reset_steps	int	1073741824	Number of steps between Fock matrix resets

PySCF

Factory name: "pyscf"

The PySCF plugin provides access to the comprehensive PySCF quantum chemistry package.

Capabilities

• Full HF support: RHF, UHF, ROHF

• Full DFT support: RKS, UKS, ROKS with extensive functional library

• Automatic spin-restricted/unrestricted selection based on multiplicity

Settings

Setting	Type	Default	Description
method	string	"hf"	Method: "hf" for Hartree-Fock, or a DFT functional name
convergence_threshold	float	1e-7	Convergence tolerance for orbital gradient norm
max_iterations	int	50	Maximum number of SCF iterations
scf_type	string	"auto"	Type of SCF calculation: "auto": Automatically detect based on spin "restricted": Force restricted calculation "unrestricted": Force unrestricted calculation

Example

from qdk_chemistry.algorithms import create  # noqa: E402
from qdk_chemistry.data import Structure  # noqa: E402

# Create and configure the PySCF solver
solver = create("scf_solver", "pyscf")
solver.settings().set("method", "b3lyp")
solver.settings().set("scf_type", "restricted")

# Run with basis set specified as input parameter
water_coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.76, 0.59], [0.0, -0.76, 0.59]])
water = Structure(water_coords, symbols=["O", "H", "H"])
energy, wfn = solver.run(water, charge=0, spin_multiplicity=1, basis_or_guess="cc-pvdz")

For more details on how to extend QDK/Chemistry with additional implementations, see the plugin system documentation.

Related classes

• Structure: Input molecular structure

• Orbitals: Output optimized molecular orbitals

Further reading

• The above examples can be downloaded as a complete Python script or C++ source file.

• Settings: Configuration settings for algorithms

• Factory Pattern: Understanding algorithm creation

• Available functionals and basis sets: Exchange-correlation functionals for DFT calculations