计算化学公社

标题: 让 openMolcas 输出真正的 Hessian 矩阵 [打印本页]

作者
Author: 王二葛 时间: 2022-12-26 14:09
标题: 让 openMolcas 输出真正的 Hessian 矩阵
# 1 - 背景

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

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

&MCKINLEY
perturbation
hessian
showhessian

复制代码

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

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

# 2 - 解决

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

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

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

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

复制代码

与您分享。

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

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

复制代码

欢迎光临计算化学公社 (http://bbs.keinsci.com/)