VTST：位移矩阵拟合与误差分析

dymfit.pl 对不同位移对应的矩阵做多项式拟合；dymanalyze.pl 比较逐步加入数据或固定参考拟合时的变化。本例先用明确给定的解析数据认识命令，再说明怎样替换成自己的有限位移计算。

1. 准备与数据来源#

下载完整操作资料包，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd vtst-matrix-tools
module load ase/3.29.0 vtstscripts/1040
bash read.sh
bash plot.sh

H0.mat 来自真实 NH₃ 质量加权动力学矩阵的 1x/1y/1z 子块。为了让拟合结果有独立已知答案，四个输入明确按下面的式子生成：

H(d) = H₀ + d²C，d = 0.01、0.02、0.03、0.04 Å。

curvature.mat 为：

2.0 0.1 0.2
0.1 3.0 0.3
0.2 0.3 4.0

这里 C 的单位与 d² 相配，使 H(d) 仍为 eV/(Ų·amu)。这些 step-*.mat 是解析构造的工具练习数据，不是四次新的 DFT 计算，也不能用来证明某个位移步长已经收敛。

2. 拟合的对象是什么#

脚本会先把输入矩阵乘以位移 d，再拟合 d×H(d) 随 d 的变化，并额外加入零点。返回一次项系数，即零位移处的导数；本例的理论答案为 H₀。

需要注意：不是把所有矩阵简单平均，也不是直接拟合 H(d) 的一次项。各矩阵应有相同原子、自定义自由度顺序、质量加权约定和单位；顺序不同先用矩阵重排工具处理。

3. 对小位移使用个人缩放副本#

本环境的旧数学库使用固定绝对消元阈值，小位移的高次幂会造成三次拟合精度不足。资料包的 setup.py 只在当前目录生成个人副本：把多项式变量改为 d/0.01，输出导数时再除以 0.01。物理位移和输入矩阵保持原值，公共脚本不变。

python3 setup.py
mkdir fit-new
cd fit-new
cp ../inputs/* .
perl ../personal/dymfit-scaled.pl 3 \
  0.01 step-0.01.mat 0.02 step-0.02.mat \
  0.03 step-0.03.mat 0.04 step-0.04.mat > fitted.mat
cat fitted.mat
cd ..

实际输出：

5.3257819300   0.0985248600  -0.0411047100
0.0985248600  35.6327578800 -12.5351893600
-0.0411047100 -12.5351893600 6.3903943900

与 H₀ 一致。3 是拟合最高阶数；增加阶数需要足够的独立位移数据，并不保证更可靠。该个人辅助程序针对已核对的 scripts 1040 源文件生成，遇到其他源文件会提示版本不匹配。

4. 逐步加入位移数据#

mkdir analysis-new-tools
cd analysis-new-tools
cp ../inputs/* .
perl "$VTSTSCRIPTS_ROOT/dymanalyze.pl" 0 0 \
  0.01 step-0.01.mat 0.02 step-0.02.mat \
  0.03 step-0.03.mat 0.04 step-0.04.mat > progressive.txt
cat progressive.txt

第一个 0 表示逐步加入位移数据，观察拟合斜率如何变化。第二个 0 是“只输出绝对值大于该阈值的矩阵元素”的筛选阈值，不是电子收敛或受力阈值。本例输出 9 个元素。

5. 固定前两个位移作为参考#

perl "$VTSTSCRIPTS_ROOT/dymanalyze.pl" 2 0 \
  0.01 step-0.01.mat 0.02 step-0.02.mat \
  0.03 step-0.03.mat 0.04 step-0.04.mat > reference.txt
cat reference.txt
cd ..

第一个 2 表示用前两个位移点与脚本加入的零点拟合一条直线，然后比较其余数据与该直线。dymanalyze 的这两种分析均为线性拟合，与上面 dymfit 的三次拟合不同。

本例固定参考分析的平均偏差依次约为 0.029%、0.007%、0.077%、0.176%；最后一列 num points=9 表示纳入统计的矩阵元素数，不是原子数或 VASP 运行次数。相对误差接近零分母时会变敏感，筛选阈值应结合矩阵实际量级选择。

6. 替换为自己的实际计算#

先固定同一平衡结构，使用多个位移大小分别收集力并生成矩阵；核对自由度顺序、单位及质量加权方式，再把命令里的 step-*.mat 换成实际输出。需要比较步长时，应保证电子收敛、ENCUT、k 点等条件一致，同时检查噪声和非谐效应。

bash run.sh
bash read.sh work-new analysis-new
bash plot.sh analysis-new

这里的快捷命令仍只执行附带的工具练习；真正的有限位移 VASP 流程见原生动力学矩阵实操。