From 69b445d625abb711f6dd806f171cb9e6babc1182 Mon Sep 17 00:00:00 2001 From: guanjihuan Date: Fri, 15 Jul 2022 14:25:27 +0800 Subject: [PATCH] update --- ...lculation of spectral function and QPI.f90 | 318 +++++++++--------- ...alculation of spectral function and QPI.py | 316 ++++++++--------- 2 files changed, 317 insertions(+), 317 deletions(-) rename academic_codes/{2019.12.30_calculation of spectral function and QPI => 2019.12.30_calculation_of_spectral_function_and_QPI}/calculation of spectral function and QPI.f90 (97%) mode change 100755 => 100644 rename academic_codes/{2019.12.30_calculation of spectral function and QPI => 2019.12.30_calculation_of_spectral_function_and_QPI}/calculation of spectral function and QPI.py (97%) mode change 100755 => 100644 diff --git a/academic_codes/2019.12.30_calculation of spectral function and QPI/calculation of spectral function and QPI.f90 b/academic_codes/2019.12.30_calculation_of_spectral_function_and_QPI/calculation of spectral function and QPI.f90 old mode 100755 new mode 100644 similarity index 97% rename from academic_codes/2019.12.30_calculation of spectral function and QPI/calculation of spectral function and QPI.f90 rename to academic_codes/2019.12.30_calculation_of_spectral_function_and_QPI/calculation of spectral function and QPI.f90 index 5c6ff86..19bc974 --- a/academic_codes/2019.12.30_calculation of spectral function and QPI/calculation of spectral function and QPI.f90 +++ b/academic_codes/2019.12.30_calculation_of_spectral_function_and_QPI/calculation of spectral function and QPI.f90 @@ -1,160 +1,160 @@ -! This code is supported by the website: https://www.guanjihuan.com -! The newest version of this code is on the web page: https://www.guanjihuan.com/archives/3785 - -module global - implicit none - double precision sqrt3,Pi - parameter(sqrt3=1.7320508075688773d0,Pi=3.14159265358979324d0) -end module global - - - -program QPI !QPI主程序 - use blas95 - use lapack95,only:GETRF,GETRI - use global - implicit none - integer i,j,info,index_0(4) - double precision omega,kx,ky,Eigenvalues(4),eta,V0,kx1,kx2,ky1,ky2,qx,qy,time_begin,time_end - parameter(eta=0.005) - complex*16 H0(4,4),green_0(4,4),green_1(4,4),green_0_k1(4,4),green_0_k2(4,4),A_spectral,V(4,4),gamma_0(4,4),Temp_0(4,4),T(4,4),g_1,rho_1 - character(len=*):: Flname - parameter(Flname='') !可以写上输出文件路径，也可以不写，输出存在当前文件的路径 - - omega=0.070d0 - open(unit=10,file=Flname//'Spectral function_w=0.07.txt') - open(unit=20,file=Flname//'QPI_intra_nonmag_w=0.07.txt') - call CPU_TIME(time_begin) - - !计算谱函数A(kx,ky) - write(10,"(f20.10,x)",advance='no') 0 - do ky=-Pi,Pi,0.01d0 !谱函数图案的精度 - write(10,"(f20.10,x)",advance='no') ky - enddo - write(10,"(a)",advance='yes') ' ' - do kx=-Pi,Pi,0.01d0 !谱函数图案的精度 - write(10,"(f20.10,x)",advance='no') kx - do ky=-Pi,Pi,0.01d0 !谱函数图案的精度 - call Greenfunction_clean(kx,ky,eta,omega,green_0) - A_spectral=-(green_0(1,1)+green_0(3,3))/Pi - write(10,"(f20.10)",advance='no') imag(A_spectral) - enddo - write(10,"(a)",advance='yes') ' ' - enddo - - !计算QPI(qx,qy) - V0=0.4d0 - V=0.d0 - V(1,1)=V0 - V(2,2)=-V0 - V(3,3)=V0 - V(4,4)=-V0 - gamma_0=0.d0 - do kx=-Pi,Pi,0.01 - do ky=-Pi,Pi,0.01 - call Greenfunction_clean(kx,ky,eta,omega,green_0) - do i=1,4 - do j=1,4 - gamma_0(i,j)=gamma_0(i,j)+green_0(i,j)*0.01*0.01 - enddo - enddo - enddo - enddo - gamma_0=gamma_0/(2*Pi)/(2*Pi) - call gemm(V,gamma_0,Temp_0) - do i=1,4 - Temp_0(i,i)=1-Temp_0(i,i) - enddo - call GETRF( Temp_0,index_0,info ); call GETRI( Temp_0,index_0,info) !求逆 - call gemm(Temp_0,V,T) !矩阵乘积 - write(20,"(f20.10,x)",advance='no') 0 - do qy=-Pi,Pi,0.01 !QPI图案的精度 - write(20,"(f20.10,x)",advance='no') qy - enddo - write(20,"(a)",advance='yes') ' ' - do qx=-Pi,Pi,0.01 !QPI图案的精度 - write(*,"(a)",advance='no') 'qx=' - write(*,*) qx !屏幕输出可以实时查看计算进度 - write(20,"(f20.10)",advance='no') qx - do qy=-Pi,Pi,0.01 !QPI图案的精度 - rho_1=0.d0 - do kx1=-Pi,Pi,0.06 !积分的精度 - kx2=kx1+qx - do ky1=-Pi,Pi,0.06 !积分的精度 - ky2=ky1+qy - call Greenfunction_clean(kx1,ky1,eta,omega,green_0_k1) - call Greenfunction_clean(kx2,ky2,eta,omega,green_0_k2) - call gemm(green_0_k1,T,Temp_0) - call gemm(Temp_0, green_0_k2, green_1) - g_1=green_1(1,1)-dconjg(green_1(1,1))+green_1(3,3)-dconjg(green_1(3,3)) - rho_1=rho_1+g_1*0.06*0.06 - enddo - enddo - rho_1=rho_1/(2*Pi)/(2*Pi)/(2*Pi)*(0.d0,1.d0) - write(20,"(f20.10,x,f20.10)",advance='no') real(rho_1) - enddo - write(20,"(a)",advance='yes') ' ' - enddo - - call CPU_TIME(time_end) - write(*,"(a)",advance='no') 'The running time of this task=' - write (*,*) time_end-time_begin !屏幕输出总的计算时间，单位为秒（按照当前步长的精度，在个人计算机上运算大概需要4个小时） -end program - - - -subroutine Greenfunction_clean(kx,ky,eta,omega,green_0) !干净体系的格林函数 -use blas95 -use lapack95,only:GETRF,GETRI -use global -integer info,index_0(4) -double precision, intent(in):: kx,ky,eta,omega -complex*16 H0(4,4) -complex*16,intent(out):: green_0(4,4) -call Hamiltonian(kx,ky,H0) -green_0=H0 -do i=1,4 - green_0(i,i)=omega+(0.d0,1.d0)*eta-green_0(i,i) -enddo -call GETRF( green_0,index_0,info ); call GETRI( green_0,index_0,info ); -end subroutine Greenfunction_clean - - - -subroutine Hamiltonian(kx,ky,Matrix) !哈密顿量 -use global -implicit none -integer i,j -double precision t1,t2,t3,t4,mu,epsilon_x,epsilon_y,epsilon_xy,delta_1,delta_2,delta_0 -double precision, intent(in):: kx,ky -complex*16,intent(out):: Matrix(4,4) - -t1=-1;t2=1.3;t3=-0.85;t4=-0.85;delta_0=0.1;mu=1.54 -Matrix=(0.d0,0.d0) - -epsilon_x=-2*t1*dcos(kx)-2*t2*dcos(ky)-4*t3*dcos(kx)*dcos(ky) -epsilon_y=-2*t1*dcos(ky)-2*t2*dcos(kx)-4*t3*dcos(kx)*dcos(ky) -epsilon_xy=-4*t4*dsin(kx)*dsin(ky) -delta_1=delta_0*dcos(kx)*dcos(ky) -delta_2=delta_1 - -Matrix(1,1)=epsilon_x-mu -Matrix(2,2)=-epsilon_x+mu -Matrix(3,3)=epsilon_y-mu -Matrix(4,4)=-epsilon_y+mu - -Matrix(1,2)=delta_1 -Matrix(2,1)=delta_1 -Matrix(1,3)=epsilon_xy -Matrix(3,1)=epsilon_xy -Matrix(1,4)=0.d0 -Matrix(4,1)=0.d0 - -Matrix(2,3)=0.d0 -Matrix(3,2)=0.d0 -Matrix(2,4)=-epsilon_xy -Matrix(4,2)=-epsilon_xy - -Matrix(3,4)=delta_2 -Matrix(4,3)=delta_2 +! This code is supported by the website: https://www.guanjihuan.com +! The newest version of this code is on the web page: https://www.guanjihuan.com/archives/3785 + +module global + implicit none + double precision sqrt3,Pi + parameter(sqrt3=1.7320508075688773d0,Pi=3.14159265358979324d0) +end module global + + + +program QPI !QPI主程序 + use blas95 + use lapack95,only:GETRF,GETRI + use global + implicit none + integer i,j,info,index_0(4) + double precision omega,kx,ky,Eigenvalues(4),eta,V0,kx1,kx2,ky1,ky2,qx,qy,time_begin,time_end + parameter(eta=0.005) + complex*16 H0(4,4),green_0(4,4),green_1(4,4),green_0_k1(4,4),green_0_k2(4,4),A_spectral,V(4,4),gamma_0(4,4),Temp_0(4,4),T(4,4),g_1,rho_1 + character(len=*):: Flname + parameter(Flname='') !可以写上输出文件路径，也可以不写，输出存在当前文件的路径 + + omega=0.070d0 + open(unit=10,file=Flname//'Spectral function_w=0.07.txt') + open(unit=20,file=Flname//'QPI_intra_nonmag_w=0.07.txt') + call CPU_TIME(time_begin) + + !计算谱函数A(kx,ky) + write(10,"(f20.10,x)",advance='no') 0 + do ky=-Pi,Pi,0.01d0 !谱函数图案的精度 + write(10,"(f20.10,x)",advance='no') ky + enddo + write(10,"(a)",advance='yes') ' ' + do kx=-Pi,Pi,0.01d0 !谱函数图案的精度 + write(10,"(f20.10,x)",advance='no') kx + do ky=-Pi,Pi,0.01d0 !谱函数图案的精度 + call Greenfunction_clean(kx,ky,eta,omega,green_0) + A_spectral=-(green_0(1,1)+green_0(3,3))/Pi + write(10,"(f20.10)",advance='no') imag(A_spectral) + enddo + write(10,"(a)",advance='yes') ' ' + enddo + + !计算QPI(qx,qy) + V0=0.4d0 + V=0.d0 + V(1,1)=V0 + V(2,2)=-V0 + V(3,3)=V0 + V(4,4)=-V0 + gamma_0=0.d0 + do kx=-Pi,Pi,0.01 + do ky=-Pi,Pi,0.01 + call Greenfunction_clean(kx,ky,eta,omega,green_0) + do i=1,4 + do j=1,4 + gamma_0(i,j)=gamma_0(i,j)+green_0(i,j)*0.01*0.01 + enddo + enddo + enddo + enddo + gamma_0=gamma_0/(2*Pi)/(2*Pi) + call gemm(V,gamma_0,Temp_0) + do i=1,4 + Temp_0(i,i)=1-Temp_0(i,i) + enddo + call GETRF( Temp_0,index_0,info ); call GETRI( Temp_0,index_0,info) !求逆 + call gemm(Temp_0,V,T) !矩阵乘积 + write(20,"(f20.10,x)",advance='no') 0 + do qy=-Pi,Pi,0.01 !QPI图案的精度 + write(20,"(f20.10,x)",advance='no') qy + enddo + write(20,"(a)",advance='yes') ' ' + do qx=-Pi,Pi,0.01 !QPI图案的精度 + write(*,"(a)",advance='no') 'qx=' + write(*,*) qx !屏幕输出可以实时查看计算进度 + write(20,"(f20.10)",advance='no') qx + do qy=-Pi,Pi,0.01 !QPI图案的精度 + rho_1=0.d0 + do kx1=-Pi,Pi,0.06 !积分的精度 + kx2=kx1+qx + do ky1=-Pi,Pi,0.06 !积分的精度 + ky2=ky1+qy + call Greenfunction_clean(kx1,ky1,eta,omega,green_0_k1) + call Greenfunction_clean(kx2,ky2,eta,omega,green_0_k2) + call gemm(green_0_k1,T,Temp_0) + call gemm(Temp_0, green_0_k2, green_1) + g_1=green_1(1,1)-dconjg(green_1(1,1))+green_1(3,3)-dconjg(green_1(3,3)) + rho_1=rho_1+g_1*0.06*0.06 + enddo + enddo + rho_1=rho_1/(2*Pi)/(2*Pi)/(2*Pi)*(0.d0,1.d0) + write(20,"(f20.10,x,f20.10)",advance='no') real(rho_1) + enddo + write(20,"(a)",advance='yes') ' ' + enddo + + call CPU_TIME(time_end) + write(*,"(a)",advance='no') 'The running time of this task=' + write (*,*) time_end-time_begin !屏幕输出总的计算时间，单位为秒（按照当前步长的精度，在个人计算机上运算大概需要4个小时） +end program + + + +subroutine Greenfunction_clean(kx,ky,eta,omega,green_0) !干净体系的格林函数 +use blas95 +use lapack95,only:GETRF,GETRI +use global +integer info,index_0(4) +double precision, intent(in):: kx,ky,eta,omega +complex*16 H0(4,4) +complex*16,intent(out):: green_0(4,4) +call Hamiltonian(kx,ky,H0) +green_0=H0 +do i=1,4 + green_0(i,i)=omega+(0.d0,1.d0)*eta-green_0(i,i) +enddo +call GETRF( green_0,index_0,info ); call GETRI( green_0,index_0,info ); +end subroutine Greenfunction_clean + + + +subroutine Hamiltonian(kx,ky,Matrix) !哈密顿量 +use global +implicit none +integer i,j +double precision t1,t2,t3,t4,mu,epsilon_x,epsilon_y,epsilon_xy,delta_1,delta_2,delta_0 +double precision, intent(in):: kx,ky +complex*16,intent(out):: Matrix(4,4) + +t1=-1;t2=1.3;t3=-0.85;t4=-0.85;delta_0=0.1;mu=1.54 +Matrix=(0.d0,0.d0) + +epsilon_x=-2*t1*dcos(kx)-2*t2*dcos(ky)-4*t3*dcos(kx)*dcos(ky) +epsilon_y=-2*t1*dcos(ky)-2*t2*dcos(kx)-4*t3*dcos(kx)*dcos(ky) +epsilon_xy=-4*t4*dsin(kx)*dsin(ky) +delta_1=delta_0*dcos(kx)*dcos(ky) +delta_2=delta_1 + +Matrix(1,1)=epsilon_x-mu +Matrix(2,2)=-epsilon_x+mu +Matrix(3,3)=epsilon_y-mu +Matrix(4,4)=-epsilon_y+mu + +Matrix(1,2)=delta_1 +Matrix(2,1)=delta_1 +Matrix(1,3)=epsilon_xy +Matrix(3,1)=epsilon_xy +Matrix(1,4)=0.d0 +Matrix(4,1)=0.d0 + +Matrix(2,3)=0.d0 +Matrix(3,2)=0.d0 +Matrix(2,4)=-epsilon_xy +Matrix(4,2)=-epsilon_xy + +Matrix(3,4)=delta_2 +Matrix(4,3)=delta_2 end subroutine Hamiltonian \ No newline at end of file diff --git a/academic_codes/2019.12.30_calculation of spectral function and QPI/calculation of spectral function and QPI.py b/academic_codes/2019.12.30_calculation_of_spectral_function_and_QPI/calculation of spectral function and QPI.py old mode 100755 new mode 100644 similarity index 97% rename from academic_codes/2019.12.30_calculation of spectral function and QPI/calculation of spectral function and QPI.py rename to academic_codes/2019.12.30_calculation_of_spectral_function_and_QPI/calculation of spectral function and QPI.py index 12fdd64..b95151f --- a/academic_codes/2019.12.30_calculation of spectral function and QPI/calculation of spectral function and QPI.py +++ b/academic_codes/2019.12.30_calculation_of_spectral_function_and_QPI/calculation of spectral function and QPI.py @@ -1,158 +1,158 @@ -""" -This code is supported by the website: https://www.guanjihuan.com -The newest version of this code is on the web page: https://www.guanjihuan.com/archives/3785 -""" - -import numpy as np -from math import * -import matplotlib.pyplot as plt -from matplotlib.colors import ListedColormap -import time - - -def green_function(fermi_energy, k1, k2, hamiltonian): # 计算格林函数 - matrix0 = hamiltonian(k1, k2) - dim = matrix0.shape[0] - green = np.linalg.inv(fermi_energy * np.identity(dim) - matrix0) - return green - - -def spectral_function(fermi_energy, k1, k2, hamiltonian): # 计算谱函数 - dim1 = k1.shape[0] - dim2 = k2.shape[0] - spectrum = np.zeros((dim1, dim2)) - i0 = 0 - for k10 in k1: - j0 = 0 - for k20 in k2: - green = green_function(fermi_energy, k10, k20, hamiltonian) - spectrum[i0, j0] = (np.imag(green[0,0])+np.imag(green[2,2]))/(-pi) - j0 += 1 - i0 += 1 - # print(spectrum) - print() - print('Spectral function显示的网格点 =', k1.shape[0], '*', k1.shape[0], '; 步长 =', k1[1] - k1[0]) - print() - return spectrum - - -def qpi(fermi_energy, q1, q2, hamiltonian, potential_i): # 计算QPI - dim = hamiltonian(0, 0).shape[0] - ki1 = np.arange(-pi, pi, 0.01) # 计算gamma_0时，k的积分密度 - ki2 = np.arange(-pi, pi, 0.01) - print('gamma_0的积分网格点 =', ki1.shape[0], '*', ki1.shape[0], '; 步长 =', ki1[1] - ki1[0]) - gamma_0 = integral_of_green(fermi_energy, ki1, ki2, hamiltonian)/np.square(2*pi) - t_matrix = np.dot(np.linalg.inv(np.identity(dim)-np.dot(potential_i, gamma_0)), potential_i) - ki1 = np.arange(-pi, pi, 0.06) # 计算induced_local_density时，k的积分密度 - ki2 = np.arange(-pi, pi, 0.06) - print('局域态密度变化的积分网格点 =', ki1.shape[0], '*', ki1.shape[0], '; 步长 =', ki1[1] - ki1[0]) - print('QPI显示的网格点 =', q1.shape[0], '*', q1.shape[0], '; 步长 =', q1[1] - q1[0]) - step_length = ki1[1] - ki1[0] - induced_local_density = np.zeros((q1.shape[0], q2.shape[0]))*(1+0j) - print() - i0 = 0 - for q10 in q1: - print('i0=', i0) - j0 = 0 - for q20 in q2: - for ki10 in ki1: - for ki20 in ki2: - green_01 = green_function(fermi_energy, ki10, ki20, hamiltonian) - green_02 = green_function(fermi_energy, ki10+q10, ki20+q20, hamiltonian) - induced_green = np.dot(np.dot(green_01, t_matrix), green_02) - temp = induced_green[0, 0]-induced_green[0, 0].conj()+induced_green[2, 2]-induced_green[2, 2].conj() - induced_local_density[i0, j0] = induced_local_density[i0, j0]+temp*np.square(step_length) - j0 += 1 - i0 += 1 - write_matrix_k1_k2(q1, q2, np.real(induced_local_density*1j/np.square(2*pi)/(2*pi)), 'QPI') # 数据写入文件（临时写入，会被多次替代） - induced_local_density = np.real(induced_local_density*1j/np.square(2*pi)/(2*pi)) - return induced_local_density - - -def integral_of_green(fermi_energy, ki1, ki2, hamiltonian): # 在计算QPI时需要对格林函数积分 - dim = hamiltonian(0, 0).shape[0] - integral_value = np.zeros((dim, dim))*(1+0j) - step_length = ki1[1]-ki1[0] - for ki10 in ki1: - for ki20 in ki2: - green = green_function(fermi_energy, ki10, ki20, hamiltonian) - integral_value = integral_value+green*np.square(step_length) - return integral_value - - -def write_matrix_k1_k2(x1, x2, value, filename='matrix_k1_k2'): # 把矩阵数据写入文件（格式化输出） - with open(filename+'.txt', 'w') as f: - np.set_printoptions(suppress=True) # 取消输出科学记数法 - f.write('0 ') - for x10 in x1: - f.write(str(x10)+' ') - f.write('\n') - i0 = 0 - for x20 in x2: - f.write(str(x20)) - for j0 in range(x1.shape[0]): - f.write(' '+str(value[i0, j0])+' ') - f.write('\n') - i0 += 1 - - -def plot_contour(x1, x2, value, filename='contour'): # 直接画出contour图像（保存图像） - plt.contourf(x1, x2, value) #, cmap=plt.cm.hot) - plt.savefig(filename+'.eps') - # plt.show() - - -def hamiltonian(kx, ky): # 体系的哈密顿量 - t1 = -1; t2 = 1.3; t3 = -0.85; t4 = -0.85; delta_0 = 0.1; mu = 1.54 - epsilon_x = -2*t1*cos(kx)-2*t2*cos(ky)-4*t3*cos(kx)*cos(ky) - epsilon_y = -2*t1*cos(ky)-2*t2*cos(kx)-4*t3*cos(kx)*cos(ky) - epsilon_xy = -4*t4*sin(kx)*sin(ky) - delta_1 = delta_0*cos(kx)*cos(ky) - delta_2 = delta_0*cos(kx)*cos(ky) - h = np.zeros((4, 4)) - h[0, 0] = epsilon_x-mu - h[1, 1] = -epsilon_x+mu - h[2, 2] = epsilon_y-mu - h[3, 3] = -epsilon_y+mu - - h[0, 1] = delta_1 - h[1, 0] = delta_1 - h[0, 2] = epsilon_xy - h[2, 0] = epsilon_xy - h[0, 3] = 0 - h[3, 0] = 0 - - h[1, 2] = 0 - h[2, 1] = 0 - h[1, 3] = -epsilon_xy - h[3, 1] = -epsilon_xy - - h[2, 3] = delta_2 - h[3, 2] = delta_2 - return h - - -def main(): # 主程序 - start_clock = time.perf_counter() - fermi_energy = 0.07 # 费米能 - energy_broadening_width = 0.005 # 展宽 - k1 = np.arange(-pi, pi, 0.01) # 谱函数的图像精度 - k2 = np.arange(-pi, pi, 0.01) - spectrum = spectral_function(fermi_energy+energy_broadening_width*1j, k1, k2, hamiltonian) # 调用谱函数子程序 - write_matrix_k1_k2(k1, k2, spectrum, 'Spectral_function') # 把谱函数的数据写入文件 - # plot_contour(k1, k2, spectrum, 'Spectral_function') # 直接显示谱函数的图像（保存图像） - - q1 = np.arange(-pi, pi, 0.01) # QPI数的图像精度 - q2 = np.arange(-pi, pi, 0.01) - potential_i = (0.4+0j)*np.identity(hamiltonian(0, 0).shape[0]) # 杂质势 - potential_i[1, 1] = - potential_i[1, 1] # for nonmagnetic - potential_i[3, 3] = - potential_i[3, 3] - induced_local_density = qpi(fermi_energy+energy_broadening_width*1j, q1, q2, hamiltonian, potential_i) # 调用QPI子程序 - write_matrix_k1_k2(q1, q2, induced_local_density, 'QPI') # 把QPI数据写入文件（这里用的方法是计算结束后一次性把数据写入） - # plot_contour(q1, q2, induced_local_density, 'QPI') # 直接显示QPI图像（保存图像） - end_clock = time.perf_counter() - print('CPU执行时间=', end_clock - start_clock) - - -if __name__ == '__main__': - main() +""" +This code is supported by the website: https://www.guanjihuan.com +The newest version of this code is on the web page: https://www.guanjihuan.com/archives/3785 +""" + +import numpy as np +from math import * +import matplotlib.pyplot as plt +from matplotlib.colors import ListedColormap +import time + + +def green_function(fermi_energy, k1, k2, hamiltonian): # 计算格林函数 + matrix0 = hamiltonian(k1, k2) + dim = matrix0.shape[0] + green = np.linalg.inv(fermi_energy * np.identity(dim) - matrix0) + return green + + +def spectral_function(fermi_energy, k1, k2, hamiltonian): # 计算谱函数 + dim1 = k1.shape[0] + dim2 = k2.shape[0] + spectrum = np.zeros((dim1, dim2)) + i0 = 0 + for k10 in k1: + j0 = 0 + for k20 in k2: + green = green_function(fermi_energy, k10, k20, hamiltonian) + spectrum[i0, j0] = (np.imag(green[0,0])+np.imag(green[2,2]))/(-pi) + j0 += 1 + i0 += 1 + # print(spectrum) + print() + print('Spectral function显示的网格点 =', k1.shape[0], '*', k1.shape[0], '; 步长 =', k1[1] - k1[0]) + print() + return spectrum + + +def qpi(fermi_energy, q1, q2, hamiltonian, potential_i): # 计算QPI + dim = hamiltonian(0, 0).shape[0] + ki1 = np.arange(-pi, pi, 0.01) # 计算gamma_0时，k的积分密度 + ki2 = np.arange(-pi, pi, 0.01) + print('gamma_0的积分网格点 =', ki1.shape[0], '*', ki1.shape[0], '; 步长 =', ki1[1] - ki1[0]) + gamma_0 = integral_of_green(fermi_energy, ki1, ki2, hamiltonian)/np.square(2*pi) + t_matrix = np.dot(np.linalg.inv(np.identity(dim)-np.dot(potential_i, gamma_0)), potential_i) + ki1 = np.arange(-pi, pi, 0.06) # 计算induced_local_density时，k的积分密度 + ki2 = np.arange(-pi, pi, 0.06) + print('局域态密度变化的积分网格点 =', ki1.shape[0], '*', ki1.shape[0], '; 步长 =', ki1[1] - ki1[0]) + print('QPI显示的网格点 =', q1.shape[0], '*', q1.shape[0], '; 步长 =', q1[1] - q1[0]) + step_length = ki1[1] - ki1[0] + induced_local_density = np.zeros((q1.shape[0], q2.shape[0]))*(1+0j) + print() + i0 = 0 + for q10 in q1: + print('i0=', i0) + j0 = 0 + for q20 in q2: + for ki10 in ki1: + for ki20 in ki2: + green_01 = green_function(fermi_energy, ki10, ki20, hamiltonian) + green_02 = green_function(fermi_energy, ki10+q10, ki20+q20, hamiltonian) + induced_green = np.dot(np.dot(green_01, t_matrix), green_02) + temp = induced_green[0, 0]-induced_green[0, 0].conj()+induced_green[2, 2]-induced_green[2, 2].conj() + induced_local_density[i0, j0] = induced_local_density[i0, j0]+temp*np.square(step_length) + j0 += 1 + i0 += 1 + write_matrix_k1_k2(q1, q2, np.real(induced_local_density*1j/np.square(2*pi)/(2*pi)), 'QPI') # 数据写入文件（临时写入，会被多次替代） + induced_local_density = np.real(induced_local_density*1j/np.square(2*pi)/(2*pi)) + return induced_local_density + + +def integral_of_green(fermi_energy, ki1, ki2, hamiltonian): # 在计算QPI时需要对格林函数积分 + dim = hamiltonian(0, 0).shape[0] + integral_value = np.zeros((dim, dim))*(1+0j) + step_length = ki1[1]-ki1[0] + for ki10 in ki1: + for ki20 in ki2: + green = green_function(fermi_energy, ki10, ki20, hamiltonian) + integral_value = integral_value+green*np.square(step_length) + return integral_value + + +def write_matrix_k1_k2(x1, x2, value, filename='matrix_k1_k2'): # 把矩阵数据写入文件（格式化输出） + with open(filename+'.txt', 'w') as f: + np.set_printoptions(suppress=True) # 取消输出科学记数法 + f.write('0 ') + for x10 in x1: + f.write(str(x10)+' ') + f.write('\n') + i0 = 0 + for x20 in x2: + f.write(str(x20)) + for j0 in range(x1.shape[0]): + f.write(' '+str(value[i0, j0])+' ') + f.write('\n') + i0 += 1 + + +def plot_contour(x1, x2, value, filename='contour'): # 直接画出contour图像（保存图像） + plt.contourf(x1, x2, value) #, cmap=plt.cm.hot) + plt.savefig(filename+'.eps') + # plt.show() + + +def hamiltonian(kx, ky): # 体系的哈密顿量 + t1 = -1; t2 = 1.3; t3 = -0.85; t4 = -0.85; delta_0 = 0.1; mu = 1.54 + epsilon_x = -2*t1*cos(kx)-2*t2*cos(ky)-4*t3*cos(kx)*cos(ky) + epsilon_y = -2*t1*cos(ky)-2*t2*cos(kx)-4*t3*cos(kx)*cos(ky) + epsilon_xy = -4*t4*sin(kx)*sin(ky) + delta_1 = delta_0*cos(kx)*cos(ky) + delta_2 = delta_0*cos(kx)*cos(ky) + h = np.zeros((4, 4)) + h[0, 0] = epsilon_x-mu + h[1, 1] = -epsilon_x+mu + h[2, 2] = epsilon_y-mu + h[3, 3] = -epsilon_y+mu + + h[0, 1] = delta_1 + h[1, 0] = delta_1 + h[0, 2] = epsilon_xy + h[2, 0] = epsilon_xy + h[0, 3] = 0 + h[3, 0] = 0 + + h[1, 2] = 0 + h[2, 1] = 0 + h[1, 3] = -epsilon_xy + h[3, 1] = -epsilon_xy + + h[2, 3] = delta_2 + h[3, 2] = delta_2 + return h + + +def main(): # 主程序 + start_clock = time.perf_counter() + fermi_energy = 0.07 # 费米能 + energy_broadening_width = 0.005 # 展宽 + k1 = np.arange(-pi, pi, 0.01) # 谱函数的图像精度 + k2 = np.arange(-pi, pi, 0.01) + spectrum = spectral_function(fermi_energy+energy_broadening_width*1j, k1, k2, hamiltonian) # 调用谱函数子程序 + write_matrix_k1_k2(k1, k2, spectrum, 'Spectral_function') # 把谱函数的数据写入文件 + # plot_contour(k1, k2, spectrum, 'Spectral_function') # 直接显示谱函数的图像（保存图像） + + q1 = np.arange(-pi, pi, 0.01) # QPI数的图像精度 + q2 = np.arange(-pi, pi, 0.01) + potential_i = (0.4+0j)*np.identity(hamiltonian(0, 0).shape[0]) # 杂质势 + potential_i[1, 1] = - potential_i[1, 1] # for nonmagnetic + potential_i[3, 3] = - potential_i[3, 3] + induced_local_density = qpi(fermi_energy+energy_broadening_width*1j, q1, q2, hamiltonian, potential_i) # 调用QPI子程序 + write_matrix_k1_k2(q1, q2, induced_local_density, 'QPI') # 把QPI数据写入文件（这里用的方法是计算结束后一次性把数据写入） + # plot_contour(q1, q2, induced_local_density, 'QPI') # 直接显示QPI图像（保存图像） + end_clock = time.perf_counter() + print('CPU执行时间=', end_clock - start_clock) + + +if __name__ == '__main__': + main()