About IBM Quantum Challenge Fall 2022

The last IBM IBM Quantum Challenge Fall 2022 was really great. Participants dove into different quantum computing topics as a space-ship Captain.

We got through a spectacular space story which, I think, got everyone engaged.

The First Lab was about Quantum Primitives. Primitives are meant to serve as an foundational, elementary, building blocks for users to perform quantum computations, developers to implement quantum algorithms, and researchers to solve complex problems and deliver new applications.[1] Examples of them are Samplers (for sampling probability distributions) and Estimators (for estimating values). We solved 6 graded exercises.

The Second Lab showed up with Quantum Kernel Learning. Topics covered were data encoding, error mitigation, estimating fidelity states and Quantum Support Vector Machines (QSVM). We had to deal with 5 graded problems.

From where I stand, the third lab was the most difficult. Solving the Travel SaleMan Problem (TSP) with parameterized quantum circuits was really a Quantum Optimization Challenge. This problem derives its importance from its "hardness" and ubiquitous equivalence to other relevant combinatorics optimization problems that arise in practice. [2] The largest one, 9 graded exercises to solve.

The last one, Quantum Chemistry, was an opportunity to know what kind of real problems quantum computing can handle. Chemical Reactions, Dipole Moments and Complex Molecular Geometries were the tasks to overcome. Six problems to deal with: 5 graded exercises and a final scored challenge. Let's talk about it.

The final challenge was to compute the energy reaction of cyclopropenylidene (C3H2). Cyclopropenylidene, or c-C3H2, is a partially aromatic molecule belonging to a highly reactive class of organic molecules known as carbenes. On Earth, cyclopropenylidene is only seen in the laboratory due to its reactivity. However, cyclopropenylidene is found in significant concentrations in the interstellar medium (ISM) and on Saturn's moon Titan.[3] Reaction is showed in Figure 1. Atomic Carbon (C) reacts with Acetylene (C2H2) to obtain Cyclopropenylidene. This reaction process is known as a key reaction process as it provides a starting point to a mechanism for the growth of carbon chains. Also, this process occurs with almost zero potential barrier at the room temperature, and it is believed that this reaction also occurs fast in low temperature environments like in interstellar space[4]. Let's see the reference value that IBM Quantum provided.

Let's see how we can approach to this problem.

From a classical point of view, the reaction is breaking one C-C bond and building up two C-C bonds in products. Thus, the net outcome is one C-C bond. As bonds in products releases energy, this amount of energy is negative. From tables, we can get the average bond value as 348 kJ/mol (3,607 eV). However, let's dive into Qiskit to show how we can handle this reaction.

First of all, let's import all modules and libraries that we are going to use:

# Import necessary libraries and packages

import math
import matplotlib.pyplot as plt
import numpy as np

import warnings
warnings.filterwarnings('ignore')

from qiskit.primitives import Estimator, BackendEstimator
from qiskit.providers.aer import StatevectorSimulator
from qiskit.utils import QuantumInstance
from qiskit.circuit.library import EfficientSU2

from qiskit.tools.jupyter import *
from qiskit.visualization import *

# Import Qiskit libraries for VQE
from qiskit.algorithms import MinimumEigensolverResult, VQE
from qiskit.algorithms.optimizers import SLSQP, SPSA


# Import Qiskit Nature libraries
from qiskit_nature.algorithms import GroundStateEigensolver
from qiskit_nature.circuit.library import HartreeFock
from qiskit_nature.drivers import Molecule
from qiskit_nature.drivers.second_quantization import ElectronicStructureDriverType, ElectronicStructureMoleculeDriver, ElectronicStructureProblem, ActiveSpaceTransformer, QubitConverter, ParityMapper

# Import Utils
from qiskit.utils import algorithm_globals
from qiskit.providers.fake_provider import FakeLagos
from qiskit.providers.aer.noise import NoiseModel
from qiskit.utils.mitigation import CompleteMeasFitter
from qiskit.opflow.primitive_ops import Z2Symmetries

# Prototype-zne / Noisy Simulator
!pip install prototype-zne --quie
from qiskit_nature.settings import settings
from zne import ZNEStrategy
from zne.extrapolation import PolynomialExtrapolator, LinearExtrapolator
from zne.noise_amplification import LocalFoldingAmplifier, GlobalFoldingAmplifier
settings.dict_aux_operators = True        

Let's define our seed:

algorithm_globals.random_seed = 1024        

Let's define a function to diagonalize Hamiltonians:

def exact_diagonalizer(problem, converter)
  solver = NumPyMinimumEigensolverFactory()
  calc = GroundStateEigensolver(converter, solver)
  result = calc.solve(problem)
  return result        

Let's define a function to plot Ground State Energies:

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()        

Right now, we have to construct our problem. We will re factor the function provided in the lab. We are going to introduce molecular geometries to transform them to Fermionic Operators and then, by a mapper, we get the number of qubits and observables of every chemical compound involved in the reaction. Furthermore, we will leverage on Z2 Symmetries to taper qubits:

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


# Molecule Problem Definition

        molecule = Molecule(geometry=geometry, charge=charge, multiplicity=multiplicity)
        ?? ?
        driver = ElectronicStructureMoleculeDriver(molecule, basis=basis, driver_type=ElectronicStructureDriverType.PYSCF)
        
        properties = driver.run()
        ?????? ?
        trafo = ActiveSpaceTransformer(num_electrons=num_electrons, num_molecular_orbitals=num_molecular_orbitals)
        ?????? ?
        problem_reduced = ElectronicStructureProblem(driver, transformers=[trafo])

# Qubit Tapering

parity_mapper = ParityMapper()
?????? ???????? 
parity_converter = QubitConverter(parity_mapper, z2symmetry_reduction='auto')
???? ???????? 
second_q_ops_reduced = problem_reduced.second_q_ops()?

qubit_op_parity = parity_converter.convert(second_q_ops_reduced.get('ElectronicEnergy'), num_particles=problem_reduced.num_particles, sector_locator=problem_reduced.symmetry_sector_locator)?
????? ???????? 
return parity_converter, qubit_op_parity, problem_reduced        

After constructing our problem, we have to deal with the different ansatz's. Every compound needs a 'guesstimate', so let's code a function to get an ansatz for reactants and products. Also, let's compose the ansatz with a HartreeFock Initial State and then transpile the EfficientSU2 gates with basic gates like U and CNOT's. Furthermore, a Variational Quantum Eigensolver (VQE) with some optimization parameters is needed, to pass it to the Ground State EigenSolver to determine the real solution of the problem. Also, we will call the function exact_diagonalize() to get the exact eigenvalues:

def ansatz_effSU2(parity, ops, problem):
    #Ansatz Initial State 
    init_state = HartreeFock(problem.num_spin_orbitals, problem.num_particles, parity)
    ansatz_temp = EfficientSU2(ops.num_qubits, reps=1, entanglement='linear')
    ansatz_temp.compose(init_state, front=True, inplace=True)
        _ansatz = transpile(ansatz_temp, basis_gates=['u', 'cx'], optimization_level=3))
    # VQE 
    vqe_solver = aVQE(ansatz=_ansatz, optimizer=SLSQP(), quantum_instance=StatevectorSimulator(seed_simulator=42), include_custom=True, expectation=AerPauliExpectation())
    # GSE Solver
    solver = GroundStateEigensolver(parity, vqe_solver)
    real_solution = solver.solve(problem)
    # Final Ansatz
    ansatz = vqe_solver.ansatz
    # Diagonalized Hamiltonian
    result_exact = exact_diagonalizer(problem, parity)
    exact_energy = np.real(result_exact.eigenenergies[0])
          
    return ansatz, real_solution, exact_energyy        

Now, we will call the function custom_vqe() provided by IBM to leverage the optimal points:

def custom_vqe(estimator, ansatz, ops, problem, maxiter_for_opt)

??? # 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([ansatz], ops, 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))

??? # Define initial point.
??? initial_point = np.random.random(ansatz.num_parameters)

??? # Define optimizer and pass callback function
??? optimizer = SPSA(maxiter=maxiter_for_opt, 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.interpret(sol).total_energies[0]
??????? vqe_interpret.append(sol)

??? return vqe_interpret, job_list, result        

After that, we are going to define the geometry for products and reactants:

# Define geometry for products and reactant

#Products
# https://atct.anl.gov/Thermochemical%20Data/version%201.122/species/?species_number=442

Cyclopropenylidene =? [["C", [2.2883,??? 0.6993,??? 0.3468 ]],
?????????? ["C",[??? 1.9543,??? 2.0133,??? 0.7806]],
?????????? ["C",[??? 1.0108,??? 0.9522,??? 0.6802]],
?????????? ["H",[??? 3.0291,??? 0.0000,??? 0.0000]],
?????????? ["H",[??? 0.0000,??? 0.5997,??? 0.7904]]]


#Reactants
carbon = [["C",[0.0,0.0,0.0]]]

# https://webbook.nist.gov/cgi/cbook.cgi?Name=acetylene&Units=SI
acetylene = [["C", [0.0000,??? 0.0000,?? -0.6025]],
??????????? ["H",[??? 0.0000,??? 0.0000,?? -1.6691]],
??????????? ["C",[??? 0.0000,??? 0.0000,??? 0.6025]],
??????????? ["H",[??? 0.0_000,??? 0.0000,??? 1.6691]]]s        

Now, we will call the construct_problem_refactored() function to get the parity, the observables and reduced problem variables. Basis Set is going to be 'aug-ccpvdz'.

# Sample defininition Cyclopropenylidene
parity_cy, ops_cy, problem_reduced_cy = construct_problem_refactored(geometry=Cyclopropenylidene, charge=0, multiplicity=1, basis="aug-ccpvdz", num_electrons=2, num_molecular_orbitals=2)

# Sample definition Carbon
parity_c, ops_c, problem_reduced_c = construct_problem_refactored(geometry=carbon, charge=0, multiplicity=3, basis="aug-ccpvdz", num_electrons=4, num_molecular_orbitals=4)

# Sample definition Acetylene
parity_ac, ops_ac, problem_reduced_ac = construct_problem_refactored(geometry=acetylene, charge=0, multiplicity=1, basis='aug-ccpvdz', num_electrons=4, num_molecular_orbitals=4)        

Let's move on an call the ansatz function to calculate them:

ansatz_cy, real_solution_cy, exact_energy_cy = ansatz_effSU2(parity_cy, ops_cy, problem_reduced_cy
ansatz_c, real_solution_c, exact_energy_c = ansatz_effSU2(parity_c, ops_c, problem_reduced_c)
ansatz_ac, real_solution_ac, exact_energy_ac = ansatz_effSU2(parity_ac, ops_ac, problem_reduced_ac))        

Now, let's analyze the outputs:

print(ansatz_cy.num_parameters, ansatz_c.num_parameters, ansatz_ac.num_parameters)        

8, 12, 12

print(len(ops_cy), len(ops_c), len(ops_ac))        

9, 8, 20

Let's check if we have leveraged all the Z2 Symmetries:

print(Z2Symmetries.find_Z2_symmetries(ops_cy)
print(Z2Symmetries.find_Z2_symmetries(ops_c))
print(Z2Symmetries.find_Z2_symmetries(ops_ac))        

Cyclopropendylene Z2 Symmetries: Symmetries: Single-Qubit Pauli X: Cliffords: Qubit index:[] Tapering values: - Possible values: []

Carbon Z2 Symmetries: Symmetries:ZIZ IXI Single-Qubit Pauli X:IIXIZI Cliffords: 0.7071067811865475 * ZIZ+ 0.7071067811865475 * IIX0.7071067811865475 * IXI+ 0.7071067811865475 * IZI Qubit index:[0, 1] Tapering values: - Possible values: [1, 1], [1, -1], [-1, 1], [-1, -1]

Acetylene Z2 Symmetries: Symmetries: Single-Qubit Pauli X: Cliffords: Qubit index:[] Tapering values: - Possible values: []

Oops! Carbon seems to have another sector to leverage which the 'auto' parameter doesn't seem to cover. Maybe Sensei Yukio Kawashima could enlight this. Maybe it is a misconception.

Now, let's see the ansatz produced for Cyclo, Carbon and Acetylene:

We have moved too far. Now, let's see what the real solution states for every compound:

print(real_solution_cy
print(real_solution_c)
print(real_solution_ac)

CYCLOPROPENDYLENE
=== GROUND STATE ENERGY ==
 
* Electronic ground state energy (Hartree): -166.994973814548
  - computed part:      -1.204561706092
  - ActiveSpaceTransformer extracted energy part: -165.790412108456
~ Nuclear repulsion energy (Hartree): 52.358783483228
> Total ground state energy (Hartree): -114.63619033132
 
=== MEASURED OBSERVABLES ===
 
  0:  # Particles: 2.000 S: 0.000 S^2: 0.000 M: 0.000
 
=== DIPOLE MOMENTS ===
 
~ Nuclear dipole moment (a.u.): [65.28909274  42.68607856  21.98885319]
 
  0: 
  * Electronic dipole moment (a.u.): [65.63571639  44.03668228  22.29263495]
    - computed part:      [7.32297796  7.36272718  2.89479567]
    - ActiveSpaceTransformer extracted energy part: [58.31273842  36.67395511  19.39783928]
  > Dipole moment (a.u.): [-0.34662365  -1.35060372  -0.30378176]  Total: 1.42708154
                 (debye): [-0.88102935  -3.43289192  -0.77213614]  Total: 3.62727913

CARBON 
=== GROUND STATE ENERGY ===
 
* Electronic ground state energy (Hartree): -37.686791640847
  - computed part:      -5.327610251933
  - ActiveSpaceTransformer extracted energy part: -32.359181388913
~ Nuclear repulsion energy (Hartree): 0.0
> Total ground state energy (Hartree): -37.686791640847
 
=== MEASURED OBSERVABLES ===
 
  0:  # Particles: 4.000 S: 1.000 S^2: 2.000 M: 1.000
 
=== DIPOLE MOMENTS ===
 
~ Nuclear dipole moment (a.u.): [0.0  0.0  0.0]
 
  0: 
  * Electronic dipole moment (a.u.): [None  None  None]
    - computed part:      [None  None  None]
    - ActiveSpaceTransformer extracted energy part: [0.0  0.0  0.0]
  > Dipole moment (a.u.): [None  None  None]  Total: None
                 (debye): [None  None  None]  Total: None
 
ACETYLENE
=== GROUND STATE ENERGY ===
 
* Electronic ground state energy (Hartree): -101.545359183881
  - computed part:      -4.280823588566
  - ActiveSpaceTransformer extracted energy part: -97.264535595315
~ Nuclear repulsion energy (Hartree): 24.717023304133
> Total ground state energy (Hartree): -76.828335879747
 
=== MEASURED OBSERVABLES ===
 
  0:  # Particles: 4.000 S: 0.000 S^2: 0.000 M: 0.000
 
=== DIPOLE MOMENTS ===
 
~ Nuclear dipole moment (a.u.): [0.0  0.0  0.0]
 
  0: 
  * Electronic dipole moment (a.u.): [None  None  None]
    - computed part:      [None  None  None]
    - ActiveSpaceTransformer extracted energy part: [0.0  0.0  0.0]
  > Dipole moment (a.u.): [None  None  None]  Total: None
                 (debye): [None  None  None]  Total: None        

Let's calculate the Real Energy Reaction from the VQE in eV:

print("Reaction Real Energy VQE eV", round(27.2114*float(np.real(real_solution_cy.total_energies[0]-(real_solution_c.total_energies[0]+real_solution_ac.total_energies[0]))),10)," eV","\n")        

Reaction Real Energy VQE eV -3.2942885678 eV

Different from the value provided, but let's remember that observables are strongly affected with the Basis Set used, but it seems close to the classical energy value that we have calculated above. Let's move further and call the function custom_vqe() to get the optimal results to pass to the grader and get the score and the accuracy.

#Optimizers
max_iter_opt_cy = 200
max_iter_opt_c = 200
max_iter_opt_ac = 200


Energy_H_cy,_,jobs_cy = custom_vqe(estimator=Estimator(), ansatz=ansatz_cy, ops=ops_cy, maxiter_for_opt=max_iter_opt_cy, problem=problem_reduced_cy)

Energy_H_c,_,jobs_c = custom_vqe(estimator=Estimator(), ansatz=ansatz_c, ops=ops_c, maxiter_for_opt=max_iter_opt_c, problem=problem_reduced_c)

Energy_H_ac,_,jobs_ac = custom_vqe(estimator=Estimator(), ansatz=ansatz_ac, ops=ops_ac, maxiter_for_opt=max_iter_opt_ac, problem=problem_reduced_ac))        

Let's see if the calculated value converges to the real one:

react_vqe_ev = 27.2114*((Energy_H_cy[max_iter_opt_cy-1])-(Energy_H_c[max_iter_opt_cy-1] + Energy_H_ac[max_iter_opt_cy-1])

print("Reaction Energy VQE Estimator eV", round(np.real(react_vqe_ev),3)," eV","\n"))        

Reaction Energey VQE Estimator eV -3.446 eV

Let's plot VQE Energy for every compound calling the plot_graph() function:

plot_graph(Energy_H_cy, real_solution_cy.total_energies, "C3H2",color = "tab:green"
plot_graph(Energy_H_c, real_solution_c.total_energies, "C",color = "tab:blue")
plot_graph(Energy_H_ac, real_solution_ac.total_energies, "C2H2",color = "tab:red"))        

So, the custom VQE for every compound is converging. Now, let's take the optimal parameters from the jobs variable for Cyclo, Carbon and Acetylene:

initial_point_cy = [-3.34445129e-04, -1.77371722e+00,? 1.39224347e+00, -5.52903906e-04
??????? 2.25388734e-03,? 1.36890719e+00, -2.51063383e-01,? 1.40781414e+00]

initial_point_c = [-0.03238777, -0.01407733, -0.15615178,? 0.14009419,? 0.95310953,
??????? 0.22246143, -0.00216823,? 0.11353584,? 0.15452344, -0.53998153,
??????? 1.16457453,? 0.49118004]

initial_point_ac = [-0.00373967,? 0.58638724, -1.42120467,? 0.81087586,? 0.01428446,
?????? -0.00472041,? 0.00659752,? 2.54566948,? 1.73374402,? 1.39459684,
??????? 0.54686524, -0.21798029]        

We are moving closer to the end. Let's define the ultimate variables that the grader requires to compute the energy reaction and call it.

opt_cy = SPSA(maxiter=1
opt_c = SPSA(maxiter=1)
opt_ac = SPSA(maxiter=1) 
?
ansatz_list = [ansatz_cy, ansatz_c, ansatz_ac] # List of ansatz circuits

ops_list = [ops_cy, ops_c, ops_ac] # List of operators

problem_reduced_list = [problem_reduced_cy,? problem_reduced_c,? problem_reduced_ac] # List of ElectronicStrucutreProblem

initial_point_list = [initial_point_cy, initial_point_c, initial_point_ac]
optimizer_list= [opt_cy, opt_c, opt_ac] # List of optimal initial values

## Grade and submit your solution
from qc_grader.challenges.fall_2022 import grade_lab4_final

grade_lab4_final(ansatz_list, ops_list, problem_reduced_list, initial_point_list, optimizer_list))        

The output was:

Your computed reaction energy: -4.623832677235864 eV?

Your total score is 52125

Good result. Finally, we are going to define a Noise Reduction Strategy and pass it again to the grader to see what happens.

# ZNE Strategy
# Define Extrapolator
extrapolator = PolynomialExtrapolator(degree=2)
# Define Amplifier
noise_amplifier = LocalFoldingAmplifier(gates_to_fold=2)

# Define strategy
zne_strategy = ZNEStrategy(
??????? noise_factors=[1, 3, 5],
??????? noise_amplifier = noise_amplifier,
??????? extrapolator=extrapolator)        

Let's call again the grader:

grade_lab4_final(ansatz_list, ops_list, problem_reduced_list, initial_point_list, optimizer_list, zne_strategy)        

And the output is:

Again cool score, but no so accurate. Sure that the noisy strategy needs to be tuned.

Well, that's the end. I expect you find this article valuable. I may fall into some flaws due to I am not a Quantum Chemistry Expert, and it is a large and complex article, but as far as I can see, it was a great exercise of what Qiskit can do with Quantum Chemistry and the power of quantum computing.

Finally, it was a great experience to share knowledge with valuable professionals like Jacob Cybulski , Alain Chancé and Yan Wang. Also, mentors like Vishal Bajpe where always available to support and assist you when you were unable to move further. In the same way, assisting others participants like Mikhaela-Paige Early , Sai Ganesh Manda , Viktorija Bezganovi?, Jayakumar Vaithiyashankar , PhD and Louis Chen to pass the labs was fulfilling too.

If you find it suitable, let's engage in a debate in comments about the problem. Drop your ideas and suggestions.

References:

[1] Lab 1 - IBM Quantum Challenge Fall 2022

[2] Lab 2 - IBM Quantum Challenge Fall 2022

[3] Cyclopropendylene https://en.wikipedia.org/wiki/Cyclopropenylidene

[4] Lab 4 - IBM Quantum Challenge Fall 2022 (https://pubs.acs.org/doi/10.1021/jp020310z)