学术, 量子输运

局域电流的计算(附Python代码)

一、使用格林函数算法进行计算

局域电流 (Local current) 公式见参考资料[1,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/3888
"""

import numpy as np
import matplotlib.pyplot as plt
import time


def matrix_00(width=10):  # 电极元胞内跃迁,width不赋值时默认为10  
    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=10):  # 电极元胞间跃迁,width不赋值时默认为10
    h01 = np.identity(width)
    return h01


def matrix_LC(width=10, length=300):  # 左电极跳到中心区
    h_LC = np.zeros((width, width*length))
    for width0 in range(width):
        h_LC[width0, width0] = 1
    return h_LC


def matrix_CR(width=10, length=300):  # 中心区跳到右电极
    h_CR = np.zeros((width*length, width))
    for width0 in range(width):
        h_CR[width*(length-1)+width0, width0] = 1
    return h_CR


def matrix_center(width=10, length=300):  # 中心区哈密顿量
    hamiltonian = np.zeros((width*length, width*length))
    for length0 in range(length-1):
        for width0 in range(width):
            hamiltonian[width*length0+width0, width*(length0+1)+width0] = 1  # 长度方向跃迁
            hamiltonian[width*(length0+1)+width0, width*length0+width0] = 1 
    for length0 in range(length):
        for width0 in range(width-1):
            hamiltonian[width*length0+width0, width*length0+width0+1] = 1  # 宽度方向跃迁
            hamiltonian[width*length0+width0+1, width*length0+width0] = 1
    # 中间加势垒
    for j0 in range(6):
        for i0 in range(6): 
            hamiltonian[width*(int(length/2)-3+j0)+int(width/2)-3+i0, width*(int(length/2)-3+j0)+int(width/2)-3+i0]= 1e8
    return hamiltonian


def main():
    start_time = time.time()
    fermi_energy = 0
    width = 60
    length = 100
    h00 = matrix_00(width)
    h01 = matrix_01(width)
    G_n = Green_n(fermi_energy+1j*1e-9, h00, h01, width, length)
    # 下面是提取数据并画图
    direction_x = np.zeros((width, length))
    direction_y = np.zeros((width, length))
    for length0 in range(length-1):
        for width0 in range(width):
            direction_x[width0, length0] = G_n[width*length0+width0, width*(length0+1)+width0]
    for length0 in range(length):
        for width0 in range(width-1):
            direction_y[width0, length0] = G_n[width*length0+width0, width*length0+width0+1]
    # print(direction_x)
    # print(direction_y)
    X, Y = np.meshgrid(range(length), range(width))
    plt.quiver(X, Y, direction_x, direction_y)
    plt.show()
    end_time = time.time()
    print('运行时间=', end_time-start_time)



def transfer_matrix(fermi_energy, h00, h01, dim):   # 转移矩阵T。dim是传递矩阵h00和h01的维度
    transfer = np.zeros((2*dim, 2*dim), dtype=complex)
    transfer[0:dim, 0:dim] = np.dot(np.linalg.inv(h01), fermi_energy*np.identity(dim)-h00)   # np.dot()等效于np.matmul()
    transfer[0:dim, dim:2*dim] = np.dot(-1*np.linalg.inv(h01), h01.transpose().conj())
    transfer[dim:2*dim, 0:dim] = np.identity(dim)
    transfer[dim:2*dim, dim:2*dim] = 0  # a:b代表 a <= x < b
    return transfer  # 返回转移矩阵


def green_function_lead(fermi_energy, h00, h01, dim):  # 电极的表面格林函数
    transfer = transfer_matrix(fermi_energy, h00, h01, dim)
    eigenvalue, eigenvector = np.linalg.eig(transfer)
    ind = np.argsort(np.abs(eigenvalue))
    temp = np.zeros((2*dim, 2*dim))*(1+0j)
    i0 = 0
    for ind0 in ind:
        temp[:, i0] = eigenvector[:, ind0]
        i0 += 1
    s1 = temp[dim:2*dim, 0:dim]
    s2 = temp[0:dim, 0:dim]
    s3 = temp[dim:2*dim, dim:2*dim]
    s4 = temp[0:dim, dim:2*dim]
    right_lead_surface = np.linalg.inv(fermi_energy*np.identity(dim)-h00-np.dot(np.dot(h01, s2), np.linalg.inv(s1)))
    left_lead_surface = np.linalg.inv(fermi_energy*np.identity(dim)-h00-np.dot(np.dot(h01.transpose().conj(), s3), np.linalg.inv(s4)))
    return right_lead_surface, left_lead_surface  # 返回右电极的表面格林函数和左电极的表面格林函数


def self_energy_lead(fermi_energy, h00, h01, width, length):  # 电极的自能
    h_LC = matrix_LC(width, length)
    h_CR = matrix_CR(width, length)
    right_lead_surface, left_lead_surface = green_function_lead(fermi_energy, h00, h01, width)
    right_self_energy = np.dot(np.dot(h_CR, right_lead_surface), h_CR.transpose().conj())
    left_self_energy = np.dot(np.dot(h_LC.transpose().conj(), left_lead_surface), h_LC)
    return right_self_energy, left_self_energy  # 返回右电极的自能和左电极的自能


def Green_n(fermi_energy, h00, h01, width, length):  # 计算G_n
    right_self_energy, left_self_energy = self_energy_lead(fermi_energy, h00, h01, width, length)
    hamiltonian = matrix_center(width, length)
    green = np.linalg.inv(fermi_energy*np.identity(width*length)-hamiltonian-left_self_energy-right_self_energy)
    right_self_energy = (right_self_energy - right_self_energy.transpose().conj())*1j
    left_self_energy = (left_self_energy - left_self_energy.transpose().conj())*1j
    G_n = np.imag(np.dot(np.dot(green, left_self_energy), green.transpose().conj()))
    return G_n


if __name__ == '__main__':
    main()

当费米能fermi_energy=0时,计算结果:

二、使用Guan软件包计算

Guan官网:https://py.guanjihuan.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/3888
"""

import numpy as np
import matplotlib.pyplot as plt
import time
import guan

def matrix_00(width=10):  # 电极元胞内跃迁,width不赋值时默认为10  
    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=10):  # 电极元胞间跃迁,width不赋值时默认为10
    h01 = np.identity(width)
    return h01


def matrix_LC(width=10, length=300):  # 左电极跳到中心区
    h_LC = np.zeros((width, width*length))
    for width0 in range(width):
        h_LC[width0, width0] = 1
    return h_LC


def matrix_CR(width=10, length=300):  # 中心区跳到右电极
    h_CR = np.zeros((width*length, width))
    for width0 in range(width):
        h_CR[width*(length-1)+width0, width0] = 1
    return h_CR


def matrix_center(width=10, length=300):  # 中心区哈密顿量
    hamiltonian = np.zeros((width*length, width*length))
    for length0 in range(length-1):
        for width0 in range(width):
            hamiltonian[width*length0+width0, width*(length0+1)+width0] = 1  # 长度方向跃迁
            hamiltonian[width*(length0+1)+width0, width*length0+width0] = 1 
    for length0 in range(length):
        for width0 in range(width-1):
            hamiltonian[width*length0+width0, width*length0+width0+1] = 1  # 宽度方向跃迁
            hamiltonian[width*length0+width0+1, width*length0+width0] = 1
    # 中间加势垒
    for j0 in range(6):
        for i0 in range(6): 
            hamiltonian[width*(int(length/2)-3+j0)+int(width/2)-3+i0, width*(int(length/2)-3+j0)+int(width/2)-3+i0]= 1e8
    return hamiltonian


def main():
    start_time = time.time()
    fermi_energy = 0
    width = 60
    length = 100
    h00 = matrix_00(width)
    h01 = matrix_01(width)
    h_LC = matrix_LC(width, length)
    h_CR = matrix_CR(width, length)
    hamiltonian = matrix_center(width, length) 
    
    G_n = guan.electron_correlation_function_green_n_for_local_current(fermi_energy, h00, h01, h_LC, h_CR, hamiltonian)

    direction_x = np.zeros((width, length))
    direction_y = np.zeros((width, length))
    for length0 in range(length-1):
        for width0 in range(width):
            direction_x[width0, length0] = G_n[width*length0+width0, width*(length0+1)+width0]
    for length0 in range(length):
        for width0 in range(width-1):
            direction_y[width0, length0] = G_n[width*length0+width0, width*length0+width0+1]

    X, Y = np.meshgrid(range(length), range(width))
    plt.quiver(X, Y, direction_x, direction_y)
    plt.show()
    end_time = time.time()
    print('运行时间=', end_time-start_time)


if __name__ == '__main__':
    main()

三、使用Kwant软件包计算

Kwant官网:https://kwant-project.org/

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

import kwant


def make_system():
    a = 1
    lat = kwant.lattice.square(a, norbs=1)  # 创建晶格,方格子
    syst = kwant.Builder() 
    t = 1.0
    W = 60  # 中心体系宽度
    L = 100  # 中心体系长度
    # 给中心体系赋值
    for i in range(L):
        for j in range(W):
            syst[lat(i, j)] = 0
            if j > 0:
                syst[lat(i, j), lat(i, j-1)] = -t  # hopping in y-direction
            if i > 0:
                syst[lat(i, j), lat(i-1, j)] = -t  # hopping in x-direction
            if 47<=i<53 and 27<=j<33:  # 势垒
                syst[lat(i, j)] = 1e8
    # 电极
    lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
    lead[(lat(0, j) for j in range(W))] = 0
    lead[lat.neighbors()] = -t  # 用neighbors()方法
    syst.attach_lead(lead)  # 左电极
    syst.attach_lead(lead.reversed())   # 用reversed()方法得到右电极
    # 制作结束
    kwant.plot(syst) 
    syst = syst.finalized()  
    return syst


def main():
    syst = make_system()
    psi = kwant.wave_function(syst)(0)[29]
    current = kwant.operator.Current(syst)(psi)
    kwant.plotter.current(syst, current)


if __name__ == '__main__':
    main()

选择中间能级的波函数:psi = kwant.wave_function(syst)(0)[29],计算结果:

参考资料:

[1] Quantum blockade and loop currents in graphene with topological defects

[2] Numerical study of the topological Anderson insulator in HgTe/CdTe quantum wells

2,121 次浏览

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

评论说明:
(1)在保留浏览器缓存的前提下,目前支持72小时自主修改或删除个人评论。如果自己无法修改或删除评论,可再次评论或联系我。如有发现广告留言,请勿点击链接,博主会不定期删除。
(2)评论支持Latex公式。把latexpage作为标签放在任何位置,评论中的公式可正常编译,示例:
$Latex formula$  [latexpage]

4 thoughts on “局域电流的计算(附Python代码)”

  1. 您好,请问利用格林函数求解局域电流有参考文献嘛 或者编程前的公式?

    1. 公式我是根据本篇博文底部给出的参考资料。可以在此基础上继续搜相关的文献或公式的原始文献。

发表回复

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