Hubbard模型平均场demo代码

import numpy as np
import matplotlib.pyplot as plt
from numpy import cos, pi, exp
from numpy.linalg import eigh


t = 1
la = 1.4
e = t * (la**2 - 2)
nk = 256
delta = 0.01
U = 3
mu = 0
T = 0.01


def Hamiltonian0(k):
    H = np.zeros((2, 2))
    H[0, 0] = 2 * t * cos(k) + 2 * t
    H[0, 1] = 2 * t * la * cos(k/2)
    H[1, 0] = 2 * t * la * cos(k/2)
    H[1, 1] = e  + delta + 2 * t
    return H

def Fermi(E, mu, T):
    return 1/(exp((E - mu)/T) + 1)



# Hatree_Fock 自洽计算过程
def Hatree_Fock(ks, Nup_avg, Ndown_avg, mu, ncc = 200):
    m = 2
    nk = ks.size
    for nc in range(ncc):
        Nup_avg_tmp = np.zeros((m,m), dtype='double')
        Ndown_avg_tmp = np.zeros((m,m), dtype='double')
        for i in range(nk):
            Hk0 = Hamiltonian0(ks[i])
            Hk_up = Hk0 + U*(Ndown_avg - Nup_avg * Ndown_avg)
            Hk_down = Hk0 + U*(Nup_avg - Nup_avg * Ndown_avg)
            Ek_up, Ak_up = eigh(Hk_up)
            Ek_down, Ak_down = eigh(Hk_down)
            fermi_up = Fermi(Ek_up,mu, T).reshape(1,m)
            fermi_down = Fermi(Ek_down,mu,T).reshape(1,m)
            Nk_up = (Ak_up*Ak_up.conj()*fermi_up).sum(axis=1,dtype='double')
            Nk_down = (Ak_down*Ak_down.conj()*fermi_down).sum(axis=1,dtype='double')
            Nup_avg_tmp += np.diag(Nk_up)
            Ndown_avg_tmp += np.diag(Nk_down)
        Nup_avg = Nup_avg_tmp/nk
        Ndown_avg = Ndown_avg_tmp/nk
        if nc%5==0:
            print(f"已经完成{nc}次迭代")
        
    return Nup_avg, Ndown_avg


def plot_band(Nup_avg, Ndown_avg):
    ks = np.linspace(0, 2*pi, nk)
    Ek_up = np.zeros((len(ks), 2))
    Ek_down = np.zeros((len(ks), 2))
    for i in range(len(ks)):
        Hk0 = Hamiltonian0(ks[i])
        Ek_up[i, :], _ = eigh(Hk0 + U*(Ndown_avg - Nup_avg * Ndown_avg))
        Ek_down[i, :], _ = eigh(Hk0 + U*(Nup_avg - Nup_avg * Ndown_avg))
    plt.plot(ks, Ek_up, 'r')
    plt.plot(ks, Ek_down, 'b')
    plt.show()
    

def Hatree_Fork_band():
    Nup_avg = 0.25 * np.eye(2)
    Ndown_avg = np.zeros((2, 2))
    ks = np.linspace(0, 2*pi, nk)
    Nup_avg, Ndown_avg = Hatree_Fock(ks, Nup_avg, Ndown_avg, mu)
    print(Nup_avg)
    print(Ndown_avg)
    print(np.trace(Nup_avg + Ndown_avg))
    plot_band(Nup_avg, Ndown_avg)

Hatree_Fork_band()