使用EzReson进行化学共振分析（3）：实例：1,3-丁二烯的共振分析

在前两篇推送中，我们介绍了定量的共振理论和EzReson的安装，从本篇开始，我们通过实例展示用EzReson程序进行化学共振分析的方法。

一. 进行WFRT分析的一般步骤

 (1) 对给定分子进行DFT单点计算（同时要调用NBO程序），以得到包含分子轨道等信息的fchk文件和包含NAO系数矩阵的33文件。如果是基于简单的Hückel理论分析，则无需DFT计算，自然也无需fchk和33文件。

(2) 根据所要进行的共振分析任务类型，准备EzReson输入文件。

(3) 运行EzReson，得到输出结果。

具体的操作见下文的实例部分。

二、具体实例

这里我们介绍几个例子：分析1,3-丁二烯的共振结构、用共振论研究一下Diels–Alder环加成反应、用共振分析研究氢键的强度。这些实例的所有输入和输出文件都能在EzReson软件包的tests文件夹里找到。限于篇幅，本篇先介绍第一个应用。

1. 1,3-丁二烯的共振分析（简单Hückel框架下）

我们首先看一个非常简单的例子，基于简单Hückel分子轨道（HMO）理论对1,3-丁二烯分子进行共振分析。用GaussView等软件搭建1,3-丁二烯分子的结构，导出为gjf原子坐标文件（需保存成直角坐标）。不妨取名该文件为butadiene.gjf，其具体内容如下：

# B3LYP/6-31G* OPT

1,3-butadiene

0 1
 C   -0.44433094   -1.18998977    0.00000000
 C    0.23094337   -0.01501247    0.00000000
 C    1.77094337   -0.01501247   -0.00000000
 C    2.45318677    1.15593201   -0.00000000
 H    0.08883281   -2.11769469    0.00000000
 H   -1.51433094   -1.18998977    0.00000000
 H   -0.30222038    0.91269245    0.00000000
 H    2.30410712   -0.94271739    0.00000000
 H    1.92554430    2.08678830    0.00000000
 H    3.52316789    1.14957471   -0.00000000

需要注意的是，由于简单HMO计算只考虑碳原子，EzReson输出结果中将不包含氢原子，所以最好把所有碳原子都放到最前面来排列，这样才能保证EzReson输出的原子序号与gjf文件中的碳原子序号相一致，否则就不太方便解读结果。

然后准备一个基于HMO的WFRT分析的输入文件，如下列名为butadiene_wfrt-huckel.in的文件：

File = butadiene
Job = WFRT
Atoms = 1 2 3 4
Huckel = TRUE

其中，File关键词指定文件的关键名为butadiene（注意不用加引号），这样EzReson就知道该从butadiene.gjf文件载入分子的原子坐标。Job关键词指定当前的计算任务类型是进行WFRT分析。Atoms关键词指定参与共振分析的原子是序号为1、2、3和4的碳原子。Huckel关键词指定当前的分析是基于简单HMO理论的，当Huckel取值为FALSE或默认不指定Huckel关键词时，EzReson将进行基于DFT的共振分析。

现在一切准备就绪了。在终端执行：

ezreson butadiene_wfrt-huckel.in > butadiene_wfrt-huckel.out

大约不到1秒的时间，计算就结束了。计算结果都输出在文件butadiene_wfrt-huckel.out中，我们来解读一下其基本内容：

这部分说明程序识别出该共振体系为4c-4e体系，共有21个Lewis共振结构。

在这里，程序确定出所有10个Lewis单元轨道（LUOs），其中前6个对应于共价键（BD，最后两个数字是形成该键的两个原子的序号），后4个LUO对应于孤对电子（LP，最后一个数字是其所在原子的序号）。

这里列出了程序穷举出的所有可能的21个Lewis共振结构，第二列Projection是各共振结构波函数在真实分子的波函数上的投影{p_i}，第三列给出了构成各共振结构波函数行列式的LUO序号，最后一列信息是各共振结构的具体Lewis结构式，采用的是一套独创的表示法：数字是各原子序号，冒号（:）表示孤对电子，横线（-）表示共价键，各共价键和/或孤对电子之间用空格隔开，请对照图1自行解读一下。

图1. 1,3-丁二烯分子的主要共振结构及所占比例（HMO和DFT框架下的结果分别用非粗体和粗体给出）。最后一个共振结构是违反Rumer规则的结构。

这里程序确定出第16号Lewis结构违反了Rumer规则（见图1中最后一个结构），因此将把它排除出去，只用剩下的20个Lewis结构来展开体系的真实波函数。

以上显示了最终的WFRT求解结果，其中第3列Coefficient给出的是各共振结构对应的展开系数{c_i}；第4列RE是各共振结构的参考能量值，即各共振结构的能量E_i对共振杂化体的能量E（见之前的理论介绍部分）的相对值，在这里是以简单HMO理论中的|β|为单位的；后面几列给出的是各种方案下的共振结构占比{w_i}，如Mulli.指的是Chirgwin–Coulson权重，可以看出整个结果列表是按PWSO权重进行排序的；最后一列是各共振式的具体Lewis结构表示。

如图1所示，含两个C=C的纯共价结构占了77.44%，但其它共振结构也或多或少地有所贡献，这说明了1,3-丁二烯中π电子的共轭离域效应。其中，2,3号碳原子形成双键的Lewis结构（4:/2-3和1:/2-3）是第二重要的共振结构。

最后可以看到，该WFRT分析中，真实波函数的重现度为100%，即波函数的展开是绝对精确的。共轭体系的真实HMO能量是–4.4721 |β|，而共振杂化体的能量也是–4.4721 |β|，说明该共振杂化体就能完全表示真实的共轭体系。该体系的共振能是0.4721 |β|，也就是最稳定的共振结构（1号结构）与真实体系之间的能量差。

*注意事项：

(1)  关于原子序号的排列：

Atoms关键词中指定的原子序号的顺序是有讲究的，因为应用Rumer规则时把所有原子放在一个圆环上来判断是否存在相交的共价键。对于本例的顺式1,3-丁二烯，分子中的原子排布顺序与把它们放在圆环上的排布顺序是一致的，所以1-3/2-4是共价键相交的结构，想要排除掉这个结构可以设Atoms=1,2,3,4；但对于反式1,3-丁二烯，分子中的原子排布顺序与把它们放在圆环上的排布顺序不一致，并且1-4/2-3才是共价键相交的结构，所以应设置Atoms=1,2,4,3才能把它排除掉。更多详细解释见文献DOI:10.1039/D0CP06207C。

(2)  EzReson输入文件中的关键词不区分大小写。也可加注释，放在#号后即可。

2. 1,3-丁二烯的共振分析（DFT框架下）

现在，我们通过更可靠的DFT计算来对1,3-丁二烯进行共振分析。我们首先用Gaussian软件计算上述butadiene.gjf，得到优化后的分子结构。由于目前EzReson只支持fchk文件格式，所以需要使用Gaussian软件来计算。然后，取优化好的几何构型再用Gaussian做以下的单点计算，不妨取文件名为butadiene-sgp.gjf，内容如下：

%chk=butadiene-sgp
# B3LYP/6-31G* POP=NBO6Read

1,3-butadiene  single-point calc.

0 1
  6    1.339919   -0.922552   -0.000000
  6    0.000000   -0.926478    0.000000
  6   -0.894356    0.241893    0.000000
  6   -0.540716    1.534311   -0.000000
  1    1.921683   -0.004033   -0.000000
  1    1.904522   -1.850053   -0.000000
  1   -0.506630   -1.891193    0.000000
  1   -1.957882    0.004699    0.000000
  1    0.497843    1.856036   -0.000000
  1   -1.288622    2.321503    0.000000

$NBO
  NOBOND
  AONAO=W
$END

对这个Gaussian输入文件解释一下：

(1) 第一行是让Gaussian生成chk文件，注意取与gjf文件相同的关键名，即“butadiene-sgp”。

(2) POP=NBO6Read用于调用NBO 6或NBO 7程序。如果没有购买NBO程序，可使用Gaussian自带的NBO 3.1程序（用POP=NBORead关键词），不过这个老版本程序给出的NAO可能没有NBO 5.x以上版本的结果可靠。

(3) 末尾的$NBO…$END中，NOBOND指令是不让NBO程序来寻找自然键轨道（Natural Bond Orbitals，NBOs），这样做一方面可节省时间，另一方面避免了对于某些体系NBO程序会出现死循环的bug；AONAO=W指令是让NBO程序把原始基组AO→NAO的变换矩阵写入扩展名为33的外部文件中。

计算完成后，用Gaussian的formchk工具将得到的butadiene-sgp.chk文件转为butadiene-sgp.fchk文件。至此，让我们再检查一下下列3个文件是否齐全：

• butadiene-sgp.out

• butadiene-sgp.fchk

• butadiene-sgp.33

若在Linux下计算，一般输出文件后缀为.log，此时要将后缀改为.out。对于.33文件，如果使用Gaussian自带的NBO 3.1，则是在工作目录下生成FILE.33文件，需要手动改为butadiene-sgp.33；如果使用NBO 7程序，则默认情况下.33文件是生成在Gaussian的临时文件目录中，名字类似Gau-XXX.33，此时需要将其移动到工作目录下，并命名成上面的形式，也可以在$NBO段落中加入archive file=butadiene-sgp关键词，则可在当前目录下生成butadiene-sgp.33文件。

接着再准备EzReson的输入文件butadiene_wfrt.in如下：

File = butadiene-sgp
Job = WFRT 
LMOs = 14 15
Atoms = 1 2 3 4

这里设置了LMOs为第14和15条轨道来定义共振子系统，因为1,3-丁二烯共有4个π电子，对应2条占据的π轨道，由于π轨道能量高于σ轨道，所以这2条π占据轨道是最高占据轨道，因此我们应该取HOMO–1 和HOMO两条LMO，分别是第14和15条（可从butadiene-sgp.out给出的体系总电子数来确定）。

然后使用ezreson程序运行butadiene_wfrt.in得到输出文件butadiene_wfrt.out。计算得出的各共振式所占比例见图1中粗体所示数值。可见，简单HMO方法给出的共振结构比例与DFT结果还是基本一致的，虽然具体数值有所不同，但各共振结构的重要性顺序都是完全一致的。

另外，可以看到B3LYP/6-31G*级别下WFRT的波函数重现度高达99.89%，表明真实波函数的展开还是相当精确的。