拓扑不变量, 学术

空间反演对称性破缺的石墨烯的贝里曲率分布(附Python代码)

参考文献为:Di Xiao, Wang Yao, and Qian Niu, “Valley-Contrasting Physics in Graphene: Magnetic Moment and Topological Transport”, Phys. Rev. Lett. 99, 236809 (2007)。在文献中考虑的石墨烯模型是在狄拉克点附近的线性模型,而本篇考虑的石墨烯的紧束缚模型。

本篇采用的是贝里曲率的定义公式,其他的贝里曲率计算方法阅读这篇:贝里曲率的计算(附Python代码)。本篇的参考博文为:

石墨烯在倒空间的示意图:

1. 考虑k_y=0对称轴上的贝里曲率的分布

代码为:

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

import numpy as np
import matplotlib.pyplot as plt
from math import *  
import cmath
import time


def hamiltonian(k1, k2, t1=2.82, a=1/sqrt(3)):  # 石墨烯哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
    h = np.zeros((2, 2), dtype=complex)
    h[0, 0] = 0.28/2
    h[1, 1] = -0.28/2
    h[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a))
    h[0, 1] = h[1, 0].conj()
    return h

def main():
    start_time = time.time()
    n = 2000  # 取点密度
    delta = 1e-9  # 求导的偏离量
    for band in range(2):
        F_all = []  # 贝里曲率
        for kx in np.linspace(-2*pi, 2*pi, n):
            for ky in [0]: # 这里只考虑ky=0对称轴上的情况 # np.linspace(-pi, pi, n):
                H = hamiltonian(kx, ky)
                eigenvalue, eigenvector = np.linalg.eig(H)
                vector = eigenvector[:, np.argsort(np.real(eigenvalue))[band]]  # 价带波函数
                # print(np.argsort(np.real(eigenvalue))[0])  # 排序索引(从小到大)
                # print(eigenvalue)  # 排序前的本征值
                # print(np.sort(np.real(eigenvalue)))  # 排序后的本征值(从小到大)
            
                H_delta_kx = hamiltonian(kx+delta, ky) 
                eigenvalue, eigenvector = np.linalg.eig(H_delta_kx)
                vector_delta_kx = eigenvector[:, np.argsort(np.real(eigenvalue))[band]]   # 略偏离kx的波函数
                # vector_delta_kx = find_vector_with_the_same_gauge(vector_delta_kx, vector)  # 如果波函数不连续需要使用这个

                H_delta_ky = hamiltonian(kx, ky+delta)  
                eigenvalue, eigenvector = np.linalg.eig(H_delta_ky)
                vector_delta_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[band]]  # 略偏离ky的波函数
                # vector_delta_ky = find_vector_with_the_same_gauge(vector_delta_ky, vector)  # 如果波函数不连续需要使用这个

                H_delta_kx_ky = hamiltonian(kx+delta, ky+delta)  
                eigenvalue, eigenvector = np.linalg.eig(H_delta_kx_ky)
                vector_delta_kx_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[band]]  # 略偏离kx和ky的波函数
                # vector_delta_kx_ky = find_vector_with_the_same_gauge(vector_delta_kx_ky, vector)  # 如果波函数不连续需要使用这个

                # 价带的波函数的贝里联络(berry connection) # 求导后内积
                A_x = np.dot(vector.transpose().conj(), (vector_delta_kx-vector)/delta)   # 贝里联络Ax(x分量)
                A_y = np.dot(vector.transpose().conj(), (vector_delta_ky-vector)/delta)   # 贝里联络Ay(y分量)
                
                A_x_delta_ky = np.dot(vector_delta_ky.transpose().conj(), (vector_delta_kx_ky-vector_delta_ky)/delta)  # 略偏离ky的贝里联络Ax
                A_y_delta_kx = np.dot(vector_delta_kx.transpose().conj(), (vector_delta_kx_ky-vector_delta_kx)/delta)  # 略偏离kx的贝里联络Ay

                # 贝里曲率(berry curvature)
                F = ((A_y_delta_kx-A_y)/delta-(A_x_delta_ky-A_x)/delta)*1j
                # print(F)
                F_all = np.append(F_all,[F], axis=0) 
        plt.plot(np.linspace(-2*pi, 2*pi, n)/pi, np.real(F_all))
        plt.xlabel('k_x (pi)')
        plt.ylabel('Berry curvature')
        if band==0:
            plt.title('Valence Band')
        else:
            plt.title('Conductance Band')
        plt.show()
    end_time = time.time()
    print('运行时间(min)=', (end_time-start_time)/60)

def find_vector_with_the_same_gauge(vector_1, vector_0):
    # 寻找近似的同一的规范
    phase_1_pre = 0
    phase_2_pre = pi
    n_test = 10001
    for i0 in range(n_test):
        test_1 = np.sum(np.abs(vector_1*cmath.exp(1j*phase_1_pre) - vector_0))
        test_2 = np.sum(np.abs(vector_1*cmath.exp(1j*phase_2_pre) - vector_0))
        if test_1 < 1e-9:
            phase = phase_1_pre
            # print('Done with i0=', i0)
            break
        if i0 == n_test-1:
            phase = phase_1_pre
            print('Gauge Not Found with i0=', i0)
        if test_1 < test_2:
            if i0 == 0:
                phase_1 = phase_1_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2
            else:
                phase_1 = phase_1_pre
                phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2
        else:
            if i0 == 0:
                phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_2_pre+(phase_2_pre-phase_1_pre)/2
            else:
                phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_2_pre 
        phase_1_pre = phase_1
        phase_2_pre = phase_2
           
    vector_1 = vector_1*cmath.exp(1j*phase)
    # print('二分查找找到的规范=', phase)   
    return vector_1

if __name__ == '__main__':
    main()

计算结果为:

2. 考虑二维k空间的贝里曲率分布

代码为:

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

import numpy as np
import matplotlib.pyplot as plt
from math import *  
import cmath
import time


def hamiltonian(k1, k2, t1=2.82, a=1/sqrt(3)):  # 石墨烯哈密顿量(a为原子间距,不赋值的话默认为1/sqrt(3))
    h = np.zeros((2, 2), dtype=complex)
    h[0, 0] = 0.28/2
    h[1, 1] = -0.28/2
    h[1, 0] = t1*(cmath.exp(1j*k2*a)+cmath.exp(1j*sqrt(3)/2*k1*a-1j/2*k2*a)+cmath.exp(-1j*sqrt(3)/2*k1*a-1j/2*k2*a))
    h[0, 1] = h[1, 0].conj()
    return h


def main():
    start_time = time.time()
    n = 400  # 取点密度
    delta = 1e-10  # 求导的偏离量
    kx_array = np.linspace(-2*pi, 2*pi, n)
    ky_array = np.linspace(-2*pi, 2*pi, n)
    for band in range(2):
        F_all = np.zeros((ky_array.shape[0], ky_array.shape[0]))  # 贝里曲率
        j0 = 0
        for kx in kx_array:
            print(kx)
            i0 = 0
            for ky in ky_array:
                H = hamiltonian(kx, ky)
                eigenvalue, eigenvector = np.linalg.eig(H)
                vector = eigenvector[:, np.argsort(np.real(eigenvalue))[band]]  # 价带波函数
                # print(np.argsort(np.real(eigenvalue))[0])  # 排序索引(从小到大)
                # print(eigenvalue)  # 排序前的本征值
                # print(np.sort(np.real(eigenvalue)))  # 排序后的本征值(从小到大)
            
                H_delta_kx = hamiltonian(kx+delta, ky) 
                eigenvalue, eigenvector = np.linalg.eig(H_delta_kx)
                vector_delta_kx = eigenvector[:, np.argsort(np.real(eigenvalue))[band]]   # 略偏离kx的波函数
                # vector_delta_kx = find_vector_with_the_same_gauge(vector_delta_kx, vector)  # 如果波函数不连续需要使用这个

                H_delta_ky = hamiltonian(kx, ky+delta)  
                eigenvalue, eigenvector = np.linalg.eig(H_delta_ky)
                vector_delta_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[band]]  # 略偏离ky的波函数
                # vector_delta_ky = find_vector_with_the_same_gauge(vector_delta_ky, vector)  # 如果波函数不连续需要使用这个

                H_delta_kx_ky = hamiltonian(kx+delta, ky+delta)  
                eigenvalue, eigenvector = np.linalg.eig(H_delta_kx_ky)
                vector_delta_kx_ky = eigenvector[:, np.argsort(np.real(eigenvalue))[band]]  # 略偏离kx和ky的波函数
                # vector_delta_kx_ky = find_vector_with_the_same_gauge(vector_delta_kx_ky, vector)  # 如果波函数不连续需要使用这个

                # 价带的波函数的贝里联络(berry connection) # 求导后内积
                A_x = np.dot(vector.transpose().conj(), (vector_delta_kx-vector)/delta)   # 贝里联络Ax(x分量)
                A_y = np.dot(vector.transpose().conj(), (vector_delta_ky-vector)/delta)   # 贝里联络Ay(y分量)
                
                A_x_delta_ky = np.dot(vector_delta_ky.transpose().conj(), (vector_delta_kx_ky-vector_delta_ky)/delta)  # 略偏离ky的贝里联络Ax
                A_y_delta_kx = np.dot(vector_delta_kx.transpose().conj(), (vector_delta_kx_ky-vector_delta_kx)/delta)  # 略偏离kx的贝里联络Ay

                # 贝里曲率(berry curvature)
                F = ((A_y_delta_kx-A_y)/delta-(A_x_delta_ky-A_x)/delta)*1j
                # print(F)
                F_all[i0, j0] = np.real(F)
                i0 += 1
            j0 += 1

        if band==0:
            plot_3d_surface(kx_array/pi, ky_array/pi, F_all, xlabel='k_x (pi)', ylabel='k_y (pi)', zlabel='Berry curvature', title='Valence Band', rcount=300, ccount=300)
        else:
            plot_3d_surface(kx_array/pi, ky_array/pi, F_all, xlabel='k_x (pi)', ylabel='k_y (pi)', zlabel='Berry curvature', title='Conductance Band', rcount=300, ccount=300)
        # import guan
        # if band==0:
        #     guan.plot_3d_surface(kx_array/pi, ky_array/pi, F_all, xlabel='kx', ylabel='ky', zlabel='Berry curvature', title='Valence Band', rcount=300, ccount=300)
        # else:
        #     guan.plot_3d_surface(kx_array/pi, ky_array/pi, F_all, xlabel='kx', ylabel='ky', zlabel='Berry curvature', title='Conductance Band', rcount=300, ccount=300)
    end_time = time.time()
    print('运行时间(min)=', (end_time-start_time)/60)

def find_vector_with_the_same_gauge(vector_1, vector_0):
    # 寻找近似的同一的规范
    phase_1_pre = 0
    phase_2_pre = pi
    n_test = 10001
    for i0 in range(n_test):
        test_1 = np.sum(np.abs(vector_1*cmath.exp(1j*phase_1_pre) - vector_0))
        test_2 = np.sum(np.abs(vector_1*cmath.exp(1j*phase_2_pre) - vector_0))
        if test_1 < 1e-9:
            phase = phase_1_pre
            # print('Done with i0=', i0)
            break
        if i0 == n_test-1:
            phase = phase_1_pre
            print('Gauge Not Found with i0=', i0)
        if test_1 < test_2:
            if i0 == 0:
                phase_1 = phase_1_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2
            else:
                phase_1 = phase_1_pre
                phase_2 = phase_1_pre+(phase_2_pre-phase_1_pre)/2
        else:
            if i0 == 0:
                phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_2_pre+(phase_2_pre-phase_1_pre)/2
            else:
                phase_1 = phase_2_pre-(phase_2_pre-phase_1_pre)/2
                phase_2 = phase_2_pre 
        phase_1_pre = phase_1
        phase_2_pre = phase_2
           
    vector_1 = vector_1*cmath.exp(1j*phase)
    # print('二分查找找到的规范=', phase)   
    return vector_1

def plot_3d_surface(x_array, y_array, matrix, xlabel='x', ylabel='y', zlabel='z', title='', fontsize=20, labelsize=15, show=1, save=0, filename='a', file_format='.jpg', dpi=300, z_min=None, z_max=None, rcount=100, ccount=100): 
    import matplotlib.pyplot as plt
    from matplotlib import cm
    from matplotlib.ticker import LinearLocator
    matrix = np.array(matrix)
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
    plt.subplots_adjust(bottom=0.1, right=0.65) 
    x_array, y_array = np.meshgrid(x_array, y_array)
    if len(matrix.shape) == 2:
        surf = ax.plot_surface(x_array, y_array, matrix, rcount=rcount, ccount=ccount, cmap=cm.coolwarm, linewidth=0, antialiased=False) 
    elif len(matrix.shape) == 3:
        for i0 in range(matrix.shape[2]):
            surf = ax.plot_surface(x_array, y_array, matrix[:,:,i0], rcount=rcount, ccount=ccount, cmap=cm.coolwarm, linewidth=0, antialiased=False) 
    ax.set_title(title, fontsize=fontsize, fontfamily='Times New Roman')
    ax.set_xlabel(xlabel, fontsize=fontsize, fontfamily='Times New Roman') 
    ax.set_ylabel(ylabel, fontsize=fontsize, fontfamily='Times New Roman') 
    ax.set_zlabel(zlabel, fontsize=fontsize, fontfamily='Times New Roman') 
    ax.zaxis.set_major_locator(LinearLocator(5)) 
    ax.zaxis.set_major_formatter('{x:.2f}')  
    if z_min!=None or z_max!=None:
        if z_min==None:
            z_min=matrix.min()
        if z_max==None:
            z_max=matrix.max()
        ax.set_zlim(z_min, z_max)
    ax.tick_params(labelsize=labelsize) 
    labels = ax.get_xticklabels() + ax.get_yticklabels() + ax.get_zticklabels()
    [label.set_fontname('Times New Roman') for label in labels] 
    cax = plt.axes([0.8, 0.1, 0.05, 0.8]) 
    cbar = fig.colorbar(surf, cax=cax)  
    cbar.ax.tick_params(labelsize=labelsize)
    for l in cbar.ax.yaxis.get_ticklabels():
        l.set_family('Times New Roman')
    if save == 1:
        plt.savefig(filename+file_format, dpi=dpi) 
    if show == 1:
        plt.show()
    plt.close('all')

if __name__ == '__main__':
    main()

计算结果为:

2,070 次浏览

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

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

10 thoughts on “空间反演对称性破缺的石墨烯的贝里曲率分布(附Python代码)”

  1. 关老师您好,请问贝利曲率有单位吗,我看文献上单位都写的是Bohr^2,但是通过贝利曲率的公式感觉是对波函数的求导以及内积,没看出有单位是什么,您能帮忙解答一下吗

    1. 贝里曲率的单位确实是长度的平方,只是平时会比较少去提及。有两个角度可以理解:(1)陈数是没有单位的,陈数为贝里曲率对布里渊区的面积分,因为倒空间坐标的单位是长度的倒数,所以贝里曲率的单位就是长度的平方。需要注意的是:用上波尔长度的平方(Bohr^2),数值上是需要进行调整的。(2)不妨把波函数想象为平面波exp(ika),贝里联络中需要波函数对k求导,于是会出现一个长度单位;贝里曲率需要对贝里联络进行求导,还会跑出一个长度单位,所以贝里曲率的单位是长度的平方。

  2. 博主你好,我很好奇,Valley-Contrasting Physics in Graphene: Magnetic Moment and Topological Transport文章哈密顿量中的赝自旋系数是怎样作用在Hamilton
    量上的。可以看出,数值计算的模型应该是一个加了质量想的石墨烯模型,并没有添加t_z项

    1. 嗯,就只是加了一个质量项。你说的t_z项是什么,石墨烯是一个二维的体系,没有z方向的跃迁。

    1. 暂时没有,语义差不多的,可以自己翻译过去。或者根据定义直接写。

  3. 博主你好,我想问下,理论和计算凝聚态物理真的需要学python吗?是不是代码都是用python写的?

    1. 不是呀,其实凝聚态物理中很多代码都是用Fortran和Matlab写的,也有用C语言的。主要还是看写的内容,跟编程语言无关,想用哪个语言都可以。Fortran和C语言的运行效率会比较高一些。Python+NumPy和Matlab差不多的,Matlab是商业软件,Python是开源的。Python运行效率没有Fortran高,我主要是图书写方便不用定义变量,各种库比较多方便调用,又是开源的,而且最近因为机器学习Python语言也比较火,然后我就用了。

发表回复

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