H₂ Molecule Homework Assignment
Quantum Software Development Journey: From Theory to Application with Classiq  Part 3

Similarly to what we have done in class, in this exercise we will implement the VQE on H2 molecule.

This time instead of using the builtin methods and functions (such as
Molecule
andMoleculeProblem
) to difne and solve the problem, you will be provided with a two qubits Hamiltonian.
Submission

Submit the completed Jupyter notebook and report via GitHub. Ensure all files are correctly named and organized.

Use the Typeform link provided in the submission folder to confirm your submission.
Additional Resources

The notebook from live session #3
Important Dates

Assignment Release: 22.5.2024

Submission Deadline: 3.6.2024 (7 A.M GMT+3)
Happy coding and good luck!
Part 1
Given the following Hamiltonian:
Complete the following code
from typing import List
from classiq import *
# Define Hamiltonian
HAMILTONIAN = QConstant(
"HAMILTONIAN",
List[PauliTerm],
[
PauliTerm(pauli=[Pauli.I, Pauli.I], coefficient=1.0523),
PauliTerm(pauli=[Pauli.I, Pauli.Z], coefficient=0.3979),
PauliTerm(pauli=[Pauli.Z, Pauli.I], coefficient=0.3979),
PauliTerm(pauli=[Pauli.Z, Pauli.Z], coefficient=0.0112),
PauliTerm(pauli=[Pauli.X, Pauli.X], coefficient=0.1809),
],
)
@qfunc
def main(q: Output[QArray[QBit]], angles: CArray[CReal, 3]) > None:
# Create an ansatz which allows each qubit to have arbitrary rotation
allocate(2, q)
U(angles[0], angles[1], angles[2], 0, q[0])
U(angles[0], angles[1], angles[2], 0, q[1])
@cfunc
def cmain() > None:
res = vqe(
HAMILTONIAN,
False,
[],
optimizer=Optimizer.COBYLA,
max_iteration=1000,
tolerance=0.001,
step_size=0,
skip_compute_variance=False,
alpha_cvar=1.0,
)
save({"result": res})
qmod = create_model(
main, classical_execution_function=cmain
) # complete the line, use classical_execution_function
qprog = synthesize(qmod)
show(qprog)
Opening: https://platform.classiq.io/circuit/5897752d79124b31bca898d6869c39de?version=0.42.1
execution = execute(qprog)
res = execution.result()
# execution.open_in_ide()
vqe_result = res[0].value
print(f"Optimal energy: {vqe_result.energy}")
print(f"Optimal parameters: {vqe_result.optimal_parameters}")
print(f"Eigenstate: {vqe_result.eigenstate}")
Optimal energy: 1.063205078125
Optimal parameters: {'angles_0': 6.181820728010731, 'angles_1': 4.1946685536401125, 'angles_2': 1.5078000440044168}
Eigenstate: {'01': (0.05412658773652741+0j), '10': (0.05412658773652741+0j), '00': (0.9970660083464885+0j)}
Does it similar to the optimal energy
we calculated in class? \
Does it similar to the total energy
we calculated in class?
**Does it similar to the `optimal energy` we calculated in class?**: No, the `optimal energy` is not similar to the one calculated in class (with the help of the `Molecule` and `MoleculeProblem` classes), since we are solving the Hamiltonian with two qubits separately, that is, the information between the qubits is not correlated in any way; therefore, there is no parameter setting that considers both qubits as a single system.
**Does it similar to the `total energy` we calculated in class?**: No, it is not similar to the `total energy` either because the `optimal energy` is not similar, so when calculating the `total energy` we cannot obtain a value similar to the one calculated in the class.\ (In the answers to _Part 2_, I calculate the `total energy`)
Optimal energy [Ha]  Total energy [Ha]  
From class (session 3)  1.857821819076969  1.137852824627989 
From this HW3 (Part 1)  1.064230957031249  ((na)) 
**Plot**: In the next cell, I show the plot of the energy optimization process of the VQE algorithm used. We can see that there is a convergence, but it reaches a plateau relatively quickly, where it remains "stuck" and is still unable to find a better value.
import base64
import io
from PIL import Image
image_bytes = base64.b64decode(vqe_result.convergence_graph_str)
# Usa io.BytesIO para convertir los bytes en un objeto de archivo
image_file = io.BytesIO(image_bytes)
# Abre la imagen con PIL
image = Image.open(image_file)
# Muestra la imagen en el notebook
# display(image) # uncomment to display
@qfunc
def main(q: Output[QArray[QBit]], angles: CArray[CReal, 3]) > None:
# Create an ansatz which allows each qubit to have arbitrary rotation
allocate(2, q)
U(angles[0], angles[1], angles[2], 0, q[0])
U(angles[0], angles[1], angles[2], 0, q[1])
CX(q[0], q[1])
@cfunc
def cmain() > None:
res = vqe(
HAMILTONIAN,
False,
[],
optimizer=Optimizer.COBYLA,
max_iteration=1000,
tolerance=0.001,
step_size=0,
skip_compute_variance=False,
alpha_cvar=1.0,
)
save({"result": res})
qmod = create_model(
main, classical_execution_function=cmain
) # complete the line, use classical_execution_function
qprog = synthesize(qmod)
# show(qprog)
execution = execute(qprog)
res = execution.result()
# execution.open_in_ide()
vqe_result = res[0].value
print(f"Optimal energy: {vqe_result.energy}")
print(f"Optimal parameters: {vqe_result.optimal_parameters}")
print(f"Eigenstate: {vqe_result.eigenstate}")
**Does it similar to the `optimal energy` we calculated in class?**: Yes, the `optimal energy` is similar to the one calculated in class (with the help of the `Molecule` and `MoleculeProblem` classes), this is because the calculation in this Part 2 considers the two correlated qubits (we added a CNOT gate between them) so that now the complete system can be considered formed by two qubits, and when calculating the values of the parameters (the angles) these are subject to working for both qubits at the same time.
**Does it similar to the `total energy` we calculated in class?**: Everything seems to indicate that this will be the case, but to confirm it we need to determine the value of the `nuclear repulsion energy`. The nuclear repulsion energy in a molecule refers to the energy associated with the repulsion between the nuclei of the atoms that compose it. In a molecule, the nuclei of atoms are positively charged and, therefore, repel each other due to the electrostatic force. To calculate this energy, we need to know the internuclear distance $R$. We could determine the distance between the two protons of the $H_2$ molecule from the Hamiltonian provided. To do this, we need to remember that the Hamiltonian coefficients are related to integrals of two electrons and other terms that depend on the internuclear distance $R$. But this is not simple to calculate. In fact, this is an inverse problem that is not solved directly from the Hamiltonian in its given form. Instead, we will consider that the hydrogen molecule is in its ground state, so we can take the average distance between the two nuclei, which would be equal to approximately $74 \ \text{pm} $ or $ 0.74 \ \text{Å} $ or $0.74 \times 10^{10} \ \text{m}$. To calculate the nuclear repulsion energy for the $H_2$ molecule, the formula for the electrostatic potential energy between two charges $q_1$ and $q_2$ can be used: $$E_{\text{repulsion}} = \frac{k \cdot q_1 \cdot q_2}{R}$$ where:  $ E_{\text{repulsion}} $ is the the nuclear repulsion energy  $ k $ is the Coulomb constant ($8.987 \times 10^9 \ \text{N·m}^2/\text{C}^2 $)  $ q_1 $ and $ q_2 $ are the charges of the nuclei (in the case of $H_2$, both are $ +1e $, where $ e $ is the charge of the proton, $ 1.602 \times 10^{19} \ \text{C} $)  $ R $ is the distance between the nuclei Therefore, \begin{align} E_{\text{repulsion}} &= \frac{(1.602 \times 10^{19} \ \text{C})^2}{4 \pi (8.854 \times 10^{12} \ \text{F/m}) \times 0.74 \times 10^{10} \ \text{m}}\\\\ &= \frac{2.566 \times 10^{38}}{1.112 \times 10^{10} \times 0.74 \times 10^{10}}\\\\ &= \frac{2.566 \times 10^{38}}{8.237 \times 10^{21}}\\\\ &= 3.115 \times 10^{18} \ \text{J} \end{align} We convert to Hartree: $$E_{\text{repulsion},(\text{Hartree})} = \frac{3.115 \times 10^{18} \ \text{J}}{4.3597447222071 \times 10^{18} \ \text{J/Hartree}} = 0.714 \ \text{Hartree}$$
Then, the `total energy` of the molecule is obtained by adding this nuclear repulsion energy to the electronic energy calculated from the given Hamiltonian: $$E_{\text{total}}= 1.8541486328124999 + 0.714 = 1.1401486328124999 \ \text{Hartree}$$
Therefore, we can conclude that the total energy is similar to that obtained in class with the help of the `Molecule` and `MoleculeProblem` classes, since a VQE algorithm for a twoqubit system works to determine the total energy of the hydrogen molecule $H_2$, as long as the qubits are correlated with each other.
Optimal energy [Ha]  Total energy [Ha]  
From class (session 3)  1.857821819076969  1.137852824627989 
From this HW3 (Part 2)  1.854148632812499  1.140148632812499 
**What can we learn about the provided form this result Hamiltonian?**: It is essential to consider the complete system as a whole; in this case, the system was composed of two qubits, and these must have a relationship between them to complete the representation of the Hamiltonian in a quantum circuit. We can find the value of the total energy of an appropriate way, the system must be internally correlated, in this case, the two qubits that make it up must have some quantum gate that correlates them.
**Plot**: In the next cell, I show the plot of the energy optimization process of the VQE algorithm used. We can see that there is a convergence, but it reaches a plateau relatively quickly, where it remains "stuck" and is still unable to find a better value.
import base64
import io
from PIL import Image
image_bytes = base64.b64decode(vqe_result.convergence_graph_str)
# Usa io.BytesIO para convertir los bytes en un objeto de archivo
image_file = io.BytesIO(image_bytes)
# Abre la imagen con PIL
image = Image.open(image_file)
# Muestra la imagen en el notebook
# display(image) # uncomment to display