Variational Quantum Eigensolver (VQE) on Molecules

Overview

VQE is a variational algorithm for calculating the ground state of some given hamiltonian H which we call it \(\psi_g\) that satisfied \(H \left|\psi_g\right> =E_g\left|\psi_g\right>\). For an arbitrary normalized wavefunction \(\psi_f\), the expectation value \(\left<\psi_f|H|\psi_f \right>\) is always not lower than the ground state energy unless \(\psi_f = \psi_g\) to some phase factor (here we assume there is no degeneracy in ground state). Based on that fact, if we use a parameterized wavefunction \(\psi_\theta\), e.g. given by a parameterized quantum circuit (PQC) with parameters \(\theta\), we can give an approximation for the ground state enery and wavefunction by minimizing the expectation value of \(H\). In practical quantum hardware, this algorithm can be realized in a quantum-neural hybrid paradigm with the gradient calculated using finite difference or paremeter shift in quantum hardware and the optimization using gradient descent method in classical computer. While in a numerical simulation, we can just calculate the gradients using automatic differentiation.

Calculating the ground state energy for molecules is often important for quantum chemistry tasks since it can be used to find out the atom structure of the molecules. In the simulation of molecules, we do not consider the motion of nuclei which means we fix the nuclear coordinates of its constituent atoms. We only consider the electrons in the molecules since the nuclei are way heavier than the electrons and thus the energy carried by phonons is negligible or can be reconsidered using Born-Oppenheimer approximation. Strictly speaking, the eletrons lives in continuous space, thus the Hilbert space is of infinite dimensions. To conduct a practical calculation, we only preserve some important single-particle basis, e.g. the low energy atomic orbitals. In the second quantization formalism, we can represent these atomic orbitals as \(c_i^\dagger|0>\). By considering the interactions of nuclei and electrons as background and also the electron-electron interaction, a molecules hamiltonian can in generally be represented as \(H = \sum_{i, j} h_{i,j} c_i^\dagger c_j + \sum_{i, j, k, l} \alpha_{i, j, k, l} c_i^\dagger c_j^\dagger c_k c_l\). Notice that the spin index is also absorbed into the orbital index. There are many softwares that can give these parameters in H such as pyscf which we will use later in this tutorial. Now we have a fermionic description for moleculars. By using a mapping from fermionic operators to spin operators such as Jordan-Wigner transformation or Bravyi-Kitaev transformation, we can map the fermionic hamiltonian to a spin hamiltonian which is more compatible with quantum computer. For a spin hamiltonian, we can easily use a PQC to construct a trail wavefunction and conduct the VQE algorithm. In the following part of this tutorial, we will demonstrate a complete example of how to use TensorCircuit to simulate VQE algorithm on Molecules.

Setup

We should first pip install openfermion openfermionpyscf to generate fermionic and qubit Hamiltonian of H2O molecule based on quantum chemistry computation provided by openfermion and pyscf.

import numpy as np
from openfermion.chem import MolecularData
from openfermion.transforms import (
    get_fermion_operator,
    jordan_wigner,
    binary_code_transform,
    checksum_code,
    reorder,
)
from openfermion.chem import geometry_from_pubchem
from openfermion.utils import up_then_down
from openfermion.linalg import LinearQubitOperator
from openfermionpyscf import run_pyscf
import tensorflow as tf

import tensorcircuit as tc

K = tc.set_backend("tensorflow")

Generate Hamiltonian

• Get molecule energy info and molecule orbitals

multiplicity = 1
basis = "sto-3g"
# 14 spin orbitals for H2O
geometry = geometry_from_pubchem("h2o")
description = "h2o"
molecule = MolecularData(geometry, basis, multiplicity, description=description)
molecule = run_pyscf(molecule, run_mp2=True, run_cisd=True, run_ccsd=True, run_fci=True)
print(molecule.fci_energy, molecule.ccsd_energy, molecule.hf_energy)

-75.0155301894916 -75.01540899923558 -74.96444758276998

• Get Fermionic Hamiltonian

mh = molecule.get_molecular_hamiltonian()

fh = get_fermion_operator(mh)

print(fh.terms[((0, 1), (0, 0))])  # coefficient of C0^\dagger C_0

-32.68991541360029

• Transform into qubit Hamiltonian

# The normal transformation such as JW or BK requires 14 qubits for H2O's 14 orbitals

a = jordan_wigner(fh)
LinearQubitOperator(a).n_qubits

14

We can use binary code to save two further qubits, as the number of spin up and spin down filling is both 5 (5/odd electrons in 7 orbitals)

b = binary_code_transform(reorder(fh, up_then_down), 2 * checksum_code(7, 1))
# 7 is 7 spin polarized orbitals, and 1 is for odd occupation
LinearQubitOperator(b).n_qubits

12

print(b.terms[((0, "Z"),)])  # coefficient of Z_0 Pauli-string

12.412562749393349

• Transform the qubit Hamiltonian in openfermion to the format in TensorCircuit

lsb, wb = tc.templates.chems.get_ps(b, 12)
lsa, wa = tc.templates.chems.get_ps(a, 14)

• Inspect Hamiltonian in matrix form

ma = tc.quantum.PauliStringSum2COO_numpy(lsa, wa)

mb = tc.quantum.PauliStringSum2COO_numpy(lsb, wb)

mad, mbd = ma.todense(), mb.todense()

The corresponding Hartree Fock product state in these two types of Hamiltonian

bin(np.argmin(np.diag(mad)))

'0b11111111110000'

bin(np.argmin(np.diag(mbd)))

'0b111110111110'

VQE Setup

We can in principle evaluate each Pauli string of the Hamiltonian as an expectation measurement, but it costs lots of simulation time, instead we fuse them as a Hamiltonian matrix as shown above to run the VQE.

• Using dense matrix expectation

n = 12
depth = 4
mbd_tf = tc.array_to_tensor(mbd)


def vqe(param):
    c = tc.Circuit(n)
    for i in [0, 1, 2, 3, 4, 6, 7, 8, 9, 10]:
        c.X(i)
    for j in range(depth):
        for i in range(n - 1):
            c.exp1(i, i + 1, unitary=tc.gates._xx_matrix, theta=param[j, i, 0])
        for i in range(n):
            c.rx(i, theta=param[j, i, 1])
        for i in range(n):
            c.ry(i, theta=param[j, i, 2])
        for i in range(n):
            c.rx(i, theta=param[j, i, 3])
    return tc.templates.measurements.operator_expectation(c, mbd_tf)

vags = tc.backend.jit(tc.backend.value_and_grad(vqe))
lr = tf.keras.optimizers.schedules.ExponentialDecay(
    decay_rate=0.5, decay_steps=300, initial_learning_rate=0.5e-2
)
opt = tc.backend.optimizer(tf.keras.optimizers.Adam(lr))

param = tc.backend.implicit_randn(shape=[depth, n, 4], stddev=0.02, dtype="float32")
for i in range(600):
    e, g = vags(param)
    param = opt.update(g, param)
    if i % 100 == 0:
        print(e)

tf.Tensor(-74.76671, shape=(), dtype=float32)
tf.Tensor(-74.95493, shape=(), dtype=float32)
tf.Tensor(-74.95319, shape=(), dtype=float32)
tf.Tensor(-74.954315, shape=(), dtype=float32)
tf.Tensor(-74.956116, shape=(), dtype=float32)
tf.Tensor(-74.95809, shape=(), dtype=float32)

• Using sparse matrix expectation

We can also use the sparse Hamiltonian matrix for circuit expectation evaluation, the only difference is to replace mbd_tf with mb_tf

mb_tf = tc.backend.coo_sparse_matrix(
    np.transpose(np.stack([mb.row, mb.col])), mb.data, shape=(2**n, 2**n)
)

A micro-benchmark between sparse matrix evaluation and dense matrix evaluation for expectation in terms of time, sparse always wins in terms of space, of course.

def dense_expt(param):
    c = tc.Circuit(n)
    for i in range(n):
        c.H(i)
        c.rx(i, theta=param[i])
    return tc.templates.measurements.operator_expectation(c, mbd_tf)


def sparse_expt(param):
    c = tc.Circuit(n)
    for i in range(n):
        c.H(i)
        c.rx(i, theta=param[i])
    return tc.templates.measurements.operator_expectation(c, mb_tf)

dense_vag = tc.backend.jit(tc.backend.value_and_grad(dense_expt))
sparse_vag = tc.backend.jit(tc.backend.value_and_grad(sparse_expt))

v0, g0 = dense_vag(tc.backend.ones([n]))
v1, g1 = sparse_vag(tc.backend.ones([n]))

# consistency check

np.testing.assert_allclose(v0, v1, atol=1e-5)
np.testing.assert_allclose(g0, g1, atol=1e-5)

%timeit dense_vag(tc.backend.ones([n]))

30.7 ms ± 1.45 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit sparse_vag(tc.backend.ones([n]))

3.6 ms ± 63 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Therefore, sparse matrix evaluation also saves time apart from space, which is always recommended.