使用KST48进行交叉点的频率计算

Accelerator

交叉点频率的处理一直以来缺乏方便的程序。由于交叉点并非3N-6维坐标下的驻点，直接用Gaussian取其中一个电子态进行频率计算的结果严格来说是没有意义的。为了解决这一问题，可以将它的Hessian矩阵投影到3N-7维子空间去，即如下的过程：

其中H1和H2是Gaussian直接得到的笛卡尔力常数矩阵，P是投影算符，投影掉了6个平动转动自由度和1个g1-g2方向的自由度。g1和g2是两个态的梯度。

KST48的KST48Freq模块（https://github.com/RimoAccelerator/KST48/）可以用于这项工作。目前KST48Freq只支持Gaussian；需要输入2个态下频率计算的fchk文件。运行方法非常简单，只需执行python3 kst48_freq.py，接着会以交互式的方式运行。以下是一个亚甲基卡宾的例子：

1. Freq fchk file for state A?a.fchk
   
2. Freq fchk file for state B?b.fchk
   
3. Cosine for the angle formed by f1 and f2 is 0.9999592481491033
   
4. Frequencies for the HEff of CP:
   
5. [-2.13821824e-05 -2.13821824e-05 -4.41518653e-06  7.24372107e-06
   
6.   7.24372107e-06  1.51873184e-05  2.94043963e-02  2.68187257e+03
   
7.   3.04107649e+03]
   
8. Now it's time to output something.
   
9. Give me a freq log file so that I could replace the vibration information
   
10. and output a new file to kst48_freq.out.
    
11. Note that only the frequencies and modes are replaced,
    
12. and all the other things are remained.
    
13. It can be read by GoodVibes to obtain free energy correction.
    
14. Your freq .log file?5_A_freq.log

复制代码

程序会让你输入2个fchk文件的位置，接着会要求一个频率计算的log文件，它起到模板的作用，你可以使用任何一个态的文件。程序会把处理后的交叉点的频率和模式替换进这个文件中，随后输出到kst48_freq.out里。这个文件可以用GaussView可视化，也可以用GoodVibes读取并得到热力学校正量。Shermo无法读取该文件，可能是由于由于其中输出的空格数量等与Gaussian标准格式有差别。
理想情况下，应该会有3N-7个模式。在KST48Freq的输出过程中，自动过滤掉了波数小于5 cm-1的模式（在屏幕的输出中，Frequencies for the HEff of CP后面的列表包含了全部模式的频率，7个被投影掉的也有，由于数值误差并不会为0。单位是cm-1），因此在少数情况下可能会发现输出的模式数量少一两个。
以亚甲基卡宾为例，打开输出文件，会发现记录了3N-7 = 2个振动模式：

用Goodvibes可以得到热力学校正：

1. o  kst48_freq                                 -38.619793   0.013038    -38.602979   0.022260   0.022260    -38.625238    -38.625238
   

复制代码

当然也可以借助此功能通过虚频数量对是否正确找到了MECP进行验证。另外，由于交叉点必须满足g1与g2共线，当两者不共线时在运行过程中也会进行警告。

sobereva

ORCA有个SurfCrossNumFreq关键词可以在MECP搜索最后做3N-7维空间里的振动分析，和这个方法之间的异同如何？

Accelerator

sobereva 发表于 2024-10-17 10:32
ORCA有个SurfCrossNumFreq关键词可以在MECP搜索最后做3N-7维空间里的振动分析，和这个方法之间的异同如何？

思想是相同的，不过ORCA手册上没有给如何组合HEff的公式细节，不确定实现过程是否完全相同。实测ORCA与KST48Freq的输出结果有微小的差异，如亚甲基卡宾的例子（B3LYP/G STO-3G），ORCA给出的两个频率是2633和3052 cm-1，而振动模式则相同。
ORCA进行此计算需要共计2次数值频率，导致在较大体系没有使用价值，并且只能用于ORCA的MECP模块支持的方法组合（主要是2个不同自旋多重度的情形，以及部分支持CASSCF）。
而KST48Freq只要能拿到两个态的频率计算输出文件就能使用，结合KST48本体高度灵活的交叉点搜索功能，对于激发态之间的交叉、不同透热态的交叉等，都可以方便地得到结果。KST48Freq本身在计算耗时上是免费的，经过实测，对于100多原子的体系也可以瞬间完成，结合Gaussian的解析频率计算，使用wB97xD/def2-SVP水平只需要几个小时就能完成全部计算，而使用ORCA的数值频率简直要不知道算到什么时候。