Phonopy 位移准备与力数据收集

下载位移准备、真实力数据与QHA资料包。在已分配CPU资源内，用普通账号解压到新目录运行：

unzip phonopy-preparation-qha-wave242-*.zip -d phonopy-qha-practice
cd phonopy-qha-practice
bash run.sh
ls preparation thermal figures

脚本使用1个物理CPU，无需GPU。先实际生成位移和晶格采样，再从已提供的真实数值力与能量计算热膨胀，最后生成四张高清图；不会提交新的VASP计算。输出目录已存在时请另解压一份。

Phonopy为4.8.3，Python3.11.16、C/OpenMP由GCC13.3编译、默认phonors0.6.0后端，symfc为1.7.3。源数据来自本集群已经完成的VASP6.6.1 GPU计算：NVHPC25.3、CUDA12.8、HPC-X2.22.1、MKL2023.2。电子计算与本包的CPU后处理是两个阶段。

1. 按什么顺序做#

优化并确认结构 → 生成超胞位移 → 每个位移结构单独计算原子力 → 按原顺序收集输出 → 生成力常数 → 计算声子和热性质。

本例用两原子Si原胞，2×2×2超胞共16个原子；正、负位移各一个，位移幅度0.02 Å。复杂结构所需位移数由对称性决定，不能直接套用“只算两个”。完整声子能带和DOS流程见硅有限位移声子。

2. 原胞POSCAR#

平衡参考常规晶格常数为5.4689259126 Å；这里使用两原子的FCC原胞：

Si primitive reference
1.0
0.000000000000000 2.734462956298492 2.734462956298492
2.734462956298492 0.000000000000000 2.734462956298492
2.734462956298492 2.734462956298492 0.000000000000000
Si
2
Direct
0.00 0.00 0.00
0.25 0.25 0.25

在自己的新目录保存为POSCAR，然后执行：

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

本例生成POSCAR-001、POSCAR-002、未位移的SPOSCAR及phonopy_disp.yaml。保留原始YAML，它记录超胞、原子顺序和每次位移；后面收集力要使用同一份文件。

3. 位移结构怎样计算力#

将每个POSCAR-00*分别作为一个VASP单点输入。不要在位移结构上继续优化，否则原子会离开指定位置。源算例完整INCAR为：

SYSTEM = Si displaced supercell
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

NSW=0表示只做单点受力，这时EDIFFG不用于离子停止判断。ISYM=0用于保留实际位移后的受力。16原子位移超胞的KPOINTS：

Gamma mesh
0
Gamma
4 4 4
0 0 0

电子k网格需要按自己的体系检查收敛。资料包inputs/DFT-inputs/保留五个体积、15个已计算条件的完整INCAR/POSCAR/KPOINTS；势文件按自己的VASP授权准备。实际VASP作业通过Slurm提交，模块与提交方法见VASP环境和入口。资料包不包含完整电子计算原始输出，run.sh仅复现数值处理。

4. 收集两个真实输出#

假设自己的结果分别位于disp-001/vasprun.xml和disp-002/vasprun.xml，在保留原始YAML的目录执行：

phonopy-collect phonopy_disp.yaml \
  disp-001/vasprun.xml disp-002/vasprun.xml
head -12 FORCE_SETS

输出FORCE_SETS包含位移及逐原子受力。文件顺序必须与YAML中的位移顺序一致；两个目录不要仅凭名字相近就互换。本集群五个体积的10份真实位移输出，已按此方式逐项收集，并与原始数值力比较一致。

希望把结构、位移和力放在一份YAML中时，加--sp：

phonopy-collect --sp phonopy_disp.yaml \
  disp-001/vasprun.xml disp-002/vasprun.xml
ls phonopy_params.yaml

4.8.3加--sp后写出phonopy_params.yaml；不加时写出FORCE_SETS。 不要运行第一种形式却只寻找第二种文件。phonopy-init --sp -f ...也可以执行相关准备，phonopy-collect将收集步骤单独列成入口。

图中受力来自实际计算；右侧的三个晶格点只是新生成的几何采样，不代表三个点都完成了新DFT。资料包提供数值FORCE_SETS与带力YAML，可直接继续声子计算，不需要重新计算同一组力。

5. 读入收集结果#

python - <<'PYCODE'
import phonopy
p = phonopy.load('inputs/v100/phonopy_disp.yaml',
                 force_sets_filename='inputs/v100/FORCE_SETS',
                 produce_fc=False, is_nac=False)
print('unit cell atoms:', len(p.unitcell))
print('supercell atoms:', len(p.supercell))
print('forces:', p.forces.shape)
PYCODE

在解压目录运行，得到原胞2原子、超胞16原子、受力数组(2,16,3)。数组第一维是两个位移计算，第二维是超胞原子，第三维是Fx/Fy/Fz。

6. 随机位移与有限温度位移要分开#

只生成随机方向、固定幅度的四个结构，可以在另一个新目录执行：

phonopy-init -c POSCAR --rd 4 --random-seed 24231 \
  --amplitude .01 --dim 2 2 2

本例每个原子位移模长0.01 Å；固定种子得到相同的结构。它们是待计算的几何，不自带真实原子力，也不是有限温度MD轨迹。有限温度谐近似快照需要已有声子信息，使用模式与有限温度随机位移教程的后处理流程。

4.8.3把结构准备放在phonopy-init，后处理放在phonopy。在线官方准备命令说明可查流程；命令仍以当前安装版本帮助为准。