学术, 电子态密度

使用Dyson方程迭代方法计算态密度（附Python代码）

发布时间：

最近修改时间：2024年1月25日

补充说明：除了本篇的计算方法，还有另外一种比本篇更加简单的算法（由zhangjiayan同学提供），思想是把每一层左右两侧都当成自能。该算法的代码可以参考这篇：使用Dyson方程计算格林函数的对角分块矩阵（第二种方法）。公式参考书籍“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}

此外也可以阅读这几篇：

在计算电导时，我们只需要知道格林函数右上角的一个分块矩阵G_{1n}。在计算态密度公式中，我们需要知道格林函数对角线上的元素，因此需要计算出格林函数所有的对角分块矩阵G_{11}，G_{22}，G_{33}……

本篇通过不断套用Dyson方程，计算格林函数的对角分块矩阵，进而计算出态密度。考虑的体系是有限长宽的方格子。用上Dyson方程后计算的矩阵维度可以大幅度降低，矩阵维度从width*length下降到width。

戴森方程（Dyson equation） 为[1]：

计算G_{1n}具体公式可以写为：

\begin{aligned}G_{11}^{(1)}&=[(E+i\eta)I-H_{11}]^{-1}\\G_{12}^{(2)}&=G_{11}^{(1)}V_{12}G_{22}^{(2)}\\G_{13}^{(3)}&=G_{12}^{(2)}V_{23}G_{33}^{(3)}\\&\quad\quad......\end{aligned}

其中，

\begin{aligned}G_{22}^{(2)}&=[(E+i\eta)I-H_{22}-V_{12}^{\dagger}G_{11}^{(1)}V_{12}]^{-1}\\G_{33}^{(3)}&=[(E+i\eta)I-H_{33}-V_{23}^{\dagger}G_{22}^{(2)}V_{23}]^{-1}\\&\quad\quad\quad\quad\quad\quad......\end{aligned}

通常，H_{11}=H_{22}=H_{33}，V_{12}=V_{23}。

此外，如果考虑了左电极的自能\Sigma_{L}，则G_{11}^{(1)}=[(E+i\eta)I-H_{11}-\Sigma_{L}]^{-1}

需要说明的是：格林函数G^{(x)}的右上角括号(x)表示当前的体系大小。因此，这里的G_{11}^{(1)}，G_{22}^{(2)}，G_{33}^{(3)}……不是计算态密度所需要的对角分块矩阵。我们需要的是G_{11}^{(n)}，G_{22}^{(n)}，G_{33}^{(n)}……其中，(n)表示需要计算的体系大小。

下面给出具体的表达式：

\begin{aligned}G_{11}^{(1)}&=[(E+i\eta)I-H_{11}]^{-1}\\G_{11}^{(2)}&=G_{11}^{(1)}+G_{11}^{(1)}V_{12}G_{22}^{(2)}V_{12}^{\dagger}G_{11}^{(1)}\\G_{11}^{(3)}&=G_{11}^{(2)}+G_{12}^{(2)}V_{23}G_{33}^{(3)}V_{23}^{\dagger}G_{21}^{(2)}\\G_{11}^{(4)}&=G_{11}^{(3)}+G_{13}^{(3)}V_{34}G_{44}^{(4)}V_{34}^{\dagger}G_{31}^{(3)}\\&\quad\quad\quad\quad......\end{aligned}

\begin{aligned}G_{22}^{(2)}&=[(E+i\eta)I-H_{22}-V_{12}^{\dagger}G_{11}^{(1)}V_{12}]^{-1}\\G_{22}^{(3)}&=G_{22}^{(2)}+G_{22}^{(2)}V_{23}G_{33}^{(3)}V_{23}^{\dagger}G_{22}^{(2)}\\G_{22}^{(4)}&=G_{22}^{(3)}+G_{23}^{(3)}V_{34}G_{44}^{(4)}V_{34}^{\dagger}G_{32}^{(3)}\\G_{22}^{(5)}&=G_{22}^{(4)}+G_{24}^{(4)}V_{45}G_{55}^{(5)}V_{45}^{\dagger}G_{42}^{(3)}\\&\quad\quad\quad\quad\quad......\end{aligned}

其中，

\begin{aligned}G_{21}^{(2)}&=G_{22}^{(2)}V_{12}^{\dagger}G_{11}^{(1)}\\G_{31}^{(3)}&=G_{33}^{(3)}V_{23}^{\dagger}G_{21}^{(2)}\\&\quad......\end{aligned}

\begin{aligned}G_{23}^{(3)}&=G_{22}^{(2)}V_{23}G_{33}^{(3)}\\G_{24}^{(4)}&=G_{23}^{(3)}V_{34}G_{44}^{(4)}\\&\quad......\end{aligned}

\begin{aligned}G_{32}^{(3)}&=G_{33}^{(3)}V_{23}^{\dagger}G_{22}^{(2)}\\G_{42}^{(4)}&=G_{44}^{(4)}V_{34}^{\dagger}G_{32}^{(3)}\\&\quad......\end{aligned}

同理有G_{33}^{(n)}，G_{44}^{(n)}，G_{55}^{(n)}等。上面的表达式看起来虽然复杂，但都只是在套用Dyson方程的表达式。

1. 方格子

Python代码如下：

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

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(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(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):  # i为格林函数的右下指标
        # 初始化开始
        G_nn_n_minus = G_11_1
        G_in_n_minus = G_11_1
        G_ni_n_minus = G_11_1
        G_ii_n_minus = G_11_1
        for i0 in range(i):
            G_nn_n = Green_nn_n(E, eta, h00, h01, G_nn_n_minus)
            G_nn_n_minus = G_nn_n
        if i!=0:
            G_in_n_minus = G_nn_n
            G_ni_n_minus = G_nn_n
            G_ii_n_minus = G_nn_n
        # 初始化结束
        for j0 in range(length-1-i): # j0为格林函数的右上指标，表示当前体系大小，即G^{(j0)}
            G_nn_n = Green_nn_n(E, eta, h00, h01, G_nn_n_minus)
            G_nn_n_minus = G_nn_n

            G_ii_n = Green_ii_n(G_ii_n_minus, G_in_n_minus, h01, G_nn_n, G_ni_n_minus)  # 需要求的对角分块矩阵
            G_ii_n_minus = G_ii_n

            G_in_n = Green_in_n(G_in_n_minus, h01, G_nn_n)
            G_in_n_minus = G_in_n

            G_ni_n = Green_ni_n(G_nn_n, h01, G_ni_n_minus)
            G_ni_n_minus = G_ni_n
        G_ii_n_array[i, :, :] = G_ii_n_minus
    return G_ii_n_array

def Green_nn_n(E, eta, H00, V, G_nn_n_minus): # n>=2
    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

def Green_in_n(G_in_n_minus, V, G_nn_n):  # n>=2
    G_in_n = np.dot(np.dot(G_in_n_minus, V), G_nn_n)
    return G_in_n

def Green_ni_n(G_nn_n, V, G_ni_n_minus): # n>=2
    G_ni_n = np.dot(np.dot(G_nn_n, V.transpose().conj()), G_ni_n_minus)
    return G_ni_n

def Green_ii_n(G_ii_n_minus, G_in_n_minus, V, G_nn_n, G_ni_n_minus):  # n>=i
    G_ii_n = G_ii_n_minus+np.dot(np.dot(np.dot(np.dot(G_in_n_minus, V), G_nn_n), V.transpose().conj()),G_ni_n_minus)
    return G_ii_n

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

import numpy as np

def hamiltonian(width, length):  
    h = np.zeros((width*length, width*length))
    for i0 in range(length):
        for j0 in range(width-1):
            h[i0*width+j0, i0*width+j0+1] = 1
            h[i0*width+j0+1, i0*width+j0] = 1
    for i0 in range(length-1):
        for j0 in range(width):
            h[i0*width+j0, (i0+1)*width+j0] = 1
            h[(i0+1)*width+j0, i0*width+j0] = 1
    return h

def main():
    width = 2
    length = 3
    h = hamiltonian(width, length)
    E = 0
    green = np.linalg.inv((E+1e-2j)*np.eye(width*length)-h)  
    for i0 in range(length):
        # print('G_{'+str(i0+1)+','+str(i0+1)+'}^{('+str(length)+')}:')
        # print(green[i0*width+0: i0*width+width, i0*width+0: i0*width+width], '\n') 
        print('x=', i0+1, ':')
        for j0 in range(width):
            print('     y=', j0+1, ' ', -np.imag(green[i0*width+j0, i0*width+j0])/np.pi) 

if __name__ == "__main__":
    main()

计算结果为：

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

2. 立方格子

把每一层切片当成一个单元，进行Dyson方程的迭代。代码为：

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

import numpy as np

def matrix_00(width, length):  
    h00 = np.zeros((width*length, width*length))
    for x in range(length):
        for y in range(width-1):
            h00[x*width+y, x*width+y+1] = 1
            h00[x*width+y+1, x*width+y] = 1
    for x in range(length-1):
        for y in range(width):
            h00[x*width+y, (x+1)*width+y] = 1
            h00[(x+1)*width+y, x*width+y] = 1
    return h00

def matrix_01(width, length): 
    h01 = np.identity(width*length)
    return h01

def main():
    height = 2  # z
    width = 3   # y
    length = 4  # x
    eta = 1e-2
    E = 0
    h00 = matrix_00(width, length)
    h01 = matrix_01(width, length)
    G_ii_n_array = G_ii_n_with_Dyson_equation(width, length, height, E, eta, h00, h01)
    for i0 in range(height):
        print('z=', i0+1, ':')
        for j0 in range(width):
            print('      y=', j0+1, ':')
            for k0 in range(length):
                print('             x=', k0+1, ' ', -np.imag(G_ii_n_array[i0, k0*width+j0, k0*width+j0])/np.pi)   # 态密度

def G_ii_n_with_Dyson_equation(width, length, height, E, eta, h00, h01):
    dim = length*width
    G_ii_n_array = np.zeros((height, dim, dim), dtype=complex)
    G_11_1 = np.linalg.inv((E+eta*1j)*np.identity(dim)-h00)
    for i in range(height):  # i为格林函数的右下指标
        # 初始化开始
        G_nn_n_minus = G_11_1
        G_in_n_minus = G_11_1
        G_ni_n_minus = G_11_1
        G_ii_n_minus = G_11_1
        for i0 in range(i):
            G_nn_n = Green_nn_n(E, eta, h00, h01, G_nn_n_minus)
            G_nn_n_minus = G_nn_n
        if i!=0:
            G_in_n_minus = G_nn_n
            G_ni_n_minus = G_nn_n
            G_ii_n_minus = G_nn_n
        # 初始化结束
        for j0 in range(height-1-i): # j0为格林函数的右上指标，表示当前体系大小，即G^{(j0)}
            G_nn_n = Green_nn_n(E, eta, h00, h01, G_nn_n_minus)
            G_nn_n_minus = G_nn_n

            G_ii_n = Green_ii_n(G_ii_n_minus, G_in_n_minus, h01, G_nn_n, G_ni_n_minus)  # 需要求的对角分块矩阵
            G_ii_n_minus = G_ii_n

            G_in_n = Green_in_n(G_in_n_minus, h01, G_nn_n)
            G_in_n_minus = G_in_n

            G_ni_n = Green_ni_n(G_nn_n, h01, G_ni_n_minus)
            G_ni_n_minus = G_ni_n
        G_ii_n_array[i, :, :] = G_ii_n_minus
    return G_ii_n_array

def Green_nn_n(E, eta, H00, V, G_nn_n_minus): # n>=2
    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

def Green_in_n(G_in_n_minus, V, G_nn_n):  # n>=2
    G_in_n = np.dot(np.dot(G_in_n_minus, V), G_nn_n)
    return G_in_n

def Green_ni_n(G_nn_n, V, G_ni_n_minus): # n>=2
    G_ni_n = np.dot(np.dot(G_nn_n, V.transpose().conj()), G_ni_n_minus)
    return G_ni_n

def Green_ii_n(G_ii_n_minus, G_in_n_minus, V, G_nn_n, G_ni_n_minus):  # n>=i
    G_ii_n = G_ii_n_minus+np.dot(np.dot(np.dot(np.dot(G_in_n_minus, V), G_nn_n), V.transpose().conj()),G_ni_n_minus)
    return G_ii_n

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

import numpy as np

def hamiltonian(width, length, height):  
    h = np.zeros((width*length*height, width*length*height))
    for i0 in range(length):
        for j0 in range(width):
            for k0 in range(height-1):
                h[k0*width*length+i0*width+j0, (k0+1)*width*length+i0*width+j0] = 1
                h[(k0+1)*width*length+i0*width+j0, k0*width*length+i0*width+j0] = 1
    for i0 in range(length):
        for j0 in range(width-1):
            for k0 in range(height):
                h[k0*width*length+i0*width+j0, k0*width*length+i0*width+j0+1] = 1
                h[k0*width*length+i0*width+j0+1, k0*width*length+i0*width+j0] = 1
    for i0 in range(length-1):
        for j0 in range(width):
            for k0 in range(height):
                h[k0*width*length+i0*width+j0, k0*width*length+(i0+1)*width+j0] = 1
                h[k0*width*length+(i0+1)*width+j0, k0*width*length+i0*width+j0] = 1
    return h

def main():
    height = 2  # z
    width = 3  # y
    length = 4  # x
    h = hamiltonian(width, length, height)
    E = 0
    green = np.linalg.inv((E+1e-2j)*np.eye(width*length*height)-h)  
    for k0 in range(height):
        print('z=', k0+1, ':')
        for j0 in range(width):
            print('      y=', j0+1, ':')
            for i0 in range(length):
                print('             x=', i0+1, ' ', -np.imag(green[k0*width*length+i0*width+j0, k0*width*length+i0*width+j0])/np.pi)   # 态密度

if __name__ == "__main__":
    main()

计算结果为：

结果在误差范围内是一致的。

另外，下面把利用Dyson方程迭代方法计算态密度（立方格子）的代码做了一些细节上的修改，这里称为第二版本“version II"，主要是删除G_ii_n_array这个变量，因为这个变量太占内存了。经过处理后，在有限的内存下，能够计算的体系大小会更大一些，消耗的内存降低为1/height倍的量级。

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

import numpy as np

def matrix_00(width, length):  
    h00 = np.zeros((width*length, width*length))
    for x in range(length):
        for y in range(width-1):
            h00[x*width+y, x*width+y+1] = 1
            h00[x*width+y+1, x*width+y] = 1
    for x in range(length-1):
        for y in range(width):
            h00[x*width+y, (x+1)*width+y] = 1
            h00[(x+1)*width+y, x*width+y] = 1
    return h00

def matrix_01(width, length): 
    h01 = np.identity(width*length)
    return h01
    
def main():
    height = 2  # z
    width = 3   # y
    length = 4  # x
    eta = 1e-2
    E = 0
    h00 = matrix_00(width, length)
    h01 = matrix_01(width, length)
    G_ii_n_with_Dyson_equation_version_II(width, length, height, E, eta, h00, h01)

def G_ii_n_with_Dyson_equation_version_II(width, length, height, E, eta, h00, h01):
    dim = length*width
    G_11_1 = np.linalg.inv((E+eta*1j)*np.identity(dim)-h00)
    for i in range(height):  # i为格林函数的右下指标
        # 初始化开始
        G_nn_n_minus = G_11_1
        G_in_n_minus = G_11_1
        G_ni_n_minus = G_11_1
        G_ii_n_minus = G_11_1
        for i0 in range(i):
            G_nn_n = Green_nn_n(E, eta, h00, h01, G_nn_n_minus)
            G_nn_n_minus = G_nn_n
        if i!=0:
            G_in_n_minus = G_nn_n
            G_ni_n_minus = G_nn_n
            G_ii_n_minus = G_nn_n
        # 初始化结束
        for j0 in range(height-1-i): # j0为格林函数的右上指标，表示当前体系大小，即G^{(j0)}
            G_nn_n = Green_nn_n(E, eta, h00, h01, G_nn_n_minus)
            G_nn_n_minus = G_nn_n

            G_ii_n = Green_ii_n(G_ii_n_minus, G_in_n_minus, h01, G_nn_n, G_ni_n_minus)  # 需要求的对角分块矩阵
            G_ii_n_minus = G_ii_n

            G_in_n = Green_in_n(G_in_n_minus, h01, G_nn_n)
            G_in_n_minus = G_in_n

            G_ni_n = Green_ni_n(G_nn_n, h01, G_ni_n_minus)
            G_ni_n_minus = G_ni_n
        # 输出
        print('z=', i+1, ':')
        for j0 in range(width):
            print('      y=', j0+1, ':')
            for k0 in range(length):
                print('             x=', k0+1, ' ', -np.imag(G_ii_n_minus[k0*width+j0, k0*width+j0])/np.pi)   # 态密度

def Green_nn_n(E, eta, H00, V, G_nn_n_minus): # n>=2
    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

def Green_in_n(G_in_n_minus, V, G_nn_n):  # n>=2
    G_in_n = np.dot(np.dot(G_in_n_minus, V), G_nn_n)
    return G_in_n

def Green_ni_n(G_nn_n, V, G_ni_n_minus): # n>=2
    G_ni_n = np.dot(np.dot(G_nn_n, V.transpose().conj()), G_ni_n_minus)
    return G_ni_n

def Green_ii_n(G_ii_n_minus, G_in_n_minus, V, G_nn_n, G_ni_n_minus):  # n>=i
    G_ii_n = G_ii_n_minus+np.dot(np.dot(np.dot(np.dot(G_in_n_minus, V), G_nn_n), V.transpose().conj()),G_ni_n_minus)
    return G_ii_n

if __name__ == '__main__': 
    main()

参考资料：

[1] MacKinnon, "A. The calculation of transport properties and density of states of disordered solids",  Z. Physik B - Condensed Matter 59, 385–390 (1985).

3,042 次浏览

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

头像

Published by guanjihuan

发表评论

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

Captcha Code