【Qbics】用Qbics计算二体和多体能量分解（EDA）分析

coolrainbow

1 Qbics介绍

大家好，我是coolrainbow。我非常高兴能在这里向大家推广一下我们开发的多尺度电子结构模拟软件Qbics。Qbics是我们在深圳湾实验室开发的新一代电子结构软件，其开发重点是利用密度泛函理论和QM/MM方法进行复杂体系的能量和动力学模拟。除了通用功能外，还包括大量原创的理论。Qbics提供全部源代码和二进制文件，欢迎大家下载学习、使用！除了电子结构程序Qbics外，我们还提供了可视化程序Qbics-MolStar和分子拓扑建模软件pdbtop，欢迎大家使用！

Qbics下载（源代码与二进制文件）：Qbics — ZHANG Jun's Website documentation
Qbics安装与使用手册：Qbics Manual — Qbics Manual
Qbics已引用案例：Qbics Applications — ZHANG Jun's Website documentation
Qbics输入文件案例：Input Library — Qbics Manual
可视化程序Qbics-MolStar: Mol* Viewer
拓扑建模软件pdbtop：pdbtop — ZHANG Jun's Website documentation

Qbics已经实现的功能，其中加粗的为原创理论。还有大量功能在开发中。

能量与梯度计算引擎：

• 密度泛函理论DFT
• 半经验方法xTB
• 半经验方法NDDO(AM1, PM3等)
• 经典力场（CHARMM，AMBER）
• QMMM：DFT/MM，xTB/MM，NDDO/MM, PHO
• 激发态方法：TSO-DFT， TSO-DFT/MM，BLW-DFT, BLW-DFT/MM
• 含时密度泛函方法：TDDFT
• PCM溶剂模型（IEFPCM，DPCM，COSMO）
  

势能面探索

• 几何优化（可以包含约束）
• 过渡态搜索（可以包含约束，NEB，DIMER，QST2）
• MECP搜索
• 分子动力学（NVE， NVT， NPT）
• 增强采样分子动力学（Plumed）
• 经典自由能微扰（FEP）
• 二体与多体能量分解
  

量子动力学：

• 路径积分方法（PI-FEP）
• 非正交组态相互作用
  

Qbics的开发坚持面向用户，让用户用最简单的输入执行最复杂的功能，可直接输出波函数文件mwfn、轨迹文件xyz、光谱文件txt等配合Multiwfn、Qbics-MolStar等用于分析。我们将在论坛上陆续提供相关教程。


本文的更多详情请见：Energy Decomposition Analysis — Qbics Manual



2 二体能量分解分析（EDA）

能量分解是量子化学中一个非常有用的功能。在Qbics中，该功能可以以一个关键词eda实现。下面举一个例子，以 (H2SO4)(NH3) 为例：

1. basis
   
2.     def2-svp
   
3. end
   
4. 
   
5. scf
   
6.     charge     0
   
7.     spin2p1    1
   
8.     type       U # This is needed for EDA.
   
9. end
   
10. 
    
11. grimmedisp
    
12.     type bj       # Options: none, zero, bj.
    
13. end
    
14. 
    
15. eda
    
16.     type      tso
    
17.     frag  0 1 1-7
    
18.     frag  0 1 8-11
    
19. end
    
20. 
    
21. mol
    
22.     S                 -0.62161000   -0.11181000    0.09167200
    
23.     O                  0.03138700   -0.14316200    1.35576500
    
24.     O                 -1.82120800   -0.82692400   -0.15912300
    
25.     O                  0.38678100   -0.46814500   -1.02324300
    
26.     O                 -0.91128600    1.43385100   -0.17168300
    
27.     H                  1.34667100   -0.25497400   -0.71202200
    
28.     H                 -1.58729700    1.50230600   -0.85583700
    
29.     N                  2.78190500    0.03915200   -0.05091700
    
30.     H                  3.20716100    0.92239700   -0.30232700
    
31.     H                  3.46881800   -0.69082400   -0.19060100
    
32.     H                  2.55167500    0.07103700    0.93671700
    
33. end
    
34. 
    
35. task
    
36.     eda b3lyp
    
37. end

复制代码

输入文件的大部分内容比较好懂：

• basis...end: 基组。可以用任意常见基组，也可以自定义。
• scf...end：定义SCF计算。charge和spin2p1分别是电荷和自旋多重度。type U表示使用unrestricted SCF。
• grimmedisp...end：使用Grimme色散校正。
• mol...end：分子坐标。
• eda...edn：定义能量分解的计算。
• type tso: EDA的类型。有两种：
  

• TSO: 使用target state optimization (TSO)方法进行能量分解（推荐）。该方法主要由coolrainbow开发和实现。请引用：J. Chem. Theory Comput. 2023, 19, 1777
• GKS：使用generalized Kohn-Sham (GKS) 进行能量分解。该方法由厦大苏培峰老师开发，由唐振实现。请引用：J. Phys. Chem. A 2014, 118, 2531；WIREs Comput. Mol. Sci. 2020, 10, e1460；Phys. Chem. Chem. Phys. 2024, 26, 17549
• frag 0 1 1-7 将原子1-7设为一个碎片，电荷和自选多重度为0和1。
• frag 0 1 8-11 将原子8-11设为一个碎片，电荷和自选多重度为0和1。
  

碎片的定义见图片：



将该文件存为eda.inp，用Qbics执行(-n 24 表示用24核并行)：

1. qbics eda.inp -n 24 > eda.out

复制代码

结果可在eda.out找到：

1. WITHOUT BSSE correction:
   
2. Electrostatic interaction energy:                 -30.95 kcal/mol
   
3. Exchange-correlation interaction energy:           27.55 kcal/mol
   
4. Polarization interaction energy:                   -5.44 kcal/mol
   
5. Charge transfer interaction energy:               -15.25 kcal/mol
   
6. Grimme's dispersion interaction:                   -2.22 kcal/mol
   
7. ----------------------------------------------------------------
   
8. Total interaction energy:                         -26.31 kcal/mol
   
9. 
   
10. WITH BSSE correction:
    
11. Electrostatic interaction energy:                 -30.95 kcal/mol
    
12. Exchange-correlation interaction energy:           27.55 kcal/mol
    
13. Polarization interaction energy:                   -5.44 kcal/mol
    
14. Charge transfer interaction energy:               -10.44 kcal/mol
    
15. Grimme's dispersion interaction:                   -2.22 kcal/mol
    
16. ----------------------------------------------------------------
    
17. Total interaction energy:                         -21.50 kcal/mol

复制代码

可以找到：

• 静电能（Electrostatic interaction energy）：纯静电相互作用；
• 交换能，或者叫排斥能（Exchange-correlation interaction energy）：由Pauli不相容原理导致的排斥能（比如空间位阻）；
• 极化能（Polarization interaction energy）：一个分子的电子云受到另一个分子的影响产生的变化所消耗的能量；
• 电荷转移能（Charge transfer interaction energy）：一个分子的电子云转移到另一个分子的变化所消耗的能量。这个能量，可以包含或者不包含BSSE，两个结果都已给出。
• 色散能（Grimme's dispersion interaction）：顾名思义。
• 总相互作用能：A+B->AB的能量变化。这里可以看到，Qbics输出了不包含BSSE的结果（-26.31 kcal/mol）和包含BSSE的结果（-21.50 kcal/mol）。
  

由此可见，Qbics可以非常简单明了定义分子片，然后计算相互作用能量及其分量，非常方便。

3 多体能量分解分析（MB-EDA）

多体能量分解是Qbics中一个独有的功能。Qbics除了可以对A+B->AB进行能量分解外，还可以对任意


A+B+C+ ... + X -> ABC...X

进行能量分解。该方法由唐振和coolrainbow开发并实现。详细原理请参见和引用：

• J. Chem. Theory Comput. 2023, 19, 1777

• Phys. Chem. Chem. Phys. 2024, 26, 17549

  
  

这里以(Ca2+)(SO42-)(H2O)5为例：

1. basis
   
2.     def2-svp
   
3. end
   
4. 
   
5. scf
   
6.     charge     0
   
7.     spin2p1    1
   
8.     type       U # This is needed for EDA.
   
9. end
   
10. 
    
11. grimmedisp
    
12.     type bj
    
13. end
    
14. 
    
15. eda
    
16.     type     mb_tso
    
17.     mb_level 4  # The many-body truncation level.
    
18.     frag  +2 1 1
    
19.     frag  -2 1 2-6
    
20.     frag   0 1 7-9
    
21.     frag   0 1 10-12
    
22.     frag   0 1 13-15
    
23.     frag   0 1 16-18
    
24.     frag   0 1 19-21
    
25. end
    
26. 
    
27. mol
    
28.    Ca      1.98604937      1.24582501      1.89905383
    
29.     S      0.69103489      3.45226578      0.98237460
    
30.     O      2.08708072      3.44779165      1.50308184
    
31.     O      0.24826049      2.05102931      0.73622996
    
32.     O     -0.21032250      4.08940001      1.98323165
    
33.     O      0.63912085      4.22084216     -0.29304506
    
34.     O      3.02175028      0.83377391     -0.05175855
    
35.     H      3.13194669      0.09601115     -0.65158745
    
36.     H      2.90415260      1.59241194     -0.62348488
    
37.     O      2.51307423      2.92493798     -1.72944079
    
38.     H      1.80535522      3.44397361     -1.34738364
    
39.     H      2.91500306      3.50585612     -2.37536585
    
40.     O      0.60914213     -0.51387571      1.58896733
    
41.     H      0.07645358     -1.26917200      1.83796991
    
42.     H      0.01571971      0.04244367      1.08439778
    
43.     O      0.73702072      1.93958625      3.70546127
    
44.     H      0.73393089      2.03509947      4.65787899
    
45.     H      0.20818330      2.67203308      3.38909995
    
46.     O      3.87639195      1.93942131      2.91131853
    
47.     H      4.78820918      1.95286631      3.20224881
    
48.     H      3.69756418      2.83995160      2.64058309
    
49. end
    
50. 
    
51. task
    
52.     eda b3lyp
    
53. end

复制代码

这里的输入和之前的几乎一样，唯一的区别就在eda：

• type mb_tso 这里用多体版的TSO-EDA（推荐）。也可以使用mb_gks。
• mb_level 4 多体截断的级别。这里表示最多考虑到4体相互作用。也可以只到3。注意，更高的数字计算量将会非常大，需要(C(n, 1)+C(n,2)+C(n,3)+...+C(n, mb_level))*6个计算。
• frag +2 1 1 将原子1设为一个碎片，电荷和自选多重度为+2和1。
• frag -2 1 2-6 将原子2-6设为一个碎片，电荷和自选多重度为-2和1。
• frag 0 1 7-9 将原子7-9设为一个碎片，电荷和自选多重度为0和1。
• frag 0 1 10-12 类似
• frag 0 1 13-15 类似
• frag 0 1 16-18 类似
• frag 0 1 19-21 类似
  

碎片的定义见图片：




将该文件存为eda.inp，用Qbics执行(-n 24 表示用24核并行)：

1. qbics eda.inp -n 24 > eda.out

复制代码

结果可在eda.out找到：

1. Table 5. Summary (kcal/mol).
   
2. ---------------------------------------------------------------------------------------------------------------------------------
   
3.     Interactions         delE_el         delE_xc         delE_pl         delE_ct       delE_bsse       delE_disp        delE_tot
   
4. ---------------------------------------------------------------------------------------------------------------------------------
   
5.    SUM of 2-body -6.73148230E+02  1.47696370E+02 -1.11395126E+02 -2.04547116E+02  1.03684957E+02 -2.06082308E+01 -7.58317375E+02
   
6.    SUM of 3-body  1.42501290E-08 -2.67700179E+00  6.09855572E+01  1.30716087E+02 -4.86421459E+01  1.07484123E+01  1.51130908E+02
   
7.    SUM of 4-body -2.82504918E-08  5.22296932E-01 -4.66678754E+00 -4.78862109E+01  1.47634652E+01 -1.52639538E+01 -5.25311901E+01
   
8.           Remain  1.63009297E-08 -9.27025377E-02  2.21610464E-01  1.02616142E+01 -2.89029982E+00  7.47122033E+00  1.49714426E+01
   
9. ---------------------------------------------------------------------------------------------------------------------------------
   
10.              SUM -6.73148230E+02  1.45448963E+02 -5.48547462E+01 -1.11455626E+02  6.69159769E+01 -1.76525519E+01 -6.44746214E+02

复制代码

我们可以看到，总相互作用能为-644.75 kcal/mol，分解为二体、三体、四体和余项。二体项是最重要的项 （-758.32 kcal/mol），而三体项也很重要，但anti-cooperative（破坏团簇稳定）（+151.13 kcal/mol）。四体项小（-52.53 kcal/mol，略微cooperative）。与低阶项相比，剩余项（5、6 和 7 体之和）非常小（+14.97 kcal/mol）。

我们还可以看到，静电能和交换能是高度可加（additive）的（超过3体几乎为零），而极化和电荷转移能是非可加（non-additive）的：



对于不同类型的团簇，三体效应（多体相互作用）可能大不相同。



简单说一下原理。对于多体相互作用，


最重要的就是三体效应，其效果可以分解为三种(Phys. Chem. Chem. Phys. 2024, 26, 17549):

• delta E(3) < 0：表示团簇中单体的协同cooperative效应。多体相互作用正在稳定团簇。这经常出现在氢键团簇中，例如水簇。
• delta E(3) > 0：表示团簇中单体的反协同anti-cooperative作用。多体相互作用正在破坏团簇的稳定性。这通常出现在带电物质团簇中，例如离子液体团簇。
• delta E(3) = 0：表示团簇中单体的非协同non-cooperative效应。集群中几乎没有多体相互作用。这通常出现在没有电荷或氢键的分子簇中。
  

由此可见，在Qbics中，MB-EDA方法可以几乎一键实现，没有任何难度。该方法已在J. Comput. Chem. 2025, 46, e70114 使用，用于分析团簇的生长规律。