abinit学习日记二十七——tbs_2.abi

前一个教程中，我们创建了 BS 方程需要的两个 WFK 文件和 SCR 文件，本节旨在演示如何在 Tamm-Dancoff 近似（TDA）中使用 Haydock 迭代技术进行标准的激子计算。

在运行计算之前，我们需要将这次计算与之前在 tbs_1.abi 中生成的结果联系起来。使用 Unix 命令，

ln -s 444_shifted_WFK tbs_2i_WFK
ln -s 444_SCR tbs_2i_SCR

创建两个符号链接，分别指向偏移的 WFK 文件和 SCR 文件。这样做的原因将在我们讨论输入文件时变得清晰。

串行启动计算程序会花费几分钟的时间，因此我使用并行运行命令

mpirun -n 8 abinit tbs_2.abi >& log_2

花费了几十秒的时间。接着看看输入文件

# Crystalline silicon
# BS run: Tamm-Dancoff approximation solved with the Haydock algorithm.# BSE run with Haydock iterative method (only resonant + W + v)
optdriver 99  # BS calculation
# 设置99来调用 BS 计算的子程序
irdwfk  1     # Read the WFK file produced in tbs_1
# 读取 tbs_1 中生成的 WFK 文件
irdscr  1     # Read the SCR file produced in tbs_1
# 读取 tbs_1 中生成的 SCR 文件bs_calctype       1       # L0 is constructed with KS orbitals and energies.
# 使用 KS 轨道和能量构建 L0
mbpt_sciss        0.8 eV  # Scissors operator used to correct the KS band structure.
# 用于修正 KS 能带结构的剪刀算符
bs_exchange_term  1       # Exchange term included.
# 包含交换项
bs_coulomb_term   11      # Coulomb term included using the full matrix W_GG'
# 使用完整的矩阵 W_GG' 包含库仑项
bs_coupling       0       # Tamm-Dancoff approximation.
# Tamm-Dancoff 近似bs_loband         2
nband             8
# 跃迁空间的基集包括了所有带隙附近的带，其索引介于 2 和 8 之间bs_freq_mesh 0 6 0.02 eV   # Frequency mesh.
# 频率网格bs_algorithm        2      # Haydock method.
# Haydock 方法（这是默认值）
bs_haydock_niter   200     # Max number of iterations for the Haydock method.
# Haydock 方法的最大迭代次数
bs_haydock_tol     0.05 0  # Tolerance for the iterative method.
# 迭代方法的容忍度
zcut               0.15 eV # complex shift to avoid divergences in the continued fraction.
# 复杂位移，用于避免连分数中的发散# Definition of the k-point grid
# MUST be equal to the grid used for generating the WFK file.
kptopt 1                  # Option for the automatic generation of k points,
# 自动生成 k 点的选项
ngkpt  4 4 4              # This mesh is too coarse for optical properties.
# 这个网格对光学性质来说太粗糙了
nshiftk 1
shiftk    0.11 0.21 0.31  # This shift breaks the symmetry of the k-mesh.
# 这个偏移破坏了 k 网格的对称性
chksymbreak 0             # Mandatory for using symmetry-breaking k-meshes in the BS code.
# 对于使用对称性破坏的 k 网格的 BS 代码来说是必需的ecutwfn 12.0              # Cutoff for the wavefunction.
# 波函数的截断
ecuteps 2.0               # Cutoff for W and /bare v used to calculate the BS matrix elements.
# 用于计算 BS 矩阵元的 W 和 /bare v 的截断
inclvkb 2                 # The commutator for the optical limit is correctly evaluated.
# 正确评估光学极限的交换子
gw_icutcoul 3               # Legacy value# VARIABLES COMMON TO THE DIFFERENT DATASETS
# 其他共同参数# Definition of the unit cell: fcc
acell 3*10.217         # This is equivalent to   10.217 10.217 10.217
rprim  0.0  0.5  0.5   # FCC primitive vectors (to be scaled by acell)0.5  0.0  0.50.5  0.5  0.0# Definition of the atom types
ntypat 1         # There is only one type of atom
znucl 14         # The keyword "zatnum" refers to the atomic number of the# possible type(s) of atom. The pseudopotential(s)# mentioned in the "files" file must correspond# to the type(s) of atom. Here, the only type is Silicon.# Definition of the atoms
natom 2           # There are two atoms
typat 1 1         # They both are of type 1, that is, Silicon.
xred              # Reduced coordinate of atoms0.0  0.0  0.00.25 0.25 0.25# Definition of the planewave basis set
ecut 12        # Maximal kinetic energy cut-off, in Hartreeistwfk *1
nstep  50      # Maximal number of SCF cycles
diemac 12.0pp_dirpath "$ABI_PSPDIR"pseudos "Psdj_nc_sr_04_pw_std_psp8/Si.psp8"

bs_calctype = 1 指定使用 Kohn-Sham 轨道和能量构建独立粒子极化率。为了模拟自能修正，KS 能量通过一个剪刀算符进行了修正，能量为 mbpt_sciss = 0.8 eV。这使得我们可以避免对过渡空间中包含的每个状态进行繁琐的 GW 计算。剪刀算符对硅是一个合理的近似，但在更复杂的系统中可能会失效，因为 GW 修正不能用简单的刚性平移来表示。

bs_exchange_term = 1 告诉代码计算核的交换部分，因此这个计算包括了局域场效应。
bs_coulomb_term 用于选择库仑项的不同选项（请花时间阅读该变量的描述以及 Bethe-Salpeter 笔记中的相关方程）。bs_coupling = 0 指定忽略非对角耦合块（Tamm-Dancoff 近似）。这种特定的参数组合对应于包含局域场效应的 Tamm-Dancoff 近似中的 Bethe-Salpeter 计算。

bs_freq_mesh 0 6 0.02 eV 这个三元组定义了一个覆盖 [0, 6] eV 范围的线性网格，步长为 0.02 eV。需要注意的是，网格中频率点的数量对计算时间没有显著影响，但是跃迁空间中包含的带数，结合 k 点的数量，定义了可以描述的频率范围。因此，bs_loband 和 nband 应该进行严格的收敛性研究。

bs_algorithm 指定了用于计算宏观介电函数的算法。这里我们使用了 Haydock 迭代技术，其最大迭代次数由 bs_haydock_niter 给出。迭代算法在连续两次光谱计算结果的差异小于 bs_haydock_tol 时停止。输入变量 zcut 给出了避免连分数发散的复位移。从物理角度看，这个参数模拟了吸收峰的实验展宽。在这个测试中，由于 k 网格的粗糙，我们需要使用一个比默认值（0.1 eV）稍大的值，以促进 Haydock 算法的收敛。理想情况下，应该通过减少 zcut 的值，同时增加 k 点的数量，来进行收敛性研究。

kptopt、ngkpt、nshiftk 和 shiftk 的值必须与用于指定 WFK 文件网格的值相等。chksymbreak = 0 用于绕过对称性破坏的检查，否则代码会停止运行。

其他细节可以参见文章

接下来我们看看输出文件。输出文件 tbs_2.abo 报告了计算的基本参数以及如果迭代方法不收敛时发出的任何警告。请花一些时间理解其结构。在开始阅读之前给出了三个问题：

• 跃迁空间中包含了多少跃迁？
• 用于评估光学极限的方向有多少？
• 在连分数中使用的洛伦兹展宽的值是多少？

最重要的结果存储在五个不同的文件中：

• tbs_2o_BSR：该二进制文件存储了共振块的上三角部分（矩阵是厄米的，因此只计算并保存非冗余部分）。BSR 文件可用于使用变量 getbsreso 或 irdbsreso 从之前的计算中重新启动运行。这种重新启动能力对于在未收敛时重新启动 Haydock 方法或以不同的 zcut 值执行 Haydock 计算非常有用。getbsreso 和 irdbsreso 也很方便，如果你想在预先存在的 TDA 计算基础上包含耦合，因为代码使用两个不同的文件来存储共振和耦合块（BSC 是存储耦合项的文件的前缀）。
• tbs_2o_HAYDR_SAVE：这是包含 Haydock 方法结果的二进制文件：三对角矩阵的系数和迭代算法中使用的三个向量。通常用于在未收敛时重新启动算法（参见相关的输入变量 gethaydock 和 irdhaydock）。
• tbs_2o_RPA_NLF_MDF 和 tbs_2o_GW_NLF_MDF：分别是使用 Kohn-Sham (KS) 能量和模拟的 GW 能量得到的不含局域场效应的随机相位近似（RPA）光谱（缩写：NLF 表示 No Local Field，MDF 表示 Macroscopic Dielectric Function）。
• tbs_2o_EXC_MDF：格式化的文件，报告包含激子效应的宏观介电函数。

EXC_MDF 文件包含了我们计算的最重要的结果，因此值得花一些时间来讨论它的格式。

# Macroscopic dielectric function obtained with the BS equation.
#  RPA L0 with KS energies and KS wavefunctions     LOCAL FIELD EFFECTS INCLUDED
# RESONANT-ONLY calculation
# Coulomb term constructed with full W(G1,G2)
# Scissor operator energy =  0.8000 [eV]
# Tolerance =  0.0500 0.0000
# npweps  = 27
# npwwfn  = 531
# nbands  = 8
# loband  = 2
# nkibz   = 64
# nkbz    = 64
# Lorentzian broadening =  0.1500 [eV]
# 头部信息：报告了计算的基本参数
# List of q-points for the optical limit:
# q =  0.938821, 0.000000, 0.000000, [Reduced coords]
# q =  0.000000, 0.938821, 0.000000, [Reduced coords]
# q =  0.000000, 0.000000, 0.938821, [Reduced coords]
# q =  0.000000, 0.813043, 0.813043, [Reduced coords]
# q =  0.813043, 0.000000, 0.813043, [Reduced coords]
# q =  0.813043, 0.813043, 0.000000, [Reduced coords]
# q 点列表：列出了用于光学极限的 q 点方向
# omega [eV]    RE(eps(q=1)) IM(eps(q=1) RE(eps(q=2) ) ...
# 频率网格：报告了不同方向上的频率网格以及对应的宏观介电函数的实部和虚部0.000   17.9912    0.0000   17.9578    0.0000   14.2584    0.0000   13.9627    0.0000   17.0848    0.0000   17.0422    0.00000.020   17.9917    0.0072   17.9583    0.0072   14.2587    0.0051   13.9630    0.0048   17.0853    0.0068   17.0426    0.00670.040   17.9931    0.0145   17.9597    0.0145   14.2597    0.0101   13.9640    0.0097   17.0866    0.0136   17.0439    0.01340.060   17.9955    0.0217   17.9621    0.0217   14.2614    0.0152   13.9656    0.0145   17.0889    0.0203   17.0461    0.02010.080   17.9989    0.0290   17.9655    0.0289   14.2638    0.0203   13.9678    0.0194   17.0920    0.0271   17.0492    0.02670.100   18.0032    0.0362   17.9698    0.0362   14.2668    0.0254   13.9707    0.0243   17.0961    0.0339   17.0532    0.03350.120   18.0085    0.0435   17.9751    0.0434   14.2705    0.0305   13.9742    0.0291   17.1010    0.0407   17.0581    0.04020.140   18.0147    0.0508   17.9813    0.0507   14.2749    0.0356   13.9784    0.0340   17.1069    0.0476   17.0639    0.0469
...................

现在回答上述问题：

• 包含的跃迁数量：跃迁空间中包含的跃迁数量由 bs_loband 和 nband 决定，即从第2带到第8带的所有跃迁。
• 评估光学极限的方向数量：默认情况下，代码考虑了 q 空间的六个不同方向（三个倒格子基向量和三个笛卡尔轴）。可以通过输入变量 gw_nqlwl 和 gw_qlwl 指定自定义方向。
• 洛伦兹展宽的值：在连分数中使用的洛伦兹展宽值由 zcut 指定，为 0.15 eV。

可以使用以下命令实现文件 tbs_2o_MDF.nc 结果可视化

(abipy_env) (abinit_env)[hanjinzhuo@mu012 Work_bs$] abiopen.py tbs_2o_MDF.nc --expose --seaborn
================================= File Info =================================
Name: tbs_2o_MDF.nc
Directory: /public/home/hanjinzhuo/local/src/abinit-10.4.7/tests/tutorial/Input/Work_bs
Size: 177.07 kB
Access Time: Tue Nov  4 10:12:54 2025
Modification Time: Mon Nov  3 15:50:32 2025
Change Time: Mon Nov  3 15:50:32 2025================================= Structure =================================
Full Formula (Si2)
Reduced Formula: Si
abc   :   3.823046   3.823046   3.823046
angles:  60.000000  60.000000  60.000000
pbc   :       True       True       True
Sites (2)#  SP       a     b     c
---  ----  ----  ----  ----0  Si    0     0     01  Si    0.25  0.25  0.25Abinit Spacegroup: spgid: 0, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True
================================== Q-points ==================================
0) [+0.939, +0.000, +0.000]
1) [+0.000, +0.939, +0.000]
2) [+0.000, +0.000, +0.939]
3) [+0.000, +0.813, +0.813]
4) [+0.813, +0.000, +0.813]
5) [+0.813, +0.813, +0.000]Loading all matplotlib figures before showing them. It may take some time...
All figures in memory, elapsed time: 5.824 s

这两幅图绘制的是宏观介电函数（Macroscopic Dielectric Function, DF）随频率变化的曲线。宏观介电函数的实部和虚部分别提供了材料对电场的极化和能量损耗的信息。实部与材料的折射率相关，而虚部与材料的吸收系数相关。通过分析这些曲线，可以了解材料在不同频率下的光学和电子特性。

收敛激子光谱需要仔细分析许多不同的参数：

• bs_loband
• nband
• ecutwfn
• ecuteps
• ngkpt
• nshiftk
• shiftk

由于内存需求与布里渊区中 k 点的数量、价带数量和跃迁空间中包含的导带数量的乘积成二次方关系，因此找到精度与计算效率之间的良好折衷非常重要。

首先，应选择感兴趣的频率范围，因为这个选择对需要包含在跃迁空间中的价带和导带数量有重要影响。光学光谱在带的数量上的收敛速度通常比 GW 修正更快，因为只有那些能量“接近”所研究频率范围的跃迁才会贡献。

ecutwfn 通常起次要作用，因为它只影响振子矩阵元的精度。建议避免对初始基集进行任何截断，方法是将 ecutwfn 设置为略大于生成 WFK 文件时使用的 ecut 值。只有在遇到内存问题时，才应该考虑截断初始平面波基集。不过，这种问题通常可以通过增加处理器数量或适当选择 gwmem 来解决。

ecuteps 的值影响库仑项矩阵元的精度，这是产生激子的基本项。因此，ecuteps 应当进行严格的收敛性测试。一般来说，可以将 ecuteps 选择为与收敛 GW 修正所需的值相等或有时甚至更小。

正如已经提到的：光学光谱对布里渊区采样的收敛性很差，因此 k 点数量的收敛性测试是收敛性研究中最重要且繁琐的部分。因此，应在找到其他参数的收敛值后，再进行这个测试。