A General Case

Author(s): Kyle Godbey

Maintainer: Kyle Godbey

To take this particular problem to its natural conclusion, we will be now devising a way to automatically generate our Hamiltonian and ansatz circuits based on what basis dimension we want to calculate.

To do this in the most straightforward manner, I will first implement the Hamiltonian in second quantization before mapping it to a Pauli representation. To help out, I’ll be importing some tools from qiskit as well as pennylane, proving that frameworks can indeed work together!

%matplotlib inline

import matplotlib.pyplot as plt
from pennylane import numpy as np
import pennylane as qml
from qiskit_nature.operators.second_quantization import FermionicOp
from qiskit_nature.mappers.second_quantization import JordanWignerMapper, ParityMapper
import time
from functools import partial
import warnings
warnings.filterwarnings('ignore')

plt.style.use(['science','notebook'])

We’ll also need some helper functions to parse Pauli strings and construct our pennylane Hamiltonian, so I’ll just dump them here.

def pauli_token_to_operator(token):
    qubit_terms = []

    for term in range(len(token)):
        # Special case of identity
        if token[term] == "I":
            pass
        else:
            #pauli, qubit_idx = term, term
            if token[term] == "X":
                qubit_terms.append(qml.PauliX(int(term)))
            elif token[term] == "Y":
                qubit_terms.append(qml.PauliY(int(term)))
            elif token[term] == "Z":
                qubit_terms.append(qml.PauliZ(int(term)))
            else:
                print("Invalid input.")
    if(qubit_terms==[]):
            qubit_terms.append(qml.Identity(0))
    full_term = qubit_terms[0]
    for term in qubit_terms[1:]:
        full_term = full_term @ term

    return full_term


def parse_hamiltonian_input(input_data):
    # Get the input
    coeffs = []
    pauli_terms = []
    chunks = input_data.split("\n")
    # Go through line by line and build up the Hamiltonian
    for line in chunks:
        #line = line.strip()
        tokens = line.split(" ")
        # Parse coefficients
        sign, value = tokens[0][0], tokens[1]

        coeff = float(value)
        if sign == "-":
            coeff *= -1
        coeffs.append(coeff)

        # Parse Pauli component
        pauli = tokens[3][::-1]
        
        pauli_terms.append(pauli_token_to_operator(pauli))

    return qml.Hamiltonian(coeffs, pauli_terms)

Something else that will be very important is a function to compute the matrix elements \(\langle n' | (\hat{T} + \hat{V}) | n \rangle\) needed for our Hamiltonian.

def kron(i,j):
    if(i==j):
        return 1
    else:
        return 0

def matrix_element(i,j):
    
    ele = 0.0
    
    ele = ((2.0*i + 1.5)*kron(i,j) - np.sqrt(i*(i+0.5))*kron(i,j+1) \
    - np.sqrt((i+1)*(i+1.5))*kron(i,j-1)) * 3.5
    
    ele += -5.68658111 * kron(i,0) * kron(i,j)
    
    
    return ele

Lastly, we’ll make a function that returns a pennylane compatible Hamiltonian given a requested basis dimension, N.

def deuteron_ham(N,mapper=JordanWignerMapper):
    # Start out by zeroing what will be our fermionic operator
    op = 0
    for i in range(N):
        for j in range(N):
            # Construct the terms of the Hamiltonian in terms of creation/annihilation operators 
            op += matrix_element(i,j) * \
            FermionicOp([([("+", i),("-", j)], 1.0)])

    hamstr = "+ "+str(mapper().map(second_q_op=op))

    hamiltonian = parse_hamiltonian_input(hamstr)

    return hamiltonian 

Let’s take a look at our new Hamiltonian generator in the next code block. We’ll start by setting the basis dimension to 3 so we can compare our two Hamiltonians: handcrafted and computer generated.

dim = 3

# Building our Hamiltonian for the N=3 case, as before

coeffs = [15.531709 ,-2.1433, -2.1433, 0.21829, -6.125, -9.625, -3.91, -3.91]
obs = [qml.Identity(0),qml.PauliX(0) @ qml.PauliX(1), qml.PauliY(0) @ qml.PauliY(1),qml.PauliZ(0),qml.PauliZ(1),qml.PauliZ(2),qml.PauliX(1) @ qml.PauliX(2),qml.PauliY(1) @ qml.PauliY(2)]

ham = qml.Hamiltonian(coeffs, obs)

# Let's print it out

print("Original Hamiltonian: \n",ham)

# Let's now use our new function!

H = deuteron_ham(dim)

print("Generated Hamiltonian: \n",H)

Original Hamiltonian: 
   (-9.625) [Z2]
+ (-6.125) [Z1]
+ (0.21829) [Z0]
+ (15.531709) [I0]
+ (-3.91) [X1 X2]
+ (-3.91) [Y1 Y2]
+ (-2.1433) [X0 X1]
+ (-2.1433) [Y0 Y1]
Generated Hamiltonian: 
   (-9.625) [Z2]
+ (-6.125) [Z1]
+ (0.21829055499999983) [Z0]
+ (15.531709445) [I0]
+ (-3.913118960624632) [Y1 Y2]
+ (-3.913118960624632) [X1 X2]
+ (-2.1433035249352805) [Y0 Y1]
+ (-2.1433035249352805) [X0 X1]

We still need to define our device and circuit, so we’ll let them be general as well.

# Set the order you'd like to go to

dim = 5

dev = qml.device("default.qubit", wires=dim)

# Define a general ansatz for arbitrary numbers of dimensions

particles = 1

ref_state = qml.qchem.hf_state(particles, dim)

ansatz = partial(qml.ParticleConservingU2, init_state=ref_state)

layers = dim - 2

@qml.qnode(dev)
def draw_circuit(params):
    ansatz(params,wires=dev.wires).decomposition()
    return qml.expval(H)

# Defining Hamiltonian

H = deuteron_ham(dim)

print(H)

cost_fn = qml.ExpvalCost(ansatz, H, dev)

qml.drawer.use_style('default')
print(qml.draw_mpl(draw_circuit)(init_params))

  (-16.625) [Z4]
+ (-13.125) [Z3]
+ (-9.625) [Z2]
+ (-6.125) [Z1]
+ (0.21829055499999983) [Z0]
+ (45.281709445000004) [I0]
+ (-7.424621202458749) [Y3 Y4]
+ (-7.424621202458749) [X3 X4]
+ (-5.670648111106878) [Y2 Y3]
+ (-5.670648111106878) [X2 X3]
+ (-3.913118960624632) [Y1 Y2]
+ (-3.913118960624632) [X1 X2]
+ (-2.1433035249352805) [Y0 Y1]
+ (-2.1433035249352805) [X0 X1]
(<Figure size 2600x600 with 1 Axes>, <Axes:>)

Quintuple Amazing! That’s yet another good looking Hamiltonian, and the circuit isn’t bad either…

Now we’ll set up what we’ll need for the VQE procedure, this time doing a loop over how many parameters we’ll need.

# Our parameter array

init_params = np.random.uniform(low=-np.pi / 2, high=np.pi / 2, size=qml.ParticleConservingU2.shape(n_layers=layers, n_wires=dim))

# Convergence information and step size

max_iterations = 500
conv_tol = 1e-05
step_size = 0.1

Finally, the VQE block. We’re still using the standard gradient descent optimizer since it worked so well before, but soon the time will come to shop around for better options. I have also added a quick and dirty variable step size just to speed things up slightly.

opt = qml.GradientDescentOptimizer(stepsize=step_size)

params = init_params

gd_param_history = [params]
gd_cost_history = []

accel = 0
prev_conv = -1.0

start = time.time()

for n in range(max_iterations):
    fac = (1.0)
    opt = qml.GradientDescentOptimizer(stepsize=step_size)

    # Take a step in parameter space and record your energy
    params, prev_energy = opt.step_and_cost(cost_fn, params)

    # This keeps track of our energy for plotting at comparisons
    gd_param_history.append(params)
    gd_cost_history.append(prev_energy)

    # Here we see what the energy of our system is with the new parameters
    energy = cost_fn(params)

    # Calculate difference between new and old energies
    conv = np.abs(energy - prev_energy)

    if(energy - prev_energy > 0.0 and step_size > 0.001):
        #print("Lowering!")
        accel = 0
        step_size = 0.5*step_size

    if(conv < prev_conv): accel += 1
    prev_conv = conv
    
    if(accel > 10 and step_size < 1.0):
        #print("Accelerating!")
        step_size = 1.1*step_size
    end = time.time()

    if n % 10 == 0:
        print(
            "It = {:},  Energy = {:.8f} MeV,  Conv = {"
            ":.8f} MeV, Time Elapsed = {:.3f} s".format(n, energy, conv,end-start)
        )
        start = time.time()


    if conv <= conv_tol:
        break

print()
print("Final value of the energy = {:.8f} MeV".format(energy))
print("Number of iterations = ", n)

It = 0,  Energy = 25.43598895 MeV,  Conv = 3.74111148 MeV, Time Elapsed = 0.811 s
It = 10,  Energy = -1.84990698 MeV,  Conv = 0.01447123 MeV, Time Elapsed = 8.448 s
It = 20,  Energy = -1.98298605 MeV,  Conv = 0.01369684 MeV, Time Elapsed = 8.302 s
It = 30,  Energy = -2.14395469 MeV,  Conv = 0.01093344 MeV, Time Elapsed = 8.459 s
It = 40,  Energy = -2.16828413 MeV,  Conv = 0.00097018 MeV, Time Elapsed = 8.162 s
It = 50,  Energy = -2.17146313 MeV,  Conv = 0.00033026 MeV, Time Elapsed = 8.317 s
It = 60,  Energy = -2.17539857 MeV,  Conv = 0.00686297 MeV, Time Elapsed = 8.247 s
It = 70,  Energy = -2.17669224 MeV,  Conv = 0.00011480 MeV, Time Elapsed = 8.408 s
It = 80,  Energy = -2.17847872 MeV,  Conv = 0.00023824 MeV, Time Elapsed = 8.537 s
It = 90,  Energy = -2.18142522 MeV,  Conv = 0.00030970 MeV, Time Elapsed = 8.183 s
It = 100,  Energy = -2.18264553 MeV,  Conv = 0.00060788 MeV, Time Elapsed = 8.361 s
It = 110,  Energy = -2.18294148 MeV,  Conv = 0.00002063 MeV, Time Elapsed = 8.373 s
It = 120,  Energy = -2.18321824 MeV,  Conv = 0.00003227 MeV, Time Elapsed = 8.444 s

Final value of the energy = -2.18336450 MeV
Number of iterations =  129

If you choose a bigger basis, you do a better job! This is to be expected, but notice that we’ve taken quite the performance hit by going to a more general representation. We’ll explore how to do better at this in the optimization section, but for now let’s check the convergence.

plt.plot(gd_cost_history, "b", label="Gradient descent")

plt.ylabel("Energy (MeV)")
plt.xlabel("Optimization steps")
plt.legend()
plt.show()

Since we are now allowed to move in higher dimensions, making a general PES plotter gets a little tricky, but feel free to play around and see if you can think of clever ways to plot how your solution explores the space!

This pretty much concludes our exploration of quantum computing for the deuteron from the point of view of Hamiltonian and state preparation, but there’s plenty of other things that can be tried to improve performance! One approach would involve optimizing the ansatz instead of using a generic form, which could be a fun challenge. Another approach is to speed up the VQE algorithm itself by being clever with how we iterate which is what we’ll briefly explore next.