如何计算p-C6H4分子的自旋极化单重态？

AndrewX

p-C6H4分子即苯环的两个对位氢去掉后的结构，即双自由基结构。我在阅读Sob老师的博文 http://sobereva.com/82 后，了解到了计算自旋极化单重态的方法。我的具体做法：
第一步，将之前计算好的闭壳层单重态的p-C6H4结构（事先用stable检测过，有从RHF>UHF的不稳定性）另存为输入文件，写上关键词# stable=opt guess=mix b3lyp/6-311++G** nosymm 来寻找稳定的基态波函数，此时的计算文件名为p-C6H4-brokensymm；
第二步，直接读取上一步的chk文件信息进行结构优化加频率计算，关键词为#  opt freq ub3lyp/6-311++g(d,p) guess=read geom=check nosymm，此时的计算文件名为p-C6H4-openshelledsinglet；
第三步，用stable检测了上一步的波函数稳定性，此时的计算文件名为p-C6H4-openshelledsinglet-stable，提示“The wavefunction has an internal instability”。

我跟导师说我算出来开壳层单重态是p-C6H4的基态，但是他觉得我算的不对，他提出的观点详见附件中的图片。
输入和输出也都在附件中，麻烦老师看看我的操作是否哪里有错误？并且我有以下疑问：
第一步中，我忘记给b3lyp加u了，但是照样找到了对称破缺的波函数，所以“加了stable=opt guess=mix不加u”是否只会增加计算耗时而不会导致无法找到对称破缺波函数？
第二步中，读取对称破缺函数进行opt+freq计算时，是否一定要加u？我试了不加u，照样能优化得到开壳层单重态，但是怕我这个只是个特例，所以有此一问。
第三步中，我还算了其他分子的一些自旋极化单重态，检测波函数稳定性也都提示有internal instability，所以是自旋极化单重态都这样还是我算的有问题？
还有一个疑问：自旋极化单重态和开壳层单重态这两个概念指向的是否是同一个东西？

Ps：从Sob老师的博文已知6-311++G**基组不适合算有机体系，这里是想和之前的工作保持一致的计算水平，方便比较。

wzkchem5

这个体系有显著的多参考态特征，所以自旋极化单重态结果也不准，只不过比闭壳层单重态准得多而已，在个别你确认不需要高精度的场合、或者更高精度的方法算不动的情况下可以用。如果需要更高的精度，必须用多参考态方法，至少用NEVPT2起步。或者用SF-TDDFT以三重态双自由基为参考态来算单重态基态也可以。
自旋极化单重态不排除可能有internal instability，但是最稳定的自旋极化单重态不可能有任何instability。你这里的问题在于最后一步stable计算读波函数没读对，导致程序从默认初猜重新做scf得到了一个波函数，检验稳定性也是检验的这个重新计算的波函数，不是你chk文件里的波函数。
自旋极化单重态和开壳层单重态确实是一个意思。
6-311++G**不是不适合算有机体系，这个基组的主要问题在于对氢加了弥散，所以适合用来算有的氢带负电的体系，比如你要算LiAlH4还原醛酮就挺合适的。如果体系没有带负电的氢，那么这个基组不仅不适合算有机体系，也不适合算无机体系。是不是有机体系根本不是这里的重点

sobereva

我的博文从来没说过什么“6-311++G**基组不适合算有机体系”这种谬论
认真看此文，别乱加弥散函数
谈谈弥散函数和“月份”基组
http://sobereva.com/119

用stable=opt的时候没必要刻意加U。Gaussian如果发现闭壳层波函数不稳定，搜索稳定波函数的时候会自动尝试以U方式计算得到更低能量。
自旋极化单重态和开壳层单重态语义侧重点不一样，但你可以认为在实际中指的电子结构状态是相同的，都是alpha和beta电子结构不匹配。

AndrewX

wzkchem5 发表于 2022-4-13 19:01
这个体系有显著的多参考态特征，所以自旋极化单重态结果也不准，只不过比闭壳层单重态准得多而已，在个别你 ...

谢谢老师解答！所以我这里用计算普通自旋极化单重态的方法得到的结果是不是误差小到可以用在文章发表中？并且，用闭壳层单重态方法计算得到的结果是不是误差已经大到不能用在发表文章里？我的文章里C6H4只是解离产物之一，目的是计算反应的ΔE，用来解释气相中的解离反应实验现象的。这种情况是否是老师说的“不需要高精度的场合”？

wzkchem5

AndrewX 发表于 2022-4-13 17:54
谢谢老师解答！所以我这里用计算普通自旋极化单重态的方法得到的结果是不是误差小到可以用在文章发表中？ ...

闭壳层单重态的结果肯定是不能发表的。自旋极化单重态的结果当且仅当你用高精度方法算不动的时候才可以发表。如果你发表的话，审稿人可能会叫你补NEVPT2甚至MRCI计算，那么就会有一个问题，C6H4做多参考态计算肯定是能算得动的，你没有理由不补这个计算。那么与其这样，不如现在就确保结果是靠谱的，除了实际做一个多参考态计算以外，还有第二种方法，因为C6H4这个体系研究的人很多，你查已有的benchmark文章，确认一下你用的泛函可以把这个体系算准，然后引在文章里，就说你知道这类体系DFT一般算不准，但是对于你这个特定的体系，文献证明你用的泛函是能算准的。

AndrewX

sobereva 发表于 2022-4-13 23:46
我的博文从来没说过什么“6-311++G**基组不适合算有机体系”这种谬论
认真看此文，别乱加弥散函数
谈谈弥 ...

抱歉老师，我脑子抽了，本意是想说从博文中得知b3lyp不适合有机体系。

AndrewX

wzkchem5 发表于 2022-4-14 03:24
闭壳层单重态的结果肯定是不能发表的。自旋极化单重态的结果当且仅当你用高精度方法算不动的时候才可以发 ...

好的，我按您说的先找找相关的benchmark文章，谢谢老师！

wzkchem5

AndrewX 发表于 2022-4-14 04:50
抱歉老师，我脑子抽了，本意是想说从博文中得知b3lyp不适合有机体系。

这个也不对，B3LYP只是不适合一部分有机体系，比如特别大的共轭体系（对于小的共轭体系，比如萘这种大小的，也是适合的）。对于大部分有机体系，B3LYP是适合的。反倒是对于过渡金属体系，B3LYP不适合的概率比较大。
看博文一定要把原结论精确地记下来，不要只看个大概、后来回想博文里举的不适合的体系都是有机体系，就以为B3LYP不适合有机体系。

橘垚垚

wzkchem5 发表于 2022-4-14 03:24
闭壳层单重态的结果肯定是不能发表的。自旋极化单重态的结果当且仅当你用高精度方法算不动的时候才可以发 ...

非常感谢老师的回答，我想请问一下老师在利用诸如CASPT2的多参考方法对这种双自由基体系进行计算时（molpro软件），要形成对称破缺波函数需要怎么设置关键词呢？另外由于开壳层的计算得到的是自旋轨道，那在后续的CASSCF计算中确定活性空间该怎样指定轨道呢？因为说明书里活性空间的指定好像都是针对可以容纳两个电子的轨道而言的。谢谢老师啦

wzkchem5

橘垚垚 发表于 2022-5-9 04:38
非常感谢老师的回答，我想请问一下老师在利用诸如CASPT2的多参考方法对这种双自由基体系进行计算时（molp ...

对称破缺波函数是只对单参考方法才存在的概念，多参考方法不用且不能用对称破缺波函数，不需要对称破缺也能得到正确的结果。选活性轨道也是同理，对于不复杂的体系，用RHF轨道就行了。

橘垚垚

wzkchem5 发表于 2022-5-9 14:46
对称破缺波函数是只对单参考方法才存在的概念，多参考方法不用且不能用对称破缺波函数，不需要对称破缺也 ...

非常感谢老师的回复，您的回复很有帮助，确实CASSCF本身就是为了解决这个问题的存在。