考察分子磁感生電流的程序GIMIC 2.0的使用（含24分鐘演示視頻）

考察分子磁感生電流的程序GIMIC 2.0的使用

First release: 2019-Jun-22 Last update: 2020-Sep-18

1 前言

GIMIC的主要用途是基于量子化學程序的輸出信息來計算分子的磁感生電流，這對于考察芳香性尤為重要，見《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）中的相關介紹。GIMIC可以給出電流矢量、電流矢量的模，還可以對穿越特定截面的電流密度進行積分來定量考察電流的強度。

之前筆者在《使用GIMIC計算和分析磁感應電流密度》（http://www.shanxitv.org/155）中已經非常詳細地介紹了如何使用GIMIC考察磁感生電流密度，對相關原理也予以了說明。那篇文章里介紹的GIMIC是老版本，一個關鍵缺點是沒法結合最主流的量子化學程序Gaussian來使用。如今GIMIC已經有了2.0版本，已可以結合Gaussian使用了，用法也有很大變化。本文的目的就是介紹一下怎么將GIMIC 2.0和Gaussian結合做磁感生電流密度的分析。凡是前面那篇文章里已經講過的內容本文就不再怎么提了，所以讀完本文后別忘了把前面那篇文章也看一看。

GIMIC 2.0程序的在線手冊是https://gimic.readthedocs.io，里面有用法說明。但寫得很混亂，好多地方語焉不詳令人糊涂，有不少地方還是錯的。讀者讀過本文后基本就沒必要看那個手冊了。

筆者將此文的方法應用于了電子結構十分新穎的18碳環體系，得到了很有意義的結果，見此文介紹的筆者的ChemRxiv上的論文：《一篇最全面、系統的研究新穎獨特的18碳環的理論文章》（http://www.shanxitv.org/524），后來相關部分的研究筆者發表在了Carbon, 165, 468 (2020)中。此論文是本文介紹的方法的很好的應用實例，十分建議一讀，也非常歡迎引用。論文的補充材料中還給出了18碳環在外磁場下環電流的動態動畫，效果十分非常酷炫，強烈建議一看。

2 GIMIC 2.0的獲取和安裝

去GIMIC 2.0的官網https://github.com/qmcurrents/gimic，點擊Clone or download按鈕，點擊Download ZIP，即可把源代碼包下載下來。筆者在2019年6月18日下載的程序包可以在這里下：http://www.shanxitv.org/attach/491/gimic_2019-Jun-18.zip，本文所述都是對這個版本而言。

筆者使用的是CentOS 7.6系統，root賬戶，安裝過程與《在VMware 15中安裝CentOS 7.6的完整過程視頻演示》（http://www.shanxitv.org/454）里演示的一致。如果你的情況和本文的不同，本文的做法可能無法安裝成功，請根據提示隨機應變。

在系統里依次運行以下命令，安裝EPEL源，并且裝上cmake和conda。GIMIC 2.0是用Fortran和Python混合寫的，所以需要安裝conda來管理Python的包。
yum install epel-release
yum install cmake
yum install conda

然后運行下面的命令創建一個名為myconda的環境，并且初始化conda
conda create --name myconda
conda init bash
在~/.bashrc文件的末尾加入一行conda activate myconda使得每次進入終端后自動激活myconda環境。

退出終端，然后重新進入終端，輸入以下命令安裝GIMIC 2.0運行時依賴的一些Python庫
conda install cython numpy pyparsing

將GIMIC壓縮包解壓，比如筆者解壓到了/sob/gimic下面，進入此目錄后運行
./setup --omp
cd build
make install

編譯完成后，在~/.bashrc里末尾加入
export PATH=$PATH:/sob/gimic/build/bin

之后重新進入終端，程序就可以用了。

注：在上面的步驟中運行了conda init bash之后，會自動在~/.bashrc里添加一串內容，從# >>> conda initialize >>>開始到# <<< conda initialize <<<結束，這可能導致進入終端后命令提示符需要按一下回車之后才能顯示出來。只要把這段內容從.bashrc里刪掉就沒這個問題了，不影響gimic的正常運行。

3 使用方法

運行以下命令即可啟動GIMIC 2.0并讀取test.inp輸入文件里的設置進行運算，信息輸出到gimic.out里也同時顯示在屏幕上。
gimic test.inp |tee gimic.out
如果當前目錄下有名為gimic.inp的文件，直接輸入gimic就可以計算，程序會將它當做默認的輸入文件。

輸入文件里要指定XDENS和MOL文件的路徑，前者包含分子坐標和基組信息，后者包含密度矩陣和密度矩陣對磁場的導數。對于不同的量子化學程序，XDENS和MOL文件有不同的得到方法。下面是對于Gaussian的用法（筆者只測試了Gaussian09的情況）。

運行以下任務
%chk=test.chk
# B3LYP/6-311G* Int=NoBasisTransform IOp(10/33=2) NMR
[空行]
test
[空行]
0 1
[坐標...]

然后將chk轉換為fchk，把gimic\tools\g092gimic目錄里的BasisSet.py和Gaussian2gimic.py拷到當前目錄下，之后運行./Gaussian2gimic.py --input=test.fchk，當前目錄下就出現了XDENS和MOL。

如果是開殼層體系，關鍵詞應當用gfprint pop=full Int=NoBasisTransform IOp(10/33=2) NMR，令輸出文件后綴設為log，然后用Gaussian2gimic.py --input=test.log。

注意Gaussian計算時會自動把結構擺到標準朝向下，可能導致最終分子不是處于你預期的朝向。加上nosymm關鍵詞即可避免這點，詳見《談談Gaussian中的對稱性與nosymm關鍵詞的使用》（http://www.shanxitv.org/297）。

4 關于ParaView

如果你想把GIMIC 2.0得到的感生電流圖繪制出來，一般需要用ParaView可視化程序。這是一個非常強大、知名的體數據可視化程序。此程序可以在https://www.paraview.org/download/免費下載。筆者使用的是Windows 5.6.1版，大家就下載比如ParaView-5.6.1-Windows-msvc2015-64bit.exe就行了。

安裝之后需要設置一下。啟動ParaView，進入Edit - Setting，點擊右上角齒輪圖標，然后選上Auto Apply和Enable Auto MPI，Auto MPILimit設成實際物理核數。然后在Camera標簽頁里把Camera 2D Manipulators和Camera 3D Manipulators下面把左鍵、中鍵、右鍵分別設為Rotate、Pan、Roll，這樣視角操作就常用的GaussView一樣了，筆者感覺比較順手（在ParaView界面里，可以在圖形窗口上方點擊3D/2D按鈕來在透視視圖和正交視圖間切換）。

在主菜單的Tools - Manage Plugins里分別展開SurfaceLIC、StreamLinesRepresentation、StreamingParticles條目，將其中的Auto Load開啟。重啟ParaView后這些顯示方式將能夠使用。

5 實例：圖形化考察苯分子的磁感生電流

本節和下一節涉及的文件都可以在此處下載：http://www.shanxitv.org/attach/491/file.rar。

文件包里有benzene.gjf，使用Gaussian計算之，用formchk將chk轉換為fchk，再用./Gaussian2gimic.py --input=benzene.fchk命令處理，就得到了文件包中的MOL和XDENS。

5.0 產生ParaView繪圖需要的文件

在用ParaView繪圖之前，首先需要產生vti文件，里面可以記錄標量場和向量場，可以被ParaView直接可視化。

寫一個GIMIC的輸入文件gimic.inp，內容如下。此文件在本文的文件包里的3D子目錄下可以找到。//后面是筆者寫的說明。文件里面#后頭的內容都是注釋，會被程序忽略

calc=cdens //任務類型。cdens代表產生格點數據，integral代表計算截面積分值
basis="../MOL" //MOL文件的路徑。../代表上一級目錄
xdens="../XDENS" //XDENS文件的路徑
openshell=false //如果是開殼層體系，需要設成true
magnet=[0,0,1] //外磁場的方向，這里要求外磁場方向是Z軸正方向。當前體系里的苯環是在XY平面上

Grid(base) {
    type=even //格點是均勻分布的
    origin=[-8.0, -8.0, -8.0]   //格子原點
    ivec=[1.0, 0.0, 0.0]   //格子有ijk三個方向，這里定義i方向矢量
    jvec=[ 0.0, 1.0, 0.0]     //定義j定義方向矢量。k的方向無需定義，程序自動將i與j做叉乘得到
    lengths=[16.0, 16.0, 16.0]   //i,j,k三個方向的長度，單位是Bohr
#   spacing=[0.5, 0.5, 0.5]     //i,j,k三個方向的格點間距，單位是Bohr
    grid_points=[50, 50, 50]     //i,j,k三個方向的格點數。spacing和grid_points二者只能定義一個，因為是彼此矛盾的，所以此例把spacing注釋掉了
}

Advanced {
    spherical=off //這個別管它，但必須寫，否則結果異常
    diamag=on      //on代表考慮diamagnetic部分的貢獻
    paramag=on     //on代表考慮paramagnetic部分的貢獻
}

Essential {
acid=on //on代表輸出acid.vti文件
jmod=on //on代表輸出jmod.vti文件
}

將MOL和XDENS文件放在上一級目錄下，把gimic.inp放在當前目錄下，然后直接運行gimic命令就開始進行運算，很快就能算完。顯然，grid_points越大，或者spacing越小，計算耗時越高，但圖像也越精細。

算完后當前目錄下出現了以下文件，在本文文件包里3D目錄下也可以看到。
? mol.xyz：分子結構文件。不懂xyz格式的話看《談談記錄化學體系結構的xyz文件》（http://www.shanxitv.org/477）
? grid.xyz：也是分子結構文件，但是末尾還多了8個X原子，對應盒子的8個頂點。末尾還有一個Be原子，無意義，可刪掉
? acid.vti：ACID函數的格點數據。關于這個函數的介紹看《使用AICD程序研究電子離域性和磁感應電流密度》（http://www.shanxitv.org/147），它本質上反映的是相應位置電子對磁場感應的各項同性強度
? jmod.vti：電流矢量的模的格點數據（標量場），體現各個位置電流的大小
? jvec.vti：電流矢量的格點數據（矢量場）

為了能讓ParaView正確地顯示出分子結構，繪圖之前還需要產生分子的cml格式的文件。先在機子里安裝OpenBabel程序，運行yum install openbabel即可。然后將本文文件包里的xyz2cml.sh拷到當前目錄，運行./xyz2cml.sh，此腳本就調用OpenBabel把mol.xyz轉換成了mol.cml，然后再自動用awk工具把此文件里坐標單位從埃轉換成Bohr并輸出為mol-bohr.cml。由于vti文件是用原子單位記錄的，所以ParaView載入的也必須是以Bohr記錄坐標的cml文件。

現在相關文件都已經準備好了，可以繪制各種圖了。下面筆者依次介紹各種常見類型的圖怎么繪制。由于操作步驟很多，所以過程筆者都錄成了視頻而不以文字和圖片去介紹了。繪制各種圖的演示視頻見：https://www.bilibili.com/video/av56441881（全長23分半！）。

值得一提的是在ParaView里選項很多，往往需要做很多操作才能得到期望的圖像，設好之后建議用File - Save State將當前的狀態保存為pvsm文件，這樣以后隨時可以用Load State載入它來恢復之前的作圖狀態。

5.1 繪制感生電流的模的等值面

基于jmod.vti，效仿本文的視頻第1部分即可繪制出來，效果如下。分子外圍藍色部分是diatropic（滿足左手定則的）電流的模的等值面，等值面數值為0.05；內部紅色部分是paratropic電流的模的等值面，等值面數值為-0.05。

5.2 用VMD繪制感生電流的模的等值面

如《在VMD里將cube文件瞬間繪制成效果極佳的等值面圖的方法》（http://www.shanxitv.org/483）所示，用VMD可以得到效果絕佳的等值面圖，而用ParaView明顯達不到這個效果，然而VMD又不支持直接載入vti文件，怎么辦？2019-Jun-21及之后發布的Multiwfn都支持載入vti文件（只支持標量數據形式）。用Multiwfn載入jmod.vti后，選擇主功能13，再選擇導出cube文件，之后就可以用VMD按照博文里的做法來繪制了。Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載。

完整操作過程在本文視頻第2部分做了演示。繪制效果如下，可見非常漂亮。

5.3 繪制感生電流的模的截面著色圖

按照本文視頻第3部分即可繪制出感生電流的模的截面著色圖。在ParaView里可以設置無數多個截面，每個截面可以通過著色方式展示出標量函數的數值或者向量函數的模。下圖把分子平面和分子平面上方1.5 Bohr處的環電流的模同時通過顏色進行了展示。為了避免遮住下方的平面，上方的平面設置了透明效果。

5.4 繪制ACID等值面圖+向量場圖

按照本文視頻第4部分即可繪制出來向量場圖，可以展現整個三維空間中各個位置的電流方向。但由于倒處都是箭頭，為了看起來更清楚，還同時把ACID等值面也展現了出來。從此圖上也可以看出在分子平面上方、下方以及外圍，電流都是diatropic的。

5.5 繪制感生電流三維流線圖

基于jvec.vti，效仿本文的視頻第5部分即可繪制出來感生電流的三維流線圖，效果如下。

流線圖可以清晰地展現出環電流的路徑。繪制流線圖的時候，流線的起點是用戶通過設置一個球形范圍來定義的，球形位置可以用鼠標挪動，球形半徑、流線的長度和方向都可調，默認是前后兩個方向都產生流線。

還可以顯示成流線箭頭圖

由上圖可見，在苯環上方形成了顯著的環繞苯環的環電流。圖中流線是根據矢量的模來著色的，圖上箭頭都普遍偏橙黃，體現出環電流強度較大。本例磁場是往Z軸正方向加的，從箭頭方向上也看出在苯環上方形成的是滿足左手定則的diatropic電流。

繪制的流線圖也可以只往正方向或逆方向生成，比如把Integration Direction從默認的BOTH切換為FORWARD之后，圖像如下所示。磁場是從屏幕內側往外側加，可見在分子外圍區域都是diatropic環電流（順時針方向），而且在C-C sigma鍵區域也形成了顯著的局部環流，體現出了這個區域電子較高的定域性。

5.6 繪制感生電流的模的平面著色流線圖+箭頭圖

按照本文視頻第6部分即可繪制出來感生電流的模的平面著色流線圖，如下所示。此圖選的截面正好在分子平面上，圖中越白的地方說明電流大小越大，可見繞著每個sigma鍵都形成了顯著環電流。

下圖是把平面挪到距離分子平面1.5 Bohr的位置，用來反映pi電子形成的環電流的情況。為了看得更清楚，把色彩刻度上限調小到了0.2。可見，在碳環上方形成了一圈白色，體現出在這個環狀區域環電流很大。

之后可以再增加一個截面用來顯示箭頭，位置也設成距離平面1.5 Bohr處，這樣就把特定截面上電流的大小和方向同時清楚地展現出來了，如下所示。

5.7 繪制全空間動態流線場圖

這種圖的動畫效果極其酷炫，保證每個看到這個圖第一眼的人都會被驚到！這種圖以極度直觀的方式將外磁場下環電流的流動體現了出來。按照本文視頻第7部分即可輕易繪制出來，大家一定要看一下。這種動態圖的截圖如下。

6 實例：對穿越苯分子C-C鍵截面的電流進行積分

這一節講一下怎么計算穿越苯分子C-C鍵截面的環電流的積分值，這可以從定量上考察環電流的大小、考察芳香性的強弱。做這種計算需要定義一個積分截面，對于當前例子如下所示

可見需要設置兩個原子來定義鍵軸，定義截面距離第一個原子的距離，以及定義鍵軸與截面交點向四周的延展距離。為了唯一地定義截面區域，還需要設置一個固定點（fixed point），這個點的位置就取為相同環平面上的另一個原子即可，并且應當取的是相對于鍵軸靠環中心那一側的原子。

本例相關文件都在本文文件包里的int目錄下。其中輸入文件gimic.inp的內容如下，有必要說明的地方我都進行了注釋。

calc=integral
basis="../MOL"
xdens="../XDENS"
openshell=false
magnet=[0,0,1]

Grid(bond) {
    type=gauss   //在截面上分布用于高斯積分的格點
    gauss_order=9   //高斯積分階數。當前值總是恰當的，不需要改
    bond=[1,2]    //定義鍵軸的兩個原子序號為1和2
#   coord1=[0.0, 0.0, 2.145166]   //也可以通過兩個坐標來定義連線方向。注意單位是Bohr
#   coord2=[0.0, 0.0, -2.145166]
    distance=1.32 //截面距離第一個原子的垂直距離，單位為Bohr。1.32 Bohr是苯的C-C鍵的一半，通常都取為鍵的一半
    fixpoint=4 //固定點對應的原子的序號（4號原子位置見下文的圖）
#   fixcoord=[0.0, 0.0, 0.0] //也可以直接提供固定點的坐標
    grid_points=[30, 30, 0]   //只需要定義前兩個值，用來設定截面上兩個方向的點數。越大積分越準確，通常用30*30就夠了
    height=[-5.0, 5.0]   //對應于上圖的down和up數值。由于down對應下方，所以寫成負值。當前數值大小一般來說足夠大了
    width=[-2.2, 5.0]    //對應于上圖的in和out數值。由于in對應內側，所以寫成負值，其絕對值對應于苯環中心與鍵軸的垂直距離。當前out數值一般來說足夠大了
}

Advanced {
    lip_order=5     //拉格朗日插值多項式的階數。不用改
    spherical=off
    diamag=on
    paramag=on
}

確保MOL和XDENS在上一級目錄下，運行gimic命令執行此目錄下的gimic.inp。由于做截面積分需要算的點數很少，瞬間計算就完成了，主要輸出信息如下

*** Integrating current
Magnetic field <x,y,z> = 0.00000 0.00000 1.00000

************************************************************
    Induced current (au)    :     0.428903
       Positive contribution:     0.604900 ( 17.045730 )
       Negative contribution:    -0.175997 ( -4.959491 )

Induced current (nA/T) : 12.086239
(conversion factor) : 28.179409
************************************************************

即穿越C-C鍵截面的感生電流數值為0.4289 a.u.，還把diatropic和paratropic電流分別貢獻的正值和負值部分都給出了，括號里是常用的nA/T為單位的數值，即a.u.的數值乘上28.179409。

當前目錄下還出現了grid.xyz，其中的X原子對應截面的四個頂點位置。將其中最后一個意義不明的Be原子刪掉，然后在VMD里進行恰當設置，可得下圖。可見確實截面位置設置得正確。建議大家每次總是這么通過繪圖檢查一下截面設置到底對不對，避免得到無意義的結果。（具體來說繪制過程是用《在VMD中顯示原子序號的方法》http://www.shanxitv.org/197里的做法顯示出序號，然后在Graphics - Representation里新增Rep，選區設成name X，然后顯示風格設成DynamicBonds，然后把Distance Cutoff逐漸設大直到四個X之間連成上圖的框）

注意環電流積分的正負號是取決于鍵軸定義的方向的。當前是用1,2兩個原子定義鍵軸，如果反過來，寫成bond=[2,1]，則積分值就成了-0.4289了。使用bond=[1,2]時，順著1→2原子方向的電流對積分值是負貢獻，順著2→1原子方向的電流對積分值是正貢獻，由于此例總積分值為正，2→1方向的電流占主導，因此diatropic電流比paratropic電流整體要大，體現出苯分子的芳香性很強。

7 其它

目前研究環電流的程序主要有兩個，一個是本文說的GIMIC，一個是AICD，后者在《使用AICD程序研究電子離域性和磁感應電流密度》（http://www.shanxitv.org/147）和《使用AICD 2.0繪制磁感應電流圖》（http://www.shanxitv.org/294）都詳細介紹了。這兩個程序各有利弊。相對于AICD而言，GIMIC的優點在于可以給出三維空間中各個位置的環電流矢量和模，并且結合ParaView可以以不同方式作出很漂亮的圖進行考察，繪圖設定非常靈活。GIMIC的缺點在于沒法將環電流分解為不同軌道的貢獻，因此沒法實現sigma和pi環電流的分離考察，而且也沒法得到AICD程序給出的那種在AICD等值面上標記電場方向箭頭的圖（這種圖非常直觀，而且對于展現非平面體系的環電流特征尤為方便和好用）。所以GIMIC和AICD之間沒有替代關系，二者在展現環電流的角度上是高度互補的，二者都應當會用，并且根據實際情況判斷怎么展現最清楚、最能說明問題。

實際研究的環往往不是正好在XY、YZ或XZ平面上的，這給加磁場的設定、設置作圖平面帶來麻煩。對于這種情況，筆者建議用《調節平面分子間距的方法》（http://www.shanxitv.org/178）文中提供的筆者寫的VMD腳本alignplane來把要考察的環弄成平行于XY平面的狀態。之后Gaussian計算時必須加nosymm關鍵詞，否則朝向會被Gaussian自動改變。

我估計大部分讀者對ParaView都不熟悉，本文的視頻里已經充分體現了用ParaView可視化格點數據涉及的各種基本操作，希望讀者舉一反三，靈活運用，結合實際情況得到又酷炫又能很好說明問題的圖。另外，目前官網上最新的Multiwfn已經可以在計算或載入格點數據后利用主功能100的子功能2導出vti文件和以Bohr為單位的cml文件，因此各種實空間函數、cube文件也都可以在ParaView里可視化。

可能有讀者不知道怎么恰當設置盒子的原點（origin）和邊長（lengths），其實利用Multiwfn，很容易確定這兩個參數。啟動Multiwfn（這里用的是3.6版），載入fch/fchk文件，然后依次輸入
5
100
1
然后在屏幕上就可以看到這些信息
Coordinate of origin in X,Y,Z is     -10.078470 -10.709411   -6.000000 Bohr
Coordinate of end point in X,Y,Z is    9.993583   10.746921    6.112446 Bohr
Grid spacing in X,Y,Z is    0.346070    0.346070    0.346070 Bohr
這說明對于當前體系可以設origin=[-10.1, -10.1, -6.0]。X方向邊長適合設9.993583-(-10.078470)=20.1 Bohr，對Y、Z方向也類似，因此設lengths=[20.1, 21.4, 12.1]是合適的。