计算化学公社

标题: Inv3 failed in PCMMkU (L502/L508 报错) 的成因与解决方法 [打印本页]
作者
Author: liyuanhe211    时间: 2017-9-6 02:00
标题: Inv3 failed in PCMMkU (L502/L508 报错) 的成因与解决方法
本帖最后由 liyuanhe211 于 2017-9-7 14:56 编辑

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

  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.

  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）提供的三种溶质空腔示意图：
(, 下载次数 Times of downloads: 167)


关键词：

更换溶剂表面的类型用 scrf=read 关键词启用溶剂化控制的段落，并在相应段落声明 surface=xxx，如 surface=SAS。
示例如下：
  1. #p pbe1pbe/genecp opt freq
  2. empiricaldispersion=gd3bj
  3. scrf=(smd,solvent=water,read)

  5. title

  7. 0 1
  8. C                   1.629576   -0.586781    1.988383
  9. ....
  10. H                   3.604870   -0.932068   -2.049232

  12. [Basis set definition]
  13. ...

  15. [Pseudopotential definition]
  16. ...

  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
复制代码






作者
Author: sobereva    时间: 2017-9-6 05:38
附上一张图便于读者理解
(, 下载次数 Times of downloads: 169)
作者
Author: kyuu    时间: 2017-9-6 06:01
sobereva 发表于 2017-9-6 05:38
附上一张图便于读者理解

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


稍有判断能力的一看就知道怎么回事，没必要注明
作者
Author: sobereva    时间: 2017-9-6 06:44
貌似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都存在。可以向官方反映下。
作者
Author: liyuanhe211    时间: 2017-9-6 14:31
sobereva 发表于 2017-9-6 06:44
貌似surface=SES没用，依然显示的是默认的Cavity type          : VdW (van der Waals Surface) (Alpha=1.0 ...

向高斯help反应了 SES/vdW 的问题，和 G16 的 scrf=G09Defaults 的问题
作者
Author: liyuanhe211    时间: 2017-9-7 14:59
sobereva 发表于 2017-9-6 06:44
貌似surface=SES没用，依然显示的是默认的Cavity type          : VdW (van der Waals Surface) (Alpha=1.0 ...

Gaussian Help 回复 Surface=SES 还需要加上 AddSph 选项程序才会用 SES 空腔。经测确实如此。
(, 下载次数 Times of downloads: 146)




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

噫，这手册真是不明不白
作者
Author: liyuanhe211    时间: 2017-9-8 08:33
sobereva 发表于 2017-9-8 03:18
噫，这手册真是不明不白

手册都没说这词→_→溶剂化整个一部分的手册都不像高斯写的似得。
作者
Author: lei234    时间: 2018-3-31 03:24
老师您好，以下是我理解该问题的解决办法，请您过目，如有不对，麻烦予以改正，在此谢过！
1）进行sas 优化，关键词 scrf=(smd,solvent=water,read) ，在最后加上surface=sas；
2）SES单点计算，关键词同1），最后加上Surface=SES  AddSph。
作者
Author: 赵云跳槽    时间: 2018-3-31 10:50
lei234 发表于 2018-3-31 03:24
老师您好，以下是我理解该问题的解决办法，请您过目，如有不对，麻烦予以改正，在此谢过！
1）进行sas 优 ...

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

(, 下载次数 Times of downloads: 172)
10.1021/acs.jpcb.7b06693

作者
Author: lei234    时间: 2018-4-1 01:59
赵云跳槽 发表于 2018-3-31 10:50
SAS计算，关键词 scrf=read，在最后加上surface=sas；
SES计算，关键词 scrf=read，最后加上addsph, rmi ...

谢谢！
这两个参数，rmin和ofac，是对水而言是吗？那溶剂是别的怎么办呢？这段话的意思是不是溶剂是水的时候可以不加这两个？
作者
Author: 赵云跳槽    时间: 2018-4-1 08:51
lei234 发表于 2018-4-1 01:59
谢谢！
这两个参数，rmin和ofac，是对水而言是吗？那溶剂是别的怎么办呢？这段话的意思是不是溶剂是水的 ...

那就不清楚了，只能自己尝试了
作者
Author: leebo    时间: 2019-1-27 21:56
sobereva 发表于 2017-9-6 06:44
貌似surface=SES没用，依然显示的是默认的Cavity type          : VdW (van der Waals Surface) (Alpha=1.0 ...

G16 B.01也有这个问题
作者
Author: Aridea    时间: 2019-10-7 15:49
刚刚也遇到同样报错问题，
感谢前辈的解决办法，我去试一试，不行我会回来告诉大家的，成功就不过来了，哈哈
作者
Author: Aridea    时间: 2019-10-9 12:45
前辈，用surface=sas 果然能够opt完成，但是这个结构能直接用吗——我想接下来算freq得到其Gibbs自由能校正项，随后计算高精度气相单点能，会不会不准。还请前辈赐教
作者
Author: liyuanhe211    时间: 2019-10-9 13:46
Aridea 发表于 2019-10-9 12:45
前辈，用surface=sas 果然能够opt完成，但是这个结构能直接用吗——我想接下来算freq得到其Gibbs自由能校正 ...

你拿几个两种都能运行完的比较一下就知道了
作者
Author: Aridea    时间: 2019-10-9 19:44
liyuanhe211 发表于 2019-10-9 13:46
你拿几个两种都能运行完的比较一下就知道了

嗯额，谢谢前辈
作者
Author: sunjk    时间: 2021-5-11 15:40
请问老师，我在优化结构时，加surface=sas的和不加的最终优化构型，金属键长差了超过0.1nm。也就是说，加了surface=sas后，对结构优化的结果是有影响。这样加了surface=sas后的优化结果，其构型，能量等是不是就不能和其它的没加的相比较了？
作者
Author: 量化小菜鸡    时间: 2021-5-11 16:22
sunjk 发表于 2021-5-11 15:40
请问老师，我在优化结构时，加surface=sas的和不加的最终优化构型，金属键长差了超过0.1nm。也就是说，加了 ...

你比较能量和自由能矫正了吗，我也遇到这种问题。我最后优化清一色全改IEFPCM去优化了，最后SMD做单点。
作者
Author: liyuanhe211    时间: 2021-5-11 17:28
sunjk 发表于 2021-5-11 15:40
请问老师，我在优化结构时，加surface=sas的和不加的最终优化构型，金属键长差了超过0.1nm。也就是说，加了 ...

键长一共多长，差0.1nm还有键吗？
作者
Author: CCCC    时间: 2021-11-22 09:17
g09 D.01也报了同样的错误，来试试前辈们的解决办法
作者
Author: pzl    时间: 2022-2-20 02:12
用bp86(d)/def2svp 算了一个乙烷的，发现在sas孔洞模型的时候，能量高一点，可能是体积大了，相互极化的小一点。然后ses和vdw差不多，因为体积几乎没差
作者
Author: ionexchangeC    时间: 2022-3-19 22:04
我在计算中发现溶剂空腔问题会导致原本高对称性的结构出现不合理的虚频。（B3LYP/def2-TZVP级别）在气相下优化[Ca(H2O)6]2+水合钙离子的结构，得到的结构为Th点群，振动分析后没有虚频；但改在溶剂模型下优化，得到的Th点群分子出现了6个虚频（使用SMD溶剂模型优化时则出现了21个虚频）。破坏对称性后继续优化，结构震荡难以收敛，尝试各种手段时偶然出现该报错，于是用GaussView打开输出文件观看溶剂化表面，发现取的空腔明显偏小（SMD则更小）。重新用高对称性的结构加上surface=SAS在溶剂模型优化，优化后Th点群保持并且没有虚频。也顺便试了下加上关键词surface=SES优化，结果在第一步优化的SCF迭代过程中就提示“Density matrix breaks symmetry”。 (, 下载次数 Times of downloads: 40)
作者
Author: GoldenBaby    时间: 2022-7-24 12:07
量化小菜鸡 发表于 2021-5-11 16:22
你比较能量和自由能矫正了吗，我也遇到这种问题。我最后优化清一色全改IEFPCM去优化了，最后SMD做单点。

好像IEFPCM也会出现这种情况，目前个人觉得比较好的解决方案还是用能算的模型跳过不行的结构再接着切回来吧
作者
Author: laplacesuanzi    时间: 2022-10-2 22:00
本帖最后由 laplacesuanzi 于 2022-10-2 22:15 编辑

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


作者
Author: Winnie1998    时间: 2022-10-16 11:34
各位老师好，我在结构优化的时候出现了上述报错，结构优化已经进行到第58帧，请问是要选取最后一帧的结构修改关键词后继续优化吗？还是要重新开始？
作者
Author: JCenter    时间: 2023-8-10 11:11
老师，我有一个疑问。当我进行结构优化时出现这个报错时，我在使用surface=sas完成结构优化时，还需要使用ses表面计算单点吗。我直接用sas优化的结构在更高的计算级别单点能，使用默认的表面不可以吗。
作者
Author: S221100847    时间: 2023-12-28 17:21
您好，想请教下，我这样写这个功能可以被使用吗
作者
Author: zhang123    时间: 2024-5-9 10:00
您好，请问您的问题解决了嘛？我是用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表面，不知道为啥还会报错呢？求大佬指点一下
作者
Author: Freeman    时间: 2024-5-10 14:52
本帖最后由 Freeman 于 2024-5-10 14:55 编辑

奇怪啊。手册上说SAS可以在SCRF关键词里指定，只能算单点能，跟高斯客服说法不一样啊。
(, 下载次数 Times of downloads: 12)
而且我发现有时结构和相对能量，加不加SAS差别很大，到了影响定性分析的程度。
作者
Author: 日暮知心    时间: 2024-7-26 10:18
请问这种问题能不能通过scrf=iterative解决？我试过前面提到的方法但还是报同样的错
作者
Author: pingju    时间: 2024-8-30 21:07
今天我也遇到了这个问题，博主您这个问题怎么解决的？请教
作者
Author: pingju    时间: 2024-8-31 10:47
lei234 发表于 2018-3-31 03:24
老师您好，以下是我理解该问题的解决办法，请您过目，如有不对，麻烦予以改正，在此谢过！
1）进行sas 优 ...

我用您这个方法，不行



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3