GPUMD 模态本征矢与结构准备

GKMA和HNEMA都需要先知道体系的振动模态。本篇先完成“结构→谐振动本征矢→模态计算输入”这一步，再分别进入GKMA或HNEMA。

1 下载与查看#

下载Si模态分析输入、真实结果和绘图脚本，解压后先查看已有结果：

unzip gpumd-Si216-GKMA-HNEMA-wave229*.zip -d modal-example
cd modal-example
sha256sum -c MANIFEST.sha256
bash view.sh
ls analysis/figures

view.sh只处理附件中的参考结果，不提交新计算。新计算使用普通账号，一任务一张GPU、一个CPU核。本例环境为GPUMD5.9.1、GCC13.3、CUDA12.9；绘图使用ASE3.29.0、NumPy和Matplotlib。附件提供P100与V100 16G模板，套餐用户可用 yskit 取得匹配账户的分区模板，保留其中的程序启动命令。

2 本例结构是什么#

Si金刚石常规胞重复3×3×3，共216个原子，晶格常数5.431 Å，立方盒长16.293 Å，三个方向周期。使用附件里的 Si_Fan_2019 Tersoff-mini势；固定晶胞，先最小化原子位置。

左侧是从原生本征矢读取的一组模态方向，箭头作了可视化缩放，不能当成热振动的实际振幅；右侧来自真实MD末帧坐标。该216原子超胞在自身Γ点有3N=648个模态，是原胞各波矢折叠后的超胞模态，不是原胞Γ点有648条声子支。

3 先提交本征矢计算#

python3 prepare.py
cd cases/basis-p100-d0.005
cat job.slurm
sbatch job.slurm
squeue -u "$USER"

模板中的 run_native.py 生成完整 model.xyz、势文件、kpoints.in 和 run.in。真实输入如下：

replicate 1 1 1
potential potential.txt
minimize fire 0.000001 5000
compute_phonon 0.005

kpoints.in 本例仅包含超胞Γ点：

0 0 0 G

查看结构头部及前几个原子：

head -8 model.xyz
cat run.in
ls -lh eigenvector.out

216
Lattice="16.293 0 0 0 16.293 0 0 0 16.293" Properties=species:S:1:pos:R:3:group:I:1 pbc="T T T"
Si 0 0 0 0
Si 1.35775 1.35775 1.35775 0
Si 0 2.7155 2.7155 0
Si 1.35775 4.07325 4.07325 0
Si 2.7155 0 2.7155 0
Si 4.07325 1.35775 4.07325 0

compute_phonon 0.005 中的0.005是有限位移量，单位Å，不是MD时间步。附件另有0.01 Å对照。对多体势，力常数截断可大于势的直接截断；本例程序打印的力常数截断为8 Å，16.293 Å盒长满足两倍截断的要求。修改材料或势时应重新核对盒长，不要机械复用这个最小尺寸。有限位移与盒长要求

4 检查原生二进制文件#

本例Γ点路径实际写出 eigenvector.out，文件大小1,682,208字节。它是二进制float32数据，不要用文本编辑器修改，也不要期待每种声子计算路径都输出同样的文件名组合。

module load ase/3.29.0
python3 - <<'PY'
import numpy as np
x = np.fromfile("eigenvector.out", dtype=np.float32)
n = 648
assert x.size == n + n*n
w2 = x[:n]
frequency = np.sign(w2)*np.sqrt(abs(w2))/(2*np.pi)
print("modes:", n)
print("frequency range (THz):", frequency.min(), frequency.max())
PY

这几个接近零的声学频率保留了数值误差的符号，不能仅凭它们判定结构出现明显软模。本征矢正交误差约1.5×10⁻⁸。位移对照适合检查本征矢输入的稳定性，不能代替势函数本身的适用性研究。

5 将本征矢传给MD#

手动运行模态输入时，把对应文件按二进制原样复制：

cp ../basis-p100-d0.005/eigenvector.out eigenvector.in
sha256sum eigenvector.in ../basis-p100-d0.005/eigenvector.out

上面命令应在匹配的MD案例目录执行。附件runner会自动完成这一步。参考结果为节省空间，只在四个 reference-results/basis-* 中保存共享本征矢；各MD目录不重复存放同一份文件。

必须保持原子数、元素质量、原子排序、晶胞与势文件相匹配。不能先算本征矢，再对MD结构排序、删原子或改盒长。新结构应重新计算本征矢。

最后得到 eigenvector.out，供后续两篇的 eigenvector.in 使用。本征矢文件大小随原子数平方增长，先用小结构跑通，再估算显存、内存与输出空间。