Post

Study notes：芳香性的计算

Posted Updated
By Bane Dysta
9 min read
Study notes：芳香性的计算

根据sob老师的博文衡量芳香性的方法以及在Multiwfn中的计算，计算芳香性的方法有很多种。笔者最近因为课题需要把大大小小的方法基本都实测了一遍，在此记录一下几个笔者发现的比较靠谱的方法。文中计算均以芴正离子为例。

一维NICS曲线

按照使用Multiwfn绘制一维NICS曲线并通过积分衡量芳香性操作：

1
2
3
4
5
6
7
8
9
10

25   //离域化和芳香性分析
13   //一维NICS扫描和积分
2    //通过一批原子定义环，将扫描的两个端点置于环中心上方和下方一定距离处，且连线垂直穿越环中心
5-9   //先看五元环部分。用芴环的5-9号碳原子定义平面。之后屏幕上会显示通过最小二乘法拟合得到的环平面的法矢量
[回车]  //用环上的原子的几何中心作为环中心。此处也可以自己输入其它方式得到的中心坐标，比如按照http://sobereva.com/108用Multiwfn做电子密度拓扑分析得到的环临界点（rcp）
[回车]   //一个端点位于环下方10埃处
[回车]   //另一个端点位于环上方10埃处。之后会从屏幕上看到扫描的起点和终点分别为0,0,-10和0,0,10
[回车]   //扫描的点数200。此例相当于每隔约0.1埃一个点，足够精细了。点数越多计算越耗时
1   //产生Gaussian输入文件
D:\program\Multiwfn_3.8_dev_bin_Win64\examples\NICS_scan\template_NMR.gjf   //载入预设模板

产生NICS_1D.gjf，提交该计算。(记得改电荷！)算完之后，载入文件绘图：

1
2
3

2
NICS_1D.log
1  //绘图

可以看到五元环部分是有反芳香性的。NICS曲线积分值285.89，这个挺高了。

NICS_ZZ (NICS-2D scan)

按照使用Multiwfn巨方便地绘制二维NICS平面图考察芳香性操作：

1
2
3
4
5
6
7
8
9
10

25   //离域性与芳香性分析
14   //绘制NICS二维平面图
1   //填色图
[回车]  //用默认的格点数100*100
8   //绘制某些原子构成的拟合平面的上方或下方的平面
5-9   //环中的原子。这一步输出的向量1.00000000   -0.00000000   -0.00000000记住，有用。
1   //绘制的是与拟合平面相平行而在它上方1埃的平面。如果输入负值代表在下方
15   //绘图平面的边长为15埃
1  //输出gjf文件
D:\program\Multiwfn_3.8_dev_bin_Win64\examples\NICS_scan\template_NMR.gjf //载入预设模板

产生NICS_2D.gjf，提交该计算。(记得改电荷！)算完之后，载入文件绘图：

1
2
3
4

2
NICS_2D.log
0  //取特定矢量投影
1.00000000   -0.00000000   -0.00000000 //刚刚打印的矢量

随后调整图像。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

4 //显示绿色标签
2
8 //显示棕色化学键
14
-8 //改单位
1 //改色彩范围
-50,50
-2 //改坐标轴刻度
2,2,10
19 //改配色
8
17 //改显示标签阈值
3
2 //等值线
3
8
-20,15,6 //根据载入时打印的minimum和maximum确定
y
1

最终效果：

看到三个环在S0其实都有一定反芳香性。

再来看看S1的。流程与上述完全一致，只是生成的NICS_2D.gjf要手动补上guess=read，读取ΔSCF波函数：

发现激发态下芳香性发生了反转，原本的反芳香环变成了芳香性的，这一点从S0激发能从1.4558eV减到0.8734eV上也可以侧面印证。

ICSS

载入文件和，输入：

1
2
3
4
5
6

25  //离域性和芳香性分析
3   //ICSS分析
-10   //设定格点数据计算区域的延展距离
12    //延展距离设为12 Bohr，也就是相对于边界原子向四周扩12 Bohr。延展距离不能太小，否则绘制出来的ICSS等值面会被截断。通常12 Bohr够大了
1     //低质量格点。Gaussian接下来将总共计算约125000个位置的磁屏蔽张量。如果用更高质量的格点计算耗时将会很高
n    //不跳过Gaussian输入文件生成的步骤

此时程序请求一个gjf文件，这个gjf文件只要与载入的那个一致即可，别的就按照默认的来：

1
2
3
4
5
6
7
8
9
10
11

# b3lyp/6-31+g(d)

title

1 1
 C                  0.00000000    2.51309068    1.23211960
 C                  0.00000000    3.45934787    0.20266322
 C                  0.00000000    3.04461175   -1.12685120
 C                  0.00000000    1.67397844   -1.49236328
...

注意：如果要并行计算，这里不要有chk行，而且生成输入文件后要手动移除所有guess=read。 计算完成后，下载所有log，然后：

1
2
3
4
5

2  //载入log文件
path\to\NICS  //载入path\to路径下以NICS起头的NICS0001.log、NICS0002.log...
5  //查看ICSS-ZZ，或者1查看ICSS
1  //查看等值面
2  //导出格点数据

导出的格点数据用VMD作图，可以看到ICSS等值面全貌。

有了ICSS，NICS-2D和NICS曲线都可以从这里载入，而不必重新交一遍相应任务，因为ICSS相当于把格点范围内每个点的磁属性都扫出来了。

1
2
3
4
5
6
7
8
9

1000
2
-1  // 使自定义函数通过格点数据插值获得数值
4
100
1  // 填色图
[enter]
3  // YZ
1a // 作图距离环上1埃

这样就弹出来了一个YZ平面上方1埃的等值面图，可以根据2D绘图流程修改作图显示风格。

同样的，一维曲线也可以：

1
2
3
4
5
6
7
8
9
10
11
12
13

1000
2
-1  // 使自定义函数通过格点数据插值获得数值
100
21
5-9  // 计算五元环几何中心，找到Geometry center为0.00000000    0.00000000    0.49697297
1-6  // 再计算一下六元环几何中心，0.00000000    2.09835945   -0.13196727 
q
0
3   //绘制函数在一条直线上的变化曲线
100   //自定义函数
2    //自定义这条线的两端坐标
0,0,0.49697297,0,2.09835945,-0.13196727 A    //始端为五元环中心(0,0,0.49697297），末端为六元环中心(0,2.09835945,-0.13196727)，单位埃

作图如下：

因为曲线不太平滑，看了一眼，我算的这个芴正离子是在XZ平面上的，不同于卢老师给的例子在XY平面，当时应该导出ICSS-YY，对应垂直于环平面的投影。

AICD

靠谱分析方法的全都是磁属性😸

note
note Aromaticity
This post is licensed under CC BY 4.0 by the author.
Total hits!