Cirq: Variational Quantum Eigensolver for \(H_2\)

Author: Veiko Palge (University of Tartu)

In this notebook, we are going to implement the variational quantum eigensolver (VQE) algorithm to compute the ground energy of molecular hydrogen \(H_2\). This is one of the smallest systems one can imagine to start experimenting with in VQE. It has also been implemented on an actual quantum computer as reported in O’Malley etal 2015, Scalable Quantum Simulation of Molecular Energies, arXiv:1512.06860v2. The following figure given in the article summarizes the entire setup:

We are going to use Google’s Cirq library and implement the VQE algorithm for \(H_2\) from scratch using the circuit in the figure and the data given in the paper.

Synopsis

Before we begin, let us summarize the task ahead of us, and what exactly we are going to do. Remember that VQE starts with classical preparations, this involves classical computational chemistry computations. As a result of these preparations we have the qubit Hamiltonian and the ansatz circuit that represents the parameterized unitary. We assume that all this work has been completed: the result is the ansatz circuit and the Hamiltonian shown in the figure above.

Let us summarize the computation shown in the figure.

The algorithm consists of two parts. Firstly, there is a computation carried out on a quantum computer. This is divided into three parts in the figure, (1) Prepare initial state, (2) Apply Parameterized Ansatz, (3) Measure Expectation values. The entire quantum computation is performed using the circuit shown at the bottom part, ‘Software’. Note the circuit contains only two qubits.

Secondly, there is the classical part, comprising (1) Calculate Energy and (2) Classical Optimizer Suggests New Parameters \(\theta\). The qubit Hamiltonian of \(H_2\) is also shown on the right.

Note that the top part of the figure also shows a rough outline of physical pulses with timings that are needed to run the circuit.

We will implement both the classical and quantum parts in Python using Google’s Cirq. We will implement the algorithm from scratch without using the more advanced Cirq libraries. We will only rely circuits and measurements. In principle, we will show how to implement the algorithm in any language that contains circuits and measurements.

The rationale of showing this simple, low level implementation is to give a taste of what is going on behind the scenes of the more advanced libraries.

Implementation

import cirq
import numpy as np
import sympy
from cirq.circuits import InsertStrategy

Let us start by creating a routine for the circuit in the figure.

##########################################################
#
#  Routines
#
##########################################################
def create_ansatz_circuit(qbits, parameter):
    """
    Create H2 ansatz circuit given in O'Malley etal 2015.

    Args:
        qbits          - qubits
        parameter      - parameter for parameterized circuit

    Returns:
        ansatz_circuit - parameterized circuit
    """
    ansatz_circuit = cirq.Circuit()

    # H2 circuit
    # 1 - create HF reference state
    ansatz_circuit.append(cirq.rx(np.pi)(qbits[0]), strategy=InsertStrategy.NEW)

    # 2 - ansatz: U(theta)
    ansatz_circuit.append(cirq.rx(-np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    # 3
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 4 - parameterized Z-rotation
    ansatz_circuit.append(cirq.rz(parameter)(qbits[0]), strategy=InsertStrategy.NEW)

    # 5
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 6
    ansatz_circuit.append(cirq.rx(np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(-np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    return ansatz_circuit


##########################################################
#
#  Main program
#
##########################################################

# variables
theta = sympy.Symbol("θ")
n_qubits = 2

qubits = cirq.LineQubit.range(n_qubits)
ansatz_circuit = create_ansatz_circuit(qubits, theta)

print(ansatz_circuit)

0: ───Rx(π)───Rx(-0.5π)──────────────X───Rz(θ)───X───Rx(0.5π)───────────────
                                     │           │
1: ───────────────────────Ry(0.5π)───@───────────@──────────────Ry(-0.5π)───

The new feature here is a parameter that is inserted in the rz() rotation in L39. It means every time when we update the parameter, we do not need to create a new circuit from scratch. This saves a lot of time in general. VQE circuits can be quite large. Imagine a circuit containing \(10\)k gates and \(300\) parameters. Creating fixed circuits each time one parameter needs to be updated will take much longer than updating a single value in the circuit. For this reason, when implementing VQE, it’s always best to use parameterized circuits if they are available in the programming language.

At this point we have accomplished the first two quantum components, (1) Prepare initial state, (2) Apply Parameterized Ansatz.

Measure Expectation values

Next, let us focus on (3) Measure Expectation values. The qubit Hamiltonian we need to measure is as follows

\[ H = g_0 \mathbb{1} + g_1 Z_0 + g_2 Z_1 + g_3 Z_0 Z_1 + g_4 Y_0 Y_1 + g_5 X_0 X_1 \]

where \(g_i\) are (real valued) coefficients that have been computed classically. They are given in the Table 1 in the appendix of O’Malley etal 2015. The values of \(g_i\) are functions of hydrogen-hydrogen bond length \(R\). We will take the coefficients at bond length \(R = 0.75\) (in units \(10^{-10}\) m) where the energy is lowest (the actual bond length where the energy is lowest is \(0.74\), which is close).

There is a problem though. We can perform measurements in the \(Z\) basis using the measure() routine. However, we cannot directly measure the Pauli \(X\) or \(Y\) operators of the terms \(X_0 X_1\) and \(Y_0 Y_1\). In order to measure the latter, we have to use the following trick:

• to measure \(X\), perform basis transformation from the \(X\) and \(Z\)

• then measure the qubits as usual (in the \(Z\) basis) with measure().

The same applies to \(Y\) measurements, in this case we need to perform basis transformation from \(Y\) to \(Z\).

In order to compute the average value \(\langle H \rangle\) of the (total) Hamiltonian we will use the fact that the total average is equal to the sum of averages of its terms,

\[ \langle H\rangle = g_0 \langle \mathbb{1}\rangle + g_1 \langle Z_0\rangle + g_2 \langle Z_1\rangle + g_3 \langle Z_0 Z_1\rangle + g_4 \langle Y_0 Y_1\rangle + g_5 \langle X_0 X_1\rangle. \]

So we need to compute the expectation value of each term, then multiply each term with its coefficient \(g_i\), and finally add averages of all terms. Note that \(\langle\mathbb{1}\rangle = 1\), so we can just add the coefficient \(g_0\).

Pseudo code wise, this means:

• Compute avg of H:

  • Loop over terms, compute avg of each term

Computing the average of a term like \(X_0 X_1\) means that we need to add the respective measurement mini-circuit at the end of the ansatz. This is how the part ‘(3) Measure Expectation values’ is realized in the circuit above in the figure, see the final light grey box.

In terms of implementation, it means that instead of looping over terms we can loop over measurement mini-circuits. And, in order to be more efficient, we will create the mini-circuits prior to running the algorithm.

So let us add a routine that takes in a term (e.g. \(X_0 X_1\)) and creates the corresponding mini-circuit, which we can add at the end of the main ansatz circuit.

##########################################################
#
#  Routines
#
##########################################################
def create_ansatz_circuit(qbits, parameter):
    """
    Create H2 ansatz circuit given in O'Malley etal 2015.

    Args:
        qbits          - qubits
        parameter      - parameter for parameterized circuit

    Returns:
        ansatz_circuit - parameterized circuit
    """
    ansatz_circuit = cirq.Circuit()

    # H2 circuit
    # 1 - create HF reference state
    ansatz_circuit.append(cirq.rx(np.pi)(qbits[0]), strategy=InsertStrategy.NEW)

    # 2 - ansatz: U(theta)
    ansatz_circuit.append(cirq.rx(-np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    # 3
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 4 - parameterized Z-rotation
    ansatz_circuit.append(cirq.rz(parameter)(qbits[0]), strategy=InsertStrategy.NEW)

    # 5
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 6
    ansatz_circuit.append(cirq.rx(np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(-np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    return ansatz_circuit


def paulis_to_measure_circ(qbits, hamlist):
    """
    Generate measure circuits from hamiltonian pauli list
    for measurements in different bases.

    Args:
        qbits     - qubits
        hamlist   - total hamiltonian, i.e. pauli strings

    Returns:
        circuits - list of circuits that can be used to average over
    """
    circuits = []

    for elem in hamlist:
        minicirc = cirq.Circuit()

        for qubitno, qubitop in elem:
            if qubitop == "Z":
                pass
            elif qubitop == "X":
                minicirc.append(cirq.H(qbits[qubitno]))
            elif qubitop == "Y":
                minicirc.append(cirq.S(qbits[qubitno]) ** -1)
                minicirc.append(cirq.H(qbits[qubitno]))
            else:
                assert False, "Error: INVALID qubit operation"

        minicirc.append(cirq.measure(qbits, key="m"))
        circuits.append(minicirc)

    return circuits


##########################################################
#
#  Main program
#
##########################################################

# variables
theta = sympy.Symbol("θ")
n_qubits = 2
hamiltonian = (
    ((0, "Z"),),
    ((1, "Z"),),
    ((0, "Z"), (1, "Z")),
    ((0, "X"), (1, "X")),
    ((0, "Y"), (1, "Y")),
)

qubits = cirq.LineQubit.range(n_qubits)
ansatz_circuit = create_ansatz_circuit(qubits, theta)

ham_mini_circuits = paulis_to_measure_circ(qubits, hamiltonian)

# print result
[print(circ) for circ in ham_mini_circuits]

0: ───M('m')───
      │
1: ───M────────
0: ───M('m')───
      │
1: ───M────────
0: ───M('m')───
      │
1: ───M────────
0: ───H───M('m')───
          │
1: ───H───M────────
0: ───S^-1───H───M('m')───
                 │
1: ───S^-1───H───M────────

[None, None, None, None, None]

The new routine paulis_to_measure_circ() takes in the total qubit Hamiltonian and outputs circuits which correspond to the terms of the Hamiltonian. We also needed to represent the Hamiltonian, to this end we use a tuple of tuples in L95. The smallest element has the form (qubit_number, Pauli_matrix), where qubit_number is the index of the qubit for the Pauli observable.

Let’s summarize the computation of the average of the Hamiltonian. We need to compute the averages of terms. We do this by looping over terms, which are represented as mini-circuits. For each mini-circuit, we will take the ansatz, add to it the mini-circuit, and run the resulting circuit, collect measurement statistics and compute the average.

Scan the parameter interval

We will not bother with classical optimization in this small implementation. The circuit has only one parameter—the angle \(\theta\) of the \(R_z\) rotation—which ranges between \([-\pi, \pi]\). Instead of searching for the minimum value of \(\langle H\rangle\), we will scan this interval and record the minimum energy.

Let us add a loop to do this.

##########################################################
#
#  Routines
#
##########################################################
def create_ansatz_circuit(qbits, parameter):
    """
    Create H2 ansatz circuit given in O'Malley etal 2015.

    Args:
        qbits          - qubits
        parameter      - parameter for parameterized circuit

    Returns:
        ansatz_circuit - parameterized circuit
    """
    ansatz_circuit = cirq.Circuit()

    # H2 circuit
    # 1 - create HF reference state
    ansatz_circuit.append(cirq.rx(np.pi)(qbits[0]), strategy=InsertStrategy.NEW)

    # 2 - ansatz: U(theta)
    ansatz_circuit.append(cirq.rx(-np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    # 3
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 4 - parameterized Z-rotation
    ansatz_circuit.append(cirq.rz(parameter)(qbits[0]), strategy=InsertStrategy.NEW)

    # 5
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 6
    ansatz_circuit.append(cirq.rx(np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(-np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    return ansatz_circuit


def paulis_to_measure_circ(qbits, hamlist):
    """
    Generate measure circuits from hamiltonian pauli list
    for measurements in different bases.

    Args:
        qbits     - qubits
        hamlist   - total hamiltonian, i.e. pauli strings

    Returns:
        circuits - list of circuits that can be used to average over
    """
    circuits = []

    for elem in hamlist:
        minicirc = cirq.Circuit()

        for qubitno, qubitop in elem:
            if qubitop == "Z":
                pass
            elif qubitop == "X":
                minicirc.append(cirq.H(qbits[qubitno]))
            elif qubitop == "Y":
                minicirc.append(cirq.S(qbits[qubitno]) ** -1)
                minicirc.append(cirq.H(qbits[qubitno]))
            else:
                assert False, "Error: INVALID qubit operation"

        minicirc.append(cirq.measure(qbits, key="m"))
        circuits.append(minicirc)

    return circuits


def get_ham_avg(params, ansatz, coeffs, ham_mini_circuits, obs, shots):
    """
    Compute average of the total qubit Hamiltonian.

    Args:
        params           - parameter value for the ansatz circuit
        ansatz           - ansatz circuit
        ham_mini_ciruits - pre-built circuits to measure Hamiltonian terms
        obs              - data for observables
        shots            - number of shots

    Returns:
        avg   - average of Hamiltonian
    """
    avg: float = 0.0

    # TODO

    return avg


##########################################################
#
#  Main program
#
##########################################################

# variables
theta = sympy.Symbol("θ")
n_qubits = 2
hamiltonian = (
    ((0, "Z"),),
    ((1, "Z"),),
    ((0, "Z"), (1, "Z")),
    ((0, "X"), (1, "X")),
    ((0, "Y"), (1, "Y")),
)
shots: int = 1000
no_of_points: int = 40
param_range = np.arange(-np.pi, np.pi, np.pi / no_of_points)
# g_i coeffs for R = 0.75
g_coeff = [0.3435, -0.4347, 0.5716, 0.091, 0.091]
g_coeff_id = 0.2252
# graph
graph = []
#
obs = []

qubits = cirq.LineQubit.range(n_qubits)
ansatz_circuit = create_ansatz_circuit(qubits, theta)

ham_mini_circuits = paulis_to_measure_circ(qubits, hamiltonian)

# record min val of hamiltonian
h_min: float = 100.0

# scan theta domain [-pi, pi]
for x in param_range:
    h_avg = get_ham_avg(
        {theta: x}, ansatz_circuit, g_coeff, ham_mini_circuits, obs, shots
    )
    h_avg += g_coeff_id  # add coeff of id term
    graph.append(h_avg)  # collect data for plotting
    if h_avg < h_min:  # record min value
        h_min = h_avg

# result
print("\nH_min = {}\n".format(h_min))

H_min = 0.2252

The loop in lines L141-147 walks through the interval \([-\pi, \pi]\) and finds the minimum energy. Of course, the answer is wrong—the value \(0.2252\) is just the value of the constant term—since we have not implemented the actual computation of the average.

Mapping outcome states to values

To actually compute the average of a Hamiltonian term, we need to know how the outcomes of the measurement correspond to the actual real valued outcomes \(+1, -1\) for a particular observable like \(Z_0\) or \(Z_0 Z_1\) we measure. Without discussing this at length here, we are going to introduce dictionaries that record the mapping for each observable. Of course, there are other ways to implement this.

Let us fill in the generic routine for computing the average.

##########################################################
#
#  Routines
#
##########################################################
def create_ansatz_circuit(qbits, parameter):
    """
    Create H2 ansatz circuit given in O'Malley etal 2015.

    Args:
        qbits          - qubits
        parameter      - parameter for parameterized circuit

    Returns:
        ansatz_circuit - parameterized circuit
    """
    ansatz_circuit = cirq.Circuit()

    # H2 circuit
    # 1 - create HF reference state
    ansatz_circuit.append(cirq.rx(np.pi)(qbits[0]), strategy=InsertStrategy.NEW)

    # 2 - ansatz: U(theta)
    ansatz_circuit.append(cirq.rx(-np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    # 3
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 4 - parameterized Z-rotation
    ansatz_circuit.append(cirq.rz(parameter)(qbits[0]), strategy=InsertStrategy.NEW)

    # 5
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 6
    ansatz_circuit.append(cirq.rx(np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(-np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    return ansatz_circuit


def paulis_to_measure_circ(qbits, hamlist):
    """
    Generate measure circuits from hamiltonian pauli list
    for measurements in different bases.

    Args:
        qbits     - qubits
        hamlist   - total hamiltonian, i.e. pauli strings

    Returns:
        circuits - list of circuits that can be used to average over
    """
    circuits = []

    for elem in hamlist:
        minicirc = cirq.Circuit()

        for qubitno, qubitop in elem:
            if qubitop == "Z":
                pass
            elif qubitop == "X":
                minicirc.append(cirq.H(qbits[qubitno]))
            elif qubitop == "Y":
                minicirc.append(cirq.S(qbits[qubitno]) ** -1)
                minicirc.append(cirq.H(qbits[qubitno]))
            else:
                assert False, "Error: INVALID qubit operation"

        minicirc.append(cirq.measure(qbits, key="m"))
        circuits.append(minicirc)

    return circuits


def get_paulistr_avg(trials, obs, shots):
    """
    Calculate average of a Pauli string (e.g. XXY).

    Args:
        trials      - result of simulation
        obs         - data for calculating probabilities
        shots       - number of repetitions

    Returns:
        avg         - the average value
    """
    avg: float = 0.0

    # TODO

    return avg


def get_ham_avg(params, ansatz, coeffs, ham_mini_circuits, obs, shots):
    """
    Compute average of the total qubit Hamiltonian.

    Args:
        params           - parameter value for the ansatz circuit
        ansatz           - ansatz circuit
        ham_mini_ciruits - pre-built circuits to measure Hamiltonian terms
        obs              - data for observables
        shots            - number of shots

    Returns:
        avg   - average of Hamiltonian
    """
    avg: float = 0.0
    sim = cirq.Simulator()

    for coeff, measure_circuit, obsval in zip(coeffs, ham_mini_circuits, obs):
        total_circuit = ansatz + measure_circuit
        cirq_trials = sim.run_sweep(
            program=total_circuit, params=params, repetitions=shots
        )
        paulistring_avg = get_paulistr_avg(cirq_trials, obsval, shots)
        avg += coeff * paulistring_avg

    return avg


##########################################################
#
#  Main program
#
##########################################################

# variables
theta = sympy.Symbol("θ")
n_qubits = 2
hamiltonian = (
    ((0, "Z"),),
    ((1, "Z"),),
    ((0, "Z"), (1, "Z")),
    ((0, "X"), (1, "X")),
    ((0, "Y"), (1, "Y")),
)
shots: int = 1000
no_of_points: int = 40
param_range = np.arange(-np.pi, np.pi, np.pi / no_of_points)
# g_i coeffs for R = 0.75
g_coeff = [0.3435, -0.4347, 0.5716, 0.091, 0.091]
g_coeff_id = 0.2252
# graph
graph = []
#
obs = []
# map outcome states to outcome values
obs = [
    {"00": 1, "10": 1, "01": -1, "11": -1},  # 1 x Z
    {"00": 1, "10": -1, "01": 1, "11": -1},  # Z x 1
    {"00": 1, "10": -1, "01": -1, "11": 1},  # Z x Z
    {"00": 1, "10": -1, "01": -1, "11": 1},  # Y x Y
    {"00": 1, "10": -1, "01": -1, "11": 1},
]  # X x X

qubits = cirq.LineQubit.range(n_qubits)
ansatz_circuit = create_ansatz_circuit(qubits, theta)

ham_mini_circuits = paulis_to_measure_circ(qubits, hamiltonian)

# record min val of hamiltonian
h_min: float = 100.0

# scan theta domain [-pi, pi]
for x in param_range:
    h_avg = get_ham_avg(
        {theta: x}, ansatz_circuit, g_coeff, ham_mini_circuits, obs, shots
    )
    h_avg += g_coeff_id  # add coeff of id term
    graph.append(h_avg)  # collect data for plotting
    if h_avg < h_min:  # record min value
        h_min = h_avg

# result
print("\nH_min = {}\n".format(h_min))

H_min = 0.2252

The routine get_ham_avg() in L104, and more specifically the loop in L121 implements the idea we outlined above. It iterates over all mini-circuits, combining each mini-circuit with the ansatz circuit, running it, collecting the statistics and finally computing the average for the term.

The answer is still wrong because the very last step, computing the average of a single Pauli string in get_paulistr_avg() is still not implemented. Let us fill in that too.

##########################################################
#
#  Routines
#
##########################################################
def create_ansatz_circuit(qbits, parameter):
    """
    Create H2 ansatz circuit given in O'Malley etal 2015.

    Args:
        qbits          - qubits
        parameter      - parameter for parameterized circuit

    Returns:
        ansatz_circuit - parameterized circuit
    """
    ansatz_circuit = cirq.Circuit()

    # H2 circuit
    # 1 - create HF reference state
    ansatz_circuit.append(cirq.rx(np.pi)(qbits[0]), strategy=InsertStrategy.NEW)

    # 2 - ansatz: U(theta)
    ansatz_circuit.append(cirq.rx(-np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    # 3
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 4 - parameterized Z-rotation
    ansatz_circuit.append(cirq.rz(parameter)(qbits[0]), strategy=InsertStrategy.NEW)

    # 5
    ansatz_circuit.append(cirq.CNOT(qbits[1], qbits[0]), strategy=InsertStrategy.NEW)

    # 6
    ansatz_circuit.append(cirq.rx(np.pi / 2)(qbits[0]), strategy=InsertStrategy.NEW)
    ansatz_circuit.append(cirq.ry(-np.pi / 2)(qbits[1]), strategy=InsertStrategy.NEW)

    return ansatz_circuit


def paulis_to_measure_circ(qbits, hamlist):
    """
    Generate measure circuits from hamiltonian pauli list
    for measurements in different bases.

    Args:
        qbits     - qubits
        hamlist   - total hamiltonian, i.e. pauli strings

    Returns:
        circuits - list of circuits that can be used to average over
    """
    circuits = []

    for elem in hamlist:
        minicirc = cirq.Circuit()

        for qubitno, qubitop in elem:
            if qubitop == "Z":
                pass
            elif qubitop == "X":
                minicirc.append(cirq.H(qbits[qubitno]))
            elif qubitop == "Y":
                minicirc.append(cirq.S(qbits[qubitno]) ** -1)
                minicirc.append(cirq.H(qbits[qubitno]))
            else:
                assert False, "Error: INVALID qubit operation"

        minicirc.append(cirq.measure(qbits, key="m"))
        circuits.append(minicirc)

    return circuits


def get_paulistr_avg(trials, obs, shots):
    """
    Calculate average of a Pauli string (e.g. XXY).

    Args:
        trials      - result of simulation: {'00': 45, '01': 34, . . .}
        obs         - data for calculating probabilities
        shots       - number of repetitions

    Returns:
        avg         - the average value
    """
    avg: float = 0.0

    for result_dec_number in trials[0].histogram(key="m"):
        prob = trials[0].histogram(key="m")[result_dec_number] / shots
        binkey = format(result_dec_number, "02b")
        avg += obs[binkey] * prob

    return avg


def get_ham_avg(params, ansatz, coeffs, ham_mini_circuits, obs, shots):
    """
    Compute average of the total qubit Hamiltonian.

    Args:
        params           - parameter value for the ansatz circuit
        ansatz           - ansatz circuit
        ham_mini_ciruits - pre-built circuits to measure Hamiltonian terms
        obs              - data for observables
        shots            - number of shots

    Returns:
        avg   - average of Hamiltonian
    """
    avg: float = 0.0
    sim = cirq.Simulator()

    for coeff, measure_circuit, obsval in zip(coeffs, ham_mini_circuits, obs):
        total_circuit = ansatz + measure_circuit
        cirq_trials = sim.run_sweep(
            program=total_circuit, params=params, repetitions=shots
        )
        paulistring_avg = get_paulistr_avg(cirq_trials, obsval, shots)
        avg += coeff * paulistring_avg

    return avg


##########################################################
#
#  Main program
#
##########################################################

# variables
theta = sympy.Symbol("θ")
n_qubits = 2
hamiltonian = (
    ((0, "Z"),),
    ((1, "Z"),),
    ((0, "Z"), (1, "Z")),
    ((0, "X"), (1, "X")),
    ((0, "Y"), (1, "Y")),
)
shots: int = 1000
no_of_points: int = 40
param_range = np.arange(-np.pi, np.pi, np.pi / no_of_points)
# g_i coeffs for R = 0.75
g_coeff = [0.3435, -0.4347, 0.5716, 0.091, 0.091]
g_coeff_id = 0.2252
# graph
graph = []
#
obs = []
# map outcome states to outcome values
obs = [
    {"00": 1, "10": 1, "01": -1, "11": -1},  # 1 x Z
    {"00": 1, "10": -1, "01": 1, "11": -1},  # Z x 1
    {"00": 1, "10": -1, "01": -1, "11": 1},  # Z x Z
    {"00": 1, "10": -1, "01": -1, "11": 1},  # Y x Y
    {"00": 1, "10": -1, "01": -1, "11": 1},
]  # X x X

qubits = cirq.LineQubit.range(n_qubits)
ansatz_circuit = create_ansatz_circuit(qubits, theta)

ham_mini_circuits = paulis_to_measure_circ(qubits, hamiltonian)

# record min val of hamiltonian
h_min: float = 100.0

# scan theta domain [-pi, pi]
for x in param_range:
    h_avg = get_ham_avg(
        {theta: x}, ansatz_circuit, g_coeff, ham_mini_circuits, obs, shots
    )
    h_avg += g_coeff_id  # add coeff of id term
    graph.append(h_avg)  # collect data for plotting
    if h_avg < h_min:  # record min value
        h_min = h_avg

# result
print("\nH_min = {:.6f}\n".format(h_min))

H_min = -1.151544

This result is now close to the result obtained in arXiv:1512.06860v2 where the circuit was executed on a real, physical quantum computer—compare it with the result in the paper in Fig. 3 at bond length \(R = 0.75\).

Further problems

The article arXiv:1512.06860v2 contains several figures that visualize the results in more detail. Using the above code, you can reproduce them with a little more work. This will give a better understanding of how VQE works.

Exercise 1. Record the average of each Hamiltonian term at \(R = 0.75\) and plot the averages. The resulting graph should be similar to Fig. 2a in arXiv:1512.06860v2.

Exercise 2. Extend the program to find the minimum energies at all bond lengths given in Table 1 in the Appendix. Plot the resulting graph. It should look like Fig. 3a in arXiv:1512.06860v2.