Gaussian在频率计算后的驻点测试疑似存在错误

beefly

Gaussian在频率计算后会打印 Eigenvalues （位于Hessian矩阵和梯度收敛检测之间），Eigenvalues 个数通常就是振动频率个数。它们的来源很好理解：Gaussian 用 Newton-Raphson （NR）方法估算下一步步长时，会求解振动久期方程，Eigenvalues就是方程的特征值。示例输出如下：

1.                           Z2        X3        Y3        Z3
   
2.            Z2           0.17587
   
3.            X3          -0.00000  -0.00003
   
4.            Y3          -0.03718  -0.00000   0.31336
   
5.            Z3           0.00604  -0.00000   0.17987   0.17587
   
6. ...
   
7.      Eigenvalues ---    0.16571   0.60781   0.80113
   
8. ...
   
9. Variable       Old X    -DE/DX   Delta X   Delta X   Delta X     New X
   
10.                                  (Linear)    (Quad)   (Total)
    
11.     X1        0.00000  -0.00000   0.00000  -0.00000  -0.00000   0.00000
    
12.     Y1        0.00000  -0.00000   0.00000  -0.00000  -0.00000   0.00000
    
13.     Z1        0.23890  -0.00002   0.00000  -0.00017  -0.00027   0.23862

复制代码

很久以前我就发现，Eigenvalues与Hessian矩阵真正的特征值（也就是原子单位下的振动频率平方）一般不存在对应关系。这些天经过大量测试，我终于找到了其中的原因。以下为水分子测试，坐标是在B3LYP/3-21G级别优化的：

1. O      0.000000000000      0.000000000000    0.126418473599
   
2. H      0.000000000000      0.784764729383   -0.487838361800
   
3. H      0.000000000000     -0.784764729383   -0.487838361800

复制代码

首先排除原子量的影响。将氢、氧的原子量都设为1.0（坐标部分改为O(iso=1.0)、H(iso=1.0)），然后计算频率。这时候Eigenvalues与频率平方是一样的，证明二者之间确实存在对应关系。

接下来改回默认原子量，发现二者就完全不一样了。因为结构是充分优化过的，可以排除振-转模式是否进行投影造成的影响。很多计算都显示Eigenvalues总是更大一些，这说明对于水分子，采用的O、H“原子量”可能是小一点的数。首先想到的可能是核电荷数。把O的原子量改为8（O(iso=8.0)）试试。计算后，发现Eigenvalues三个数中最大的值0.80113与B2振动模式的频率平方一致，然而前两个数还是对不上。即便如此，这也足以证明Gaussian的 NR 方法估算步长时，用了核电荷数而非原子量。

为什么要用核电荷数而不是原子量呢？这大概是为了避免结构优化受同位素影响。例如重水HDO，用原子量估算步长会在结构优化中导致H-O、D-O键长存在轻微差别，破坏波函数的C2v对称性。由于大部分量子化学计算基于B-O近似（电子波函数和原子框架可以有不一样的点群对称性），对称破缺的结构不利于分析。

怎么得到前两个数呢？水分子有3个振动模式：2个A1模式，1个B2，其中B2和A1模式是永远正交的，而两个A1模式需要对角化才能保证正交。经过很多失败的尝试后，我做了一个最坏的设想：Gaussian用错了公式，导致两个A1模式的Eigenvalues不正交！于是我用matlab写了下面的测试程序进行验证。

1. % H2O
   
2. %  O      0.000000000000     -0.000000000000    0.126418473599
   
3. %  H     -0.000000000000      0.784764729383   -0.487838361800
   
4. %  H      0.000000000000     -0.784764729383   -0.487838361800
   
5. %
   
6. Natm=3;
   
7. % atomic masses
   
8. m=[8 1 1];
   
9. % force constant matrix computed at the B3LYP/3-21G level of theory
   
10. f=[
    
11.   1.94853383E-05 -9.30582821E-16 -8.22729703E-16 -9.74266831E-06 -5.04562860E-17  4.85890549E-17 -9.74266830E-06  8.24775901E-16  6.13401158E-16
    
12. -9.30582821E-16  5.54616379E-01  3.72245637E-14 -2.35078770E-17 -2.77308189E-01  2.17053899E-01  8.46266535E-16 -2.77308189E-01 -2.17053899E-01
    
13. -8.22729703E-16  3.72245637E-14  3.63828365E-01 -3.97504322E-17  1.42689337E-01 -1.81914183E-01  7.37686754E-16 -1.42689337E-01 -1.81914183E-01
    
14. -9.74266831E-06 -2.35078770E-17 -3.97504322E-17 -2.82354345E-05  2.59983521E-18  4.72518175E-18  3.79781028E-05  5.71487206E-17  6.55811015E-17
    
15. -5.04562860E-17 -2.77308189E-01  1.42689337E-01  2.59983521E-18  3.13360609E-01 -1.79871618E-01  3.66031172E-17 -3.60524196E-02  3.71822811E-02
    
16.   4.85890549E-17  2.17053899E-01 -1.81914183E-01  4.72518175E-18 -1.79871618E-01  1.75870959E-01 -4.07033663E-17 -3.71822811E-02  6.04322314E-03
    
17. -9.74266830E-06  8.46266535E-16  7.37686754E-16  3.79781028E-05  3.66031172E-17 -4.07033663E-17 -2.82354345E-05 -7.62847125E-16 -5.66799750E-16
    
18.   8.24775901E-16 -2.77308189E-01 -1.42689337E-01  5.71487206E-17 -3.60524196E-02 -3.71822811E-02 -7.62847125E-16  3.13360609E-01  1.79871618E-01
    
19.   6.13401158E-16 -2.17053899E-01 -1.81914183E-01  6.55811015E-17  3.71822811E-02  6.04322314E-03 -5.66799750E-16  1.79871618E-01  1.75870959E-01
    
20. ];
    
21. mm=zeros(Natm*3,1);
    
22. for i=1:Natm
    
23.   for j=3*i-2 : 3*i
    
24.     mm(j)=m(i);
    
25.   end
    
26. end
    
27. ma=1./sqrt(mm);
    
28. % mass-weighted force constant matrix
    
29. ff=diag(ma)*f*diag(ma);
    
30. [p,q]=eig(ff);q=diag(q);
    
31. % mass-weighted vib. modes
    
32. vec=p(:,7:9);
    
33. 
    
34. % check results
    
35. diag(vec'*ff*vec)        %  0.11760, 0.45920, 0.51651
    
36. 
    
37. eig(vec'*f*vec)          %  0.16571, 0.60781, 0.80115 (G16: 0.16571   0.60781   0.80113)

复制代码

对公式的解释：
振动久期方程为
F L = M L E
其中F是直角坐标系下的Hessian矩阵（力常数矩阵），M是对角的原子量矩阵（为了重复Gaussian结果用了核电荷数），L、E分别是特征矢量和特征值（也就是频率平方和Eigenvalues）。
以上方程是赝-特征值方程，在实际计算中人们习惯把它写成标准特征值方程的形式：
FF LL = LL E
其中，FF = M^(-1/2) F M^(-1/2)，称为质量加权的Hessian矩阵，LL = M^(1/2) L。根据特征值方程的特性，有如下关系：
LL^T FF LL = E
但是以下关系一般不成立（除非原子量都是1）：
LL^T F LL = E

在最后一行代码中，由于两个A1模式不正交，导致E矩阵并不是对角的，因此需要额外执行一次对角化。

matlab测试程序的结果表明，Gaussian 确实是用了最后一个不成立的公式（或者等价的算法）来计算 Eigenvalues ，再通过 NR 方法估算步长。由于 F 矩阵是未经过质量加权的（没有除以原子量矩阵），导致 Eigenvalues 偏大，进而把 NR 方法估算的 Displacement 也降低了。这可能带来的后果是（纯猜测，未经验证）：即便优化质量不高的几何结构，也可能会被误判为驻点结构，含有重元素的分子尤甚。

大家都不喜欢出现 No 的情况，所以这个错误不改也挺好。