计算化学公社

标题: TRAH-SCF无法收敛的体系一例——高自旋钐Sm配合物 [打印本页]

作者
Author: Freeman 时间: 2023-6-6 15:00
标题: TRAH-SCF无法收敛的体系一例——高自旋钐Sm配合物
ORCA的TRAH方法很厉害，许多难以SCF收敛的体系，用上TRAH，问题都能迎刃而解。虽然TRAH远比常规SCF耗时，但它毕竟提供了一个可靠的最终解决方案。另一方面，稀土元素的SCF一直是很困难的，因为稀土元素有很高的轨道角动量和（往往）很高的自旋态。那么TRAH对高自旋态稀土元素的表现如何呢？

用PBE0-D4/def2-SVP算SmI2(THF)5七重态单点能，结构是从https://pubs.acs.org/doi/10.1021/acs.inorgchem.2c03752取的，优化理论水平为B3LYP-D3/6-311G(d,p)+SDD。算高自旋稀土元素时，最好先用G16算个片段组合波函数作为初猜，再用Multiwfn将chk转换为gbw文件，让ORCA读取gbw初猜波函数，方便SCF最终收敛到我们想要的电子组态上。

1、默认的SCF策略
不加任何多余的关键词，使用ORCA5默认的SCF策略。在第20次普通SCF迭代后自动进入TRAH，又进行了1458次迭代，最终不收敛。见sm1.orca.inp和sm1.orca.out。

2、不用TRAH
使用

! notrah

复制代码

关键词，表示不用TRAH。20次普通SCF迭代后收敛，而且结果显示Sm上有6个多余自旋，与初猜一致。这个策略在我算过的其他体系里，有时导致虽然收敛了，但轨道梯度[F,P]依然很大，表示其实能量还有进一步下降的空间。所以这样做的话，收敛后最好查看一下轨道梯度。见sm2.orca.inp和sm2.orca.out。

3、避免NR步
观察1中的TRAH过程

grep 2413 sm1.orca.out

复制代码

我们发现TRAH初期使能量下降，但下降到一定程度，就会出现NR步，而NR步使能量持续上升，上升到一定程度就又触发TRAH了。如此循环往复，最终导致不收敛。所以我们可以试着避免触发NR，一直TRAH到收敛。要规避NR，使用以下关键词。

%scf
trah quadregionstart 1.e-5 end # 默认值为1.e-3
end

复制代码

结果显示能量确实一直在下降，但是下降到一个点就下不去了，TRAH步一直在被reject，最后还是不收敛。但是我在其他Sm体系中应用此策略，有时奏效。见sm3.orca.inp和sm3.orca.out。

我在论坛上看到说TRAH在数学上可证明是令所有体系收敛的，但是我对这个Sm体系测试结论是，TRAH可能还没有普通SCF奏效。这是怎么回事呢？

(, 下载次数 Times of downloads: 22)

作者
Author: zjxitcc 时间: 2023-6-6 15:15
我问一下哈，你是已经用其他办法收敛了、只是想探究一下有无更好的收敛策略？还是尝试过的办法都没法收敛，现在想知道怎么将其收敛？

作者
Author: wzkchem5 时间: 2023-6-6 15:16
格点不够大。尝试defgrid3甚至refgrid，应该能收敛。
之前也有人在论坛上测过某个4f8构型的稀土离子，也是TRAH不收敛，但加refgrid就收敛了。
一般SCF收敛算法都要求Fock矩阵元等于能量对密度矩阵元的导数，但TRAH对这个条件的依赖比其他收敛算法要强，但凡因为数值误差导致这个条件被破坏，就很容易不收敛，因为程序会发现，本来按Fock矩阵判断会让能量降低的轨道旋转，实际反而会让能量变高，于是程序就不断降低trust radius，直至完全停滞不前。

作者
Author: Freeman 时间: 2023-6-6 15:32

zjxitcc 发表于 2023-6-6 15:15
我问一下哈，你是已经用其他办法收敛了、只是想探究一下有无更好的收敛策略？还是尝试过的办法都没法收敛， ...

已经用常规SCF收敛完了。想搞清楚为什么TRAH反而不管用。

作者
Author: Freeman 时间: 2023-6-6 15:35

wzkchem5 发表于 2023-6-6 15:16
格点不够大。尝试defgrid3甚至refgrid，应该能收敛。
之前也有人在论坛上测过某个4f8构型的稀土离子，也是 ...

谢谢解答。那么为什么不含稀土元素的体系就（大概）没有此问题呢？积分格点为什么对稀土元素影响那么大呢？

作者
Author: wxsxjj 时间: 2023-6-6 17:16
与题无关，我前几天刚算过类似的Sm体系，没猜错你要算SmI2的还原反应？
我用Gaussian16 A.03 SCF=XQC没遇到问题。同样ORCA的TRAH收敛有问题，后来干脆就都用Gaussian了

作者
Author: Freeman 时间: 2023-6-6 18:53

wxsxjj 发表于 2023-6-6 17:16
与题无关，我前几天刚算过类似的Sm体系，没猜错你要算SmI2的还原反应？
我用Gaussian16 A.03 SCF=XQC没遇 ...

不是还原反应。我也这么尝试过，但有时scf=xqc也不收敛，而且g16的scf速度比不上orca的rijcosx。

作者
Author: wzkchem5 时间: 2023-6-6 19:28

Freeman 发表于 2023-6-6 11:53
不是还原反应。我也这么尝试过，但有时scf=xqc也不收敛，而且g16的scf速度比不上orca的rijcosx。

话说，我觉得最好能普及一下把orca波函数转成chk格式、作为高斯初猜的workflow，现在sob老师写了把高斯波函数传给orca的工具、有很多用户在用，虽然原则上把orca波函数传给高斯也一样容易，但是没有人把这种方法给大家demonstrate出来，导致的后果就是orca不收敛、大家都想着去尝试高斯，如果发现高斯也不收敛就是体系太复杂，如果高斯收敛了就是orca不如高斯；但是如果高斯scf不收敛，没什么人想着尝试orca，这样即便orca收敛概率比高斯更高，大家也发现不了。以前orca的SCF收敛确实做得不如高斯，但是现在真的未必如此了（尤其如果不看SCF迭代次数，而是看总SCF时间，因为RIJCOSX等算法经常会轻微牺牲迭代次数换取加速），我们orca团队的人又因为没有高斯license，没法做benchmark，做了也因为众所周知的原因不能发出来。

作者
Author: wzkchem5 时间: 2023-6-6 19:32
本帖最后由 wzkchem5 于 2023-6-6 12:34 编辑

Freeman 发表于 2023-6-6 08:35
谢谢解答。那么为什么不含稀土元素的体系就（大概）没有此问题呢？积分格点为什么对稀土元素影响那么大呢 ...

因为4f6、4f8电子构型，多余的那个f electron或f hole不管在哪个f轨道，能量都极其接近，但因为配位环境不对称以及格点误差，能量又不完全一样。所以极小的能量或Fock矩阵误差就会使得波函数有极大的变动。这种特殊电子结构对误差的放大作用，和TRAH对数值误差要求严格这个问题叠加在一起，等于雪上加霜，导致实际发现不收敛。
其实我在想，反正这种体系高精度计算也要用多参考态方法，如果用DFT优化结构的话，或许可以考虑Fermi smearing，在有限电子温度下做计算？因为稀土的f能级极其接近，应该在室温下f轨道占据数就已经分布得非常均匀了，会极大有利于SCF收敛，而且又因为你设的电子温度就是实际反应温度，你又没有引入任何近似，甚至可以说比电子温度0K的计算还靠谱一点。不过这个只是从理论上觉得可行，我不知道有没有文献先例支持

作者
Author: zjxitcc 时间: 2023-6-6 20:23
本帖最后由 zjxitcc 于 2023-6-6 20:25 编辑

wzkchem5 发表于 2023-6-6 19:28
话说，我觉得最好能普及一下把orca波函数转成chk格式、作为高斯初猜的workflow，现在sob老师写了把高斯波 ...

MOKIT的mkl2gjf小程序就可以把ORCA的轨道传给Gaussian，有些用户使用过。不过，由于mkl文件里没有赝势数据，所以如果用户使用赝势，产生的gjf文件就不含赝势（非mkl2gjf所能解决），需要用户手动添加，全电子基组无此问题。所以要想推广ORCA -> Gaussian，我的建议是ORCA开发者要扩展一下mkl文件，加进赝势数据。

作者
Author: chrinide 时间: 2023-6-6 21:43

zjxitcc 发表于 2023-6-6 20:23
MOKIT的mkl2gjf小程序就可以把ORCA的轨道传给Gaussian，有些用户使用过。不过，由于mkl文件里没有赝势数 ...

mkl这个格式设计得..........，体系稍微大一点那转化速度简直不要忒感人了

作者
Author: wzkchem5 时间: 2023-6-6 21:56

zjxitcc 发表于 2023-6-6 13:23
MOKIT的mkl2gjf小程序就可以把ORCA的轨道传给Gaussian，有些用户使用过。不过，由于mkl文件里没有赝势数 ...

我也是觉得molden比mkl方便一些，印象中molden的有效数字好像也比默认导出的mkl多？
加赝势数据的问题你可以在orca论坛上提一下

作者
Author: zjxitcc 时间: 2023-6-6 22:57
本帖最后由 zjxitcc 于 2023-6-6 22:59 编辑

chrinide 发表于 2023-6-6 21:43
mkl这个格式设计得..........，体系稍微大一点那转化速度简直不要忒感人了

你是指orca_2mkl对大文件转换比较慢是吧，不是指mkl2gjf

作者
Author: wxsxjj 时间: 2023-6-6 23:05

wzkchem5 发表于 2023-6-6 19:28
话说，我觉得最好能普及一下把orca波函数转成chk格式、作为高斯初猜的workflow，现在sob老师写了把高斯波 ...

multiwfn里example/scripts里面有gbw2chk脚本应该可以，然后guess=read。不知道您是不是这个意思

作者
Author: wzkchem5 时间: 2023-6-6 23:49

wxsxjj 发表于 2023-6-6 16:05
multiwfn里example/scripts里面有gbw2chk脚本应该可以，然后guess=read。不知道您是不是这个意思

对，我就是这个意思。我感觉这种工具应该普及，不然的话高斯不收敛但orca收敛的体系大家根本没机会发现，别人一搜咱们论坛，全都是orca不收敛、用高斯收敛以后传给orca做初猜的例子。甚至很多人把高斯预先收敛当作默认workflow了，万一自己的体系用orca直接收敛比高斯收敛得还快，都不知道。别的用户一看这么多人用高斯收敛scf以后把波函数传给orca，心想肯定是orca的scf收敛做得不好，其实未必如此。

作者
Author: chrinide 时间: 2023-6-7 06:30

zjxitcc 发表于 2023-6-6 22:57
你是指orca_2mkl对大文件转换比较慢是吧，不是指mkl2gjf

对头

作者
Author: biogon 时间: 2023-6-7 08:19
这个体系可能不用trah，用veryslowconv反而能收敛

作者
Author: ionexchangeC 时间: 2023-6-7 14:40
我记得multiwfn可以把ORCA收敛的波函数转换为Gaussian格式啊

作者
Author: sobereva 时间: 2023-6-7 14:44

ionexchangeC 发表于 2023-6-7 14:40
我记得multiwfn可以把ORCA收敛的波函数转换为Gaussian格式啊

是的。molden经过Multiwfn转成mkl，再经过orca_2mkl转成gbw。10行左右shell脚本的事

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