让 openMolcas 输出真正的 Hessian 矩阵

王二葛

# 1 - 背景

最近在用多参考方法，尝试了一众量化软件，先后放弃了 Gaussian、Molpro。听闻 Molcas 是业界多参考的标杆，尝试用开源版 openMolcas 做多参考。

由于下一步的工作需要精确 Hessian，在输入文件中加入了

1. &MCKINLEY
   
2. perturbation
   
3. hessian
   
4. showhessian

复制代码

它的确让程序打印出了许多的 3N*3N 矩阵，可是这些矩阵都与其它量化程序对不上。

openMolcas 做多参考很快，特别是相比 Molpro 惨不忍睹的并行效率。几天后，我又将目光转向了 openMolcas。



# 2 - 解决

经对比，发现其输出的频率、本征向量竟然跟 Gaussian 一模一样（精确到 0.1 厘米波数）。这令人心生激动！如果频率正确，那么 Hessian 矩阵一定也正确。怎么找到正确的 Hessian 矩阵呢？

翻看了 openMolcas 的源码，大致了解了 Hessian 的计算流程。的确如同手册所说，用了某种微扰方法，在微扰之后才有了正确的 Hessian 矩阵，可惜没有关键词可以输出它。

根据输出文件中的文本，定位到了正确的 Hessian 矩阵第一次出现在 src/mclr/output_mclr.f ~530 行，因此只要在调用 FreqAnal 之前输出 Hessian 矩阵即可。

1. write(6,*) 'Real Hessian by wb'
   
2. write(6,'(3E18.10)') Hess

复制代码

与您分享。

Freeman

用vitables打开频率任务正常结束的输出文件***.slapaf.h5，里面有HESSIAN标签页，格式是https://gaussian.com/external/的FILE FORMATS最后一段所描述的那样。
也可以用python批量读取：

1. import h5py
   
2. import numpy as np
   
3. h5=h5py.File('***.slapaf.h5','r')
   
4. np.savetxt('hessian.txt',h5['HESSIAN'][:])
   
5. h5.close()

复制代码