CREST：独立数值Hessian与振动频率

CREST的独立 numhess 入口直接对给定结构做数值Hessian。本篇用已优化的NH₃说明输入、提交、矩阵文件和频率读取；它与“先优化再求Hessian”的 ohess 是不同入口。

1. 下载并读取实际结果#

下载完整输入、实际输出、读取和绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd crest-numhess-sort
bash read.sh
bash plot.sh

默认读取附带结果，CSV写入 analysis-read，图片写入 figures。这些命令不会重新提交计算。

2. 先准备优化后的结构#

cat inputs/ammonia-numhess/input.xyz

4
 energy: -4.426244036894 gnorm: 0.000000570487 xtb: 6.7.1 (edcfbbe)
N            0.00000000001755        0.00000023323951        0.10039380387376
H           -0.00000000006857        0.94388664568029       -0.26644287397807
H            0.81742959044634       -0.47194293951421       -0.26644296494362
H           -0.81742959039533       -0.47194293940559       -0.26644296495207

本例结构由xTB 6.7.1的GFN2模型优化得到，原始优化命令为：

xtb input.xyz --gfn 2 --chrg 0 --uhf 0 --opt vtight --acc 0.1 --parallel 2

优化应放在独立目录，以 xtbopt.xyz 作为后续Hessian的 input.xyz。资料包已提供该优化结果，不需要为读取示例再次优化。

3. NUMHESS输入与提交#

input = "input.xyz"
runtype = "numhess"
threads = 2
[calculation]
elog = "energies.log"
[[calculation.level]]
method = "gfn2"
chrg = 0
uhf = 0
accuracy = 0.01
rddip = true

accuracy 控制电子结构计算精度；threads=2 是本例线程设置。rddip 不代表输出的每列都已经是物理红外强度，见第5步。

本例使用未改源码的个人CREST 3.0.3构建，GNU Fortran 13、MKL，并启用内置GFN2计算接口；外部xTB 6.7.1用于独立单点。首次编译应在已分配的CPU资源内执行：

bash build-private.sh

程序位于 ~/.local/ys-crest-tools/bin/crest。核心命令：

~/.local/ys-crest-tools/bin/crest input.toml > native.log 2>&1

整包的新任务入口为：

bash submit.sh

此资料包也包含独立结构排序例子，会在互相隔离的子目录依次运行。

4. 找到矩阵和频率文件#

head -12 reference-results/ammonia-numhess/numhess1
cat reference-results/ammonia-numhess/vibspectrum1

4原子体系的Cartesian Hessian为12×12，单位Eh/Bohr²。非线性分子扣除平移和转动后有6个内部振动模式。本次内部频率依次为：

模式序号按原生频率表排序。本例六个内部模式均为正；没有做实验谱或高精度从头算频率标定。

5. 为什么不能直接把1000画成红外强度#

本入口的 vibspectrum1 给内部模式填入了统一的1000占位值。这些值不能用于红外峰强比较，即使旁边存在偶极梯度文件，也不能把占位列当成真实谱强度。

需要实际红外强度，可采用已测试的构象红外属性流程。本篇只展示数值Hessian与振动频率。

6. 用ASE读取原始Hessian#

资料包 read.py 将Eh/Bohr²转换为eV/Ų，按原子质量进行质量加权，再对角化：

from ase import units
from ase.vibrations import VibrationsData
vibrations = VibrationsData.from_2d(atoms, hessian * units.Ha / units.Bohr**2)
frequencies = vibrations.get_frequencies()

独立质量加权结果与原生六个内部频率的最大差约0.105cm⁻¹，相对差约3.05×10⁻⁵。原子质量常数及刚体投影处理存在实现差别，不要求每一位一致；原始未投影矩阵的低频刚体残余也不应误当作内部振动。

bash read.sh
cat analysis-read/CREST-NH3-numhess-frequencies.csv

下载频率对照表。自己的新任务完成后，改用 bash read.sh work-new，再运行 bash plot.sh 即可绘图。