FourPhonon：极性ZnS的非解析修正

极性材料的长程静电响应需要与结构匹配的Born有效电荷和高频介电张量，不能只把Si例子的nonanalytic改成true。本篇使用本集群已有真实ZnS受力、Born电荷和介电数据，计算谐波与三声子相空间，不重新进行DFT。

ZnS为两原子闪锌矿原胞，立方晶格常数5.45 Å，2×2×2超胞16原子。元素顺序为S、Zn，质量分别32.065与65.38 u；Born与力常数的原子顺序也必须一致。当前模型3³ q网格，是快速操作例。

下载实际输入、参考数值和终端脚本。

unzip fourphonon-ZnS-NAC-wave670-*.zip
cd fourphonon-ZnS-NAC
bash results.sh
bash plot.sh

以上只读参考表并绘图，不启动DFT。新的原生计算以普通账号sbatch job.slurm提交，脚本自动准备干净目录、计算并提取结果；reference-data/是已完成参考，results/才是新执行的提取结果。CPU申请cpu-8173m独占整节点1任务×56核，本小例内部1MPI×1线程，计费仍按整节点。程序为FourPhonon1.3-fix20260924 CPU版，oneAPI2024.2、Intel MPI2021.13与顺序MKL，适用AVX512 CPU；本例不报告CPU节点性能对照。

1. 原始数据与实际CONTROL#

cat source-inputs/BORN
head -18 inputs/cases/NAC-on-T300/FORCE_CONSTANTS_2ND
cat inputs/cases/NAC-on-T300/CONTROL

本例实际Born电荷为S约−2.06281946 e、Zn约+2.06281946 e，张量总和接近0；各向同性高频介电常数约6.13221705。二阶常数由同一真实受力模型通过Phonopy4.8.3生成，并按元素和周期位置对应到FourPhonon超胞索引，完整16×16个3×3块保留。

&allocations
nelements=2,natoms=2,ngrid(:)=3 3 3
&end
&crystal
lfactor=0.1
lattvec(:,1)=0.0,2.725,2.725
lattvec(:,2)=2.725,0.0,2.725
lattvec(:,3)=2.725,2.725,0.0
elements="S","Zn",types=1 2
masses(1)=32.065,masses(2)=65.38
gfactors=0 0
scell(:)=2 2 2
positions(:,1)=0.25,0.25,0.25
positions(:,2)=0.0,0.0,-0.0
epsilon(:,1)=6.132217050000005,-2.1729834716729258e-32,-2.172983471672925e-32
epsilon(:,2)=-1.3799192150979052e-32,6.132217050000005,1.864834680570967e-32
epsilon(:,3)=1.3799192150979054e-32,1.3799192150979054e-32,6.132217050000005
born(:,1,1)=-2.0628194599999996,-9.872503362906365e-17,-2.65018587562681e-32
born(:,2,1)=-9.872503362906363e-17,-2.0628194599999996,9.872503362906362e-17
born(:,3,1)=2.4659469448956835e-32,9.872503362906365e-17,-2.0628194599999996
born(:,1,2)=2.0628194599999996,9.872503362906365e-17,2.65018587562681e-32
born(:,2,2)=9.872503362906363e-17,2.0628194599999996,-9.872503362906362e-17
born(:,3,2)=-2.4659469448956835e-32,-9.872503362906365e-17,2.0628194599999996
&end
&parameters
T=300,scalebroad=1
&end
&flags
nonanalytic=.true.,espresso=.false.,tdep=.false.
isotopes=.false.,autoisotopes=.false.,nanowires=.false.
onlyharmonic=.true.,four_phonon=.false.,convergence=.false.
&end

lfactor=0.1把晶格输入的Å尺度转换到程序内部nm；born和epsilon填写完整张量。autoisotopes=.false.保留指定质量，isotopes=.false.不把同位素散射混入这个谐波例。这里是通常二阶力常数路径，tdep=.false.，不要额外填写只对极性TDEP路径适用的Ewald。

2. 四个独立目录计算#

sbatch job.slurm
ls cases/NAC-off-T300 cases/NAC-on-T300
head cases/NAC-off-T300/BTE.omega
head cases/NAC-on-T300/BTE.omega

脚本还计算两个600 K目录，各温度独立保存，具体原因见谐波多温度的目录安排。每个目录只含二阶常数，onlyharmonic=true、four_phonon=false，不需要三、四阶常数，也不计算完整弛豫时间热导。

3. 比较频率先换单位#

原生BTE.omega为rad/ps，转THz除以2π。本例关闭NAC的全部不可约q点频率与独立Phonopy结果最大差7.23×10⁻⁶ rad/ps，确认了原子映射、质量和二阶数据的对应。

启用NAC后，有限q点的光学频率确有变化，本次最大排序频率差约0.42974117 THz。图中分支按每个q点的频率顺序编号；简并或交叉附近不能因此认为编号相同就是相同本征矢。

本例精确Γ点没有指定趋近方向，三条Γ光学频率在开关前后相同。这不表示极性材料没有LO–TO劈裂；方向相关的q→0极限与精确Γ点、粗网格有限q点是不同对象。不同程序的有限q长程修正方案也不能在未核对方法时直接逐值等同。

head reference-data/ZnS-FourPhonon-NAC-frequencies.csv

这篇没有把完整二阶常数直接改名充作极性TDEP的裸短程常数；后者必须与长程减除和Ewald定义匹配，否则容易重复加入静电项。