新建
2024秋计算材料学-上机练习：ABACUS能带和态密度计算
stanfordbshan
FermiNoodles
推荐镜像 :comput-mater:0512v2
推荐机型 :c2_m4_cpu
赞 3
14
4
目录
2024秋计算材料学-上机练习：ABACUS能带和态密度计算
代码
文本
©️ Copyright 2024 @ Authors
作者：Jinghang Wang (AISI电子结构团队实习生), 斯坦福大厨 📨
日期：2024-05-28
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的 开始连接 按钮，选择 comput-mater:0512v2 镜像及 c2_m4_cpu 节点配置即可开始运行。
上机练习：为了使你更深入地理解计算的流程并体会动手完成计算的快乐，在上机练习中我们设置了一些互动环节，部分单元格需要你来填写。
恭喜您已经发现了这份神奇的计算材料学课件！这份课件是我在熬夜吃掉不计其数的披萨和咖啡后创作出来的，配套的教材是由单斌、陈征征、陈蓉合著的《计算材料学--从算法原理到代码实现》。学习资料合集您可以在这个网址找到：www.materialssimulation.com/book，您也可以跟着up主无人问津晦涩难懂的B站视频一起进行学习。希望它能帮您在计算材料学的道路上摔得不那么痛。
就像您尊重那些一边烘焙披萨一边写代码的大厨一样，当您使用这份课件时，请：
- 记得告诉大家这份课件是斯坦福大厨写的，并且他在华中科技大学微纳中心工作
- 别用它去赚大钱，这个课件是用来学习的，不是用来买披萨的
- 保持开放共享的精神
如果你有关于计算材料学的想法，或者你只是想和我讨论最好吃的披萨口味，欢迎通过邮件 bshan@mail.hust.edu.cn 联系我。
代码
文本
代码
文本
1. 引言
代码
文本
这里是2024秋计算材料学-上机练习第一课——ABACUS能带和态密度计算。相比于以往的课程，我们在该系列中设置了更多互动环节，使你在更深入地理解相应的概念的同时又不会感觉到枯燥。
代码
文本
我们引用保罗·狄拉克的一段话作为计算材料学系列课程的开端：
The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. It therefore becomes desirable that approximate practical methods of applying quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computation.
|
这段话首先指出物理和化学领域的数学理论所需的物理定律已经完备，困难在于求解这些复杂的方程。接着提到有必要发展量子力学的近似实用方法，这些方法在不需要太多计算的情况下，可以解释复杂原子系统的主要特征。
代码
文本
密度泛函理论就是应用最为广泛的量子力学近似方法之一。密度泛函理论是不依赖经验参数的第一性原理计算方法，最终目的是精确求解多体薛定谔方程，从而获取材料的各种物理性质。
代码
文本
ABACUS（Atomic-orbtial Based Ab-initio Computation at UStc，中文名原子算筹）是一款国产开源密度泛函理论软件，由中国科学技术大学量子信息与超级计算中心重点实验室-中国科学院计算机网络与信息中心开发，旨在从第一性原理出发进行大规模电子结构模拟。ABACUS目前与DeepModeling社区合作，在北京科学智能研究院、物理所、北京大学等高校和研究院所均有团队开发人员。
代码
文本
在本节课程中，我们将通过石墨烯和硅这两种材料案例，简单介绍晶体结构、倒空间以及电子能带这些概念，并且使用ABACUS完成相关第一性原理计算。
代码
文本
2. 石墨烯的能带和态密度计算
代码
文本
石墨烯是一种由单层碳原子密集排列成二维蜂窝状晶格的材料。2004年，英国曼彻斯特大学的安德烈·海姆 (Andre Geim)和康斯坦丁·诺沃肖洛夫 (Konstantin Novoselov)使用透明胶带对石墨进行了分离，得到了在之前被认为是不可能存在的石墨烯。他们因这一开创性的工作而在2010年获得了诺贝尔物理学奖。他们的发表以及他们所描述的令人惊讶的简单制备方法引发了“石墨烯淘金热”。研究扩展到许多不同的子领域，探索材料的不同特殊特性——量子力学、电学、化学、机械、光学、磁学等。
石墨烯的特点之一便是其电子的能量和动量具有线性的色散关系，因此其费米能级附近的电子能带结构呈现出上下两个锥体，分别代表电子和空穴。两个锥体的顶端刚好相连，形成“零带隙”的半金属相，这显著影响了其电子属性，包括高电导率和异常的载流子迁移率。态密度（DOS）在石墨烯中也具有独特性。它在狄拉克点处消失，并且从这些点向外线性增加。这影响了电子和空穴（电子空位）对石墨烯导电性的贡献。
代码
文本
使用紧束缚模型可以得到蜂巢格子的色散关系：
在 K/K' 附近，色散关系 ε(q) 变成：
即在 K 点附近，色散关系是线性的，因此态密度 ρ(ε) ∝ ε，也是线性的。
图1 石墨烯的晶体结构、能带结构以及狄拉克点附近的色散关系
（来源: Understanding the Electrical Conductivity of Graphene ）
（来源: Understanding the Electrical Conductivity of Graphene ）
代码
文本
接下来我们就使用第一性原理软件ABACUS计算石墨稀的能带结构以及态密度，并可视化输出结果。
代码
文本
2.1 石墨烯的原子结构
代码
文本
对于任何材料体系的计算，建模都是第一步。对于任何材料体系的计算，建模都是第一步。石墨烯的碳原子排列与石墨的单原子层相同，是碳原子以sp2杂化轨道呈蜂巢晶格排列构成的单层二维晶体。但你并不需要真的构建一个包含上万个原子的石墨烯构型，因为凝聚态领域的第一性原理计算软件（例如VASP，ABACUS等）都会使用周期性边界条件。因此，下面你只需要构建一个包含两个碳原子的原胞。
代码
文本
(a)
(b)
图 2 石墨烯的蜂窝晶格，不同的颜色用于表示两个子晶格。基向量用黑色箭头表示
(a)基矢夹角为60° (b)基矢夹角为120°
(a)基矢夹角为60° (b)基矢夹角为120°
代码
文本
🔖接下来，为了更深刻地理解石墨烯的晶格结构，请你跟着提示做：
我们的目标是构建一个包含两个碳原子的原胞，其中原子的位置是相对于基矢的。但是石墨烯的基矢取法不唯一。因此，在图2中我们给出了两种基矢取法及对应的基向量，基矢的原点均在六方晶格的中心，夹角分别为60°和120°。在下方的代码单元格中，你可以通过'basis'参数选择一组基矢，然后填写该基矢下的A，B原子的相对坐标。
请在下方单元格中填写，确认后请点击单元格左上角的➡️运行单元格。
代码
文本
[1]
import os
import py3Dmol
from ase import Atoms
import numpy as np
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
# Fill in the cooridnate of two C atoms within the unit cell
# example
# atom1=[0.5, 0.4]
sqrt3 = np.sqrt(3)
# lattice parameter
a = 2.46
# choose a basis, basis = 'a'/'b'
basis = ''
# position of first atom, in fractional coordinates
atomA_frac = [ , ]
# position of second atom, in fractional coordinates
atomB_frac = [ , ]
代码
文本
填充完上一格中的晶胞基矢以及原子坐标后，请运行下面的单元格可视化石墨烯的结构，确保原子坐标的正确性。如果一切顺利，你将看到一个5x5的六方蜂窝型石墨烯结构。如果没有得到理想的结果，我们在下面给出了答案。
- 按住鼠标左键拖拽，可以旋转观察二维石墨烯结构
- 按住鼠标右键上下移动，可以进行视角缩放
代码
文本
[2]
import py3Dmol
from ase import Atoms
import numpy as np
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
# lattice parameter
a = 2.46
# create cell
if basis == 'a':
cell = [
[a*np.sqrt(3)/2, a*(1/2), 0],
[a*np.sqrt(3)/2, a*(-1/2), 0],
[0, 0, 25] # vacuum layer
]
elif basis == 'b':
cell = [
[-a*np.sqrt(3)/2, a*(1/2), 0],
[a*np.sqrt(3)/2, a*(1/2), 0],
[0, 0, 25] # vacuum layer
]
else:
print("Invalid basis.")
# Position of two carbon atoms in graphene (fractional coordinates)
atomA_frac_3 = [None, None, 0]
atomA_frac_3 = atomA_frac[:2] + atomA_frac_3[2:]
atomB_frac_3 = [None, None, 0]
atomB_frac_3 = atomB_frac[:2] + atomB_frac_3[2:]
# frac 2 cart
atomA = np.dot(atomA_frac_3, cell)
atomB = np.dot(atomB_frac_3, cell)
# build
graphene = Atoms('C2', positions=[atomA, atomB], cell=cell, pbc=True)
# extend
extended_graphene = graphene * (5, 5, 1)
# show
view = py3Dmol.view(width=800, height=400)
xyz_data = extended_graphene.get_positions()
symbols = extended_graphene.get_chemical_symbols()
for i, (x, y, z) in enumerate(xyz_data):
view.addSphere({'center': {'x': x, 'y': y, 'z': z}, 'radius': 0.4, 'color': 'gray'})
view.zoomTo()
view.show()
代码
文本
下一个单元格是答案
代码
文本
答案在此，点击展开
- basis = 'a'
atomA_frac = [1/3, 1/3]
atomB_frac = [2/3, 2/3]
- basis = 'b'
atomA_frac = [2/3, 1/3]
atomB_frac = [1/3, 2/3]
代码
文本
运行下方单元格中的代码可以将你建模的石墨烯结构转变为ABACUS可识别的STRU文件，从而使用你的结构完成2.4.2和2.4.3节的计算。你也可以跳过这一步，直接使用我们准备好的STRU文件。
代码
文本
已隐藏单元格
代码
文本
2.2 石墨烯的布里渊区
- 布里渊区：在倒格子空间中以任意一个倒格点为原点，做原点和其它所有倒格点连线的中垂面(或中垂线)，这些中垂面(或中垂线)将倒格子空间分割成许多区域，这些区域称为布里渊区(Brillouin zone)；把连接两个倒格点连线之间的垂直平分面称为布拉格平面。在倒格子空间中，以一个倒格点为原点，从原点出发，不经过任何布拉格平面所能到达的所有点的集合，称为第一布里渊区(first Brillouin zone)，也叫简约布里渊区。显然，它是围绕原点的最小闭合区域。由布里渊区的定义我们容易看出，第一布里渊区和前面所讲的魏格纳-塞茨(Wigner-Seitz)原胞的取法一样，所以通常人们把第一布里渊区定义为倒格子空间中的魏格纳-塞茨原胞。
根据布洛赫定理，简约波矢 k 表示原胞之间电子波函数位相的变化，如果 k 改变一个倒格矢量，它们所标志的原胞之间波函数位相的变化是相同的，也就是说E(k)是 k 空间的周期函数，其周期等于倒格矢。我们利用这种平移对称性可以将第二Brillouin区的每一块各自平移一个倒格矢而与第一Brillouin区重合。同理，更高的Brillouin区也可通过适当的平移与第一区重合，因此我们可以把注意力仅限制在第一区内，它包含了晶体能带的所有必要信息。
上面我们通过晶体的平移对称性将研究对象从整个倒空间约化为第一布里渊区。可以证明，在 k 空间中 En(k) 具有与晶体点群完全相同的对称性，这样就可以把第一布里渊区再分成若干个等价的小区域，即不可约布里渊区 (Irreducible Brillouin-Zone)，只取其中一个就足够了。区域大小为第一布里渊区的 1/f，f 为晶体点群对称操作元素数。 - 布里渊区中的高对称点：对于一般 k 点，简约区中对称相关的波矢量数就等于点群的阶数。但若 k 点在简约区中的某些特殊位置（对称点、对称轴或对称面）上，简约区中等价波矢量数就少于点群的阶数。这些 k 点就是能带计算中常用的高对称点（high symmetry points），又称临界点（critical points）。布里渊区的高对称点和线通常对应于由晶体对称性导致电子能带结构简化，这些点和线也通常是出现能带极值(最大值和最小值)的临界点，这对于理解电子性质(如带隙、有效质量和光学性质)至关重要。因此在实际能带计算中，我们会在布里渊区内用一条k点路径连接所有的高对称点作为自变量，得到的E-k关系就是我们熟悉的能带结构图。
代码
文本
🔖接下来，为了更深刻的理解布里渊区，请你跟着提示做：
- 二维的倒格矢的表达式我们已经给出。请根据上面的你所选择基矢，求出对应的倒格矢'b1'和'b2'，填写在下方的代码单元格中。（如果你想尝试不同基矢下的布里渊区，也可以在单元格中重新设定basis参数，然后求出对应的基矢及高对称点继续后续操作）
- 二维六方晶格的高对称点包括：G：布里渊区中心；K, K'：布里渊区的角点；M：布里渊区边界的中点。请根据你求出的倒格矢，在下方的代码单元格填写它们的相对坐标。
- 填写并运行完毕后，你可以运行再下一个代码单元格可视化布里渊区。
请在下方单元格中填写，确认后请点击单元格左上角的➡️运行单元格。
代码
文本
代码
文本
[19]
# basis = ''
# Define the reciprocal lattice vectors
b1 = 2 * np.pi / a * np.array([ , ])
b2 = 2 * np.pi / a * np.array([ , ])
# Define the high-symmetry points in the Brillouin zone
K_frac = [ , ]
K_prime_frac = [, ]
M_frac = [, ]
G_frac = [ , ]
代码
文本
填写并运行完毕后，你可以运行下方的代码单元格可视化布里渊区及对应的高对称点。如果一切顺利，你将看到一个六方型的格子，并且红色的高对称点位于其中心、角点以及边界的中点。如果没有得到理想的结果，我们在下面给出了答案。
代码
文本
[24]
import matplotlib.pyplot as plt
import numpy as np
# Plot the Brillouin zone
fig, ax = plt.subplots()
ax.plot([b1[0], b1[0] + b2[0]], [b1[1], b1[1] + b2[1]], 'k--')
ax.plot([b2[0], b1[0] + b2[0]], [b2[1], b1[1] + b2[1]], 'k--')
ax.quiver(0, 0, b1[0], b1[1], angles='xy', scale_units='xy', scale=1, color='k')
ax.quiver(0, 0, b2[0], b2[1], angles='xy', scale_units='xy', scale=1, color='k')
ax.text(b1[0]-0.2, b1[1]-0.2, 'b1', fontsize=12, color='k', ha='right')
ax.text(b2[0]-0.2, b2[1]+0.1, 'b2', fontsize=12, color='k', ha='right')
# frac 2 cart
K = K_frac[0]*b1 + K_frac[1]*b2
K_prime = K_prime_frac[0]*b1 + K_prime_frac[1]*b2
M = M_frac[0]*b1 + M_frac[1]*b2
G = G_frac[0]*b1 + G_frac[1]*b2
# Function to plot the Brillouin zone boundaries
if basis == 'a':
def plot_bz(ax, b1, b2):
BZ_K = (b1 + 2*b2) / 3
BZ_K_prime = (2*b1 + b2) / 3
BZ_M = (b1 + b2) / 2
BZ_G = np.array([0, 0])
corners = [
BZ_K_prime,
BZ_K,
(BZ_K-BZ_M)*2,
-BZ_K_prime,
-BZ_K,
-(BZ_K-BZ_M)*2,
BZ_K_prime
]
corners.append(corners[0]) # close the loop
# Extract x and y coordinates of corners
x = [corner[0] for corner in corners]
y = [corner[1] for corner in corners]
# Plot the Brillouin zone
ax.plot(x, y, 'k-')
elif basis == 'b':
def plot_bz(ax, b1, b2):
BZ_K = (b1 + b2) / 3
BZ_K_prime = ((-1)*b1 + 2*b2) / 3
BZ_M = b2 / 2
BZ_G = np.array([0, 0])
corners = [
BZ_K_prime,
BZ_K,
(BZ_K-BZ_M)*2,
-BZ_K_prime,
-BZ_K,
-(BZ_K-BZ_M)*2,
BZ_K_prime
]
corners.append(corners[0]) # close the loop
# Extract x and y coordinates of corners
x = [corner[0] for corner in corners]
y = [corner[1] for corner in corners]
# Plot the Brillouin zone
ax.plot(x, y, 'k-')
else:
print("Invalid basis.")
# Plot first Brillouin zone
plot_bz(ax, b1, b2)
# Plot high-symmetry points
ax.plot(G[0], G[1], 'ro')
ax.text(G[0], G[1], 'G', fontsize=12, ha='right')
ax.plot(K[0], K[1], 'ro')
ax.text(K[0], K[1], 'K', fontsize=12, ha='right')
ax.plot(K_prime[0], K_prime[1], 'ro')
ax.text(K_prime[0], K_prime[1], 'K\'', fontsize=12, ha='right')
ax.plot(M[0], M[1], 'ro')
ax.text(M[0], M[1], 'M', fontsize=12, ha='right')
ax.set_aspect('equal')
ax.set_xlabel('kx')
ax.set_ylabel('ky')
plt.title('Brillouin Zone of Graphene')
plt.grid(True)
plt.show()
代码
文本
下一个单元格是答案
代码
文本
答案在此，点击展开
basis = 'a'
b1 = 2 * np.pi / a * np.array([np.sqrt(3)/3, 1])
b2 = 2 * np.pi / a * np.array([np.sqrt(3)/3, -1])
K_frac = [ 1/3, 2/3 ]
K_prime_frac = [ 2/3 , 1/3 ]
M_frac = [ 1/2 , 1/2]
G_frac = [ 0, 0 ]basis = 'b'
b1 = 2 * np.pi / a * np.array([-np.sqrt(3)/3, 1])
b2 = 2 * np.pi / a * np.array([np.sqrt(3)/3, 1])
K_frac = [ 1/3, 1/3 ]
K_prime_frac = [ -1/3 , 2/3 ]
M_frac = [ 0 , 1/2]
G_frac = [ 0, 0 ]
代码
文本
2.3 紧束缚模型石墨烯色散关系
代码
文本
- 紧束缚模型：我们假定原子实对电子的束缚作用很强，因此，当电子距某个原子实比较近时，电子的运动主要受该原子势场的影响，受其它原子势场影响很弱。因此固体中电子的行为同孤立原子中电子的行为更为相似。这时可将孤立原子看成零级近似，而将其他原子势场的影响看成小的微扰，由此可以给出电子的原子能级和晶体能带之间的相互联系。
代码
文本
🔖接下来，为了得到紧束缚模型下的石墨烯色散关系，请你跟着提示做：
对于本节开始时给出紧束缚模型下的石墨烯色散关系，ε0=0且为正值时的表达式如下。请你将表达式转化为python代码，填写在下方单元格中的energy_experssion后。
请在下方单元格中填写，确认后请点击单元格左上角的➡️运行单元格。
代码
文本
代码
文本
[16]
# Constants
sqrt3 = np.sqrt(3)
a = 1.42 # bond length in Angstroms
t = 2.7 # hopping parameter in eV
# Dispersion relation of graphene in the tight binding model
def E(kx, ky):
energy_expression =
return energy_expression
代码
文本
填写并运行完毕后，你可以运行下方的代码单元格可视化紧束缚模型下的石墨烯色散关系。如果一切顺利，你将看到红色和蓝色的两条能带，分别对应能量为正值和负值，并且在K点处有唯一的交点。如果没有得到理想的结果，我们在下面给出了答案。
代码
文本
[33]
# PLOT 2D
import numpy as np
import matplotlib.pyplot as plt
# High symmetry points in reciprocal space
Gamma = np.array([0, 0])
M = np.array([2 * np.pi / (sqrt3 * a), 0])
K = np.array([2 * np.pi / (sqrt3 * a), 2 * np.pi / (3 * a)])
# Generate k-points along the high symmetry path Gamma -> K -> M -> Gamma
num_points = 100
k_points = np.concatenate([
np.linspace(Gamma, K, num_points, endpoint=False),
np.linspace(K, M, num_points, endpoint=False),
np.linspace(M, Gamma, num_points)
])
# Calculate energy bands
E_plus = []
E_minus = []
for k in k_points:
kx, ky = k
energy = E(kx, ky)
E_plus.append(energy)
E_minus.append(-energy)
# Plot energy bands
fig, ax = plt.subplots(figsize=(8, 6))
k_dist = np.linspace(0, 3 * num_points, 3 * num_points)
ax.plot(k_dist, E_plus, 'r-', label='$E_+(\mathbf{k})$')
ax.plot(k_dist, E_minus, 'b-', label='$E_-(\mathbf{k})$')
ax.axhline(0, color='green', linestyle='--') # Zero energy line
# Adding high symmetry points
ax.axvline(num_points, color='k', linestyle='--')
ax.axvline(2 * num_points, color='k', linestyle='--')
# Set x-axis limits and labels
ax.set_xlim([0, 3 * num_points])
ax.set_xticks([0, num_points, 2 * num_points, 3 * num_points])
ax.set_xticklabels(['$\\Gamma$', '$K$', '$M$', '$\\Gamma$'])
# Labels and title
ax.set_ylabel('Energy [$t$]')
ax.set_title('Dispersion Relation for Graphene in Tight-Binding Model')
ax.legend()
plt.tight_layout()
plt.show()
代码
文本
下一个单元格是答案
代码
文本
答案在此，点击展开
energy_expression = t * np.sqrt(1 + 4 * np.cos(ky * a / 2)**2 + 4 * np.cos(ky * a / 2) * np.cos(sqrt3 * kx * a / 2))
代码
文本
我们还可以画出三维的石墨烯的色散关系。你看到图中的狄拉克锥了吗？
代码
文本
- 狄拉克锥（Dirac cone）：描述狄拉克费米子在晶体中的电子能带结构的一种图形。这种结构最著名的例子是在石墨烯中发现的。电子的能量和动量具有线性的色散关系，因此其费米能级附近的电子能带结构呈现出上下两个锥体，分别代表电子和空穴。由于是锥型，电传导可以用无质量费米子的电荷载流子来描述。两个锥体的顶端刚好相连，交点即“狄拉克点”。在狄拉克点处，能带之间没有带隙。除了石墨烯，狄拉克锥结构也在其他二维材料中被发现，例如过渡金属二硫族化物（如MoS2）和某些拓扑绝缘体。理解和研究狄拉克锥有助于揭示材料中的新奇物理现象，例如高电子迁移率、自旋电子学、和量子霍尔效应等。这些研究不仅在基础物理学上有重要意义，还可能在电子器件、传感器和其他高科技应用中带来革命性的变化。
代码
文本
已隐藏单元格
代码
文本
2.4 使用ABACUS计算石墨烯能带
代码
文本
2.4.1 计算前的准备
代码
文本
在使用ABACUS进行计算时，所有操作都需要在指定的文件夹目录下进行。因此我们首先进入/root目录，并创建两个路径Graphene和Si，并且在两种材料的目录下创建band和dos子目录，分别用于进行石墨烯和硅的能带及态密度计算。
代码
文本
[1]
import os
# 设置工作目录
work_dir = "/root"
os.chdir(work_dir)
# 创建 Graphene 和 Si 目录
Graphene_dir = os.path.join(work_dir, "Graphene")
Si_dir = os.path.join(work_dir, "Si")
os.makedirs(Graphene_dir, exist_ok=True)
os.makedirs(Si_dir, exist_ok=True)
# 在 Graphene 目录下创建子目录
Graphene_band_dir = os.path.join(Graphene_dir, "band")
Graphene_strained_band_dir = os.path.join(Graphene_dir, "strained_band")
Graphene_bilayer_band_dir = os.path.join(Graphene_dir, "bilayer_band")
Graphene_bilayer_band_noelec_dir = os.path.join(Graphene_bilayer_band_dir, "noelec")
Graphene_bilayer_band_elec_dir = os.path.join(Graphene_bilayer_band_dir, "elec")
Graphene_dos_dir = os.path.join(Graphene_dir, "dos")
os.makedirs(Graphene_band_dir, exist_ok=True)
os.makedirs(Graphene_strained_band_dir, exist_ok=True)
os.makedirs(Graphene_bilayer_band_dir, exist_ok=True)
os.makedirs(Graphene_bilayer_band_elec_dir, exist_ok=True)
os.makedirs(Graphene_bilayer_band_noelec_dir, exist_ok=True)
os.makedirs(Graphene_dos_dir, exist_ok=True)
# 在 Si 目录下创建 band 和 dos 子目录
Si_band_dir = os.path.join(Si_dir, "band")
Si_dos_dir = os.path.join(Si_dir, "dos")
os.makedirs(Si_band_dir, exist_ok=True)
os.makedirs(Si_dos_dir, exist_ok=True)
print(os.getcwd())
print(f"files in {Graphene_dir}:", os.listdir(Graphene_dir))
print(f"files in {Graphene_bilayer_band_dir}:", os.listdir(Graphene_bilayer_band_dir))
print(f"files in {Si_dir}:", os.listdir(Si_dir))
/root files in /root/Graphene: ['dos', 'band', 'strained_band', 'bilayer_band'] files in /root/Graphene/bilayer_band: ['noelec', 'elec'] files in /root/Si: ['dos', 'band']
代码
文本
计算开始前的输入文件：要成功运行ABACUS，需要准备三个必要的输入文件，分别是INPUT、KPT、STRU。
- INPUT: 可类比VASP的INCAR文件。INPUT文件是ABACUS的核心控制文件，它指定了ABACUS需要完成的任务和所采用的方法等关键信息。类似VASP，ABACUS的INPUT文件中的指令也采用“tag=value”的形式，其中tag是ABACUS规定的关键词，value是用户输入的符合ABACUS规定的值。但是区别在于：1.ABACUS中不能使用等号来赋值，因为ABACUS中参数的读取实际上使用的是C++里的getline函数，所以不能有“=”号 2.如果同一个tag在INPUT中出现多次，VASP只会识别第一次出现的值，而ABAUCS只会识别最后一次出现的值。
- KPT: 可类比VASP的KPOINTS文件。KPT文件用于描述在倒易空间布里渊区中k点的分布。
- STRU: 可类比VASP的POSCAR以及POTCAR文件。STRU文件同时包含了计算对象的结构信息、赝势以及原子轨道信息。
接下来我们将在每个子目录下准备这三个输入文件并详细解释文件格式以及参数设置。
代码
文本
计算结束后的输出文件：计算结束后通常会在当前目录下产生time.json文件和OUT.xxx文件夹。time.json文件是ABACUS开发团队使用rapidjson做一些格式化的输出，主要用来监控软件的耗时和效率。OUT.xxx文件夹包含ABACUS产生的所有输出文件。其中xxx为后缀，默认名为OUT.ABACUS，可通过在INPUT文件中设置suffix参数修改。OUT.xxx文件夹通常包含：
- INPUT：与输入文件INPUT同名，但其内容更多，因为其中包含了运行前输入文件INPUT中没有设置的变量的默认值。
- STRU_READIN_ADJUST.cif：输入文件STRU的cif格式，可以用VESTA等可视化软件打开。
- istate.info：包含电子能级和占据数。
- kpoints：包含k点位置和权重。
- running_scf.log：可类比VASP的OUTCAR文件。最重要的输出文件，记录了ABACUS详尽的计算过程。
- warning.log：ABACUS执行中的stderr。
以及由INPUT中的各种out_xxx参数控制的输出文件，例如本教程将会使用到的out_chg、out_band、out_bandgap以及out_dos，分别控制是否输出电荷、能带、带隙以及态密度信息。
代码
文本
2.4.2 输入文件及关键参数
代码
文本
STRU文件：
按ATOMIC_SPECIES，NUMERICAL_ORBITAL等关键字分为不同的section。
- ATOMIC_SPECIES 下的每一行对应一种元素，每种元素需要：元素类型 相对原子质量，赝势文件名，赝势文件类型。
- NUMERICAL_ORBITAL下为每种的元素的轨道。关于轨道的命名和使用详见：NAO
- LATTICE_CONSTANT 下为晶胞的缩放因子。注意这里的单位为Bohr。 1.8897259886 Bohr = 1.0 Angstrom。
- LATTICE_VECTORS 下为晶胞的三个基矢，它们乘上面的缩放因子得到晶胞真正的大小。
- ATOMIC_POSITIONS 下为原子坐标。首先第一行为原子坐标的类型，有分数坐标“Direct”和直角坐标“Cartesian”。再下面开始以“元素类型 磁矩 原子数 原子位置”为一组，多种原子就有多个组。原子位置一行，前三个数为原子位置，后三个数决定原子在弛豫时能否移动。0：不能移动；1：可以移动。
代码
文本
由于石墨烯是二维材料，因此我们在Z方向上加了 25Å 的真空层。
代码
文本
[34]
#如果你在2.1节最后运行了格式转换的代码，那么这里你就可以选择跳过该单元格，从而使用你的结构完成2.4.2和2.4.3节的计算；如果没有，请运行该单元格生成后续计算所需的STRU文件
stru = f'''ATOMIC_SPECIES
C 12.011 C_ONCV_PBE-1.0.upf upf201
NUMERICAL_ORBITAL
C_gga_8au_100Ry_2s2p1d.orb
LATTICE_CONSTANT
1.8897259886
LATTICE_VECTORS
2.4677236000 0.0000000000 0.0000000000
1.2338618000 2.1371113000 0.0000000000
0.0000000000 0.0000000000 25.0000000000
ATOMIC_POSITIONS
Direct
C #label
0 #magnetism
2 #number of atoms
0.3333333300 0.3333333300 0.0000000000 m 1 1 1
0.6666666700 0.6666666700 0.0000000000 m 1 1 1
'''
with open(os.path.join(Graphene_band_dir,"STRU"),"w") as f: f.write(stru)
with open(os.path.join(Graphene_dos_dir,"STRU"),"w") as f: f.write(stru)
代码
文本
这样我们就准备好了石墨烯的结构文件STRU。
代码
文本
INPUT文件：
- 第一行为关键词，不能有任何修改，否则ABACUS无法识别该文件。并且任何在该关键词之前的内容都会被忽略。
- ABACUS中使用
#当做注释符。 pseudo_dir和orbital_dir为分别为原子的赝势和轨道所在的文件夹，对于平面波 (PW) 计算，只需要指定赝势；如果选用数值原子轨道作为基组 (LCAO) 进行计算，还需要额外指定轨道文件。大部分ABAUCS相关镜像都会在系统盘下自带/abacus-develop文件夹，即ABACUS源代码，在其中的/abacus-develop/tests/PP_ORB文件夹中就包含赝势和轨道文件。ecutwfc为平面波函数的截断能。首先，这里能量的单位为Ry；其次，ecutwfc的默认值为50；最后，即便使用数值原子轨道算法，也需要设置该值，因为ABACUS中局部赝势部分和相关力是由平面波基集计算的。详情见官方说明：ecutwfc。scf_thr决定自洽收敛的精度。basis_type决定使用的基组类型。ABACUS目前同时支持pw和lcao基组。calculation决定计算类型：不同于VASP使用ISTART，ICHARG，IBRION等参数控制计算类型，ABACUS直接使用关键词控制。scf-自洽计算，nscf-非自洽计算，relax-离子弛豫等。详情见：calculationsymmetry决定是否打开体系的对称性，1是打开对称性；0关闭晶体对称性，但保留时间反演对称性；-1关闭所有对称性。init_chg决定初始化电荷的方式。atomic：原子电荷密度叠加；file：从已有的电荷密度文件读取。- 输出文件相关参数
out_chg、out_band、out_bandgap以及out_dos：分别控制是否输出电荷、能带、带隙以及态密度信息。默认均为关闭。
当然，由于我们这里计算的体系非常简单，大部分参数可以不设置（使用默认值），所以不可能碰到所有的输出参数。全部的输入参数列表和说明详见官方文档：Full List of INPUT Keywords
代码
文本
🔖接下来，为了更深入地理解上述参数的作用，请你跟着提示做：
请完成下方单元格中scf计算和能带计算的输入文件中的calculation，inti_chg参数的设置。
请在下方单元格中填写，确认后请点击单元格左上角的➡️运行单元格。
代码
文本
[35]
input_scf = '''
INPUT_PARAMETERS
symmetry 1
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
init_chg
calculation
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_chg 1
'''
with open(os.path.join(Graphene_band_dir,"INPUT_scf"),"w") as f: f.write(input_scf)
with open(os.path.join(Graphene_dos_dir,"INPUT_scf"),"w") as f: f.write(input_scf)
input_nscf_band = '''
INPUT_PARAMETERS
symmetry 0
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
init_chg
calculation
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_band 1
'''
with open(os.path.join(Graphene_band_dir,"INPUT_nscf"),"w") as f: f.write(input_nscf_band)
input_nscf_dos = '''
INPUT_PARAMETERS
symmetry 0
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
init_chg file
calculation nscf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_dos 1
'''
with open(os.path.join(Graphene_dos_dir,"INPUT_nscf"),"w") as f: f.write(input_nscf_dos)
代码
文本
下一个单元格是答案
代码
文本
答案在此，点击展开
input_scf:
init_chg atomic
calculation scfinput_band:
init_chg file
calculation nscf
可以看到，相比于VASP，ABACUS的参数设置更为直观。当然，如果你本身使用的是VASP，习惯ABACUS的参数设置以及单位制等仍然需要一些时间。
代码
文本
这样我们就准备好了石墨烯的自洽、能带以及态密度计算的三个INPUT文件。
代码
文本
KPT文件：
常规K-mesh：
- 第一行为关键词。
- 第二行K点数量。手动设置K点时需最后修改该值为K点总数，若自动生成K点则设置该值为0.
- 第三行为产生K点的方式。对于K-mesh，有Gamma和Monkhorst Pack (MP) 两种。Gamma方式保证k点一定包含倒空间原点，而MP方式则不一定。
- 第四行为由空格分开的6个整数。前三个数字表示沿三个倒易基矢方向的k点数目，后三个数字表示k点网格的平移。
对于能带计算：
- 第一行也为关键词。
- 第二行为高对称K点数量。
- Line。
- 第四行开始为高对称点的坐标，以及每段高对称路径中间插入的K点数。
代码
文本
[98]
kpt_scf = '''
K_POINTS
0
Gamma
12 12 1 0 0 0
'''
with open(os.path.join(Graphene_band_dir,"KPT_scf"),"w") as f: f.write(kpt_scf)
with open(os.path.join(Graphene_dos_dir,"KPT_scf"),"w") as f: f.write(kpt_scf)
kpt_nscf_band = '''
K_POINTS
4
Line
0.000 0.000 0.000 50 # G
0.333 0.666 0.000 50 # K
0.500 0.500 0.000 50 # M
0.000 0.000 0.000 1 # G
'''
with open(os.path.join(Graphene_band_dir,"KPT_nscf"),"w") as f: f.write(kpt_nscf_band)
kpt_nscf_dos = '''
K_POINTS
0
Gamma
40 40 1 0 0 0
'''
with open(os.path.join(Graphene_dos_dir,"KPT_nscf"),"w") as f: f.write(kpt_nscf_dos)
代码
文本
这样我们就准备好了石墨烯的自洽、能带以及态密度计算的三个KPT文件。
注意：
- 由于在Z方向上加了真空层，所以在该方向上设置K点为1即可。
- 由于石墨烯的能带结构为狄拉克锥，因此在做态密度计算时必须设置非常密的K点，从而保证在狄拉克点附近有足够的采样精度。
代码
文本
2.4.3 能带计算
代码
文本
准备好了输入文件，我们就可以开始进行能带和态密度的计算。
代码
文本
[99]
# band calculation
os.chdir(Graphene_band_dir)
os.system("cp INPUT_scf INPUT && cp KPT_scf KPT")
os.system("abacus | tee run_scf.log")
os.system("cp INPUT_nscf INPUT && cp KPT_nscf KPT")
os.system("abacus | tee run_nscf.log")
已隐藏输出
代码
文本
我们使用ABACUS提供的python包abacus-plot绘制能带图。
代码
文本
[100]
# band plot
import os
import matplotlib.pyplot as plt
import shutil
from PIL import Image
def gen_band_config(efermi, energy_r,jsonf = "config.json"):
import json
config = {
"bandfile": "BANDS_1.dat",
"efermi": efermi,
"energy_range": energy_r,
"bandfig": "band.png",
"kptfile": "KPT",
"dpi": 300
}
json.dump(config,open(jsonf,"w"),indent=4)
def get_efermi(logfile):
with open(logfile) as f:
for line in f.readlines()[::-1]:
if " EFERMI = " in line:
return float(line.split()[-2])
# cp KPT into OUT
current_dir = os.getcwd()
target_dir = os.path.join(current_dir, "OUT.ABACUS")
source_path = os.path.join(current_dir, "KPT")
target_path = os.path.join(target_dir, "KPT")
shutil.copy(source_path, target_path)
# change to OUT
os.chdir("OUT.ABACUS")
# plot band
efermi = get_efermi("running_scf.log")
config = gen_band_config(efermi, [-4,4] ,"config.json")
os.system("abacus-plot -b")
# plot show
if os.path.isfile("band.png"):
plt.figure(dpi=300)
plt.axis('off')
plt.imshow(Image.open('band.png'))
plt.show()
else:
print("找不到band.png文件！")
代码
文本
可以看到，通过第一性原理计算，我们也成功得到了石墨烯的能带结构。
代码
文本
2.4.4 态密度计算
代码
文本
- 态密度：在凝聚态物理学中，态密度为某一能量附近每单位能量区间里微观状态的数目。具有同一能量的微观状态被称为简并的，简并态的个数叫做简并数。在离散能级处，简并数就是相应能量的态密度。在连续和准连续能态处，设为态密度，则处在能量和区间的态的个数为。
代码
文本
[8]
# dos calculation
os.chdir(Graphene_dos_dir)
os.system("cp INPUT_scf INPUT && cp KPT_scf KPT")
os.system("abacus | tee scf.log")
os.system("cp INPUT_nscf INPUT && cp KPT_nscf KPT")
os.system("abacus | tee nscf.log")
已隐藏输出
代码
文本
我们使用ABACUS提供的python包abacus-plot绘制态密度图。
代码
文本
[9]
# dos plot
import os
import matplotlib.pyplot as plt
import shutil
from PIL import Image
def gen_tdos_config(efermi, energy_r,dos_r,jsonf = "config.json"):
import json
config = {
"tdosfile": "TDOS",
"efermi": efermi,
"energy_range": energy_r,
"dos_range": dos_r,
"tdosfig": "tdos.png",
"dpi": 300
}
json.dump(config,open(jsonf,"w"),indent=4)
def get_efermi(logfile):
with open(logfile) as f:
for line in f.readlines()[::-1]:
if " EFERMI = " in line:
return float(line.split()[-2])
# change to OUT
os.chdir("OUT.ABACUS")
# plot dos
efermi = get_efermi("running_scf.log")
config = gen_tdos_config(efermi, [-4,4], [0,1], "config.json")
os.system("abacus-plot -d")
# plot show
if os.path.isfile("tdos.png"):
plt.figure(dpi=300)
plt.axis('off')
plt.imshow(Image.open('tdos.png'))
plt.show()
else:
print("找不到tdos.png文件！")
/usr/local/lib/python3.10/dist-packages/matplotlib/projections/__init__.py:63: UserWarning: Unable to import Axes3D. This may be due to multiple versions of Matplotlib being installed (e.g. as a system package and as a pip package). As a result, the 3D projection is not available.
warnings.warn("Unable to import Axes3D. This may be due to multiple versions of "
代码
文本
图中蓝色虚线代表费米能级。可以看到，石墨烯仅在费米能级处态密度为0，对应导带和价带的唯一交点的能量值，再次验证了石墨烯是零能隙材料。
代码
文本
拓展内容：
现在你已经掌握了ABACUS的基本输入参数及操作，那么让我们来做一点更有趣的事情吧！😆
代码
文本
石墨烯是碳的二维结晶形式，其非凡的电子迁移率和其他独特特性为纳米级电子学和光子学带来了巨大的前景。但有一个问题：石墨烯没有带隙！没有带隙极大地限制了石墨烯在电子领域的应用。例如，你可以用石墨烯构建场效应晶体管，但如果没有带隙，你就无法将它们关闭！然而，如果你能实现石墨烯带隙，你应该能够制造出非常好的晶体管。
代码
文本
🔖那么，请你思考一下，用什么方式可以打开石墨烯的带隙呢？
（提示：破坏石墨烯的对称性）
下一个单元格是答案
代码
文本
答案在此，点击展开
📖打开石墨烯的带隙
- 化学掺杂，包括表面吸附掺杂和替位式掺杂。前者不会破坏石墨烯的结构，通过表面电荷转移实现n、p型掺杂；后者直接使用硼、氮等原子替换石墨烯晶格中的碳原子，破坏石墨烯的空间反演对称性，从而打开带隙。
- 量子限制效应（零维：石墨烯量子点；一维：石墨烯纳米带），使得原本连续的能级离散化，从而引入带隙。
- 双层石墨烯电场调控（如图2所示）。无外电场时两层中的电子态是相同的，外加电场后破坏了垂直石墨烯表面方向上的反演对称性，从而打开带隙。
- 非对称应变（如图3所示），直接改变石墨烯的晶格结构，破坏石墨烯本身的空间反演对称性。
- 衬底诱导（如图4所示）。通过石墨烯与衬底的相互作用，有希望破坏其空间反演对称性从而打开能隙。
- 异质结构（如图4所示）。例如，通过石墨烯/锂/石墨烯异质结将相反手性的狄拉克费米子耦合起来，就可以实现手征对称性的破缺，从而打开带隙。
除此之外，还有更为复杂的转角石墨烯等。也可以将上述方式排列组合，例如双层-单层石墨烯外加电场、石墨烯纳米带-掺杂石墨烯堆叠等等。
代码
文本
接下来，我们将以“应变石墨烯”和“电场调控双层石墨烯”作为例子，使用ABACUS软件计算这两种情况下的石墨烯的能带结构。
代码
文本
图3 单层石墨烯(top)没有带隙。对称的双层石墨烯(middle)也缺少带隙。电场(箭头)将不对称引入双层结构(bottom)，产生可以选择性调节的带隙(Δ)。
（图源: Bilayer Graphene Gets a Bandgap)
（图源: Bilayer Graphene Gets a Bandgap)
代码
文本
图4 (a) 沿C-C键方向的非对称应变分布的石墨烯体系的能带结构、PDOS和ELF (b) Lx随Ly的变化情况 (c) 能隙宽度随Ly的变化情况
（doi: 10.1103/PhysRevB.78.075435 ）
（doi: 10.1103/PhysRevB.78.075435 ）
代码
文本
图5 (a)，(b)碳化硅上外延生长石墨烯的结构，其电子结构在K点打开带隙；(c)，(d)石墨烯/氮化硼异质结的结构，在K点原始狄拉克锥与二级狄拉克锥均打开能隙；(e)，(f)锂插层石墨烯的结构，其能带表现为石墨烯K点与K′点的狄拉克锥复制到Γ点并打开能隙
（来源: 二维材料的新奇物理及异质结的能带调控 ）
（来源: 二维材料的新奇物理及异质结的能带调控 ）
代码
文本
2.5 应变石墨烯能带计算
代码
文本
研究发现，具有对称应变分布的石墨烯始终是零带隙半导体，其赝带隙在弹性范围内随应变强度线性减小。然而，施加一定程度的非对称应变可以打开石墨烯的带隙。下面我们通过对石墨烯平行于C-C键方向的非对称应变来尝试打开单层石墨烯的带隙。
代码
文本
2.5.1 应变石墨烯的原子结构
代码
文本
我们已经准备好了建模的文件，直接运行即可。
注意⚠️：这里我们省略了两个步骤：
- 1.固定Ly，计算对应的平衡Lx。由于ABACUS还无法进行固定x/y/z轴的结构优化，这里我们通过计算一系列Lx值对应的总能，选择总能最低对应的Lx作为平衡Lx。计算脚本放在下方隐藏单元格中。最终得到平衡Lx值为1.218Å。
- 2.离子位置优化。在确定Lx后，还需要再做一次固定基矢的离子位置优化，从而得到最终的平衡构型。
代码
文本
[34]
import py3Dmol
from ase import Atoms
import numpy as np
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
# lattice parameter
Ly = 2.396
Lx = 1.218
# create cell
cell = [
[-Lx, Ly, 0],
[Lx, Ly, 0],
[0, 0, 25] # vaccum layer
]
# Position of two carbon atoms in graphene (fractional coordinates)
atom1_frac = [1/3, 1/3, 0.0]
atom2_frac = [2/3, 2/3, 0.0]
# frac 2 cart
atom1 = np.dot(atom1_frac, cell)
atom2 = np.dot(atom2_frac, cell)
# build
graphene = Atoms('C2', positions=[atom1, atom2], cell=cell, pbc=True)
# extend
extended_graphene = graphene * (5, 5, 1)
#show
view = py3Dmol.view(width=800, height=400)
xyz_data = extended_graphene.get_positions()
symbols = extended_graphene.get_chemical_symbols()
for i, (x, y, z) in enumerate(xyz_data):
view.addSphere({'center': {'x': x, 'y': y, 'z': z}, 'radius': 0.4, 'color': 'gray'})
view.zoomTo()
view.show()
代码
文本
请根据上面的数据，计算平行于C-C键方向的应变。
代码
文本
答案在此，点击展开
代码
文本
2.5.2 输入文件准备
代码
文本
[35]
stru = f'''
ATOMIC_SPECIES
C 12.011 C_ONCV_PBE-1.0.upf upf201
NUMERICAL_ORBITAL
C_gga_8au_100Ry_2s2p1d.orb
LATTICE_CONSTANT
1.8897259886
LATTICE_VECTORS
-1.2180000000 2.3960000000 0.0000000000 #latvec1
1.2180000000 2.3960000000 0.0000000000 #latvec2
0.0000000000 0.0000000000 25.0000000000 #latvec3
ATOMIC_POSITIONS
Direct
C #label
0 #magnetism
2 #number of atoms
0.3323413789 0.3323413789 0.0000000000 m 1 1 1
0.6676586211 0.6676586211 0.0000000000 m 1 1 1
'''
with open(os.path.join(Graphene_strained_band_dir,"STRU"),"w") as f: f.write(stru)
代码
文本
已隐藏单元格
代码
文本
[36]
input_scf = '''
INPUT_PARAMETERS
symmetry 0
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
calculation scf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_chg 1
'''
with open(os.path.join(Graphene_strained_band_dir,"INPUT_scf"),"w") as f: f.write(input_scf)
input_nscf = '''
INPUT_PARAMETERS
symmetry 0
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
init_chg file
calculation nscf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_band 1
out_bandgap 1
'''
with open(os.path.join(Graphene_strained_band_dir,"INPUT_nscf"),"w") as f: f.write(input_nscf)
代码
文本
[41]
kpt_scf = '''
K_POINTS
0
Gamma
12 12 1 0 0 0
'''
with open(os.path.join(Graphene_strained_band_dir,"KPT_scf"),"w") as f: f.write(kpt_scf)
kpt_nscf_band = '''
K_POINTS
4
Line
0.000 0.000 0.00 50 # G
-0.3146042305 0.3146042305 0.000 50 # K
0.000 0.500 0.000 50 # M
0.000 0.000 0.000 1 # G
'''
with open(os.path.join(Graphene_strained_band_dir,"KPT_nscf"),"w") as f: f.write(kpt_nscf_band)
代码
文本
2.5.3 应变石墨烯能带计算
代码
文本
[42]
# band calculation
os.chdir(Graphene_strained_band_dir)
os.system("cp INPUT_scf INPUT && cp KPT_scf KPT")
os.system("abacus | tee run_scf.log")
os.system("cp INPUT_nscf INPUT && cp KPT_nscf KPT")
os.system("abacus | tee run_nscf.log")
已隐藏输出
代码
文本
[43]
# band plot
import os
import matplotlib.pyplot as plt
import shutil
from PIL import Image
def gen_band_config(efermi, energy_r,jsonf = "config.json"):
import json
config = {
"bandfile": "BANDS_1.dat",
"efermi": efermi,
"energy_range": energy_r,
"bandfig": "band.png",
"kptfile": "KPT",
"dpi": 300
}
json.dump(config,open(jsonf,"w"),indent=4)
def get_efermi(logfile):
with open(logfile) as f:
for line in f.readlines()[::-1]:
if " EFERMI = " in line:
return float(line.split()[-2])
# cp KPT into OUT
current_dir = os.getcwd()
target_dir = os.path.join(current_dir, "OUT.ABACUS")
source_path = os.path.join(current_dir, "KPT")
target_path = os.path.join(target_dir, "KPT")
shutil.copy(source_path, target_path)
# change to OUT
os.chdir("OUT.ABACUS")
# plot band
efermi = get_efermi("running_scf.log")
config = gen_band_config(efermi, [-10,10] ,"config.json")
os.system("abacus-plot -b")
# plot show
if os.path.isfile("band.png"):
plt.figure(dpi=300)
plt.axis('off')
plt.imshow(Image.open('band.png'))
plt.show()
else:
print("找不到band.png文件！")
代码
文本
[44]
os.system("grep -i -n 'bandgap' running_nscf.log")
3171: E_bandgap 0.102018 eV
0
代码
文本
可以看到，施加沿C-C方向的非对称应变后成功地在石墨烯中产生了带隙。第一性原理计算和实验表明，施加的单轴应变可以在石墨烯中引入0.1eV左右的带隙（低于这个值时仅表现为狄拉克点位置在布里渊区的移动），这与我们的结果相符。并且相比于armchari方向（垂直于C-C键方向），沿zigzag方向（平行于C-C键方向）的应变更容易打开带隙。
代码
文本
2.6 电场调控双层石墨烯能带计算
代码
文本
2.6.1 双层石墨烯的原子结构
代码
文本
🔖接下来，为了更深刻地理解双层石墨烯的晶格结构，请你跟着提示做：
双层石墨烯的AB堆叠（又称伯纳尔堆叠形式，Bernal-stacked form）相对较为稳定，其中一半原子直接位于下层石墨烯六边形的中心上方，另一半原子直接位于六边形角的一个原子上方。请填写两层原子的相对坐标在下方单元格中。你还是可以通过'basis'选择不同的基组，然后填写对应的第一层原子的相对坐标，在此基础上再填写第二层原子的即可。
请在下方单元格中填写，确认后请点击单元格左上角的➡️运行单元格。
代码
文本
[3]
import numpy as np
# Fill in the cooridnate of four C atoms within the unit cell
# example
# atom1=[0.5, 0.4, 0.0]
# lattice parameter
a = 2.46
# lattice vector
cell = [
[a*(-1/2), a*np.sqrt(3)/2, 0],
[a*(1/2), a*np.sqrt(3)/2, 0],
[0, 0, 25] # vaccum layer
]
# interlayer distance
b = 3.35
b_frac = b/25
# choose a basis, basis = 'a'/'b'
basis = ''
# position of C atoms in the first layer, in fractional coordinates
atom1_frac= [ , ]
atom2_frac= [ , ]
# position of C atoms in the second layer, in fractional coordinates
atom3_frac= [ , ]
atom4_frac= [ , ]
代码
文本
填充完上一格中两层原子原子坐标后，请运行下面的单元格可视化石墨烯的结构，确保原子坐标的正确性。如果一切顺利，你将看到一个5x5的六方蜂窝型双层石墨烯结构，注意检查两层原子的相对位置是否正确。如果没有得到理想的结果，我们在下面给出了答案。
- 按住鼠标左键拖拽，可以旋转观察二维石墨烯结构
- 按住鼠标右键上下移动，可以进行视角缩放
代码
文本
[4]
import py3Dmol
from ase import Atoms
import numpy as np
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
# create cell
if basis == 'a':
cell = [
[a*np.sqrt(3)/2, a*(1/2), 0],
[a*np.sqrt(3)/2, a*(-1/2), 0],
[0, 0, 25] # vacuum layer
]
elif basis == 'b':
cell = [
[-a*np.sqrt(3)/2, a*(1/2), 0],
[a*np.sqrt(3)/2, a*(1/2), 0],
[0, 0, 25] # vacuum layer
]
else:
print("Invalid basis.")
# position of C atoms in the first layer, in fractional coordinates
atom1_frac_3 = [None, None, 0]
atom1_frac_3 = atom1_frac[:2] + atom1_frac_3[2:]
atom2_frac_3 = [None, None, 0]
atom2_frac_3 = atom2_frac[:2] + atom2_frac_3[2:]
# position of C atoms in the second layer, in fractional coordinates
atom3_frac_3 = [None, None, b_frac]
atom3_frac_3 = atom3_frac[:2] + atom3_frac_3[2:]
atom4_frac_3 = [None, None, b_frac]
atom4_frac_3 = atom4_frac[:2] + atom4_frac_3[2:]
# frac 2 cart
atom1 = np.dot(atom1_frac_3, cell)
atom2 = np.dot(atom2_frac_3, cell)
atom3 = np.dot(atom3_frac_3, cell)
atom4 = np.dot(atom4_frac_3, cell)
# build
graphene = Atoms('C4', positions=[atom1, atom2, atom3, atom4], cell=cell, pbc=True)
# extend
extended_graphene = graphene * (5, 5, 1)
# show
view = py3Dmol.view(width=800, height=400)
xyz_data = extended_graphene.get_positions()
symbols = extended_graphene.get_chemical_symbols()
for i, (x, y, z) in enumerate(xyz_data):
view.addSphere({'center': {'x': x, 'y': y, 'z': z}, 'radius': 0.4, 'color': 'gray'})
view.zoomTo()
view.show()
代码
文本
下一个单元格是答案
代码
文本
答案在此，点击展开
- basis = 'a'
atom1_frac = [1/3, 1/3]
atom2_frac = [2/3, 2/3]
atom3_frac = [ 1/3,1/3 ]
atom4_frac = [ 0,0 ]
- basis = 'b'
atom1_frac = [2/3, 1/3]
atom2_frac = [1/3, 2/3]
atom3_frac = [ 2/3,1/3 ]
atom4_frac = [ 0,0 ]
代码
文本
2.6.2 输入文件准备
代码
文本
[5]
stru = f'''ATOMIC_SPECIES
C 12.011 C_ONCV_PBE-1.0.upf upf201
NUMERICAL_ORBITAL
C_gga_8au_100Ry_2s2p1d.orb
LATTICE_CONSTANT
1.8897259886
LATTICE_VECTORS
2.4677236000 0.0000000000 0.0000000000
1.2338618000 2.1371113000 0.0000000000
0.0000000000 0.0000000000 25.0000000000
ATOMIC_POSITIONS
Direct
C #label
0 #magnetism
4 #number of atoms
0.3333333300 0.3333333300 0.0000000000 m 1 1 1
0.6666666700 0.6666666700 0.0000000000 m 1 1 1
0.3333333300 0.3333333300 0.1340000000 m 1 1 1
0.0000000000 0.0000000000 0.1340000000 m 1 1 1
'''
with open(os.path.join(Graphene_bilayer_band_noelec_dir,"STRU"),"w") as f: f.write(stru)
with open(os.path.join(Graphene_bilayer_band_elec_dir,"STRU"),"w") as f: f.write(stru)
代码
文本
[6]
input_scf_noelec = '''
INPUT_PARAMETERS
symmetry 1
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
calculation scf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_chg 1
'''
with open(os.path.join(Graphene_bilayer_band_noelec_dir,"INPUT_scf"),"w") as f: f.write(input_scf_noelec)
input_scf_elec = '''
INPUT_PARAMETERS
symmetry 0
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
calculation scf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
efield_flag 1
dip_cor_flag 1
efield_amp 0.0019447 # 1 eV/nm
out_chg 1
'''
with open(os.path.join(Graphene_bilayer_band_elec_dir,"INPUT_scf"),"w") as f: f.write(input_scf_elec)
input_nscf = '''
INPUT_PARAMETERS
symmetry 0
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
init_chg file
calculation nscf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_band 1
out_bandgap 1
'''
with open(os.path.join(Graphene_bilayer_band_noelec_dir,"INPUT_nscf"),"w") as f: f.write(input_nscf)
with open(os.path.join(Graphene_bilayer_band_elec_dir,"INPUT_nscf"),"w") as f: f.write(input_nscf)
代码
文本
INPUT文件外加电场相关参数介绍：
efield_flag：是否打开电场。dip_cor_flag：是否添加偶极矫正。Default: 2。efield_dir：外加电场方向。0：x；1：y；2：z。efield_pos_max：最大电势位置。Default: center of vacuum - width of vacuum / 20。efield_pos_dec：电势降低区域。Default: vacuum / 10。efield_amp：电场的振幅。在efield_pos_max+efield_pos_dec-1到efield_pos_max的区域内，电势随斜率efield_amp增加，然后在efield_pos_max到efield_pos_max+efield_pos_dec的区域内减小，方向沿efield_dir。单位：a.u.。0.0019447 a.u. = 1 V/nm
代码
文本
[7]
kpt_scf = '''
K_POINTS
0
Gamma
12 12 1 0 0 0
'''
with open(os.path.join(Graphene_bilayer_band_noelec_dir,"KPT_scf"),"w") as f: f.write(kpt_scf)
with open(os.path.join(Graphene_bilayer_band_elec_dir,"KPT_scf"),"w") as f: f.write(kpt_scf)
kpt_nscf_band = '''
K_POINTS
4
Line
0.000 0.000 0.000 50 # G
0.333 0.666 0.000 50 # K
0.500 0.500 0.000 50 # M
0.000 0.000 0.000 1 # G
'''
with open(os.path.join(Graphene_bilayer_band_noelec_dir,"KPT_nscf"),"w") as f: f.write(kpt_nscf_band)
with open(os.path.join(Graphene_bilayer_band_elec_dir,"KPT_nscf"),"w") as f: f.write(kpt_nscf_band)
代码
文本
2.6.3 无外电场时双层石墨烯的能带结构
代码
文本
[8]
# band calculation
os.chdir(Graphene_bilayer_band_noelec_dir)
os.system("cp INPUT_scf INPUT && cp KPT_scf KPT")
os.system("abacus | tee run_scf_noelec.log")
os.system("cp INPUT_nscf INPUT && cp KPT_nscf KPT")
os.system("abacus | tee run_nscf_noelec.log")
已隐藏输出
代码
文本
[23]
# band plot
# cp KPT into OUT
current_dir = os.getcwd()
target_dir = os.path.join(current_dir, "OUT.ABACUS")
source_path = os.path.join(current_dir, "KPT")
target_path = os.path.join(target_dir, "KPT")
shutil.copy(source_path, target_path)
# change to OUT
os.chdir("OUT.ABACUS")
# plot band
efermi = get_efermi("running_scf.log")
config = gen_band_config(efermi, [-4,4] ,"config.json")
os.system("abacus-plot -b")
# plot show
if os.path.isfile("band.png"):
plt.figure(dpi=300)
plt.axis('off')
plt.imshow(Image.open('band.png'))
plt.show()
else:
print("找不到band.png文件！")
代码
文本
可以看到，无外电场时，双层石墨烯的带隙是闭合的。
代码
文本
2.6.4 有外电场时双层石墨烯的能带结构
代码
文本
[24]
# band calculation
os.chdir(Graphene_bilayer_band_elec_dir)
os.system("cp INPUT_scf INPUT && cp KPT_scf KPT")
os.system("abacus | tee run_scf_elec.log")
os.system("cp INPUT_nscf INPUT && cp KPT_nscf KPT")
os.system("abacus | tee run_nscf_elec.log")
已隐藏输出
代码
文本
[25]
# band plot
# cp KPT into OUT
current_dir = os.getcwd()
target_dir = os.path.join(current_dir, "OUT.ABACUS")
source_path = os.path.join(current_dir, "KPT")
target_path = os.path.join(target_dir, "KPT")
shutil.copy(source_path, target_path)
# change to OUT
os.chdir("OUT.ABACUS")
# plot band
efermi = get_efermi("running_scf.log")
config = gen_band_config(efermi, [-4,4] ,"config.json")
os.system("abacus-plot -b")
# plot show
if os.path.isfile("band.png"):
plt.figure(dpi=300)
plt.axis('off')
plt.imshow(Image.open('band.png'))
plt.show()
else:
print("找不到band.png文件！")
代码
文本
[26]
os.system("grep -i -n 'bandgap' running_nscf.log")
3779: E_bandgap 0.181449 eV
0
代码
文本
理论和计算的表明，双层石墨烯中打开带隙的阈值为1eV/nm，带隙值约0.1eV，继续增大电场可以增加带隙值至0.25~0.30eV。可以看到，我们的计算结果表明，1eV/nm的电场已经可以在双层石墨烯中产生带隙。接下来，让我们再增加电场强度，看看能带会有什么样的变化。
代码
文本
[27]
# band calculation
os.chdir(Graphene_bilayer_band_elec_dir)
os.system("sed -i '/efield_amp/s/.*/efield_amp 0.0077784/' INPUT_scf") # 4 eV/nm
os.system("cp INPUT_scf INPUT && cp KPT_scf KPT")
os.system("abacus | tee run_scf_elec.log")
os.system("cp INPUT_nscf INPUT && cp KPT_nscf KPT")
os.system("abacus | tee run_nscf_elec.log")
已隐藏输出
代码
文本
[28]
# band plot
# cp KPT into OUT
current_dir = os.getcwd()
target_dir = os.path.join(current_dir, "OUT.ABACUS")
source_path = os.path.join(current_dir, "KPT")
target_path = os.path.join(target_dir, "KPT")
shutil.copy(source_path, target_path)
# change to OUT
os.chdir("OUT.ABACUS")
# plot band
efermi = get_efermi("running_scf.log")
config = gen_band_config(efermi, [-4,4] ,"config.json")
os.system("abacus-plot -b")
# plot show
if os.path.isfile("band.png"):
plt.figure(dpi=300)
plt.axis('off')
plt.imshow(Image.open('band.png'))
plt.show()
else:
print("找不到band.png文件！")
代码
文本
[29]
os.system("grep -i -n 'bandgap' running_nscf.log")
3779: E_bandgap 0.333814 eV
0
代码
文本
可以看到，此时带隙进一步增大，为0.3eV，与理论和计算结果相符。并且在K点附近出现了“墨西哥帽”形状的能带，导带底和价带顶均发生了偏移，不再位于K点处。该能带结构也可以通过紧束缚模型求解，感兴趣的同学可以自行尝试。
代码
文本
3. Si的能带和态密度计算
代码
文本
硅是地壳中含量最丰富的元素之一，它在自然界中以二氧化硅（SiO2）的形式存在，同时也以硅石、石英、石英砂等矿物形式存在。硅同时是一种半导体材料，其在集成电路（IC）、太阳能电池、光电器件等领域有广泛的应用，其优良的电子性质使得硅成为现代电子工业的基石。
图6 集成电路中无处不在的硅
（来源: Britannica ）
（来源: Britannica ）
代码
文本
与石墨烯不同，硅晶体的结构为三维金刚石结构，这是一种具有高对称性的立方晶格结构，在该结构中，每个硅原子都与其他四个硅原子共享电子，形成共价键。并且，硅存在（间接）带隙。硅晶体的能带结构是其电子性质的重要决定因素。通过第一性原理计算，我们可以确定硅的禁带宽度以及导带和价带的形状和对称性，从而更深刻地理解硅电子性质。
图7 硅的能带结构及布里渊区
（来源: tuwien ）
（来源: tuwien ）
代码
文本
3.1. 硅的原子结构
代码
文本
图8 硅的晶体结构。左侧为参考坐标系，图中绿色为基向量，原胞中处于不同位点的两个原子A、B用黄色字体标记。
代码
文本
通过前面的内容，你应该已经掌握了如何通过设置基矢和原子的相对坐标来构建二维材料的结构。现在就让我们来尝试构建三维金刚石结构硅的结构。
代码
文本
🔖硅晶体建模：
金刚石结构可以看作是由两套面心立方晶格套构而成，图8中的A、B原子就分别对应两套晶格的顶点处的原子。硅原胞就包含A、B原子，原胞基矢在图中给出。请在下方单元格中填写三个单位基向量的坐标以及A，B原子的相对位置，完成Si晶体的建模。
请在下方单元格中填写，确认后请点击单元格左上角的➡️运行单元格。
代码
文本
[1]
import os
import py3Dmol
from ase import Atoms
import numpy as np
from ase.visualize.plot import plot_atoms
import matplotlib.pyplot as plt
# lattice parameter
a = 5.43
# Fill in the cell and the cooridnate of two Si atoms within the unit cell
cell = [
[ , , ], #a1
[ , , ], #a2
[ , , ] #a3
]
# position of first atom, in fractional coordinates
atomA_frac = [ , , ]
# position of second atom, in fractional coordinates
atomB_frac = [ , , ]
代码
文本
填充完上一格中的晶胞基矢以及原子坐标后，请运行下面的单元格可视化硅的结构，确保原子坐标的正确性。
晶体的原胞是晶格的最小周期性单元，是能带计算的首选。但由于不能反映晶格的对称性，因此在可视化结构时不如单胞那样便于分析。你可以尝试寻找其中存在的正四面体结构，如果可以找到，那么你的建模大概率是成功的。如果没有，或者你不确定自己的结果，我们在下方准备了答案。
代码
文本
[4]
# frac 2 cart
atomA = np.dot(atomA_frac, cell)
atomB = np.dot(atomB_frac, cell)
# build
silicon = Atoms('Si2', positions=[atomA, atomB], cell=cell, pbc=True)
# extend
extended_silicon = silicon * (5, 5, 5)
# show
view = py3Dmol.view(width=800, height=400)
xyz_data = extended_silicon.get_positions()
symbols = extended_silicon.get_chemical_symbols()
for i, (x, y, z) in enumerate(xyz_data):
view.addSphere({'center': {'x': x, 'y': y, 'z': z}, 'radius': 0.1, 'color': 'blue'})
view.zoomTo()
view.show()
代码
文本
下一个单元格是答案
代码
文本
答案在此，点击展开
cell = [
[ 0.5, 0.5, 0 ], #a1
[ 0, 0.5, 0.5 ], #a2
[ 0.5, 0, 0.5 ] #a3
]
atomA_frac = [ 0, 0, 0 ]
atomB_frac = [ 0.25, 0.25, 0.25 ]
代码
文本
3.2. 计算前的准备
代码
文本
我们首先准备Si原胞的结构文件。你可以尝试将你刚才的建模和下面STRU格式文件中LATTICE_VECTORS以及ATOMIC_POSITIONS部分对应。
代码
文本
[10]
stru = f'''ATOMIC_SPECIES
Si 28 Si_ONCV_PBE-1.0.upf upf201
NUMERICAL_ORBITAL
Si_gga_8au_60Ry_2s2p1d.orb
LATTICE_CONSTANT
10.2
LATTICE_VECTORS
0.5000000000 0.5000000000 0.0000000000
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.0000000000 0.5000000000
ATOMIC_POSITIONS
Direct
Si #label
0 #magnetism
2 #number of atoms
0.0000000000 0.0000000000 0.0000000000 m 1 1 1
0.2500000000 0.2500000000 0.2500000000 m 1 1 1
'''
with open(os.path.join(Si_band_dir,"STRU"),"w") as f: f.write(stru)
with open(os.path.join(Si_dos_dir,"STRU"),"w") as f: f.write(stru)
代码
文本
然后准备输入文件INPUT。自洽和态密度计算的INPUT和石墨烯的完全一致；能带计算的INPUT中我们设置了输出带隙，其余参数和石墨烯一致。
代码
文本
[11]
input_scf = '''
INPUT_PARAMETERS
symmetry 1
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
calculation scf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_chg 1
'''
with open(os.path.join(Si_band_dir,"INPUT_scf"),"w") as f: f.write(input_scf)
with open(os.path.join(Si_dos_dir,"INPUT_scf"),"w") as f: f.write(input_scf)
input_nscf_band = '''
INPUT_PARAMETERS
symmetry 0
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
init_chg file
calculation nscf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_band 1
out_bandgap 1
'''
with open(os.path.join(Si_band_dir,"INPUT_nscf"),"w") as f: f.write(input_nscf_band)
input_nscf_dos = '''
INPUT_PARAMETERS
symmetry 0
pseudo_dir /abacus-develop/tests/PP_ORB
orbital_dir /abacus-develop/tests/PP_ORB
init_chg file
calculation nscf
ecutwfc 50
scf_thr 1e-7
scf_nmax 300
basis_type lcao
out_dos 1
'''
with open(os.path.join(Si_dos_dir,"INPUT_nscf"),"w") as f: f.write(input_nscf_dos)
代码
文本
最后准备输入文件KPT。由于石墨烯是二维蜂窝状材料，而硅是三维金刚石结构材料，对应的布里渊区维度以及高对称点坐标不同，因此非自洽计算时KPT文件也不相同。
代码
文本
[12]
kpt_scf = '''
K_POINTS
0
Gamma
9 9 9 0 0 0
'''
with open(os.path.join(Si_band_dir,"KPT_scf"),"w") as f: f.write(kpt_scf)
with open(os.path.join(Si_dos_dir,"KPT_scf"),"w") as f: f.write(kpt_scf)
kpt_nscf_band = '''
K_POINTS
8
Line
0.000 0.000 0.000 30 # G
0.500 0.000 0.500 20 # X
0.625 0.250 0.625 20 # U
0.375 0.375 0.750 50 # K
0.000 0.000 0.000 50 # G
0.500 0.500 0.500 20 # L
0.500 0.250 0.750 20 # W
0.500 0.000 0.500 1 # X
'''
with open(os.path.join(Si_band_dir,"KPT_nscf"),"w") as f: f.write(kpt_nscf_band)
kpt_nscf_dos = '''
K_POINTS
0
Gamma
15 15 15 0 0 0
'''
with open(os.path.join(Si_dos_dir,"KPT_nscf"),"w") as f: f.write(kpt_nscf_dos)
代码
文本
3.2 能带计算
代码
文本
准备好了输入文件，我们就可以开始进行能带和态密度的计算。
代码
文本
[13]
# band calculation
os.chdir(Si_band_dir)
os.system("cp INPUT_scf INPUT && cp KPT_scf KPT")
os.system("abacus | tee scf.log")
os.system("cp INPUT_nscf INPUT && cp KPT_nscf KPT")
os.system("abacus | tee nscf.log")
已隐藏输出
代码
文本
[14]
# 读取输出文件中的带隙信息
os.system('grep -i -n "bandgap" OUT.ABACUS/running_nscf.log')
4309: E_bandgap 0.823445 eV
0
代码
文本
[15]
# band plot
import os
import matplotlib.pyplot as plt
import shutil
from PIL import Image
def gen_band_config(efermi, energy_r,jsonf = "config.json"):
import json
config = {
"bandfile": "BANDS_1.dat",
"efermi": efermi,
"energy_range": energy_r,
"bandfig": "band.png",
"kptfile": "KPT",
"dpi": 300
}
json.dump(config,open(jsonf,"w"),indent=4)
def get_efermi(logfile):
with open(logfile) as f:
for line in f.readlines()[::-1]:
if " EFERMI = " in line:
return float(line.split()[-2])
# cp KPT into OUT
current_dir = os.getcwd()
target_dir = os.path.join(current_dir, "OUT.ABACUS")
source_path = os.path.join(current_dir, "KPT")
target_path = os.path.join(target_dir, "KPT")
shutil.copy(source_path, target_path)
# change to OUT
os.chdir("OUT.ABACUS")
# plot band
efermi = get_efermi("running_scf.log")
config = gen_band_config(efermi, [-4,4] ,"config.json")
os.system("abacus-plot -b")
# plot show
if os.path.isfile("band.png"):
plt.figure(dpi=300)
plt.axis('off')
plt.imshow(Image.open('band.png'))
plt.show()
else:
print("找不到band.png文件！")
/usr/local/lib/python3.10/dist-packages/matplotlib/projections/__init__.py:63: UserWarning: Unable to import Axes3D. This may be due to multiple versions of Matplotlib being installed (e.g. as a system package and as a pip package). As a result, the 3D projection is not available.
warnings.warn("Unable to import Axes3D. This may be due to multiple versions of "
代码
文本
🔖Note：
- 从能带图中我们可以看到，不同于石墨烯，硅的能带结构中存在一个能量范围，在该能量范围内不存在任何电子能级，这便是“禁带”。禁带的上限为“导带底”，下限为“价带顶”。从能带图中可以看出，硅的价带顶位于倒空间的Gamma点，而价带顶位于Gamma点和X点之间。这说明：硅是一种间接带隙半导体。
- 上面我们读取到的带隙为0.82eV，比实际的1.17eV小，造成这种差异的主要原因是：上面计算中默认的交换关联泛函为PBE泛函，其对相互作用较强的轨道电子的描述不够准确，在计算半导体带隙时往往存在“带隙低估”的问题。为了获得更精确的带隙结果，我们可以使用杂化泛函，但计算的成本会显著提高（约为原来的100倍）。
如果你想了解杂化泛函的结果并学习如何使用ABACUS完成杂化泛函计算，可以继续浏览notebook：ABACUS计算模拟实例 | VIII. 基于HSE06的态密度与能带计算
。
代码
文本
3.3 态密度计算
代码
文本
[16]
# dos calculation
os.chdir(Si_dos_dir)
os.system("cp INPUT_scf INPUT && cp KPT_scf KPT")
os.system("abacus | tee scf.log")
os.system("cp INPUT_nscf INPUT && cp KPT_nscf KPT")
os.system("abacus | tee nscf.log")
已隐藏输出
代码
文本
[19]
# dos plot
import os
import matplotlib.pyplot as plt
import shutil
from PIL import Image
def gen_tdos_config(efermi, energy_r,dos_r,jsonf = "config.json"):
import json
config = {
"tdosfile": "TDOS",
"efermi": efermi,
"energy_range": energy_r,
"dos_range": dos_r,
"tdosfig": "tdos.png",
"dpi": 300
}
json.dump(config,open(jsonf,"w"),indent=4)
def get_efermi(logfile):
with open(logfile) as f:
for line in f.readlines()[::-1]:
if " EFERMI = " in line:
return float(line.split()[-2])
# change to OUT
os.chdir("OUT.ABACUS")
# plot dos
efermi = get_efermi("running_scf.log")
config = gen_tdos_config(efermi, [-4,4], [0,5], "config.json")
os.system("abacus-plot -d")
# plot show
if os.path.isfile("tdos.png"):
plt.figure(dpi=300)
plt.axis('off')
plt.imshow(Image.open('tdos.png'))
plt.show()
else:
print("找不到tdos.png文件！")
/usr/local/lib/python3.10/dist-packages/matplotlib/projections/__init__.py:63: UserWarning: Unable to import Axes3D. This may be due to multiple versions of Matplotlib being installed (e.g. as a system package and as a pip package). As a result, the 3D projection is not available.
warnings.warn("Unable to import Axes3D. This may be due to multiple versions of "
代码
文本
从硅的态密度图中可以更直观地看到，不同于石墨烯仅在费米能级处态密度为0，硅在费米能级附近的一段范围内的态密度均为0，这段能量范围不存在电子可以跃迁的能级，即禁带。
代码
文本
恭喜你完成了2024秋计算材料学-上机练习的第一课！🎉
代码
文本
🔚结语
通过这一课，你应该了解了石墨烯以及硅的晶体结构，倒空间、布里渊区及高对称点等概念，紧束缚模型下石墨烯的色散关系，掌握了如何使用第一性原理软件ABACUS完成自洽、能带及态密度计算。如果你对ABACUS这款软件感兴趣，并且还想尝试更多的计算实例，可以点击下方了解ABACUS计算模拟合集链接： ABACUS计算模拟实例 | 概述 。
代码
文本
已赞3
本文被以下合集收录
ABACUS@计算材料学 | 计算模拟实例
weiqingzhou@whu.edu.cn
更新于 2024-08-09
13 篇8 人关注
ABACUS计算模拟实例
FermiNoodles
更新于 2024-07-15
13 篇1 人关注
推荐阅读
公开
2024秋计算材料学-上机练习：ABACUS能带和态密度计算副本lujilin@dp.tech
更新于 2024-08-07
公开
2024秋计算材料学-上机练习：经验势与机器学习势stanfordbshan
更新于 2024-09-13