Quickstart

These quickstart instructions provide a high-level overview of the typical workflow when using the QDK/Chemistry library. More comprehensive documentation can be found in the In-depth user guide. A list of features can be found in the Features, methods and dependencies document.

Installation

To install QDK/Chemistry, please see the installation instructions.

End-to-end example

This document is intended to provide a brief introduction to the QDK/Chemistry library by walking through a minimal end-to-end example for ground state energy estimation with state preparation and measurement. The emphasis of this example is optimization: reducing the resources required for the quantum computer to run a simple chemistry application. The example starts with a molecular structure and ends with an energy estimation computed by simulating a quantum circuit. The focus of this example is on high-level concepts and common coding patterns that can be extended to other applications.

You can also view a related code sample in a Jupyter notebook format, along with several other examples, in the examples folder of the GitHub repository.

Create a Structure object

The Structure class represents a molecular structure, i.e. the coordinates of its atoms. Structure objects can be constructed manually, or via deserialization from a file. QDK/Chemistry supports multiple serialization formats for Structure objects, including the standard XYZ file format, as well as QDK/Chemistry-specific JSON and HDF5 serialization schemes. Internally to QDK/Chemistry, coordinates are always stored in Bohr/atomic units; however, when reading or writing to files, the units follow the file format convention (Angstrom for XYZ) or can be specified (JSON, HDF5). See below for language specific examples of creating and serializing Structure objects.

  // Load para-benzyne structure from XYZ file
  auto structure = std::make_shared<data::Structure>(
      data::Structure::from_xyz_file("../data/para_benzyne.structure.xyz"));

  std::cout << "Created structure with " << structure->get_num_atoms()
            << " atoms" << std::endl;
  std::cout << "Elements: ";
  for (const auto& elem : structure->get_elements()) {
    std::cout << data::Structure::element_to_symbol(elem) << " ";
  }
  std::cout << std::endl;
from pathlib import Path
import numpy as np
from qdk_chemistry.algorithms import create
from qdk_chemistry.data import Structure
from qdk_chemistry.data.qubit_hamiltonian import (
    filter_and_group_pauli_ops_from_wavefunction,
)
from qdk_chemistry.utils.wavefunction import get_top_determinants

# Load para-benzyne structure from XYZ file
structure = Structure.from_xyz_file(
    Path(__file__).parent / "../data/para_benzyne.structure.xyz"
)

print(f"Created structure with {structure.get_num_atoms()} atoms")
print(f"Elements: {structure.get_elements()}")

Run a self-consistent field (SCF) calculation

Once a Structure is created, an SCF calculation can be performed to produce an initial Wavefunction as well as an SCF energy. QDK/Chemistry performs SCF calculations via instantiations of a Self-consistent field (SCF) solver algorithm, and is the first instance of the separation Data classes and Algorithm classes most will encounter in QDK/Chemistry. See the design principles documentation for more information on this pattern and how data flow is generally treated in QDK/Chemistry. Instantiations of the Self-consistent field (SCF) solver algorithm (and all other Algorithm classes) are managed by a factory. See the Factory pattern documentation for more information on how it is used in the code base.

The inputs for an SCF calculation are a Structure object, the total charge and spin multiplicity of the molecular system, and information about the single-particle basis to be used. Optionally, Settings specific to the particular Self-consistent field (SCF) solver can be configured to control the execution of the SCF algorithm (e.g. convergence tolerances, etc) by accessing the settings() method. The basis for the SCF calculation can be set via a string input (specifying one of the Available basis sets), a custom BasisSet or initial Orbitals can also be provided.

  // Perform a SCF calculation
  auto scf_solver = algorithms::ScfSolverFactory::create();
  auto [E_hf, wfn_hf] = scf_solver->run(structure, 0, 1, "cc-pvdz");
  std::cout << "SCF energy is " << E_hf << " Hartree" << std::endl;

  // Display a summary of the molecular orbitals
  std::cout << "SCF Orbitals:\n"
            << wfn_hf->get_orbitals()->get_summary() << std::endl;
# Perform an SCF calculation, returning the energy and wavefunction
scf_solver = create("scf_solver")
E_hf, wfn_hf = scf_solver.run(
    structure, charge=0, spin_multiplicity=1, basis_or_guess="cc-pvdz"
)
print(f"SCF energy is {E_hf:.3f} Hartree")

# Display a summary of the molecular orbitals obtained from the SCF calculation
print("SCF Orbitals:\n", wfn_hf.get_orbitals().get_summary())

Select an active space

While a full set of SCF orbitals are useful for many applications, they are often not the optimal set for further post-SCF calculations, including algorithms intended for a quantum computer. For this reason, both orbital localization and active space selection algorithms are provided within QDK/Chemistry.

QDK/Chemistry offers many methods for the selection of active spaces to reduce the problem size: accurately modeling the quantum many-body problem while avoiding the prohibitive computational scaling of full configuration interaction. See the Active space selection documentation for a list of supported methods, along with their associated Settings, which accompany the standard QDK/Chemistry distribution.

The following are language-specific examples of how to select a so-called “valence” active space containing a subset of only those orbitals surrounding the Fermi level - in this case 6 electrons to be distributed in 6 orbitals (6e, 6o).

  // Select active space (6 electrons in 6 orbitals)
  auto active_space_selector =
      algorithms::ActiveSpaceSelectorFactory::create("qdk_valence");
  active_space_selector->settings().set("num_active_electrons", 6);
  active_space_selector->settings().set("num_active_orbitals", 6);

  auto active_wfn = active_space_selector->run(wfn_hf);
  auto active_orbitals = active_wfn->get_orbitals();

  // Print a summary of the active space orbitals
  std::cout << "Active Space Orbitals:\n"
            << active_orbitals->get_summary() << std::endl;
# Select active space (6 electrons in 6 orbitals for para-benzyne)
#   to choose most chemically relevant orbitals
active_space_selector = create(
    "active_space_selector",
    algorithm_name="qdk_valence",
    num_active_electrons=6,
    num_active_orbitals=6,
)
active_wfn = active_space_selector.run(wfn_hf)
active_orbitals = active_wfn.get_orbitals()

# Print a summary of the active space orbitals
print("Active Space Orbitals:\n", active_orbitals.get_summary())

Calculate the Hamiltonian

Once an active space has been selected, the electronic Hamiltonian can be computed within that active space to describe the energetic interactions between electrons. QDK/Chemistry provides flexible Hamiltonian construction capabilities through the Hamiltonian construction algorithm. The Hamiltonian constructor generates the one- and two-electron integrals (or factorizations thereof) needed for subsequent quantum many-body calculations and quantum algorithms.

  // Construct Hamiltonian in the active space
  auto hamiltonian_constructor =
      algorithms::HamiltonianConstructorFactory::create();
  auto hamiltonian = hamiltonian_constructor->run(active_orbitals);
  std::cout << "Active Space Hamiltonian:\n"
            << hamiltonian->get_summary() << std::endl;
# Construct Hamiltonian in the active space and print its summary
hamiltonian_constructor = create("hamiltonian_constructor")
hamiltonian = hamiltonian_constructor.run(active_orbitals)
print("Active Space Hamiltonian:\n", hamiltonian.get_summary())

Compute a multi-configurational wavefunction for the active space

With the active space Hamiltonian constructed, quantum many-body calculations can be performed to obtain multi-configurational wavefunctions that go beyond the single-determinant SCF approximation. QDK/Chemistry supports various Multi-Configuration (MC) methods including Complete Active Space Configuration Interaction (CASCI) and selected CI approaches. MC calculations are performed via instantiations of the Multi-configuration calculations algorithm, which takes as input an instance of an active space Hamiltonian, and the number of alpha and beta electrons, and produces as output a Wavefunction representing the multi-configurational state as well as its associated energy.

While multi-configurational methods provide more accurate energy estimates than SCF, their primary role in the quantum applications workflow is to generate high-quality initial states for quantum algorithms. On scaled fault-tolerant quantum computers, these classically-computed wavefunctions serve as the foundation for state preparation circuits, enabling algorithms such as quantum phase estimation to achieve chemical accuracy for systems that remain intractable for purely classical methods. These methods also serve as a critical analysis tool, allowing users to better understand the electronic structure of their systems of interest and the potential for quantum algorithms to provide meaningful utility over classical state of the art.

In the following example, as the aforementioned (6e, 6o) active space is relatively small, we perform a CASCI calculation to obtain the exact ground state wavefunction within the active space.

  // Perform CASCI calculation
  auto mc = algorithms::MultiConfigurationCalculatorFactory::create();
  auto [E_cas, wfn_cas] = mc->run(hamiltonian, 3, 3);
  std::cout << "CASCI energy is " << E_cas
            << " Hartree, and the electron correlation energy is "
            << (E_cas - E_hf) << " Hartree" << std::endl;
# Perform CASCI calculation to get the wavefunction and exact energy for the active space
mc = create("multi_configuration_calculator")
E_cas, wfn_cas = mc.run(
    hamiltonian, n_active_alpha_electrons=3, n_active_beta_electrons=3
)
print(
    f"CASCI energy is {E_cas:.3f} Hartree, and the electron correlation energy is {E_cas - E_hf:.3f} Hartree"
)

Select important configurations

For large active spaces, the multi-configuration Wavefunction may contain thousands or millions of configurations, but often only a small subset contributes significantly to the overall state. By identifying and retaining only the dominant configurations (those with the largest amplitudes), we can create a sparse wavefunction that maintains high fidelity with respect to the original Wavefunction while dramatically reducing resource requirements for quantum state preparation.

/**
 * @brief Get top N determinants from a wavefunction sorted by coefficient
 * magnitude
 * @param wfn The wavefunction to extract determinants from
 * @param max_determinants Maximum number of determinants to return
 * @return Vector of top configurations
 */
std::vector<data::Configuration> get_top_determinants(
    std::shared_ptr<data::Wavefunction> wfn, size_t max_determinants) {
  const auto& determinants = wfn->get_active_determinants();
  const auto& coeffs = wfn->get_coefficients();

  // Extract coefficients as doubles (handle real case)
  std::vector<double> coeff_values;
  const auto& real_coeffs = std::get<Eigen::VectorXd>(coeffs);
  coeff_values.resize(real_coeffs.size());
  for (int i = 0; i < real_coeffs.size(); ++i) {
    coeff_values[i] = std::abs(real_coeffs[i]);
  }

  // Create indices and sort by coefficient magnitude
  std::vector<size_t> indices(determinants.size());
  std::iota(indices.begin(), indices.end(), 0);
  std::sort(indices.begin(), indices.end(),
            [&coeff_values](size_t i1, size_t i2) {
              return coeff_values[i1] > coeff_values[i2];
            });

  // Extract top determinants
  std::vector<data::Configuration> top_dets;
  size_t n = std::min(max_determinants, determinants.size());
  top_dets.reserve(n);
  for (size_t i = 0; i < n; ++i) {
    top_dets.push_back(determinants[indices[i]]);
  }

  return top_dets;
}
  // Get top 2 determinants from the CASCI wavefunction
  auto top_configurations = get_top_determinants(wfn_cas, 2);

  // Perform PMC calculation with selected configurations
  auto pmc = algorithms::ProjectedMultiConfigurationCalculatorFactory::create();
  auto [E_pmc, wfn_pmc] = pmc->run(hamiltonian, top_configurations);
  std::cout << "Reference energy for top 2 determinants is " << E_pmc
            << " Hartree" << std::endl;
# Get top 2 determinants from the CASCI wavefunction to form a sparse wavefunction
top_configurations = get_top_determinants(wfn_cas, max_determinants=2)

# Compute the reference energy of the sparse wavefunction
pmc_calculator = create("projected_multi_configuration_calculator")
E_sparse, wfn_sparse = pmc_calculator.run(hamiltonian, list(top_configurations.keys()))

print(f"Reference energy for top 2 determinants is {E_sparse:.6f} Hartree")

Preparing a qubit representation of the Hamiltonian

To execute quantum algorithms on actual quantum hardware, the Fermionic Hamiltonian must be mapped to a qubit Hamiltonian composed of Pauli operators. QDK/Chemistry supports multiple encoding schemes including the Jordan-Wigner and Bravyi-Kitaev transformations. The resulting qubit Hamiltonian is represented as a sum of Pauli strings (tensor products of Pauli \(X, Y, Z\), and identity operators), each with an associated coefficient.

For efficient energy estimation, the qubit Hamiltonian can be optimized in two ways:

  1. Filter out terms that have negligible expectation values given the sparse Wavefunction, pre-computing their classical contributions.

  2. Group the remaining Pauli operators into commuting sets that can be measured simultaneously, significantly reducing the number of measurement circuits required on the quantum device.

// This step is currently only available in the Python API.
# Prepare qubit Hamiltonian
qubit_mapper = create("qubit_mapper", algorithm_name="qiskit", encoding="jordan-wigner")
qubit_hamiltonian = qubit_mapper.run(hamiltonian)

# Print the number of Pauli strings in the full Hamiltonian
print(
    f"Number of Pauli strings in the Hamiltonian: {len(qubit_hamiltonian.pauli_strings)}"
)

# Filter and group Pauli operators based on the wavefunction
filtered_hamiltonian_ops, classical_coeffs = (
    filter_and_group_pauli_ops_from_wavefunction(qubit_hamiltonian, wfn_sparse)
)
print(
    f"Filtered and grouped qubit Hamiltonian contains {len(filtered_hamiltonian_ops)} groups:"
)
for igroup, group in enumerate(filtered_hamiltonian_ops):
    print(f"Group {igroup + 1}: {[group.pauli_strings]}")
print(f"Number of classical coefficients: {len(classical_coeffs)}")

Generate the state preparation circuit

Given the classical representation of the sparse multi-configurational Wavefunction, a quantum circuit can be generated to prepare this state on a quantum computer. This can be done in many ways, including via Isometry encoding [ICK+16], linear combinations of unitaries, and tensor product methods. However, when the wavefunction is very sparse, these methods can be inefficient. In QDK/Chemistry, we provide a specialized method for generating state preparation circuits for sparse wavefunctions based on the construction of sparse isometries over GF(2) with X gates, provided as a State preparation algorithm. See State preparation for more details on the other algorithms provided for state preparation in QDK/Chemistry.

// This step is currently only available in the Python API.

# Generate state preparation circuit for the sparse state via sparse isometry (GF2 + X)
state_prep = create("state_prep", "sparse_isometry_gf2x")
sparse_isometry_circuit = state_prep.run(wfn_sparse)

Estimate the ground state energy by sampling

The final step combines the state preparation circuit with the measurement circuits derived from the grouped qubit Hamiltonian to estimate the ground state energy. Each measurement circuit appends specific Pauli basis rotations to the state preparation circuit, followed by computational basis measurements. The measurement outcomes are repeated many times (shots) to build up statistics for each observable group.

The energy expectation value is computed by combining the measurement statistics from each group with their corresponding Hamiltonian coefficients and the pre-computed classical contributions. The statistical nature of quantum measurements introduces variance in the energy estimate, which decreases as the number of shots is increased.

// This step is currently only available in the Python API.
# Estimate energy using the optimized circuit and filtered Hamiltonian operators
estimator = create("energy_estimator", algorithm_name="qdk_base_simulator")
energy_results, simulation_data = estimator.run(
    circuit=sparse_isometry_circuit,
    qubit_hamiltonians=filtered_hamiltonian_ops,
    total_shots=250000,
    classical_coeffs=classical_coeffs,
)

for i, results in enumerate(simulation_data.bitstring_counts):
    print(
        f"Measurement Results for Hamiltonian Group {i + 1}: {simulation_data.hamiltonians[i].pauli_strings}"
    )

# Print statistics for the measured energy
energy_mean = energy_results.energy_expectation_value + hamiltonian.get_core_energy()
energy_stddev = np.sqrt(energy_results.energy_variance)
print(
    f"Estimated energy from quantum circuit: {energy_mean:.3f} ± {energy_stddev:.3f} Hartree"
)

# Print comparison with reference energy
print(f"Difference from reference energy: {energy_mean - E_sparse} Hartree")

Additional examples

For more information, see:

  • Above examples as complete C++ and Python scripts.

  • Additional examples of workflows using QDK/Chemistry in the examples directory of the source repository.

  • In-depth user guide for links to more detailed documentation on specific components of QDK/Chemistry.