VASP与Phonopy：ZnS的LO–TO分裂

极性晶体在Γ点附近的纵光学LO和横光学TO频率可以不同。Born有效电荷与电子介电张量提供长程电场项；把它加入动力学矩阵，才是本篇实际执行的非解析修正NAC。这里复用鸢算科技集群已经完成的ZnS原胞DFPT，依次读取数据、运行Phonopy、检查方向和模式，再画结果。

本例只有两原子原胞的Γ点力常数，用于Γ方向极限；它不是整条布里渊区色散。 需要一般q点色散时，仍应准备足够大的超胞力常数或匹配的完整q网格计算。

1. 下载后先看参考图#

下载ZnS原生响应数值与NAC实操包，解压到新目录：

python3 -m zipfile -e vasp-ZnS-Gamma-NAC-wave303-*.zip .
cd vasp-ZnS-Gamma-NAC-wave303
bash preview.sh
bash run.sh

preview/重画两张参考图；run.sh调用真实Phonopy4.8.3完成三个方向的Wang Γ点NAC，再独立构造Born/介电修正矩阵核对。results/保存完整6×6动力学矩阵、复数本征矢、CSV、BORN、FORCE_CONSTANTS和Phonopy YAML，plots/保存汇总图。只需1个CPU后处理，不重新运行DFT；再次运行请换全新解压目录。

包内reference/ZnS-DFPT.npz是从实际DFPT导出的完整自有数值，不含原始完整XML、OUTCAR、POTCAR、波函数或厂商程序。实际VASP输入位于inputs/ZnS-DFPT/。

2. 从头计算时准备什么#

本例闪锌矿ZnS，常规晶格常数5.45 Å，两原子fcc原胞，固定晶胞与坐标。S在前、Zn在后；POTCAR必须按相同顺序合法准备。

cd inputs/ZnS-DFPT
ls INCAR POSCAR KPOINTS run.slurm
nano INCAR

完整POSCAR：

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-response
ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
EDIFFG = -0.02
IBRION = 8
NSW = 1
ISIF = 3
ISMEAR = -5
SIGMA = 0.05
LREAL = .FALSE.
LASPH = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
NELM = 150
ALGO = Normal
ISYM = 2
LCALCEPS = .FALSE.
LEPSILON = .TRUE.
NFREE = 2
POTIM = 0.015

完整KPOINTS：

Gamma mesh
0
Gamma
7 7 7
0 0 0

IBRION=8是DFPT响应，不是结构弛豫；LEPSILON请求电子介电与Born响应，NSW=1用于进入相应响应任务。要研究材料平衡结构，先独立优化，再保持响应设置相容。EDIFF=1E-6是本例阈值，微小频差还需针对自己的体系检查电子、网格和k点精度。

本例原生DFPT使用VASP6.6.1 CPU程序、Intel oneAPI2024.2、Intel MPI2021.13与MKL2024.2。附件模板在可用V10016分区调用vasp/6.6.1-8173m的CPU vasp_std，一个MPI进程、一个线程；GPU随该模板分配，并没有执行这条CPU程序的DFPT。不要因加载某个GPU分区就把结果称为GPU加速。

ls POTCAR
sbatch run.slurm
squeue -u "$USER"

结束后，检查实际响应输出：

grep 'EDIFF is reached' OUTCAR | tail -n 3
grep -n 'BORN EFFECTIVE CHARGES' OUTCAR
grep -n 'DIELECTRIC TENSOR' OUTCAR
ls -lh vasprun.xml

3. 从自己的VASP结果提取#

回到资料包根目录：

cd ../..
module purge
module load phonopy/4.8.3
python read_vasp_response.py inputs/ZnS-DFPT --output my-response
python analyze_nac.py --data my-response --output my-results
python plot.py --data my-results --output my-plots

第一条Python命令调用安装环境中Phonopy的VASP读取器，并提取原生Hessian、Born、电子/离子介电、结构与质量；后两条做NAC和画图。目标目录已有结果时，新建目录保存，避免覆盖。

VASP6.6.1 XML中的Hessian单位为THz²且已质量加权。不能把它直接写成eV/Ų的FORCE_CONSTANTS；读取器会恢复符号、原子质量与单位。本例按原生PAW记录采用S=32.066、Zn=65.390 amu，没有擅自改成其他同位素或软件默认质量。

当前脚本专门对应本例两原子ZnS与已有检查条件。更换材料、原胞、原子顺序或质量时，需要一起修改模型与相应检查，不能只把一个文件名替换掉。

4. 方向必须用正确坐标#

NAC的长波极限取决于趋近Γ的q方向。Phonopy使用实际原胞的倒易分数坐标；它与立方常规笛卡尔方向不是同一组数字。

脚本用实际原胞矩阵转换方向，随后执行：

ph.nac_params = {
    'born': born,
    'dielectric': epsilon_electronic,
    'factor': 14.39965172592227,
    'method': 'wang',
}
result = ph.run_qpoints(
    [[0, 0, 0]],
    nac_q_direction=q_direction_fractional,
    with_eigenvectors=True,
    with_dynamical_matrices=True,
)

这里的转换系数对应当前VASP/Phonopy的eV、Å与电子电荷单位；born形状为原胞原子数×3×3，dielectric是电子介电张量ε∞，不能把含离子响应的静态总ε₀代入同一个NAC公式，否则会重复计入晶格屏蔽。

5. 本次实际结果#

立方ZnS三个方向的光学频率一致，但LO位移向量会随q方向旋转。两个TO模式保持横向；LO模式中S与Zn沿q反向运动。图中箭头来自实际本征矢，统一缩放方便观察，不是真实温度下的振幅。

Born对角元为S −2.06281946 e、Zn +2.06281946 e；它们的和为零。电子介电对角元约6.13221705。原子投影电荷、名义化合价与Born有效电荷定义不同，不能互相替代。

用同一次DFPT的电子与离子介电响应作Lyddane–Sachs–Teller关系对照，得到LO约326.416837 cm⁻¹，与实际NAC矩阵结果差约0.083356 cm⁻¹，即0.0255%。保留这两个有限精度结果，不把其中一个改成另一个。

完整矩阵还用Born/介电的秩一长波修正独立重建，与Phonopy矩阵相对差约1.4×10⁻¹⁶；这是同一输入数据的实现核对，不代表电子结构、超胞或材料参数已经充分收敛。

6. 近零声学模和后续使用#

本次未施加声学和规则ASR，也未修正Born电荷；原生声学数值约−0.498、−0.498、−0.396 cm⁻¹。它们靠近数值零，不能只看负号就宣布材料失稳。先检查电子/力常数精度、对称性和声学模形状，再决定是否在明确记录的前提下施加ASR；不要隐藏原值。

cat results/frequencies.csv
cat results/BORN
head results/FORCE_CONSTANTS
ls results/full-dynamical-matrices.npz results/phonopy_params.yaml

frequencies.csv列出方向、TO、LO和LST对照；NPZ保留三组完整复数矩阵与本征矢。Γ点DFPT数据可以直接教学演示LO–TO极限，但不能用这个1×1×1原胞模型宣称整条声子色散、声子DOS、热导或动力学稳定性。

继续阅读有限位移晶体声子、介电与Born响应、红外强度及声子投影与单位。方法参照Phonopy DFPT流程和非解析修正公式。