关于团簇结构与有机分子间吸附能计算时的BSSE处理问题

monoyu

各位老师，请教下关于石墨烯团簇与单个苯酚分子之间吸附能计算以及BSSE的处理问题。软件：M$2018-Dmol3
处理步骤如下：1.建立石墨烯团簇模型（23个碳六元环，边缘加H原子作饱和处理，共计66个C原子，22个H原子），然后做结构优化（GGA, PBE, DFT-D, DNP基组，Global orbital cutoff 5A，use COSMO water作隐式溶剂化处理）；
2.建立单个苯酚结构模型，然后做结构优化（GGA, PBE, DFT-D, DNP基组，Global orbital cutoff 5A，use COSMO water作隐式溶剂化处理）；
3.将上述优化后的两个结构放在一起组成复合物，调整相对位置，进行结构优化（GGA, PBE, DFT-D, DNP基组，Global orbital cutoff 5A，use COSMO water作隐式溶剂化处理），得到总能量 E(复合结构）；4.按照步骤3优化后的结构进行复合物的单点能计算（考虑BSSE，GGA, PBE, DFT-D, DNP基组，Global orbital cutoff 5A），得到校正能量 E(Bsse) ；
5. 按照步骤3优化后的结构分别计算石墨烯团簇和苯酚分子的单点能（没有考虑分子内BSSE，GGA, PBE, DFT-D, DNP基组，Global orbital cutoff 5A，use COSMO water作隐式溶剂化处理），得到E(石墨烯）和E(苯酚）。

那么石墨烯团簇与单个苯酚分子之间的作用能 Eab= E(复合结构）- E(石墨烯）- E(苯酚）- E(Bsse) 是否正确和合理呢？

卡开发发

DMol3的BSSE我个人很少使用，有几个保守一些的观点不妨讨论一下：
1、按照DMol3给的BSSE公式，这里Eab的计算公式应该是没问题的，不过BSSE一般是在气相计算下面做（也就是只有第五步不考虑溶剂化），考虑溶剂化的情况结果可能会有异常，然后再把得到的E_BSSE校正到溶剂化下Eab= E(复合结构）- E(石墨烯）- E(苯酚）的结果上。
2、DMol3的数值原子基组一般BSSE是不太大的（可以参考Inada Y.; Orita H. J. Comput Chem., 29, 225 (2008).），因此不进行高精度计算的情况可能不太重要，并且DMol3也并不太适合做小体系的高精度计算，毕竟所支持的方法不够好。

sobereva

E(BSSE)是正值，应当加到E(复合结构）- E(石墨烯）- E(苯酚）上去。看此文
谈谈BSSE校正与Gaussian对它的处理
http://sobereva.com/46
算BSSE校正能时不能带溶剂模型，看上文

而且为了精确，应当用TZ级别的而非DNP。

PS：用Dmol3算这种问题远远不如用免费、精度好得多的ORCA。PBE-D2算这个太糙了，ORCA的话用PWPB95-D3(BJ)/def2-QZVP算这个精度强得多得多还花不了多少时间，Multiwfn直接就能产生输入文件。Dmol3在算弱相互作用方面已全面过时，碰见内行的审稿人恐怕都不能忍PBE-D2。

monoyu

卡开发发 发表于 2020-3-14 16:53
DMol3的BSSE我个人很少使用，有几个保守一些的观点不妨讨论一下：
1、按照DMol3给的BSSE公式，这里Eab的计 ...

谢谢卡开发发老师回复。
BSSE就是在气相计算下面做的，不支持溶剂化的情况。输出文件中得到的 “BSSE  correction energy ” 应该就是E_BSSE了吧？不过怎么BSSE  correction energy 都是负的？跟社长提到的不一样，所以有点迷惑。

monoyu

sobereva 发表于 2020-3-14 18:49
E(BSSE)是正值，应当加到E(复合结构）- E(石墨烯）- E(苯酚）上去。看此文
谈谈BSSE校正与Gaussian对它的 ...

谢谢社长的回复和指点。
BSSE计算后输出文件中得到的 “BSSE  correction energy ” 应该就是E_BSSE了吧？不过怎么BSSE  correction energy 都是负的？看了几篇文献中貌似也是Eab=E(A+B）-E(A )-E(B) - E(BSSE)的公式。所以有点迷惑。目前也在学习ORCA使用，希望可以像您说的用更合适的工具去处理这些问题。

卡开发发

monoyu 发表于 2020-3-14 21:53
谢谢卡开发发老师回复。
BSSE就是在气相计算下面做的，不支持溶剂化的情况。输出文件中得到的 “BSSE  c ...

因为DMol3的BSSE定义是反的：E(BSSE)=[E(A,AB)-E(A)]+[E(B,AB)-E(B)]。

monoyu

卡开发发 发表于 2020-3-14 22:01
因为DMol3的BSSE定义是反的：E(BSSE)=[E(A,AB)-E(A)]+[E(B,AB)-E(B)]。

谢谢卡开发发老师。十几个初猜的吸附构型，计算得到的E(BSSE）基本都在0.0005-0.00015Ha范围。
如社长的建议，还是学习下ORCA，再去解决此类问题。

卡开发发

monoyu 发表于 2020-3-14 22:21
谢谢卡开发发老师。十几个初猜的吸附构型，计算得到的E(BSSE）基本都在0.0005-0.00015Ha范围。
如社长的 ...

这个很正常，DMol3的基组设计的时候BSSE就不大，上面给了参考文献。
ORCA是个好的选择。当然DMol3或者说其他数值轨道程序并非完全不值得用，对于比较大的体系尤其是包含周期性边界的情况，效率上还是有不小的优势，应当根据你所研究的问题和预期的精度选择合适的计算程序。

telunsu

卡开发发 发表于 2020-3-14 22:53
这个很正常，DMol3的基组设计的时候BSSE就不大，上面给了参考文献。
ORCA是个好的选择。当然DMol3或者说 ...

因为Dmol3计算的BSSE是负值，那么这么计算是否是对的呢？
Eab= E(复合结构）- E(石墨烯）- E(苯酚）- E(Bsse)

卡开发发

telunsu 发表于 2021-8-19 11:36
因为Dmol3计算的BSSE是负值，那么这么计算是否是对的呢？
Eab= E(复合结构）- E(石墨烯）- E(苯酚）- E ...

按照定义应当是Eab= E(复合结构）- E(石墨烯）- E(苯酚）+ E(Bsse)。

telunsu

卡开发发 发表于 2021-8-19 13:40
按照定义应当是Eab= E(复合结构）- E(石墨烯）- E(苯酚）+ E(Bsse)。

老师的意思是需要＋一个负值吗？，还是要取这个负值的绝对值相加？因为我现在Eab很大，如果再加上一个负值，那吸附能会变得异常大

DMOL3 定义的是 ,而sobereva 在http://sobereva.com/46写道是E_BSSE，即(E_A - E_A,bAB) + (E_B - E_B,bAB)，这两个公式是互为相反数的，所以如果使用Dmol计算两个物质的吸附能，E1=10 ev，E2=110ev E1+2=100ev  BSSE=-5ev   那么吸附能BE= E1+2 - E1 -E2 +（BSSE的绝对值)=100-10-110+(5ev)=-15ev吗？

wzkchem5

telunsu 发表于 2021-8-19 07:13
老师的意思是需要＋一个负值吗？，还是要取这个负值的绝对值相加？因为我现在Eab很大，如果再加上一个负 ...

首先检查能量有没有减错，石墨烯和苯酚之间的吸附是纯粹的非共价吸附，吸附能应该在1eV以内，BSSE的绝对值也应该在1eV以内。不可能算出10eV量级的作用能。检查一下吸附前后的结构有没有优化，吸附前后石墨烯的结构是否定性相同，石墨烯是否有不饱和键，吸附前后体系波函数是否定性一致。
BSSE全称是basis set superposition error，也就是基组重叠误差。有人用BSSE指代这个误差本身，也有人指代校正掉这个误差需要做的那个校正，这两者互为相反数，这就是为什么你看到两种定义。至于该取什么符号，很好记，就是BSSE总是有利于结合的，BSSE校正总是不利于结合的，所以应该取让结合变得更不利的那个符号。

telunsu

卡开发发 发表于 2021-8-19 13:40
按照定义应当是Eab= E(复合结构）- E(石墨烯）- E(苯酚）+ E(Bsse)。

DMol3 BSSE Energy Calculations (Ha)
Total energy A(B) = -232.219935
Total energy B(A) = -1756.972062
Total energy A    = -232.205937
Total energy B    = -1756.970660
BSSE correction energy = -0.015400
Total energy BSSE corrected = -1989.452953
Counterpoise corrected interaction energy = -0.276356

DMOL 的输出如上，我在计算吸附能的时候是否可以 使用 Total energy BSSE corrected = -1989.452953代替E_1+2和E_BSSE 这两项，计算总能量？

即Eab= E(复合结构）- E(石墨烯）- E(苯酚）+ E(Bsse)=Total energy BSSE corrected--E(石墨烯）- E(苯酚）

卡开发发

telunsu 发表于 2021-8-19 14:13
老师的意思是需要＋一个负值吗？，还是要取这个负值的绝对值相加？因为我现在Eab很大，如果再加上一个负 ...

他这边BSSE correction energy是按照

E=E_A(B)-E_A+E_B(A)-E_B

来算的，
所以应该取

Eab= E(复合结构）- E(石墨烯）- E(苯酚）- BSSE correction energy

能量没减错的话，可能出现的问题是几个体系截断半径的值等导致的，首先截断半径值要设置足够，然后每个体系都按这个取，不要只是设置fine或medium这样的等级，默认对轻元素截断半径可能给的有点小。

相互作用能确实有点大，你看下最终吸附能是多少。

喵星大佬

卡开发发 发表于 2021-8-19 19:08
他这边BSSE correction energy是按照

来算的，

我觉得出现这种问题应该先检查模型和SCF的问题