Exploring Qiskit (23) — Quantum Chemistry Calculations with the Latest Qiskit

Purpose

The IBM Quantum Challenge Fall 2022 event was held in the winter of 2022, and it included:

Lab 4: Quantum Chemistry

There was a problem concerning the hydrogen atom H, which roughly involves solving $\mathcal{H} \psi = E \psi$ for the Hamiltonian:

$$\begin{align*} \mathcal{H} = - \frac{\hbar^2}{2 \mu} \Delta - \frac{Z}{|x|} \end{align*}$$

In particular, if we are interested in the ground state eigenvalue $E_1$:

• Second-quantize $\mathcal{H}$,
• Change it to a spin operator form through something like the Jordan-Wigner transformation (we'll call this $\mathcal{H}$ as well),
• Use an appropriate ansatz (this time using the UCCSD ansatz) $\ket{\psi (\theta)} = U(\theta) \ket{0}^{\otimes n}$,

By minimizing:

$$\begin{align*} f(\theta) = \braket{\psi (\theta) | \mathcal{H} | \psi (\theta)} \geq E_1 \end{align*}$$

it is expected that $\inf \braket{\psi (\theta) | \mathcal{H} | \psi (\theta)} = E_1$ will approximately hold.

While this is quite a fascinating topic, the purpose here is to get it working again, as the previous code no longer runs with the latest Qiskit.

Previous Code

A notebook containing the problem creator's solutions is available at solutions-by-authors/lab-4, so I'll gratefully make use of it. It seems to be under the Apache-2.0 license.

%%bash

pip install "qiskit==1.0.2" "qiskit[visualization]==1.0.2"
pip install "qiskit_aer==0.13.3" "qiskit_algorithms==0.3.0" "qiskit_nature==0.7.2"
pip install "pyscf==2.5.0"
pip list | grep qiskit

qiskit 1.0.2
qiskit-aer 0.13.3
qiskit-algorithms 0.3.0
qiskit-nature 0.7.2

Under this environment, several APIs and modules have been deprecated and will not work.

I rewrote the code by referring to migration guides such as the (Minimal) EigensolverFactory Migration Guide and Problem Transformers.

Most of the changes are located within construct_problem.

Updated Code

For now, I've got it working as follows.

First, it seems enough to import something like this. Many components have been consolidated under qiskit_nature.second_q, so I followed that.

# Import necessary libraries and packages
import math
import matplotlib.pyplot as plt
import numpy as np

import warnings
warnings.filterwarnings('ignore')

from qiskit import QuantumCircuit
from qiskit.primitives import Estimator
from qiskit_aer import StatevectorSimulator

# Import Qiskit libraries for VQE
from qiskit_algorithms import MinimumEigensolverResult, VQE
from qiskit_algorithms.optimizers import SLSQP, SPSA
from qiskit_algorithms.utils import algorithm_globals
algorithm_globals.random_seed = 1024

# Import Qiskit Nature libraries
from qiskit_nature.second_q.algorithms import GroundStateEigensolver
from qiskit_nature.second_q.algorithms.initial_points import HFInitialPoint
from qiskit_nature.second_q.circuit.library import UCC, UCCSD, HartreeFock
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import (
    BravyiKitaevMapper, JordanWignerMapper, ParityMapper
)
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer

from qiskit_nature.settings import settings

settings.dict_aux_operators = True

1. Helper Function for constructing problem

In the new approach, it seems PySCFDriver is used directly, so I rewrote it by referring to the migration guide and how parameters were passed in the original code. Originally, ParityMapper was used, but I tried using JordanWignerMapper instead.

def construct_problem(
    geometry, charge, multiplicity, basis, num_electrons, num_molecular_orbitals
):

    # Convert to a format like 'H 0.0 0.0 0.0; H 0.0 0.0 0.735'
    atom = "; ".join(
        [" ".join([atm[0]] + [str(v) for v in atm[1]])
            for atm in geometry]
    )
    # https://github.com/qiskit-community/qiskit-nature/blob/0.5.2/qiskit_nature/drivers/second_quantization/pyscfd/pyscfdriver.py#L335
    spin = multiplicity - 1
    driver = PySCFDriver(atom=atom, charge=charge, spin=spin, basis=basis)

    # Run the preliminary quantum chemistry calculation
    problem = driver.run()

    # Set the active space
    active_space_trafo = ActiveSpaceTransformer(
        num_electrons=num_electrons,
        num_spatial_orbitals=num_molecular_orbitals
    )

    # Now you can get the reduced electronic structure problem
    problem_reduced = active_space_trafo.transform(problem)

    # The second quantized Hamiltonian of the reduce problem
    second_q_ops_reduced = problem_reduced.hamiltonian.second_q_op()

    # Set the mapper to qubits
    # Tried using JordanWignerMapper
    mapper = JordanWignerMapper()

    # Compute the Hamitonian in qubit form
    qubit_op = mapper.map(second_q_ops_reduced)
    aux_ops = {}
    aux_ops.update(
        mapper.map(problem_reduced.properties.particle_number.second_q_ops())
    )
    aux_ops.update(
        mapper.map(problem_reduced.properties.angular_momentum.second_q_ops())
    )

    ansatz = UCCSD(
        problem_reduced.num_spatial_orbitals,
        problem_reduced.num_particles,
        mapper,
        initial_state=HartreeFock(
            problem_reduced.num_spatial_orbitals,
            problem_reduced.num_particles,
            mapper,
        ),
    )

    initial_point = HFInitialPoint()
    initial_point.ansatz = ansatz
    initial_point.problem = problem_reduced

    solver = VQE(Estimator(), ansatz, SLSQP())
    solver.initial_point = initial_point.to_numpy_array()
    result = solver.compute_minimum_eigenvalue(qubit_op, aux_ops)
    real_solution = result.optimal_value

    return ansatz, qubit_op, real_solution, problem_reduced

2. Helper function for Running VQE

This is the standard VQE optimization part, and since there was no particular need to change it, I used it as is. It calculates expectation values using the modern Estimator and optimizes parameters using SPSA. While it might be related to running in a cloud environment, it seems EstimatorV2 is being introduced, so that might affect things again in the future.

def custom_vqe(
    estimator, ansatz, ops, problem_reduced, optimizer = None,
    initial_point=None,
):

    # Define convergence list
    convergence = []

    # Keep track of jobs (Do-not-modify)
    job_list = []

    # Define evaluate_expectation function
    def evaluate_expectation(x):
        x = list(x)

        # Define estimator run parameters
        job = estimator.run(
            circuits=[ansatz], observables=[ops], parameter_values=[x]
        ).result()
        results = job.values[0]
        job_list.append(job)

        # Pass results back to callback function
        return np.real(results)

    # Call back function
    def callback(x,fx,ax,tx,nx):
        # Callback function to get a view on internal states and statistics of the optimizer for visualization
        convergence.append(evaluate_expectation(fx))

    np.random.seed(10)

    # Define initial point. We shall define a random point here based on the number of parameters in our ansatz
    if initial_point is None:
        initial_point = np.random.random(ansatz.num_parameters)

    # Define optimizer and pass callback function
    if optimizer == None:
        optimizer = SPSA(maxiter=100, callback=callback)

    # Define minimize function
    result =  optimizer.minimize(evaluate_expectation, x0=initial_point)

    vqe_interpret = []
    for i in range(len(convergence)):
        sol = MinimumEigensolverResult()
        sol.eigenvalue = convergence[i]
        sol = problem_reduced.interpret(sol).total_energies[0]
        vqe_interpret.append(sol)

    return vqe_interpret, job_list, result

3. Custom Function to Plot Graphs

This function visualizes the optimization history, and I used this as is as well.

import matplotlib.pyplot as plt

def plot_graph(energy, real_solution, molecule, color="tab:blue"):

    plt.rcParams["font.size"] = 14

    # plot loss and reference value
    plt.figure(figsize=(12, 6), facecolor='white')
    plt.plot(energy, label="Estimator VQE {}".format(molecule),color = color)
    plt.axhline(y=real_solution.real, color="tab:red", ls="--", label="Target")

    plt.legend(loc="best")
    plt.xlabel("Iteration")
    plt.ylabel("Energy [H]")
    plt.title("VQE energy")
    plt.show()

1. Compute the Ground State Energy of Each Ingredient

Define Geometry

It simply places a proton at the origin of a 3D space. I used this as is.

# Constructing H
hydrogen_a = [["H", [0.0 ,0.0, 0.0]]]

Construct Problem

H - VQE Run

Computation using VQE. I didn't change much here either, but since the eigenvalues and other values are in numpy.complex128, I took the real part for appearance's sake.

algorithm_globals.random_seed = 1024
# For H
ansatz_a, ops_a, real_solution_a, problem_reduced_a = \
    construct_problem(geometry=hydrogen_a, charge=0, multiplicity=2,
                      basis="ccpvdz", num_electrons=(1,0),
                      num_molecular_orbitals=2)

# Estimator VQE for H
Energy_H_a,_,jobs = custom_vqe(estimator=Estimator(), ansatz=ansatz_a,
                               ops=ops_a, problem_reduced=problem_reduced_a)


Energy_H_a = [v.real for v in Energy_H_a]
real_solution_a = real_solution_a.real

# Plot Graph H
plot_graph(Energy_H_a, real_solution_a, "H",color = "tab:purple")

Viewing the Results

The ground state energy obtained via VQE is calculated in units of Hartree energy. According to the CODATA Value: Hartree energy in eV, we can use the following conversion:

27.211 386 245 988(53) eV

So:

# Energy of the ground state of the hydrogen atom.: -13.6 [eV]
E = -13.6

# 1 hartree-energy = 27.211386245988 eV
print(np.min(Energy_H_a) * 27.211386245988)

-13.586057479730476

This is the result. While I am not sure about the exact level of precision, it seems to be generally consistent.

Summary

Since I haven't tested other cases, I'm not entirely sure if the code migration has been handled perfectly, but at least the hydrogen atom case is now working.

The parts that failed were the steps up to converting the Hamiltonian into spin operator form. This was caused by changes in the qiskit-nature API, stemming from updates to what was formerly known as qiskit-terra as well as its own internal changes. Since I am not an expert in this area at all, it was a struggle to tediously replace the code while consulting the migration guide.

For general experimentation, it might be easier to just use the versions and coding styles that were functional at the time.