ABACUS：DFPT声子与有限位移对照
以两原子Si为例，准备PW输入、提交Γ点DFPT、读取声子，并与同条件有限位移受力比较。
阅读本篇目录
想用线性响应计算声子，先明确两件事：电子基态要收敛，DFPT响应本身也要收敛。本例用两原子Si原胞在Γ点完成这条流程，并用同一赝势、晶胞和电子网格下的正负位移受力作独立参考。
unzip abacus-dfpt-si-wave684-*.zip
cd abacus-dfpt-si
bash results.sh
bash plot.sh
以上只读取参考数据并重绘，不启动新计算。重新计算请使用新解压的目录。 本例使用ABACUS实验功能个人环境，固定3.11.0-beta10源码。公共3.10.1模块仍可用于原有计算；不能直接用旧模块运行本篇新版功能。
bash use-toolchain.sh ../abacus-ri-env/toolchain
bash prepare.sh
sbatch job.slurm
squeue -u "$USER"
脚本申请cpu-8173m的56核独占整节点，本例内部使用一个MPI进程和一个线程，仍按整节点计费。CPU资源暂不可用时可以先读取参考数据。这里的小网格用于学习输入输出，实际研究需分别收敛k点、截断能、基组或空带数。
1. 准备结构与k点#
晶格常数5.43 Å。LATTICE_CONSTANT采用每Å的Bohr数，因此下列晶格矢量数值可以按Å理解：
ATOMIC_SPECIES
Si 28.0855 Si_ONCV_PBE-1.0.upf
LATTICE_CONSTANT
1.889726877774355
LATTICE_VECTORS
0 2.715 2.715
2.715 0 2.715
2.715 2.715 0
ATOMIC_POSITIONS
Direct
Si
0
2
0.000000000000000 0.000000000000000 0.000000000000000 1 1 1
0.25 0.25 0.25 1 1 1
KPT：
K_POINTS
0
Gamma
1 1 1 0 0 0
这里只使用Γ点电子采样；它适合快速掌握流程，不能代表声子已经对电子k网格收敛。
2. 写DFPT输入#
INPUT_PARAMETERS
suffix dfpt
calculation scf
esolver_type dfpt
basis_type pw
ntype 1
nspin 1
nelec 8
nbands 8
ecutwfc 60
scf_thr 1e-6
scf_nmax 160
dft_functional pbe
symmetry 0
ks_solver cg
smearing_method fixed
mixing_type broyden
mixing_beta 0.3
cal_force 1
cal_stress 0
out_chg 0
dfpt_qmesh 1 1 1
dfpt_compute_q0 0
dfpt_loto 0
dfpt_conv_thr 1e-6
dfpt_max_iter 160
dfpt_mix_beta 0.2
esolver_type dfpt选择本版实验性PW响应求解器。scf_thr 1e-6是电子基态阈值，dfpt_conv_thr 1e-6是响应阈值，两者各自控制不同迭代。dfpt_qmesh 1 1 1只求Γ点，不能据此画整条声子色散。- 本篇
dfpt_compute_q0 0、dfpt_loto 0；极性晶体的介电、Born电荷与LO–TO修正见后续SiC教程。 symmetry 0保留本例的显式响应流程。不要为了消除虚频而任意改动结果。
3. 提交与读取原生结果#
sbatch job.slurm
# 将下面的作业号换成返回的数字
tail -f job-作业号.log
作业脚本依次计算DFPT与两个有限位移SCF。结束后：
grep -n 'SCF IS CONVERGED' work/dfpt/OUT.dfpt/running_scf.log
grep -n -E 'DFPT phonon frequencies|mode.*cm' work/dfpt/run.log
python3 extract.py
cat my-results.csv
声子频率写在标准输出run.log中，不能只在running_scf.log里找。本例六个原始模式为：
| 模式 | 原始频率（cm⁻¹） |
|---|---|
| 1–3，声学 | −15.220736、−13.980168、−13.319948 |
| 4–6，光学 | 743.508997、743.517606、743.520895 |
原始声学模有小负值，通常以虚频显示。本例保留这些值；先收敛电子k网格、基组与响应参数，再判断实际材料稳定性，不把小体系教学值直接解释成结构失稳或完全稳定。
4. 用有限位移受力交叉查看#
两个输入目录force-plus/与force-minus/将第一个Si原子沿笛卡尔x方向移动±0.005 Å，其他电子条件相同：
diff -u inputs/force-minus/STRU inputs/force-plus/STRU
grep -A5 'TOTAL-FORCE' work/force-plus/OUT.dfpt/running_scf.log
grep -A5 'TOTAL-FORCE' work/force-minus/OUT.dfpt/running_scf.log
cat reference-data/Si-finite-force-column.csv
力常数列按Hᵢₓ,₀ₓ = −[Fᵢₓ(+h)−Fᵢₓ(−h)]/(2h)计算。本例再利用相同Si子晶格与立方对称性，得到一个光学频率参考742.261986 cm⁻¹。它是一个受力方向的对照，不是完整的所有原子、所有方向有限位移矩阵。
DFPT光学频率相对参考的最大差约0.17%。两种方法相互接近，只说明当前同条件小例相符，仍需针对实际材料收敛网格。
5. 最后会得到什么#
work/dfpt/run.log：原生Γ点频率。work/force-*/OUT.dfpt/running_scf.log：正负位移受力。my-results.csv：本次运行提取值。reference-data/：本集群已完成算例的频率和受力数据，便于比较读取方式。
要计算完整声子色散，应进一步设计q点或超胞有限位移流程，不能把本例的六个Γ点频率沿路径重复绘制。