FourPhonon：从ZnS实算力到极性TDEP力常数
准备54原子ZnS位移结构、读取四次实际DFT力，使用原生TDEP进行极性力常数拟合。
阅读本篇目录
下载ZnS极性TDEP与FourPhonon示例包，保存为example.zip：
python3 -m zipfile -e example.zip .
cd polar-tdep-fourphonon-examples
bash results.sh
bash read.sh
bash plot.sh
以上命令读取本例的实算结果并画图，不会重算VASP。需要自己重新拟合时，运行sbatch run.slurm，然后用bash read.sh work读取新的原生输出。CPU队列按整节点申请56核，本例小程序实际使用1个MPI进程和1个线程；申请核数不能用来缩减整节点计费。绘图使用ase-3.29.0的Python，HDF5读取使用phonopy-4.8.3的Python，脚本已写好解释器路径。
本例使用TDEP25.03（d38f435）与FourPhonon1.3-fix20260924（5b814b3）的CPU程序，Intel oneAPI2024.2、Intel MPI2021.13、MKL顺序库。四个位移结构的力由VASP6.6.1实算获得：NVHPC25.3、CUDA12.8、HPC-X2.22.1、MKL2023.2，P100/V10016G单卡。
1. 从什么结构开始#
使用闪锌矿ZnS，常规晶格常数5.45 Å，原胞顺序固定为S、Zn。3×3×3原胞超胞包含54个原子。对独立位移方向取±0.02 Å，共四个位移结构。本例是静态有限位移的谐力常数演示，不是有限温度分子动力学拟合。
| 输入 | 内容与作用 |
|---|---|
inputs/fit/infile.ucposcar |
2原子平衡原胞，确定原子顺序 |
infile.ssposcar |
54原子平衡超胞 |
infile.positions |
四个实际计算构型的全精度分数坐标 |
infile.forces |
按构型、按原子排列的实际DFT力，eV/Å |
infile.meta |
原子数54、构型数4；后两项0 |
infile.lotosplitting |
同一组Born有效电荷与介电张量 |
infile.isotopes |
与原胞顺序一致的32.065、65.38 u |
infile.stat |
实际DFT能量；未用的温度/时间/压力等列为0 |
不能将infile.stat中的零温度、零时间列理解成一条热轨迹。换体系时，应一起更新结构、力、物种顺序、质量、Born电荷及介电张量。
2. 四个VASP输入如何准备#
实际POSCAR较长，完整文件在optional-DFT-inputs/disp-000至disp-003。可以直接检查：
head -n 12 optional-DFT-inputs/disp-000/POSCAR
cat optional-DFT-inputs/disp-000/INCAR
cat optional-DFT-inputs/disp-000/KPOINTS
实际INCAR：
SYSTEM = ZnS finite displacement
ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
NELM = 100
ALGO = Normal
ISMEAR = 0
SIGMA = 0.05
IBRION = -1
NSW = 0
ISYM = 0
LREAL = .FALSE.
LASPH = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
实际KPOINTS：
ZnS54 matched primitive6cube density
0
Gamma
2 2 2
0 0 0
这里2×2×2的超胞网格保持与原2×2×2超胞、3×3×3网格相同的原胞k点密度。POTCAR使用S、Zn顺序的PBE势，需要自行准备有授权的文件；资料包不分发POTCAR。每个结构单独计算，不能把前一个构型的原子排序带进下一个构型。提取vasprun.xml中的力时保留构型和原子顺序。XML坐标只打印8位小数，本例用全精度实际POSCAR构造TDEP位置文件；不要据打印舍入重新修改实际结构。
3. 用原生TDEP拟合#
以下命令应写入计算脚本，在申请的CPU资源中执行；完整顺序已放进run.slurm：
source /etc/profile.d/lmod.sh
module load tdep/25.03-8173m
cd inputs/fit
extract_forceconstants --polar -rc2 5.5 > fit.log
--polar使用极性长程分离，-rc2 5.5是请求的二阶相互作用截断，单位Å。实际可用截断仍受超胞大小限制；本例实际4.51890 Å、每个原子29个邻居，不应把请求5.5 Å直接写成实际截断。
sed -n '1,12p' outfile.forceconstant
sed -n '/FORCES (eV\/A)/,$p' fit.log
得到outfile.forceconstant后，把它作为FourPhonon的infile.forceconstant。这个文件尾部同时包含极性信息，不能只复制前面的二阶矩阵。
4. 实际结果说明什么#
原16原子超胞的实际截断只到2.35992 Å，力残差RMS为0.009423 eV/Å；改为54原子、同类独立位移后，较长截断拟合的RMS为0.001116 eV/Å。本例采用0.005 eV/Å作流程比较条件，两个超胞的结果都保留。
图中结构显示平衡超胞，实际四个位移POSCAR随包提供。该比较只说明本例有限位移数据与截断的匹配情况；正式有限温度TDEP需要代表目标温度的足够独立构型，并继续检查超胞、采样、截断和力误差。
可下载拟合比较CSV与四次DFT条件CSV。下一步见极性Ewald输入。TDEP输入文件含义可查官方文件说明。