对化合物 Cyclobuta[1,2,3,4]dicyclopentene 的理论分析 (2)

前言

笔者在刚开始学量化时，曾对化合物 Cyclobuta[1,2,3,4]dicyclopentene 的一些性质进行了简要的分析，其中包括了对键合性质、离域性质、磁性质的分析，感兴趣的可以浏览笔者的第一篇博文《对化合物 Cyclobuta[1,2,3,4]dicyclopentene 的理论分析 (1)》。不过笔者在对该结构芳香性的考察时，没有涉及到该结构的 NICS，这属实是一种遗憾。因此，本文将着重介绍 NICS 以及通过 NICS 衡量该结构的芳香性。本文主要参考了 Sob 老师的以下几个博文：

• 《衡量芳香性的方法以及在 Multiwfn 中的计算》
• 《使用 Multiwfn 绘制一维 NICS 曲线并通过积分衡量芳香性》
• 《使用 Multiwfn 巨方便地绘制二维 NICS 平面图考察芳香性》

什么是 NICS

NICS（核独立化学位移）是目前最流行、可接受度最高的衡量芳香性的方法。NICS 原始的定义是：环上重原子的几何中心处的各向同性化学屏蔽值的负值，以 ppm 为单位。对于芳香性体系，一般这个屏蔽值是正值，即 NICS 为负，本质是因为外磁场导致的共轭环上的环电流产生的感生磁场会一定程度对外磁场有抵消（屏蔽）作用。对于反芳香性体系，环电流的感生磁场和外磁场方向相同，故会加强外磁场，所以这个屏蔽值为负，NICS 即为正。

最初的 NICS 定义，也被称为 NICS(0)，随着 NICS(0) 被用于各种不同的体系，其已经被指出在不少体系都不适用。目前推荐使用 NICS(1)_ZZ 来衡量芳香性，NICS(1)_ZZ 就是指的在环平面上方 1 埃处（具体来说，是以环中心为起点向垂直于环平面方向挪1埃），垂直于环平面方向的屏蔽张量分量值的负值。由于在环上方 1 埃处 sigma 轨道效应影响很小，所以 NICS(1)_ZZ 衡量的主要是 pi 芳香性，而绝大多数结构的芳香性直接来源就是 pi 芳香性，因此这中衡量方法是非常有效的。

在NICS的基础上有人还做了进一步扩展，比如不仅仅考虑某个特定的点的屏蔽值，而是将磁屏蔽值作为一个实空间函数来考虑，通过绘制曲线图、平面图、等值面图来研究以获得更全面的信息，称为 ICSS 或 NICS scan。其中 ICSS 已经在笔者的第一篇博文中使用到了，这篇博文主要使用 NICS scan 来对该结构的芳香性进行衡量。

一维 NICS 曲线图

Multiwfn 能够很方便的完成一维的 NICS scan，在第一篇博文中，笔者已经将该分子的结构优化好了，并且确定是一个具有 D2h 对称性的平面结构。体系处于 Z=0 的 XY 平面上。此例对它绘制从环中心下方 10 埃到环中心上方 10 埃的NICS曲线，扫描路径垂直穿越环中心，并且考察的是垂直于环方向的 NICS 分量。启动 Multiwfn 输入优化好的带有笛卡尔坐标的结构文件，这里假设叫 NICS_opt.gjf，接着输入以下内容。

25   		// 离域化和芳香性分析
13   		// 一维NICS扫描和积分
2    		// 通过一批原子定义环，将扫描的两个端点置于环中心上方和下方一定距离处，且连线垂直穿越环中心
1-6   		// 用于1-6号碳原子定义环。之后屏幕上会显示通过最小二乘法拟合得到的环平面的法矢量
[Enter]  	// 用环上的原子的几何中心作为环中心。
10   		// 一个端点位于环下方10埃处
10   		// 另一个端点位于环上方10埃处。之后会从屏幕上看到扫描的起点和终点分别为0,0,-10和0,0,10
200   		// 扫描的点数。此例相当于每隔约0.1埃一个点，足够精细了。点数越多计算越耗时
1   		// 产生Gaussian输入文件
template_NMR.gjf 	//这是一个 Gaussian 做 NMR 计算的模板文件。打开它便知计算是在 B3LYP/6-31+G* 下进行的，和普通输入文件唯一的差异是坐标部分用 [geometry] 代替

现在 Multiwfn 在当前目录下产生了 NICS_1D.gjf 感兴趣的同学可以用 GaussView 打开这个文件，这里就不展示了。这里展示一下 template_NMR.gjf 的内容。

#p b3lyp/6-31+G* NMR

template file

  0  1
[geometry]

使用 Gaussian 对 NICS_1D.gjf 进行计算，该任务计算起来很快，比优化快得多，假设得到的输出文件为 NICS_1D_all.out。接着在 Multiwfn 中输入 2，然后输入 NICS_1D_all.out 的绝对路径。在之后的界面里可以选 1 和 2 分别绘制和保存 NICS 曲线图像，选 3 可以导出数据以便在 Origin 等程序里更灵活地绘制（导出的文本文件里每一列是什么含义在屏幕上提示得超级清楚，别无视）。这里我们选 3 得到数据，保存到名为 NICS_1D_all.txt 的文本文件中。

接着计算计算垂直于环的 NICS-sigma 和 NICS-pi 曲线，将上述的 template_NMR.gjf 分别改为 template_NMR_sigma.gjf 和 template_NMR_pi.gjf 文件。

#p b3lyp/6-31+G* nmr=csgt iop(10/93=2)

template pi file

  0  1
[geometry]

AICD.txt

22,28,31-33

其中 nmr=csgt iop(10/93=2) 必须写，AICD.txt 是产生的给 AICD 程序用的文件，当前情况不用管，文件名随意，产生后可以删掉。22,28,31-33 是波函数 pi 轨道的序号，可以通过 Multiwfn 轻松找到。同样如果要考察 sigma 轨道，只用将最后的轨道序号修改为 sigma 轨道序号就可以了。

#p b3lyp/6-31+G* nmr=csgt iop(10/93=2)

template file

  0  1
[geometry]

AICD.txt

1-21,23-27,29-30

重复上述操作，就能得到整体 NICS 曲线、sigma 贡献的 NICS 曲线以及 pi 贡献的 NICS 曲线的数据。再通过 Origin 或者其他绘图软件或程序绘图就可以了。笔者这里使用 Python 中的 matplotlib、proplot 库编写了一个绘制曲线的脚本，绘制的结果如下所示。

如图所示，X=0 的位置是环中心位置，在这个位置整体的 NICS 非常正，如果只考虑整体的 NICS 很容易将其归为反芳香性体系。但是将其分解为 sigma 贡献的和 pi 贡献的之后，可以发现 pi 贡献的 NICS 曲线都是负的，并且在离平面上下 1 埃的地方 NICS 最负。也就是说分子上方和下方大约 1 埃处对垂直于环平面方向的磁场的屏蔽最强。

二维 NICS 平面图

上一节介绍了如何使用 Multiwfn 计算并生成 NICS 曲线图，本节介绍如何使用 Multiwfn 计算并生成 NICS 平面填色图。计算的流程和绘制 NICS 曲线图的流程相差不大，启动 Multiwfn 后输入 NICS_opt.gjf，接着输入以下内容。

25   		// 离域性与芳香性分析
14   		// 绘制NICS二维平面图
1   		// 填色图
[Enter]  	// 用默认的格点数，即两个方向都是100个点，因此要计算100*100=10000个Bq
0   		// 设置延展距离
1			// 1 Bohr。
1   		// XY平面
1a   		// Z=1埃
1   		// 产生Gaussian的NICS二维扫描的输入文件
template_NMR.gjf   //这是Gaussian做NMR计算的模板文件，里面原子坐标部分用[geometry]代替，会被自动替换

这里的 template_NMR.gjf 和上一节的提到的是一样的。同样，Multiwfn 会在当前文件夹下输出一个名为 NICS_2D.gjf 文件，直接用 Gaussian 计算这个文件。假设得到的输出文件名为 NICS_2D.out，然后在 Multiwfn 界面里输入 2 选择载入 Gaussian 输出文件，然后输入其绝对路径。接着程序问你要获得哪种 NICS，可以选各项同性值、各向异性值、平行于笛卡尔 X 或 Y 或 Z 方向的分量值、顺着特定矢量的分量值。这里选择 5，即平行于 Z 方向的值。由于当前体系平行于 XY 平面，因此这么选得到的对应于一般意义的 NICS_ZZ。

关闭图像，输入以下命令修改作图效果：

4   		//显示原子标签
1  			//红色
8    		//显示化学键
14  		//棕色
17   		//设置显示标签的距离阈值
5   		//距离作图平面5 Bohr以内的原子标签才显示出来
y   		//更远的原子用细体字显示标签
1   		//修改色彩刻度范围
-45,45
-8   		//坐标轴改为以埃为单位
-2   		//修改坐标轴刻度
2,2,10
2   		//显示出等值线
3   		//修改等值线设置
8   		//按等差数列生成等值线数值
-50,5,21  	//起始值，步长，步数
y			//替换原有的等值线数值。之后等值线数值为-50,-45,-40...略...40,45,50
1   		//保存并返回
-1  		//重新作图

得到的图像如下图所示，此图颜色越深说明 NICS_ZZ 越负，即对垂直于体系平面方向的磁场屏蔽效应越强。此图中两个五元环的颜色比中间四元环的颜色要深体现出五元环确实是存在芳香性的。

绘制曲线的 Python 脚本

Python 版本为 3.9、matplotlib 版本为 3.4.3、proplot 版本为 0.9.5、pandas 版本随意。

import pandas as pd
import proplot as pplt
import matplotlib.pyplot as plt

from proplot import rc

# 读取 Multiwfn 生成的 NICS-1D-all、sigma、pi 的数据
nics_pi = pd.read_csv("NICS_1D_pi.txt", delimiter="\s+")
nics_sigma = pd.read_csv("NICS_1D_sigma.txt", delimiter="\s+")
nics_all = pd.read_csv("NICS_1D_all.txt", delimiter="\s+")

# 设置绘图的默认参数，如字体、字号等
rc['font.family'] = "Arial"
rc['tick.width'] = 1.3
rc['meta.width'] = 1.3
rc['font.size'] = 9.5
rc['label.size'] = 12.5
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

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

colors = ['off yellow', 'cherry red', 'true blue']

# 绘制整体贡献的 NICS 曲线
ax.plot(nics_all.iloc[:, -2], nics_all.iloc[:, -1], label='total', linewidth=1.2, color=colors[0])
# 绘制 sigma 体系贡献的 NICS 曲线
ax.plot(nics_sigma.iloc[:, -2], nics_sigma.iloc[:, -1], label='sigma', linewidth=1.2, color=colors[1])
# 绘制 pi 体系贡献的 NICS 曲线
ax.plot(nics_pi.iloc[:, -2], nics_pi.iloc[:, -1], label='pi', linewidth=1.2, color=colors[2])

# 在 y=0 处绘制一条虚线
axhline = ax.axhline(y=0, color='black', linewidth=1.2)
# 设置为最底层
axhline.set_zorder(0)

# 设置图例
ax.legend(loc='best', ncols=1, fontweight='bold', fontsize='12.5', frame=False, bbox_to_anchor=(0.95, 0.96))

# 格式化图像
fig.format(
    grid=False, ylabel='Shielding (in ppm)', xlabel='Position (in Å)',
    xlim=(-10, 10), xminorlocator=1, xlocator=2, ylim=(-30, 40), yminorlocator=5, ylocator=10
)

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