详谈使用Gaussian做势能面扫描

详谈使用Gaussian做势能面扫描

文/Sobereva@北京科音  2019-Apr-1


目录
1 基本知识
2 用Gaussian做刚性扫描的方法和实例
    2.1 基础知识
    2.2 刚性扫描实例1:扫描乙醇的羟基的二面角
    2.3 刚性扫描实例2:扫描Li+...苯之间的距离
    2.4 刚性扫描实例3:扫描氟化氢键长获得解离曲线
    2.5 刚性扫描实例4:同时扫描水分子的键长与键角
3 使用dimerscan结合Gaussian做二聚体的刚性扫描
    3.1 实例:苯酚二聚体的氢键距离扫描
    3.2 实例:主-客体复合物的接触距离扫描
4 使用gentor结合Gaussian做二面角的刚性扫描
5 用Gaussian做柔性扫描的方法和实例
    5.1 基本用法
    5.2 柔性扫描的一些相关问题
    5.3 柔性扫描实例1:柔性扫描C2H4ClBr的二面角
    5.4 柔性扫描实例2:乙酰和甲胺封闭的丙氨酸(ACE-ALA-NME)的构象搜索
    5.5 柔性扫描实例3:辅助寻找乙醇脱水的过渡态
    5.6 柔性扫描实例4:把环丁烷拉成丁烷双自由基
6 通过广义化内坐标(GIC)进行柔性扫描
    6.1 GIC扫描实例1:水+氮气二聚体的几何中心距离扫描
    6.2 GIC扫描实例2:令1,3-丁二烯两个亚甲基同步旋转的扫描
7 总结&其它

本文将介绍势能面扫描的概念,通过最常用的量子化学程序Gaussian演示怎么实现各种类型的势能面扫描,并且介绍利用笔者开发的gentor和dimerscan程序结合Gaussian来实现一些只靠Gaussian实现不了的特殊的扫描。相信笔者读完本文后会对势能面扫描有十分全面的了解。在北京科音(http://www.keinsci.com)开办的初级/基础量子化学培训班里会对这个主题讲得更充分并给出更多例子。


1 基本知识

势能面扫描用来考察体系能量随着一个或多个几何变量的改变而发生的改变。势能面扫描有很多实际用处,比如
• 构造完整势能面或者势能面的子空间。之后基于此可以跑量子动力学、计算反应速率常数时考虑多维隧道效应等等
• 产生力场参数。一般通过势能面拟合实现
• 帮助确定搜索过渡态适合的初猜结构
• 求解振动问题(振动能级、振动平均结构等)。相关信息参考《Molcas的计算双原子分子光谱常数的模块vibrot使用简介》(http://sobereva.com/372)、《谈谈温度、压力、同位素设定对量子化学计算结果产生的影响(http://sobereva.com/423)。
• 寻找能量较低或最低的构象。参考《gentor:扫描方式做分子构象搜索的便捷工具》(http://bbs.keinsci.com/thread-2388-1-1.html
• 帮助弄清楚反应机理、键的解离等化学上感兴趣的过程
• 考察体系电子结构随几何结构特定方式的改变而发生的变化。见比如《制作动画分析电子结构特征》(http://sobereva.com/86)。

扫描分为刚性扫描(rigid scan)和柔性扫描(relaxed scan)。刚性扫描指的是让程序按照你指定的扫描设定依次改变体系几何结构,结构每变一次就算一次单点能,而不被扫描的那些几何变量还是保持在初始结构的值。而在柔性扫描中,每次按照你的要求而改变特定几何参数时,做的不是单点计算而是限制性优化,即没有被扫描的变量都会被优化以使得体系能量尽可能低,相当于允许这些变量自发地弛豫(relax)。我们平时说扫描的时候一般默认指前者。因此如果你想告诉别人你做的是柔性扫描的话,那么绝对不要简单地说成“扫描”或scan,免得造成误会。

下面,笔者就结合最常用的量子化学程序Gaussian对势能面扫描的操作进行介绍,并且举各种势能面扫描的例子。所有计算使用Gaussian16 A.03版完成,GaussView用的是6.0.16版,同时会利用到的xyz2QC和gentor程序来自于molclus 1.8版程序包(可在此免费获取http://www.keinsci.com/research/molclus.html)。

本文例子的输入文件都可以在此文件包里找到:http://sobereva.com/attach/474/file.rar。文件名在例子中都提到了。输出文件也给了。


2 用Gaussian做刚性扫描的方法和实例

2.1 基础知识

在Gaussian里做刚性扫描很简单,关键词里写上scan,然后对要扫的变量进行定义即可。键长、键角、二面角都可以扫描,被扫描的几何坐标必须以变量形式表达。比如某个键长变量原本定义是B1= 1.3,那么把这个变量设置改为B1= 1.3 10 0.1,就代表以1.3埃为初始值扫描10步,每步增加0.1埃,最终扫到1.3+10*0.1=2.3埃。由于初始结构也会被计算一次,所以总共会计算11个单点能。算完之后,可以从输出文件末尾读取信息汇总,也可以把输出文件用GaussView打开观看扫描过程每一步的结构变化,还可以在Result - Scan里图形化显示能量在扫描过程中的变化。

还值得一提的是,扫描的每一步用的初猜波函数会自动用上一步收敛的波函数,主要考虑是相邻两步之间结构变化不是特别大,所以这样获得初猜波函数比重新产生初猜波函数通常更容易令SCF收敛,从而节约时间,通常还能使得上一步的电子结构能够承接下去。如果你做扫描的时候,发现一开始能收敛,但是扫描中途有的点SCF不收敛,或者收敛到奇怪的情况(体现在能量出现突然波动),不妨将扫描步长改小试试。

下面我们来看一些例子,被考察的问题由简单到复杂。


2.2 刚性扫描实例1:扫描甲醇的羟基的二面角

这个简单例子的目的是考察一下甲醇的羟基二面角旋转过程的能量变化。做扫描之前,我们先在B3LYP/6-31G*下对甲醇做一下几何优化,优化后是下图的结构,两张图是不同的视角。可见当前结构是Cs点群:

用GaussView把优化后的结构保存为methanol.gjf文件,在保存的界面里记得把Write Cartesians复选框取消掉,以使得保存出的坐标是以内坐标方式记录的。

我们打算在B3LYP/6-31G*下做扫描,想让C1-O5键转180度,使得上图的右图中H6从H2的反方向转到与H2重合。因此,我们得从methanol.gjf里去看内坐标的定义,看让哪个二面角变量的改变可以使得H6绕着C1-O5轴旋转。在gjf里可见
 C             
 H                  1            B1
 H                  1            B2    2            A1
 H                  1            B3    2            A2    3            D1
 O                  1            B4    2            A3    3            D2
 H                  5            B5    1            A4    2            D3
显然,D3对应H6-O5-C1-H2二面角,因此如果让D3变化,即绕着C1-O5转,就可以达到目的了。我们让这个变量每10度变化一次,因此总步数应当是180/10=18步。遂把gjf文件改成下面这样,对应methanol_scan.gjf:

# b3lyp/6-31g(d) scan nosymm
[空行]
Title Card Required
[空行]
0 1
 C             
 H                  1            B1
 H                  1            B2    2            A1
 H                  1            B3    2            A2    3            D1    0
 O                  1            B4    2            A3    3            D2    0
 H                  5            B5    1            A4    2            D3    0
[空行]
   B1             1.09337920
   B2             1.10125884
   B3             1.10125884
   B4             1.41863900
   B5             0.96872260
   A1           108.05956655
   A2           108.05956655
   A3           106.69693538
   A4           107.66742719
   D1           117.09085412
   D2          -121.45457294
   D3           180.00000000 18 10.

注意,凡是程序需要读入浮点数的地方,绝对不能写成整数,因此此例步长必须写10.或者10.0而不能写10。另外,做刚性扫描的过程中,如果涉及到点群的改变,不加nosymm关键词的话有时会出问题,而且看到的扫描轨迹可能不连续。一开始结构是Cs点群,二面角改变一下就成C1点群了,所以我们这里用了nosymm。关于nosymm的更多说明看《谈谈Gaussian中的对称性与nosymm关键词的使用》(http://sobereva.com/297)。

这个体系小、计算级别低,很快就算完了。计算过程中途输出的信息不用管,直接把输出文件拉到末尾,可以看到下面的信息,是扫描过程的汇总。D3是指被扫描的那个变量,当前用的B3LYP属于SCF一类方法,所以每个点的能量显示在SCF那一列下面,单位是Hartree。
 Summary of the potential surface scan:
   N      D3          SCF   
 ----  ---------  -----------
    1   180.0000   -115.71441
    2   190.0000   -115.71425
    3   200.0000   -115.71382
    4   210.0000   -115.71320
...略

把methanol_scan.out拖到GaussView里,在窗口左上角可以看到帧号。当前扫描18步,初始结构也算一个结构,因此一共19帧,第一帧对应输入文件里的结构。选择Results - Scan,可以看到扫描过程的能量变化,如下图所示。点击其中一个点,图形窗口就会切换到相应的帧。在窗口下方可以看到这个点对应的能量和被扫描的坐标的当前值。

从扫描曲线上,可以看到有两个极小点(对于这条曲线而非整个势能面而言),一个是D3=180度,正好对应初始结构,这显然是能量最低的结构。另一个局部极小点在D3大约300度的位置,结构如上图左边所示。而当H6转到与H2沿着C-O键轴看正好重合的时候,能量达到了全局最大。

通过这个图中的极大点和相邻极小点位置,我们可以粗略估计旋转势垒。但是要想精确计算的话,应当用这个图中极小点和极大点的坐标分别作为初猜结构去准确优化极小点和过渡态结构(找过渡态参考《简谈Gaussian里找过渡态的关键词opt=TS和QST2、3》http://sobereva.com/460),然后再用更高级别算单点能再求差。

做刚性扫描比较烦人的一个情况是想扫描的那个坐标在GaussView直接保存的gjf文件没有恰好对应的变量。此时要么想办法修改体系坐标书写方式使得被扫描的坐标能通过某个几何变量表示;要么手动在GaussView里每修改一次坐标就保存一个输入文件,扫多少步就产生多少个输入文件,然后批量运行、批量提取单点能。Linux下批量执行的方法看《使用Gaussian时的几个实用脚本和命令》(http://sobereva.com/258)。


2.3 刚性扫描实例2:扫描Li+...苯之间的距离

这个例子中我们要扫描Li+与苯的苯环中心在垂直于苯环方向的距离。这个目的有两种实现方式,第一种是基于笛卡尔坐标来设定扫描方式,下面说一下。我们先在GaussView里画一个苯,做对称化成为D6h点群,用便宜的B3LYP/6-31G*优化一下,然后以笛卡尔坐标方式保存为gjf文件。此文件里,苯正好在Z=0的XY笛卡尔平面上,苯环中心就是(0,0,0)位置。假设我们想扫描Li+,使之从距离苯环中心5埃处逐渐接近苯环中心,共10步,每一步移动0.3埃,最终移动到相距2埃的位置,我们应当把输入文件写成这样(Ben_Li_Cart.gjf):

# M062X/6-311g(d) scan pop=always
[空行]
Title Card Required
[空行]
1 1
 C                  0.00000000    1.39650157    0.00000000
 C                  1.20940584    0.69825078   -0.00000000
 C                  1.20940584   -0.69825078   -0.00000000
 C                  0.00000000   -1.39650157    0.00000000
 C                 -1.20940584   -0.69825078    0.00000000
 C                 -1.20940584    0.69825078    0.00000000
 H                  0.00000000    2.48320512   -0.00000000
 H                  2.15051871    1.24160256   -0.00000000
 H                  2.15051871   -1.24160256   -0.00000000
 H                 -0.00000000   -2.48320512    0.00000000
 H                 -2.15051871   -1.24160256    0.00000000
 H                 -2.15051871    1.24160256    0.00000000
 Li                0. 0. Z
[空行]
Z= 5.0 10 -0.3

可见我们手动添加了一个Li原子,X和Y都为0,Z坐标以变量表示,初值为5.0,最后一步距离苯环中心将是5.0+10*(-0.3)=2.0埃。此例为了让结果达到基本定量合理,用了算弱相互作用不错的M06-2X结合比前例稍微大一些的基组。由于扫描过程始终是C6v点群,所以此例我们不用写nosymm关键词,如果写了这个关键词的话由于没法利用对称性加速计算会多花很多时间。

有量子化学常识的人都知道Gaussian里只能设定体系总电荷,没法直接指定Li带的电荷。但Li在实际体系中到底带多少电荷?随着扫描的进行其电荷量如何变化?想了解这个问题,最简单的做法就是在扫描的时候像本例这样同时写pop=always关键词,这样扫描的每一步都会做一次布居分析从而得到原子电荷。

执行上面这个任务,把输出文件用GaussView打开,扫描曲线如下。图应该从右往左看

可见随着Li逐渐接近苯环中心,体系能量显著下降。我们可以直接在GaussView里观看原子电荷的变化,做法是在Scan界面里选择Plots,选Plot Molecular Property,再选Atomic Charge,然后输入13,点OK。之后在能量变化图的下方就可以看到序号为13的Li原子的Mulliken电荷随扫描的变化了,如下所示

由图可见,当Li距离苯很远的时候,Li带的原子电荷接近1.0,是个典型的阳离子。而随着距离接近,苯上的电子就逐渐往Li+上转移,导致在距离较近的时候Li带的原子电荷明显达不到1.0了。注意,Mulliken电荷只是一个比较low的原子电荷模型,存在容易低估离子性、怕弥散函数等问题,相关讨论参见笔者的《原子电荷计算方法的对比》(http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml)一文,在实际研究中用ADCH或NPA电荷说事更好。如果想对电荷转移特征了解更多,还可以用Multiwfn做更细致的布居分析,看Multiwfn手册4.7.0节的例子,以及绘制密度差图,看《使用Multiwfn作电子密度差图》(http://sobereva.com/113)。

Plot Molecular Property界面里还可以绘制别的,比如绘制几何变量随扫描的改变。使用pop=always后,还可以在这个界面里选择绘制偶极矩矢量或其总大小。而如果不用pop=always的话,程序只会在初始结构计算时做布居分析。

这个例子也可以在内坐标下实现扫描,但是需要用到虚原子。做法是在GaussView里显示出苯之后,在苯环中心添加一个虚原子X(把苯环上的六个碳都选中成为黄色,在Builder窗口里选X原子,然后点右键选Builder - Place Fragment at Centroid...),然后再在虚原子上加个氢,把X-H的键轴调成与苯环垂直,然后把氢替换成Li原子。之后保存内坐标形式的gjf文件,会发现其中有一个变量正好对应虚原子和Li原子之间的距离,因此对这个距离扫描即可。输入文件是Ben_Li_Zmat.gjf。

实际研究中,可能环是斜着的,要扫描某个原子与它的环中心的距离。此时既可以用上面提的内坐标方式扫描,也可以用笛卡尔坐标方式扫描。对于后者,由于Gaussian里在刚性扫描的时候没法让某个原子按照某个自定义矢量来移动,因此需要先旋转体系,让环平面恰好平行于某个笛卡尔平面,比如平行于XY,这样才能通过扫描Z坐标实现目的。按照《调节平面分子间距的方法》(http://sobereva.com/178)里的“方法二”,借助笔者自写的VMD程序的脚本,可以轻易让某个环平面恰好平行XY平面。

值得一提的是,在扫描过程中如果能进行深入细致的电子结构分析,可以获得远比结构+能量丰富得多的信息,让研究文章充实起来。想实现这点,可以用《产生Gaussian的IRC和SCAN任务每个点的波函数文件的工具》(http://sobereva.com/199)一文介绍的笔者写的SCANsplit工具载入刚性扫描的输出文件,将扫描过程中每个结构转化为单点任务输入文件,批量执行后就有了记录每个点波函数的fch或wfn文件,之后再通过脚本和一些Linux下的命令就可以自动调用Multiwfn考察电荷分布、成键等电子结构方面的信息随扫描过程的变化,详见《通过键级曲线和ELF/LOL/RDG等值面动画研究化学反应过程》(http://sobereva.com/200),可以讨论的手段见《Multiwfn支持的分析化学键的方法一览》(http://sobereva.com/471)。


2.4 刚性扫描实例3:扫描氟化氢键长获得解离曲线

这一节我们通过DFT方法研究氟化氢分子的H-F键断裂过程的势能曲线。这个曲线我之前在《浅谈为什么优化和振动分析不需要用大基组》(http://sobereva.com/387)给出过,当时对不同级别扫描出的曲线进行了对比。H-F键的断裂过程中体系会由闭壳层分子逐渐变成H自由基和F自由基。当键长在平衡距离附近的时候,体系是闭壳层状态,计算没有什么特殊的,而距离拉远后,超过所谓的“不稳定点”对应的键长,整个体系就成了双自由基状态,属于自旋极化单重态,这时对于DFT计算而言需要做对称破缺计算。如果对这点不懂,一定要看《谈谈片段组合波函数与自旋极化单重态》(http://sobereva.com/82)。因此,扫描这种共价键均裂解离曲线不是看起来那么简单。而如果你扫描的是NaCl这样的离子键异裂成Na+和Cl-的解离曲线,那就不需要考虑这些问题了,因为解离后变成的Na+和Cl-都是闭壳层状态。

对于这个体系的扫描,我的建议是这样:
(1)先在H和F距离很远、键完全断裂的情况下产生对称破缺波函数,对应双自由基状态
(2)用非限制性开壳层形式(U)做扫描,让H-F距离逐渐减小,第一步的初猜波函数直接读取(1)收敛的波函数。由于之后每一步都会用上一步收敛的波函数做初猜,所以一开始的波函数的对称破缺状态会一直维持下去。当距离减小到已经小于不稳定点的时候,由于此时基态波函数是闭壳层波函数,因此从此开始每一步的结果和用限制性闭壳层(R)形式计算的结果将会完全一样。

第1步的输入文件是下面这样(HF_broken.gjf),对应键长4埃。用U明确指定做非限制性计算,用guess=mix是为了打破初猜波函数的对称性,但是直接用guess=mix未必能收敛到真正的对应当前结构下基态的对称破缺波函数,因此又加上了stable=opt来对收敛的波函数进行检测,如果发现不稳定,会自动尝试优化出最稳定的波函数。

%chk=C:\HF.chk
# UB3LYP/def2TZVP guess=mix stable=opt
[空行]
Title Card Required
[空行]
0 1
 F
 H      1 B1
[空行]
B1 4.0

第2步的输入文件如下(HF_scan.gjf),从4.0埃扫到4.0-0.06*55=0.7埃。涉及对称破缺的计算加上nosymm会比较稳妥,所以这里用了nosymm。

%oldchk=C:\HF.chk
# UB3LYP/def2TZVP scan guess=read nosymm
[空行]
Title Card Required
[空行]
0 1
 F
 H      1 B1
[空行]
B1 4.0 55 -0.06

最终扫出来的曲线如下,非常理想。如果扫描出的这种曲线有一些明显的突跃,通常是那个点恰好没有收敛到与周围的点特征基本相同的波函数所致,此时可以尝试把扫描步长减小再试,这样有助于保持波函数的连续性,毕竟每个点的初猜自动用上一个点的收敛的波函数,二者结构相差越小则波函数的状态越容易继承下去。

扫描共价键的解离曲线还有个做法,只需要一步,就是使用比如# UM062X/TZVP guess(always,mix) scan nosymm这种关键词,此时让键长从大变小来扫描,还是让键长从小变大来扫描,其实都一样,因为这里用了always,代表每一步重新产生初猜波函数;而且由于用了mix来使得每一步都试图得到对称破缺态,如果对称破缺态更稳定,那么这个点就收敛到对称破缺态,而如果闭壳层态更稳定,那么就自动收敛到对应闭壳层的波函数。大家可以尝试对乙烷的C-C键用这个关键词扫描,会发现可以得到正确的解离曲线,扫到最后相当于两个甲基自由基。但是这套关键词用于扫描上面的氟化氢体系的话会发现行不通,会扫得乱七八糟,因为对这个体系的很多点,光靠guess=mix产生的初猜波函数收敛不到实际基态波函数上。不信的话可以看一下HF_scan_mix_always.out的曲线,就是以这种方式算的,其最后一个点SCF没有收敛,只看其它点的话,会看到能量曲线不合理,因为拉远后能量变化没有趋于水平。虽然从输出文件中看到超过不稳定点的结构的<S**2>确实不为0,因此确实得到了对称破缺态,但是和之前我们扫描出来的输出文件里相应的点<S**2>不符(明显偏小)。有兴趣的话读者可以对诸如键长为3.5埃的时候用UB3LYP/def2TZVP guess(mix,always) nosymm产生fch文件,用Multiwfn看看自旋密度分布了解是怎么回事(参见《谈谈自旋密度、自旋布居以及在Multiwfn中的绘制和计算》http://sobereva.com/353),你会发现,并不是如期望的在H和F上各有单电子且自旋相反,而是alpha和beta单电子同时出现在了F上面!

可以计算共价键解离曲线的方法极多。使用恰当的泛函,通过对称破缺方式计算通常就可以给出不错的解离曲线。但如果要求更高,可以用UCCSD(T)。而对于牵扯多重键断裂的问题,情况复杂、静态相关很强,通常需要考虑用较复杂但是普适性强的多参考方法如CASPT2以确保得到靠谱的结果。


2.5 刚性扫描实例4:同时扫描水分子的键长与键角

这个例子演示二维扫描。对水分子,令两个O-H键长同时从0.92扫到1.02埃,每步0.01埃,共10步。与此同时,让键角从95度扫到115度,每步2埃,共10步。因此扫描任务中要算的单点能的次数是(10+1)*(10+1)=121个。显然,随着被扫描的变量数增加,刚性扫描的耗时呈几何式增加。对某个刚性扫描任务要想估计能不能扫得动,你可以算一个单点看看耗时,然后乘以要算的点数。

当前这个任务的输入文件如下(H2O_2Dscan.gjf)
# B3LYP/def2SVP scan
[空行]
Title Card Required
[空行]
0 1
 O             
 H                  1    B1
 H                  1    B1    2  A1
[空行]
B1 0.92 10 0.01
A1 95.0 10 2.0

算完之后,可以把输出文件末尾的汇总数据用Sigmaplot、Surfer、Origin之类绘制成曲面图或者填色图或者等值线图便于考察。如果将输出文件载入GaussView,也可以看到有121帧,并可以观看扫描轨迹。如果进入Results - Scan,会看到下图,曲面图的纵坐标对应电子能量,在下方还有投影的等值线图。

点击曲面上的小白点使之成为绿色,可以使状态栏下方显示相应的点对应的几何变量的数值和能量,图形窗口里的结构也会切换到对应的帧。由上图可见,所有扫描的点里能量最低的是B1=0.97埃、A1=103度的那个点,比起其它点更接近于水分子的势能面最低点。

类似地,我们还可以用Gaussian同时扫描更多的坐标,但超过两个坐标时,结果就没法直接用GaussView绘图直观展现了。


3 使用dimerscan结合Gaussian做二聚体的刚性扫描

借助笔者写的dimerscan和xyz2QC程序,我们可以实现很多没法直接用Gaussian的scan设定实现的二聚体的刚性扫描。下面通过两个例子来体现这点。dimerscan可以从《考察SAPT能量分解的能量项随分子二聚体间距变化的简单方法》(http://sobereva.com/469)页面里下载,xyz2QC是笔者开发的团簇构型与分子构象搜索molclus程序包中带的工具,可以在molclus主页下载:http://www.keinsci.com/research/molclus.html

3.1 实例:苯酚二聚体的氢键距离扫描

经常研究分子间相互作用的人,往往会遇到一个问题就是想对二聚体之间做刚性扫描,但是GaussView产生的gjf文件里没有正好对应于想扫描的坐标的几何变量。比如下面的苯酚二聚体

我们想考察体系能量随氢键距离变化而发生的变化,从而得到不同距离下两个单体间的弱相互作用能,因此我们应该扫描H9-O20或者O7-O20键长。但不幸的是,GaussView里保存的内坐标形式的输入文件中(phenoldimer.gjf),O20的键长项是参考C14来定义的,因此只使用Gaussian的话没法达到我们的扫描目的。虽然如后文所述,用柔性扫描的话想怎么定义都可以,但是柔性扫描不仅昂贵,而且扫描过程中其它坐标都可能显著变化,和我们期望的研究目的不符。像这种情况,我们要扫描氢键距离,就得用dimerscan结合xyz2QC程序来实现。下面提到的文件都在本文文件包里rigid\phenoldimer目录下。

我们先用GaussView打开phenoldimer.gjf,保存成笛卡尔坐标形式的gjf。将gjf里原子坐标以外的部分都删掉,在坐标前面插入两行,第一行是两个单体各自的原子数,第二行是计算两个单体时分别用的电荷和自旋多重度(对当前研究无影响),然后把这个文件保存为dimer.txt。此文件当前内容的形式如下

13 13
0 1 0 1
[苯酚单体1的坐标]
[苯酚单体2的坐标]

假设我们想扫描H9-O20,从1.7埃扫20步,每步增加0.1埃。将dimer.txt放到dimerscan.rar解压后的目录中,启动dimerscan,输入dimer.txt的路径,然后依次输入
9,20
1.7   //初始距离
20    //扫20步
0.1   //每步伸长0.1埃
最后按一次回车退出。此时当前目录下出现了一堆.inp文件,不用管,那是给PSI4程序做SAPT用的。当前目录下还出现了scan.xyz,里面每一帧对应扫描过程的一个结构。如果想检验下一生成的对不对,可以把这个文件拖到VMD程序(http://www.ks.uiuc.edu/Research/vmd/)里播放动画确认一下,由动画可见确实生成得没错:

把molclus程序包中的template.gjf里的关键词改为要计算每个点用的级别,比如我们想用比较便宜的PM6-D3方法实现苯酚二聚体扫描,就把这个文件里的关键词设为# PM6D3。然后启动molclus程序包中的xyz2QC程序,选择1 Generate multi-step Gaussian input file,输入scan.xyz的路径,然后直接按回车代表考虑所有帧,再按回车退出。当前目录下马上出现了Gaussian.gjf文件,此文件是Gaussian的多步任务文件,每一步是对每个结构算一次单点。

用Gaussian运行Gaussian.gjf,用ultraedit或者Linux下的grep命令将里面所有带有E(RPM6D3)字样的行都提取到一个文本文件里,会发现一共有21行(初始结构+扫描的20个结构)。然后把多余的列都删掉只保留能量的列,导入到Origin里,把扫描的距离信息也插入到里面作为一列,然后用Origin作图,如下所示(原始文件是phenoldimer.opj),完全与我们预期的相符。

注:当二聚体结构中不同单体的原子顺序存在相互交错的时候,需要先处理一下,确保单体里的原子序号是连着的,否则dimerscan没法用。处理过程是用GaussView打开二聚体结构文件,在其中一个单体的原子上点右键选Select Fragment of ...使这个单体变成黄色,然后按Ctrl+X复制到剪切板,创建一个新窗口,再按Ctrl+Shift+V将之粘贴进去,保存成笛卡尔坐标形式的输入文件1.gjf。之前的窗口还剩另一个单体,保存为2.gjf。最后将2.gjf里的坐标拼接到1.gjf的坐标后头去。


3.2 实例:主-客体复合物的接触距离扫描

此例我们想对下图所示的主-客体复合物在红线方向上进行扫描,看这个富勒烯脱离主体分子过程的能量变化曲线。复合物结构是本文文件包里host-guest目录下的host-guest.gjf。

做这个体系的刚性扫描也没法直接用Gaussian实现,不仅被扫描的两端恰好没有原子出现,而且即便有原子,一般也不会恰好有对应的变量可供扫描。对这种情况,我们还是可以用dimerscan+xyz2QC的组合来实现,只不过用之前我们先得在扫描的两端位置添加虚原子。

我们先在GaussView里将客体分子剪切并粘贴到新窗口里,然后将主体和客体分子分别保存成host.gjf和guest.gjf。然后打开主体分子,在扫描的位点增加一个虚原子,然后保存gjf。对客体分子也这样增加虚原子并保存。虚原子位置如下所示

然后把host.gjf和guest.gjf里的坐标合并在一起,改成下面的格式,并保存为dimer.txt。目前包含虚原子X在内主体和客体分子分别有89和61个原子。
89 61
0 1 0 1
[主体分子的坐标]
[客体分子的坐标]

将dimer.txt拷到dimerscan目录下,启动dimerscan程序并输入
dimer.txt
89,150  //新添加的两个虚原子在整体中的序号,都是每个单体最后一个原子,因此第二个虚原子序号是89+61=150
2.8  //初始距离
30   //30步
0.15  //每步增加0.15埃

用VMD播放一下新产生的scan.xyz,会看到下面的动画

还是按照上一节的方法,通过xyz2QC程序将scan.xyz转换成多任务的Gaussian输入文件Gaussian.gjf,使其中每个子任务都对应于用PM6-D3计算单点。用Gaussian执行此任务,然后提取SCF Done的数据进行作图,会看到下图的结果:

值得一提的是,一般绘制这种图都应该用极小点结构作为纵轴零点,令每个点的能量都减去这个能量,使得图上纵坐标体现的是相对于极小点的能量变化,此时纵坐标才有化学意义。


4 使用gentor结合Gaussian做二面角的刚性扫描

往往我们想刚性扫描二面角,并且扫描的时候让整个基团连带地转动。比如下面这个体系,我们想看看绕