Learn OpenFermion-FQE
Simulating the 1D Fermi-Hubbard model in OpenFermion-FQE.
In [1]:
try:
import fqe
except:
!pip install git+https://github.com/quantumlib/OpenFermion-FQE --quiet
In [2]:
import numpy as np
from scipy import sparse
from scipy.linalg import expm
import openfermion as of
import fqe
import matplotlib.pyplot as plt
try:
plt.style.use('publish')
except:
pass
In [36]:
import copy
In [42]:
from fqe.algorithm.low_rank import evolve_fqe_givens, evolve_fqe_charge_charge_alpha_beta
from fqe.algorithm.low_rank import evolve_fqe_givens_sector
In [3]:
class model(dict):
def __getattr__(self, name):
try:
return self[name]
except KeyError as e:
raise AttributeError(name) from e
__setattr__ = dict.__setitem__
__delattr__ = dict.__delitem__
In [4]:
params = model(nsites=4, U = 2.0, t = 1.0)
hubbard = of.fermi_hubbard(1, params.nsites,
tunneling=params.t,
coulomb=params.U,
periodic=False)
In [5]:
nelectrons, Sz = 2, 0;
init_wavefunction = fqe.Wavefunction([[nelectrons,
Sz,
params.nsites]])
init_wavefunction.set_wfn(strategy='random')
In [13]:
%%time
evolve_time = 1.0;
true_evolved_fqe = init_wavefunction.time_evolve(evolve_time, hubbard);
CPU times: user 21.4 ms, sys: 2.8 ms, total: 24.2 ms Wall time: 23.7 ms
In [14]:
%%time
init_cirq_wfn=fqe.to_cirq(init_wavefunction)
# U = expm(-i*t*H)
H_matrix = of.get_sparse_operator(hubbard).todense()
unitary = expm(-1j*evolve_time*H_matrix)
true_evolved_cirq = unitary@init_cirq_wfn
CPU times: user 132 ms, sys: 10.5 ms, total: 142 ms Wall time: 109 ms
In [16]:
fidelty=abs(fqe.vdot(true_evolved_fqe,
fqe.from_cirq(true_evolved_cirq,
thresh=1e-12)))**2
print('F-1 = ', abs(fidelty-1.0))
F-1 = 4.440892098500626e-16
In [23]:
nqubits = 2*params.nsites
one_body = [ op + of.hermitian_conjugated(op)
for op in (of.FermionOperator(((i,1), (i+2,0)),
coefficient= -1*params.t) for i in range(nqubits-2))]
two_body = [of.FermionOperator(((i,1),(i,0), (i+1,1),(i+1,0)), coefficient=params.U)
for i in range(0, nqubits, 2)]
In [25]:
hopping_matrix = of.get_sparse_operator(sum(one_body))
unitary = expm(-1j*evolve_time*hopping_matrix.todense())
evolved_cirq_wfn = unitary @ init_cirq_wfn
In [27]:
hopping_matrix = -params.t*(np.diag([1]*(params.nsites-1), k=1)
+np.diag([1]*(params.nsites-1), k=-1))
In [28]:
umat = expm(-1j*evolve_time*hopping_matrix)
evolved_wavefunction=evolve_fqe_givens(init_wavefunction, umat)
In [30]:
fidelity=abs(fqe.vdot(evolved_wavefunction,
fqe.from_cirq(evolved_cirq_wfn,thresh=1e-6)))**2
print('F-1 = ', abs(fidelity-1))
F-1 = 0.0
In [31]:
coulomb_matrix=of.get_sparse_operator(sum(two_body))
unitary=expm(-1j*evolve_time*coulomb_matrix.todense())
evolved_cirq_wfn=unitary@init_cirq_wfn
In [32]:
coulomb_matrix=np.diag([params.U]*params.nsites)
In [34]:
evolved_fqe_wfn=evolve_fqe_charge_charge_alpha_beta(init_wavefunction,
coulomb_matrix,
evolve_time)
In [35]:
fidelity = abs(fqe.vdot(evolved_fqe_wfn, fqe.from_cirq(evolved_cirq_wfn, thresh=1e-6)))**2
print('F-1 = ', abs(fidelty-1))
F-1 = 4.440892098500626e-16
In [37]:
steps=10
umat = expm(-1j*hopping_matrix*(evolve_time/steps))
current_wfn=copy.deepcopy(init_wavefunction)
for st in range(steps):
current_wfn=evolve_fqe_givens(current_wfn, u=umat)
current_wfn=evolve_fqe_charge_charge_alpha_beta(
current_wfn,
coulomb_matrix,
evolve_time/steps
)
In [38]:
trotterized_fidelity = abs(fqe.vdot(true_evolved_fqe, current_wfn))**2
print(trotterized_fidelity)
0.9988722601807759
In [68]:
params = model(nsites = 8,
l_up = 4,
m_up = 4.5,
sigma_up = 1,
U = 3.0
)
site_index = np.arange(1, params.nsites+1)
Hσ_up = np.diag([-1.]*(params.nsites-1), k=1)+np.diag([-1.] *(params.nsites - 1), k=-1)
Hσ_up += np.diag(-params.l_up*np.exp(-0.5*(site_index-params.m_up)**2)/params.sigma_up**2)
Hσ_dn = np.diag([-1.]*(params.nsites-1), k=1)+np.diag([-1.] *(params.nsites - 1), k=-1)
In [69]:
init_wfn = fqe.Wavefunction([[params.nsites//2, 0, params.nsites]])
init_wfn.set_wfn(strategy='hartree-fock')
init_wfn.print_wfn()
Sector N = 4 : S_z = 0 a'00000011'b'00000011' (1+0j)
In [70]:
# find the eigenbasis
vup = np.linalg.eigh(Hσ_up)[1]
vdn = np.linalg.eigh(Hσ_dn)[1]
In [71]:
init_wfn = evolve_fqe_givens_sector(init_wfn, vup, sector='alpha')
init_wfn = evolve_fqe_givens_sector(init_wfn, vdn, sector='beta')
In [94]:
opdm_a, opdm_b = init_wfn.sector((params.nsites // 2, 0)).get_spin_opdm()
dt = 0.3
steps = 55
hopping_u = expm(-1j*dt*Hσ_dn)
charge_u = np.diag([params.U]*params.nsites)
opdms_alpha=[opdm_a]
opdms_beta=[opdm_b]
times=[0.0]
current_time = 0.0
final_wfn = init_wfn
for t in range(steps):
final_wfn=evolve_fqe_givens(final_wfn, hopping_u)
final_wfn=evolve_fqe_charge_charge_alpha_beta(final_wfn, charge_u, dt)
current_time += dt
opdm_a_current, opdm_b_current = final_wfn.sector((params.nsites//2,0)).get_spin_opdm()
opdms_alpha.append(opdm_a_current)
opdms_beta.append(opdm_b_current)
times.append(current_time)
In [116]:
select_times=[0.0, 1.2, 1.8, 3.]
time_indices=[int(t/dt) for t in select_times]
charge = [np.diagonal(opdms_alpha[idx]+opdms_beta[idx]).real
for idx in time_indices]
spin = [np.diagonal(opdms_alpha[idx]-opdms_beta[idx]).real
for idx in time_indices]
In [138]:
fig,ax = plt.subplots(2,len(select_times), sharex=True, figsize=(12, 5))
for i,c,s in zip(range(len(select_times)), charge, spin):
ax[0,i].plot(site_index, c, 'o-', mfc='none', color='tab:blue')
ax[1,i].plot(site_index, s, 'o-', mfc='none', color='tab:red')
ax[0,i].set_ylim(0, 1.1)
ax[1,i].set_ylim(-0.6, 0.6)
ax[1,i].set_xticks(site_index)
ax[1,i].set_xlabel('site index')
if i != 0: ax[0, i].set_yticks([]);
if i != 0: ax[1, i].set_yticks([]);
ax[0,i].set_title(r't = {:1.1f} ($\hbar/t$)'.format(select_times[i]))
ax[0,0].set_ylabel(r'charge density $\rho^{+}$')
ax[1,0].set_ylabel(r'spin density $\rho^{-}$')
plt.savefig('snapshot.pdf')
In [150]:
nup_all_time = np.array([np.diagonal(alpha).real for (alpha, beta) in zip(opdms_alpha, opdms_beta)])
ndn_all_time = np.array([np.diagonal(beta).real for (alpha, beta) in zip(opdms_alpha, opdms_beta)])
#spin_all_time = np.array([np.diagonal(alpha-beta).real for (alpha, beta) in zip(opdms_alpha, opdms_beta)])
In [153]:
nup_all_time
Out[153]:
array([[0.002055 , 0.02196527, 0.18623579, 0.78974394, 0.78974394,
0.18623579, 0.02196527, 0.002055 ],
[0.00223427, 0.02967753, 0.24682716, 0.72126104, 0.72126104,
0.24682716, 0.02967753, 0.00223427],
[0.00524075, 0.06088448, 0.34364209, 0.59023269, 0.59023269,
0.34364209, 0.06088448, 0.00524075],
[0.01840544, 0.12211441, 0.39455177, 0.46492837, 0.46492837,
0.39455177, 0.12211441, 0.01840544],
[0.05298825, 0.18757967, 0.39686223, 0.36256986, 0.36256986,
0.39686223, 0.18757967, 0.05298825],
[0.114017 , 0.22367181, 0.37639258, 0.28591861, 0.28591861,
0.37639258, 0.22367181, 0.114017 ],
[0.19230246, 0.2344619 , 0.32605593, 0.24717971, 0.24717971,
0.32605593, 0.2344619 , 0.19230246],
[0.26898591, 0.2363266 , 0.25235402, 0.24233347, 0.24233347,
0.25235402, 0.2363266 , 0.26898591],
[0.32883715, 0.22461168, 0.20179626, 0.24475491, 0.24475491,
0.20179626, 0.22461168, 0.32883715],
[0.36169906, 0.2098958 , 0.1912606 , 0.23714455, 0.23714455,
0.1912606 , 0.2098958 , 0.36169906],
[0.35555807, 0.22302191, 0.19745801, 0.22396201, 0.22396201,
0.19745801, 0.22302191, 0.35555807],
[0.31039463, 0.26080961, 0.21957591, 0.20921985, 0.20921985,
0.21957591, 0.26080961, 0.31039463],
[0.24980318, 0.29650348, 0.25336209, 0.20033124, 0.20033124,
0.25336209, 0.29650348, 0.24980318],
[0.19662206, 0.32812046, 0.25929787, 0.21595961, 0.21595961,
0.25929787, 0.32812046, 0.19662206],
[0.15573637, 0.35619217, 0.24107294, 0.24699853, 0.24699853,
0.24107294, 0.35619217, 0.15573637],
[0.13071431, 0.36068206, 0.24999333, 0.2586103 , 0.2586103 ,
0.24999333, 0.36068206, 0.13071431],
[0.12747294, 0.33884156, 0.2856799 , 0.2480056 , 0.2480056 ,
0.2856799 , 0.33884156, 0.12747294],
[0.14011813, 0.31375264, 0.31129835, 0.23483088, 0.23483088,
0.31129835, 0.31375264, 0.14011813],
[0.15952901, 0.30071896, 0.3262153 , 0.21353673, 0.21353673,
0.3262153 , 0.30071896, 0.15952901],
[0.18395151, 0.30388074, 0.33836872, 0.17379903, 0.17379903,
0.33836872, 0.30388074, 0.18395151],
[0.2073382 , 0.32930919, 0.32817918, 0.13517343, 0.13517343,
0.32817918, 0.32930919, 0.2073382 ],
[0.22887904, 0.37033366, 0.28113719, 0.11965011, 0.11965011,
0.28113719, 0.37033366, 0.22887904],
[0.27272756, 0.38723288, 0.21747348, 0.12256608, 0.12256608,
0.21747348, 0.38723288, 0.27272756],
[0.34922851, 0.34690639, 0.17463427, 0.12923084, 0.12923084,
0.17463427, 0.34690639, 0.34922851],
[0.41347239, 0.28435121, 0.1706201 , 0.1315563 , 0.1315563 ,
0.1706201 , 0.28435121, 0.41347239],
[0.41835912, 0.25992641, 0.19428957, 0.1274249 , 0.1274249 ,
0.19428957, 0.25992641, 0.41835912],
[0.37505109, 0.27496079, 0.22587388, 0.12411424, 0.12411424,
0.22587388, 0.27496079, 0.37505109],
[0.31281579, 0.30180458, 0.24817686, 0.13720278, 0.13720278,
0.24817686, 0.30180458, 0.31281579],
[0.23469907, 0.34613848, 0.24872891, 0.17043354, 0.17043354,
0.24872891, 0.34613848, 0.23469907],
[0.15141908, 0.39566343, 0.23928078, 0.21363671, 0.21363671,
0.23928078, 0.39566343, 0.15141908],
[0.0985753 , 0.38888043, 0.24923712, 0.26330715, 0.26330715,
0.24923712, 0.38888043, 0.0985753 ],
[0.09288278, 0.31233079, 0.27935238, 0.31543405, 0.31543405,
0.27935238, 0.31233079, 0.09288278],
[0.1140928 , 0.22808886, 0.29386725, 0.36395109, 0.36395109,
0.29386725, 0.22808886, 0.1140928 ],
[0.13753763, 0.17906144, 0.27031277, 0.41308816, 0.41308816,
0.27031277, 0.17906144, 0.13753763],
[0.15189908, 0.16737755, 0.23434065, 0.44638272, 0.44638272,
0.23434065, 0.16737755, 0.15189908],
[0.15118929, 0.19377334, 0.22516054, 0.42987683, 0.42987683,
0.22516054, 0.19377334, 0.15118929],
[0.14238955, 0.23854862, 0.24167833, 0.3773835 , 0.3773835 ,
0.24167833, 0.23854862, 0.14238955],
[0.1483039 , 0.25893564, 0.26329003, 0.32947044, 0.32947044,
0.26329003, 0.25893564, 0.1483039 ],
[0.17505961, 0.24311495, 0.28587807, 0.29594737, 0.29594737,
0.28587807, 0.24311495, 0.17505961],
[0.19992018, 0.22273316, 0.2967273 , 0.28061936, 0.28061936,
0.2967273 , 0.22273316, 0.19992018],
[0.20791167, 0.22045732, 0.27844606, 0.29318495, 0.29318495,
0.27844606, 0.22045732, 0.20791167],
[0.20652158, 0.2223829 , 0.25110829, 0.31998723, 0.31998723,
0.25110829, 0.2223829 , 0.20652158],
[0.19739406, 0.22221787, 0.24446542, 0.33592265, 0.33592265,
0.24446542, 0.22221787, 0.19739406],
[0.17154275, 0.23700171, 0.25605553, 0.33540001, 0.33540001,
0.25605553, 0.23700171, 0.17154275],
[0.13515694, 0.25676317, 0.28243583, 0.32564406, 0.32564406,
0.28243583, 0.25676317, 0.13515694],
[0.10824455, 0.25131245, 0.3312888 , 0.30915419, 0.30915419,
0.3312888 , 0.25131245, 0.10824455],
[0.09862348, 0.22959338, 0.37490813, 0.296875 , 0.296875 ,
0.37490813, 0.22959338, 0.09862348],
[0.09924836, 0.22161757, 0.37235462, 0.30677945, 0.30677945,
0.37235462, 0.22161757, 0.09924836],
[0.10680358, 0.22666258, 0.33976988, 0.32676396, 0.32676396,
0.33976988, 0.22666258, 0.10680358],
[0.12427806, 0.23230582, 0.3160371 , 0.32737902, 0.32737902,
0.3160371 , 0.23230582, 0.12427806],
[0.14753597, 0.24368196, 0.29319855, 0.31558352, 0.31558352,
0.29319855, 0.24368196, 0.14753597],
[0.16914123, 0.26366122, 0.25601447, 0.31118307, 0.31118307,
0.25601447, 0.26366122, 0.16914123],
[0.19482008, 0.27530047, 0.23387412, 0.29600532, 0.29600532,
0.23387412, 0.27530047, 0.19482008],
[0.2353433 , 0.26007741, 0.24609469, 0.2584846 , 0.2584846 ,
0.24609469, 0.26007741, 0.2353433 ],
[0.28573776, 0.22556153, 0.25776838, 0.23093234, 0.23093234,
0.25776838, 0.22556153, 0.28573776],
[0.33071472, 0.19715894, 0.24110484, 0.2310215 , 0.2310215 ,
0.24110484, 0.19715894, 0.33071472]]) In [152]:
fig, ax = plt.subplots(1, 2, sharey=True)
ax[0].imshow(nup_all_time, cmap='Greens', aspect='auto', origin='lower')
ax[1].imshow(ndn_all_time, cmap='Blues', aspect='auto', origin='lower')
Out[152]:
<matplotlib.image.AxesImage at 0x7fbf8225b2e0>
In [99]:
spread= np.abs(np.arange(1, params.nsites+1)-(params.nsites+1)/2)
charge_spread=[np.multiply(spread,np.diagonal(alpha+beta).real).sum()
for (alpha,beta) in zip(opdms_alpha, opdms_beta)]
spin_spread=[np.multiply(spread,np.diagonal(alpha-beta).real).sum()
for (alpha, beta) in zip(opdms_alpha, opdms_beta)]
In [137]:
fig,ax=plt.subplots(2,1, sharex=True, figsize=(4,6))
stop=12
ax[0].plot(times[:stop], charge_spread[:stop], 'o-', mfc='none', color='tab:blue', label='charge')
ax2=ax[0].twinx()
ax[0].set_ylim(5, 10)
ax[0].plot(times[:stop], spin_spread[:stop], 'o-', mfc='none', color='tab:red', label='spin')
#ax[1].set_xticks(list(range(0, 3, 1)))
ax2.plot(times[:stop], spin_spread[:stop], 'o-', mfc='none', color='tab:red', label='spin')
#ax[1].set_xticks(list(range(0, 3, 1)))
ax[-1].set_xlabel(r'Time ($\hbar/t$)')
ax[0].set_ylabel(r'Charge spread $\kappa^{+}$'); ax2.set_ylabel(r'Spin spread $\kappa^{-}$')
ax2.set_yticks([-2,0,2])
ax[1].plot(times[:stop], np.gradient(charge_spread[:stop]), 'o-', mfc='none', label='charge')
#ax2=ax[1].twinx()
ax[1].plot(times[:stop], np.gradient(spin_spread[:stop]), 'o-', mfc='none', color='tab:red', label='spin')
ax[1].set_ylabel(r'Spread rate $(t/\hbar)$')
plt.subplots_adjust(hspace=0.3)
ax[1].legend(); ax[0].legend()
ax[1].set_ylim(-1, 2)
plt.savefig('spread.pdf')
