py4vasp光学谱与振动模式读取

下载本篇输入、完整脚本、数据与高清图。解压后在新目录运行，reference保存鸢算科技实际运行的参考结果。

1 光学谱需要先计算相应输出#

py4vasp不能从普通SCF文件自动补出未计算的光学谱。本次Si计算在VASP 6.6.1 HDF5版本中设置：

LH5 = .TRUE.
LOPTICS = .TRUE.
NBANDS = 32
NEDOS = 1001
CSHIFT = 0.1

其余结构、ENCUT=550 eV、电子收敛1E-6和11×11×11 k点见下载包输入。读取复介电张量：

import py4vasp
c = py4vasp.Calculation.from_path('input/si-optics-hdf5')
d = c.dielectric_function.read()
energy = d['energies']
epsilon = d['dielectric_function']
print(epsilon.shape)

数组形状为(3,3,1001)，最后一维为光子能量网格。epsilon[0,0].real和.imag分别为xx分量实部、虚部。

零频实部εxx约 13.82101；立方结构的三个对角分量最大差小于8×10⁻¹⁵。该谱来自PBE独立粒子计算，没有GW带隙修正或BSE激子效应；空带数、k点和展宽也需按研究目标收敛。

2 Γ点振动模式#

本次Si原胞采用IBRION=8、NSW=1进行DFPT，并保存HDF5。读取方式：

c = py4vasp.Calculation.from_path('input/si-dfpt-hdf5')
modes = c.phonon.mode.read()
fc = c.force_constant.read()['force_constants']
print(modes['frequencies'], modes['eigenvectors'])

2个原子给出6个模式。py4vasp的模式频率数组以eV表示，乘1000为meV；虚频以虚部表示，不能直接丢弃虚部后宣称没有虚频。

三个光学模式均约 62.31207 meV，其余为数值零模；最大虚部小于2×10⁻¹³ meV。本例仅计算Γ点，不能据此判断全布里渊区动力学稳定性，也不是完整声子色散。

3 力常数矩阵的符号和求和规则#

本例原生6×6矩阵对称，按相邻原子求和后接近零，符合平移不改变能量的要求。需要留意读取接口返回的原生符号约定：本地force_constant.eigenvectors()实现使用-force_constants进行对角化。不要只因原生矩阵出现负本征值，就将其解释为虚频；应结合质量加权、符号约定以及phonon.mode输出判断。

完整模式、原生矩阵和检查结果在下载包，plots.py可重新生成两幅图。更换体系时，应先优化结构，再选择所需q点、超胞或DFPT方案。

参考：py4vasp官方物理量目录、力常数接口。