使用AICD 2.0绘制磁感应电流图

sobereva

使用AICD 2.0绘制磁感应电流图

Drawing magnetic induction current map using AICD 2.0

文/Sobereva @北京科音

First release: 2015-Jun-13  Last update: 2024-Jan-26

以前写过一篇帖子《使用AICD程序研究电子离域性和磁感应电流密度》（http://sobereva.com/147），那个帖子当时用的是1.5.7.1版。AICD本身用起来简单，但是使用前必须自行编译高斯，这使得绝大部分用户都无缘使用AICD。好在从Gaussian09 D.01开始直接提供了IOp(10/93)选项来输出AICD所需的电流密度张量数据，并且AICD 2.0版开始可以直接读取这些数据绘图，比原先实在方便太多了，这使得AICD能够广泛普及开来。这个帖子就说下AICD 2.0结合G09 D.01的使用方法。老版本的AICD的输出文件和各个选项在2.0版里完全没变，2.0也没有什么新增的选项，所以本帖仅仅用薁分子做个演示，更多内容请看上面的帖子。

1 安装AICD 2.0

注：由于原版的AICD 2.0与Gaussian 16兼容性有问题，而且与比较新的gcc编译器的兼容也有问题（像CentOS 6.x这种老系统带的gcc没问题，而CentOS 7.x的就不行了），因此这里给的是http://bbs.keinsci.com/thread-13577-1-1.html里提供的修改版AICD包，安装方法也是针对这个修改版而言的。

在这里下载修改版AICD 2.0：http://sobereva.com/soft/AICD-2.0.0-modified.tar.gz。

将这个AICD 2.0包解压到任意一处，比如/sob/AICD-2.0.0，然后进入此目录，输入make。然后在用户目录的.bashrc文件里添加alias AICD=/sob/AICD-2.0.0/AICD。重新登录终端或输入bash使之生效。现在就可以在任意目录下直接通过AICD命令调用AICD了。

如果运行时提示无权限的问题，运行chmod +x /sob/AICD-2.0.0/* -R给此目录下所有文件加可执行权限。


2 计算薁

输入文件azulene.gjf内容如下。NMR=CSGT必须写，IOp(10/93=1)代表输出电流密度张量到坐标后面的那个文件里，此例是test.txt。

1. %chk=./azulene.chk
   
2. # b3lyp/6-31g(d) nmr=csgt iop(10/93=1)
   
3. 
   
4. b3lyp/6-31g(d) opted
   
5. 
   
6. 0 1
   
7. C 0.00000000 0.00000000 2.50418951
   
8. C 0.00000000 1.26637425 1.91156618
   
9. C 0.00000000 1.59417912 0.55241750
   
10. C 0.00000000 -1.26637425 1.91156618
    
11. C 0.00000000 0.75015229 -0.55452268
    
12. C 0.00000000 -1.59417912 0.55241750
    
13. C 0.00000000 -0.75015229 -0.55452268
    
14. H 0.00000000 0.00000000 3.59316833
    
15. H 0.00000000 2.10900376 2.59929497
    
16. H 0.00000000 2.66005032 0.32418052
    
17. H 0.00000000 -2.10900376 2.59929497
    
18. H 0.00000000 -2.66005032 0.32418052
    
19. C 0.00000000 1.15011348 -1.90182297
    
20. H 0.00000000 2.17701495 -2.24803140
    
21. C 0.00000000 0.00000000 -2.70850570
    
22. H 0.00000000 0.00000000 -3.79360290
    
23. C 0.00000000 -1.15011348 -1.90182297
    
24. H 0.00000000 -2.17701495 -2.24803140
    
25. 
    
26. test.txt
    
27. 
    

复制代码

用g09 < azulene.gjf > azulene.out运行此文件，在当前目录下得到azulene.out和test.txt。


3 产生感应电流图

我们把外磁场垂直于环平面加。若打开azulene.out会看到垂直于分子平面的是x轴，即外磁场矢量应为1 0 0。我们输入
AICD -m 4 -b 1 0 0 -pov azulene.out
此时AICD就会根据azulene.out和test.txt里的信息生成AICD格点数据，并且产生文件名以azulene_40000_0.050_1_0_0_Aniso_4.4为开头的5个文件。将这5个文件拷到Windows系统下，安装好Povray渲染器，然后打开RenderMich.pov后缀的那个，渲染，即得到下图，包含多个视角：

可见是有明显环电流产生的。


4 只获得特定轨道的贡献

也可以只获得特定轨道对AICD的贡献，写上IOp(10/93=2)，并且在输入文件末尾写上要考虑的轨道编号即可。

比如薁这个分子，我们想得到所有pi轨道对AICD的贡献。首先我们先找出pi轨道编号，虽然可以挨个看轨道图来找，但对于轨道比较多的体系会比较劳神。对于薁这样纯平面的体系找pi轨道的最方便的方法是用Multiwfn (http://sobereva.com/multiwfn)。将薁这个体系的fch文件载入Multiwfn，依次输入100、22、2（因为当前体系平行于YZ平面）。程序马上就给出了所有pi MO的编号，我们把占据数为2.0的那5个复制到前面azulene.gjf的末尾去（从命令行窗口直接复制数据的操作见Multiwfn手册5.4节）。此时azulene.gjf内容如下

1. %chk=./azulene.chk
   
2. # b3lyp/6-31g(d) nmr=csgt IOp(10/93=2)
   
3. 
   
4. b3lyp/6-31g(d) opted
   
5. 
   
6. 0 1
   
7. [坐标，略]
   
8. 
   
9. test.txt
   
10. 
    
11. 26
    
12. 31
    
13. 32
    
14. 33
    
15. 34
    
16. 
    

复制代码

为了简便，上面也可以写成26 31-34。用前述方法绘图，结果如下

十分值得一提的是，用上面的做法不仅可以得到不同轨道对AICD的贡献，也可以得到不同轨道对磁屏蔽值、磁化率的贡献（即它们可以严格分解成轨道的贡献）。比如某体系有13个轨道，用上面的关键词，分别选择1-5号轨道，以及6-13号轨道，二者得到的磁屏蔽值或者磁化率的加和和直接算整体时的结果是一样的。因此，可以利用这个方法考察NICS-sigma和NICS-pi。关于这点，我专门写了篇文章介绍：《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（http://sobereva.com/670）。

对于非限制性开壳层（U）形式的计算，IOp手册里没写清楚怎么指定要考虑的轨道。我这里根据我读源代码以及实际测试做出的理解介绍一下。比如计算三重态水分子，原本一共有6个占据的alpha轨道和4个占据的beta轨道，在你设置要考虑的轨道时，记住beta占据轨道的序号是排在alpha占据轨道后面的，而不用管alpha和beta的空轨道。比如，如果你要让原先的那10个占据轨道都被考虑，输入文件末尾就写1-10。如果你只想考虑alpha占据轨道当中后三个、beta占据轨道中第2个，就应该写4-6 8。之所以是8，是因为6+2=8，即原本的alpha占据轨道数目加上要考虑的这个beta的序号。空轨道不用管，比如你把11、12等序号超过原先占据轨道总数的轨道也纳入的话，不会影响结果。

对于非平面体系，没法严格定义哪些是sigma轨道、哪些是pi轨道，或者说判断pi轨道的方式有很大含糊性、人为任意性。但是也没关系，强大的Multiwfn可以给出各个轨道的pi成份，超过特定阈值（比如80%）的轨道你都可以姑且当做pi轨道。关于怎么用Multiwfn计算pi成份，看《Multiwfn已支持计算任意轨道的pi成份》（http://bbs.keinsci.com/thread-13110-1-1.html）。



5 其它

如果对本文有任何不清楚的地方，或者想弄明白AICD的各种选项、细节以及相关理论，请看《使用AICD程序研究电子离域性和磁感应电流密度》（http://sobereva.com/147）。

笔者通过此文的方法考察了电子结构非常特殊的18碳环体系，得到了漂亮的图像和很有意义的结果，充分揭示了此体系由in-plane和out-of-plane电子产生的双芳香性特征。文章介绍见《一篇最全面、系统的研究新颖独特的18碳环的理论文章》（http://sobereva.com/524），其中与成键和芳香性有关的研究部分后来发表在了Carbon, 165, 468 (2020)，非常建议大家阅读和引用，这是AICD很好的应用实例。下图左边是18碳环的in-plane pi电子的环电流，右边是out-of-plane pi电子的环电流。

在《深入揭示18碳环的重要衍生物C18-(CO)n的电子结构和光学特性》（http://sobereva.com/640）中介绍的笔者研究C18-(CO)n的电子离域和芳香性特征的文章Chem. Eur. J., 28, e202103815 (2022)中，以及在《不寻常的环[18]碳前驱体C18Br6的电子结构和芳香性》（http://sobereva.com/664）中介绍的笔者研究C18-Br6的电子结构的文章Chem. Eur. J., e202300348 (2023) DOI: 10.1002/chem.202300348中，以及在《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（http://sobereva.com/696）介绍的笔者研究B6N6C6芳香性并与B9N9对比的Inorg. Chem., 62, 19986 (2023)文章中，以及在《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（http://sobereva.com/725）介绍的笔者研究18氮环的芳香性的文章中，都使用了AICD来直观展现体系的感生电流，是AICD的很好的应用范例，建议阅读里面的讨论，非常推荐作为AICD的应用范例引用。

有另外两个也非常知名、常用的考察磁感生电流的程序GIMIC和SYSMOIC，与AICD有很强的互补性，强烈建议读者仔细看《考察分子磁感生电流的程序GIMIC 2.0的使用》（http://sobereva.com/491）和《使用SYSMOIC程序绘制磁感生电流图和计算键电流强度》（http://sobereva.com/702）。

有人问我除了考察感生电流外还有没有其它可视化研究芳香性的方法，实际上在Multiwfn中支持很多，如ICSS、LOL-pi、AdNDP，分别参看
通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性
http://sobereva.com/216
使用Multiwfn巨方便地绘制二维NICS平面图考察芳香性
http://sobereva.com/682（http://bbs.keinsci.com/thread-39020-1-1.html）
在Multiwfn中单独考察pi电子结构特征
http://sobereva.com/432（http://bbs.keinsci.com/thread-10610-1-1.html）
使用AdNDP方法以及ELF/LOL、多中心键级研究多中心键
http://sobereva.com/138
Multiwfn还支持非常丰富的定量考察芳香性的方法，参见《衡量芳香性的方法以及在Multiwfn中的计算》（http://sobereva.com/176）。

aqhuangry

谢谢版主分享。能否提供一下AICD 2.0的下载连接？谢谢！

sobereva

aqhuangry 发表于 2015-6-14 15:28
谢谢版主分享。能否提供一下AICD 2.0的下载连接？谢谢！

补上了，见第1节

aqhuangry

好的。十分感谢！

moonking

多谢分享。最近尝试自己分析下自己的体系，前面所有运行都没有问题，只有最后Povray渲染时，双击RenderMich.pov后缀的那个，然后点击Run，可是为什么出来的是全黑的图，没有环电流……求指点

sobereva

moonking 发表于 2016-2-19 05:26
多谢分享。最近尝试自己分析下自己的体系，前面所有运行都没有问题，只有最后Povray渲染时，双击RenderMich ...

绘制文中的例子能得到正常结果不？

moonking

sobereva 发表于 2016-2-19 14:54
绘制文中的例子能得到正常结果不？

谢谢您的建议，折腾了好几天，还是搞不懂，只能求助了……
我严格重复文中例子，可以顺利得到azulene.out和test.txt，运行AICD也没有提示错误，可以生成5个对应的pov文件，可是RenderMich.pov的大小仍然是0kb。。。双击打开，点run后还是一片黑，想上传附件来着，可是提示不能上传0KB文件，可见确实是这一步出了问题，求指导，谢谢了……

我怀疑问题出在如下步骤，运行AICD -m 4 -b 1 0 0 -pov azulene.out时最后要输出pov文件部分时提示如下

****************************************************************************************
modifying azulene.pov to azulene_40000_0.050_1_0_0_Aniso_4.4.Molekuel.pov
/home/wjn/AICD-2.0.0/AICD: 行 656:  5676 Aborted                 (核心已转储)$An                                                                                                                imation $inputbasename.pdb $Atom1 $Atom2 $Atom3 < "$RenderMich" > $Pov_RenderMic                                                                                                                h_pov
Replacing Camera Postion to \<0, 0, -250\>
Inserting #version directive in RenderMich.pov
calling Replace_Pov_Filenames azulene_40000_0.050_1_0_0_Aniso_4.4.Isoober.inc
calling Replace_Pov_Filenames azulene_40000_0.050_1_0_0_Aniso_4.4.Molekuel.inc
calling Replace_Pov_Filenames azulene_40000_0.050_1_0_0_Aniso_4.4.Molekuel.pov
calling Replace_Pov_Filenames azulene_40000_0.050_1_0_0_Aniso_4.4.RenderMich.pov
calling Replace_Pov_Filenames azulene_40000_0.050_1_0_0_Aniso_4.4.Rotate.inc
calling AICD-isocut -m -1 -M -1 -r 10000000 < azulene_40000_0.050_1_0_0_Aniso_4.                                                                                                                4.Isoober.inc.noncut > azulene_40000_0.050_1_0_0_Aniso_4.4.Isoober.inc
This is AICD-isocut....
  arrow_maxlength = -1
  arrow_minlength = -1
  Point of origin:  0 0 0
  Maximum distance: 1e+07
Number of deleted arrows: 0


Call PovRay manually with the following commandline:

povray +FN +iazulene_40000_0.050_1_0_0_Aniso_4.4.RenderMich.pov +oazulene_40000_                                                                                                                0.050_1_0_0_Aniso_4.4.png /home/wjn/AICD-2.0.0/povray-AICD-templates/AICD-pov.in                                                                                                                i -P -W1024 -H768

sobereva

moonking 发表于 2016-2-24 08:50
谢谢您的建议，折腾了好几天，还是搞不懂，只能求助了……
我严格重复文中例子，可以顺利得到azulene.ou ...

/home/wjn/AICD-2.0.0/AICD: 行 656:  5676 Aborted                 (核心已转储)$An
这句已经提示可能出现错误了

什么原因我也不好说

moonking

换了台Linux的电脑重新来了一遍，有了一个更详细的提示

“Error in `AICD-rotate_mol': free(): invalid next size (fast): 0x099340d8 ”

经过不同参数测试输出，只有要求 -pov 时才会报错，所以是我的AICD-rotate编译出问题了？

因为是做合成出身的，也不太精通计算机，发的文章中包含的一点点理论分析还都是用Multiwfn做的，所以对这个问题没辙了……

我给作者写了封邮件咨询，看看他们怎么回复。

如果可能，可否把你编译好的文件也打包传一份，我试试看？

多谢啦

sobereva

moonking 发表于 2016-2-27 01:05
换了台Linux的电脑重新来了一遍，有了一个更详细的提示

“Error in `AICD-rotate_mol': free(): invalid ...

这是我编译的

moonking

sobereva 发表于 2016-2-27 02:58
这是我编译的

今天终于搞定了，多谢啦，我发现原来我用的是64位的Linux系统，换成32位的之后马上就没问题了……

MAKE

moonking 发表于 2016-3-4 13:01
今天终于搞定了，多谢啦，我发现原来我用的是64位的Linux系统，换成32位的之后马上就没问题了……

根据你之前的那个报错应该是gcc 和 g++的版本导致，换用4.4 或4.5就应该能解决，不需要换32位系统

sunchunlin

我也发现这个AICD软件和用的系统有关系，好几台系统用的centos7的就是会出现，渲染文件大小为0的问题，而同样的操作换了一台红帽子系统，就正常了。编译器也没有报错。而且都是最新的编译器。

ggdh

有人问显示不全怎么办，或者是要换方向看怎么办，可以直接用文本编辑工具打开RenderMich.pov文件改文件开头下面的参数（写完后才发现sob的原帖中有更多选项的说明，可以直接过去看。。。）
camera{
  location <-1000, 0, 0>   //调整摄像头的位置，这里放在x轴的-1000的位置，如果显示不全可以放的更远
  direction < 2, 0, 0 >      //这个是摄像头看的方向，这里看向x的正方向，也就是看向yz平面
  sky < 1,0,0>
  up < 0, 1, 0 >              //色相头上下的广角，显示不全也可以把这个数调大。这里从x方向看yz平面，0,1,0表示上下是y轴
  right < 0, 0, 1.33 >      //色相头左右的广角，上下和左右的广角最好按比例调，否者图像会变形
  orthographic
}

light_source { < -160, 60, -60 > color rgb < 0.35, 0.35, 0.35 > }
light_source { < -160, 60, 60 > color rgb < 0.35, 0.35, 0.35 > }
light_source { < -160, -60, -60 > color rgb < 0.35, 0.35, 0.35 > }
light_source { < -160, -60, 60 > color rgb < 0.35, 0.35, 0.35 > }//这4个是加了4个光源，左边矢量控制光源的位置，右边控制光源的颜色，如果换方向看就要调整光源的位置。注意要把光源放在摄像头的同侧，不然是逆光。。

ggdh

今天有人问：对于开壳层体系，发现β轨道无法正常运算，
研究一番之后发现，β轨道的编号并不能用正常的编号
正常的β轨道的编号是α占据+α空+β编号。
而这里β轨道编号应该是α占据+β编号。
举个例子：
甲醛阳离子自由基分子，α轨道34个β轨道34个，αHOMO 是α轨道中的第8号轨道，βHOMO是β轨道中的第7号。
Multiwfn会显示βHOMO是34+7=41号轨道，而这这里输入βHOMO应该输入8+7=15号轨道。