FQE vs OpenFermion vs Cirq: Quadratic Hamiltonian Evolution

View on QuantumAI Run in Google Colab View source on GitHub Download notebook

This notebook demonstrates how evolve a state under a quadratic Hamiltonian with FQE, Cirq, and OpenFermion.

We first generate a random quadratic hamiltonian using OpenFermion. Then we set up a full 2n×2n operator and a computational basis state corresponding to a Slater determinant with the zero and one indexed alpha spin-orbitals occupied. We then compare the full space time evolutions to the LU-decomposition evolution implemented in FQE. Lastly, we using the optimal_givens_decomposition to implement the quadratic Hamiltonian time evolution in Cirq. The 1-RDMs are compared for correctness. We purposefully picked a high spin state so the spin-summed 1-RDM from FQE is equivalent to the spinless 1-RDM.

try:
    import fqe
except ImportError:
    !pip install fqe --quiet

from typing import Union
from itertools import product
import fqe
from fqe.hamiltonians.restricted_hamiltonian import RestrictedHamiltonian

import numpy as np
import cirq

import openfermion as of
from openfermion.testing.testing_utils import random_quadratic_hamiltonian
from openfermion.ops import QuadraticHamiltonian
from openfermion.linalg.givens_rotations import givens_decomposition_square
from openfermion.circuits.primitives import optimal_givens_decomposition

from scipy.sparse import csc_matrix
from scipy.linalg import expm
from scipy.special import binom
import scipy as sp

# set up the number of orbs and a random quadratic Hamiltonian
norbs = 4
ikappa = random_quadratic_hamiltonian(norbs, conserves_particle_number=True, real=True, expand_spin=False, seed=2)
ikappa_mat = of.get_sparse_operator(of.get_fermion_operator(ikappa), n_qubits=norbs)

# set up initial full Hilbert space wavefunction
wfn = np.zeros((2**(norbs), 1), dtype=np.complex128)
wfn[int("".join([str(x) for x in [1] * (norbs//2) + [0] * (norbs//2)]), 2), 0] = 1  # alpha1-alpha2 Hartree-Fock state

# full space unitary
time = np.random.random()
hamiltonian_evolution = expm(-1j * ikappa_mat.todense() * time)

# build a FQE operator corresponding to the Hamiltonian evolution
fqe_ham = RestrictedHamiltonian((ikappa.n_body_tensors[1, 0],))
assert fqe_ham.quadratic()

# initialize the FQE wavefunction [n-elec, sz, norbs]
fqe_wfn = fqe.Wavefunction([[norbs // 2, norbs//2, norbs]])
hf_wf = np.zeros((int(binom(norbs, norbs // 2)), 1), dtype=np.complex128)
hf_wf[0, 0] = 1  # right most bit is zero orbital.
fqe_wfn.set_wfn(strategy='from_data',
                raw_data={(norbs//2, norbs//2): hf_wf})
fqe_wfn.print_wfn()

Sector N = 2 : S_z = 2
a'0011'b'0000' (1+0j)

Aside

The FQE and cirq/OpenFermion have reverse ordering of the occupation number vectors. Give a bitstring representing occupations 01001 one commonly wants to determine the index of the orbitals used in the 2-orbital Slater determinant. In OpenFermion the orbitals used are 1 and 4. In FQE the orbitals used are 0 and 3. The difference is that OpenFermion starts indexing from the left most bit whereas FQE starts indexing from the right most bit. 

# initialize wavefunction [n-elec, sz, norbs]
fqe_wfn = fqe.Wavefunction([[norbs // 2, norbs//2, norbs]])
hf_wf = np.zeros((int(binom(norbs, norbs // 2)), 1), dtype=np.complex128)
hf_wf[0, 0] = 1  # right most bit is zero orbital. 2-alpha electron HF is 0011.  This is the first element of our particle conserved Hilbert space
fqe_wfn.set_wfn(strategy='from_data',
                raw_data={(norbs//2, norbs//2): hf_wf})
fqe_wfn.print_wfn()

# OpenFermion indexes from the left as zero.  The RHF state is 1100
wfn[int("".join([str(x) for x in [1] * (norbs//2) + [0] * (norbs//2)]), 2), 0] = 1

print("\nHF-Vector-OpenFermion ", bin(int("".join([str(x) for x in [1] * (norbs//2) + [0] * (norbs//2)]), 2)))

Sector N = 2 : S_z = 2
a'0011'b'0000' (1+0j)

HF-Vector-OpenFermion  0b1100

fqe_wfn.print_wfn()  # before
fqe_wfn = fqe_wfn.time_evolve(time, fqe_ham)
fqe_wfn.print_wfn()  # after

Sector N = 2 : S_z = 2
a'0011'b'0000' (1+0j)
Sector N = 2 : S_z = 2
a'0011'b'0000' (0.7640756559841199+0.20777158789940361j)
a'0101'b'0000' (-0.17313823898087422+0.2638365952059404j)
a'1001'b'0000' (0.40785856832720757+0.06983556964686519j)
a'0110'b'0000' (-0.10456038050286236-0.10035394746174761j)
a'1010'b'0000' (-0.1978764858966419+0.19803152903956048j)
a'1100'b'0000' (0.008829529249272671-0.05247669185209847j)

final_wfn = hamiltonian_evolution @ wfn  # Full space time-evolution

Compare 1-RDMs

It is easy to compare the 1-RDMs of the wavefunctions to confirm they are representing the same state. We will compute the 1-RDM for the wfn and fqe_wfn and compare to the opdm computed from fqe.to_cirq.

def get_opdm(wfn: Union[sp.sparse.coo_matrix, sp.sparse.csr_matrix,
                        sp.sparse.csc_matrix], nso: int) -> np.ndarray:
    """
    A slow way to get 1-RDMs from density matrices or states
    """
    if isinstance(wfn, (sp.sparse.coo_matrix, sp.sparse.csc_matrix)):
        wfn = wfn.tocsr()
        wfn = wfn.toarray()
    opdm = np.zeros((nso, nso), dtype=wfn.dtype)
    a_ops = [
        of.get_sparse_operator(of.jordan_wigner(of.FermionOperator(((p, 0)))),
                               n_qubits=nso)
        for p in range(nso)
    ]
    for p, q in product(range(nso), repeat=2):
        operator = a_ops[p].conj().T @ a_ops[q]
        if wfn.shape[0] == wfn.shape[1]:  # we are working with a density matrix
            val = np.trace(wfn @ operator)
        else:
            val = wfn.conj().transpose() @ operator @ wfn
        # print((p, q, r, s), val.toarray()[0, 0])
        if isinstance(val, (float, complex, int)):
            opdm[p, q] = val
        else:
            opdm[p, q] = val[0, 0]

    return opdm

fqe_opdm = fqe_wfn.rdm('i^ j') # get the FQE-opdm

# get the exactly evolved opdm
initial_opdm = np.diag([1] * (norbs//2) + [0] * (norbs//2))
u = expm(1j * ikappa.n_body_tensors[1, 0] * time)
final_opdm = u @ initial_opdm @ u.conj().T
print(final_opdm)
print()

# contract the wavefunction opdm
test_opdm = get_opdm(final_wfn, norbs)
print(test_opdm)

assert np.allclose(test_opdm.real, final_opdm.real)
assert np.allclose(test_opdm.imag, final_opdm.imag)
assert np.allclose(final_opdm, fqe_opdm)

cirq_wfn = fqe.to_cirq(fqe_wfn).reshape((-1, 1))

cfqe_opdm = get_opdm(cirq_wfn, 2 * norbs)
print()
print(cfqe_opdm[::2, ::2])
assert np.allclose(final_opdm, cfqe_opdm[::2, ::2])

# cirq evolution
qubits = cirq.LineQubit.range(norbs)
prep = cirq.Moment([cirq.X.on(qubits[0]), cirq.X.on(qubits[1])])
rotation = cirq.Circuit(optimal_givens_decomposition(qubits, u.conj().copy()))
circuit = prep + rotation

final_state = circuit.final_state_vector().reshape((-1, 1))
cirq_opdm = get_opdm(final_state, norbs)
print()
print(cirq_opdm)

assert np.allclose(final_opdm, cirq_opdm)

[[ 0.89779286+0.j         -0.07524962+0.13954963j  0.1006792 +0.03293365j
   0.12542128-0.19918035j]
 [-0.07524962-0.13954963j  0.72635602+0.j         -0.08961216+0.2461997j
   0.32180163-0.03775482j]
 [ 0.1006792 -0.03293365j -0.08961216-0.2461997j   0.12342215+0.j
  -0.05137394-0.16026316j]
 [ 0.12542128+0.19918035j  0.32180163+0.03775482j -0.05137394+0.16026316j
   0.25242897+0.j        ]]

[[ 0.89779286+0.j         -0.07524962+0.13954963j  0.1006792 +0.03293365j
   0.12542128-0.19918035j]
 [-0.07524962-0.13954963j  0.72635602+0.j         -0.08961216+0.2461997j
   0.32180163-0.03775482j]
 [ 0.1006792 -0.03293365j -0.08961216-0.2461997j   0.12342215+0.j
  -0.05137394-0.16026316j]
 [ 0.12542128+0.19918035j  0.32180163+0.03775482j -0.05137394+0.16026316j
   0.25242897+0.j        ]]

[[ 0.89779286+0.j         -0.07524962+0.13954963j  0.1006792 +0.03293365j
   0.12542128-0.19918035j]
 [-0.07524962-0.13954963j  0.72635602+0.j         -0.08961216+0.2461997j
   0.32180163-0.03775482j]
 [ 0.1006792 -0.03293365j -0.08961216-0.2461997j   0.12342215+0.j
  -0.05137394-0.16026316j]
 [ 0.12542128+0.19918035j  0.32180163+0.03775482j -0.05137394+0.16026316j
   0.25242897+0.j        ]]

[[ 0.89779286+0.j         -0.07524962+0.13954963j  0.1006792 +0.03293365j
   0.12542128-0.19918035j]
 [-0.07524962-0.13954963j  0.72635602+0.j         -0.08961216+0.2461997j
   0.32180163-0.03775482j]
 [ 0.1006792 -0.03293365j -0.08961216-0.2461997j   0.12342215+0.j
  -0.05137394-0.16026316j]
 [ 0.12542128+0.19918035j  0.32180163+0.03775482j -0.05137394+0.16026316j
   0.25242897+0.j        ]]

Givens rotations with FQE

We can use the same circuit generated by OpenFermion except apply the Givens rotations to the FQE wavefunction. Each Givens rotation gate sequence is mapped to

e−iθ(a†paq−a†qap)eiϕa†pap

which can be applies sequentially to the state.

# Get the basis rotation circuit
rotations, diagonal = givens_decomposition_square(u.conj().T)

# Reinitialize the wavefunction to Hartree-Fock
fqe_wfn = fqe.Wavefunction([[norbs // 2, norbs // 2, norbs]])
hf_wf = np.zeros((int(binom(norbs, norbs // 2)), 1),
                 dtype=np.complex128)
hf_wf[0, 0] = 1  # right most bit is zero orbital.
fqe_wfn.set_wfn(strategy='from_data',
                raw_data={(norbs // 2, norbs // 2): hf_wf})

# Iterate through each layer and time evolve by the appropriate fermion operators
for layer in rotations:
    for givens in layer:
        i, j, theta, phi = givens
        op = of.FermionOperator(((2 * j, 1), (2 * j, 0)), coefficient=-phi)
        fqe_wfn = fqe_wfn.time_evolve(1.0, op)

        op = of.FermionOperator(((2 * i, 1), (2 * j, 0)), coefficient=-1j * theta) + of.FermionOperator(((2 * j, 1), (2 * i, 0)), coefficient=1j * theta)
        fqe_wfn = fqe_wfn.time_evolve(1.0, op)

# evolve the last diagonal phases
for idx, final_phase in enumerate(diagonal):
    if not np.isclose(final_phase, 1.0):
        op = of.FermionOperator(((2 * idx, 1), (2 * idx, 0)), -np.angle(final_phase))
        fqe_wfn = fqe_wfn.time_evolve(1.0, op)

# map the wavefunction to a cirq wavefunction and get the opdm
crq_wfn = fqe.to_cirq(fqe_wfn).reshape((-1, 1))
crq_opdm = get_opdm(crq_wfn, 2 * norbs)
print(crq_opdm[::2, ::2])

# check if it is the same as our other 1-RDMs
assert np.allclose(final_opdm, crq_opdm[::2, ::2])

[[ 0.89779286+0.j         -0.07524962+0.13954963j  0.1006792 +0.03293365j
   0.12542128-0.19918035j]
 [-0.07524962-0.13954963j  0.72635602+0.j         -0.08961216+0.2461997j
   0.32180163-0.03775482j]
 [ 0.1006792 -0.03293365j -0.08961216-0.2461997j   0.12342215+0.j
  -0.05137394-0.16026316j]
 [ 0.12542128+0.19918035j  0.32180163+0.03775482j -0.05137394+0.16026316j
   0.25242897+0.j        ]]