将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析

将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析 

文/Sobereva @北京科音

First release: 2018-May-29  Last update: 2018-Aug-1


0 前言

在《盘点Grimme迄今对理论化学的贡献》(http://sobereva.com/388)一文中笔者曾简单提到GFN-xTB方法,说白了就是类似于DFTB那种半经验意味的DFT,精度不错普适性也好,对与之耗时相仿佛的半经验方法和DFTB带来了极大的冲击。Grimme在其主页上也给出了名为xtb的专门做GFN-xTB计算的程序,目前可按照此网页上的说明申请免费下载:https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/xtb/xtb。xtb程序用过的人都说好,速度很快,普适性挺好,推出不久已经开始有很多人使用了,在《2018年度计算化学公社杯最常用的量子化学程序和DFT泛函投票结果统计》(http://sobereva.com/420)里已经有一些得票率了,而且一些第三方程序已经支持xtb了,比如Multiwfn (http://sobereva.com/multiwfn)可以读取xtb的输出来绘制红外光谱,molclus程序可以结合xtb来做构型和构象搜索(见http://www.keinsci.com/research/molclus.html),另外计算化学公社上fhh2626还写了NAMD与xtb结合做QM/MM的界面(http://bbs.keinsci.com/thread-7583-1-1.html)。

xtb程序目前可以做单点、优化、振动分析等任务,但是对于一般计算化学研究者来说,还希望能够找过渡态、产生IRC,并且希望振动分析的结果(特别是虚频模式)可以可视化,但这些xtb程序目前还做不到。好在Gaussian从09开始加入了external关键词,在进行极小点/过渡态优化、IRC、振动分析等任务时,可以从外部文件直接读入能量、受力、Hessian,而外部文件的这些信息可以用任意程序来产生,当然也包括xtb,不过需要自己写个接口才能实现。因此Gaussian可以被当做一个“optimizer”来使用,这种用法似于ASE(atomic simulation environment)程序。

本文的目的就是介绍如何将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC,并且使得xtb振动分析结果能够被gview来可视化,从而弥补xtb的不足、极大地扩展xtb的实用价值。本文Gaussian使用G09 E.01版,使用的xtb程序是2018-Jul-13发布的版本,这个版本实际上默认用的是在未来才要发表的GFN-xTB二代方法的尝鲜版(带了DFT-D4校正)。本文的做法对于以后版本也许会有不兼容的问题,请届时根据本文的说明随机应变恰当修改Fortran源码和脚本(本文的脚本至少对于更早期的xtb不适合,因为输出的文件格式不同)。

本文涉及的各种文件都可以从这里下载:gau_xtb.zip(xtb可执行程序不在里面,需自行申请)。

笔者还另写了一篇文章,《Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析》(http://sobereva.com/422),将Gaussian和ORCA联用也带来很多好处,有兴趣可以看看。


1 xtb使用简介

借本文的机会简单介绍一下xtb的使用。向Grimme发送邮件申请xtb后,会给你下载链接,只能下载到Linux下的预编译版。

下载xtb包后解压,然后把这个目录加入$PATH和$XTBHOME环境变量,然后就可以在任何位置直接输入比如xtb test.xyz -opt > test.out命令进行计算了。此目录中有很多.开头的文件,这在Linux中是隐藏文件。其中.xtbrc里含有$set ... $end字段,用于设定各种任务的详细执行参数,比如几何优化收敛限,一般不需要改。此目录中还有很多.param开头的,是参数文件。

如果需要并行计算,计算前应运行以下命令,N是机子的CPU物理核心数(建议把以下内容加入到用户目录下的.bashrc,使得每次进入终端后自动生效):
export OMP_NUM_THREADS=N
export MKL_NUM_THREADS=N
export OMP_STACKSIZE=1000m
ulimit -s unlimited

xtb的输入文件就是一个xyz文件,这是最常用的记录分子结构格式之一,很多程序都可以产生。用Multiwfn产生也可以,可以把fch、pdb、mol、wfn、wfx、molden等Multiwfn能认的格式载入Multiwfn,然后进入主功能100的子功能2,选择输出xyz文件。由于这个文件格式非常简单,比如自己把Gaussian的.gjf文件编辑一下产生也可以。

xtb的详细使用说明见自带的HOWTO文件,也可以直接输入xtb命令查看常见选项的使用提示。几个比较常用的选项如下
-chrg:设定体系净电荷
-uhf:设定alpha电子数减beta电子数(相当于自旋多重度减1)
-sp:计算单点
-grad:计算梯度
-opt:优化
-hess:计算Hessian并做振动分析
-ohess:优化后计算Hessian并做振动分析
-omd:优化并做分子动力学
-gbsa:使用隐式溶剂模型。目前支持的溶剂有toluene、thf、methanol、h2o、ether、chcl3、acetonitrile、acetone、cs2
例如:
对yoshiko.xyz做真空中的单点计算。电荷为1,自旋多重度为2:xtb yoshiko.xyz -chrg 1 -uhf 1 -sp
对yohane.xyz做甲苯溶剂下优化和振动分析。体系是默认的中性单重态:xtb yohane.xyz -ohess -gbsa toluene

xtb运行时一方面会在屏幕上输出信息,同时也会在当前目录下产生一大堆文件。这些文件的含义在HOWTO末尾有说明。

xtb目前有解析梯度,但只支持数值Hessian。-hess任务做完会输出g98.out和g98_canmode.out。前者是模仿高斯freq输出格式来输出频率、红外强度、正则坐标。后者没用。

-opt任务产生的xtbopt.coord是最后结构的xyz坐标文件,xtbopt.log是含有优化过程每一帧的多帧xyz文件,后缀改为.xyz后就可以拖入VMD查看优化轨迹。

.xtbrc里$set ... $end字段设的是一般用途的参数,如果对于计算的某些体系需要用特殊参数,可以自行编辑xyz文件,在后头写$set ... $end字段来设定。比如写
$set
chrg 1
uhf 1
$end
则此体系在计算的时候是当做带一个正电荷的二重态来计算,命令行里将就不用写-chrg 1 -uhf 1了。

优化时若要对比如1-30号原子冻结,则在xyz文件末尾写上
$set
fix 1-30
fixfc 0
$end
如果只是限制而不冻结的话,fixfc可以设具体的限制势力常数。

其它的xtb细节看自带的HOWTO,在北京科音(http://www.keinsci.com)的高级量子化学培训班里会对GFN-xTB的原理和xtb的用法做深入全面讲授。


2 Gaussian的external功能简介

关于Gaussian的external功能的使用,详见http://sobereva.com/g09/k_external.htm。简单来说,Gaussian的输入文件里写上比如external='./xtb.sh',在计算时候就会以这样的方式调用当前目录下的xtb.sh脚本:
xtb.sh layer InputFile OutputFile MsgFile FChkFile MatElFile
各个参数的含义可以看手册,后五个是文件名,这里我们主要关心的是其中第二个参数InputFile和第三个参数OutputFile。如果查看Gaussian输出文件,会发现这样的提示
 Running external command "./xtb.sh R"
         input file       "/sob/g09/scratch/Gau-28355.EIn"
         output file      "/sob/g09/scratch/Gau-28355.EOu"
         message file     "/sob/g09/scratch/Gau-28355.EMs"
         fchk file        "/sob/g09/scratch/Gau-28355.EFC"
         mat. el file     "/sob/g09/scratch/Gau-28355.EUF"
其中"/sob/g09/scratch/Gau-28355.EIn"和"/sob/g09/scratch/Gau-28355.EOu"正是分别传递给xtb.sh的InputFile和OutputFile的文件名。

InputFile文件是Gaussian产生的,记录了当前步的信息,格式为:
原子数  需要的导数  电荷  自旋多重度
原子1元素序号  X  Y  Z  MM电荷 MM原子类型
原子2元素序号  X  Y  Z  MM电荷 MM原子类型
...
原子N元素序号  X  Y  Z  MM电荷 MM原子类型
如果当前任务只需要能量信息,“需要的导数”为0;如果需要受力,则为1(比如几何优化任务);如果还需要Hessian,则为2(比如freq任务,以及优化或IRC时用了calcfc等情况)。

OutputFile是要在执行xtb.sh时候由这个脚本来生成的,里面记录能量、受力、Hessian等,按照手册的要求格式为:

其中诸如4D20.12是数据格式,稍微懂点Fortran就能明白。按照这个格式产生好OutputFile文件,则Gaussian就会从中读取当前任务需要的信息开展计算。比如如果当前任务只需要能量,则填上能量和偶极矩那一行即可,而如果比如需要Hessian,则整个文件所有信息都得填上。此文件中的偶极矩、极化率对于本文涉及的优化、走IRC、振动频率计算都是不需要的,可以直接填0(但相应地,涉及到这些信息的Gaussian输出,比如偶极矩、freq任务的红外强度等也将都为0)。

要想将Gaussian和xtb联用,关键就是要恰当编写xtb.sh,使得这个脚本可以基于InputFile里的信息去调用xtb计算出当前结构下的能量、受力、Hessian,并按照Gaussian要求的格式转化出OutputFile文件。下一节就介绍怎么实现。


3 Gaussian与xtb的接口xtb.sh的编写

笔者写好的xtb.sh文件在本文一开始提到的压缩包里有。这是bash shell脚本,本节解释一下脚本内容,需要懂得一些Linux命令和shell编程知识才能完全理解(值得一提的是,写这种脚本并非必须用bash shell,也并非必须是Linux环境才能用external功能。比如在Windows下也完全可以写成.bat脚本,例如Windows版Gaussian调用NBO6就是通过external来调用NBO6的.bat文件实现的)。

xtb.sh文件里$2、$3分别对应于xtb.sh接收到的第2、第3个参数,也即InputFile和OutputFile文件名。脚本首先用
read atoms derivs charge spin < $2
从InputFile中把原子数、需要的导数、电荷、自旋多重度分别读到atoms、derivs、charge、spin四个变量里,然后用以下命令构建一个mol.tmp文件
cat >> mol.tmp <<EOF
$atoms

$(sed -n 2,$(($atoms+1))p < $2 | cut -c 1-72)
EOF
这个mol.tmp文件是xyz文件的雏形,还需要做两个处理才能变成xyz格式文件,一方面是把从InputFile读过来的元素序号转化成元素名,另一方面是把读过来的以Bohr为单位的坐标转化成埃。为此,笔者写了个Fortran小程序genxyz.f90,编译好的可执行文件是压缩包里的genxyz。这个小程序稍微懂点编程的人都能看懂,就不解释了。xtb.sh以下两行就是调用genxyz把当前目录下的mol.tmp转化为mol.xyz
./genxyz
rm -f mol.tmp

之后,脚本根据读入的自旋多重度,将之减1,算出来-uhf后面的参数,然后根据当前任务需要的导数信息,来判断在调用xtb时是用-grad还是用-hess
uhf=`echo "$spin-1" | bc` #nalpha-nbeta
if [ $derivs == "2" ] ; then
 echo "Running: xtb mol.xyz -chrg $charge -uhf $uhf -hess > xtbout"
 xtb mol.xyz -chrg $charge -uhf $uhf -hess > xtbout
elif [ $derivs == "1" ] ; then
 echo "Running: xtb mol.xyz -chrg $charge -uhf $uhf -grad > xtbout"
 xtb mol.xyz -chrg $charge -uhf $uhf -grad > xtbout
fi

最后,xtb.sh如下调用笔者自编的extderi程序产生OutputFile文件。对应的源代码extderi.f90也很简单,就不解释了
./extderi $3 $atoms $derivs
extderi会读取xtb运行后在当前目录下产生的energy、gradient、hessian文件,分别提取能量、受力、Hessian信息,输出到OutputFile中。传递给extderi的三个参数$3、$atoms、$derivs分别告诉这个程序要产生的OutputFile文件的文件名是什么、总共多少原子、要读取/写入哪些导数信息。

xtb.sh中还用了一些rm -f命令,用来删除xtb产生的各种文件,确保Gaussian运算后当前目录不残留多余的文件。


4 Gaussian与xtb联用搜索过渡态、做振动分析、产生IRC应用示例

我们这里将Gaussian与xtb联用,搜索一下下图所示的环丙基卡宾异构化过程的过渡态

对应的Gaussian输入文件是本文压缩包里的TS.gjf,开头两行如下
%chk=mol.chk
#P opt(ts,calcfc,noeigen,nomicro) external='./xtb.sh'

可见这里用的是常用的opt=TS方式搜索过渡态,因此输入文件里的结构应当是这个与实际过渡态比较像的过渡态初猜结构。opt里必须写nomicro,否则Gaussian在优化的时候会试图调用分子力学的optimizer去搞,达不到我们的目的。这里我们刻意保留了chk文件,因为这样的话之后做IRC、freq任务就可以直接用geom=allcheck从chk文件中读取已经优化好的过渡态结构来计算了。注意对于当前情况,不能在这里直接写opt freq,必须把opt和freq拆成两步做才行,否则freq任务会出错。另外,由于当前任务的能量、导数都调用xtb来算了,因此理论方法和基组就不需要写了。

我们确保机子里已经装好Gaussian了,xtb也已经配置好了从而可以直接通过xtb命令调用了,然后把TS.gjf、genxyz和extderi都放到当前目录下,运行诸如g09 < TS.gjf |tee TS.out,就开始计算了。

这里特别强调一点,如果你的Gaussian的Default.route里设的默认是并行做Gaussian计算的话,强烈建议改为串行计算(如果不想动这个文件就在.gjf文件里设%nproc=1),否则在Gaussian通过external方式调用外部脚本期间会造成很高无意义的CPU资源消耗,导致总耗时增加。

此任务收敛很顺利,13步就收敛了。找到的过渡态精度如何?下图上半部分是将Gaussian+xtb找到的过渡态(白线)与B3LYP/TZVP下找到的过渡态(红线)放在一起进行对比,已经按照《在VMD中计算RMSD衡量两个结构间的几何偏差》(http://sobereva.com/290)文中的做法将两个结构进行了Align使之最大程度匹配。下图下半部分是过渡态结构的球棍图,便于读者看清楚结构。

 

从上面的对比来看,GFN-xTB方法当初虽然没有特意考虑过渡态问题,但是优化过渡态的结果和较好精度的B3LYP/TZVP比,误差基本可以接受,至少定性正确。大家可以用Gaussian+xtb尝试用不同初猜搜索过渡态,等搜出来一个看着基本合理的过渡态,再用DFT去进一步优化(当然,Gaussian+xtb不是干这个的唯一选择,笔者也尝试了用PM7半经验方法优化这个过渡态,结果也一样定性正确,至于和Gaussian+xtb给出的结构孰优孰劣,从相对于B3LYP/TZVP结构的RMSD偏差上看半斤八两)

接下来再做一下振动分析,看看虚频情况。虽然xtb直接就能做振动分析给出振动频率,但是不便于观看振动模式,而通过Gaussian+xtb联用,结果就可以直接用gview看了。输入文件是本文文件包里的freq.gjf,内容只有两行,为
%chk=mol.chk
#P freq geom=allcheck external='./xtb.sh'
运行之,结果是压缩包里的freq.out。虚频模式的正则矢量如下

从振动动画上看,过渡态确实找对了。虚频大小是732cm-1,而B3LYP/TZVP下是686.7cm-1,PM7下是843.1cm-1,可见Gaussian+xtb的结果合理,而且误差比PM7明显更小。

最后,我们再走一下IRC。输入文件是本文文件包里的IRC.gjf,内容为
%oldchk=mol.chk
#P IRC(maxpoints=20,calcfc) geom=allcheck external='./xtb.sh'
这里用%oldchk是避免IRC任务改写之前的chk。用Gaussian执行之。gview看到的IRC如下

可见IRC曲线很光滑,而且所有点的结构都正常,证明Gaussian和xtb联用很成功。

虽然本文只测试了一个体系,但至少证明Gaussian+xtb用来粗略研究化学反应是充分可行的,值得在实际研究中广泛使用。不过对于普通有机体系,这种做法比起直接用Gaussian自带的PM6/PM7优势不显著。但碰到略诡异体系,预感半经验方法可能连定性正确结果也给不出的时候,则十分建议改用Gaussian+xtb。

最后提醒一下,虽然GFN-xTB往往很不错,但也别以为它的普适性和精度能和一般的DFT计算抗衡。例如优化Li2,M06-2X/def2-TZVP下结果是2.7064埃,实验值是2.6729埃,但目前版本xtb(对应预发布的GFN2-xTB)优化出来是2.2991埃,误差还是不小的,尽管比PM7优化出来的1.8106埃已经强得多了。

评论已关闭