Phonopy 多项式势训练与声子

下载拟合、声子与短SSCHA操作包。在已分配的CPU资源中，以普通账号解压到新目录：

unzip phonopy-optional-alm-mlp-sscha-wave246.zip -d phonopy-optional
cd phonopy-optional
bash setup.sh
bash run.sh

setup.sh按包内实际使用的Linux x86_64版本清单创建个人.env，不修改公共模块；需要联网下载conda-forge包和足够的个人目录空间。run.sh依次拟合ALM及多项式势，计算势函数声子，运行两组短SSCHA，再提取数据和画图。计算均使用1个CPU线程，不需要GPU；结果放在fit/、native/、analysis/、figures/，已有结果目录不覆盖。

可选环境使用Phonopy4.8.3、pypolymlp0.22.1、symfc1.7.3、Python3.11.18、NumPy2.4.6；由conda-forge预编译包组成，实际C/C++运行库libgcc/libstdcxx为16.2.0。ALM的conda包版本为2.0.0_dev.2，Python元数据显示1.2.0；这里的ALM是ttadano/ALM力常数拟合库，不要安装PyPI中同名的金融建模ALM0.0.1。绘图另用集群phonopy/4.8.3模块中的Matplotlib。

本例复用实际Si位移—力数据：Si16，晶格常数5.43 Å，共65个构型，其中48个训练、8个验证、8个独立测试以及1个平衡参考。标签来自此前集群实际运行的GPUMD5.9.1 Si NEP势计算，不是DFT标签。固定原子序、位移和数据划分，测试集不参与拟合或选参数。模型输出中的dataset_type: vasp是库的默认格式字段，不改变本例标签来源。

1 先确定训练与验证数据#

cat fit/MLP-METRICS.json
head -n 24 fit/polymlp.yaml

48个训练构型用于拟合，8个验证构型用于正则化参数选择，最后8个独立测试构型只报告误差。不要把训练回算误差当作独立测试误差；正式模型还应按体积、温度、结构类型等考虑更严格的划分。

脚本实际使用的主要参数：

params = PypolymlpParams(
    cutoff=5.0,
    model_type=3,
    max_p=2,
    gtinv_order=2,
    gtinv_maxl=(2,),
    gaussian_params2=(0.0, 4.5, 6),
    reg_alpha_params=(-5.0, -1.0, 5),
    atom_energies={"Si": 0.0},
)

cutoff为Å，后面几项控制多项式不变量、径向基和岭回归候选范围。这是一组小体系教学设置，不是通用Si势的推荐参数。

2 保存模型，并读取实际误差#

cat fit/MLP-METRICS.json
ls fit/polymlp.yaml fit/predictions.npz

实际结果：

这表示小模型在当前构型分布上逼近了NEP标签，不是得到了0.00148 eV/Å的DFT准确度。fit.py还会将模型保存后重新载入，并比较能量、力以及一个坐标的能量负梯度；这样可以在换文件、换接口后发现单位或原子序不一致。

3 用模型生成新的位移受力与谐力常数#

脚本在native/harmonic/中调用模型。需要单独重做时，使用一个新目录：

mkdir manual-phonon
cp fit/polymlp.yaml manual-phonon/
cd manual-phonon
../.env/bin/phonopy ../native/geometry.yaml \
  --pypolymlp -d --noreadfc --amplitude 0.01 \
  --include-fc --writefc --writefc-format text --full-fc \
  --qpoints '0.25 0 0'
cd ..

模型文件应位于当前计算目录。-d要求生成有限位移，0.01 Å为位移尺度，随后模型求力、symfc重建力常数，再计算q点。

**只有--pypolymlp而不指定-d或--rd时，当前版本可以只加载模型后正常退出，不生成受力或力常数。**所以实际使用时检查FORCE_CONSTANTS、qpoints.yaml及生成的数据，不能只看终端没有报错。

cat native/harmonic/qpoints.yaml
ls native/harmonic/FORCE_CONSTANTS

本例q=(0.25,0,0)六个频率约为2.152355、2.152355、6.124735、13.979023、14.365756、14.365756 THz。这些是本例多项式势的谐结果，不能与另一晶格常数、另一DFT参数的频率直接当作同条件误差比较。

下一步看300 K短SSCHA流程；训练数据整理方式另见结构、能量、力与应力HDF5。软件来源见pypolymlp官方项目。