动力学矩阵法计算石墨烯声子谱

phyzj

已于 2022-02-16 20:51:30 修改

阅读量3.1k

点赞数 5

文章标签： python

于 2022-01-26 22:33:50 首次发布

本文链接：https://blog.csdn.net/dashenthk/article/details/122709500

版权

上个学期，学习固体物理，老师让写一个程序来计算石墨烯的声子谱和振动的态密度进行画图。大二的假期里自学了Python但是没有实际自己写过，决定用Python代码来实现这个任务（最重要的是自己能力不行，c语言学的不行，用Python来实现真的很方便！！！）。最终竟然用了100多行才完成任务，也是第一次写这么长的代码，记录一下。~ ~ ~也希望能给下一届的学弟学妹们一个参考。

动力学矩阵方法计算石墨烯声子谱的计算步骤参考 Physical Propertier of Carbon Nanotubes（第九章163页开始） .这里只有利用动力学矩阵法计算石墨烯声子谱部分，振动态密度部分计算时间较长，没有找到较好的解决方法，就直接删去了。计算中只考虑碳原子的四近邻的相互作用。

"""
在x-y-z坐标系中
正格子基矢:a1 = [ 3 / 2,sqrt(3) / 2 ] * a , a2 = [ 0 , sqrt(3) ] * a
倒格子基矢：b1 = [ 4 * pi / 3 / a , 0 , 0], b2 = [-2*pi/3/a , 2*pi/a/sqrt(3),0]
"""
from math import sqrt, cos, sin, pi, acos
import matplotlib.pyplot as plt
import numpy as np

a = 1.42e-10    # 晶格常数
m = 1.99e-26     # 碳原子质量

# 生成力常数的对角矩阵

f = np.diag([36.50, 24.50, 9.82, 8.80, -3.23, -0.40, 3.00, -5.25, 0.15, -1.92, 2.29, -0.58])

# 定义旋转函数


def rotate(theta, r):
    """
        theta: 旋转角度
        r: 位置坐标或者K矩阵
        旋转之后通过旋转矩阵U可以将以A类原子为中心的近邻原子的相应的K矩阵转化成以B类原子为中心的
    """
    R = []
    RB = []
    Um = np.array([[cos(theta), sin(theta), 0], [-sin(theta), cos(theta), 0], [0, 0, 1]])
    U = np.array([[cos(pi), sin(pi), 0], [-sin(pi), cos(pi), 0], [0, 0, 1]])

    for i in range(int(2*pi / theta)):    # 确定旋转次数
        if np.size(r) == 3:               # 通过判断来确定是位置坐标还是K矩阵
            r = np.matmul(np.linalg.inv(Um), r)    # 位置坐标旋转
            rb = r*(-1)                            # 中心原子改变后的位置
        else:
            r = np.linalg.inv(Um)@r@Um            # 对K矩阵的操作，同位置坐标的操作类似
            rb = np.linalg.inv(U)@r@U
        R.append(r)                       # 用一个列表储存，接下来直接求和，计算方便
        RB.append(rb)                     # 储存B原子为中心的位置和K矩阵
    return R, RB


# 以A为圆心周围的原子的坐标 以及 B的相应坐标
FA, FB = rotate(2/3*pi, np.array([[1], [0], [0]])*a)          # 第一近邻
SA, SB = rotate(1/3*pi, np.array([[3/2], [sqrt(3)/2], [0]])*a)     # 第二近邻
TA, TB = rotate(2/3*pi, np.array([[1], [sqrt(3)], [0]])*a)         # 第三近邻
LA1, LB1 = rotate(2/3*pi, np.array([[2.5], [sqrt(3)/2], [0]])*a)   # 第四近邻，两种角度
LA2, LB2 = rotate(2/3*pi, np.array([[2.5], [-sqrt(3)/2], [0]])*a)
LA = LA1+LA2                                                       # 第四近邻两种情况进行合并，可省略
LB = LB1+LB2


def K(theta, k):
    """
    构建初始的Kij，中心原子与正右侧的原子间的K矩阵为对角矩阵，通过旋转，得到中心原子与每个近邻的第一个原子间的K矩阵
    """
    U = np.array([[cos(theta), sin(theta), 0], [-sin(theta), cos(theta), 0], [0, 0, 1]])
    return np.linalg.inv(U)@k@U


# 得到初始K矩阵，角度由文献的图9.1 可知，分别得到（a）图和(b)图情况的初始K矩阵
KAB1, KBA1 = rotate(2/3*pi, f[0:3, 0:3])
KAA, KBB = rotate(1/3*pi, K(1/6*pi, f[3:6, 3:6]))
KAB3, KBA3 = rotate(2/3*pi, K(1/3*pi, f[6:9, 6:9]))
KAB4f, KBA4f = rotate(2/3*pi, K(acos(2.5/sqrt(7)), f[9:12, 9:12]))
KAB4s, KBA4s = rotate(2/3*pi, K(2*pi-acos(2.5/sqrt(7)), f[9:12, 9:12]))
KAB4 = KAB4f+KAB4s
KBA4 = KBA4f+KBA4s


# 构造D矩阵
def Dm(k):
    D = np.zeros([6, 6], dtype=complex)
    DAAs = 0
    DBBs = 0
    DBAs = 0
    DABs = 0
# 以上四个矩阵由文献的 9.7 式可知，并组合成一个完整的D矩阵
    for i in range(3):
        DAAs = DAAs+KAA[i]*np.exp(1j*np.matmul(k, -SA[i])) + KAA[i+3]*np.exp(1j*np.matmul(k, -SA[i+3]))
        DBBs = DBBs + KBB[i] * np.exp(1j * np.matmul(k, -SB[i])) + KBB[i+3]*np.exp(1j*np.matmul(k, -SB[i+3]))
        DABs = DABs + KAB1[i] * np.exp(1j * np.matmul(k, -FA[i])) + KAB3[i] * np.exp(1j * np.matmul(k, -TA[i])) +\
               KAB4[i] * np.exp(1j * np.matmul(k, -LA[i])) + KAB4[i + 3] * np.exp(1j * np.matmul(k, -LA[i + 3]))
        DBAs = DBAs + KBA1[i] * np.exp(1j * np.matmul(k, -FB[i])) + KBA3[i] * np.exp(1j * np.matmul(k, -TB[i])) + \
               KBA4[i] * np.exp(1j * np.matmul(k, -LB[i])) + KBA4[i + 3] * np.exp(1j * np.matmul(k, -LB[i + 3]))
    D[0:3, 3:6] = -DABs
    D[3:6, 0:3] = -DBAs
    D[0:3, 0:3] = sum(KAB1)+sum(KAA)+sum(KAB3)+sum(KAB4f)+sum(KAB4s)-DAAs
    D[3:6, 3:6] = sum(KAB1) + sum(KAA) + sum(KAB3) + sum(KAB4f) + sum(KAB4s) - DBBs

    return D


# 求解矩阵
n = 100         # 沿第一布里渊区的最小重复单元的边界每条边单位长度取100个k点进行计算（主要是为了使点均匀分布，同样可以每条边取n 个点进行计算）
result = np.zeros([(30 + int(sqrt(3) * 10)) * n, 6])    # 解的矩阵
for i in range((30 + int(sqrt(3) * 10)) * n):           # 在这里将sqrt(3)近似取为17，没有什么特别的意义
    if i < n * int(10 * sqrt(3)):                       # 判断i的大小确定k的取值
        kx = i*2*pi/3/a/(n * int(10 * sqrt(3)))
        ky = 0
    elif i < (10 + int(10 * sqrt(3))) * n:
        kx = 2 * pi / 3 / a
        ky = (i-n * int(10 * sqrt(3))) / (10 * n-1) * 2 * pi / 3 / a/sqrt(3)
    else:
        kx = 2 * pi / 3 / a - (i-(10 + int(10 * sqrt(3))) * n)/(n * 20-1) * 2 * pi / 3 / a
        ky = kx/sqrt(3)
    k = np.array([kx, ky, 0])                 # 得到k值，带入D矩阵
    w, t = np.linalg.eig(Dm(k))
    w = list(w)
    w.sort()
    result[i, :] = (np.real(np.sqrt(w) / sqrt(m)))      # 将本征值进行保存


# 画图
s = sqrt(3)
xk = [0, s, s+1, s+3]
kk = np.linspace(0, 4.7, num=(30 + int(sqrt(3) * 10)) * n)   # 作为x轴，使其和本征值矩阵每一列的y的值个数相同
plt.plot(kk, result, c="r")
plt.xticks(xk, ["Γ", "M", "K", "Γ"])
plt.xlim(0, 4.75)
plt.ylim(0, 1e14)
plt.ylabel("ω", fontsize=14)
plt.show()