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

复制代码

pingju

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

我用您这个方法，不行

pingju

今天我也遇到了这个问题，博主您这个问题怎么解决的？请教

日暮知心

请问这种问题能不能通过scrf=iterative解决？我试过前面提到的方法但还是报同样的错

Freeman

奇怪啊。手册上说SAS可以在SCRF关键词里指定，只能算单点能，跟高斯客服说法不一样啊。

而且我发现有时结构和相对能量，加不加SAS差别很大，到了影响定性分析的程度。

zhang123

您好，请问您的问题解决了嘛？我是用gaussian 16B.01算的单点能，输入文件如下：
%oldchk=Co_P_2_1_gas-sp.chk
%chk=Co_P_2_1_solv-sp.chk
#p ub3lyp/6-31+g* guess=read geom=allcheck scrf=g03defaults，也出现了一样问题，我看了一下G03默认的就是SES表面，不知道为啥还会报错呢？求大佬指点一下

S221100847

您好，想请教下，我这样写这个功能可以被使用吗

JCenter

老师，我有一个疑问。当我进行结构优化时出现这个报错时，我在使用surface=sas完成结构优化时，还需要使用ses表面计算单点吗。我直接用sas优化的结构在更高的计算级别单点能，使用默认的表面不可以吗。

Winnie1998

各位老师好，我在结构优化的时候出现了上述报错，结构优化已经进行到第58帧，请问是要选取最后一帧的结构修改关键词后继续优化吗？还是要重新开始？

laplacesuanzi

我最近也遇到了相同的问题，一模一样的报错。关键词最开始为“#p opt freq b3lyp scrf=(solvent=generic,read) em=gd3bj tzvp”，通过参考http://sobereva.com/61，我在关键词中又加上了“SCF=noincfock”，这个问题就解决了。感谢卢天老师。

GoldenBaby

量化小菜鸡 发表于 2021-5-11 16:22
你比较能量和自由能矫正了吗，我也遇到这种问题。我最后优化清一色全改IEFPCM去优化了，最后SMD做单点。

好像IEFPCM也会出现这种情况，目前个人觉得比较好的解决方案还是用能算的模型跳过不行的结构再接着切回来吧

ionexchangeC

我在计算中发现溶剂空腔问题会导致原本高对称性的结构出现不合理的虚频。（B3LYP/def2-TZVP级别）在气相下优化[Ca(H2O)6]2+水合钙离子的结构，得到的结构为Th点群，振动分析后没有虚频；但改在溶剂模型下优化，得到的Th点群分子出现了6个虚频（使用SMD溶剂模型优化时则出现了21个虚频）。破坏对称性后继续优化，结构震荡难以收敛，尝试各种手段时偶然出现该报错，于是用GaussView打开输出文件观看溶剂化表面，发现取的空腔明显偏小（SMD则更小）。重新用高对称性的结构加上surface=SAS在溶剂模型优化，优化后Th点群保持并且没有虚频。也顺便试了下加上关键词surface=SES优化，结果在第一步优化的SCF迭代过程中就提示“Density matrix breaks symmetry”。

pzl

用bp86(d)/def2svp 算了一个乙烷的，发现在sas孔洞模型的时候，能量高一点，可能是体积大了，相互极化的小一点。然后ses和vdw差不多，因为体积几乎没差

CCCC

g09 D.01也报了同样的错误，来试试前辈们的解决办法

liyuanhe211

sunjk 发表于 2021-5-11 15:40
请问老师，我在优化结构时，加surface=sas的和不加的最终优化构型，金属键长差了超过0.1nm。也就是说，加了 ...

键长一共多长，差0.1nm还有键吗？