Performing a Standard BEAST-VQE Convergence Study

In this example, we plot the convergence of a BEAST Variational Quantum Eigensolver (BEAST-VQE) calculation.

Full Script

The script below shows the full script, which are then explained section by section.

  1"""Demonstrate BEAST-VQE convergence plotting on LiH molecule."""
  2
  3# pyright: reportUnknownVariableType=false, reportMissingImports=false
  4# ruff: noqa
  5from pathlib import Path
  6
  7import matplotlib.pyplot as plt
  8import numpy as np
  9
 10import qrunch as qc
 11
 12
 13
 14def plot_energies(
 15    beast_energies: list[float],
 16    reference_energy: float,
 17    *,
 18    outdir: Path = Path("dist"),
 19    show: bool = False,
 20) -> None:
 21    """
 22    Plot VQE energy convergence and log-scale error vs reference energy.
 23
 24    Args:
 25        beast_energies: Sequence of total energies per VQE iteration.
 26        reference_energy: Reference total energy (e.g., FCI) to compare against.
 27        outdir: Output directory where plots are written. Defaults to "dist".
 28        show: Whether to display the figures interactively after saving.
 29
 30    """
 31    # Convert to a NumPy array for safe numeric operations and validation.
 32    energies = np.asarray(list(beast_energies), dtype=float)
 33
 34    # Basic validation to help catch silent failures early.
 35    if energies.size == 0:
 36        msg = "beast_energies must contain at least one value."
 37        raise ValueError(msg)
 38    if not np.isfinite(energies).all():
 39        msg = "beast_energies contains non-finite values."
 40        raise ValueError(msg)
 41    if not np.isfinite(reference_energy):
 42        msg = "reference_energy must be finite."
 43        raise ValueError(msg)
 44
 45    # Prepare x-axis as 1-based iteration indices for readability.
 46    iterations = np.arange(1, energies.size + 1, dtype=int)
 47
 48    # Ensure output directory exists.
 49    outdir.mkdir(parents=True, exist_ok=True)
 50
 51    # -----------------------------
 52    # Figure 1: Energies vs Iteration
 53    # -----------------------------
 54    fig1 = plt.figure()
 55    plt.plot(iterations, energies, "-*", color="darkgreen", label="BEAST-VQE")  # points to visualize steps
 56    plt.axhline(reference_energy, linestyle="--", color="black", label="Reference")
 57    plt.xlabel("Iteration number")
 58    plt.ylabel("Total energy [Hartree]")
 59    plt.title("BEAST-VQE convergence")
 60    plt.legend()
 61    plt.tight_layout()
 62
 63    convergence_plot = outdir / "vqe_convergence.png"
 64    fig1.savefig(convergence_plot, dpi=200, bbox_inches="tight")
 65
 66    # -----------------------------
 67    # Figure 2: |Energy - Reference| (log scale)
 68    # -----------------------------
 69    # Use absolute error; add a tiny epsilon to avoid log(0) if we hit the reference exactly.
 70    eps = np.finfo(float).eps
 71    abs_error = np.abs(energies - reference_energy) + eps
 72
 73    fig2 = plt.figure()
 74    plt.semilogy(iterations, abs_error, "-*", color="darkgreen", label="|E - E_ref|")
 75    plt.xlabel("Iteration number")
 76    plt.ylabel("Absolute error [Hartree]")
 77    plt.title("BEAST-VQE error vs reference")
 78    plt.legend()
 79    plt.tight_layout()
 80
 81    error_convergence_plot = outdir / "vqe_error_semilogy.png"
 82    fig2.savefig(error_convergence_plot, dpi=200, bbox_inches="tight")
 83
 84    if show:
 85        plt.show()
 86    else:
 87        # Close figures to free memory when running in batch contexts.
 88        plt.close(fig1)
 89        plt.close(fig2)
 90
 91
 92
 93
 94def main() -> list[float]:
 95    """Run BEAST-VQE on LiH as the first script."""
 96    path_to_molecule_xyz_file = Path().absolute() / "qrunch/demo_scripts/lih.xyz"
 97    # Build Lithium Hydride (LiH) molecular configuration from xyz file.
 98    molecular_configuration = qc.build_molecular_configuration(
 99        molecule=Path(path_to_molecule_xyz_file),
100        basis_set="sto3g",
101        spin_difference=0,
102        charge=0,
103        units="angstrom",
104    )
105
106    # Build the ground state problem.
107    problem_builder = qc.problem_builder_creator().ground_state().standard().create()
108    ground_state_problem = problem_builder.build_restricted(molecular_configuration)
109
110    adaptive_vqe_options = qc.options.IterativeVqeOptions(
111        max_iterations=100,  # Increase the maximum number of iterations
112        force_all_iterations=True,  # Force all iterations to run, independent of convergence.
113    )
114
115    # Build the BEAST-VQE calculator using the user-configured VQE instance
116    beast_vqe_calculator = (
117        qc.calculator_creator()  # Start creating a calculator
118        .vqe()  # Narrow: pick Variational quantum eigensolver (VQE)
119        .iterative()  # Narrow: pick the iterative VQE
120        .beast()  # Narrow: pick the iterative VQE
121        .with_options(adaptive_vqe_options)  # Use the user-defined iterative VQE options
122        .choose_minimizer()  # Start sub-choice: Choose the gate parameter minimizer
123        .quick_default()  # Perform the sub-choice selection - Here we pick the greedy and quick minimizer
124        .create()  # Create the calculator instance
125    )
126
127    result = beast_vqe_calculator.calculate(ground_state_problem)
128
129    # We can get a nice timings report.
130
131    # Extract BEAST-VQE energies for plotting
132    beast_energies = result.total_energy_per_macro_iteration_with_initial_energy_and_final_energy.values
133
134    # Build a Full Configuration Interaction (FCI) calculator
135    full_configuration_interaction_calculator = (
136        qc.calculator_creator()  # Start creating a calculator
137        .configuration_interaction()  # Narrow: pick Configuration Interaction (CI)
138        .paired_electron_approximation()  # Narrow: pick paired CI
139        .create()  # Create the calculator instance
140    )
141
142    # Calculate the paired-FCI result for reference energy
143    pfci_result = full_configuration_interaction_calculator.calculate(ground_state_problem)
144
145    # Plot the convergence and error relative to FCI reference.
146    plot_energies(beast_energies, reference_energy=pfci_result.total_energy.value, show=True)
147
148    return beast_energies
149
150
151if __name__ == "__main__":
152    main()

Explanation

  1. Import the QDK library

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

import qrunch as qc

This code snippet loads the python packages we need. import qrunch as qc loads the public and stable API. This is the only import you need for most tasks, and the only one that is guaranteed to be supported across versions. In addition, we use matplotlib, numpy, and Path.

  1. The plotting method

def plot_energies(
    beast_energies: list[float],
    reference_energy: float,
    *,
    outdir: Path = Path("dist"),
    show: bool = False,
) -> None:
    """
    Plot VQE energy convergence and log-scale error vs reference energy.

    Args:
        beast_energies: Sequence of total energies per VQE iteration.
        reference_energy: Reference total energy (e.g., FCI) to compare against.
        outdir: Output directory where plots are written. Defaults to "dist".
        show: Whether to display the figures interactively after saving.

    """
    # Convert to a NumPy array for safe numeric operations and validation.
    energies = np.asarray(list(beast_energies), dtype=float)

    # Basic validation to help catch silent failures early.
    if energies.size == 0:
        msg = "beast_energies must contain at least one value."
        raise ValueError(msg)
    if not np.isfinite(energies).all():
        msg = "beast_energies contains non-finite values."
        raise ValueError(msg)
    if not np.isfinite(reference_energy):
        msg = "reference_energy must be finite."
        raise ValueError(msg)

    # Prepare x-axis as 1-based iteration indices for readability.
    iterations = np.arange(1, energies.size + 1, dtype=int)

    # Ensure output directory exists.
    outdir.mkdir(parents=True, exist_ok=True)

    # -----------------------------
    # Figure 1: Energies vs Iteration
    # -----------------------------
    fig1 = plt.figure()
    plt.plot(iterations, energies, "-*", color="darkgreen", label="BEAST-VQE")  # points to visualize steps
    plt.axhline(reference_energy, linestyle="--", color="black", label="Reference")
    plt.xlabel("Iteration number")
    plt.ylabel("Total energy [Hartree]")
    plt.title("BEAST-VQE convergence")
    plt.legend()
    plt.tight_layout()

    convergence_plot = outdir / "vqe_convergence.png"
    fig1.savefig(convergence_plot, dpi=200, bbox_inches="tight")

    # -----------------------------
    # Figure 2: |Energy - Reference| (log scale)
    # -----------------------------
    # Use absolute error; add a tiny epsilon to avoid log(0) if we hit the reference exactly.
    eps = np.finfo(float).eps
    abs_error = np.abs(energies - reference_energy) + eps

    fig2 = plt.figure()
    plt.semilogy(iterations, abs_error, "-*", color="darkgreen", label="|E - E_ref|")
    plt.xlabel("Iteration number")
    plt.ylabel("Absolute error [Hartree]")
    plt.title("BEAST-VQE error vs reference")
    plt.legend()
    plt.tight_layout()

    error_convergence_plot = outdir / "vqe_error_semilogy.png"
    fig2.savefig(error_convergence_plot, dpi=200, bbox_inches="tight")

    if show:
        plt.show()
    else:
        # Close figures to free memory when running in batch contexts.
        plt.close(fig1)
        plt.close(fig2)

This code plots the convergence of the VQE energy and the error relative to the reference. The paired Full Configuration Interaction (pFCI) energy.

The code generates two PNG files under dist/: - vqe_convergence.png - vqe_error_semilogy.png

  1. Molecule from XYZ file

    # Build Lithium Hydride (LiH) molecular configuration from xyz file.
    molecular_configuration = qc.build_molecular_configuration(
        molecule=Path(path_to_molecule_xyz_file),
        basis_set="sto3g",
        spin_difference=0,
        charge=0,
        units="angstrom",
    )

We load LiH from an XYZ file, set the STO-3G basis, neutral charge, and zero spin difference.

  1. Build the ground-state problem

    # Build the ground state problem.
    problem_builder = qc.problem_builder_creator().ground_state().standard().create()
    ground_state_problem = problem_builder.build_restricted(molecular_configuration)

We create a standard ground-state problem using the problem builder. Here we use the restricted variant for LiH This is required for the BEAST-VQE algorithm, and the Paired-FCI reference calculation.

  1. Configure and create the BEAST-VQE calculator

    adaptive_vqe_options = qc.options.IterativeVqeOptions(
        max_iterations=100,  # Increase the maximum number of iterations
        force_all_iterations=True,  # Force all iterations to run, independent of convergence.
    )

    # Build the BEAST-VQE calculator using the user-configured VQE instance
    beast_vqe_calculator = (
        qc.calculator_creator()  # Start creating a calculator
        .vqe()  # Narrow: pick Variational quantum eigensolver (VQE)
        .iterative()  # Narrow: pick the iterative VQE
        .beast()  # Narrow: pick the iterative VQE
        .with_options(adaptive_vqe_options)  # Use the user-defined iterative VQE options
        .choose_minimizer()  # Start sub-choice: Choose the gate parameter minimizer
        .quick_default()  # Perform the sub-choice selection - Here we pick the greedy and quick minimizer
        .create()  # Create the calculator instance
    )

We specify that we want a maximum of 100 iterations, one gate per iteration, and we force all iterations to be executed, independent of any stopping criteria.

We then build the adaptive VQE with the specified options and wrap it in the BEAST-VQE calculator, where we choose a quick gate parameter optimization method, optimizing only the parameter of the last gate in every iteration.

  1. Run and extract energies

    result = beast_vqe_calculator.calculate(ground_state_problem)

We run the VQE calculation, and here we also print a timing report for profiling:

    # We can get a nice timings report.

  1. Compute a high-accuracy paired-FCI reference

    # Build a Full Configuration Interaction (FCI) calculator
    full_configuration_interaction_calculator = (
        qc.calculator_creator()  # Start creating a calculator
        .configuration_interaction()  # Narrow: pick Configuration Interaction (CI)
        .paired_electron_approximation()  # Narrow: pick paired CI
        .create()  # Create the calculator instance
    )

    # Calculate the paired-FCI result for reference energy
    pfci_result = full_configuration_interaction_calculator.calculate(ground_state_problem)

We build a paired-FCI calculator and calculate the paired-FCI energy for the same problem.

  1. Make the plots

    # Build a Full Configuration Interaction (FCI) calculator
    full_configuration_interaction_calculator = (
        qc.calculator_creator()  # Start creating a calculator
        .configuration_interaction()  # Narrow: pick Configuration Interaction (CI)
        .paired_electron_approximation()  # Narrow: pick paired CI
        .create()  # Create the calculator instance
    )

    # Calculate the paired-FCI result for reference energy
    pfci_result = full_configuration_interaction_calculator.calculate(ground_state_problem)

We plot the VQE convergence and the error relative to the paired-FCI reference, by calling the plotting function defined above.

Running the Example

After saving the script as run_vqe_convergence.py, you can run it directly from the command line:

$ python run_vqe_convergence.py

You should see 2 plots:

BEAST-VQE convergence plot BEAST-VQE convergence plot, as the error relative to paired-FCI