O1NumHess：只用常数个梯度计算数值Hessian的程序

wzkchem5

一、简介

众所周知，计算含N个原子的体系的数值Hessian，利用双向差分的方法，需要计算6N个梯度，即使采用较粗糙的单向差分方法，也需要计算3N+1个梯度。因此数值Hessian的耗时往往很高，对于不支持解析Hessian的方法（如BDF、ORCA的TDDFT Hessian计算），Hessian计算往往会成为整个计算流程的瓶颈。

最近我们发现，Hessian矩阵的离对角线较远的块呈现低秩性，因此虽然Hessian有O(N^2)个矩阵元，但只需要O(N)个参数即可准确描述任意体系的Hessian。因此，只需要计算O(1)个梯度，即可得到足以确定Hessian所需的信息（因为每个梯度包含O(N)个已知数），由此可以大大减少计算数值Hessian计算所需的梯度数目。受此启发，我们开发了O1NumHess方法，对于较小的体系相比传统双向差分数值Hessian算法节省一半的梯度计算，对于足够大的体系，所需梯度数目稳定在100个左右；频率/ZPE/焓/Gibbs自由能误差仅约为传统数值Hessian的2倍，如振动频率误差约3 cm-1，Gibbs自由能误差约1 kcal/mol；对于基态计算，内存消耗仅为解析Hessian的几十分之一，且计算速度往往较BDF的解析Hessian还要更快。因此可以预期O1NumHess主要适用于以下场景：

• TDDFT数值频率计算，以及其他不支持解析Hessian的方法的数值频率计算（对于100原子体系加速5~8倍）；
• 50个原子以上、对于自由能精度要求不高（允许有1~2 kcal/mol误差）情况下的频率计算。虽然BDF支持DFT的基态解析频率，但O1NumHess对于50个原子以上的体系，一般较解析频率更快。需要注意的是随着将来BDF解析频率代码的优化，解析频率与O1NumHess算法的crossover point可能会出现得更晚，也就是对于更大的体系O1NumHess才能表现出相比解析频率的速度优势；
• 特别大的体系（如150原子以上）的频率计算，此时解析频率因所需内存过多，需要分批求解耦合微扰方程，导致计算显著变慢，使得O1NumHess具有明显优势；
• 并行核数较多（如32个核以上）的频率计算，O1NumHess继承传统数值Hessian算法的优势，具有很高的并行效率，并行核数越多，较解析Hessian的优势越明显。
  

O1NumHess的理论与测试结果已上传至arXiv（https://arxiv.org/pdf/2508.07544），对应的Python代码已于GitHub上公开（https://github.com/ilcpm/O1NumHess、https://github.com/ilcpm/O1NumHess_QC），已实现BDF和ORCA程序的接口，可以直接按程序文档（https://github.com/ilcpm/O1NumHess_QC/blob/main/doc/doc.md）安装使用。其中O1NumHess库实现了底层的数学算法，可以用来求解任意多元函数的Hessian，而不限于量子化学里的Hessian（电子能对核坐标的Hessian）；O1NumHess_QC库是O1NumHess与BDF、ORCA的接口，实现了读取xyz文件，并调用BDF、ORCA计算梯度等功能。

二、使用方法

O1NumHess的使用方法已于文档https://github.com/ilcpm/O1NumHess_QC/blob/main/doc/doc.md中说明，此处仅进行大致的复述。

1. 安装（当使用2025年7月及以后的BDF版本时，可跳过这一步）
O1NumHess需要在Linux系统（包括WSL等虚拟Linux环境）下运行。
（1）首先确认安装了Python 3.6以上版本，方法是在命令行下输入python --version，若报错或显示的版本号低于3.6（如显示Python 2.7.0、Python 3.5.2等），则需要安装Python 3，方法参见https://wiki.python.org/moin/BeginnersGuideChinese。这一步如果遇到问题可以直接询问任意大语言模型，得到的回答一般是可靠的。
（2）下载O1NumHess库。对于不熟悉Git的用户，最简单的方式是：单击https://github.com/ilcpm/O1NumHess的绿色Code按钮，选择Download ZIP，在任意不含中文字符/空格/特殊符号的路径下解压即可。熟悉Git的用户也可用Git clone的方式下载。
（3）安装O1NumHess库。在命令行窗口中，cd到O1NumHess的根目录下，然后运行python3 setup.py install --prefix ~/.local。若当前计算机没有安装NumPy、SciPy，则程序会自动联网下载NumPy、SciPy。若因没有联网或网速过慢导致无法下载，则应按https://numpy.org/install/、https://scipy.org/install/的方法自行下载安装，再重新安装O1NumHess库。这一步如果遇到问题可以直接询问任意大语言模型，得到的回答一般是可靠的。
注意：虽然可以用Anaconda一步到位安装Python 3、NumPy和SciPy，但据我们测试，某些版本的Anaconda结合BDF进行并行的O1NumHess计算时，会因个别梯度计算失败而导致Hessian计算失败。因此除非确实有必要，否则不建议使用Anaconda。
（4）下载O1NumHess_QC库。方法与（2）相同，但链接改为https://github.com/ilcpm/O1NumHess_QC。
（5）安装O1NumHess_QC库。方法与（3）相同，但改为在O1NumHess_QC的根目录下执行操作。
（6）按https://github.com/ilcpm/O1NumHess_QC/blob/main/doc/doc.md#config的说明，将~/.O1NumHess_QC/BDF_config_example.py复制为~/.O1NumHess_QC/BDF_config.py（当打算使用BDF进行梯度计算时），或将~/.O1NumHess_QC/ORCA_config_example.py复制为~/.O1NumHess_QC/ORCA_config.py（当打算使用ORCA进行梯度计算时），并根据当前计算机的BDF、ORCA、OpenMPI等的安装路径，修改复制的.py文件中的BDFHOME、LD_LIBRARY_PATH、PATH、MPI_HOME等变量，以及最后的小写的path变量。

2. 基于2025年7月之前的BDF版本进行O1NumHess计算
首先创建一个包含待计算结构的xyz文件，例如job.xyz。然后创建一个计算梯度的BDF输入文件，注意其中原子坐标必须用"file=job.xyz"的语法从job.xyz读取，而不能将原子坐标直接写入job.inp中（参见https://bdf-manual.readthedocs.i ... d%20Output.html#id5）。运行Python 3，输入以下命令（或将以下命令保存为.py文件并用Python 3运行）：
from O1NumHess_QC import O1NumHess_QC

qc = O1NumHess_QC("path/to/job.xyz")
hessian = qc.calcHessian_BDF(
    method = "o1numhess",
    delta = 5e-3,
    core = 4,
    total_cores = 16,
    mem = "4G",
    inp = "path/to/job.inp",
)
其中delta为以Bohr为单位的数值差分步长，core为每个梯度使用的核数，total_cores为所有梯度计算一共能使用的核数，mem为每个核所能使用的内存，"path/to/"为job.inp、job.xyz所在的绝对路径。具体说明参见https://github.com/ilcpm/O1NumHess_QC/blob/main/doc/doc.md#bdf，这里仅指出一些一般的规则：

• 若忽略total_cores，则程序用当前节点的所有核进行计算。若使用排队系统，或当前节点有其他任务在运行，则必须指定total_cores。
• core太小会导致并行负载不均衡，太大会导致梯度计算的并行效率低，两者都会导致计算时间变长。core的最优取值可以按total_cores/min(20,0.5*N)估算，其中N为体系的原子数；取离计算结果较近且能整除total_cores的一个非零整数，作为core的取值。
  

计算结束后，程序打印Hessian矩阵，数据类型为np.ndarray。可自行撰写Python脚本将hessian变量以任意格式输出到文件，或进行后续数学处理。

3. 基于ORCA进行O1NumHess计算
首先创建一个包含待计算结构的xyz文件，例如job.xyz。然后创建一个计算梯度的ORCA输入文件（一般需在感叹号开头的行中添加EnGrad关键词），注意其中原子坐标必须用"* xyzfile xxx xxx job.xyz"的语法从job.xyz读取，而不能将原子坐标直接写入job.inp中（参见https://www.faccts.de/docs/orca/ ... from-external-files）。运行Python 3，输入以下命令（或将以下命令保存为.py文件并用Python 3运行）：
from O1NumHess_QC import O1NumHess_QC

qc = O1NumHess_QC("path/to/job.xyz")
hessian = qc.calcHessian_ORCA(
    method = "o1numhess",
    delta = 5e-3,
    total_cores = 16,
    inp = "path/to/job.inp",
)
各参数的意义参见上一小节。注意对于ORCA计算，每个梯度计算所需的核数不应用core关键词指定，而应该用%pal nprocs xxx end的语法直接写入ORCA输入文件中；每个核所需的内存不应用mem关键词指定，而应该用%maxcore xxx的语法直接写入ORCA输入文件中。计算结束后，程序打印Hessian矩阵，可自行以任意格式输出或进行任意后续处理。

4. 基于2025年7月及之后的BDF版本进行O1NumHess计算
自2025年7月起，BDF自带O1NumHess、O1NumHess_QC库，无需另行安装O1NumHess、O1NumHess_QC库即可使用，且用户无需Python相关知识。只需在一个BDF的Hessian输入文件中加入O1NumHess关键词，即可自动调用O1NumHess进行数值Hessian计算。更多细节参见https://bdf-manual.readthedocs.i ... .html#o1numhess-o-1。目前（2025年8月），自带O1NumHess的BDF版本还没有公开发行，但对于在鸿之微云上使用BDF的用户，可在“内蒙A区”超算上使用2025年7月初的BDF版本，计算结果可能与上述arXiv文章的结果略有区别。相关运行问题请咨询鸿之微客服。


三、未来计划

未来将进一步改善O1NumHess的程序生态，包括：

• 实现将ORCA计算的Hessian导出成.hess格式的功能，从而可以被ORCA读取，进行热化学分析、ESD等计算（但实现Herzberg-Teller效应的计算仍需较多工作）；
• 实现断点续算功能；
• 实现其他软件的接口，如Gaussian等（我们缺乏人手，欢迎有兴趣者贡献接口代码）；
• 改善对Anaconda的支持；
• 实现红外、拉曼光谱计算。
  

此外，O1NumHess的arXiv文章（https://arxiv.org/pdf/2508.07544）将于近几天投稿，欢迎大家提出宝贵的意见和建议！若意见或建议被采纳，我们将在文章中予以诚挚致谢。