Stability analysis

The StabilityChecker algorithm in QDK/Chemistry performs wavefunction stability analysis to verify that a Self-Consistent Field (SCF) solution corresponds to a true energy minimum rather than a saddle point. Following QDK/Chemistry’s algorithm design principles, it takes a Wavefunction instance as input and produces a StabilityResult object containing stability information as output. Its primary purpose is to identify directions in orbital space where the energy can be further lowered, which is crucial for ensuring the reliability of quantum chemistry calculations.

Unstable SCF solutions can occur in challenging cases such as stretched molecular geometries, open-shell systems, and molecules with near-degenerate orbitals. This iterative stability analysis of checking stability, rotating orbitals, and re-running SCF continues until a stable minimum is reached.

Overview

Stability analysis [SM91] examines the second-order response of the wavefunction energy to orbital rotations. Mathematically, this involves computing the eigenvalues of the electronic Hessian matrix. A stable wavefunction should have all positive (or non-negative) eigenvalues, indicating that the energy cannot be lowered by any orbital rotation. Negative eigenvalues signal instabilities and identify directions in which the energy can decrease.

There are two types of stability that can be checked:

Internal stability

Examines whether the wavefunction is stable with respect to rotations within the same wavefunction type. For example, restricted Hartree-Fock (RHF) internal stability checks if the energy can be lowered while maintaining the restricted formalism.

External stability

Examines whether the wavefunction is stable with respect to transitions to a different wavefunction type. The most common case is RHF external stability, which checks for RHF to unrestricted Hartree-Fock (UHF) instabilities. This is particularly important for systems that may benefit from spin symmetry breaking, such as stretched bonds or transition metal complexes.

The typical workflow for handling instabilities involves an iterative procedure:

1. Run an SCF calculation to obtain initial orbitals

2. Check stability of the resulting wavefunction

3. If unstable, extract the eigenvector corresponding to the most negative eigenvalue

4. Rotate orbitals along this eigenvector direction

5. Use rotated orbitals as initial guess for a new SCF calculation

6. Repeat until a stable solution is found

When an external instability is detected in a restricted calculation, the procedure typically switches to an unrestricted calculation type, as the instability indicates that spin symmetry breaking would lower the energy.

Running a stability check

This section demonstrates how to setup, configure, and run a stability check. The run method returns two values: a boolean indicating overall stability status and a StabilityResult object containing detailed analysis results including eigenvalues and eigenvectors.

Input requirements

The StabilityChecker requires a converged wavefunction from an SCF calculation:

Wavefunction

A Wavefunction instance containing orbital information. This is typically obtained as output from a ScfSolver calculation.

Creating a stability checker

C++ API

#include <iomanip>
#include <iostream>
#include <qdk/chemistry.hpp>
#include <qdk/chemistry/constants.hpp>
#include <qdk/chemistry/utils/orbital_rotation.hpp>
#include <string>

using namespace qdk::chemistry::algorithms;
using namespace qdk::chemistry::data;
using namespace qdk::chemistry::utils;

// Create the default StabilityChecker instance
auto stability_checker = StabilityCheckerFactory::create();

Python API

import numpy as np
from qdk_chemistry.algorithms import create
from qdk_chemistry.data import Structure
from qdk_chemistry.utils import rotate_orbitals
from qdk_chemistry.constants import ANGSTROM_TO_BOHR

# Create the default StabilityChecker instance
stability_checker = create("stability_checker", "qdk")

Configuring settings

Settings can be modified using the settings() object. See Available settings below for a complete list of options. Key settings include whether to perform internal and external stability checks, and the tolerance for determining stability.

C++ API

  // Configure stability checker settings
  stability_checker->settings().set("internal", true);
  stability_checker->settings().set(
      "external", true);  // Will be adjusted based on calculation type
  stability_checker->settings().set("stability_tolerance", -1e-4);
  stability_checker->settings().set("davidson_tolerance", 1e-4);
  stability_checker->settings().set("max_subspace", static_cast<int64_t>(30));
  stability_checker->settings().set("method", "hf");

Python API

# Configure stability checker settings
stability_checker.settings().set("internal", True)
# Will be adjusted based on calculation type
stability_checker.settings().set("external", True)
stability_checker.settings().set("stability_tolerance", -1e-4)
stability_checker.settings().set("davidson_tolerance", 1e-4)
# Maximum subspace size for Davidson solver
stability_checker.settings().set("max_subspace", 30)
stability_checker.settings().set("method", "hf")

Running the stability check and iterative workflow

The example below shows a complete iterative workflow for achieving a stable wavefunction. This includes running an initial SCF calculation, checking stability, and iteratively rotating orbitals and re-running SCF until convergence to a stable solution.

The workflow handles both internal instabilities (requiring orbital rotation within the same calculation type) and external instabilities (requiring a switch from restricted to unrestricted calculation). Note that the Davidson eigenvalue solver used in stability analysis may occasionally fail to converge; when this occurs, increasing max_subspace or adjusting davidson_tolerance can help.

C++ API

  // Create N2 molecule at stretched geometry (1.4 Angstrom)
  std::vector<Eigen::Vector3d> coords = {{0.0, 0.0, 0.0}, {1.4, 0.0, 0.0}};
  std::vector<std::string> symbols = {"N", "N"};
  for (auto& coord : coords) {
    coord *= qdk::chemistry::constants::angstrom_to_bohr;
  }
  auto n2 = std::make_shared<Structure>(coords, symbols);

  // Create and configure SCF solver with auto scf_type
  auto scf_solver = ScfSolverFactory::create();
  scf_solver->settings().set("scf_type", "auto");
  scf_solver->settings().set("method", "hf");

  // Run initial SCF calculation
  int spin_multiplicity = 1;
  auto [energy, wavefunction] =
      scf_solver->run(n2, 0, spin_multiplicity, "def2-svp");
  std::cout << "Initial SCF Energy: " << std::fixed << std::setprecision(10)
            << energy << " Hartree" << std::endl;

  // Determine if calculation is restricted and configure stability checker
  // accordingly
  bool is_restricted =
      wavefunction->get_orbitals()->is_restricted() && spin_multiplicity == 1;
  if (is_restricted) {
    stability_checker->settings().set("external", true);
  } else {
    stability_checker->settings().set("external", false);
  }

  // Iterative workflow: Check stability and rotate orbitals until convergence
  const int max_iterations = 5;
  int iteration = 0;
  bool is_stable = false;
  std::shared_ptr<StabilityResult> result;

  std::cout << "\n=== Starting Iterative Stability Workflow ===\n" << std::endl;

  while (iteration < max_iterations) {
    iteration++;
    std::cout << "--- Iteration " << iteration << " ---" << std::endl;
    std::cout << "Current Energy: " << std::fixed << std::setprecision(10)
              << energy << " Hartree" << std::endl;

    // Check stability - handle potential Davidson convergence failure
    try {
      auto [stable, stability_result] = stability_checker->run(wavefunction);
      is_stable = stable;
      result = stability_result;
    } catch (const std::runtime_error& e) {
      std::string error_msg = e.what();
      if (error_msg.find("Davidson Did Not Converge!") != std::string::npos) {
        std::cout
            << "Try increasing max_subspace or adjusting davidson_tolerance"
            << std::endl;
        throw std::runtime_error("Davidson Did Not Converge!");
      } else {
        throw std::runtime_error("Stability check failed: " + error_msg);
      }
    }

    if (is_stable) {
      std::cout << "\nConverged to stable wavefunction!" << std::endl;
      break;
    }

    // Determine rotation type based on which instability is present
    bool do_external = false;
    double smallest_eigenvalue;
    Eigen::VectorXd rotation_vector;

    if (!result->is_internal_stable()) {
      auto [eigenvalue, eigenvector] =
          result->get_smallest_internal_eigenvalue_and_vector();
      smallest_eigenvalue = eigenvalue;
      rotation_vector = eigenvector;
      std::cout << "Internal instability detected. Smallest eigenvalue: "
                << std::fixed << std::setprecision(6) << smallest_eigenvalue
                << std::endl;
    } else if (!result->is_external_stable() && result->has_external_result()) {
      auto [eigenvalue, eigenvector] =
          result->get_smallest_external_eigenvalue_and_vector();
      smallest_eigenvalue = eigenvalue;
      rotation_vector = eigenvector;
      do_external = true;
      std::cout << "External instability detected. Smallest eigenvalue: "
                << std::fixed << std::setprecision(6) << smallest_eigenvalue
                << std::endl;
    } else {
      throw std::runtime_error(
          "Unexpected state: neither internal nor external instability "
          "detected");
    }

    // Rotate orbitals along the instability direction
    auto [num_alpha, num_beta] = wavefunction->get_total_num_electrons();
    auto orbitals = wavefunction->get_orbitals();
    auto rotated_orbitals = rotate_orbitals(orbitals, rotation_vector,
                                            num_alpha, num_beta, do_external);

    // If external instability detected, switch to unrestricted calculation
    if (do_external) {
      std::cout
          << "Switching to unrestricted calculation due to external instability"
          << std::endl;

      // Get current settings values to preserve them
      std::string method = scf_solver->settings().get<std::string>("method");
      std::string stab_method =
          stability_checker->settings().get<std::string>("method");
      double davidson_tol =
          stability_checker->settings().get<double>("davidson_tolerance");
      double stability_tol =
          stability_checker->settings().get<double>("stability_tolerance");
      int64_t max_sub =
          stability_checker->settings().get<int64_t>("max_subspace");
      bool internal = stability_checker->settings().get<bool>("internal");

      // Create new solver instances
      scf_solver = ScfSolverFactory::create();
      stability_checker = StabilityCheckerFactory::create();

      // Reconfigure with updated settings
      scf_solver->settings().set("scf_type", "unrestricted");
      scf_solver->settings().set("method", method);

      stability_checker->settings().set("internal", internal);
      stability_checker->settings().set("external", false);
      stability_checker->settings().set("method", stab_method);
      stability_checker->settings().set("davidson_tolerance", davidson_tol);
      stability_checker->settings().set("stability_tolerance", stability_tol);
      stability_checker->settings().set("max_subspace", max_sub);
    }

    // Re-run SCF with rotated orbitals as initial guess
    auto [new_energy, new_wavefunction] =
        scf_solver->run(n2, 0, spin_multiplicity, rotated_orbitals);
    energy = new_energy;
    wavefunction = new_wavefunction;
    std::cout << "New Energy after rotation: " << std::fixed
              << std::setprecision(10) << energy << " Hartree" << std::endl;
    std::cout << std::endl;
  }

  std::cout << "\nFinal Energy: " << std::fixed << std::setprecision(10)
            << energy << " Hartree" << std::endl;
  std::cout << "Final stability status: " << (is_stable ? "stable" : "unstable")
            << std::endl;

Python API

# Create N2 molecule at stretched geometry (1.4 Angstrom)
symbols = ["N", "N"]
coords = np.array([[0.0, 0.0, 0.0], [1.4, 0.0, 0.0]])
coords *= ANGSTROM_TO_BOHR
n2 = Structure(symbols, coords)

# Create and configure SCF solver with auto scf_type
scf_solver = create("scf_solver", "qdk")
scf_solver.settings().set("scf_type", "auto")
scf_solver.settings().set("method", "hf")

# Run initial SCF calculation
spin_multiplicity = 1
energy, wavefunction = scf_solver.run(
    n2, charge=0, spin_multiplicity=1, basis_or_guess="def2-svp"
)
print(f"Initial SCF Energy: {energy:.10f} Hartree")

# Determine if calculation is restricted and configure stability checker accordingly
is_restricted = wavefunction.get_orbitals().is_restricted() and spin_multiplicity == 1
if is_restricted:
    stability_checker.settings().set("external", True)
else:
    stability_checker.settings().set("external", False)

# Iterative workflow: Check stability and rotate orbitals until convergence
max_iterations = 5
iteration = 0

print("\n=== Starting Iterative Stability Workflow ===\n")

while iteration < max_iterations:
    iteration += 1
    print(f"--- Iteration {iteration} ---")
    print(f"Current Energy: {energy:.10f} Hartree")

    # Check stability - handle potential Davidson convergence failure
    try:
        is_stable, result = stability_checker.run(wavefunction)
    except RuntimeError as e:
        if "Davidson Did Not Converge!" in str(e):
            print("Try increasing max_subspace or adjusting davidson_tolerance")
            raise RuntimeError(f"Davidson solver did not converge: {str(e)}")
        else:
            raise RuntimeError(f"Stability check failed: {str(e)}")

    if is_stable:
        print("\nConverged to stable wavefunction!")
        break

    # Determine rotation type based on which instability is present
    do_external = False
    if not result.is_internal_stable():
        smallest_eigenvalue, rotation_vector = (
            result.get_smallest_internal_eigenvalue_and_vector()
        )
        print(
            f"Internal instability detected. Smallest eigenvalue: {smallest_eigenvalue:.6f}"
        )
    elif not result.is_external_stable() and result.has_external_result():
        smallest_eigenvalue, rotation_vector = (
            result.get_smallest_external_eigenvalue_and_vector()
        )
        print(
            f"External instability detected. Smallest eigenvalue: {smallest_eigenvalue:.6f}"
        )
        do_external = True
    else:
        print("Unexpected state: neither internal nor external instability detected")
        break

    # Rotate orbitals along the instability direction
    num_alpha, num_beta = wavefunction.get_total_num_electrons()
    orbitals = wavefunction.get_orbitals()
    rotated_orbitals = rotate_orbitals(
        orbitals, rotation_vector, num_alpha, num_beta, do_external
    )

    # If external instability detected, switch to unrestricted calculation
    if do_external:
        print("Switching to unrestricted calculation due to external instability")
        # Create new solver instances with updated settings
        scf_solver_name = scf_solver.name()
        stability_checker_name = stability_checker.name()

        # Copy settings and update for unrestricted calculation
        scf_settings_map = scf_solver.settings().to_dict()
        scf_settings_map["scf_type"] = "unrestricted"
        new_scf_solver = create("scf_solver", scf_solver_name)
        new_scf_solver.settings().from_dict(scf_settings_map)

        stability_settings_map = stability_checker.settings().to_dict()
        stability_settings_map["external"] = False
        new_stability_checker = create("stability_checker", stability_checker_name)
        new_stability_checker.settings().from_dict(stability_settings_map)

        scf_solver = new_scf_solver
        stability_checker = new_stability_checker

    # Re-run SCF with rotated orbitals as initial guess
    energy, wavefunction = scf_solver.run(
        n2, charge=0, spin_multiplicity=1, basis_or_guess=rotated_orbitals
    )
    print(f"New Energy after rotation: {energy:.10f} Hartree")
    print()

print(f"\nFinal Energy: {energy:.10f} Hartree")
print(f"Final stability status: {is_stable}")

Available settings

The StabilityChecker accepts settings to control stability analysis behavior. All implementations share a common base set of settings:

Setting	Type	Default	Description
internal	bool	True	Check internal stability (within same wavefunction type)
external	bool	Implementation-dependent	Check external stability (RHF → UHF instabilities). Only supported for RHF wavefunctions
method	string	"hf"	The method to match the wavefunction: "hf" for Hartree-Fock, or a DFT functional name
stability_tolerance	float	-1e-4	Eigenvalue threshold for determining stability. Eigenvalues below this value indicate instability
davidson_tolerance	float	1e-8	Convergence tolerance for the Davidson eigenvalue solver
max_subspace	int	80	Maximum subspace dimension for the Davidson solver

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

Available implementations

QDK/Chemistry’s StabilityChecker provides a unified interface to stability analysis across various quantum chemistry packages. You can discover available implementations programmatically:

C++ API

  auto available = StabilityCheckerFactory::available();
  std::cout << "\nAvailable stability checker implementations: ";
  for (const auto& name : available) {
    std::cout << name << " ";
  }
  std::cout << std::endl;

Python API

from qdk_chemistry.algorithms import registry  # noqa: E402

print("available backend choices for scf_solver and stability_checker:")
print(registry.available("stability_checker"))

QDK (Native)

Factory name: "qdk" (default)

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

Capabilities

• Internal stability analysis for RHF and UHF wavefunctions

• External stability analysis for RHF wavefunctions (RHF → UHF instabilities)

• Efficient Davidson eigenvalue solver for computing lowest eigenvalues

• Support for both Hartree-Fock and DFT wavefunctions

Settings

The QDK implementation uses the common base settings listed in the Available settings section above. The default value for external is False for this implementation.

PySCF

Factory name: "pyscf"

The PySCF plugin provides access to stability analysis through the PySCF quantum chemistry package.

Capabilities

• Internal stability analysis for RHF, restricted open-shell Hartree-Fock (ROHF), and UHF wavefunctions

• External stability analysis for RHF wavefunctions (RHF → UHF instabilities)

• Support for both Hartree-Fock and DFT wavefunctions

• Point group symmetry considerations in stability analysis

• Configurable Davidson solver parameters and PySCF verbosity

• Note: ROHF wavefunctions are not supported by the QDK backend for now

Settings

The PySCF implementation supports all common base settings. The following table shows settings that differ from or extend the base implementation:

Setting	Type	Default	Description
external	bool	True	Check external stability (differs from QDK default of False)
with_symmetry	bool	False	Whether to respect point group symmetry during analysis
nroots	int	3	Number of eigenvalue roots to compute in the Davidson solver
xc_grid	int	3	DFT integration grid density level (0=coarse, 9=very fine)
pyscf_verbose	int	4	PySCF verbosity level for logging (0=silent, 4=info, 5=debug)

Note

The PySCF implementation automatically detects the wavefunction type (RHF, ROHF, or UHF) and applies the appropriate stability analysis method. External stability analysis is only supported for RHF wavefunctions and will raise an error if requested for ROHF or UHF.

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

Understanding stability results

The StabilityResult object returned by the stability checker contains detailed information about the stability analysis:

Overall stability status

The boolean return value indicates whether the wavefunction is stable overall (both internally and externally if checked).

Internal stability

The is_internal_stable() method indicates whether internal stability is satisfied. Internal instabilities suggest the energy can be lowered while maintaining the same wavefunction type.

External stability

The is_external_stable() method indicates whether external stability is satisfied (if external stability was checked). External instabilities in RHF calculations indicate that breaking spin symmetry (switching to UHF) would lower the energy.

Eigenvalues and eigenvectors

The stability result provides access to the computed eigenvalues and their corresponding eigenvectors:

• get_internal_eigenvalues() and get_internal_eigenvectors()

• get_external_eigenvalues() and get_external_eigenvectors()

• get_smallest_eigenvalue() returns the overall smallest eigenvalue

• get_smallest_internal_eigenvalue_and_vector() and get_smallest_external_eigenvalue_and_vector() provide convenient access to the most unstable direction

The eigenvector corresponding to the smallest (most negative) eigenvalue indicates the direction in orbital space that would most effectively lower the energy. This eigenvector can be used with the rotate_orbitals() utility function to generate rotated orbitals for use as an initial guess in a subsequent SCF calculation.

Related classes

• ScfSolver: Performs SCF calculations that produce wavefunctions for stability analysis

• Wavefunction: Input wavefunction to analyze

• Orbitals: Can be rotated based on stability analysis results

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

• [SM91]: Original reference for stability analysis in quantum chemistry