FourPhonon：极性TDEP与Ewald输入实操

下载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具有非零Born有效电荷。TDEP的极性拟合分开处理短程力常数与偶极长程贡献。FourPhonon读取这份短程力常数时，也必须使用匹配的长程处理；不能关掉非解析项后仍把结果当成完整的ZnS动力学矩阵。

本例原胞顺序为S、Zn，质量分别32.065、65.38 u；Born电荷为各向同性−2.06281946、+2.06281946，介电张量对角元6.13221705。电荷中性和正定介电张量同时满足。Ewald耦合参数0.786429876361969直接从实际TDEP输出读取，不是随意调出来的拟合数。

2. 完整CONTROL#

cat inputs/Ewald-on/CONTROL

&allocations
nelements=2,natoms=2,ngrid(:)=3 3 3
&end
&crystal
lfactor=.1
lattvec(:,1)=0.0000000000000000e+00,2.7250000000000001e+00,2.7250000000000001e+00
lattvec(:,2)=2.7250000000000001e+00,0.0000000000000000e+00,2.7250000000000001e+00
lattvec(:,3)=2.7250000000000001e+00,2.7250000000000001e+00,0.0000000000000000e+00
elements="S","Zn",types=1 2
masses(1)=32.065,masses(2)=65.38
gfactors=0 0
scell(:)=3 3 3
positions(:,1)=2.5000000000000000e-01,2.5000000000000000e-01,2.5000000000000000e-01
positions(:,2)=0.0000000000000000e+00,0.0000000000000000e+00,-0.0000000000000000e+00
epsilon(:,1)=6.1322170500000048e+00,-2.1729834716729258e-32,-2.1729834716729250e-32
epsilon(:,2)=-1.3799192150979052e-32,6.1322170500000048e+00,1.8648346805709670e-32
epsilon(:,3)=1.3799192150979054e-32,1.3799192150979054e-32,6.1322170500000048e+00
born(:,1,1)=-2.0628194599999996e+00,-9.8725033629063653e-17,-2.6501858756268098e-32
born(:,2,1)=-9.8725033629063628e-17,-2.0628194599999996e+00,9.8725033629063616e-17
born(:,3,1)=2.4659469448956835e-32,9.8725033629063653e-17,-2.0628194599999996e+00
born(:,1,2)=2.0628194599999996e+00,9.8725033629063653e-17,2.6501858756268098e-32
born(:,2,2)=9.8725033629063628e-17,2.0628194599999996e+00,-9.8725033629063616e-17
born(:,3,2)=-2.4659469448956835e-32,-9.8725033629063653e-17,2.0628194599999996e+00
&end
&parameters
T=300,scalebroad=1
Ewald=7.8642987636196904e-01
&end
&flags
nonanalytic=.true.,espresso=.false.,tdep=.true.
isotopes=.false.,autoisotopes=.false.,nanowires=.false.
onlyharmonic=.true.,four_phonon=.false.,convergence=.false.
&end

lfactor=.1配合上面以Å书写的晶格，按当前FourPhonon约定转换长度。tdep=.true.选择TDEP文件读入；nonanalytic=.true.恢复极性长程部分。质量、介电张量、Born电荷、原子位置与原胞顺序必须全部匹配。Ewald取本次TDEP实际输出中的值；更换拟合后重新读取。

sed -n '/Coupling parameter in Ewald/p' reference-results/fit/outfile.forceconstant
cmp reference-results/fit/outfile.forceconstant reference-results/Ewald-on/infile.forceconstant

cmp无输出表示逐字节相同。不要将已经含完整长程贡献的普通实空间IFC直接当作本例的短程IFC再加一次Ewald。

3. 提交计算并读取#

sbatch run.slurm
# 作业结束后
bash read.sh work

提交脚本先拟合、再运行TDEP参考和两组FourPhonon计算。在分配资源中，FourPhonon的核心命令是：

module load fourphonon/1.3-fix20260924-8173m
cd work/Ewald-on
mpirun -np 1 fourphonon > native.log
head BTE.omega
head BTE.qpoints_full

BTE.omega为不可约网格频率，单位rad/ps；BTE.qpoints_full给出完整网格及映射。不能直接把完整网格的27个q点序号用来索引不可约频率行。

4. 与关闭长程项的结果对照#

Ewald-off-bare与主计算使用同一份短程IFC，只将nonanalytic关闭。它用于展示漏掉长程项的后果，是另一种近似模型，不能替代完整极性结果。

左图对照同一套原胞、质量和长程参数的TDEP原生结果；26个非Γ网格点、每点6条分支的最大差为0.00698633 rad/ps。右图保留只用短程IFC的结果，显示遗漏项确实会改变频率。

本例onlyharmonic=.true.、four_phonon=.false.，只演示二阶谐频率和极性读入；没有用这四个构型计算四声子散射或热导率。完整热导计算仍需要对应的三/四阶IFC、q网格与展宽检查。下一篇见单位与Γ点。具体读取选项参考FourPhonon官方手册。