计算化学公社

标题: 请教关于FrozenCore设置引起的误差的问题 [打印本页]

作者
Author: DoorBell 时间: 2020-10-6 00:51
标题: 请教关于FrozenCore设置引起的误差的问题
如题~

用Gaussian 16 A.03在CCSD(T)/aug-cc-pVTZ级别计算金属杂硼团簇阴离子不同异构体之间的相对能量的时候，发现默认的FrozenCore设置和加上Full关键字得到的两个异构体相对能量有一定偏差（几kcal/mol）。由于本身两个异构体能量较为接近，这样的偏差导致定性判断能量高低结果不同，请问这时候应该相信谁的结果呢？

关键词均为
# CCSD(T,T1Diag)/GenECP SCF=verytight 和 # CCSD(T,T1Diag,Full)/GenECP SCF=verytight
基组均为金属aVTZ-PP，硼aVTZ

非常感谢！

作者
Author: LGC020305 时间: 2020-10-6 05:03
本帖最后由 LGC020305 于 2020-10-6 05:06 编辑

1. 给你转一段社长的答复（不针对你的具体问题）： https://www.researchgate.net/pos ... gaussian-16_package
17th Dec, 2017
Tian Lu， Beijing Kein Research Center for Natural Sciences

The number of core orbitals of each element is already known according to chemistry knowledge, since lowest lying MOs are mainly composed of core atomic orbitals, the number of lowest MOs to be frozen can be simply obtained by counting all core orbitals of atoms presented in current system.
By default, in Gaussian, all atomic orbitals except for the outermost valence ones are regarded as core atomic orbitals. This typically works reasonably. However, in rare case this treatment may cause unsatisfactory result; for example, if you optimize Li2, you will find the discrepancy between CCSD(T,full) and CCSD(T) is never small. Because Li only has one valence electron, electron correlation due to its 1s electrons also have nonnegligible contribution to chemically interesting quantites.
In Gaussian there are a few different ways of determining how many orbitals will be frozen, see http://gaussian.com/frozencore/. For example, with "FreezeInnerNobleGasCore" keyword, MOs composed of both outermost valence and sub-valence atomic orbitals will not be frozen.
If you are still confusing, you can try to calculate some elements and simple molecules, the number of frozen MOs can be found after "NFC=" label in output file.

2. Full
This specifies that all electrons be included in a correlation calculation.
FC
This specifies a “frozen core” calculation, and it implies that inner-shells are excluded from the correlation calculation. This is the default calculation mode.

3. 好比如一个是近似解，另一个不是。理论上Full 的更可靠

作者
Author: DoorBell 时间: 2020-10-6 09:14

LGC020305 发表于 2020-10-6 05:03
1. 给你转一段社长的答复（不针对你的具体问题）： https://www.researchgate.net/post/Frozen-cores_in_ga ...

感觉B也有可能有类似Li的问题啊...
不过Full就一定可靠吗？我看到社长在http://bbs.keinsci.com/thread-17401-1-1.html也提到了一个问题？

作者
Author: wzkchem5 时间: 2020-10-6 10:32

DoorBell 发表于 2020-10-6 09:14
感觉B也有可能有类似Li的问题啊...
不过Full就一定可靠吗？我看到社长在http://bbs.keinsci.com/thread- ...

你的金属是什么？一般金属core最外层轨道的影响比硼这一周期的core轨道影响大

作者
Author: beefly 时间: 2020-10-6 17:25
本帖最后由 beefly 于 2020-10-6 17:29 编辑

对于存在d占据轨道的元素来说，半芯轨道(n-1)d的关联能比较重要。原则上，full比默认的冻芯结果更准确，但是它把1s、2s、2p这些最内层轨道也考虑电子关联，一般是没有必要的（除非计算最内层电子有关的性质），白白浪费机时。最好的办法是设置FC(Window=?)，仅把(n-1)d以下的轨道进行冻结。计算条件允许的情况下也可以保留(n-1)s(n-1)p，只冻结(n-1)s以下轨道。

半芯轨道的关联能需要用特殊基组，例如cc-pCVnZ、cc-pWCVnZ系列，或者把半芯轨道对应的几个主要高斯函数做非收缩处理。一般的基组达不到目的，例如有些文献用的CCSD(T,full)/6-311++G**就是瞎算。

作者
Author: DoorBell 时间: 2020-10-6 18:22
本帖最后由 DoorBell 于 2020-10-6 19:15 编辑

beefly 发表于 2020-10-6 17:25
对于存在d占据轨道的元素来说，半芯轨道(n-1)d的关联能比较重要。原则上，full比默认的冻芯结果更准确，但 ...

金属是Yttrium，好像没有像样的全电子基组/能包括(n-1)d的基组。
譬如说cc-pwCVTZ-PP对Y有定义，但是cc-pwCVTZ没有，而且前者对Y仅有11个价电子，并不包括(n-1)d，似乎有点棘手了
另外cc-pwCVTZ只对第一过渡系有定义，这些金属也没有(n-1)d欸~

作者
Author: DoorBell 时间: 2020-10-6 18:22

wzkchem5 发表于 2020-10-6 10:32
你的金属是什么？一般金属core最外层轨道的影响比硼这一周期的core轨道影响大

Y, Yttrium

作者
Author: wzkchem5 时间: 2020-10-6 20:19

DoorBell 发表于 2020-10-6 18:22
金属是Yttrium，好像没有像样的全电子基组/能包括(n-1)d的基组。
譬如说cc-pwCVTZ-PP对Y有定义，但是cc- ...

一般来说，ECP忽略内层电子的影响和相对论效应的量级应该类似，所以如果考虑改成全电子，那么应该同时也显式考虑标量相对论，这样才比较balanced。
如果一定要用高斯，那么可以用cc-pwCVTZ-DK（需要从basis set exchange上面下载），搭配标量DKH2相对论方法。不过我没用过高斯的DKH2，不知道靠谱性。
另一种我个人比较推荐的方法是用ORCA算，可以用ZORA或DKH2，基组有cc-pwCVTZ-DK、ANO-RCC-TZP等多种选择。而且关键是ORCA的DLPNO-CCSD(T)比一般的CCSD(T)快得多，内存消耗也小得多，triple zeta基组可以算到100个原子，可以大大加快你的计算；不仅如此，还有explicitly correlated的CCSD(T)方法（F12-CCSD(T)），可以以TZ的计算量得到5Z左右的精度。

作者
Author: DoorBell 时间: 2020-10-6 20:53

wzkchem5 发表于 2020-10-6 20:19
一般来说，ECP忽略内层电子的影响和相对论效应的量级应该类似，所以如果考虑改成全电子，那么应该同时也 ...

非常感谢大神！
现在主要是需要比较不同异构体的能量，所以对即使1kcal/mol的误差都是很敏感的。
速度不是主要矛盾，可信度现在是最亟需的。

另外，我已经听好几个搞理论的老师说不相信DLPNO-CCSD(T)甚至ORCA的正则CCSD(T)结果了（超小声）

作者
Author: zjxitcc 时间: 2020-10-6 21:13
本帖最后由 zjxitcc 于 2020-10-6 21:17 编辑

DoorBell 发表于 2020-10-6 20:53
非常感谢大神！
现在主要是需要比较不同异构体的能量，所以对即使1kcal/mol的误差都是很敏感的。
速度 ...

（1）“听好几个搞理论的老师说不相信DLPNO-CCSD(T)甚至ORCA的正则CCSD(T)结果了”其实不需要听别人的这些话，因为每个人对ORCA的关键词熟悉程度是不一样的，有的学生做错计算了也不知道。到底靠不靠谱，自己拿个中小有机体系算个例子，把精度关键词（SCF收敛限，CCSD冻核数目等等）都控制好，看与高斯一不一样就行。一样的话，别人说破嘴也没用（如果不一样，可以另开一贴大家分析）。
另外DLPNO-CCSD(T)在PNO上也有NormalPNO, TightPNO等等区别，精度要高至少也得TightPNO；(T)也细分为T0和T1，精度要高得用T1。不过DLPNO-CCSD(T)是逼近CCSD(T)的，当你算的动CCSD(T)且觉得它有问题时，算DLPNO-CCSD(T)显然没用。

（2）当你的体系对于冻/不冻核->改变结论时，说明冻核可能不是关键性因素，可能是你的方法需要有所改进，例如考虑多参考特征是否较强、考虑相对论效应（相对论赝势RECP，全电子DKH2, 精确二分量x2c等等），以及二者的结合。

另外说一点基本的，你的体系是开壳层，做的是UCCSD(T)？你的UHF波函数是稳定的么？

作者
Author: wzkchem5 时间: 2020-10-6 22:37

DoorBell 发表于 2020-10-6 20:53
非常感谢大神！
现在主要是需要比较不同异构体的能量，所以对即使1kcal/mol的误差都是很敏感的。
速度 ...

1kcal/mol量级的误差的话，必须用全电子基组，显式考虑相对论，赝势就别想了。而且团簇搞不好有强相关，一定密切注意T1 diagnostic，如果太大的话CCSD(T)都不靠谱。
至于ORCA的CCSD(T)的可信度，起码正则CCSD(T)不可信是瞎扯。有的可能是用6-31G(d)跟高斯比，然后忘了高斯用笛卡尔d函数而ORCA用球谐d函数。又比如说是高斯和ORCA的HF都没有收敛到同一个解，比如这个https://orcaforum.kofo.mpg.de/viewtopic.php?f=11&t=4185，甚至高斯收敛到的解比ORCA的能量还高，他也不去分析为什么，只看到两个算出来的CCSD(T)结果有差别，就怀疑ORCA有bug，然后可能在圈子里就传开了。Neese固然没有高斯的license，写程序时不方便直接和高斯比较，但是他有Molpro等等其他支持CCSD(T)的程序可以比较，不至于连正则CCSD(T)都做错。
至于DLPNO-CCSD(T)，毕竟DLPNO方法独此一家，不看代码没人能说做的是对的，但是有一点是没错的，就是当收紧各种收敛阈值时，DLPNO-CCSD(T1)结果能无限趋近正则CCSD(T)的结果，这个测过很多了。从这个意义上讲，DLPNO-CCSD(T)如果结果不可信，那一定是阈值设得太松。这里需要注意的是Neese文章里说NormalPNO误差一般在1kcal/mol左右，有人拿NormalPNO算出来3kcal/mol误差，可能就断定DLPNO一无是处了，其实根本不是这么回事，你去看这种体系的T1 diagnostic，可能都超过0.02了。对于T1 diagnostic比较大的体系，DLPNO-CCSD(T)误差容易显著高于平均误差，因为(T)项的分母里有轨道能级差，对于强相关体系，分母小，所以会把分子的误差放大。但这种情况下CCSD(T)本身就不靠谱。有人拿CCSD(T)强行算强相关体系，然后发现DLPNO-CCSD(T)和正则CCSD(T)差别很大，这个其实说明不了什么。

作者
Author: wzkchem5 时间: 2020-10-6 22:40

zjxitcc 发表于 2020-10-6 21:13
（1）“听好几个搞理论的老师说不相信DLPNO-CCSD(T)甚至ORCA的正则CCSD(T)结果了”其实不需要听别人的这 ...

这种团簇体系可以考虑用Kohn-Sham轨道做CCSD(T)的reference，靠CCSD(T)的T1算符来起到orbital relaxation的作用。HF即使波函数稳定也可能定性错误，比如我就遇到过一个阴离子体系，HF的HOMO是个Rydberg轨道，完全不对，但B3LYP轨道是定性正确的，所以用B3LYP做reference来算CCSD(T)。

作者
Author: zjxitcc 时间: 2020-10-6 23:20
本帖最后由 zjxitcc 于 2020-10-6 23:21 编辑

wzkchem5 发表于 2020-10-6 22:40
这种团簇体系可以考虑用Kohn-Sham轨道做CCSD(T)的reference，靠CCSD(T)的T1算符来起到orbital relaxation ...

这有点风险吧，取决于CCSD(T)代码的具体实现。KS轨道相当于未收敛的HF轨道，存在非零Fock矩阵元f_{ia}，如果用的CCSD(T)代码里没有考虑f_{ia}（大部分方法文献假设是收敛的HF正则轨道，不考虑f_{ia}），那实际上相当于少了一部分贡献，结果好的话，存在误差相消的可能性。
这种情况我会选择（1）继续尝试寻找其他稳定的HF解；（2）使用多参考方法。当然了，这是个人选择了。

作者
Author: wzkchem5 时间: 2020-10-6 23:51

zjxitcc 发表于 2020-10-6 23:20
这有点风险吧，取决于CCSD(T)代码的具体实现。KS轨道相当于未收敛的HF轨道，存在非零Fock矩阵元f_{ia}， ...

你说得对，需要看程序说明书里有没有提是否假设f_{ia}=0。我只知道ORCA是不假设f_{ia}=0的

作者
Author: sobereva 时间: 2020-10-7 00:38

DoorBell 发表于 2020-10-6 09:14
感觉B也有可能有类似Li的问题啊...
不过Full就一定可靠吗？我看到社长在http://bbs.keinsci.com/thread- ...

full一定更可靠的前提是用了能体现核电子相关的基组
对于当前情况，建议先把Y用aug-cc-pwCVTZ-PP，B用aug-cc-pwCVTZ，结合full做计算看看情况，起码比你之前做的计算更可靠。想更可靠再考虑提到QZ、全电子相对论或多参考方法（参考T1诊断）方面的事。

作者
Author: DoorBell 时间: 2020-10-7 00:51

sobereva 发表于 2020-10-7 00:38
full一定更可靠的前提是用了能体现核电子相关的基组
对于当前情况，建议先把Y用aug-cc-pwCVTZ-PP，B用au ...

蟹蟹sob老师！您提到的计算已经在进行了，好慢好慢QAQ

作者
Author: beefly 时间: 2020-10-7 09:50
本帖最后由 beefly 于 2020-10-7 10:18 编辑

DoorBell 发表于 2020-10-6 18:22
金属是Yttrium，好像没有像样的全电子基组/能包括(n-1)d的基组。
譬如说cc-pwCVTZ-PP对Y有定义，但是cc- ...

前面是针对主族元素说的。如果是d区过渡元素，例如Ni，3d是价轨道，3s3p是半芯轨道，它们的电子关联比较重要。
然后查基组文献，看半芯轨道是否做了冻结处理。比如ANO-RCC基组，较重的元素做基组收缩的时候未冻结半芯轨道，而某些较轻的元素可能冻结了半芯轨道，后者就不适合做半芯轨道的电子关联计算。

如果不是用来做基组外推计算，自己改造基组也很容易。比如硼的def2-SV(P)，标准的收缩因子未考虑1s关联。先用非收缩的def2-SV(P)做B的HF或DFT计算（为了简单，算闭壳层的B+），查看需要做关联处理的1s轨道因子，发现在7.8793722710、2.4088857172两个高斯函数上的收缩因子最大（通常s，p取轨道因子最大的2-3个，d，f取因子最大的3-4个即可），这两个基函数再单独做个非收缩就行了。最后，B用来做1s芯电子关联的def2-SV(P)（为了区分，可以改名叫def2-C-SV(P)）就是

B 0
S  5  1.00
  839.31830086    -.55929201074E-02
  126.26464843    -.41565520771E-01
  28.620600763    -.18299816983
  7.8793722710    -.46540391866
  2.4088857172    -.44173884791
S  1  1.00
  .25105109036    1.0000000000
S  1  1.00
  .83648866069E-01  1.0000000000
P  3  1.00
  6.0332223619    -.35603672456E-01
  1.2499157866    -.19895775769
  .33871676350    -.50850202618
P  1  1.00
  .96415632351E-01  1.0000000000
D  1  1.00
  .50  1.00
! 以下两个函数是新加的
S  1  1.00
  7.8793722710    1.0
S  1  1.00
  2.4088857172    1.0
****

p、d半芯轨道的处理方法类似。

作者
Author: zjxitcc 时间: 2020-10-7 10:49

beefly 发表于 2020-10-7 09:50
前面是针对主族元素说的。如果是d区过渡元素，例如Ni，3d是价轨道，3s3p是半芯轨道，它们的电子关联比较 ...

借个楼请教一下，针对这个B自定义基组的例子，加完“! 以下两个函数是新加的”后，把原来S 5 1.00里第4个和第5个高斯函数删掉，然后把数字5改为3，这样做是等价的吧？我算了下看这样RHF和MP2(full)能量是不变的，应该就是相当于ExactBasisTransform与NoBasisTransform的区别

作者
Author: beefly 时间: 2020-10-11 15:57

zjxitcc 发表于 2020-10-7 10:49
借个楼请教一下，针对这个B自定义基组的例子，加完“! 以下两个函数是新加的”后，把原来S 5 1.00里第4 ...

B太轻，不明显。对于比较重的元素，两种收缩方式得到的能量、性质会有一点差异，特别是ANO类型的基组

作者
Author: zjxitcc 时间: 2020-10-11 16:26

beefly 发表于 2020-10-11 15:57
B太轻，不明显。对于比较重的元素，两种收缩方式得到的能量、性质会有一点差异，特别是ANO类型的基组

thanks!

欢迎光临计算化学公社 (http://bbs.keinsci.com/)