CP2k中BSSE的使用

丁越

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分别代表了EA和EB；BSSE-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. 
   
7. raw complexation energy : -10.61861775 kcal/mol
   
8. BSSE energy : 0.001422177
   
9. 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. 
   
7. BSSE energy:  0.00061494943495
   
8. corrected binding energy: -9.842965 kcal/mol

复制代码

可见两者计算的结合能数据符合极好。

Aridea

涨姿势

Aridea

大佬，后面如果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?这解释不通啊，还是不太理解

丁越

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 （第二个单体）

Aridea

丁越 发表于 2021-10-17 19:40
这里面CP-corrected Totalenergy = EAB+EBSSE计算出来的能量已经是包含校正BSSE的能量了，结合能显然是这 ...

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

丁越

Aridea 发表于 2021-10-17 20:40
大佬，问题是如果这么计算，最后的相互作用数值正好等于0.004686，这个数不仅是正的，而且和EBSSE一样大 ...

重新更正了内容，BSSE-free interaction energy就代表了进行BSSE校正后的A、B片段之间的结合能

Aridea

丁越 发表于 2021-10-17 22:15
不用管上面的计算数据，这是直接复制CP2k官网文件的，它的文件基本是用分子力场计算的，计算原理明白就好 ...

soga, 多谢大佬耐心解答

Aridea

丁越 发表于 2021-10-17 22:15
不用管上面的计算数据，这是直接复制CP2k官网文件的，它的文件基本是用分子力场计算的，计算原理明白就好 ...

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

丁越

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）

Aridea

丁越 发表于 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这种。不知道为啥

丁越

Aridea 发表于 2021-10-19 17:54
如果是片段和复合结构的电荷不一致，自旋多重度都为1。例如一个片段是0 1另一个是2 1，总复合结构是2 1. ...

&DFT
    CHARGE    0 #Net charge
    MULTIPLICITY  1 #Spin multiplicity
这一部分是必须要设置的啊，就和gaussian一样。

Aridea

丁越 发表于 2021-10-19 19:24
&DFT
    CHARGE    0 #Net charge
    MULTIPLICITY  1 #Spin multiplicity

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

丁越

Aridea 发表于 2021-10-19 20:57
嗯嗯，大佬这里DFT部分给01是举例子嘛，如果换成我举的例子复合结构的电荷和自旋多重度 21 这里&DFT部分 ...

是的

Aridea

丁越 发表于 2021-10-19 22:37
是的

好的，多谢多谢~

丁越

2022.3.4：重新更正了本文内容