严格投影对角化优化后的代码版本

import numpy as np
import numpy.matlib
from scipy.linalg import block_diag
from numpy import linalg as LA
import matplotlib.pyplot as plt
from math import sqrt,pi

生成三对角矩阵

# 生成块三对角矩阵函数
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

def Hamiltonian_H0_SU4(k,N,t=-1):
    """
    t 是最近邻hopping系数
    k 是Bloch波矢量
    N 是宽度
    """
    A = np.matrix([[sqrt(3),0,1,np.exp(-1j*k)],[0, sqrt(3),1,1],[1,1,sqrt(3),0],[np.exp(1j*k),1,0,sqrt(3)]])
    B = np.matrix([[0,0,0,0],[0,0,0,0],[1,0,0,0],[0,-1,0,0]])
    H = tridiag(A,B,B.H,N)
    H[0,0] = 1000
    return H

对角化

nk = 700
N = 128
k = np.linspace(0,2*pi,nk)
band = np.zeros((4*N, nk))
Ak = np.zeros((4*N,nk),dtype="complex")
for i in range(nk):
    Hk0 = Hamiltonian_H0_SU4(k[i],N)
    E, A = LA.eigh(Hk0)
    band[:,i] = E
    Ak[:,i] = A[:,N-1]

画图

for i in range(4*N):
    plt.plot(k,band[i,:])
plt.ylim((-1,1))

本征向量

Ak = Ak[:8,:]
Ap = np.power(np.absolute(Ak),2).sum(axis=1).copy()
plt.plot(Ap)

相位因子矩阵

Q = np.zeros((8,nk),dtype="complex")
Q[0,:] = 1
Q[1,:] = np.exp(sqrt(3)/2*1j)
Q[2,:] = np.exp(sqrt(3)*1j)
Q[3,:] = np.exp(sqrt(3)*3/2*1j)
Q[4,:] = 1
Q[5,:] = np.exp(sqrt(3)/2*1j)
Q[6,:] = np.exp(sqrt(3)*1j)
Q[7,:] = np.exp(sqrt(3)*3/2*1j)

生成有效哈密顿量

def Hamiltonian_Heff(nk,iq,Ak,U,N):
    H = np.zeros((nk,nk),dtype='complex')
    for j in range(nk):
        j_q = (j - iq + nk)%nk
        for i in range(j+1,nk):
            i_q = (i-iq + nk)%nk
            H[j,i] -= np.sum(Ak[:,i_q].conj()*Ak[:,i]*Ak[:,j_q]*Ak[:,j].conj())
    
    H = H + H.T
    
    for i in range(nk):
        i_q = (i-iq + nk)%nk
        H[i,i] -= np.sum(Ak[:,i_q].conj()*Ak[:,i]*Ak[:,i_q]*Ak[:,i].conj())
    
    M = H.copy()
    
    
    for i in range(nk):
        i_q = (i - iq + nk)%nk
        H[i,i] += np.sum(Ap*np.power(np.absolute(Ak[:,i]),2))
    
    
    
    return U*H/N, U*M/N

对角化

U = 1.0
band_mag = np.zeros((nk, nk))
band_M = np.zeros((nk,nk))
psi = np.zeros((nk,nk,nk),dtype='complex')
psi_M = np.zeros((nk,nk,nk),dtype='complex')
for i in range(nk):
    Heff,M = Hamiltonian_Heff(nk,i,Ak,U,N)
    E, A = LA.eigh(Heff)
    Em, Am = LA.eigh(M)
    psi[:,:,i] = A 
    band_mag[:,i] = E
    band_M[:,i] = Em
    psi_M[:,:,i] = Am

画图

for i in range(nk):
    plt.plot(k,band_mag[i,:],color='gray')
plt.ylim((0,0.6))

画

for i in range(2):
    plt.plot(k,band_M[i,:],color='gray')

辅助函数

def circshift(matrix,shiftnum1,shiftnum2):
    """
    """
    h,w=matrix.shape
    matrix=np.vstack((matrix[(h-shiftnum1):,:],matrix[:(h-shiftnum1),:]))
    matrix=np.hstack((matrix[:,(w-shiftnum2):],matrix[:,:(w-shiftnum2)]))
    return matrix

谱函数

delta = 0.002
nw = 128
w = np.linspace(0,0.7,nw)
Spectral_A = np.zeros((nw,nk),dtype='complex')
for i in range(nk):
    for j in range(nk):
        A_kq = circshift(Ak,0,j)
        S_plus = np.sum(Ak.conj()*A_kq*Q,axis=0).reshape(nk,1)
        P = np.absolute(np.dot(S_plus.T,psi[:,i,j]))**2
        #P = np.absolute(np.dot(psi_M[:,1,j].T,psi[:,i,j]))**2
        for p in range(nw):
            Spectral_A[p,j] += P/(w[p]-band_mag[i,j]+1j*delta)

画图

X,Y = np.meshgrid(k,w)
plt.figure(figsize=(5,5))
Z = -Spectral_A.imag/pi
Z = np.log(Z)
plt.pcolormesh(X, Y, Z)
plt.colorbar()
#plt.ylim(0.02, 0.05)

严格投影对角化优化后的代码版本

猜你喜欢