Inv3 failed in PCMMkU (L502/L508 报错) 的成因与解决方法

liyuanhe211

近日遇到一个比较罕见的报错，形式如下

1. (Enter /home/g16/l502.exe)
   
2. Integral symmetry usage will be decided dynamically.
   
3. Closed shell SCF:
   
4. Using DIIS extrapolation, IDIIS=  1040.
   
5. NGot= 18983288832 LenX= 18982597916 LenY= 18982263254
   
6. Requested convergence on RMS density matrix=1.00D-08 within 128 cycles.
   
7. Requested convergence on MAX density matrix=1.00D-06.
   
8. Requested convergence on             energy=1.00D-06.
   
9. No special actions if energy rises.
   
10. Fock matrices will be formed incrementally for  20 cycles.
    
11. 
    
12. Cycle   1  Pass 1  IDiag  1:
    
13. FoFJK:  IHMeth= 1 ICntrl=       0 DoSepK=F KAlg= 0 I1Cent=           0 FoldK=F
    
14. IRaf= 990000000 NMat=       1 IRICut=       1 DoRegI=T DoRafI=F ISym2E= 0 IDoP0=0 IntGTp=1.
    
15. FoFCou: FMM=T IPFlag=           0 FMFlag=      100000 FMFlg1=        2001
    
16.          NFxFlg=           0 DoJE=F BraDBF=F KetDBF=F FulRan=T
    
17.          wScrn=  0.000000 ICntrl=         0 IOpCl=  0 I1Cent=           0 NGrid=           0
    
18.          NMat0=    1 NMatS0=      1 NMatT0=    0 NMatD0=    1 NMtDS0=    0 NMtDT0=    0
    
19. Symmetry not used in FoFCou.
    
20. FMM levels:  10  Number of levels for PrismC:   9
    
21. Inv3:  Mode=1 IEnd=    55186563.
    
22. Iteration    1 A*A^-1 deviation from unit magnitude is 2.18D-11 for   1454.
    
23. Iteration    1 A*A^-1 deviation from orthogonality  is 1.36D-11 for   4245   1454.
    
24. Iteration    1 A^-1*A deviation from unit magnitude is 2.55D-11 for   1454.
    
25. Iteration    1 A^-1*A deviation from orthogonality  is 2.25D-10 for   4227   3067.
    
26. Iteration    2 A*A^-1 deviation from unit magnitude is 7.03D-07 for   4227.
    
27. Iteration    2 A*A^-1 deviation from orthogonality  is 8.89D-07 for   2879   1454.
    
28. Iteration    2 A^-1*A deviation from unit magnitude is 1.00D-06 for   1454.
    
29. Iteration    2 A^-1*A deviation from orthogonality  is 9.53D-07 for   4227   1454.
    
30. ...
    
31. ...
    
32. Iteration   10 A*A^-1 deviation from unit magnitude is 7.97D-07 for   2943.
    
33. Iteration   10 A*A^-1 deviation from orthogonality  is 1.08D-06 for   2943   1454.
    
34. Iteration   10 A^-1*A deviation from unit magnitude is 3.30D-07 for   4227.
    
35. Iteration   10 A^-1*A deviation from orthogonality  is 6.08D-07 for   4227   2061.
    
36. Inv3 failed in PCMMkU.
    
37. Error termination via Lnk1e in /home/g16/l502.exe at Mon Sep  4 11:33:20 2017.
    
38. Job cpu time:       0 days  0 hours 18 minutes 30.5 seconds.
    
39. Elapsed time:       0 days  0 hours  0 minutes 42.9 seconds.

复制代码

这一报错自 G09 D.01 至发稿时最新的 G16 A.03 均存在，发生于SCF的第一次迭代（含正常的 L502 和 QC 的 L508）。

这一问题未能在网络上找到靠谱的解决方法，调整各种常用的解决SCF不收敛问题的选项没有效果（可以理解，因为它不是SCF的问题）。故向 Gaussian Help 询问，回复很快，大意如下（译自英文）：

这个问题与溶剂模型的溶质空腔有关。高斯中为构建SMD溶质空腔的默认选项是为中等大小的溶质分子设置的（注：默认是范德华表面 surface=vdW），它是由以原子为中心、一定的原子半径的一些球构成的。其中用于构建空腔的默认原子半径较小，对较大的溶质分子会造成溶剂空腔长得像“奶酪”，有很多数学形式上暴露于溶剂中的孔洞、沟槽，但在实际体系中这些位置并不能接触溶剂。这会造成上述为 PCM 矩阵求逆迭代过程中不收敛的问题，也是你另一封邮件中提到SCF能量异常的本质原因（我同时发送了两封邮件，下文将提到）。

理想情况下，使用溶剂排斥表面（Solvent-excluding surface, SES）是更好的选择。SES 首先以原子位置为中心，构建出范德华表面；之后再用一些原点不在原子上的球面进行平滑，就消除了上述不能被溶剂够到的孔洞、沟槽，解决此问题。但其劣势在于这个方法做出的表面不随构象的改变而连续变化，造成势能面不连续、不平滑，难以进行几何优化。

另一个可能是用溶剂可及表面（Solvent-accessible surface，SAS）构建空腔。溶剂可及表面和 vdW 表面一样都是由以原子为中心的球构成的，但其半径是[原子半径+溶剂半径]（注，这个时候溶剂半径才有用，默认选项下自定义溶剂时设半径没什么用）。SAS 表面没有“奶酪问题”，且势能面连续。但其问题是半径太大，空腔扣的太多，相比 vdW 和 SES 表面会低估溶剂化效应。

一个较好的解决方法是用 SAS 表面做优化，然后用 SES 表面做单点。或者在 SAS 表面做优化后再用其结果为初猜做 SES 表面下的单点（以期待 SES 在极小点附近可能没有不连续性）。这两种做法都比在气相优化的结果好一些。

附：Sobereva（2L）提供的三种溶质空腔示意图：

关键词：

更换溶剂表面的类型用 scrf=read 关键词启用溶剂化控制的段落，并在相应段落声明 surface=xxx，如 surface=SAS。
示例如下：

1. #p pbe1pbe/genecp opt freq
   
2. empiricaldispersion=gd3bj
   
3. scrf=(smd,solvent=water,read)
   
4. 
   
5. title
   
6. 
   
7. 0 1
   
8. C                   1.629576   -0.586781    1.988383
   
9. ....
   
10. H                   3.604870   -0.932068   -2.049232
    
11. 
    
12. [Basis set definition]
    
13. ...
    
14. 
    
15. [Pseudopotential definition]
    
16. ...
    
17. 
    
18. surface=sas
    

复制代码

对 surface=SES，还需添加 AddSph 选项。即将上述文件中的 “surface=sas” 换为 “surface=SES AddSph”

如果需要，可以用 geom=allcheck 等语句接续下一步计算切换为其他表面做单点/优化。
--------------------------------------------------------------------------------------------------------------------------------------------------------------------

• 关于“另一封邮件”：
  

这个问题也有可能不报错，但造成 SCF 能量明显异常。这个错误比较隐蔽，当出现能量明显不对的时候应该检查。

例如我在一次优化中发现有两次 SCF 迭代后的能量（用>>>标出）低了好几十个 Hartree，明显是不对的，SCF 迭代所需的次数也明显较多。这也是上面所述的溶剂空腔问题。

1.         Line 68580:  SCF Done:  E(RB-P86) =  -3825.13599635     A.U. after   13 cycles
   
2.         Line 69843:  SCF Done:  E(RB-P86) =  -3825.13662431     A.U. after   12 cycles
   
3.         Line 71105:  SCF Done:  E(RB-P86) =  -3825.13636208     A.U. after   12 cycles
   
4. >>> Line 73188:  SCF Done:  E(RB-P86) =  -3845.72449495     A.U. after   54 cycles
   
5.         Line 74770:  SCF Done:  E(RB-P86) =  -3825.13665330     A.U. after   20 cycles
   
6.         Line 76261:  SCF Done:  E(RB-P86) =  -3825.13647203     A.U. after   23 cycles
   
7.         Line 77804:  SCF Done:  E(RB-P86) =  -3825.13594059     A.U. after   26 cycles
   
8.         Line 79343:  SCF Done:  E(RB-P86) =  -3825.13482223     A.U. after   26 cycles
   
9.         Line 81098:  SCF Done:  E(RB-P86) =  -3825.13408631     A.U. after   35 cycles
   
10. >>> Line 83672:  SCF Done:  E(RB-P86) =  -3888.91182114     A.U. after   73 cycles
    
11.         Line 85289:  SCF Done:  E(RB-P86) =  -3825.13665328     A.U. after   22 cycles
    
12.         Line 86819:  SCF Done:  E(RB-P86) =  -3825.13647358     A.U. after   25 cycles
    
13.         Line 88367:  SCF Done:  E(RB-P86) =  -3825.13595614     A.U. after   26 cycles
    
14.         Line 89942:  SCF Done:  E(RB-P86) =  -3825.13490799     A.U. after   27 cycles

复制代码

sobereva

附上一张图便于读者理解

kyuu

sobereva 发表于 2017-9-6 05:38
附上一张图便于读者理解

图里还可以区分一下范德华球和溶剂探针，毕竟CS和SAS是探针滚出来的

sobereva

kyuu 发表于 2017-9-6 06:01
图里还可以区分一下范德华球和溶剂探针，毕竟CS和SAS是探针滚出来的

稍有判断能力的一看就知道怎么回事，没必要注明

sobereva

貌似surface=SES没用，依然显示的是默认的Cavity type          : VdW (van der Waals Surface) (Alpha=1.000).，结果和默认时相同
surface=SAS会影响结果，显示 Cavity type          : SAS (Solvent Accessible Surface) (Alpha=1.000).
怀疑是bug。g09 E01和A16 A03都存在。可以向官方反映下。

liyuanhe211

sobereva 发表于 2017-9-6 06:44
貌似surface=SES没用，依然显示的是默认的Cavity type          : VdW (van der Waals Surface) (Alpha=1.0 ...

向高斯help反应了 SES/vdW 的问题，和 G16 的 scrf=G09Defaults 的问题

liyuanhe211

sobereva 发表于 2017-9-6 06:44
貌似surface=SES没用，依然显示的是默认的Cavity type          : VdW (van der Waals Surface) (Alpha=1.0 ...

Gaussian Help 回复 Surface=SES 还需要加上 AddSph 选项程序才会用 SES 空腔。经测确实如此。

sobereva

liyuanhe211 发表于 2017-9-7 14:59
Gaussian Help 回复 Surface=SES 还需要加上 AddSph 选项程序才会用 SES 空腔。经测确实如此。

噫，这手册真是不明不白

liyuanhe211

sobereva 发表于 2017-9-8 03:18
噫，这手册真是不明不白

手册都没说这词→_→溶剂化整个一部分的手册都不像高斯写的似得。

lei234

老师您好，以下是我理解该问题的解决办法，请您过目，如有不对，麻烦予以改正，在此谢过！
1）进行sas 优化，关键词 scrf=(smd,solvent=water,read) ，在最后加上surface=sas；
2）SES单点计算，关键词同1），最后加上Surface=SES  AddSph。

赵云跳槽

lei234 发表于 2018-3-31 03:24
老师您好，以下是我理解该问题的解决办法，请您过目，如有不对，麻烦予以改正，在此谢过！
1）进行sas 优 ...

SAS计算，关键词 scrf=read，在最后加上surface=sas；
SES计算，关键词 scrf=read，最后加上addsph, rmin=0.200 ofac=0.890
（以上溶剂是水）


10.1021/acs.jpcb.7b06693

lei234

赵云跳槽 发表于 2018-3-31 10:50
SAS计算，关键词 scrf=read，在最后加上surface=sas；
SES计算，关键词 scrf=read，最后加上addsph, rmi ...

谢谢！
这两个参数，rmin和ofac，是对水而言是吗？那溶剂是别的怎么办呢？这段话的意思是不是溶剂是水的时候可以不加这两个？

赵云跳槽

lei234 发表于 2018-4-1 01:59
谢谢！
这两个参数，rmin和ofac，是对水而言是吗？那溶剂是别的怎么办呢？这段话的意思是不是溶剂是水的 ...

那就不清楚了，只能自己尝试了

leebo

sobereva 发表于 2017-9-6 06:44
貌似surface=SES没用，依然显示的是默认的Cavity type          : VdW (van der Waals Surface) (Alpha=1.0 ...

G16 B.01也有这个问题

Aridea

刚刚也遇到同样报错问题，
感谢前辈的解决办法，我去试一试，不行我会回来告诉大家的，成功就不过来了，哈哈