 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 \(2^{n} \times 2^{n}\) 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.7640756559841229+0.20777158789940006j) a'0101'b'0000' (-0.1731382389808739+0.26383659520594166j) a'1001'b'0000' (0.40785856832720846+0.06983556964686378j) a'0110'b'0000' (-0.10456038050286318-0.10035394746174689j) a'1010'b'0000' (-0.1978764858966419+0.19803152903956175j) a'1100'b'0000' (0.008829529249272675-0.05247669185209916j)
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 \theta (a_{p}^{\dagger}a_{q} - a_{q}^{\dagger}a_{p})}e^{i \phi a_{p}^{\dagger}a_{p} } \]
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 ]]