学术, 电子态密度

使用Dyson方程计算格林函数的对角分块矩阵(第二种方法)

本篇给出另外一种计算格林函数的对角分块矩阵的方法(由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}

说明:和之前那篇的代码一样,这里不考虑电极的自能。同时,由于这里考虑的是简单方格子,因此是左右对称的,如果有更复杂的情况,您需要对代码进行改进后再使用。

代码为:

"""
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

计算结果在误差范围内是相等的。

370 次浏览

【说明:本站主要是个人的一些笔记和代码分享,内容可能会不定期修改。为了使全网显示的始终是最新版本,这里的文章未经同意请勿转载。引用请注明出处:https://www.guanjihuan.com

发表评论

您的电子邮箱地址不会被公开。 必填项已用 * 标注

Captcha Code