CREST：数值梯度与原子受力

本例通过 numgrad 对一个三原子水结构计算解析梯度，并由 CREST 原生有限差分得到数值梯度。用途是检查势能梯度接口、理解梯度和力的符号及单位；它不是结构优化，也不等于更高等级电子结构方法的准确性检验。

1. 下载并提取结果#

下载完整输入、实际输出、CSV 与画图脚本，保存为 example.zip，先读取附带结果：

python3 -m zipfile -e example.zip .
cd crest-gradient-reactor-tools
module load ase/3.29.0
bash read.sh
bash plot.sh

需要重新调用 CREST 时，先在申请到的 CPU 作业环境运行 bash build-private.sh，随后回到算例目录运行 bash submit.sh。模板申请 CPU 56 核整节点，原生工具使用 1 线程。程序独立安装到 ~/.local/ys-crest-tools/bin/crest，固定 CREST 3.0.3 维护源码、GNU Fortran 与 MKL；不替换公共 crest/3.0.2。构建脚本首次需要访问源码仓库。

新计算结果在 work-new/，读取方式是 bash read.sh work-new analysis-read。以下直接调用命令均应在已分配的 CPU 作业中执行。

2. 准备结构和完整输入#

cat inputs/gfn2-numgrad/input.xyz
cat inputs/gfn2-numgrad/input.toml

input="input.xyz"
runtype="numgrad"
threads=1
[calculation]
[[calculation.level]]
method="gfn2"
chrg=0
uhf=0
accuracy=0.0001

XYZ 的长度单位是 Å；chrg=0 表示中性，uhf=0 为闭壳层。accuracy=0.0001 用来减少电子求解噪声对有限差分的干扰。第二组输入只将 method 改成 gfnff，保持同一几何。

在两个独立目录运行：

cd inputs/gfn2-numgrad
~/.local/ys-crest-tools/bin/crest input.toml > native.log 2>&1
cd ../gfnff-numgrad
~/.local/ys-crest-tools/bin/crest input.toml > native.log 2>&1

3. 从日志查什么#

grep -A 5 'Final molecular gradient' native.log
grep -A 4 'Numerical Gradient' native.log
grep -A 4 'Gradient Difference' native.log
grep 'Total number of energy+grad calls' native.log

每组三原子有九个笛卡尔分量，中心差分需要九次正位移和九次负位移，再加原结构一次，实际各 19 次调用。

读取器统一导出 CREST-numgrad-components.csv，每行包含方法、原子号、方向、解析梯度、数值梯度、差值和换算后的力。两种方法的解析和数值梯度分量最大差均约 6×10⁻⁸ Eh/Bohr；输出打印精度会影响最后一位。

4. 梯度和力怎么换算#

梯度是 dE/dx，力是它的负值。这里日志用 Eh/Bohr，换成 eV/Å：

force_eV_A = -gradient_Eh_Bohr * 27.211386245988 / 0.529177210903

不要把正梯度直接写成正受力。也不要用两种方法能量相近来代替逐原子、逐方向的力比较。比较不同方法时还需要统一坐标、电荷、自旋、溶剂和收敛设置。

5. 最终文件#

head -12 analysis-read/CREST-numgrad-components.csv
ls figures/CREST-native-numgrad-components.png

先看全分量图确认差值分布，再根据实际研究问题决定优化或动力学输入。若要优化并计算频率，可继续阅读CREST 谐振动频率。模式语法可参考 CREST 官方源码。