# Quantum Phase Estimation (QPE) for Solving Molecular Energies¶

Quantum Phase Estimation (QPE) is a key algorithm in quantum computing, allowing one to estimate the phase (or eigenvalue) of an eigenvector of a unitary operation. The algorithm is designed so, that for a given Hamiltonian $H$, and an eigenvalue ${|\psi\rangle}$, the output that will be obtained is $\epsilon$ where

$U{|\psi\rangle} = e^{2\pi i\epsilon}{|\psi\rangle} , U = e^{2\pi iH}$ .

Therefore, by measuring the phase accumulated, the QPE algorithm allows calculating the energies relating to the chosen initial state. When using QPE for chemistry problems, it is common to search for the lowest energy of a given molecule. As the molecule can be written in the form of a Hamiltonian (Hermitian matrix which meaning is the energetic forces of the structure), usually symbolized by "H", one only need to insert the ground eigenvector in order to obtain the minimal energy value using QPE. However, obtaining the ground state is not a trivial problem. In order to overcome it, it is sufficient to use a state with big overlap with the ground state.

We define a state ${|v\rangle}$ which will be chosen as the algorithm's initial state. Let us define {$\psi_i$} to be the set of (unknown) eigenvalues of $H$. Generally, any vector can be rewritten as a superposition of any basis set, thus -

${|v\rangle} = \sum_i a_i{|\psi_i\rangle}$

and -

$U{|v\rangle} = \sum_i a_i e^{2\pi i\epsilon_i}{|\psi_i\rangle}$.

where ${\epsilon_i}$ are the eigenvalues of H, i.e. the span of energies relating to the molecule. Using execution with enough shots, one will obtain this set of $\theta_i$, i.e., a subset of the Hamiltonian's eigenvalues. As we are specifically interested in $\epsilon_0$, the ground state of H, it is important that there is a large overlap between ${\psi_0}$ and ${|v\rangle}$ so the probability to measure ${\epsilon_0}$ is high, i.e.

$P(\epsilon_0) = |\langle v|\psi_0\rangle|^2 > \epsilon $.

How large is $\epsilon$? After execution, we will obtain a set of ${E_i}$. If we have done 1000 shots of execution, and $P(\epsilon_0)>1$, it should be measured is about 10 times to get a good statistical variance.

A common choice for ${|v\rangle}$ (the initial state) is the Hartree-Fock (HF) state, which with good probability have a large overlap with the ground state. However, other guesses for the initial state are possibly good or even better fit, and choosing the right initial state is a sort of art and an active field of research.

For further reading about QPE we recommend [1].

#### What are the benefits of using QPE algorithm for finding a molecule's ground state?¶

The two most prominent methods to solve ground energy for molecules are quantum variational algorithm (VQE) and QPE. Those promise better scalability compared to classical counterparts as the molecules becomes more complex, with larger number of electrons manifesting as more degrees of freedom for the problem.

While so, the number of parameters in VQE is closely related to the number of electrons. This may create inherent difficulty achieving high-precision calculations through sampling statistical estimators, and paraphs even not converge for very large systems. While so, the number of parameters in QPE is a flexible value which is directly related to the resolution of the problem, but is not bonded with the number of electrons.

Furthermore, it is known that advanced quantum algorithms based on QPE can perform electronic structure calculations in sub-exponential time with accuracy that rivals exact diagonalization methods. This guarantee of simultaneously achieving high accuracy, efficiency, and generality is a feat that is believed to be impossible for classical algorithms. For those reasons, VQE is applicable in the near term (NISQ) era, while QPE is suited for fault-tolerant design.

#### During this tutorial, we will follow the QPE algorithms steps as following:¶

a. Define a molecule and convert it into a Hamiltonian H, manifested as a pauli list.

b. Prepare the Hamiltonian for QPE: including normalization and choosing the resolution.

c. Initializing the state for the HF state.

d. Preparing the quantum circuit.

e. Executing the circuit to find the related phases and analyzing the results to find the ground state.

f. Comparing the QPE results with exact solution for the molecule’s ground state.

## 0. Pre-requirments¶

The model is using several Classiq's libraries in addition to IBM's simulating tool.

```
import matplotlib.pyplot as plt
import numpy as np
import scipy
from numpy import linalg as LA
# for chemistry
from classiq import Model
from classiq.applications.chemistry import Molecule, MoleculeProblem
# for exponentiation
from classiq.builtin_functions import (
UCC,
Exponentiation,
HardwareEfficientAnsatz,
HartreeFock,
PhaseEstimation,
StatePreparation,
UnitaryGate,
)
from classiq.builtin_functions.exponentiation import (
ExponentiationConstraints,
ExponentiationOptimization,
PauliOperator,
)
from classiq.execution import IBMBackendPreferences
# for state preperation and phase estimation
from classiq.model import Constraints, Preferences
```

## 1. Create Your Molecule¶

```
# build your molecule
molecule_H2 = Molecule(atoms=[("H", (0.0, 0.0, 0)), ("H", (0.0, 0.0, 0.735))])
molecule_O2 = Molecule(atoms=[("O", (0.0, 0.0, 0)), ("O", (0.0, 0.0, 1.16))])
molecule_LiH = Molecule(atoms=[("H", (0.0, 0.0, 0.0)), ("Li", (0.0, 0.0, 1.596))])
molecule_H2O = Molecule(
atoms=[("O", (0.0, 0.0, 0.0)), ("H", (0, 0.586, 0.757)), ("H", (0, 0.586, -0.757))]
)
molecule_BeH2 = Molecule(
atoms=[("Be", (0.0, 0.0, 0.0)), ("H", (0, 0, 1.334)), ("H", (0, 0, -1.334))]
)
```

```
# define your molecule
molecule = molecule_H2
```

```
# define your molecule problem
gs_problem = MoleculeProblem(
molecule=molecule,
basis="sto3g",
mapping="jordan_wigner", #'bravyi_kitaev'
z2_symmetries=True,
freeze_core=True,
)
operator = gs_problem.generate_hamiltonian()
gs_problem = gs_problem.update_problem(operator.num_qubits)
print("Your Hamiltonian is", operator.show(), sep="\n")
```

```
pauli_ops = operator.pauli_list
pauli_operator = PauliOperator(pauli_list=pauli_ops)
N = pauli_operator.num_qubits
print("Your Pauli list is", pauli_ops, sep="\n")
```

## 2. Preparing the Molecule for QPE¶

### 2.1 Chose the Algorithm's Precision¶

For QPE algorithms, the precision is set by the number of qubits chosen $n$, such that the resolution is $1/{2^n}$. In case the matrix needs to be normlized, the resolution will be distorted. In case of normalization, that the span of results in the QPE is strached between the lowest and highest possible phase, thus the reslution will be mapped to $normalization-coefficient/{2^n} ~\sim 1/{((\lambda_{max}-\lambda_{min})*2^n)}$.

Here you can chose the precision in the parameter `n_qpe`

, or by setting your desired resolution. If you chose your desired resolution and set the parameter `get_recommended_n`

to be True, the number of qubits will be calculated for you accordingly.

```
# `n_qpe`- The precision of the QPE (the number of qubits to be used in the QPE. Creates resolution of 1/2^n_qpe
n_qpe = 6
```

```
# recommanded n_qpe:
get_recommanded_n = False
import math
desired_resolution = 0.02
def get_nqpe(pauli_operator, desired_resolution):
N = pauli_operator.num_qubits
A = 0
for a, b in pauli_operator.pauli_list:
A = A + abs(b)
nqpe = math.log2(2 * N * A / desired_resolution)
return math.ceil(nqpe)
if get_recommanded_n:
n_qpe = get_nqpe(pauli_operator, desired_resolution)
print("number of qubits for QPE is", n_qpe)
```

### 2.3 Normalize the Matrix¶

As QPE obtains a phase in the form $e^{2\pi i\theta}$, there is meaning only for $\theta \in [0,2\pi)$. Generally, our matrix M can have any eigenvalue, thus \theta can have any value. In order to fix this discrepancy, the values of the matrix are stretched to be rescaled. We assume $\theta \in [\lambda_{min}, \lambda_{max}]$ and use a normalization function in order to map those values into $[0, 1-1/{2^n}]$, where n is the number of qubits chosen for the QPE process in section 2.2.

We perform the normalization procedure as following:

a. We evaluate $\lambda_{min},\lambda_{max}$ (in the function ` normalization_params()`

below). In order to do so we use rough estimation of the absolute max value that can take place by adding together all the pauli coefficients and multiplying by the matrix's dimensions. That will yield us a value $\lambda$ (which is referred in the code as `normalization_coeff`

) and we now assume that the domain is $\theta \in [-\lambda, \lambda]$.
In general, one can build a more accurate assessment, which will decrease the span of solutions and thus achieve a better resolution.

b. We make sure only positive values are available by adding $\lambda*I^n$ to the pauli list. Now our evaluated span is $[0, 2*\lambda]$.

c. We normlize our matrix by multiplying all of the pauli coefficients by $(1-1/2^n)/(2*\lambda)$. Now the span of $\theta$ is $[0, 1-1/2^n]$, as required for proper QPE process.

The values of the matrix's eigenvalues should be now between $0$ to $1-(1/2^n)$. The QPE procedure will be performed on this new normalized matrix. After the phases are obtained, the original phases of the pre-normalized matrix will be gathered by performing opposite steps to this normalization procedure.

- Note that in case your matrix's eigenvalues are naturally between the values $0$ to $1-(1/2^n)$, you may not want to normalize it, as the normalization procedure may enlarge the span, thus lowering the resolution of the agorithm. In that case, you may skip those lines or change the value
`normalize`

to False.

```
# normalizing the operator
## we need to create a matrix such that its normalized version will have eigenvalues of [0,1/2^k] when k is the resolution of the QPE
normalize = True
def normalization_params(pauli_operator, N):
A = 0
for a, b in pauli_operator.pauli_list:
A = A + abs(b)
return N * A
def normalize_hamiltonian(pauli_operator, normalization_coeff, k):
new_pauli_operator = []
for a, b in pauli_operator.pauli_list:
if a == "I" * N:
new_pauli_operator.append(
(
a,
(b + normalization_coeff)
* (1 - 1 / (2**k))
/ (2 * normalization_coeff),
)
)
else:
new_pauli_operator.append(
(a, b * (1 - 1 / (2**k)) / (2 * normalization_coeff))
)
return new_pauli_operator
if normalize:
normalization_coeff = normalization_params(pauli_operator, N)
new_pauli_operator = normalize_hamiltonian(
pauli_operator, normalization_coeff, n_qpe
)
pauli_ops = new_pauli_operator
print(pauli_ops)
```

## 3. Creating the Quantum Circuit¶

We will now create a quantum circuit of the QPE algorithm using the Classiq platform. The user is able to fill in their constraints and perefences as desired.

```
constraints = Constraints()
preferences = Preferences(timeout_seconds=600)
model = Model(preferences=preferences, constraints=constraints)
```

### 3.1. Initializing the State¶

The state chosen for the initial state is HF, as it usually has good overlap with the initial state.

```
hf_params = HartreeFock(gs_problem=gs_problem)
output_dict = model.HartreeFock(params=hf_params)
sp_out = output_dict["OUT"]
```

### 3.1. Calling the Phase Estimation¶

Hereby we send the pauli list created for the molecule for creating phase estimation circuit.

```
from classiq.builtin_functions.qpe import (
ExponentiationScaling,
ExponentiationSpecification,
)
exp_max_depth = 2200
po = pauli_ops
exp_params = Exponentiation(
pauli_operator=PauliOperator(pauli_list=po), evolution_coefficient=-2 * np.pi
)
qpe_params = PhaseEstimation(
size=n_qpe,
unitary_params=exp_params,
exponentiation_specification=ExponentiationSpecification(
scaling=ExponentiationScaling(
max_depth=exp_max_depth, max_depth_scaling_factor=2**0.5
)
),
)
Qreg_name = {"IN": "IN", "OUT": "OUT"}
qpe_out = model.PhaseEstimation(params=qpe_params, in_wires={Qreg_name["IN"]: sp_out})
```

### 3.3. Defining the Registers¶

The following step is made in order to have a map of the relevant outputted qubits for analysis.

```
model.set_outputs(
{"phase_result": qpe_out["PHASE_ESTIMATION"], "out": qpe_out[Qreg_name["OUT"]]}
)
```

### 3.4. Synthesizing the Circuit¶

Hereby we syntesize the circuit and show it using the analyzer.

```
from classiq.execution import ExecutionPreferences
num_shots = 200000
model.execution_preferences = ExecutionPreferences(num_shots=num_shots)
model.sample()
qmod = model.get_model()
```

```
with open("qpe_for_molecules.qmod", "w") as f:
f.write(qmod)
```

```
from classiq import show, synthesize
qprog = synthesize(qmod)
show(qprog)
```

# 4. Measurment and Analysis¶

### 4.1. Circuit Execution¶

The circuit is now sent to execution on a default simulator

```
from classiq import execute
results = execute(qprog).result()
```

```
from classiq.execution import ExecutionDetails
results = results[0].value
```

### 4.2. Presenting the Result's Count¶

```
results.output_qubits_map
```

```
qpe_results = results.counts_of_output("phase_result")
print(qpe_results)
```

### 4.3. Chosing the Most Probable Solution¶

Hereby the user will choose the number of eigenvalues they wish to extract from the poll of results. The value `number_of_solutions`

will determine how many results out of `qpe_results`

will be analyzed. The algorithm `solution_search`

will save only the `number_of_solutions`

highest count results from the execution.

```
number_of_solutions = 1 # number of phases sought
```

```
def solution_search(result_dict, num):
def maxval(d):
v = list(d.values())
k = list(d.keys())
return k[v.index(max(v))]
sol_vec = []
for i in range(num):
A = maxval(result_dict)
sol_vec.append(A)
result_dict.pop(A)
return sol_vec
solution = solution_search(qpe_results, number_of_solutions)
print("Your qubit solution is", solution, sep="\n")
```

### 4.4. Translating into Energy value (Phase)¶

Here, the value in the vector `results`

are translated from a binary number into a full solution for the eigenvalues. ### 3.3. Translating into Eigenvalue (Phase)

At first - the binary number is translated into a decimal value.

```
# translate binary into decimal
def bin_to_dec(q_num, solution):
dec_sol = 0
for i in range(q_num):
dec_sol = dec_sol + int(solution[i]) * (2**i)
return dec_sol / (2**q_num)
dec_sol_vec = []
for i in range(number_of_solutions):
dec_solution = bin_to_dec(n_qpe, solution[i])
dec_sol_vec.append(dec_solution)
print("Your decimal solutions are", dec_sol_vec, sep="\n")
```

Secondy, the decimal value is mapped back into the original values, i.e. renormalized into it's original span.

```
# renormalize into the "real" solution -
if normalize:
solution = [
(
(value * 2 * normalization_coeff / (1 - (1 / 2**n_qpe)))
- normalization_coeff
)
for value in dec_sol_vec
]
else:
solution = dec_sol_vec
```

The final results of the Energy (matrix's eigenvlues) is therefore:

```
print(solution)
```

And the results including the error contributed from the resolution (number of qubits participating in the QPE) are:

```
if normalize:
energy_resolution = (
(1 / (2**n_qpe)) * 2 * normalization_coeff / (1 - (1 / 2**n_qpe))
)
else:
energy_resolution = 1 / (2**n_qpe)
print("the resolution of result is", energy_resolution)
for sol in solution:
print(
"the eneregy solution is beteen",
sol - energy_resolution,
"and",
sol + energy_resolution,
)
###if zero or exceed the normalization range need to add conditions
```

### 4.4. Find Exact Solution for Comparison¶

```
mat = pauli_operator.to_matrix()
w, v = np.linalg.eig(mat)
print("exact result:", np.real(min(w)))
```

### 4.4. Find the Solution's Histogram¶

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

```
import matplotlib.pyplot as plt
import numpy as np
energy_vec = []
energy_prob = []
for string, counts in qpe_results.items():
temp = bin_to_dec(n_qpe, string)
if normalize:
temp2 = (
temp * 2 * normalization_coeff / (1 - (1 / 2**n_qpe))
) - normalization_coeff
else:
temp2 = temp
energy_vec.append(temp2)
energy_prob.append(counts / num_shots)
```

```
plt.plot(energy_vec, energy_prob, ".")
plt.show()
```