计算化学公社
标题: CP2k中BSSE的使用 [打印本页]
作者Author: 丁越 时间: 2021-9-15 14:17
标题: CP2k中BSSE的使用
本帖最后由 丁越 于 2022-3-6 12:10 编辑
一、前言
BSSE来源于使用高斯型基函数的基组计算AB复合物时，A的基函数会延伸到B空间，同样的，B的基函数也会延伸到A空间，这样在计算AB复合物时相当于所用到的基组要大于单独计算A和B所用到的基组，导致计算AB的能量更负。而对于平面波基组，就不会有BSSE的问题。在cp2k中使用高斯混合平面波基组的方法，因此涉及到计算弱相互作用能的时候需要采用CP（counterpoise）方法进行BSSE校正,校正公式如下：
E_BSSE=[E(A)-EAB(A)] + [E(B)-EAB(B)]
E(A)-EAB(A)表示单独计算片段A的能量减去在AB基函数下片段A的能量，对片段B亦如此。EAB(A)的计算就要用到鬼原子，当下即将B设置为鬼原子。
二、cp2k中BSSE的设置
输入文件中需要进行三部分的设置，首先将RUN_TYPE 设置为BSSE。第二部分是&KIND部分额外添加AB中所有元素的鬼原子（鬼原子（ghost atoms）：指某个位置只带有基函数但是不带有原子电荷和自旋多重度），例如对于H2O..HF二聚体：
- &KIND O_ghost
- ELEMENT O
- BASIS_SET TZVP-MOLOPT-GTH
- POTENTIAL GTH-PBE
- GHOST
- &END KIND
- &KIND H_ghost
- ELEMENT H
- BASIS_SET TZVP-MOLOPT-GTH
- POTENTIAL GTH-PBE
- GHOST
- &END KIND
- &KIND F_ghost
- ELEMENT F
- BASIS_SET TZVP-MOLOPT-GTH
- POTENTIAL GTH-PBE
- GHOST
- &END KIND
- &KIND O
- ELEMENT O
- BASIS_SET TZVP-MOLOPT-GTH
- POTENTIAL GTH-PBE
- &END KIND
- &KIND H
- ELEMENT H
- BASIS_SET TZVP-MOLOPT-GTH
- POTENTIAL GTH-PBE
- &END KIND
- &KIND F
- ELEMENT F
- BASIS_SET TZVP-MOLOPT-GTH
- POTENTIAL GTH-PBE
- &END KIND
第三部分设置&BSSE，如下所示：
- &BSSE
- &FRAGMENT #设置第一个片段H2O
- LIST 1..3 #H2O的原子序号，注意这里不能写为1-3
- &END FRAGMENT
- &FRAGMENT #设置第二个片段HF
- LIST 4 5
- &END FRAGMENT
- &CONFIGURATION
- GLB_CONF 1 1 #该关键词指定包含的片段，如1 0指只包含片段1；0 1指只包含片段2；1 1指片同时包含段1 和2.
- SUB_CONF 1 1 #该关键词指定片段的类型（1代表实片段，0代表鬼原子片段），当前1 1就代表H2O..HF的实片段。
- CHARGE 0
- MULTIPLICITY 1
- &END CONFIGURATION
- &CONFIGURATION
- GLB_CONF 1 0 #当前代表H2O的实片段
- SUB_CONF 1 0
- CHARGE 0
- MULTIPLICITY 1
- &END CONFIGURATION
- &CONFIGURATION
- GLB_CONF 0 1 #这里表示HF的实片段
- SUB_CONF 0 1
- CHARGE 0
- MULTIPLICITY 1
- &END CONFIGURATION
- &CONFIGURATION
- GLB_CONF 1 1 #这里将HF设置为鬼原子
- SUB_CONF 1 0
- CHARGE 0
- MULTIPLICITY 1
- &END CONFIGURATION
- &CONFIGURATION
- GLB_CONF 1 1 #同上，将H2O设置为鬼原子
- SUB_CONF 0 1
- CHARGE 0
- MULTIPLICITY 1
- &END CONFIGURATION
- &END BSSE
-
CP2K在计算过程中会进行5次能量计算，分别是：
在最后计算完成后，会输出
- -------------------------------------------------------------------------------
- - -
- - BSSE RESULTS -
- - -
- - CP-corrected Total energy: -42.083993 -
- - -
- - 1-body contribution: -17.222002 -
- - 1-body contribution: -24.846305 -
- - -
- - 2-body contribution: -0.015687 -
- - BSSE-free interaction energy: -0.015687 -
其中，CP-corrected Total energy指EAB+E_BSSE；两个1-body contribution分别代表了EA和EB；BSSE-free interaction energy就代表了进行BSSE校正后的A、B片段之间的结合能，这里所有单位都是a.u.
下面将与高斯计算结果进行对比：
Gaussian：简单测试为目的，结构优化和单点都使用了b3lyp/def2tzvp计算级别，实际计算中，单点能计算对于这么简单的分子得用CCSD(T)。结果如下：
- H2O_HF: -176.970414
- H2O : -76.4633294
- HF_bq : -76.4641418
- HF : -100.4901625
- H2O_bq : -100.490772277
- raw complexation energy : -10.61861775 kcal/mol
- BSSE energy : 0.001422177
- corrected complexation energy : -9.73 kcal/mol
CP2K：直接利用高斯优化完的结构，将盒子边长大小都设置为30埃，所有元素均用PBE、TZVP-MOLOPT-GTH级别，结果如下：
- H2O: -17.22200178914723
- HF: -24.84630455020665
- HF_ghost: -17.22214957254267
- H2O_ghost: -24.84677171624616
- H2O_HF: -42.08460803369478
- BSSE energy: 0.00061494943495
- corrected binding energy: -9.842965 kcal/mol
可见两者计算的结合能数据符合极好。
作者Author: Aridea 时间: 2021-9-30 16:09
涨姿势
作者Author: Aridea 时间: 2021-10-17 19:12
大佬，后面如果CP-corrected Totalenergy = EAB+EBSSE = -34.320623；EA=1-bodycontribution=-17.162695 ; EB=1-bodycontribution=-17.162614,
最后AB的相互作用能怎么算呀，是不是Delta E = EAB +EBSSE - EA -EB = 0.004686?这解释不通啊，还是不太理解
作者Author: 丁越 时间: 2021-10-17 19:40
这里面CP-corrected Totalenergy = EAB+EBSSE计算出来的能量已经是包含校正BSSE的能量了，结合能显然是这个CP-corrected Totalenergy - 1-bodycontribution（第一个单体）- 1-bodycontribution （第二个单体）
作者Author: Aridea 时间: 2021-10-17 20:40
大佬，问题是如果这么计算，最后的相互作用数值正好等于0.004686，这个数不仅是正的，而且和EBSSE一样大小，这意味着EAB=EA+EB。
从来没见过复合结构和两个单体能量合一样的弱相互作用。所以我觉得输出的解读对应E xxx这里应该是有什么不对劲的地方，这也是不能理解的唯一原因。
作者Author: 丁越 时间: 2021-10-17 22:15
本帖最后由 丁越 于 2022-3-4 16:11 编辑
重新更正了内容，BSSE-free interaction energy就代表了进行BSSE校正后的A、B片段之间的结合能
作者Author: Aridea 时间: 2021-10-17 22:30
soga, 多谢大佬耐心解答
作者Author: Aridea 时间: 2021-10-18 20:31
本帖最后由 Aridea 于 2021-10-18 22:05 编辑
大佬，这个这里bsse里加了charge 和 multiplicity. 后面DFT部分加和不加是不是一样的呀
作者Author: 丁越 时间: 2021-10-19 12:51
本帖最后由 丁越 于 2022-3-4 16:01 编辑
不可以的，&DFT后面如果你不指定电荷和自旋多重度的话CHARGE和MULTIPLICITY的默认值都是0，显然是不正确的。（https://manual.cp2k.org/cp2k-8_2 ... FT.html#list_CHARGE）
作者Author: Aridea 时间: 2021-10-19 17:54
如果是片段和复合结构的电荷不一致，自旋多重度都为1。例如一个片段是0 1另一个是2 1，总复合结构是2 1.使用BSSE的时候在&BSSE的部分分别安装帖子定义了2个fragment,和5个configuration的电荷和自旋多重度。最后&DFT的部分里面的电荷和自旋多重度是不是还要给个复合结构的电荷和自旋多重度2 1。然后这一个文件计算所有弱相互作用需要的能量。
我按照这么操作有一批计算算出来的BSSE很离谱，全是负的而且是-0.8，-0.4这种。不知道为啥
作者Author: 丁越 时间: 2021-10-19 19:24
&DFT
CHARGE 0 #Net charge
MULTIPLICITY 1 #Spin multiplicity
这一部分是必须要设置的啊，就和gaussian一样。
作者Author: Aridea 时间: 2021-10-19 20:57
嗯嗯，大佬这里DFT部分给01是举例子嘛，如果换成我举的例子复合结构的电荷和自旋多重度 21 这里&DFT部分也给2 1 对吧
作者Author: 丁越 时间: 2021-10-19 22:37
是的
作者Author: Aridea 时间: 2021-10-19 23:05
好的，多谢多谢~
作者Author: 丁越 时间: 2022-3-4 16:04
2022.3.4：重新更正了本文内容
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |