ABACUS 有限位移声子实操

下载完整输入、42个力计算结果、6个热输运条件和绘图脚本。解压进入附件根目录。先用已有结果熟悉提取与画图：

bash analyze.sh

结果表在 reference/analysis，6张高清图在 reference/figures。这个命令只处理附件内的实际计算结果，不启动新的DFT或热输运任务。下面按准备输入、提交、检查、提取的顺序操作。新计算请使用新解压目录，已有 work 时脚本会停止，避免覆盖结果。

1 先准备平衡结构，再生成位移#

本例从晶格常数5.43 Å的金刚石硅原胞出发，原胞2个原子，扩成2×2×2超胞后是16个原子。晶格与原子位置在本教程中固定，用于演示有限位移流程；自己的材料先完成合理的结构优化，再检查应力、残余力和计算参数。

使用 Phono3py 4.8.2 生成0.03 Å位移，同时为后续三阶力常数准备结构。终端运行：

module load phono3py/4.8.2
python -B make_displacements.py

会生成 generated-inputs，并逐项核对它与附件 inputs 的39套输入一致。38套是实际位移结构，另1套是未位移超胞。内部位移数据集有114个位置，其中76个组合被距离筛选排除，不是已经计算且力为零的结构。本脚本专门复现这个Si16算例，换材料需要重新指定元素、赝势、轨道、超胞和收敛参数。

2 看懂实际INPUT、STRU和KPT#

每个目录都包含三份输入。先看第一套：

cat inputs/si16-equilibrium/INPUT
cat inputs/si16-equilibrium/STRU
cat inputs/si16-equilibrium/KPT

实际 INPUT：

INPUT_PARAMETERS
suffix tutorial
pseudo_dir ../../PP_ORB
orbital_dir ../../PP_ORB
ntype 1
calculation scf
basis_type lcao
ks_solver elpa
device gpu
precision double
nbands 48
ecutwfc 80
dft_functional pbe
nspin 1
smearing_method fixed
symmetry -1
scf_thr 1e-8
scf_nmax 100
mixing_type pulay
mixing_beta 0.4
cal_force 1
cal_stress 0
out_chg 0

calculation scf 加 cal_force 1 表示固定该位移结构，只计算原子力。不要把位移结构重新优化，否则会抹去构建力常数所需要的扰动。这里使用LCAO、PBE、80 Ry、48条能带、双精度和电子收敛1E-8；有限差分对力噪声敏感，自己的体系应检查力随参数的收敛。

ATOMIC_SPECIES
Si	28.0855	Si_ONCV_PBE-1.0.upf

NUMERICAL_ORBITAL
Si_gga_8au_100Ry_2s2p1d.orb

LATTICE_CONSTANT
1.0

LATTICE_VECTORS
0.00000000000000 10.2612169463150 10.2612169463150
10.2612169463150 0.00000000000000 10.2612169463150
10.2612169463150 10.2612169463150 0.00000000000000

ATOMIC_POSITIONS
Direct

Si
0
16
0.00000000000000 0.00000000000000 0.00000000000000 1 1 1
0.50000000000000 0.00000000000000 0.00000000000000 1 1 1
0.00000000000000 0.50000000000000 0.00000000000000 1 1 1
0.50000000000000 0.50000000000000 0.00000000000000 1 1 1
0.00000000000000 0.00000000000000 0.50000000000000 1 1 1
0.50000000000000 0.00000000000000 0.50000000000000 1 1 1
0.00000000000000 0.50000000000000 0.50000000000000 1 1 1
0.50000000000000 0.50000000000000 0.50000000000000 1 1 1
0.12500000000000 0.12500000000000 0.12500000000000 1 1 1
0.62500000000000 0.12500000000000 0.12500000000000 1 1 1
0.12500000000000 0.62500000000000 0.12500000000000 1 1 1
0.62500000000000 0.62500000000000 0.12500000000000 1 1 1
0.12500000000000 0.12500000000000 0.62500000000000 1 1 1
0.62500000000000 0.12500000000000 0.62500000000000 1 1 1
0.12500000000000 0.62500000000000 0.62500000000000 1 1 1
0.62500000000000 0.62500000000000 0.62500000000000 1 1 1

质量28.0855为实际Si有效质量；赝势为 Si_ONCV_PBE-1.0.upf，DZP数值轨道为 Si_gga_8au_100Ry_2s2p1d.orb。本例208个NAO、64个价电子。LATTICE_CONSTANT 与晶格向量共同决定真实长度，附件在工具接口中显式转换Bohr与Å；不要直接把Å晶胞交给要求Bohr的ABACUS格式写出接口。

实际KPT：

K_POINTS
0
Gamma
2 2 2 0 0 0

这是16原子超胞的电子k网格。声子能带路径、8×8×8声子DOS网格和后面的热输运q网格是不同的采样，不能相互替代。

3 提交力计算#

新解压目录直接提交：

sbatch submit.slurm

脚本在分配到的单张P100上运行39套主要输入，再运行3套仅用于比较的Gamma电子网格输入，最后自动提取力和生成声子结果。不同位移结构可以分配给不同已申请GPU并行运行，每套保持独立输出目录；附件的单GPU顺序模板便于直接使用。

环境为ABACUS3.10.1-p100，GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2、CUDA12.8、ELPA GPU，1 MPI进程、1线程、GPU双精度。不要在登录节点运行DFT。

如希望手动逐个执行，在已经分配的GPU环境中使用下列方式；它与上面的整套提交方式二选一：

module load abacus/3.10.1-p100
python3 -B setup.py
python3 -B run_case.py si16-displacement-001

其余名称见 CASES.json。仅运行一个位移还不能提取整套力常数，必须完成所有39套。Gamma诊断的3个输入见 GAMMA-CASES.json，单套加 --gamma；整套模板已经包含它们。

4 提取力，检查声子#

完整GPU脚本会自动执行后处理。若使用手动方式完成所有输入，再运行：

bash after-forces.sh

步骤依次为：读取原始力→显式去除平均力漂移及平衡力→构建二阶/受限三阶力常数→导出交换格式→计算声子能带和DOS→对照电子k采样→准备热输运输入。主要输出是：

采用Phono3py4.8.2环境内的Phonopy4.8.1、symfc1.7.3；hiPhive1.5负责力常数导出。脚本按步骤切换模块，不混用不同Python依赖环境。

实际Γ点光学频率约16.273437 THz；16原子超胞对应原胞6条声子支，DOS积分为每Si2原胞6个模式。8×8×8声子网格和路径上最小频率约−1.2×10⁻⁷ THz，属于Γ声学支数值舍入，原始有符号频率保留。DOS的0.1 THz高斯宽度只是作图/积分设置，不能当作声子寿命。

5 有虚频先查什么#

本例只把超胞电子采样改为Gamma时，最低频率约−1.991860 THz；加密成2×2×2后，明显虚频消失。几何、赝势、基组和截断能都保持一致。由此可见，负频率不总是结构本身不稳定，也可能是电子采样导致力常数不可靠。

遇到明显虚频，依次检查平衡结构、电子收敛、k网格、基组/截断能、位移幅度及超胞尺寸；别把负频率取绝对值后继续算热导率。本例Gamma-only的热输运结果不作为教程的材料结果。加密后没有明显虚频，只说明当前采样通过了这一步检查，不代表所有计算参数都已经收敛。

接口流程可对照Phonopy官方ABACUS说明；在线新版CLI与集群安装版本可能不同，本例下载包使用已实跑的Python API。