Ctrl+K

简单理解最初步的 ONIOM 能量计算

4. 简单理解最初步的 ONIOM 能量计算#

创建时间：2021-06-10

在这篇文档中，我们会使用 PySCF 实现 ONIOM 计算的最初步的结果，即单点能计算。ONIOM 方法的原始文献之一是 Dapprich, Frisch et al. [1]。

由于单点能本身是没有实际物理价值的，因此相对能量、或者光谱性质、几何构型才是真正有价值的输出量。但得到这些输出量会比较复杂，特别是大多数 (除了 Gaussian 外) 以代价较大的第一性为卖点的量化程序，对半经验方法的光谱信息的支持并不好。我们并不对这些问题作更细致的讨论。正因为 Gaussian 对 ONIOM 的各种光谱数据、结构性质和溶剂化的支持非常充分，而其它软件在没有特化程序时似乎只能计算单点能；因此这些特色使得 Gaussian 确实地成为适合药物设计的软件之一。

尽管 PySCF 确实可以给出 ONIOM 的计算结果，但为了方便地得到 ONIOM 的分块信息，我们仍然还需要 Gaussian 的支持。

from pyscf import gto, scf, mp, semiempirical
import numpy as np

HARTREE2KCAL = semiempirical.mopac_param.HARTREE2KCAL

4.1. 简单问题：三氟代乙醛的 2-layer 计算#

该问题的 Gaussian 输入卡置于 2-layer.gjf 中。

! cat 2-layer.gjf

#p ONIOM(MP2(Full)/6-31G:HF/6-31G)=InputFiles NoSymm
 
2-layer ONIOM, modified from Gaussian keyword list
 
0 1 0 1 0 1
  F     -1.041506214819     0.000000000000    -2.126109488809 L
  F     -2.033681935634    -1.142892069126    -0.412218766901 L
  F     -2.033681935634     1.142892069126    -0.412218766901 L
  C     -1.299038105677     0.000000000000    -0.750000000000 L H 5 0.7 0.8
  C      0.000000000000     0.000000000000     0.000000000000 H
  H      0.000000000000     0.000000000000     1.100000000000 H
  O      1.125833024920     0.000000000000    -0.650000000000 H

计算完毕后，最关键的 Gaussian 的输出是

ONIOM: gridpoint  1 method:  low   system:  model energy:  -113.796286020376
ONIOM: gridpoint  2 method:  high  system:  model energy:  -114.023466233349
ONIOM: gridpoint  3 method:  low   system:  real  energy:  -449.289045248409
ONIOM: extrapolated energy =    -449.516225461381

我们可以知道，最终能量的计算是通过下式完成的：

\[ E_\mathrm{ONIOM2} = E_\mathrm{low} (\mathrm{real}) + E_\mathrm{high} (\mathrm{model}) - E_\mathrm{low} (\mathrm{model}) \]

对于 2-layer 的 ONIOM 计算，它要细分成三个计算任务：

任务序号

计算级别 (method)

计算体系 (system)

能量 / Hartree

1

HF/6-31G (low)

(H)CHO (model)

-113.796286020376

2

MP2(Full)/6-31G (high)

(H)CHO (model)

-114.023466233349

3

HF/6-31G (low)

CF3-CHO (real)

-449.289045248409

注意到我们指定的模型层 (model) 是最后三个原子；但实际上，由于模型层 (model) 与全局层 (real) 之间有键相互作用 (第 4 个 real 层甲基碳原子与第 5 个 model 层醛碳原子)，因此不能简单粗暴地在此处断键。

上面的 Gaussian 输入卡会在计算模型层时引入氢原子；其引入并非直接将第 4 个碳原子替换为氢，而同时还要缩放键长度。在上述表格中，人为补上去的氢原子也由括号作标识。H 5 0.7 0.8 即指在低计算级别 (low) 缩放到 0.7 倍，而在高计算级别 (high) 缩放到 0.8 倍。

4.1.1. 任务 1：\(E_\mathrm{low} (\mathrm{model})\) 计算#

低计算级别 (low) 与高计算级别 (high) 都需要计算模型层 (model) 的分子片段。但由于氢原子添补方式的不同，使得我们需要分别定义这两个分子，及其后续计算。我们首先给出 high 级别的计算结果。

Low 级别下，引入的氢原子键长是对应 4 号碳与 5 号碳的 0.7 倍。

mol_model_low = gto.Mole()
mol_model_low.atom = """
  H     -0.90932667         0.                -0.525     
  C      0.000000000000     0.000000000000     0.000000000000 
  H      0.000000000000     0.000000000000     1.100000000000 
  O      1.125833024920     0.000000000000    -0.650000000000 
"""
mol_model_low.basis = "6-31G"
mol_model_low.verbose = 0
mol_model_low.build()

<pyscf.gto.mole.Mole at 0x7f277b9a56a0>

mf_model_low = scf.RHF(mol_model_low).run()
eng_model_low = mf_model_low.e_tot
eng_model_low

-113.79628602351522

4.1.2. 任务 2：\(E_\mathrm{high} (\mathrm{model})\) 计算#

High 级别下，引入的氢原子键长是对应 4 号碳与 5 号碳的 0.8 倍。

mol_model_high = gto.Mole()
mol_model_high.atom = """
  H     -1.03923048         0.                -0.6
  C      0.000000000000     0.000000000000     0.000000000000 
  H      0.000000000000     0.000000000000     1.100000000000 
  O      1.125833024920     0.000000000000    -0.650000000000 
"""
mol_model_high.basis = "6-31G"
mol_model_high.verbose = 0
mol_model_high.build()

<pyscf.gto.mole.Mole at 0x7f276ae590a0>

mf_model_high = mp.MP2(mol_model_high).run()
eng_model_high = mf_model_high.e_tot
eng_model_high

-114.02346625303319

4.1.3. 任务 3：\(E_\mathrm{low} (\mathrm{real})\) 计算#

mol_real = gto.Mole()
mol_real.atom = """
  F     -1.041506214819     0.000000000000    -2.126109488809 
  F     -2.033681935634    -1.142892069126    -0.412218766901 
  F     -2.033681935634     1.142892069126    -0.412218766901 
  C     -1.299038105677     0.000000000000    -0.750000000000 
  C      0.000000000000     0.000000000000     0.000000000000 
  H      0.000000000000     0.000000000000     1.100000000000 
  O      1.125833024920     0.000000000000    -0.650000000000 
"""
mol_real.basis = "6-31G"
mol_real.verbose = 0
mol_real.build()

<pyscf.gto.mole.Mole at 0x7f276ae597c0>

mf_real_low = scf.RHF(mol_real).run()
eng_real_low = mf_real_low.e_tot
eng_real_low

-449.28904521819294

4.1.4. 能量的统合#

回顾 2-layer ONIOM 的能量统合方式：

\[ E_\mathrm{ONIOM2} = E_\mathrm{low} (\mathrm{real}) + E_\mathrm{high} (\mathrm{model}) - E_\mathrm{low} (\mathrm{model}) \]
eng_real_low + eng_model_high - eng_model_low

-449.5162254477109

Gaussian 的结果是 -449.516225461381。

4.2. 较复杂问题：丙醛的 3-layer 计算#

该问题的 Gaussian 输入卡置于 3-layer.gjf 中。

! cat 3-layer.gjf

#p ONIOM(MP2(Full)/6-311g*:HF/6-31g:MINDO3)=InputFiles NoSymm

3-layer ONIOM, modified from Gaussian test job 0679

   0   1   0   1   0   1   0   1   0   1   0   1   0   1
 C  -0.006049274275    0.000000000000    0.066754956170 H
 O   0.011403425950    0.000000000000    1.308239478983 H
 H   0.944762558657    0.000000000000   -0.507359536461 H
 C  -1.307562483867    0.000000000000   -0.766510748030 M H 1 0.723886 0.723886 0.723886
 C  -1.047480751885    0.000000000000   -2.301387120377 L H 4 0.723886 0.723886 0.723886
 H  -1.903669606697   -0.885256630266   -0.468844831106 M
 H  -1.903669606697    0.885256630266   -0.468844831106 M
 H  -1.988817319373    0.000000000000   -2.842389774687 L
 H  -0.482972255230    0.881286097766   -2.591806824941 L
 H  -0.482972255230   -0.881286097766   -2.591806824941 L

计算完毕后，最关键的 Gaussian 的输出是

ONIOM: gridpoint  1 method:  low   system:  model energy:    -0.033309689782
ONIOM: gridpoint  2 method:  med   system:  model energy:  -113.805007046900
ONIOM: gridpoint  3 method:  low   system:  mid   energy:    -0.059535553265
ONIOM: gridpoint  4 method:  high  system:  model energy:  -114.255323473041
ONIOM: gridpoint  5 method:  med   system:  mid   energy:  -152.836036428501
ONIOM: gridpoint  6 method:  low   system:  real  energy:    -0.068015765875
ONIOM: extrapolated energy =    -153.294833067253

由于外推表达式比较复杂，原始文献也使用了简化的表达式作能量结果的表示 (需要注意，\(E_1 = E_\mathrm{low} (\mathrm{model})\) 没有在能量统合公式中)：

\[ E_\mathrm{ONIOM3} = E_6 - E_3 + E_5 - E_2 + E_4 \]

对于 3-layer 的 ONIOM 计算，它要细分成 5 个计算任务：

任务序号

计算级别 (method)

计算体系 (system)

能量 / Hartree

2

HF/6-31g (med)

(H)CHO (model)

-113.805007046900

3

MINDO/3 (low)

(H)CH2-CHO (mid)

-0.059535553265

4

MP2(Full)/6-311g* (high)

(H)CHO (model)

-114.255323473041

5

HF/6-31g (med)

(H)CH2-CHO (mid)

-152.836036428501

6

MINDO/3 (low)

CH3-CH2-CHO (real)

-0.068015765875

4.2.1. 使用 Gaussian 查看每个计算任务的分子构型#

一旦体系增大，模型 (model)、中间 (mid/intermediate) 与全局 (real) 层之间的断键数目增多，这时人为地确定具体的分子构型就会变得困难了。

一种比较简单粗暴的方式是使用 Gaussian ONIOM 关键词的 OnlyInputFiles 选项；这样就可以在不执行 ONIOM 具体计算的情况下，把每个需要计算的分子片段信息打印出来。如果加入关键词 InputFiles，那么可以同时打印分子片段以及计算 ONIOM。

我们以第 5 个格点为例 (中层 (mid/intermediate)，中等级算量 (med))。其在 Gaussian 中的输出是 (略去 Gaussian 推断的成键关系)

ONIOM: generating point  5 -- med level on mid system.

--------------------------------------------------------------------------------
#P Test IOp(2/15=1,5/32=2,5/38=1) HF/6-31G

3-layer ONIOM, modified from Gaussian test job 0679
Point  5 -- med level on mid system.

0,1
 C                                                -0.006049274275      0.000000000000      0.066754956170
 O                                                 0.011403425950      0.000000000000      1.308239478983
 H                                                 0.944762558657      0.000000000000     -0.507359536461
 C                                                -1.307562483867      0.000000000000     -0.766510748030
 H(Iso=12)                                        -1.119292959229      0.000000000000     -1.877586265703
 H                                                -1.903669606697     -0.885256630266     -0.468844831106
 H                                                -1.903669606697      0.885256630266     -0.468844831106
 Bq-#1(Iso=1.00782504,Spin=1,ZNuc=1.,GFac=2.792846)         -1.988817319373      0.000000000000     -2.842389774687
 Bq-#1(Iso=1.00782504,Spin=1,ZNuc=1.,GFac=2.792846)         -0.482972255230      0.881286097766     -2.591806824941
 Bq-#1(Iso=1.00782504,Spin=1,ZNuc=1.,GFac=2.792846)         -0.482972255230     -0.881286097766     -2.591806824941

我们可以看出第 5 号位的氢原子实际上是替代了低层 (low) 的碳原子。之所以要设置为氢原子 12 质量的同位素，是为了跟进一步的分子力与频率分析作准备；在单纯讨论能量时不需要考虑原子质量的问题。同时，最后的三个低层 (low) 原子都被设置为 Bq 原子，即虚原子。我们实际需要带入能量计算的原子就是前 7 个原子了。

4.2.2. 任务 2：\(E_2 = E_\mathrm{med} (\mathrm{model})\)#

mol_2 = gto.Mole()
mol_2.atom = """
C    -0.006049274275      0.000000000000      0.066754956170
O     0.011403425950      0.000000000000      1.308239478983
H     0.944762558657      0.000000000000     -0.507359536461
H    -0.948196465514      0.000000000000     -0.536434421381
"""
mol_2.basis = "6-31G"
mol_2.verbose = 0
mol_2.build()

<pyscf.gto.mole.Mole at 0x7f276ae59a30>

mf_2 = scf.RHF(mol_2).run()
eng_2 = mf_2.e_tot
eng_2

-113.80500704910236

4.2.3. 任务 3：\(E_3 = E_\mathrm{low} (\mathrm{mid})\)#

需要注意，Gaussian 在半经验方法中输出的能量并非接近于单点能，而是接近于原子化能。因此，在使用 PySCF 时尽量不要直接用半经验的 e_tot 变量作结果输出。

mol_3 = gto.Mole()
mol_3.atom = """
C     -0.006049274275      0.000000000000      0.066754956170
O      0.011403425950      0.000000000000      1.308239478983
H      0.944762558657      0.000000000000     -0.507359536461
C     -1.307562483867      0.000000000000     -0.766510748030
H     -1.119292959229      0.000000000000     -1.877586265703
H     -1.903669606697     -0.885256630266     -0.468844831106
H     -1.903669606697      0.885256630266     -0.468844831106
"""
mol_3.verbose = 0
mol_3.build()

<pyscf.gto.mole.Mole at 0x7f276ae59a60>

mf_3 = semiempirical.MINDO3(mol_3).run()
eng_3 = mf_3.e_heat_formation / HARTREE2KCAL
eng_3

-0.059543662948324035

这与 Gaussian 的结果有略微差别，但差别在 5e-3 kcal/mol 上，我们可以认为这时可以忽略的差距了。

4.2.4. 任务 4：\(E_4 = E_\mathrm{high} (\mathrm{model})\)#

需要注意，即使这个分子与 \(E_2 = E_\mathrm{med} (\mathrm{model})\) 所使用的分子相同 (因为使用了相同的断键补氢系数 \(g = 0.723886\))；但 \(E_2\) 的基组是中等级 (med) 的 6-31G，而 \(E_4\) 则是高等级 (high) 的 6-311G*。

mol_4 = gto.Mole()
mol_4.atom = """
C    -0.006049274275      0.000000000000      0.066754956170
O     0.011403425950      0.000000000000      1.308239478983
H     0.944762558657      0.000000000000     -0.507359536461
H    -0.948196465514      0.000000000000     -0.536434421381
"""
mol_4.basis = "6-311G*"
mol_4.verbose = 0
mol_4.build()

<pyscf.gto.mole.Mole at 0x7f276ae598e0>

mf_4 = mp.MP2(mol_4).run()
eng_4 = mf_4.e_tot
eng_4

-114.25532346425724

4.2.5. 任务 5：\(E_5 = E_\mathrm{med} (\mathrm{mid})\)#

mol_5 = gto.Mole()
mol_5.atom = """
C     -0.006049274275      0.000000000000      0.066754956170
O      0.011403425950      0.000000000000      1.308239478983
H      0.944762558657      0.000000000000     -0.507359536461
C     -1.307562483867      0.000000000000     -0.766510748030
H     -1.119292959229      0.000000000000     -1.877586265703
H     -1.903669606697     -0.885256630266     -0.468844831106
H     -1.903669606697      0.885256630266     -0.468844831106
"""
mol_5.basis = "6-31G"
mol_5.verbose = 0
mol_5.build()

<pyscf.gto.mole.Mole at 0x7f276ae77670>

mf_5 = scf.RHF(mol_5).run()
eng_5 = mf_5.e_tot
eng_5

-152.83603642639747

4.2.6. 任务 6：\(E_6 = E_\mathrm{low} (\mathrm{real})\)#

mol_6 = gto.Mole()
mol_6.atom = """
C     -0.006049274275      0.000000000000      0.066754956170
O      0.011403425950      0.000000000000      1.308239478983
H      0.944762558657      0.000000000000     -0.507359536461
C     -1.307562483867      0.000000000000     -0.766510748030
C     -1.047480751885      0.000000000000     -2.301387120377
H     -1.903669606697     -0.885256630266     -0.468844831106
H     -1.903669606697      0.885256630266     -0.468844831106
H     -1.988817319373      0.000000000000     -2.842389774687
H     -0.482972255230      0.881286097766     -2.591806824941
H     -0.482972255230     -0.881286097766     -2.591806824941
"""
mol_6.verbose = 0
mol_6.build()

<pyscf.gto.mole.Mole at 0x7f276ae77d90>

mf_6 = semiempirical.MINDO3(mol_6).run()
eng_6 = mf_6.e_heat_formation / HARTREE2KCAL
eng_6

-0.06801695739152386

4.2.7. 能量的统合#

回顾 3-layer ONIOM 的能量统合方式：

\[ E_\mathrm{ONIOM3} = E_6 - E_3 + E_5 - E_2 + E_4 \]
eng_6 - eng_3 + eng_5 - eng_2 + eng_4

-153.29482613599555

Gaussian 的结果是 -153.294833067253 Hartree。我们使用 PySCF 给出的计算结果与 Gaussian 相差 4e-3 kcal/mol。

文档目录