SIESTA分子振动：有限位移、Vibra与频率提取
用水分子演示原生力常数和Vibra操作，比较两种位移，保留全部频率并区分内部振动与外部转动。
阅读本篇目录
先优化结构，再用有限位移计算力常数，最后用vibra得到频率和振型。这是孤立水分子例子，和晶体的声子色散不同：这里仅处理Gamma点的9个自由度。
1. 准备输入并提交#
unzip siesta-water-wave201.zip -d water-vib
cd water-vib
sbatch job.slurm
包内包含优化、动力学和振动的完整流程；只需提交一次。结构、PBE/DZP、150Ry、12Å盒和电子收敛的设置见前一篇。环境为SIESTA5.4.2、Intel oneAPI2024.2、Intel MPI2021.13、MKL顺序LP64，配套sisl0.16.4。
2. 从优化结构计算力常数#
程序先读取cases/water-opt/h2o.XV，再为fc010和fc020写入相同的优化坐标。以下是第一组运行参数：
MD.TypeOfRun FC
FC.Displacement 0.01 Ang
FC.First 1
FC.Last 3
FC.First/Last是要位移的原子编号范围；水有3个原子，三个方向分别做正负位移。第二组把位移改成0.02Å，其他设置相同。位移过小会放大力的数值噪声，过大则可能引入非谐影响，应结合实际体系比较。
ls cases/fc010/h2o.FC
head cases/fc010/h2o.FC
h2o.FC保留原生力常数数据。包内分析直接读取它，不修改原始文件。
3. 用Vibra提取频率与振型#
Vibra的结构输入需要质量，以及它支持的坐标格式名：
SystemLabel h2o
Vibra.FC h2o.FC
NumberOfAtoms 3
LatticeConstant 1 Ang
%block LatticeVectors
12.0 0.0 0.0
0.0 12.0 0.0
0.0 0.0 12.0
%endblock LatticeVectors
AtomicCoordinatesFormat NotScaledCartesianAng
%block AtomicCoordinatesAndAtomicSpecies
5.9999996334 5.9999999996 5.9394271851 1 15.9994
6.7703988056 6.0000000261 6.5401437938 2 1.00794
5.2296015302 5.9999999943 6.5401435769 2 1.00794
%endblock AtomicCoordinatesAndAtomicSpecies
SuperCell_1 0
SuperCell_2 0
SuperCell_3 0
Eigenvectors true
BandLinesScale ReciprocalLatticeVectors
%block BandLines
1 0.0 0.0 0.0 Gamma
%endblock BandLines
上面的坐标块来自本例优化后的结构；使用自己的体系时，应换成对应优化结果，不能继续沿用这些水分子坐标。O/H质量分别为15.9994/1.00794。注意Vibra支持的格式别名与SIESTA并不完全相同，应使用这里的NotScaledCartesianAng。
以下是自动流程实际采用的后处理命令；如需重新提取，可在计算完成后运行：
module load siesta/5.4.2-8173m-mpi
cd cases/fc010
vibra < vibra.fdf > vibra.out
grep 'omega=' vibra.out
head h2o.bands
head h2o.vectors
检查确实生成9个频率和对应振型，不能只看程序退出码。h2o.bands存频率，h2o.vectors存振型，vibra.out还有零点能等输出。对于含有数值负频的分子，本例不把其零点能输出作为完整热力学结论。
4. 看懂正频、近零频和负频#
非线性三原子分子有3个内部振动自由度，其余是平移和转动。当前计算得到3个正频内部振动、3个近零平移频率，以及3个明显为负的外部模式。
这里没有把负频删掉：全部9个模式都在图和CSV中。对质量加权振型与刚体转动子空间作投影后，3个负频模式的转动重合度均超过99.99%。这说明它们主要对应当前网格和几何精度下没有准确回到零的外部转动，不能直接解读为水分子内部结构不稳定。若要得到准确外部零模，仍需更严格的几何、网格和其他参数收敛检查；不能对任意体系的负频都套用同一解释。
| 位移 | 弯曲（cm⁻¹） | 伸缩1（cm⁻¹） | 伸缩2（cm⁻¹） |
|---|---|---|---|
| .01Å | 1594.36 | 3600.70 | 3757.28 |
| .02Å | 1585.54 | 3601.95 | 3762.21 |
两种位移的三个内部频率最大相对变化为0.553%。此处只比较位移幅度，没有完成基组、盒长与网格收敛，因此不作为高精度光谱预测。
5. 读取CSV，独立核对力常数#
head analysis/fc010-all-modes.csv
head analysis/fc020-all-modes.csv
cat analysis/SUMMARY.json
包内分析按原生Vibra的平移求和规则对称化力常数，再构造质量加权动力学矩阵，独立对角化并与原生频率比较。两种位移的三个内部频率都与原生结果一致；原始负频和力常数仍保留，不用人工改值代替收敛。
官方振动与声子教程还区分分子振动与晶体声子：晶体需要超胞力常数与波矢路径，不能直接照搬本例的分子Gamma点设置。