Phonopy极性晶体：ZnS的LO–TO分裂

下载完整输入、原始输出和后处理脚本。解压后进入phonopy-workflows子目录。在已分配的CPU资源中运行bash run.sh，可重建本文全部后处理结果；这不会重新提交VASP任务。

为什么还需要Born电荷和介电张量#

极性晶体的长波光学模包含长程电场响应。有限超胞受力给出短程近似，再结合Born有效电荷和电子介电张量进行非解析修正（NAC），才能描述Γ附近的LO–TO分裂。

本例新计算ZnS 2×2×2超胞的4个位移构型，每个16原子、位移±0.02 Å，均在集群P100完成电子收敛。Born电荷和电子介电张量来自同一原胞的已有真实VASP DFPT输出。用的是电子介电响应，不能误填包含离子响应的静态介电常数。

1 结构与位移计算#

原胞POSCAR如下，元素顺序为S、Zn：

S  Zn
 1.0000000000000000
     0.0000000000000000    2.7250000000000001    2.7250000000000001
     2.7250000000000001    0.0000000000000000    2.7250000000000001
     2.7250000000000001    2.7250000000000001    0.0000000000000000
 S   Zn 
   1   1
Direct
  0.2500000000000000  0.2500000000000000  0.2499999999999999
  0.0000000000000000  0.0000000000000000 -0.0000000000000000

完整受力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.

超胞K点为3×3×3 Γ中心。四个目录zns-000至zns-003保存实际输入和输出；赝势按S、Zn的顺序准备，程序不会自动替你纠正错序。

新结构的生成命令可写成：

module load phonopy/4.8.3
phonopy-init -c POSCAR-zns --dim 2 2 2 --pa P -d --pm --amplitude 0.02

这里显式--pa P保留已选原胞基矢，便于对应Born张量与原子顺序。对其他结构不要机械照抄。

2 读取Born电荷与电子介电张量#

配套脚本从zns-born/vasprun.xml解析张量，检查电荷和规则与电子介电矩阵，并将实际受力组装为力常数：

python analyze_zns.py
cat zns-analysis/BORN
cat zns-analysis/RESULTS.json

BORN第一行为单位换算因子，第二行是3×3电子介电张量，随后每行是一个原子的3×3Born张量。行顺序必须与原胞原子一致。完整实际文件：

14.399652
6.132217050000005 -1.3799192150979052e-32 1.3799192150979054e-32 -2.1729834716729258e-32 6.132217050000005 1.3799192150979054e-32 -2.172983471672925e-32 1.864834680570967e-32 6.132217050000005
-2.0628194599999996 -9.872503362906363e-17 -6.162975822039155e-33 -9.872503362906363e-17 -2.0628194599999996 9.872503362906362e-17 -2.0543252740130515e-33 9.872503362906362e-17 -2.0628194599999996
2.0628194599999996 9.872503362906363e-17 6.162975822039155e-33 9.872503362906363e-17 2.0628194599999996 -9.872503362906362e-17 2.0543252740130515e-33 -9.872503362906362e-17 2.0628194599999996

3 指定趋近Γ的方向并比较#

脚本沿原胞倒易分数方向[0,0,1]取61个q点，从0到0.1，分别开启和关闭NAC。关键API为：

phonon.nac_params = {
    'born': born, 'dielectric': epsilon, 'factor': 14.399652
}
phonon.run_qpoints(qpoints, nac_q_direction=[0, 0, 1])

这是原胞倒易坐标方向，不等于直接把常规晶胞的笛卡尔z轴代入。

zns-analysis/gamma-line.csv给出所有q点的12列频率，便于自己画图。声学Γ模的约10⁻⁷ THz数值残差不代表本例发现结构不稳定。本例展示流程；超胞、K点、晶格常数和DFPT参数均需针对生产结果检查收敛。