Study notes：芳香性的计算

Posted Apr 6, 2025 Updated Apr 19, 2025

By Bane Dysta

15 min read

Study notes：芳香性的计算

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

一维NICS曲线

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，提交该计算。(记得改电荷！)算完之后，载入文件绘图：

2
NICS_1D.log
1  //绘图

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

NICS_ZZ (NICS-2D scan)

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

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

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

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

随后调整图像。

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

载入文件和，输入：

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

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

# 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相当于把格点范围内每个点的磁属性都扫出来了。

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

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

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

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

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

安装

参考blog31

Gaussian计算

若不需要单独考虑某些轨道的贡献，直接使用NMR=csgt和iop(10/93=1)即可

# nmr=csgt b3lyp/6-31g(d) iop(10/93=1)

 title

 [charge] [spin]
 [geometry]

 filename.txt

若要单独考虑某些轨道的贡献：

# nmr=csgt b3lyp/6-31g(d) iop(10/93=2)

 title

 [charge] [spin]
 [geometry]

 filename.txt

 [orbital1]
 [orbital2]
 ...

对于开壳层，beta列轨道号=对应alpha列轨道号+总占据轨道数。

AICD计算

下载filename.txt和高斯输出文件，在对应目录执行AICD：

AICD --povrayinput -m 4 benzene_aicd.log --runpovray

我这边看着有bug，运行结束后应该产生.png图片的，但实际上生成的是.4文件，似乎是由-m 4引起的。这边把.4文件重命名为png就可以观看渲染结果了。

参数：

--povrayinput：生成povray的输入文件。
--runpovray：并执行渲染。
-c; --povraycopy：保留povray的输入文件(指从tmp复制一份到当前目录)。
-f; --force：覆盖当前已存在的povray的输入文件。
-m [number];--view [number]：指定渲染什么视图。1=主视图；2=主视图+表面箭头；3=三视图；4=三视图+表面箭头。默认值似乎是3。
-b [x] [y] [z]; --magneticfield [x] [y] [z]：指定外加磁场方向矢量，默认为0 0 1。
-p [number]; --gridpoints [number]：设置数据点数量，默认为40000，n值越大，生成的箭头和等值面越平滑，但相应计算成本会增加。
-s; --smoothisosurface：平滑化，避免边缘棱角。相应计算成本会增加。
-maxarrowlength [number]：设置最大箭头长度，超过[number]的箭头不渲染，用于屏蔽一些经常存在的飘逸的箭头。
-l [number]; --isolimit [number]; --isovalue [number]：指定等值面值，默认0.05。通常不用调，调太小有可能引起等值面被截断。
-h：列出所有参数。

PovRay渲染

默认渲染参数很可能出现显示不全、太挤、很糊等问题，想做出高质量的图像一般都需要手动调整参数。首先渲染一个高格点质量三视图看看情况，决定怎么调参数：

AICD -p 1000000 --runpovray --povrayinput -c -f -b 0 0 1 -s -m 4 -maxarrow 1 benzene_aicd.log
for file in *.4; do mv "$file" "$file.png"; done

进入生成的.d结尾的目录，找到RenderMich.pov，这个就是渲染参数文件。

相机高度

如果渲染图像不全，修改这部分的location

camera{
  location <  0, 0, -160 >
  direction < 0, 0, 2 >
  sky < 0, 1, 0 >
  up < 0, 1, 0 >
  right < 1.333, 0, 0 >
  orthographic
}

location表示相机的位置，可以推知z越负看到的分子越全。

旋转

如果渲染的三视图看着比较模糊，可以考虑分别渲染高质量的单个视角。使用-m 2参数，然后找到RenderMich.pov文件结尾的object部分：

{ object { Molekuel }
  object 
  { Isooberflaeche
    texture
    { pigment { color Yellow }
      finish { Plastic }
    }
  }
  object { Pfeile }

  #include "Rotate.inc"
  
  rotate <0,0,90>
  translate < 0,0,0 >
  scale 6
}

修改rotate <0,0,90>的值来更改旋转方向，分别渲染三视图。

对象间距

如果清晰度还行，但是三视图叠在一起了，找到RenderMich.pov文件结尾的三视图object部分，修改translate

object 
{ MolUndMag
  translate < -45,32,0 >
}
object 
{ MolUndMag
  rotate < 90,0,0 >
  translate < -45,0,0 >
}
object 
{ MolUndMag
  rotate < 0,90,0 >
  translate < -45,-32,0 >
}


object 
{ MolUndIso
  translate < -10,29,0 >
}
object 
{ MolUndIso
  rotate < 90,0,0 >
  translate < -10,0,0 >
}
object 
{ MolUndIso
  rotate < 0,90,0 >
  translate < -10,-29,0 >
}

object 
{ MolUndIso
  rotate < 0,180,0 >
  translate < 35,29,0 >
}
object 
{ MolUndIso
  rotate < 90,0,0 >
  rotate < 0,180,0 >
  translate < 35,0,0 >
}
object 
{ MolUndIso
  rotate < 0,90,0 >
  rotate < 0,180,0 >
  translate < 35,-29,0 >
}

一共三块object，每个块控制一个视图。translate后面的坐标控制每个对象的位置，改大一些可以增加间距。(Z坐标应当保持为0)

展示一个最终效果：

外环顺时针流动为芳香性；逆时针流动为反芳香性。

note

note Aromaticity

This post is licensed under CC BY 4.0 by the author.