计算化学公社

 找回密码 Forget password
 注册 Register
Views: 5458|回复 Reply: 14

[CP2K] CP2k中BSSE的使用

[复制链接 Copy URL]

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

发表于 Post on 2021-9-15 14:17:21 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 丁越 于 2022-3-6 12:10 编辑

CP2K中BSSE的使用方法
一、前言
  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二聚体:
  1.     &KIND O_ghost
  2.       ELEMENT O
  3.       BASIS_SET TZVP-MOLOPT-GTH
  4.       POTENTIAL GTH-PBE
  5.       GHOST
  6.     &END KIND
  7.     &KIND H_ghost
  8.       ELEMENT H
  9.       BASIS_SET TZVP-MOLOPT-GTH
  10.       POTENTIAL GTH-PBE
  11.       GHOST
  12.     &END KIND
  13.     &KIND F_ghost
  14.       ELEMENT F
  15.       BASIS_SET TZVP-MOLOPT-GTH
  16.       POTENTIAL GTH-PBE
  17.       GHOST
  18.     &END KIND
  19.     &KIND O
  20.       ELEMENT O
  21.       BASIS_SET TZVP-MOLOPT-GTH
  22.       POTENTIAL GTH-PBE
  23.     &END KIND
  24.     &KIND H
  25.       ELEMENT H
  26.       BASIS_SET TZVP-MOLOPT-GTH
  27.       POTENTIAL GTH-PBE
  28.     &END KIND
  29.     &KIND F
  30.       ELEMENT F
  31.       BASIS_SET TZVP-MOLOPT-GTH
  32.       POTENTIAL GTH-PBE
  33.     &END KIND
复制代码
第三部分设置&BSSE,如下所示:
  1. &BSSE
  2.           &FRAGMENT #设置第一个片段H2O
  3.                   LIST 1..3 #H2O的原子序号,注意这里不能写为1-3
  4.           &END FRAGMENT
  5.           &FRAGMENT #设置第二个片段HF
  6.                   LIST 4 5
  7.           &END FRAGMENT
  8.           &CONFIGURATION
  9.                   GLB_CONF 1 1 #该关键词指定包含的片段,如1 0指只包含片段1;0 1指只包含片段2;1 1指片同时包含段1 和2.
  10.                   SUB_CONF 1 1 #该关键词指定片段的类型(1代表实片段,0代表鬼原子片段),当前1 1就代表H2O..HF的实片段。
  11.                   CHARGE 0
  12.                   MULTIPLICITY 1
  13.           &END CONFIGURATION
  14.           &CONFIGURATION
  15.                   GLB_CONF 1 0 #当前代表H2O的实片段
  16.                   SUB_CONF 1 0
  17.                   CHARGE 0
  18.                   MULTIPLICITY 1
  19.           &END CONFIGURATION
  20.           &CONFIGURATION
  21.                   GLB_CONF 0 1 #这里表示HF的实片段
  22.                   SUB_CONF 0 1
  23.                   CHARGE 0
  24.                   MULTIPLICITY 1
  25.           &END CONFIGURATION
  26.           &CONFIGURATION
  27.                   GLB_CONF 1 1 #这里将HF设置为鬼原子
  28.                   SUB_CONF 1 0
  29.                   CHARGE 0
  30.                   MULTIPLICITY 1
  31.           &END CONFIGURATION
  32.           &CONFIGURATION
  33.                   GLB_CONF 1 1 #同上,将H2O设置为鬼原子
  34.                   SUB_CONF 0 1
  35.                   CHARGE 0
  36.                   MULTIPLICITY 1
  37.           &END CONFIGURATION
  38.   &END BSSE
  39.                  
复制代码

CP2K在计算过程中会进行5次能量计算,分别是:

  
Calculations
  
E
与&CONFIGURATION部分对应关系
1
EA
1 0 1 0
2
EB
0 1 0 1
3
EAB(A)
1 1 1 0
4
EAB(B)
1 1 0 1
5
EAB
1 1 1 1
在最后计算完成后,会输出
  1. -------------------------------------------------------------------------------
  2. -                                                                             -
  3. -                                 BSSE RESULTS                                -
  4. -                                                                             -
  5. -                 CP-corrected Total energy:      -42.083993                  -
  6. -                                                                             -
  7. -                       1-body contribution:      -17.222002                  -
  8. -                       1-body contribution:      -24.846305                  -
  9. -                                                                             -
  10. -                       2-body contribution:       -0.015687                  -
  11. -                 BSSE-free interaction energy:       -0.015687               -
复制代码
其中,CP-corrected Total energy指EAB+E_BSSE两个1-body contribution分别代表了EAEBBSSE-free interaction energy就代表了进行BSSE校正后的A、B片段之间的结合能,这里所有单位都是a.u.

下面将与高斯计算结果进行对比:
  Gaussian:简单测试为目的,结构优化和单点都使用了b3lyp/def2tzvp计算级别,实际计算中,单点能计算对于这么简单的分子得用CCSD(T)。结果如下:
  1. H2O_HF: -176.970414
  2. H2O :  -76.4633294
  3. HF_bq : -76.4641418
  4. HF :  -100.4901625
  5. H2O_bq :  -100.490772277

  6. raw complexation energy : -10.61861775 kcal/mol
  7. BSSE energy : 0.001422177
  8. corrected complexation energy :   -9.73 kcal/mol
复制代码

  CP2K:直接利用高斯优化完的结构,将盒子边长大小都设置为30埃,所有元素均用PBE、TZVP-MOLOPT-GTH级别,结果如下:
  1. H2O:  -17.22200178914723
  2. HF:  -24.84630455020665
  3. HF_ghost:  -17.22214957254267
  4. H2O_ghost:  -24.84677171624616
  5. H2O_HF:  -42.08460803369478

  6. BSSE energy:  0.00061494943495
  7. corrected binding energy: -9.842965 kcal/mol
复制代码
可见两者计算的结合能数据符合极好。

H2O_HF.inp

4.05 KB, 下载次数 Times of downloads: 17

评分 Rate

参与人数
Participants 5
eV +30 收起 理由
Reason
tkaray + 5 GJ!
djjj148 + 5 好物!
Aridea + 5 GJ!
ChemG + 5 好物!
sobereva + 10

查看全部评分 View all ratings

自由发挥,野蛮生长

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-9-30 16:09:00 | 显示全部楼层 Show all
涨姿势

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-17 19:12:47 | 显示全部楼层 Show all
大佬,后面如果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?这解释不通啊,还是不太理解

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-10-17 19:40:14 | 显示全部楼层 Show all
Aridea 发表于 2021-10-17 19:12
大佬,后面如果CP-corrected Totalenergy = EAB+EBSSE = -34.320623;EA=1-bodycontribution=-17.162695 ;  ...

这里面CP-corrected Totalenergy = EAB+EBSSE计算出来的能量已经是包含校正BSSE的能量了,结合能显然是这个CP-corrected Totalenergy - 1-bodycontribution(第一个单体)- 1-bodycontribution (第二个单体)
自由发挥,野蛮生长

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-17 20:40:00 | 显示全部楼层 Show all
丁越 发表于 2021-10-17 19:40
这里面CP-corrected Totalenergy = EAB+EBSSE计算出来的能量已经是包含校正BSSE的能量了,结合能显然是这 ...

大佬,问题是如果这么计算,最后的相互作用数值正好等于0.004686,这个数不仅是正的,而且和EBSSE一样大小,这意味着EAB=EA+EB。
从来没见过复合结构和两个单体能量合一样的弱相互作用。所以我觉得输出的解读对应E xxx这里应该是有什么不对劲的地方,这也是不能理解的唯一原因。

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-10-17 22:15:08 | 显示全部楼层 Show all
本帖最后由 丁越 于 2022-3-4 16:11 编辑
Aridea 发表于 2021-10-17 20:40
大佬,问题是如果这么计算,最后的相互作用数值正好等于0.004686,这个数不仅是正的,而且和EBSSE一样大 ...

重新更正了内容,BSSE-free interaction energy就代表了进行BSSE校正后的A、B片段之间的结合能
自由发挥,野蛮生长

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-17 22:30:15 | 显示全部楼层 Show all
丁越 发表于 2021-10-17 22:15
不用管上面的计算数据,这是直接复制CP2k官网文件的,它的文件基本是用分子力场计算的,计算原理明白就好 ...

soga, 多谢大佬耐心解答

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-18 20:31:56 | 显示全部楼层 Show all
本帖最后由 Aridea 于 2021-10-18 22:05 编辑
丁越 发表于 2021-10-17 22:15
不用管上面的计算数据,这是直接复制CP2k官网文件的,它的文件基本是用分子力场计算的,计算原理明白就好 ...

大佬,这个这里bsse里加了charge 和 multiplicity. 后面DFT部分加和不加是不是一样的呀

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-10-19 12:51:00 | 显示全部楼层 Show all
本帖最后由 丁越 于 2022-3-4 16:01 编辑
Aridea 发表于 2021-10-18 20:31
大佬,这个这里bsse里加了charge 和 multiplicity. 后面DFT部分加和不加是不是一样的呀

不可以的,&DFT后面如果你不指定电荷和自旋多重度的话CHARGE和MULTIPLICITY的默认值都是0,显然是不正确的。(https://manual.cp2k.org/cp2k-8_2 ... FT.html#list_CHARGE

自由发挥,野蛮生长

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-19 17:54:28 | 显示全部楼层 Show all
丁越 发表于 2021-10-19 12:51
不可以的,&DFT后面如果你不指定电荷和自旋多重度的话CHARGE和MULTIPLICITY的默认值都是0,显然是不正确 ...

如果是片段和复合结构的电荷不一致,自旋多重度都为1。例如一个片段是0 1另一个是2 1,总复合结构是2 1.使用BSSE的时候在&BSSE的部分分别安装帖子定义了2个fragment,和5个configuration的电荷和自旋多重度。最后&DFT的部分里面的电荷和自旋多重度是不是还要给个复合结构的电荷和自旋多重度2 1。然后这一个文件计算所有弱相互作用需要的能量。
我按照这么操作有一批计算算出来的BSSE很离谱,全是负的而且是-0.8,-0.4这种。不知道为啥

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-10-19 19:24:55 | 显示全部楼层 Show all
Aridea 发表于 2021-10-19 17:54
如果是片段和复合结构的电荷不一致,自旋多重度都为1。例如一个片段是0 1另一个是2 1,总复合结构是2 1. ...

&DFT
    CHARGE    0 #Net charge
    MULTIPLICITY  1 #Spin multiplicity
这一部分是必须要设置的啊,就和gaussian一样。
自由发挥,野蛮生长

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-19 20:57:41 | 显示全部楼层 Show all
丁越 发表于 2021-10-19 19:24
&DFT
    CHARGE    0 #Net charge
    MULTIPLICITY  1 #Spin multiplicity

嗯嗯,大佬这里DFT部分给01是举例子嘛,如果换成我举的例子复合结构的电荷和自旋多重度 21 这里&DFT部分也给2 1 对吧

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-10-19 22:37:10 | 显示全部楼层 Show all
Aridea 发表于 2021-10-19 20:57
嗯嗯,大佬这里DFT部分给01是举例子嘛,如果换成我举的例子复合结构的电荷和自旋多重度 21 这里&DFT部分 ...

是的
自由发挥,野蛮生长

132

帖子

0

威望

701

eV
积分
833

Level 4 (黑子)

发表于 Post on 2021-10-19 23:05:02 | 显示全部楼层 Show all

好的,多谢多谢~

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-4 16:04:13 | 显示全部楼层 Show all
2022.3.4:重新更正了本文内容
自由发挥,野蛮生长

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-2 22:05 , Processed in 0.432947 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list