<jupyter_start><jupyter_text>Variational Quantum Eigensolver for Fermions: tutorialIn this notebook, we want to find a good approximation of the minimum energy ($E_0$) of a fermionic Hamiltonian $H$ of the form :$$ H=\sum_{pq}h_{pq}c_p^\dagger c_q + \frac{1}{2} \sum_{pqrs} h_{pqrs}c_p^\dagger c_q^\dagger c_r c_s $$To find this minimum, we invoke the Rayleigh-Ritz variational principle that stipulates that:$$ \frac{\langle\Psi|H|\Psi\rangle}{\langle\Psi|\Psi\rangle} \geq E_0 \; \; \forall \Psi$$The idea is to construct a family of variational wavefunctions $|\Psi(\vec{\theta})\rangle$ and find the parameter $\vec{\theta}^*$ that minimizes $\frac{\langle\Psi(\vec{\theta})|H|\Psi(\vec{\theta})\rangle}{\langle\Psi(\vec{\theta})|\Psi(\vec{\theta})\rangle} $.  To construct such states, we use a quantum circuit that generates $\Psi(\vec{\theta})$ given a list of parameters $\vec{\theta}$. We then measure the mean value of $H$, $E(\vec{\theta})$, which we feed to a classical optimizer that will find a new set of parameters, and so on until a desired accuracy is achieved.  Defining the HamiltonianHere, we define the fermionic Hamiltonian. For simplicity, we will focus on a very simple example: a **single-site Hubbard model**, defined by the Hamiltonian:$$ H = U n_\uparrow n_\downarrow - \mu \sum_{\sigma=\uparrow,\downarrow} n_\sigma $$with $n_\sigma \equiv c_\sigma^\dagger c_\sigma$.Such a Hamiltonian is already part of the QLM's ``qat.dqs`` module via the ``make_hubbard_model`` function below. This function creates the hopping and interaction integrals $h_{pq}$ and $h_{pqrs}$ that define our Hamiltonian.<jupyter_code>import numpy as np
import scipy.optimize

from qat.dqs.hamiltonians import make_hubbard_model

U = 2.0
nqbit = 2
t_mat = np.zeros((1, 1))
hamilt = make_hubbard_model(t_mat, U, mu=U/2)

print("H = ", hamilt)<jupyter_output><empty_output><jupyter_text>The first step is to transform this fermionic Hamiltonian to a Hamiltonian that acts on qubits, i.e two-level systems. This transformation is not unique. Common transformations are the Jordan-Wigner and the Bravyi-Kitaev transformation. In the following cell, we use the Jordan-Wigner transform:<jupyter_code>from qat.dqs.transforms import transform_to_jw_basis
hamilt_spin = transform_to_jw_basis(hamilt)
print("H (spin) = ", hamilt_spin)<jupyter_output><empty_output><jupyter_text>Ansatz and corresponding quantum circuitThe second step is to find a circuit with some parameters to generate our variational wavefunction $\Psi(\vec\theta)$. This ansatz has to be powerful enough to be able to come close enough to the ground state, but with as few as possible gates (to reduce the circuit depth) and parameters (to ease the task of the classical optimizer).In our example, we take a list with only two parameters and we create a function that returns a ``QRoutine``. Note that we insert a CNOT because we know that our Hamiltonian contains terms entangling two orbitals in our fermionic system.We encapsulate the QRoutine inside a Program that is itself converted into a parametric circuit:<jupyter_code>from qat.lang.AQASM import H, RX, RY, CNOT, QRoutine, Program

def simple_circuit_with_two_parameters(theta):
    """Take a parameter theta and return the corresponding circuit"""
    Qrout = QRoutine()
    Qrout.apply(H, 0)
    Qrout.apply(RY(theta[0]), 0)
    Qrout.apply(CNOT, 0, 1)
    Qrout.apply(RX(theta[1]), 1)
    return Qrout

prog = Program()
reg = prog.qalloc(hamilt_spin.nbqbits)
prog.apply(simple_circuit_with_two_parameters([prog.new_var(float, '\\theta_%s'%i)
                                               for i in range(hamilt_spin.nbqbits)]), reg)
circ = prog.to_circ()

%qatdisplay circ<jupyter_output><empty_output><jupyter_text>Job preparationWhat we now want to achieve is the minimization of the expectation value $\langle \Psi(\vec{\theta})| H |\Psi(\vec{\theta})\rangle$ of the Hamiltonian $H$ at hand on the variational wavefunction $|\Psi(\vec{\theta})\rangle$ over the two-dimensional space of variational parameters $\vec{\theta}$. Such a computational job is described by a ``Job`` object that essentially contains information about- the variational circuit: ``circ`` below- the observable (Hamiltonian) to be minimized: ``hamilt_spin``The ``nbshots`` variable indicates the number of repetitions of the state preparation and measurement that are devoted to the computation of the expectation value<jupyter_code>job = circ.to_job(job_type="OBS", observable=hamilt_spin, nbshots=0)
theta0 = [0.5, 1.23]<jupyter_output><empty_output><jupyter_text>This job is a parametric job because the parameters $\vec{\theta}$ of the circuit are not fixed. To handle such jobs, the QLM supplies special plugins whose purpose is to minimize the value of a job's observable over the job's parametric circuit.The minimization can be carried out using various classical optimizers.In the cell below, we use Scipy's Nelder-Mead optimizer as a backend for this minimization. The corresponding dedicated QLM Plugin is ``ScipyMinimizePlugin``:<jupyter_code>from qat.plugins import ScipyMinimizePlugin
from qat.qpus import LinAlg

optimizer_scipy = ScipyMinimizePlugin(method="Nelder-Mead",
                                      x0=theta0,
                                      tol=1e-3,
                                      options={"maxiter": 600})
linalg_qpu = LinAlg()
qpu = optimizer_scipy | linalg_qpu
result = qpu.submit(job)

print("Minimum VQE energy =", result.value)<jupyter_output><empty_output><jupyter_text>The ``result`` variable contains additional information about the variational procedure, such as the successive values of the variational energy over the course of the optimization, or the value of the optimal parameters $\vec{\theta}^*$:<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(eval(result.meta_data['optimization_trace']))
plt.xlabel("Steps")
plt.ylabel("Energy");

print("Theta* = ", result.meta_data["parameters"])<jupyter_output><empty_output><jupyter_text>Simultaneous Perturbation Stochastic Approximation optimizerWe now switch to another classical optimizer, SPSA:<jupyter_code>from qat.plugins import SPSAMinimizePlugin

optimizer_spsa = SPSAMinimizePlugin(theta0)

qpu = optimizer_spsa | linalg_qpu
result = qpu.submit(job)

print("Minimum energy =", result.value)
plt.plot(eval(result.meta_data['optimization_trace']))
plt.xlabel("Steps");
plt.ylabel("Energy");<jupyter_output><empty_output><jupyter_text>Particle-Swarm OptimizerSimilarly, we can use a Particle-Swarm Optimizer, PSO. We can specify the number of swarm particles used for the optimization. In the following example, we take 10.<jupyter_code>from qat.plugins import PSOMinimizePlugin

optimizer_pso = PSOMinimizePlugin(theta0, max_iter=50, swarm_number=10)
qpu = optimizer_pso | linalg_qpu
result = qpu.submit(job)

print("Minimum energy =", result.value)

plt.plot(eval(result.meta_data['optimization_trace']))
plt.xlabel("Steps");
plt.ylabel("Energy");<jupyter_output><empty_output><jupyter_text>Exact resultIt may be interesting to compare this result with an exact result. This can be done, in this simple case, by a full exact diagonalisation:<jupyter_code>eigvals, eigvecs = np.linalg.eigh(hamilt.get_matrix())
print("Exact ground-state energy =", min(eigvals))<jupyter_output><empty_output><jupyter_text>VQE on a noisy computerWe can now define a noisy hardware model and see how VQE behaves.<jupyter_code>from qat.hardware import HardwareModel, DefaultGatesSpecification
from qat.quops import ParametricAmplitudeDamping, ParametricPureDephasing
from qat.qpus import NoisyQProc

gate_durations = {"Z": 60 ,"H": 60,  "X":60, "Y":60,
                  "RY": lambda angle : 60, "RX": lambda angle : 60,
                  "CNOT":150}
T1 = 1000
T2 = 1000

amp_damping = ParametricAmplitudeDamping(T_1=T1)
pure_dephasing= ParametricPureDephasing(T_phi=1/(1/T2 - 1/(2*T1)))
hw_model = HardwareModel(DefaultGatesSpecification(gate_durations),
                         idle_noise=[amp_damping, pure_dephasing])
noisy_qproc = NoisyQProc(hardware_model=hw_model, sim_method="deterministic-vectorized")<jupyter_output><empty_output><jupyter_text>It suffices to replace an ideal qpu with a noisy one to obtain the result in the presence of noise with any of the two above-mentioned methods. One can also select various `sim_method` arguments (see commented-out line for an example with "stochastic" method).<jupyter_code>optimizer_scipy = ScipyMinimizePlugin(method="COBYLA",
                                      x0=theta0,
                                      tol=1e-3,
                                      options={"maxiter": 600})
qpu = optimizer_scipy | noisy_qproc
result = qpu.submit(job)

print("Minimum energy =", result.value)
plt.plot(eval(result.meta_data['optimization_trace']))
plt.xlabel("Steps");
plt.ylabel("Energy");<jupyter_output><empty_output><jupyter_text>Influence of the number of shotsSo far we have assumed each circuit observable was measured exactly, i.e with an infinite number of projective measurements. We can change this by setting the parameter ``nbshots`` to a finite number (it defaults to 0, i.e infinity). Choosing an insufficient number of shots leads to an inaccurate result due to shot noise.<jupyter_code>optimizer_scipy = ScipyMinimizePlugin(method="COBYLA",
                                      x0=theta0,
                                      tol=1e-3,
                                      options={"maxiter": 600})
qpu = optimizer_scipy | linalg_qpu
nbshots_list = [50, 1500, 0]
for nbshots in nbshots_list:
    job = prog.to_circ().to_job(job_type="OBS",
                                observable=hamilt_spin,
                                nbshots=nbshots)
    result = qpu.submit(job)
    print("==== nb shots =", nbshots)
    print("Minimum energy =", result.value)
    plt.plot(eval(result.meta_data['optimization_trace']),
             label=r"$n_\mathrm{shots} = %s$"%nbshots)
plt.xlabel("Steps");
plt.ylabel("Energy");
plt.legend();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10132377824226137, 'fraction_numerical': 0.009356026674629243, 'mean_word_length': 4.46383904295813, 'pattern_counts': {'":': 10, '<': 49, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '430927', 'license': 'Apache-2.0', 'n_tokens_mistral': 3187, 'n_tokens_neox': 2974, 'n_words': 965, 'path': 'misc/notebooks/tutorials/digital_quantum_simulation/vqe_fermions.ipynb', 'repo_name': 'wmkirby1/CS-VQE', 'size': '10047'}