CASEA的奇思妙想

Finding transition states with NEB-TS

发表于 更新于 分类于
本文字数： 5.8k 阅读时长 ≈ 5 分钟

使用NEB-TS方法寻找过渡态

Finding transition states with NEB-TS

使用NEB-TS方法寻找过渡态

假设你想预测乙酸甲酯水解成乙酸和甲醇 的过渡态(transition state TS)
也许你想寻找水解机理，甚至计算 ΔG 并尝试预测反应速率。
在ORCA，有一个黑箱方法 NEB-TS（from Nudged Elastic Band with TS optimization），能够从反应物和产物的几何结构中寻找到过渡态结构

定义反应物和产物

第一步是优化反应物和产物，请注意是反应式上的所有分子集合。
首先，对于反应物，在绘制出初始结构后，进行结构优化

1
!B3LYP DEF2-SVP OPT D4

亲核水与羰基形成氢键，重新排列原子，猜测产物结构，用同样的方法优化后
可以开始寻找过渡态

为了使NEB-TS工作，反应物和产物几何结构中的原子编号必须相同。反应物中的碳原子顺序和产物xyz坐标中的碳原子顺序相同。
最简单的方法是，从一个几何结构开始，复制一份，通过移动原子来画另一个结构

如果使用Avogadro，你需要禁用调整氢功能，否则氢原子在键断裂时自动添加，并且不一定是你想要的位置。

The idea is not to optimize all molecules separately and just join then in a single file, but rather already build your “reactant” and “product” adducts in a certain orientation that makes sense with respect to the expected reaction, like done in the example.

制作NEB-TS输入文件

只需在主要输入里添加NEB-TS关键词，并在%NEB选项中添加产物的坐标。

1
2
3
!B3LYP DEF2-SVP D4 NEB-TS FREQ
%NEB NEB_END_XYZFILE "product.xyz" END
* XYZfile 0 1 reactant.xyz

如果你有一个猜测的过渡态结构，可以使用如下命令

1
2
3
4
5
!B3LYP DEF2-SVP D4 NEB-TS FREQ
%NEB NEB_END_XYZFILE "product.xyz"
     NEB_TS_XYZFILE "guessTS.xyz"
END
* XYZfile 0 1 reactant.xyz

反应物和产物在默认运行情况下不会自动优化，所以确保使用的优化方法相同，如果你想重新优化他们，在 %NEB 设置PREOPT_ENDS TRUE

运行NEB

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# reactant.xyz
14

C         -3.84500       -0.17571       -0.03921
C         -2.44492        0.32082        0.15466
H         -4.02399       -1.04086        0.60409
H         -4.00401       -0.43615       -1.08861
H         -4.55059        0.61392        0.23466
O         -1.55509       -0.65188       -0.18018
C         -0.18433       -0.27497       -0.02916
H          0.05373        0.57033       -0.68208
H          0.43501       -1.12619       -0.32499
H          0.03294       -0.03494        1.01615
O         -2.15827        1.44183        0.55653
O         -0.06767        3.07881        1.16647
H         -0.50923        3.90074        1.42789
H         -0.83167        2.50085        0.95538
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# product.xyz
14 

C         -3.44596       -0.05388        0.06874
C         -2.45275        1.00891        0.40403
H         -3.34585       -0.89157        0.76374
H         -3.29683       -0.38937       -0.95986
H         -4.45598        0.35325        0.17275
O          0.13193       -1.79894       -0.30248
C          1.49893       -1.43301       -0.26662
H          1.66462       -0.57388       -0.92141
H          2.10288       -2.27456       -0.61417
H          1.78363       -1.18238        0.75829
O         -2.67152        2.20314        0.50109
O         -1.20490        0.52554        0.57750
H         -0.69584        1.33651        0.80013
H         -0.37607       -1.03530        0.03854

计算开始时，首先输出以下信息，包含一些方法的内容

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
----------------------
NEB settings
----------------------
Method type                             ....  climbing image
Threshold for climbing image            ....  2.00e-02 Eh/Bohr
Free endpoints                          ....  off
Tangent type                            ....  improved
Number of intermediate images           ....  8
Number of images free to move           ....  8
Spring type for image distribution      ....  distance between adjacent images
Spring constant                         ....  energy weighted (0.0100 -to- 0.1000) Eh/Bohr^2
Spring force perp. to the path          ....  none
Generation of initial path              ....  image dependent pair potential
Initial path via TS guess               ....  off
[...]

然后构建初始路径

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Generation of  the initial path:
Minimize RMSD between reactant and product configurations     .... done
       RMSD before minimization         ....  2.1153 Bohr
       RMSD after minimization          ....  1.6611 Bohr
Performing linear interpolation         ....  done
Interpolation using                     ....  IDPP (Image Dependent Pair Potential)
IDPP-Settings:
       Remove global transl. and rot.   ...   true
       Convergence tolerance            ....  0.0100 1/Ang.^3
       Max. numer of iterations         ....    7000
       Spring constant                  ....  1.0000 1/Ang.^4
       Time step                        ....  0.0100 fs
       Max. movement per iteration      ....  0.0500 Ang.
       Full print                       ....  false
idpp initial path generation successfully converged in 117 iterations
Displacement along initial path:             11.6511 Bohr
Writing initial trajectory to file      ....  NEBTS_initial_path_trj.xyz

使用IDPP方法创建一系列构象，从反应物到产物。这些构象会使用NEB形式一起优化。
最初的路径会被保存为basename_initial_path_trj.xyz，建议检查初始路径的合理性。
此示例将得到8个构象片段。质子在离开时从水中转移到甲醇中。
轨迹文件是包括多个结构的文件，chemcraft，molden能直接打开。
注：Avogadro需要Extensions——>Animation，文件可以被打开甚至制作成动画。

迭代矩阵开始：

1
2
3
4
5
6
7
8
Starting iterations:

Optim.  Iteration  HEI  E(HEI)-E(0)  max(|Fp|)   RMS(Fp)    dS
Switch-on CI threshold               0.020000
   LBFGS     0      5    0.297110    0.116137   0.026058  11.6511
   LBFGS     1      5    0.277441    0.092448   0.021661  11.5993
   LBFGS     2      6    0.262453    0.111904   0.022265  11.5925
   LBFGS     3      6    0.252519    0.082144   0.017474  11.6060

在收敛之后，会输出climbing image" (CI)

TS过渡态优化

在迭代之后，NEB-CI 收敛

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
                                .--------------------.
          ----------------------| CI-NEB convergence |-------------------------
          Item                value                   Tolerance       Converged
          ---------------------------------------------------------------------
          RMS(Fp)             0.0006879134            0.0100000000      YES
          MAX(|Fp|)           0.0032907070            0.0200000000      YES
          RMS(FCI)            0.0006591030            0.0010000000      YES
          MAX(|FCI|)          0.0016597640            0.0020000000      YES
          ---------------------------------------------------------------------

The elastic band and climbing image have converged successfully to a MEP in 128 iterations!

                    *********************H U R R A Y*********************
                    ***        THE NEB OPTIMIZATION HAS CONVERGED     ***
                    *****************************************************

开始过渡态优化，首先，基于先前NEB信息猜测建立Hessian矩阵并初始化鞍点搜索
如果TS全部优化，将会输出NEB-TS的完整总结

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
            ***********************HURRAY************************
            ***        THE TS OPTIMIZATION HAS CONVERGED      ***
            *****************************************************

---------------------------------------------------------------
                      PATH SUMMARY FOR NEB-TS
---------------------------------------------------------------
All forces in Eh/Bohr. Global forces for TS.

Image     E(Eh)   dE(kcal/mol)  max(|Fp|)  RMS(Fp)
  0     0.000    -344.38971      0.00       0.00014   0.00003
  1     2.024    -344.38507      2.91       0.00144   0.00053
  2     3.395    -344.38062      5.70       0.00150   0.00055
  3     4.863    -344.38237      4.61       0.00167   0.00064
  4     5.663    -344.36742     13.98       0.00329   0.00102
  5     6.043    -344.33717     32.97       0.00274   0.00097
  6     6.336    -344.32049     43.44       0.00178   0.00058 <= CI
  7     6.754    -344.34916     25.45       0.00260   0.00098
  8     7.426    -344.38096      5.49       0.00143   0.00058
  9     9.661    -344.38666      1.91       0.00023   0.00007

TS的频率会被计算来验证是否只含有一个负频

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

Scaling factor for frequencies =  1.000000000  (already applied!)

   0:         0.00 cm**-1
   1:         0.00 cm**-1
   2:         0.00 cm**-1
   3:         0.00 cm**-1
   4:         0.00 cm**-1
   5:         0.00 cm**-1
   6:     -1149.57 cm**-1 ***imaginary mode***
   7:       106.27 cm**-1
   8:       135.34 cm**-1
   9:       188.41 cm**-1

输出证明优化后的TS结构确实是一个鞍点
最终从反应物到产物的轨迹被写入basename_MEP_trj.xyz
最终的过渡态几何结构被写入basename_NEB-TS_converged.xyz
从TS结构可以看出TS具有四面体碳原子，质子协同转移到离开的甲氧基
该计算在气相中运算以简化解释，如果包含溶剂，机理会更加复杂，质子转移可能通过外部水分子发生。

Possible intermediates（可能的中间体）

在NEB-CI过程中，会输出以下信息

1
2
3
Possible intermediate minimum found at image(s): 7
  => have a look at the .interp file.
  => it is often a good idea to divide the path calculation into smaller fragments.

在沿反应路径有局部最小值，实际上可能不止一个
沿着反应最小路径的能量输出在interp文件中，可以将能量-反应路径绘制成图
此例可以清楚地看到有两个中间体，也应该优化这个结构以正确预测反应能量，并且可能用这些新的反应物和产物猜测结构重新运行 NEB-TS