计算化学公社

标题: 关于团簇结构与有机分子间吸附能计算时的BSSE处理问题 [打印本页]

作者
Author: monoyu 时间: 2020-3-14 14:01
标题: 关于团簇结构与有机分子间吸附能计算时的BSSE处理问题
本帖最后由 monoyu 于 2020-3-14 14:02 编辑

各位老师，请教下关于石墨烯团簇与单个苯酚分子之间吸附能计算以及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) 是否正确和合理呢？

作者
Author: 卡开发发 时间: 2020-3-14 16:53
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也并不太适合做小体系的高精度计算，毕竟所支持的方法不够好。

作者
Author: sobereva 时间: 2020-3-14 18:49
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。

作者
Author: monoyu 时间: 2020-3-14 21:53

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

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

作者
Author: monoyu 时间: 2020-3-14 22:01

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使用，希望可以像您说的用更合适的工具去处理这些问题。

作者
Author: 卡开发发 时间: 2020-3-14 22:01

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

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

作者
Author: monoyu 时间: 2020-3-14 22:21

卡开发发发表于 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，再去解决此类问题。

作者
Author: 卡开发发 时间: 2020-3-14 22:53

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

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

作者
Author: telunsu 时间: 2021-8-19 11:36

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

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

作者
Author: 卡开发发 时间: 2021-8-19 13:40

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

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

作者
Author: telunsu 时间: 2021-8-19 14:13
本帖最后由 telunsu 于 2021-8-19 15:15 编辑

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

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

DMOL3 定义的是 (, 下载次数 Times of downloads: 60) ,而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吗？

作者
Author: wzkchem5 时间: 2021-8-19 15:55

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

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

作者
Author: telunsu 时间: 2021-8-19 17:52
本帖最后由 telunsu 于 2021-8-19 17:53 编辑

卡开发发发表于 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(苯酚）

作者
Author: 卡开发发 时间: 2021-8-19 19:08

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这样的等级，默认对轻元素截断半径可能给的有点小。

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

作者
Author: 喵星大佬 时间: 2021-8-19 19:12

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

来算的，

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

作者
Author: 卡开发发 时间: 2021-8-19 19:26

喵星大佬发表于 2021-8-19 19:12
我觉得出现这种问题应该先检查模型和SCF的问题

确实如此，有时候我们总认为其他人会自己多检查几遍才会到论坛上来问，不过实际遇到的低级错误是不少（比如POTCAR顺序写错这种事情就发生了n次）。

作者
Author: telunsu 时间: 2021-8-19 19:43
本帖最后由 telunsu 于 2021-8-19 19:52 编辑

卡开发发发表于 2021-8-19 19:26
确实如此，有时候我们总认为其他人会自己多检查几遍才会到论坛上来问，不过实际遇到的低级错误是不少（比 ...

实在抱歉，引起大家的误会了，我只是拿这个楼主的公式举了个例子，我实际是在用DMOL3计算MIL-101与水、呋喃的相互作用
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

还是想问下： 1，是否可以Total energy BSSE corrected = -1989.452953代替E_1+2和E_BSSE 这两项，计算总能量？

因为我有时发现Total energy BSSE corrected =E_AB -BSSE correction energy,而有时相差又有点大,体系含Cr，spin为2的时候能量最低，选pw91 是参照了别人之前在JPCC的工作
2，如果不能代替，那么Total energy BSSE corrected  代表什么？十分感谢各位老师的回答！

# Task parameters
Calculate                   optimize
Opt_energy_convergence       1.0000e-005
Opt_gradient_convergence    2.0000e-003 A
Opt_displacement_convergence  5.0000e-003 A
Opt_iterations             800
Opt_max_displacement       0.3000 A

Symmetry                   off
Max_memory                   204800
File_usage                   smart
Scf_density_convergence    1.000000e-005
Scf_charge_mixing          2.000000e-001
Scf_spin_mixing             5.000000e-001
Scf_diis                   6 pulay
Scf_iterations             300

# Electronic parameters
Spin_polarization          unrestricted
Spin                         2
Charge                      0
Basis                      dnp
basis_version                basfile_v4.4
Pseudopotential             dspp
Functional                   gga(p91)
Aux_density                octupole
Dftd                         OBS
Integration_grid             medium
Occupation                   thermal 0.0020
Cutoff_Global                4.4000 angstrom

作者
Author: 卡开发发 时间: 2021-8-19 20:10

telunsu 发表于 2021-8-19 19:43
实在抱歉，引起大家的误会了，我只是拿这个楼主的公式举了个例子，我实际是在用DMOL3计算MIL-101与水、呋 ...

按照前面给的BSSE correction energy去做应该肯定没问题，这个值和上面的A(B)和A，B(A)和B之间的关系算出来是没问题的，能确保这个就可了。Total energy BSSE corrected他是咋定义的我不太清楚。

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