Gaussian + xtb 结合做过渡态搜索

前言

过渡态理论无疑是当前研究化学动力学最成功的理论之一。但是并不是每一个化学过程都存在过渡态，例如离子键断裂成为阴阳离子或阴阳离子的结合。过渡态是势能面的特征，因此只有涉及核坐标变化的化学过程才能有过渡态。单纯的得失电子、电子转移等过程并没有过渡态。

在过去，想使用量子化学的手段搜索过渡态，需要使用 Gaussian 程序进行过渡态搜索，尽管 Gaussian 在搜索过渡态时支持 qst2、qst3 算法，但是并不好用。同时如果涉及到过渡金属反应的过渡态时，作为预搜索的 PM6、PM7 等半经验方法往往都是没办法用的。如果直接上 DFT 方法，遇到大体系的话，计算时间长不说，还不一定收敛到想要的过渡态，到时候又得重算，十分浪费时间。

不过这种情况在 Grimme 推出 xtb 程序之后明显的改善了很多。xtb 程序内置了大量 Grimme 提出的半经验 DFT 方法，可以很好的替代 PM6、PM7 这种半经验方法，并且半经验 DFT 方法在处理过渡金属时的效果要比 PM6、PM7 好得多。本文主要介绍如何使用 Gaussian + xtb 结合做过渡态搜索。

本文主要参考了 Sob 老师的文章，本文所有的操作都是在 Rocky Linux 下完成的。

• 《简谈 Gaussian 里找过渡态的关键词 opt=TS 和 QST2、3》(http://sobereva.com/460)
• 《将 Gaussian 与 Grimme 的 xtb 程序联用搜索过渡态、产生 IRC、做振动分析》(http://sobereva.com/421)

安装 xtb

首先在 xtb release 页面下载已经编译好的 xtb，也就是后缀为 .tar.xz 的版本，注意别下成 source 版本了。接着将 .tar.xz 文件解压至 \home\kimariyb\xtb 文件夹下。

用 nano 或者 vim 编辑 ~/.bashrc，加入以下语句后，保存即可。

export PATH=$PATH:/home/kimariyb/xtb/bin
export XTBPATH=/home/kimariyb/xtb/share/xtb
export OMP_NUM_THREADS=8
export MKL_NUM_THREADS=8
export OMP_STACKSIZE=1000m
ulimit -s unlimited

gau_xtb 接口

由于 Gaussian 程序存在独特的 external 功能，因此可以通过这个功能调用第三方的脚本。这里就不细说 external 功能了，有兴趣的可以自行浏览 Sob 老师的文章。同时，Sob 老师贴心的为广大计算化学玩家提供了其自己实现的接口 gau_xtb，下载地址为：http://sobereva.com/soft/gau_xtb/

如果在发表文章时，使用了 Sob 的 gau_xtb 接口，需要按照如下格式进行引用。

Tian Lu, gau_xtb: A Gaussian interface for xtb code, http://sobereva.com/soft/gau_xtb (accessed month day, year)

接着将 gau_xtb 解压在某一文件夹下，在 gau_xtb 中，genxyz、extderi 和 xtb.sh 文件用来调用 xtb 计算。在调用 xtb 之前，还需要把 xtb.sh 中的 OMP_NUM_THREADS 和 MKL_NUM_THREADS 变量设置为主机的核心数。最后，就可以在 Gaussian 输入文件中写入以下四条命令，用来调用 xtb 做优化、过渡态搜索、IRC 以及频率计算。

# external='./xtb.sh' force
# opt(ts,calcfc,noeigen,nomicro) external='./xtb.sh'
# IRC(maxpoints=20,calcfc) geom=check external='./xtb.sh'
# freq geom=allcheck external='./xtb.sh'

搜索过渡态

这里随便找一个例子，例如 J. Am. Chem. Soc. 2018, 140, 18124−18131 一文探究了如下反应的过渡态。

原文是使用 Gaussian 在 B3LYP/6-31+G** 级别下结合 SMD 溶剂模型优化的。笔者并不清楚作者在优化时加弥散函数的意义，并没有看见对当前体系加弥散的好处。因此笔者不会一模一样还原原作者的基组和方法。同时为了省事，这里就只考虑气象下的反应。反应物、生成物以及过渡态和中间产物的笛卡尔坐标附在文章的 SI 里了。

在搜索过渡态时，将平常搜索过渡态的关键词修改为如下关键词，这样就可以调用 xtb 了。先将反应物摆成大概猜测的过渡态的样子，接着调用 Gaussian 计算。由于 xtb 计算的速度非常快，小体系几秒钟就算完了。

%chk=mol.chk
%nproc=1
# opt(ts,calcfc,noeigen,nomicro) external='./xtb.sh'

接着可以在 example 文件夹中找到 freq.gjf 文件，将其拖到 gau_xtb 文件夹中，用 Gaussian 运行它，就可以得到 freq 任务的输出文件。可以根据其虚频的方向初步确定是不是需要找的过渡态。

确认过渡态

在大致确认过渡态后，就可以使用 Gaussian 的 DFT 方法计算过渡态、IRC 以及能量了。这里笔者使用 M06-2X/6-311G** 做过渡态搜索，同时在相同级别下做 IRC（IRC 曲线如下所示）。最后在 revDSD-PBEP86-D3(BJ)/def2-QZVPP 级别下计算单点能。

最后可以将反应物、过渡态、中间产物和生成物的能量绘制能量曲线图。这里为了省事就不自己绘制了，直接把原文的图片搬过来。

附：Python 绘制 IRC 曲线图的代码

import proplot as pplt
import pandas as pd

from proplot import rc

# 读取 txt 文件中的数据
data = pd.read_csv('irc1_tot_ener.txt', header=None, 
                   delim_whitespace=True, skiprows=4, names=['X', 'Y'])

# 设置绘图的默认参数，如字体、字号等
rc['font.name'] = 'Arial'
rc['title.size'] = 14
rc['label.size'] = 12
rc['font.size'] = 10.5
rc['tick.width'] = 1.3
rc['meta.width'] = 1.3
rc['label.weight'] = 'bold'
rc['tick.labelweight'] = 'bold'
rc['ytick.major.size'] = 4.6
rc['ytick.minor.size'] = 2.5
rc['xtick.major.size'] = 4.6
rc['xtick.minor.size'] = 2.5

xlabel = 'Intrinsic Reaction Coordinate'
ylabel = 'Free Energy (kcal/mol)'
xlim = (0, 8, 2)
ylim = (-50, 5, 10)

# 将 DataFrame 对象中 Y 减去最小值 yMin 同时乘以 627.51
data['Y'] = (data['Y'] - data['Y'].max()) * 627.51

# 创建实例
fig, ax = pplt.subplots(figsize=(5.4 * 0.9, 4 * 0.9))

# 绘制曲线图
ax.plot(data['X'], data['Y'], linewidth=1.3, color='#FFBE7A', zorder=1)
# 绘制实际的点
ax.scatter(data['X'], data['Y'], color='#FFBE7A', zorder=2, markersize=12)

# 格式化图像
fig.format(
    grid=False, ylabel=ylabel, xlabel=xlabel,
    xlim=(xlim[0], xlim[1]), xminorlocator=(xlim[2] / 2), xlocator=xlim[2], 
    ylim=(ylim[0], ylim[1]), yminorlocator=(ylim[2] / 2), ylocator=ylim[2]
)

# 显示图像
fig.show()

# 保存图像
fig.savefig("irc.png", dpi=300, bbox_inches="tight")