这是之前的一篇：贝里曲率的计算（附Python代码）。本篇把其中的“高效法”和“Wilson loop方法”的贝里曲率计算代码改写成函数形式。此外加上这两篇文章的内容：陈数Chern number的计算（多条能带的高效法，附Python代码）、陈数Chern number的计算（多条能带的Wilson loop方法，附Python代码），多条能带同时计算可以支持能带交叉或简并的情况。
还是以这篇为例子：空间反演对称性破缺的石墨烯的贝里曲率分布（附Python代码）。计算结果和之前的相同。程度的部分函数是复制自Guan软件包的源码。为了使代码更加简洁，可以删除函数代码，直接调用Guan函数。官网：https://py.guanjihuan.com。安装方法：pip install --upgrade guan
补充：Kubo公式也可以推广到简并的情况，计算方法可以看这篇综述文献的公式（73）：First-principle calculations of the Berry curvature of Bloch states for charge and spin transport of electrons。该文献由评论区提供。这里暂时没考虑这个公式。
一、贝里曲率高效法（写成函数形式）
"""
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/24059
"""
import numpy as np
from math import *
import cmath
import math
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():
k_array, berry_curvature_array = calculate_berry_curvature_with_efficient_method(hamiltonian_function=hamiltonian, k_min=-2*math.pi, k_max=2*math.pi, precision=500, print_show=0)
plot_3d_surface(k_array, k_array, np.real(berry_curvature_array[:, :, 0]), title='Valence Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
plot_3d_surface(k_array, k_array, np.real(berry_curvature_array[:, :, 1]), title='Conductance Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
dim = berry_curvature_array.shape
plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :, 0]), title='Valence Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :, 1]), title='Conductance Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
# import guan
# k_array, berry_curvature_array = guan.calculate_berry_curvature_with_efficient_method(hamiltonian_function=hamiltonian, k_min=-2*math.pi, k_max=2*math.pi, precision=500, print_show=0)
# guan.plot_3d_surface(k_array, k_array, np.real(berry_curvature_array[:, :, 0]), title='Valence Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
# guan.plot_3d_surface(k_array, k_array, np.real(berry_curvature_array[:, :, 1]), title='Conductance Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
# dim = berry_curvature_array.shape
# guan.plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :, 0]), title='Valence Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
# guan.plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :, 1]), title='Conductance Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
def calculate_berry_curvature_with_efficient_method(hamiltonian_function, k_min=-math.pi, k_max=math.pi, precision=100, print_show=0):
if np.array(hamiltonian_function(0, 0)).shape==():
dim = 1
else:
dim = np.array(hamiltonian_function(0, 0)).shape[0]
delta = (k_max-k_min)/precision
k_array = np.arange(k_min, k_max, delta)
berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0], dim), dtype=complex)
i0 = 0
for kx in k_array:
if print_show == 1:
print(kx)
j0 = 0
for ky in k_array:
H = hamiltonian_function(kx, ky)
eigenvalue, vector = np.linalg.eigh(H)
H_delta_kx = hamiltonian_function(kx+delta, ky)
eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx)
H_delta_ky = hamiltonian_function(kx, ky+delta)
eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky)
H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta)
eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky)
for i in range(dim):
vector_i = vector[:, i]
vector_delta_kx_i = vector_delta_kx[:, i]
vector_delta_ky_i = vector_delta_ky[:, i]
vector_delta_kx_ky_i = vector_delta_kx_ky[:, i]
Ux = np.dot(np.conj(vector_i), vector_delta_kx_i)/abs(np.dot(np.conj(vector_i), vector_delta_kx_i))
Uy = np.dot(np.conj(vector_i), vector_delta_ky_i)/abs(np.dot(np.conj(vector_i), vector_delta_ky_i))
Ux_y = np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_ky_i), vector_delta_kx_ky_i))
Uy_x = np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i)/abs(np.dot(np.conj(vector_delta_kx_i), vector_delta_kx_ky_i))
berry_curvature = cmath.log(Ux*Uy_x*(1/Ux_y)*(1/Uy))/delta/delta*1j
berry_curvature_array[j0, i0, i] = berry_curvature
j0 += 1
i0 += 1
return k_array, berry_curvature_array
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')
def plot(x_array, y_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style='', y_min=None, y_max=None, linewidth=None, markersize=None, adjust_bottom=0.2, adjust_left=0.2):
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
plt.subplots_adjust(bottom=adjust_bottom, left=adjust_left)
ax.grid()
ax.tick_params(labelsize=labelsize)
labels = ax.get_xticklabels() + ax.get_yticklabels()
[label.set_fontname('Times New Roman') for label in labels]
ax.plot(x_array, y_array, style, linewidth=linewidth, markersize=markersize)
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')
if y_min!=None or y_max!=None:
if y_min==None:
y_min=min(y_array)
if y_max==None:
y_max=max(y_array)
ax.set_ylim(y_min, y_max)
if save == 1:
plt.savefig(filename+file_format, dpi=dpi)
if show == 1:
plt.show()
plt.close('all')
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/24059
"""
import numpy as np
from math import *
import cmath
import math
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():
k_array, berry_curvature_array = calculate_berry_curvature_with_efficient_method_for_degenerate_case(hamiltonian_function=hamiltonian, index_of_bands=[0], k_min=-2*math.pi, k_max=2*math.pi, precision=500)
dim = berry_curvature_array.shape
plot_3d_surface(k_array, k_array, np.real(berry_curvature_array), title='Valence Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :]), title='Valence Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
k_array, berry_curvature_array = calculate_berry_curvature_with_efficient_method_for_degenerate_case(hamiltonian_function=hamiltonian, index_of_bands=[0, 1], k_min=-2*math.pi, k_max=2*math.pi, precision=500)
dim = berry_curvature_array.shape
plot_3d_surface(k_array, k_array, np.real(berry_curvature_array), title='All Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :]), title='All Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
# import guan
# k_array, berry_curvature_array = guan.calculate_berry_curvature_with_efficient_method_for_degenerate_case(hamiltonian_function=hamiltonian, index_of_bands=[0], k_min=-2*math.pi, k_max=2*math.pi, precision=500)
# dim = berry_curvature_array.shape
# guan.plot_3d_surface(k_array, k_array, np.real(berry_curvature_array), title='Valence Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
# guan.plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :]), title='Valence Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
# k_array, berry_curvature_array = guan.calculate_berry_curvature_with_efficient_method_for_degenerate_case(hamiltonian_function=hamiltonian, index_of_bands=[0, 1], k_min=-2*math.pi, k_max=2*math.pi, precision=500)
# dim = berry_curvature_array.shape
# guan.plot_3d_surface(k_array, k_array, np.real(berry_curvature_array), title='All Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
# guan.plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :]), title='All Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
def calculate_berry_curvature_with_efficient_method_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], k_min=-math.pi, k_max=math.pi, precision=100, print_show=0):
delta = (k_max-k_min)/precision
k_array = np.arange(k_min, k_max, delta)
berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0]), dtype=complex)
i00 = 0
for kx in np.arange(k_min, k_max, delta):
if print_show == 1:
print(kx)
j00 = 0
for ky in np.arange(k_min, k_max, delta):
H = hamiltonian_function(kx, ky)
eigenvalue, vector = np.linalg.eigh(H)
H_delta_kx = hamiltonian_function(kx+delta, ky)
eigenvalue, vector_delta_kx = np.linalg.eigh(H_delta_kx)
H_delta_ky = hamiltonian_function(kx, ky+delta)
eigenvalue, vector_delta_ky = np.linalg.eigh(H_delta_ky)
H_delta_kx_ky = hamiltonian_function(kx+delta, ky+delta)
eigenvalue, vector_delta_kx_ky = np.linalg.eigh(H_delta_kx_ky)
dim = len(index_of_bands)
det_value = 1
# first dot product
dot_matrix = np.zeros((dim , dim), dtype=complex)
i0 = 0
for dim1 in index_of_bands:
j0 = 0
for dim2 in index_of_bands:
dot_matrix[i0, j0] = np.dot(np.conj(vector[:, dim1]), vector_delta_kx[:, dim2])
j0 += 1
i0 += 1
dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
det_value = det_value*dot_matrix
# second dot product
dot_matrix = np.zeros((dim , dim), dtype=complex)
i0 = 0
for dim1 in index_of_bands:
j0 = 0
for dim2 in index_of_bands:
dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx[:, dim1]), vector_delta_kx_ky[:, dim2])
j0 += 1
i0 += 1
dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
det_value = det_value*dot_matrix
# third dot product
dot_matrix = np.zeros((dim , dim), dtype=complex)
i0 = 0
for dim1 in index_of_bands:
j0 = 0
for dim2 in index_of_bands:
dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_kx_ky[:, dim1]), vector_delta_ky[:, dim2])
j0 += 1
i0 += 1
dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
det_value = det_value*dot_matrix
# four dot product
dot_matrix = np.zeros((dim , dim), dtype=complex)
i0 = 0
for dim1 in index_of_bands:
j0 = 0
for dim2 in index_of_bands:
dot_matrix[i0, j0] = np.dot(np.conj(vector_delta_ky[:, dim1]), vector[:, dim2])
j0 += 1
i0 += 1
dot_matrix = np.linalg.det(dot_matrix)/abs(np.linalg.det(dot_matrix))
det_value= det_value*dot_matrix
berry_curvature = cmath.log(det_value)/delta/delta*1j
berry_curvature_array[j00, i00] = berry_curvature
j00 += 1
i00 += 1
return k_array, berry_curvature_array
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')
def plot(x_array, y_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style='', y_min=None, y_max=None, linewidth=None, markersize=None, adjust_bottom=0.2, adjust_left=0.2):
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
plt.subplots_adjust(bottom=adjust_bottom, left=adjust_left)
ax.grid()
ax.tick_params(labelsize=labelsize)
labels = ax.get_xticklabels() + ax.get_yticklabels()
[label.set_fontname('Times New Roman') for label in labels]
ax.plot(x_array, y_array, style, linewidth=linewidth, markersize=markersize)
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')
if y_min!=None or y_max!=None:
if y_min==None:
y_min=min(y_array)
if y_max==None:
y_max=max(y_array)
ax.set_ylim(y_min, y_max)
if save == 1:
plt.savefig(filename+file_format, dpi=dpi)
if show == 1:
plt.show()
plt.close('all')
if __name__ == '__main__':
main()三、贝里曲率Wilson loop方法（写成函数形式）
"""
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/24059
"""
import numpy as np
from math import *
import cmath
import math
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():
k_array, berry_curvature_array = calculate_berry_curvature_with_wilson_loop(hamiltonian_function=hamiltonian, k_min=-2*math.pi, k_max=2*math.pi, precision_of_plaquettes=500, precision_of_wilson_loop=1)
plot_3d_surface(k_array, k_array, np.real(berry_curvature_array[:, :, 0]), title='Valence Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
plot_3d_surface(k_array, k_array, np.real(berry_curvature_array[:, :, 1]), title='Conductance Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
dim = berry_curvature_array.shape
plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :, 0]), title='Valence Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :, 1]), title='Conductance Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
# import guan
# k_array, berry_curvature_array = guan.calculate_berry_curvature_with_wilson_loop(hamiltonian_function=hamiltonian, k_min=-2*math.pi, k_max=2*math.pi, precision_of_plaquettes=500, precision_of_wilson_loop=1)
# guan.plot_3d_surface(k_array, k_array, np.real(berry_curvature_array[:, :, 0]), title='Valence Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
# guan.plot_3d_surface(k_array, k_array, np.real(berry_curvature_array[:, :, 1]), title='Conductance Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
# dim = berry_curvature_array.shape
# guan.plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :, 0]), title='Valence Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
# guan.plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :, 1]), title='Conductance Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
def calculate_berry_curvature_with_wilson_loop(hamiltonian_function, k_min=-math.pi, k_max=math.pi, precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0):
if np.array(hamiltonian_function(0, 0)).shape==():
dim = 1
else:
dim = np.array(hamiltonian_function(0, 0)).shape[0]
delta = (k_max-k_min)/precision_of_plaquettes
k_array = np.arange(k_min, k_max, delta)
berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0], dim), dtype=complex)
i00 = 0
for kx in k_array:
if print_show == 1:
print(kx)
j00 = 0
for ky in k_array:
vector_array = []
# line_1
for i0 in range(precision_of_wilson_loop):
H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky)
eigenvalue, eigenvector = np.linalg.eig(H_delta)
vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
vector_array.append(vector_delta)
# line_2
for i0 in range(precision_of_wilson_loop):
H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0)
eigenvalue, eigenvector = np.linalg.eig(H_delta)
vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
vector_array.append(vector_delta)
# line_3
for i0 in range(precision_of_wilson_loop):
H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta)
eigenvalue, eigenvector = np.linalg.eig(H_delta)
vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
vector_array.append(vector_delta)
# line_4
for i0 in range(precision_of_wilson_loop):
H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0)
eigenvalue, eigenvector = np.linalg.eig(H_delta)
vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
vector_array.append(vector_delta)
wilson_loop = 1
for i0 in range(len(vector_array)-1):
wilson_loop = wilson_loop*np.dot(vector_array[i0].transpose().conj(), vector_array[i0+1])
wilson_loop = wilson_loop*np.dot(vector_array[len(vector_array)-1].transpose().conj(), vector_array[0])
berry_curvature = np.log(np.diagonal(wilson_loop))/delta/delta*1j
berry_curvature_array[j00, i00, :]=berry_curvature
j00 += 1
i00 += 1
return k_array, berry_curvature_array
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')
def plot(x_array, y_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style='', y_min=None, y_max=None, linewidth=None, markersize=None, adjust_bottom=0.2, adjust_left=0.2):
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
plt.subplots_adjust(bottom=adjust_bottom, left=adjust_left)
ax.grid()
ax.tick_params(labelsize=labelsize)
labels = ax.get_xticklabels() + ax.get_yticklabels()
[label.set_fontname('Times New Roman') for label in labels]
ax.plot(x_array, y_array, style, linewidth=linewidth, markersize=markersize)
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')
if y_min!=None or y_max!=None:
if y_min==None:
y_min=min(y_array)
if y_max==None:
y_max=max(y_array)
ax.set_ylim(y_min, y_max)
if save == 1:
plt.savefig(filename+file_format, dpi=dpi)
if show == 1:
plt.show()
plt.close('all')
if __name__ == '__main__':
main()四、贝里曲率多条能带的Wilson loop方法（写成函数形式）
"""
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/24059
"""
import numpy as np
from math import *
import cmath
import math
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():
k_array, berry_curvature_array = calculate_berry_curvature_with_wilson_loop_for_degenerate_case(hamiltonian_function=hamiltonian, index_of_bands=[0], k_min=-2*math.pi, k_max=2*math.pi, precision_of_plaquettes=500, precision_of_wilson_loop=1)
dim = berry_curvature_array.shape
plot_3d_surface(k_array, k_array, np.real(berry_curvature_array), title='Valence Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :]), title='Valence Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
k_array, berry_curvature_array = calculate_berry_curvature_with_wilson_loop_for_degenerate_case(hamiltonian_function=hamiltonian, index_of_bands=[0, 1], k_min=-2*math.pi, k_max=2*math.pi, precision_of_plaquettes=500, precision_of_wilson_loop=1)
dim = berry_curvature_array.shape
plot_3d_surface(k_array, k_array, np.real(berry_curvature_array), title='All Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :]), title='All Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
# import guan
# k_array, berry_curvature_array = guan.calculate_berry_curvature_with_wilson_loop_for_degenerate_case(hamiltonian_function=hamiltonian, index_of_bands=[0], k_min=-2*math.pi, k_max=2*math.pi, precision_of_plaquettes=500, precision_of_wilson_loop=1)
# dim = berry_curvature_array.shape
# guan.plot_3d_surface(k_array, k_array, np.real(berry_curvature_array), title='Valence Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
# guan.plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :]), title='Valence Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
# k_array, berry_curvature_array = guan.calculate_berry_curvature_with_wilson_loop_for_degenerate_case(hamiltonian_function=hamiltonian, index_of_bands=[0, 1], k_min=-2*math.pi, k_max=2*math.pi, precision_of_plaquettes=500, precision_of_wilson_loop=1)
# dim = berry_curvature_array.shape
# guan.plot_3d_surface(k_array, k_array, np.real(berry_curvature_array), title='All Band', xlabel='kx', ylabel='ky', zlabel='Berry curvature')
# guan.plot(k_array, np.real(berry_curvature_array[int(dim[0]/2), :]), title='All Band ky=0', xlabel='kx', ylabel='Berry curvature') # ky=0
def calculate_berry_curvature_with_wilson_loop_for_degenerate_case(hamiltonian_function, index_of_bands=[0, 1], k_min=-math.pi, k_max=math.pi, precision_of_plaquettes=20, precision_of_wilson_loop=5, print_show=0):
delta = (k_max-k_min)/precision_of_plaquettes
k_array = np.arange(k_min, k_max, delta)
berry_curvature_array = np.zeros((k_array.shape[0], k_array.shape[0]), dtype=complex)
i000 = 0
for kx in k_array:
if print_show == 1:
print(kx)
j000 = 0
for ky in k_array:
vector_array = []
# line_1
for i0 in range(precision_of_wilson_loop):
H_delta = hamiltonian_function(kx+delta/precision_of_wilson_loop*i0, ky)
eigenvalue, eigenvector = np.linalg.eig(H_delta)
vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
vector_array.append(vector_delta)
# line_2
for i0 in range(precision_of_wilson_loop):
H_delta = hamiltonian_function(kx+delta, ky+delta/precision_of_wilson_loop*i0)
eigenvalue, eigenvector = np.linalg.eig(H_delta)
vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
vector_array.append(vector_delta)
# line_3
for i0 in range(precision_of_wilson_loop):
H_delta = hamiltonian_function(kx+delta-delta/precision_of_wilson_loop*i0, ky+delta)
eigenvalue, eigenvector = np.linalg.eig(H_delta)
vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
vector_array.append(vector_delta)
# line_4
for i0 in range(precision_of_wilson_loop):
H_delta = hamiltonian_function(kx, ky+delta-delta/precision_of_wilson_loop*i0)
eigenvalue, eigenvector = np.linalg.eig(H_delta)
vector_delta = eigenvector[:, np.argsort(np.real(eigenvalue))]
vector_array.append(vector_delta)
wilson_loop = 1
dim = len(index_of_bands)
for i0 in range(len(vector_array)-1):
dot_matrix = np.zeros((dim , dim), dtype=complex)
i01 = 0
for dim1 in index_of_bands:
i02 = 0
for dim2 in index_of_bands:
dot_matrix[i01, i02] = np.dot(vector_array[i0][:, dim1].transpose().conj(), vector_array[i0+1][:, dim2])
i02 += 1
i01 += 1
det_value = np.linalg.det(dot_matrix)
wilson_loop = wilson_loop*det_value
dot_matrix_plus = np.zeros((dim , dim), dtype=complex)
i01 = 0
for dim1 in index_of_bands:
i02 = 0
for dim2 in index_of_bands:
dot_matrix_plus[i01, i02] = np.dot(vector_array[len(vector_array)-1][:, dim1].transpose().conj(), vector_array[0][:, dim2])
i02 += 1
i01 += 1
det_value = np.linalg.det(dot_matrix_plus)
wilson_loop = wilson_loop*det_value
berry_curvature = np.log(wilson_loop)/delta/delta*1j
berry_curvature_array[j000, i000]=berry_curvature
j000 += 1
i000 += 1
return k_array, berry_curvature_array
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')
def plot(x_array, y_array, xlabel='x', ylabel='y', title='', fontsize=20, labelsize=20, show=1, save=0, filename='a', file_format='.jpg', dpi=300, style='', y_min=None, y_max=None, linewidth=None, markersize=None, adjust_bottom=0.2, adjust_left=0.2):
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
plt.subplots_adjust(bottom=adjust_bottom, left=adjust_left)
ax.grid()
ax.tick_params(labelsize=labelsize)
labels = ax.get_xticklabels() + ax.get_yticklabels()
[label.set_fontname('Times New Roman') for label in labels]
ax.plot(x_array, y_array, style, linewidth=linewidth, markersize=markersize)
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')
if y_min!=None or y_max!=None:
if y_min==None:
y_min=min(y_array)
if y_max==None:
y_max=max(y_array)
ax.set_ylim(y_min, y_max)
if save == 1:
plt.savefig(filename+file_format, dpi=dpi)
if show == 1:
plt.show()
plt.close('all')
if __name__ == '__main__':
main()【说明：本站主要是个人的一些笔记和代码分享，内容可能会不定期修改。为了使全网显示的始终是最新版本，这里的文章未经同意请勿转载。引用请注明出处：https://www.guanjihuan.com】
请问石墨烯的贝利曲率的积分面积为什么是[-pi, pi]？
这里没有积分，只算贝里曲率，所以只是大概算了一个区域。
博主您好，请问用kubo公式能求多能带的贝里曲率吗？
如果不是简并的情况，可以把直接多个带的贝里曲率求和。如果是简并的情况，用Kubo公式好像会出一些问题，我暂时还不知道怎么推广这个公式。
感谢回复，我想计算三维情况下的berry curvature。似乎kubo公式比较容易推广， 用高效法或wilson loop可以实现吗？
都是可以算的吧。我的理解是：在一个截面上算的贝里曲率只是三维贝里曲率的一个分量。
我在这篇文献下找到了简并情况的计算方法， 公式（73）。
https://refubium.fu-berlin.de/bitstream/handle/fub188/14194/Gradhand-et-al.pdf?sequence=1
嗯，不错，非常感谢！