本篇给出另外一种计算格林函数的对角分块矩阵的方法（由zhangjiayan同学提供），公式会比之前这篇更简洁些：使用Dyson方程迭代方法计算态密度（附Python代码）。思想是把每一层左右两侧都当成自能，公式参考书籍“2016 - Ryndyk - Theory of Quantum Transport at Nanoscale”的（3.152）式：
![G_{ii}=[E-H_{ii}-H_{i,i-1}G_{i-1, i-1}^{\rightarrow (i-1)} H_{i-1, i}-H_{i, i+1}G_{i+1, i+1}^{\leftarrow(i+1)}H_{i+1, i}]^{-1}](https://www.guanjihuan.com/wp-content/ql-cache/quicklatex.com-a5748a2c8af20dd51547939f7a5f2720_l3.png "Rendered by QuickLaTeX.com")
说明：和之前那篇的代码一样，这里不考虑电极的自能。同时，由于这里考虑的是简单方格子，因此是左右对称的，如果有更复杂的情况，您需要对代码进行改进后再使用。
代码为：
"""
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/38466
"""
import numpy as np
def matrix_00(width):
h00 = np.zeros((width, width))
for width0 in range(width-1):
h00[width0, width0+1] = 1
h00[width0+1, width0] = 1
return h00
def matrix_01(width):
h01 = np.identity(width)
return h01
def main():
width = 2
length = 3
eta = 1e-2
E = 0
h00 = matrix_00(width)
h01 = matrix_01(width)
G_ii_n_array = G_ii_n_with_Dyson_equation_2(width, length, E, eta, h00, h01)
for i0 in range(length):
# print('G_{'+str(i0+1)+','+str(i0+1)+'}^{('+str(length)+')}:')
# print(G_ii_n_array[i0, :, :],'\n')
print('x=', i0+1, ':')
for j0 in range(width):
print(' y=', j0+1, ' ', -np.imag(G_ii_n_array[i0, j0, j0])/np.pi) # 态密度
def G_ii_n_with_Dyson_equation_2(width, length, E, eta, h00, h01):
G_ii_n_array = np.zeros((length, width, width), complex)
G_11_1 = np.linalg.inv((E+eta*1j)*np.identity(width)-h00)
for i in range(length):
G_nn_n_right_minus = G_11_1
G_nn_n_left_minus = G_11_1
if i!=0:
for _ in range(i-1):
G_nn_n_right = Green_nn_n(E, eta, h00, h01, G_nn_n_right_minus)
G_nn_n_right_minus = G_nn_n_right
if i!=length-1:
for _ in range(length-i-2):
G_nn_n_left = Green_nn_n(E, eta, h00, h01, G_nn_n_left_minus)
G_nn_n_left_minus = G_nn_n_left
if i==0:
G_ii_n_array[i, :, :] = np.linalg.inv((E+eta*1j)*np.identity(width)-h00-np.dot(np.dot(h01, G_nn_n_left_minus), h01.transpose().conj()))
elif i!=0 and i!=length-1:
G_ii_n_array[i, :, :] = np.linalg.inv((E+eta*1j)*np.identity(width)-h00-np.dot(np.dot(h01.transpose().conj(), G_nn_n_right_minus), h01)-np.dot(np.dot(h01, G_nn_n_left_minus), h01.transpose().conj()))
elif i==length-1:
G_ii_n_array[i, :, :] = np.linalg.inv((E+eta*1j)*np.identity(width)-h00-np.dot(np.dot(h01.transpose().conj(), G_nn_n_right_minus), h01))
return G_ii_n_array
def Green_nn_n(E, eta, H00, V, G_nn_n_minus):
dim = H00.shape[0]
G_nn_n = np.linalg.inv((E+eta*1j)*np.identity(dim)-H00-np.dot(np.dot(V.transpose().conj(), G_nn_n_minus), V))
return G_nn_n
if __name__ == '__main__':
main()运行结果：
x= 1 :
y= 1 0.00636333452533021
y= 2 0.00636333452533021
x= 2 :
y= 1 0.00954388846688089
y= 2 0.009543888466880892
x= 3 :
y= 1 0.00636333452533021
y= 2 0.00636333452533021和之前那篇的计算结果在误差范围内是相等的。
为了进一步检查代码的正确性，不存在语义Bug，这里设置另外一组参数：width = 3，length = 6，作为简单的测试，运行结果为：
x= 1 :
y= 1 0.01587742510743398
y= 2 0.02221000774071636
y= 3 0.01587742510743398
x= 2 :
y= 1 0.019041813622148983
y= 2 0.0349014828667101
y= 3 0.019041813622148986
x= 3 :
y= 1 0.009535013625636614
y= 2 0.01270669288490087
y= 3 0.009535013625636612
x= 4 :
y= 1 0.009535013625636614
y= 2 0.012706692884900868
y= 3 0.009535013625636614
x= 5 :
y= 1 0.019041813622148983
y= 2 0.0349014828667101
y= 3 0.019041813622148986
x= 6 :
y= 1 0.01587742510743398
y= 2 0.02221000774071636
y= 3 0.01587742510743398和直接对哈密顿量求逆的算法计算结果进行比较：
x= 1 :
y= 1 0.015877425107433986
y= 2 0.022210007740714695
y= 3 0.015877425107433993
x= 2 :
y= 1 0.01904181362214891
y= 2 0.03490148286670733
y= 3 0.019041813622148913
x= 3 :
y= 1 0.009535013625636617
y= 2 0.012706692884900781
y= 3 0.009535013625636614
x= 4 :
y= 1 0.009535013625636614
y= 2 0.01270669288490078
y= 3 0.009535013625636616
x= 5 :
y= 1 0.019041813622148913
y= 2 0.034901482866707335
y= 3 0.019041813622148913
x= 6 :
y= 1 0.015877425107433993
y= 2 0.022210007740714695
y= 3 0.015877425107433993计算结果在误差范围内是相等的。
【说明：本站主要是个人的一些笔记和代码分享，内容可能会不定期修改。为了使全网显示的始终是最新版本，这里的文章未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】
