from __future__ import print_function
import numpy as np
from hubbard_1site import make_local_ops, make_matrix, hermite_conjugate
def projection(h, n):
from itertools import product
# occup = [sum(state) for state in product([0, 1], repeat=4)]
# index = [i for i, occ in enumerate(occup) if occ == n]
index = [i for i, state in enumerate(product([0, 1], repeat=4)) if sum(state)==n]
return h[index, :][:, index]
def solve_hubbard_2site(t, U, mu=0, n=None):
local_ops = make_local_ops()
cdag = local_ops['c^+']
I = local_ops['I']
F = local_ops['F']
# creation operators
Cdag = {}
Cdag['1u'] = make_matrix(I, I, I, cdag)
Cdag['1d'] = make_matrix(I, I, cdag, F)
Cdag['2u'] = make_matrix(I, cdag, F, F)
Cdag['2d'] = make_matrix(cdag, F, F, F)
C = {} # annihilation operators
N = {} # number operators
for key, cdag in Cdag.items():
C[key] = hermite_conjugate(cdag)
N[key] = cdag @ C[key]
hamil = 0
# t
for key1, key2 in [('1u', '2u'), ('1d', '2d'), ('2u', '1u'), ('2d', '1d')]:
hamil += -t * Cdag[key1] @ C[key2]
# U
for key1, key2 in [('1u', '1d'), ('2u', '2d')]:
hamil += U * N[key1] @ N[key2]
# mu
for n_op in N.values():
hamil += -mu * n_op
# print("H =\n", hamil)
if n is not None:
# projection to n-particle state
hamil = projection(hamil, n)
eigval, eigvec = np.linalg.eigh(hamil)
return eigval, eigvec
def main():
t = 1.0
U = 4.0
mu = 2
E, vec = solve_hubbard_2site(t, U, mu)
# E, vec = solve_hubbard_2site(t, U, mu, n=2)
print("\nE =\n", E)
# print("\nEigenvectors =")
# for i in range(vec.shape[1]):
# print(vec[:, i])
if __name__ == '__main__':
main()