Phonopy硅声子：从位移受力到能带与DOS

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

工作流先后顺序#

先确定平衡结构与计算精度，再生成超胞和位移；为每个位移结构运行VASP单点计算，收集受力；最后由Phonopy生成力常数、声子能带和DOS。不能只生成位移结构就开始画声子图。

本例硅原胞2原子，2×2×2超胞16原子，位移0.02 Å。复用的是鸢算科技集群此前实际完成的VASP 6.6.1 PBE PAW64受力，不是经验势或网上示例力常数。ENCUT=550 eV、EDIFF=1E-6。GPU构建使用NVHPC 25.3、CUDA 12.8、HPC-X 2.22.1和MKL 2023.2；Phonopy 4.8.3使用独立Python 3.11环境及C/OpenMP与phonors后端。

1 生成位移#

在自己的新目录准备原胞POSCAR，使用：

module purge
module load phonopy/4.8.3
phonopy-init -c POSCAR --dim 2 2 2 -d --amplitude 0.02 --pm
ls POSCAR-* phonopy_disp.yaml

--pm显式生成正负位移；结构对称性决定需要计算多少个位移。文件编号从工具实际输出读取，不要假定所有材料都只有两个。

本下载包的历史硅位移结构在si/calculations/si-phonon-v100-000和...-001中，对应的位移描述为si/si-phonon-v100.yaml。本例完整INCAR：

SYSTEM = si-phonon-v100-000
ENCUT = 550
PREC = Accurate
EDIFF = 1E-6
EDIFFG = -0.02
IBRION = -1
NSW = 0
ISIF = 2
ISMEAR = 0
SIGMA = 0.05
LREAL = .FALSE.
LASPH = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
NELM = 150
ALGO = Normal
ISYM = 0

KPOINTS：

Gamma mesh
0
Gamma
4 4 4
0 0 0

为每个位移文件建立独立VASP目录，复制相同INCAR、KPOINTS及有授权的Si赝势，将位移结构命名为POSCAR，然后用对应GPU模板提交。NSW=0保持位移，不能把位移结构重新优化回去。

2 按顺序收集受力#

从本包根目录执行以下命令可重建FORCE_SETS：

mkdir -p my-force-sets
cd my-force-sets
cp ../si/si-phonon-v100.yaml phonopy_disp.yaml
phonopy-init -f \
 ../si/calculations/si-phonon-v100-000/vasprun.xml \
 ../si/calculations/si-phonon-v100-001/vasprun.xml
head -n 12 FORCE_SETS
cd ..

文件顺序必须对应位移编号。先检查每个VASP任务正常结束并达到电子收敛，不能把中途退出的XML直接拿来组装力常数。本例配套analyze_si.py还移除每个位移构型的平均力漂移并对称化力常数；新提取的五个体积结果与旧力常数最大差异约5×10⁻¹⁶ eV/Ų。

python analyze_si.py
ls si/analysis/v100/

得到包含力常数的phonopy.yaml及FORCE_SETS。16原子超胞对长程相互作用的覆盖有限，生产计算仍要检查超胞与电子采样收敛。

3 画声子能带和DOS#

mkdir -p my-band
cd my-band
phonopy ../si/analysis/v100/phonopy.yaml \
 --band '0 0 0 0.5 0 0.5 0.5 0.25 0.75' \
 --band-points 51 -p -s
ls band.yaml band.pdf
cd ..
mkdir -p my-dos
cd my-dos
phonopy ../si/analysis/v100/phonopy.yaml --mp '12 12 12' --dos -p -s
ls total_dos.dat total_dos.pdf

上述短路径为Γ–X–W，坐标属于本例fcc原胞的倒易分数坐标。下图由配套Python脚本绘制更完整的Γ–X–W–K–Γ–L路径与16³网格DOS；不要把两套路径长度混为同一个输出。

每个q点有6支，因为原胞包含2个原子。Γ附近三支声学模应趋近零；微小数值负值要结合绝对量级判断，明显负频率则应检查结构与计算收敛。

4 看热力学数据#

phonopy si/analysis/v100/phonopy.yaml --mp '12 12 12' \
 -t --tmin 0 --tmax 800 --tstep 100 --exclude-gamma-acoustic

得到thermal_properties.yaml。自由能单位kJ/mol、熵与热容单位J/(mol·K)，归一化对象按该Phonopy原胞定义；比较不同原胞大小时需统一。继续看五体积准谐流程和模振动与随机位移。