hubbard_1site.py

hubbard_1site.py

from __future__ import print_function

import numpy as np


def make_local_ops():
    ops = {}

    # c^+
    cdag = np.zeros((2,2))
    cdag[1, 0] = 1
    ops['c^+'] = cdag

    # c
    c = cdag.transpose()
    ops['c'] = c

    # I
    ops['I'] = np.identity(2)

    # F (for fermionic anticommutation)
    ops['F'] = np.diag((1.0, -1.0))

    return ops


def make_matrix(*ops):
    r = 1.0
    for op in ops[::-1]:
        r = np.kron(op, r)
    return r


def hermite_conjugate(mat):
    return mat.conj().T


def solve_hubbard_1site(U, mu=0):

    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, cdag)  # F represents fermionic anticommutation
    Cdag['1d'] = make_matrix(cdag, F)

    C = {}  # annihilation operators
    N = {}  # number operators
    for key, cdag in Cdag.items():
        C[key] = hermite_conjugate(cdag)
        N[key] = cdag @ C[key]

    hamil = U * N['1u'] @ N['1d'] - mu * (N['1u'] + N['1d'])
    print("H =\n", hamil)

    eigval, eigvec = np.linalg.eigh(hamil)

    return eigval, eigvec


def main():
    U = 10.0
    mu = 2

    E, vec = solve_hubbard_1site(U, mu)

    print("\nE =\n", E)

    print("\nEigenvectors =")
    for i in range(vec.shape[1]):
        print(vec[:, i])


if __name__ == '__main__':
    main()