计算化学公社

标题: 关于benchmark和泛函选择问题求助 [打印本页]

作者
Author: chemicalchange 时间: 2022-12-8 01:17
标题: 关于benchmark和泛函选择问题求助
本帖最后由 chemicalchange 于 2022-12-8 01:44 编辑

各位老师好，我有一个普通有机体系（含C、H、O、N、S，存在多个组分间氢键，反应涉及这些氢键相关的质子转移，最多含93个原子）的机理研究，我用的程序是ORCA 5.0.3，结构在B3LYP-D4-CPCM(toluene)/def2-SVP下进行了优化。这个机理中存在两条竞争的途径，得到了两种产物，在此问题中产物的比例相应于两个过渡态的能量差。实验得到相应两个产物比例为2.5:1，反应在110 ℃下进行，这一比例相应活化自由能差约为0.7 kcal/mol，这一数值较小，用DFT计算这样的差异比较tricky，于是我在高精度单点计算结合def2-TZVPP基组和SMD(toluene)隐式溶剂模型对比了几个不同的泛函以求寻找较符合实验事实的结果（热力学校正量由Shermo 2.3.4获得，低频采用quasi-RRHO做了校正，用了自己测试的ZPE校正因子，SCF收敛限都是tightSCF，积分格点默认），结果如下：
ωB97M-V (2.4 kcal/mol)
M06-2X (5.3 kcal/mol，这里用了defGrid3)
ωB97X-D4 (2.0 kcal/mol)
PWPB95-D4(6.4 kcal/mol)
revDSD-PBEP86-D4 (5.6 kcal/mol)

可见，我所了解的一些比较合适的计算有机体系能量的泛函在这一体系中并不理想，普遍偏高，体系中有可能带局部负电荷的地方加了弥散（ma-def2-TZVPP）对能量差几乎没有改善，换成CPCM溶剂模型基本会缓和0.5 kcal/mol，仍然不够理想。由于使用orca，没有刻意去尝试一些pople型的基组。此时我尝试了三个研究这类问题不那么有优势的泛函，结果如下：
B3LYP-D4 (3.8 kcal/mol，我知道B3LYP算能量不行，但是在优化级别下的能量差给了我一点点希望好像可能与实验比较相符，但是结合大基组后发现并不可靠)
PBE0-D4 (0.6 kcal/mol，惊奇地发现PBE0-D4在这个体系下给出的结果非常理想，但是一些测试里它算有机体系能量还是比较平庸的)
r2SCAN0-D4 (1.5 kcal/mol，看到今年JCP上Grimme的文章测试指出他弄的这个比PBE0-D4算能量问题更好，就试了一下，照搬了Grimme文章SI给的关键词，只改了基组，积分格点应该也是defGrid3)

课题组是做有机合成的，我的计算也都是跟sob老师培训学的副业，我还是个本科生，我没有这方面的投稿经验。想问一下各位老师，此时我有几个抉择：
1. 就采用PBE0-D4-SMD(toluene)/def2-TZVPP算能量（其他还有很多个中间体和反应过程），足够充分。
2. 用r2SCAN0-D4-CPCM(toluene)/def2-TZVPP算能量，由于换了溶剂模型，这个时候算出来是1.0 kcal/mol，也可以接受，但挺像凑数据的。
3. 有测试表明PBE0-DH算势垒挺不错的，我测试了一下PBE0-DH-SMD(toluene)/def2-TZVPP得到的结果是0.4 kcal/mol，我或许可以通过这个数据去印证PBE0-D4的合理性？（虽然看起来有点狡猾）
4. 既然PBE0-DH行就都用它算得了，虽然PBE0-DH算有机反应能平庸了一些，但如果不影响一些结论也无所谓。

我知道DFT想达到这个数值的精度是很困难的，这里面很多合成结合计算的文章只是在凑数据，再做杂七杂八的测试可能意义也不太大，然后投稿也遇不到什么真的很懂计算的审稿，但我还是希望能尽可能的good result due to correct reason，想请问一下各位老师上面这四种方式站在做计算的角度有什么倾向，或者有其他更好的建议吗？谢谢。

作者
Author: sobereva 时间: 2022-12-8 04:27
有条件的话，上DLPNO-CCSD(T)/def2-TZVPP，对于这样大小的体系（还可以考虑恰当简化模型来节约耗时，比如不重要区域的甲基替换成氢，等等），有个比较好的双路服务器还是能算得动的。能算出来也就没什么可质疑的了，若和实验不符那也应该再从其它角度去考虑导致差异的原因。

PBE0-DH结果恰好较满足期望，极有可能只是误差巧合抵消而已或者偶然性的运气。毕竟大规模测试没体现出比wB97M-V有什么优势

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

作者
Author: chemicalchange 时间: 2022-12-8 13:13

sobereva 发表于 2022-12-8 04:27
有条件的话，上DLPNO-CCSD(T)/def2-TZVPP，对于这样大小的体系（还可以考虑恰当简化模型来节约耗时，比如不 ...

好的老师，我尝试一下看看能不能算得动。如果算不起DLPNO-CCSD(T)，因为大规模测试表现较好的且ORCA能直接使用的泛函似乎都没什么希望给出一个与实验数据比较相符的结果。PBE0也不是一个很小众的泛函，我基于误差抵消的角度在此选用PBE0-D4算能量这样的做法能够被接受吗

作者
Author: wzkchem5 时间: 2022-12-8 15:44
有没有做构象搜索？有没有考虑隧穿效应？

作者
Author: sobereva 时间: 2022-12-8 16:29

chemicalchange 发表于 2022-12-8 13:13
好的老师，我尝试一下看看能不能算得动。如果算不起DLPNO-CCSD(T)，因为大规模测试表现较好的且ORCA能直 ...

不能
众所周知PBE0算有机体系热力学量早就过时了。而DFT-D3/D4对势垒计算精度也没任何改进。不适当的方法即便算出来好的结果也只是巧合或者误差抵消罢了，什么也说明不了。

作者
Author: gauss98 时间: 2022-12-8 17:15
个人觉得
还是在tzvp下优化一下并做频率
另外
这么大的体系
随便一个甲基取向都会有接近一个千卡的能量差距

作者
Author: chemicalchange 时间: 2022-12-8 17:34

wzkchem5 发表于 2022-12-8 15:44
有没有做构象搜索？有没有考虑隧穿效应？

这是个相当刚性的结构，做过构象搜索，没有得到其他相关的构象。隧穿确实值得考虑，这里高能量的一个过渡态是一个对碳正离子的β位质子的攫取，虚频有五六百i，我尝试一下，谢谢。

作者
Author: chemicalchange 时间: 2022-12-8 18:54

sobereva 发表于 2022-12-8 16:29
不能
众所周知PBE0算有机体系热力学量早就过时了。而DFT-D3/D4对势垒计算精度也没任何改进。不适当的方 ...

老师，我确实算不动DLPNO-CCSD(T)，考虑这个数值也确实都在各种级别的误差窗口之外了，我也只是试图在做一些误差抵消的尝试，这个问题的核心就是反应的选择性，差太远了很难讨论，我可以说是根据测试的结果，我尝试了wB97M-V不理想于是我退而求其次测试了PBE0-DH然后得到了理想的结果，这样勉强说得过去吗。

作者
Author: scf 时间: 2022-12-9 02:11
本帖最后由 scf 于 2022-12-9 02:12 编辑

问下实验上能在多大程度上确定这个反应是热力学控制,动力学控制,还是都有的？

作者
Author: sobereva 时间: 2022-12-9 05:27

chemicalchange 发表于 2022-12-8 18:54
老师，我确实算不动DLPNO-CCSD(T)，考虑这个数值也确实都在各种级别的误差窗口之外了，我也只是试图在做 ...

说不过去

恰当简化体系到能算得动的程度，也可以恰当用混合基组。或者找机子好的人帮你算。

不把计算级别问题弄确切了，就老得惦记这部分的显著误差的可能性，在讨论其它导致误差的原因时也不放心。而且PBE0-DH本来也不是流行的泛函，认可程度还很有限，到时候大概率审稿人还会让你再用别的方法验证，更何况你已经发现其它靠谱的泛函与之结果差异不小，到底PBE0-DH对此反应准不准自己都心里没数。

作者
Author: gauss98 时间: 2022-12-9 10:39

chemicalchange 发表于 2022-12-8 18:54
老师，我确实算不动DLPNO-CCSD(T)，考虑这个数值也确实都在各种级别的误差窗口之外了，我也只是试图在做 ...

在这个精度上想要定量的跟实验值吻合，不光是单点计算的问题了
下面各种因素都有0.5-1kcal左右甚至以上的误差（相对值），对你选择性都有显著影响。

频率计算的基组 svp-tzvp
在（不同）溶剂模型下算频率
简谐近似的处理方法

溶剂模型的方法（高斯默认，SMD直接算，SMD用站长推荐的M052X算）

有了这么多因素可以调节

其实，只要心里认为定性准确
在算不动dlpno-ccsd（t）的情况下（其实tzvpp、qzvpp影响都很大）情况下，找个公认的泛函较大的基组算个单点，根据上面的方法就可以凑个相对定量准确的数据写文章了
完全靠计算来肯定定量准确是不可能的。ccsd（t）也不行

作者
Author: chemicalchange 时间: 2022-12-9 11:08

gauss98 发表于 2022-12-9 10:39
在这个精度上想要定量的跟实验值吻合，不光是单点计算的问题了
下面各种因素都有0.5-1kcal左右甚至以上 ...

情况是这样的，二三十个中间体过渡态都已经这么优化了，调单点和溶剂模型目前是最划算的，溶剂模型其实前面说了换成CPCM会缓和一些但是原理上又欠佳。只是一些数据和相关的讨论希望看起来自然一点，毕竟有的反应说差了2个kcal/mol，就几乎是得到单一的产物了，这里要再又谈2个kcal/mol其实也没啥有点牵强。

作者
Author: chemicalchange 时间: 2022-12-9 11:12

scf 发表于 2022-12-9 02:11
问下实验上能在多大程度上确定这个反应是热力学控制,动力学控制,还是都有的？

实验上只能确定这个机理大概的历程是什么样子的以及分子结构中的一些官能团保护与否和氧化态高低对反应结果会有什么样的影响。这个选择性到底是哪里来的怎么控制的，我们希望通过计算给出一个解释。这个110度是甲苯回流的温度，低了不反应的，要再高就换溶剂了就不是一个体系了。

作者
Author: scf 时间: 2022-12-9 11:33
如果是热力学控制的，并不对应不同路径的过渡态能量差

作者
Author: chemicalchange 时间: 2022-12-9 11:47

scf 发表于 2022-12-9 11:33
如果是热力学控制的，并不对应不同路径的过渡态能量差

应该不是热力学控制的，很多中间过程都是不可逆的。分离到单一产物后在相同体系里也不会发生平衡的转化。至于有没有curtin-hammett那种情况确实不太清楚，是想通过计算来讨论的。现在遇到的难题是。如果是2 kcal/mol的差异，其实是可以通过curtin-hammett讨论的（昨天跑IRC反应这个两个过渡态对应的起始结构并不一致，这两个起始结构间转化的过渡态还没有找到），如果是0.6 kcal/mol是可以通过动力学控制讨论的，如果是5-6 kcal/mol是没法讨论的。所以才特别纠结这个过程。

作者
Author: scf 时间: 2022-12-9 11:55
本帖最后由 scf 于 2022-12-9 12:12 编辑

体系中有多参考态特征吗？比如symmetry breaking Hartree Fock的S^2。

分离到单一产物后在相同体系里也不会发生平衡的转化。

分离产物是在常温下吗？如果110°C反应产生了热力学控制，常温下不发生转换未必说明110°C时没热力学参与

作者
Author: gauss98 时间: 2022-12-9 12:16

chemicalchange 发表于 2022-12-9 11:08
情况是这样的，二三十个中间体过渡态都已经这么优化了，调单点和溶剂模型目前是最划算的，溶剂模型其实前 ...

其实在已知各个结构的情况下，重算一下也没多大工作量
实在不愿意重算，调调谐振近似方法也是可以的，谐振，截断100近似，截断50近似，插值近似都差不少
这些方法没有绝对的对错

说句实话，就是好看而已

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