Python numpy实现Hatree-Fork自洽计算程序

还是numpy香，几行代码就实现了向量化的Hatree-Fork计算过程

# 生成块三对角矩阵函数, 工具函数
def tridiag(c, u, d, N): 
    # c, u, d are center, upper and lower blocks, repeat N times
    cc = block_diag(*([c]*N)) 
    shift = c.shape[1]
    uu = block_diag(*([u]*N)) 
    uu = np.hstack((np.zeros((uu.shape[0], shift)), uu[:,:-shift]))
    dd = block_diag(*([d]*N)) 
    dd = np.hstack((dd[:,shift:],np.zeros((uu.shape[0], shift))))
    return cc+uu+dd

H0 

# Hamiltonian_H0函数
def Hamiltonian_H0(k,N,t=-1):
    """
    t 是最近邻hopping系数
    k 是Bloch波矢量
    N 是宽度
    """
    T = t*np.matrix([[0,1+np.exp(-1j*k)],[1+np.exp(1j*k),0]])
    D = t*np.matrix([[0,0],[1,0]])
    H = tridiag(T,D.H,D,N)
    return H

向量化的费米函数 

# 费米函数
def fermi(E, mu, T):
    return 1/(np.exp((E-mu)/T)+1) 

Hatree-Fork过程

# Hatree_Fork 自洽计算过程
def Hatree_Fock(ks, Nup_avg, Ndown_avg, U, mu,T, N, ncc):
    m = 2*N
    nk = ks.size
    for i 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 = Hamiltonian_H0(ks[i],N)
            Hk_up = Hk0 + U*(Ndown_avg -0.5*np.eye(m))
            Hk_down = Hk0 + U*(Nup_avg - 0.5*np.eye(m))
            Ek_up, Ak_up = LA.eigh(Hk_up)
            Ek_down, Ak_down = LA.eigh(Hk_down)
            fermi_up = np.matlib.repmat(fermi(Ek_up,mu,T),m,1)
            fermi_down = np.matlib.repmat(fermi(Ek_down,mu,T),m,1)
            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 i%5==0:
            print("已经完成%d次迭代"%i)
        
        return Nup_avg, Ndown_avg

能带计算过程

nk = 1024
ks = np.linspace(0,2*pi,nk)
U = 2
mu = 0
T = 0.001
N = 64
m = 2*N
ncc = 20
Nup_avg = 0.5*np.eye(m)
Ndown_avg = 0.5*np.eye(m)
Nup_avg[0,0] = 0.7
Nup_avg[-1,-1] = 0.7
Ndown_avg[0,0] = 0.3
Ndown_avg[-1,-1] = 0.3
Nup_avg, Ndown_avg = Hatree_Fock(ks,Nup_avg,Ndown_avg,U,mu,T,N,ncc)
band_up = np.zeros((nk,m))
band_down = np.zeros((nk,m))
for i in range(nk):
    Hk0 = Hamiltonian_H0(ks[i],N)
    Hk_up = Hk0 + U*(Ndown_avg -0.5*np.eye(m))
    Hk_down = Hk0 + U*(Nup_avg - 0.5*np.eye(m))
    Ek_up, _ = LA.eigh(Hk_up)
    Ek_down, _ = LA.eigh(Hk_down)
    band_up[i,:] = Ek_up
    band_down[i,:] = Ek_down
for i in range(m):
    plt.plot(ks,band_up[:,i],'r')
    plt.plot(ks,band_down[:,i],'b')
plt.ylim((-1,1))
plt.show()

计算谱函数

nw = 400
delta = 0.002
nk = 400
ks = np.linspace(0,2*pi,nk)
omega = np.linspace(-1,1,nw) 
Ak_up = np.zeros((nw,nk))
Ak_down = np.zeros((nw,nk))
N = 32
for i in range(nw):
    for j in range(nk):
        Hk0 = Hamiltonian_H0(ks[j],N)
        Hk_up = Hk0 + U*(Ndown_avg -0.5*np.eye(m))
        Hk_down = Hk0 + U*(Nup_avg - 0.5*np.eye(m))
        Gk_up = inv((omega[i]+1j*delta)*np.eye(m)-Hk_up)
        Gk_down = inv((omega[i]+1j*delta)*np.eye(m)-Hk_down)
        Ak_up[i,j] = np.imag(Gk_up.trace())
        Ak_down[i,j] = np.imag(Gk_down.trace())
X,Y = np.meshgrid(ks,omega)
plt.pcolormesh(X, Y, Ak_up)

接下来就是实现双粒子格林函数的计算