Get started with Qunova HI-VQE Chemistry

See the overview and the API reference

• Package versions

  The code on this page was developed using the following requirements. We recommend using these versions or newer.

  qiskit-ibm-runtime~=0.45.0
  

Set up and load the function

Ways to use HI-VQE

The instance you authenticate with needs Quantum Functions access enabled, or the function will not appear in the catalog.

To get a license, request one from Qunova; IBM then enables the function on your IBM Cloud account for the term of your license. With trial access, runs are limited to function size S; see Function sizes.

To ask about a license, you can also email [email protected]. A licensed run needs no license token.

Load the function

Save your IBM Quantum® API key to your local environment if you have not already, then select the quantum function as follows:

import reprlib
from qiskit_ibm_catalog import FunctionsCatalog

catalog = FunctionsCatalog(channel="ibm_quantum_platform")

# Verify that you have access to the function
catalog.list()

Output:

[QiskitFunction(qunova/hivqe-chemistry),
 QiskitFunction(global-data-quantum/quantum-portfolio-optimizer),
 QiskitFunction(algorithmiq/tem),
 QiskitFunction(qedma/qesem),
 QiskitFunction(multiverse/singularity),
 QiskitFunction(ibm/circuit-function),
 QiskitFunction(q-ctrl/optimization-solver),
 QiskitFunction(colibritd/quick-pde),
 QiskitFunction(q-ctrl/performance-management),
 QiskitFunction(kipu-quantum/iskay-quantum-optimizer)]

# Load the function
function = catalog.load("qunova/hivqe-chemistry")

# Use the least busy QPU available to you
backend_name = catalog.least_busy(operational=True).name

Function sizes

Every job runs at one size ΓÇö S, M or L. Too small and the job is rejected; larger than needed and it runs on a bigger compute profile than the work requires.

The size depends on the number of qubits the molecule needs and the max_states you run with. HI-VQE uses two qubits for each orbital of the active space, one for each electron spin:

By default, every orbital in the basis is active. You can choose a smaller active space with molecule_options.

A problem must fit both bounds of a size, so either its qubits or its max_states can raise the size it needs. Runs with trial access are limited to S; a problem that needs a larger size is rejected with error code 4500.

The default is M, used when a job does not request a size.

You can also read the mapping off the function itself with function.sizes_map, and the size used when a job does not ask for one with function.default_size.

Calculate the size before submitting

Run this on your own inputs. It reads them and nothing else: no job is submitted and nothing is reserved, so you can call it as often as you like. It needs only PySCF (pip install pyscf). It is the same calculation the function runs when a job arrives, so a job sized by it is not rejected for being too small.

Open the code below, then run it in a cell, or save it as hivqe_chemistry_size_calculator.py and import it with from hivqe_chemistry_size_calculator import calculate_function_size. You only need to call calculate_function_size, with the same max_states you will pass to run(); the rest of the code is how the size is worked out.

• Size calculator code (run this before calculate_function_size)

  # Size calculator for Qunova HI-VQE Chemistry. Requires PySCF only.
  #
  # calculate_function_size(geometry, max_states=..., molecule_options=...) returns "S", "M"
  # or "L" from two things: the number of qubits the problem needs, and the
  # max_states it is run with. A problem is the smallest size whose bounds it fits
  # on both; above any bound of "M" it is "L".
  #
  # The qubits are two per active orbital. The active space is the orbitals AVAS
  # selects if avas_aolabels is set, else active_orbitals if set, else every
  # orbital in the basis. With fcidump_file_name, the path of an FCIDUMP file, it
  # is every orbital of the file, and geometry and molecule_options are ignored,
  # as the function ignores them. max_states is the function's argument of that
  # name. Other inputs are accepted and ignored.
  import gzip
  import re
  from typing import Any, Literal
  
  from pyscf import gto, scf
  from pyscf.mcscf import avas
  
  
  Size = Literal["S", "M", "L"]
  _ORDER = {"S": 0, "M": 1, "L": 2}
  
  # Upper bounds (inclusive) for each size, on the number of qubits and on
  # max_states. A problem above either bound of a size does not fit it; above
  # every bound here it is "L".
  _BOUNDS = (
      ("S", 20, 10**3),
      ("M", 48, 10**7),
  )
  
  # The function's defaults for the molecule_options that decide the active space.
  _DEFAULTS = {
      "active_orbitals": None,
      "avas_aolabels": None,
      "avas_threshold": 0.2,
      "basis": "sto-3g",
      "charge": 0,
      "symmetry": False,
      "symmetry_subgroup": None,
      "unit": "angstrom",
      "nucmod": None,
      "pseudo": None,
      "cart": False,
      "magmom": None,
  }
  
  # AVAS weights below this are numerical noise; counting them would make the size vary
  # from one call to the next.
  _AVAS_WEIGHT_TOLERANCE = 1e-10
  
  
  def calculate_function_size(problem: Any, **kwargs: Any) -> Size:
      """Return the size ("S", "M" or "L") required to run the function."""
      metrics = _measure(problem, **kwargs)
      metrics["max_states"] = _max_states(kwargs)
      return size_for(metrics)
  
  
  def size_for(metrics: dict) -> Size:
      """The smallest size whose qubit and max_states bounds both fit the metrics."""
      for size, max_qubits, max_states in _BOUNDS:
          if (
              metrics["num_qubits"] <= max_qubits
              and metrics["max_states"] <= max_states
          ):
              return size
      return "L"
  
  
  def _max_states(kwargs: dict) -> int:
      """max_states as the function takes it: a positive integer, which it requires."""
      value = kwargs.get("max_states")
      # JSON Schema counts 1000.0 as an integer, so the function can receive one.
      if isinstance(value, float) and value.is_integer():
          value = int(value)
      if isinstance(value, bool) or not isinstance(value, int) or value < 1:
          raise ValueError("max_states must be given as a positive integer.")
      return value
  
  
  def _measure(problem: Any, **kwargs: Any) -> dict:
      """Return the number of qubits the run uses for this geometry and these options."""
      if kwargs.get("fcidump_file_name"):
          return {
              "num_qubits": 2
              * _fcidump_num_orbitals(kwargs["fcidump_file_name"])
          }
      options = {**_DEFAULTS, **(kwargs.get("molecule_options") or {})}
      mol = gto.M(
          atom=problem,
          basis=options["basis"],
          spin=0,  # the function runs singlets only
          charge=options["charge"],
          symmetry=options["symmetry"],
          symmetry_subgroup=options["symmetry_subgroup"],
          unit=options["unit"],
          nucmod=options["nucmod"],
          pseudo=options["pseudo"],
          cart=options["cart"],
          magmom=options["magmom"],
          verbose=0,
      )
  
      if options["avas_aolabels"] is None:
          # The active space is the orbitals given, or the whole basis; building the molecule
          # is enough to know it.
          if options["active_orbitals"]:
              return {"num_qubits": 2 * len(options["active_orbitals"])}
          return {"num_qubits": 2 * mol.nao_nr()}
  
      # AVAS picks the active space from an RHF, so it has to run one.
      rhf = scf.RHF(mol)
      rhf.verbose = 0
      rhf.kernel()
      selection = avas.AVAS(
          mf=rhf,
          aolabels=options["avas_aolabels"],
          threshold=max(options["avas_threshold"], _AVAS_WEIGHT_TOLERANCE),
      )
      selection.verbose = 0
      selection.kernel()
      return {"num_qubits": 2 * selection.ncas}
  
  
  # An FCIDUMP header is read up to its end (&END or /). It is a few lines; its
  # only part that grows with the problem is ORBSYM, a few characters per orbital.
  _FCIDUMP_MAX_HEADER_LINES = 1000
  
  
  def _fcidump_num_orbitals(path: str) -> int:
      """NORB from the header of an FCIDUMP file, plain or gzip-compressed,
      without reading the integrals."""
      with open(path, "rb") as raw:
          gzipped = raw.read(2) == b"\x1f\x8b"
      opener = gzip.open if gzipped else open
      header = []
      try:
          with opener(path, "rt") as in_file:
              for line_number, line in enumerate(in_file):
                  header.append(line.upper())
                  if "&END" in header[-1] or "/" in line:
                      break
                  if line_number >= _FCIDUMP_MAX_HEADER_LINES:
                      raise ValueError(
                          "No end of header (&END or /) was found in the FCIDUMP file."
                      )
              else:
                  raise ValueError(
                      "No end of header (&END or /) was found in the FCIDUMP file."
                  )
      except (UnicodeDecodeError, OSError, EOFError) as error:
          raise ValueError(
              "The FCIDUMP file could not be read as text or as gzip-compressed text."
          ) from error
      match = re.search(r"NORB\s*=\s*(\d+)", "".join(header))
      if match is None or int(match.group(1)) < 1:
          raise ValueError(
              "No valid NORB entry was found in the FCIDUMP file header."
          )
      return int(match.group(1))

For example, H2 in the STO-3G basis has 2 orbitals, so 4 qubits:

# 4 qubits and max_states=1000 fit "S"
print(calculate_function_size("H 0 0 0; H 0 0 0.74", max_states=1000))

Pass the result to run() as function_size, as the Example does. If the size is too small, the job fails with error code 4600 and the message names the size to use instead. If it is larger than needed, the job runs normally and records a 4601 advisory event.

Example: NH3 from a geometry

The first example shows how to compute the ground state energy for an NH3 molecule using the HI-VQE algorithm.

Define the molecule

The molecular geometry of NH3 is provided with Cartesian coordinates separated with ";" for each atom.

# Define the molecule geometry
geometry = """
N         -0.85188       -0.02741        0.03141;
H          0.16545        0.00593       -0.01648;
H         -1.16348       -0.39357       -0.86702;
H         -1.16348        0.94228        0.06281;
"""

Additional options can be defined and provided for molecular system in the following dictionary format. For every option you can set, with its type, default and valid range, see molecule_options and hivqe_options in the API reference.

# Configure some options for the job.
molecule_options = {"basis": "sto3g"}
hivqe_options = {"shots": 100, "max_iter": 20}

Check the arguments

Before running, you can check the arguments with validate_arguments(). It checks their names, types and ranges against the function's argument schema without creating a job. It returns {"valid": True}, or raises an exception describing what failed. function.arguments_schema shows every argument the function accepts, with its type, range and default.

This check does not cover the function size; see Function sizes.

function.validate_arguments(
    {
        "geometry": geometry,
        "backend_name": backend_name,
        "max_states": 1000,
        "max_expansion_states": 10,
        "molecule_options": molecule_options,
        "hivqe_options": hivqe_options,
    }
)

Run the function

Run the function with the geometry and options. NH3 has 8 orbitals in the STO-3G basis, so the problem uses 16 qubits. Run with max_states=1000, it fits S, so it also runs with trial access. Calculate the size, then pass it to run().

# 16 qubits and max_states=1000 fit "S"
size = calculate_function_size(
    geometry, max_states=1000, molecule_options=molecule_options
)
print(size)

# Run HI-VQE
job = function.run(
    geometry=geometry,
    # `backend_name` is a backend you have access to,
    # for example, "ibm_marrakesh".
    backend_name=backend_name,
    # `function_size` sets the compute profile the job runs on.
    # See "Function sizes" above.
    function_size=size,
    max_states=1000,
    max_expansion_states=10,
    molecule_options=molecule_options,
    hivqe_options=hivqe_options,
)

It is a good idea to print the Function job ID so that it can be provided in support requests if something goes wrong.

print("Job ID:", job.job_id)

Output:

Job ID: 140a0768-5863-4727-8074-aaf8a15aedc1

Get the results

Check your quantum function workload's status or return results as follows:

print(job.status())

Output:

INITIALIZING

result = job.result()

# Output can be long, so we display a shortened representation
shortened_result = reprlib.repr(result)
print(shortened_result)

Output:

{'energy': -55.52089378327241, 'energy_history': [-55.51854305148325, -55.51916155061359, -55.52001948770051, -55.52010922102457, -55.52019547601691, -55.52089891983981, ...], 'energy_variance': 3.26661973062643e-10, 'function_version': '5.1.0', ...}

To access the ground state energy, use the "energy" key. The "results_file_name" key provides the name of the h5 results file that can be downloaded with catalog.file_download().

The response also carries:

• "energy_variance": the energy variance of the final trial state, as a single number,
• "energy_history": the energies computed along the way, two per iteration,
• "function_version": the version of the function that ran the job, which is worth quoting in a support request, and
• "metadata": a "resource_usage" dictionary giving the time spent in each stage of the job, including QPU time.

The results file contains:

• "final_energy": the energy of the final trial state,
• "final_energy_variance": the energy variance of the final trial state,
• "energy_history": the energy history over HI-VQE's training,
• "eigenvector": the CI coefficients of the final trial state's wave function, and
• "states": the corresponding electron configurations comprising the final trial state wave function.

fci_energy = -55.521148034704126  # the exact energy using FCI method
hivqe_energy = result["energy"]
print(
    f"|Exact Energy - HI-VQE Energy|: "
    f"{abs(fci_energy - hivqe_energy) * 1000} mHa"
)

Output:

|Exact Energy - HI-VQE Energy|: 0.25425143171275977 mHa

Download the results file

Download the results file for the calculation with catalog.file_download(), then open it as shown below. Reading it requires the h5py package, which you can install with pip install h5py. See Manage files across runs for how files are moved to and from a function.

import h5py

# Save the file in the current directory under its own name
results_file_name = result["results_file_name"]
catalog.file_download(
    results_file_name, function, target_name=results_file_name
)

# Read data from the results file for a single calculation
with h5py.File(results_file_name, "r") as f:
    states = f["states"].asstr()[:]
    eigenvector = f["eigenvector"][:]
    energy_history = f["energy_history"][:]
    final_energy = f["final_energy"][0]
    final_energy_variance = f["final_energy_variance"][0]

Read the job log

job.logs() gives the job log, which this function writes one record per line: the time since the run started, the stage of the run, and the record. The stages are start, mapping (setting up the problem), iter N+ and iter N- (the two energy evaluations of SPSA iteration N), iter N (its convergence check), final and done.

print(job.logs())

Each evaluation gives two lines: its energies, then its time by step. On a QPU, the time line also gives how long the Qiskit Runtime job queued and ran, and the QPU time it used; the done line gives the totals for the run. This excerpt is from a run of the example on a simulator:

00:00:00.0  start    HI-VQE Quantum Function version 5.1.0
00:00:00.0  start    random seed 1586140065; to have this run repeated, give it to Qunova support together with your input
00:00:00.2  mapping  16 qubits, 10 active electrons, 3136 states in the symmetry space
00:00:00.7  iter 1+  sample -55.4729814691 Ha (144 states), subspace -55.4729814691 Ha (144 states); ...
00:00:00.7  iter 1+  time 0.37 s = optimizing 0.22 + waiting 0.01 + executing 0.06 + post-processing 0.08 + status updates 0.00; CPU 0.53 s
00:00:00.7  iter 1   energy -55.4729814691 Ha
...
00:00:01.4  final    energy -55.5166284474 Ha, variance 3.11e-08, 676 states; ...
00:00:01.4  done     time 1.27 s = mapping 0.17 + optimizing 0.42 + waiting 0.05 + executing 0.15 + post-processing 0.48 + status updates 0.00; CPU 2.89 s; ...

Results of a stopped job

If you stop a job with job.cancel(), its result holds what the run had computed by then, so job.result() still returns it:

• "status": "canceled",
• "canceled_during": the job status the run was in when it was stopped,
• "iterations_completed" and "energy_history": the iterations finished and the energies computed in them,
• "lowest_energy": the lowest of those energies, or None if no iteration finished,
• "runtime_jobs": the Qiskit Runtime jobs the run made, and what became of each, and
• "function_version" and "metadata", as for a completed job.

A stopped job writes no results file.

Example: Run from an FCIDUMP file

Instead of a geometry, you can give the problem as an FCIDUMP file of its integrals, for example one written by PySCF or by another quantum chemistry program. Upload the file to the function, then pass its name as fcidump_file_name:

• The integrals are used as they are. No mean-field calculation is run, and geometry and molecule_options are ignored if you pass them as well.
• Every orbital in the file is active, so the problem uses 2 ├ù NORB qubits.
• The file can be plain text or gzip-compressed. Only closed-shell singlets are supported: MS2=0 and an even NELEC.
• The first NELEC/2 orbitals are taken as the occupied orbitals of the reference state. If the orbitals of the file are not in order of energy, the job log says so, as the run can converge slowly from that state.

This writes the NH3 problem of the example to a file with PySCF, which you can install with pip install pyscf:

import pyscf
from pyscf.tools import fcidump

mol = pyscf.M(atom=geometry, basis="sto3g")
fcidump.from_scf(mol.RHF().run(), "nh3.FCIDUMP")

Calculate the size from the file, then upload it and run. Pass the file name alone, without a directory. The file stays uploaded until you remove it with catalog.file_delete().

# nh3.FCIDUMP has 8 orbitals, so 16 qubits; with max_states=1000 it fits "S"
size = calculate_function_size(
    None, fcidump_file_name="nh3.FCIDUMP", max_states=1000
)
print(size)

catalog.file_upload("nh3.FCIDUMP", function)

job = function.run(
    fcidump_file_name="nh3.FCIDUMP",
    backend_name=backend_name,
    function_size=size,
    max_states=1000,
    max_expansion_states=10,
    hivqe_options=hivqe_options,
)

Fetch error messages

A failed job's error is fetched as described in Fetch error messages. For example, this run fails after it starts:

job = function.run(
    geometry="invalid-geometry",  # This will cause an error
    backend_name=backend_name,
    max_states=1000,
    max_expansion_states=15,
    molecule_options=molecule_options,
    hivqe_options=hivqe_options,
)

job.result()

Output:

---------------------------------------------------------------------------
QiskitServerlessException                 Traceback (most recent call last)
Cell In[12], line 10
      1 job = function.run(
      2     geometry="invalid-geometry",  # This will cause an error
      3     backend_name=backend_name,
   (...)
      7     hivqe_options=hivqe_options,
      8 )
---> 10 job.result()

File ~/work/documentation/documentation/.tox/py311/lib/python3.11/site-packages/qiskit_serverless/core/job.py:253, in Job.result(self, wait, cadence, verbose, maxwait)
    251 if self.status() == "ERROR":
    252     if results:
--> 253         raise QiskitServerlessException(results)
    255     raise QiskitServerlessException(self.filtered_logs(include=r"(?i)error|exception"))
    257 if isinstance(results, str):

QiskitServerlessException: ["runner.UnknownRuntimeError: 'An unexpected error occurred during job execution. Please make sure that your inputs are valid. If you are still experiencing problems, you can contact the Qunova Computing support service at [email protected] and provide the Function job ID of this job for more assistance. -- https://docs.quantum.ibm.com/errors#1500'\n"]

job.status()

Output:

'ERROR'

job.error_message() returns the error code and the message of the failed job. When the input was not accepted, the job log ends with a line pointing to it.

print(job.error_message())

Every failure reports an error code, a message describing what went wrong, and a solution saying what to do about it. The codes follow IBM's error code registry, and the ones this function returns are:

When a Qiskit Runtime job made by the run fails with an error code of its own, such as 1305, that code is reported instead; look it up in the registry.

Several codes cover more than one cause, so read the message as well as the code. If you need help with a failure, send the Function job ID to [email protected].

A call with a wrong argument name, type or range does not get this far: run() refuses it before a job is created; see Check the arguments.

Next steps