用内坐标做结构优化和频率，可能会得到错误结果

beefly

由于内坐标使用不当，在结构优化和振动频率计算中可能会隐藏着错误，对研究构成误导。以下是最近发现的一个文献中的错误，并且被以讹传讹地传播了数十年。

PF5分子的“旋转门”过渡态。1975年，Russegger和Brickmann提出PF5存在“旋转门”（Turnstile）转动模式，PF5的五个氟原子分成2+3两组，分别位于磷原子上下的两个平面上。通过两个平面的相对转动，PF5发生分子重排反应，氟原子重新排布。“旋转门”过程的过渡态具有Cs对称性。参见Russegger and Brickmann, Chem. Phys. Lett., 30, 276, 1975。以下是PF5旋转门模式的示意图，来自C. D. Montgomery, J. Chem. Edu. 78(6), 844, (2001)。注意：图中有个错误，四个键角APD, BPC, BPE, CPE被固定在了90度。按照原始定义，APD与BPC=BPE=CPE应该是被优化的。



这个猜测一直被学术界接受，只不过认为此过渡态的能量远远高于Berry赝旋转的过渡态（对称性C4v），因此实验上观测不到。果真如此吗？

以下是“旋转门”过渡态优化的输入文件。结构参数已经优化过了，因此很快就能收敛。这个过渡态的结构优化无法用直角坐标实现，因此必须用Z矩阵定义结构。优化过程也使用Z矩阵，因为无论是直角坐标优化还是Gaussian默认的冗余内坐标优化，都会破坏结构参数的比例关系，导致能量更低的C4v过渡态。

1. %mem=4GB
   
2. %NProcShared=4
   
3. %chk=pf5ts
   
4. # b3lyp/cc-pvtz int(grid=ultraFine) opt(z-mat,vtight)
   
5. 
   
6. PF5 Cs TS of turnstile pseudorotation
   
7. 
   
8. 0 1
   
9. P
   
10. X 1 1.0
    
11. F 1 r1 2 a1
    
12. F 1 r1 2 a1 3 120.
    
13. F 1 r1 2 a1 4 120.
    
14. F 1 r2 2 a2 3  90.
    
15. F 1 r2 2 a2 3 -90.
    
16. 
    
17. r1          1.571170
    
18. r2          1.595396
    
19. a1        121.7419
    
20. a2         42.1664

复制代码

优化结束后，接下来计算振动频率，验证上面优化的结构是否是过渡态。一种方法是使用优化好的Z矩阵参数，做一个独立的振动频率计算。当然也可以把opt和freq合并为一步。如我们所预期的，有一个虚频（i*171 cm-1, A"模式）。顺便检查频率计算之后的梯度，仍然是收敛的：

1.                   Item               Value     Threshold  Converged?
   
2. Maximum Force            0.000153     0.000450     YES  
   
3. RMS     Force            0.000098     0.000300     YES
   
4. Maximum Displacement     0.000158     0.001800     YES
   
5. RMS     Displacement     0.000109     0.001200     YES
   
6. Predicted change in Energy=-1.860418D-08

复制代码

如此看来，如果忽略DFT不适用的个别情况，“旋转门”过渡态看起来应该是真实存在的。

但是，如果我们用优化后的直角坐标（从checkpoint文件读取）做频率计算，就会得到不同结果。

1. %mem=6GB
   
2. %NProcShared=4
   
3. %chk=pf5ts
   
4. # b3lyp/cc-pvtz int(grid=ultrafine) geom=checkpoint freq
   
5. 
   
6. PF5 Cs TS of turnstile pseudorotation
   
7. 
   
8. 0 1

复制代码

通过比较发现，电子总能量，振动频率，各种热化学量，和前一个频率计算的结果都是一样的。唯一的区别是频率计算之后的梯度部分，显示直角坐标系下的梯度并未收敛，表明这并不是过渡态。

1.                   Item               Value     Threshold  Converged?
   
2. Maximum Force            0.034138     0.000450     NO
   
3. RMS     Force            0.013785     0.000300     NO
   
4. Maximum Displacement     0.192385     0.001800     NO
   
5. RMS     Displacement     0.070820     0.001200     NO
   
6. Predicted change in Energy=-1.315610D-02

复制代码

为什么会有截然相反的结果呢？

与分子对称性密切相关的是分子自由度，也就是分子独立坐标的个数。搜索Gaussian输出文件中的“Deg. of freedom”，发现数值是7，也就是说，要优化具有Cs对称性的“旋转门”过渡态，必须要用7个独立的内坐标变量来描述。7个之外的内坐标参数，要么是固定的常数，要么存在线性依赖关系（比如二面角-B1依赖于B1）。而上面的Z矩阵坐标只用到4个内坐标参数，这就会导致严重问题：当把直角坐标的梯度转化成内坐标梯度的时候，一些重要的信息丢失了。

在用内坐标做优化或频率计算的时候，如果注意以下两点，就可以避免或及时发现上面的问题。

1，用Z矩阵做OPT+FREQ或者FREQ计算的时候，一定要检查独立坐标变量个数，是否等于（实际上也可以大于）分子自由度。

2，全部坐标参数都写成坐标符号的形式，避免用数值。例如上面的坐标，要写成：

1. P
   
2. X 1 1.0
   
3. F 1 r1 2 a1
   
4. F 1 r1 2 a1 3  b1
   
5. F 1 r1 2 a1 4  b1
   
6. F 1 r2 2 a2 3  b2
   
7. F 1 r2 2 a2 3 -b2
   
8.     Variables:
   
9. r1          1.571170
   
10. r2          1.595396
    
11. a1        121.7419
    
12. a2         42.1664
    
13.     Constants:
    
14. b1        120.0
    
15. b2         90.0

复制代码

在OPT或者FREQ计算之后，不要忘记检查常数坐标参数的梯度。下面的信息显示b1和b2的梯度很大，显然这并不是过渡态。

1.                        ----------------------------
   
2.                        !   Optimized Parameters   !
   
3.                        ! (Angstroms and Degrees)  !
   
4. ----------------------                            ----------------------
   
5. !      Name          Value   Derivative information (Atomic Units)     !
   
6. ------------------------------------------------------------------------
   
7. !       r1          1.5711   -DE/DX =    0.0002                        !
   
8. !       r2          1.5954   -DE/DX =   -0.0001                        !
   
9. !       a1        121.7385   -DE/DX =   -0.0001                        !
   
10. !       a2         42.1615   -DE/DX =    0.0                           !
    
11. !       b1        120.0      -DE/DX =   -0.0509                        !
    
12. !       b2         90.0      -DE/DX =   -0.0742                        !
    
13. ------------------------------------------------------------------------
    

复制代码

ChemiAndy

谢谢，很有意思，就是没看太明白。配上图会好一些。

1.“这个过渡态结构无法用直角坐标定义，因此必须用Z矩阵。” 这句话是不是不妥？
2. 第一次优化没有把b1, b2列入变量列表，这对于opt=Zmat应该算是一个明显的错误吧？文献犯的是这个错误吗？

smutao

ChemiAndy 发表于 2015-7-12 02:05
谢谢，很有意思，就是没看太明白。配上图会好一些。

1.“这个过渡态结构无法用直角坐标定义，因此必须用 ...

beefly

ChemiAndy 发表于 2015-7-12 02:05
谢谢，很有意思，就是没看太明白。配上图会好一些。

1.“这个过渡态结构无法用直角坐标定义，因此必须用 ...

因为文章太老，电子版的分辨率太低，无法提供原图。只能用一个后来发表的图，但是存在错误。可以当作示意图来看。

1.是表述错误。应该是过渡态结构的优化。已修改。

2.Gaussian程序允许把固定的结构参数直接写在坐标中。文献用的是Pople的半经验程序，可能是Gaussian的前身，具体怎么写的就不知道了，只能从文章正文叙述以及图、表格中得知有四个坐标变量。这个错误依赖程序，幸好Gaussian会打印固定参数b1、b2的梯度。如果用其他程序（比如CFour，Molpro），这个错误可能永远都发现不了。

ChemiAndy

beefly 发表于 2015-7-13 06:03
因为文章太老，电子版的分辨率太低，无法提供原图。只能用一个后来发表的图，但是存在错误。可以当作示意 ...

Much better, thanks.

smutao

ChemiAndy 发表于 2015-7-13 10:36
Much better, thanks.

ChemiAndy 你知道原因吗？
其实很简单，Z-matrix最多提供3N-6个结构参数，当指定了opt=z-matrix的时候，B矩阵中有且仅有这3N-6个结构参数，当G'=BB+矩阵做对角化时，我们发现居然有为0的特征值出现了！这说明Z-matrix的3N-6个参数中有冗余性！这个时候，后面就不对了，更别说后续的把3N维度下的力分解到3N-6个维度去了。
而高斯默认的冗余内坐标会提供尽可能多的结构参数，当G'对角化时，我们明显可以看到有3N-6个非0特征值，其余都是0，这个时候力分解过去是对的（虽然这个过程中涉及到B矩阵的pseudo-inverse，不同的人可能有不同的定义）。

smutao

详细的公式可以参考 J. Chem. Phys., Vol. 96, No. 4, 1992
The Journal of Chemical Physics 96, 2856 (1992); doi: 10.1063/1.462844