VTST：动力学矩阵抽取与自由度重排

dymextract.pl 可以从较大的动力学矩阵中选出部分自由度；dymreorder.pl 用于把来自不同位移顺序的矩阵整理到共同顺序。本例沿用 NH₃ 原生动力学矩阵 的实际结果，不需要重新计算 VASP。

1. 解压并先看结果#

下载输入、实际输出、CSV 与绘图脚本，保存为 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

本例调用 VTST scripts 1040 的 Perl 工具，图表使用 ASE 3.29.0、NumPy 和 Matplotlib。analysis-read/matrix.csv 保存原矩阵、抽取结果和重排结果；figures/ 保存图片。

2. 先明确矩阵的自由度顺序#

原结构的原子顺序为 H、H、H、N，完整矩阵为 12×12。inputs/DISPLACECAR 为：

0.015 0.015 0.015
0.015 0.015 0.015
0.015 0.015 0.015
0.015 0.015 0.015

矩阵顺序按“原子逐行、每行 x/y/z”展开：1x、1y、1z、2x……4z。freq.mat 是质量加权矩阵，单位为 eV/(Ų·amu)，不要与未质量加权的力常数矩阵混用。

3. 抽取第一个原子的三个方向#

inputs/DISPLACECAR.subset 为：

0.01 0.01 0.01
0    0    0
0    0    0
0    0    0

第一个位移文件描述原矩阵的自由度，第二个描述要保留的自由度：

mkdir extract-new
cd extract-new
cp ../inputs/* .
perl "$VTSTSCRIPTS_ROOT/dymextract.pl" \
  DISPLACECAR DISPLACECAR.subset freq.mat > selected.mat
cat selected.mat
cd ..

实际输出：

   5.3258   0.0985  -0.0411
   0.0985  35.6328 -12.5352
  -0.0411 -12.5352   6.3904

脚本按 4 位小数输出，因此与源文件最后几位不同是输出精度所致。这个 3×3 块是“其余自由度固定时的子空间矩阵”，不是整个 NH₃ 的完整振动谱，也不是消去其余原子响应后的有效矩阵。换体系时，两个 DISPLACECAR 都保留全部原子行。

4. 把 z/x/y 排序恢复为 x/y/z#

资料包把同一 3×3 块的行列同时置换成 z、x、y，得到 new-order.mat。其位移列表分成三个文件：DISPLACECAR.0 仅 1z 非零，.1 仅 1x 非零，.2 仅 1y 非零；每个文件仍有四行。

mkdir reorder-new
cd reorder-new
cp ../inputs/* .
perl "$VTSTSCRIPTS_ROOT/dymreorder.pl" \
  1 DISPLACECAR.xyz \
  3 DISPLACECAR.0 DISPLACECAR.1 DISPLACECAR.2 \
  new-order.mat > reordered.mat
cat reordered.mat
cd ..

第一个 1 是参考位移文件数量，后面的 3 是待整理位移文件数量。脚本依据相同的“原子编号＋方向”建立映射，同时变换行和列。输出与上面的 selected.mat 一致。

两个集合必须包含相同自由度。原子数相同但顺序不同，或者一边是质量加权矩阵、一边不是，都不能靠重排自动修正。

5. 一次执行本例所有工具#

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

run.sh 只处理本资料包的小矩阵和结构，不启动 VASP；结果写入新的 work-new/，已有同名目录时停止，避免覆盖旧结果。后续可阅读位移拟合与误差分析和多模式动画。