使用SYSMOIC程序繪制磁感生電流圖和計算鍵電流強度
使用SYSMOIC程序繪制磁感生電流圖和計算鍵電流強度
文/Sobereva@北京科音 2024-Feb-29
0 前言
筆者在《衡量芳香性的方法以及在Multiwfn中的計算》（http://www.shanxitv.org/176）中介紹過，磁感生電流是一種重要、直觀的考察分子芳香性的方式。SYSMOIC是一個基于continuous transformation of the origin of the current density (CTOCD)方法考察化學體系的磁感生電流相關問題的程序。本文將結合苯分子為例，介紹這個程序的最關鍵用法。SYSMOIC程序和相關原理的詳細介紹可以看原文J. Chem. Inf. Model., 61, 270 (2021)。
AICD和GIMIC也是兩個重要的考察磁感生電流的程序，筆者在下文有專門的介紹。這些文章里介紹的背景知識對理解本文很有幫助
? 使用AICD 2.0繪制磁感應電流圖（http://www.shanxitv.org/294）
? 使用AICD程序研究電子離域性和磁感應電流密度（http://www.shanxitv.org/147）
? 考察分子磁感生電流的程序GIMIC 2.0的使用（http://www.shanxitv.org/491）
? Utilizing AICD and GIMIC programs to study magnetically induced current density（http://www.shanxitv.org/148）
之所以筆者寫此文介紹SYSMOIC，在于它對以上兩個程序有明顯互補性。AICD程序是將感生電流以箭頭方式繪制在AICD的等值面上，而等值面以外區域的感生電流就無法在圖上體現出來了，而且AICD沒法給出感生電流的鍵截面積分值（鍵電流強度）來定量考察電流大小。SYSMOIC則能夠給出三維空間或二維平面上的感生電流向量場圖，從而能展現出AICD圖表現不出來的信息，而且SYSMOIC還能非常方便地給出鍵電流強度。GIMIC程序雖然能做感生電流積分、結合Paraview程序能給出漂亮的向量場圖，但是GIMIC卻無法將感生電流分解成不同分子軌道的貢獻，而SYSMOIC則能夠方便地指定計算感生電流時只考慮或不考慮哪些分子軌道。可見SYSMOIC有獨特的價值。
本文介紹的情況對應SYSMOIC 2024.1版，對其它版本不一定適用，請結合相應版本程序提示和手冊隨機應變。
1 SYSMOIC的基本特征
SYSMOIC是免費而不開源的程序，Windows、Linux、Mac版都有，可以在http://SYSMOIC.chem.unisa.it/DISTRIB下載。諸如STABIN_LNX-2024.1.tar.gz就是Linux的2024.1版可執行文件包。Linux版需要用較新的系統，筆者用Rocky Linux 9可以正常使用，CentOS 7.4用不了。
SYSMOIC是由一堆子程序組成的，對應壓縮包里不同可執行文件，Windows版有.exe后綴，其它平臺的沒有。例如JBMAP子程序用于計算感生電流格點數據，BOCUST子程序用于計算鍵電流強度，等等。v3d子程序用來對其它子程序計算出來的.3d文件進行可視化，如繪制感生電流矢量圖。各個程序的用法在SYSMOIC的在線手冊http://sysmoic.chem.unisa.it/MANUAL/里都有說明。
SYSMOIC是解壓即用的程序。解壓后只需要把SYSMOIC目錄加入到操作系統的PATH環境變量里（不會的話自行google），然后就可以在任何目錄下輸入子程序名直接使用了。
SYSMOIC有不少官方例子文件，是http://sysmoic.chem.unisa.it/WORKED_EXAMPLES/里的WORKEXA.tar.gz。例子按照體系/理論方法/基組對目錄層級進行命名，里面RUN_TESTS文件是對例子體系做各種分析的shell腳本，可以從中看到都是以什么命令調用程序的，以及傳入的命令文件是什么。REFERENCES目錄下的文件是運行RUN_TESTS腳本后產生的.3d文件和程序輸出信息，對前者可以直接用v3d作圖。
SYSMOIC必須結合最主流的量子化學程序Gaussian使用，目前只支持閉殼層形式的HF和DFT。只需要用Gaussian做NMR計算的同時產生wfx文件，再用SYSMOIC自帶的unpackwfx將之轉化成fort.3、fort.11、fort.23、fort.28這四個二進制文件，之后就可以用JBMAP、BOCUST等子程序做各種計算了，算完了可以用v3d子程序直接作圖。
2 產生苯的SYSMOIC輸入文件
下面以最典型的芳香性分子苯為例，講解SYSMOIC繪制感生電流矢量圖和計算鍵電流強度的方法。其它的用法我覺得在一般研究中用到的機會很少，所以就不在本文介紹了，感興趣者請看程序手冊、原文并結合程序自帶的例子文件摸索。下文涉及的各種輸入輸出文件都可以在http://www.shanxitv.org/attach/702/file.rar中獲得。
首先寫一個在B3LYP/6-31G*級別下做NMR計算的Gaussian輸入文件benzene.gjf，內容如下，在本文的文件包里已經提供了。其中的坐標是在當前級別下已經優化過的。NMR=CSGT要求以CSGT方法做NMR計算，out(wfx,CSGT)要求產生.wfx波函數文件，并且把CSGT計算過程中產生的擾動的軌道展開系數寫入到wfx文件中，輸出位置為末尾指定的D:\benzene.wfx。
# b3lyp/6-31g(d) NMR=CSGT out(wfx,CSGT)
[空行]
benzene
[空行]
0 1
C 0.00000000 1.39621600 0.00000000
C 1.20915900 0.69810800 0.00000000
C 1.20915900 -0.69810800 0.00000000
C 0.00000000 -1.39621600 0.00000000
C -1.20915900 -0.69810800 0.00000000
C -1.20915900 0.69810800 0.00000000
H 0.00000000 2.48297800 0.00000000
H 2.15032200 1.24148900 0.00000000
H 2.15032200 -1.24148900 0.00000000
H 0.00000000 -2.48297800 0.00000000
H -2.15032200 -1.24148900 0.00000000
H -2.15032200 1.24148900 0.00000000
[空行]
D:\benzene.wfx
[空行]
[空行]
算完后就有了benzene.wfx，確保它在當前目錄下，并且已經把SYSMOIC目錄加入到了PATH環境變量中，然后運行unpackwfx benzene命令，當前目錄下立刻就出現了fort.3、fort.11、fort.23、fort.28，它們都已提供在了本文的文件包里。后文說的所有操作都是對于這四個fort文件都處于當前目錄下的情況而言的。
注意當前的苯分子是處在Z=0的XY平面上，用Multiwfn載入Gaussian的out文件然后進主功能0，或者用GaussView載入out文件并顯示笛卡爾軸，都可以確認這一點。后面我們要考察的都是外磁場方向垂直于苯環所導致的感生電流的情況。
3 繪制磁感生電流向量場圖的例子
3.1 默認方式繪制
先來看怎么產生一個最簡單的磁感生電流量場圖，之后再解釋細節。直接在命令行窗口輸入JBMAP運行之，屏幕上會出現好多默認用的參數，比如下面這行，就是告訴你默認的外磁場是(0 0 1)矢量，顯然正垂直于當前苯環，和我們期望的一致
magnetic field...........B 0.0000 0.0000 1.0000
程序問你當前顯示出來的參數是否ok，這里就直接按回車代表用默認的y（yes）。然后再依次問你三個問題：
? external grid points y/n? [n]：即是否讀取外部格點，這里直接按回車代表用默認的n（如果要計算的點不是均勻分布在一個矩形區域，而是比如分布在一個球層表面，就需要輸入y并且提供格點定義文件）
? reference arrow y/n? [n]：含義不明。直接按回車用默認的n
? B direction y/n? [n]：即是否把外磁場箭頭畫在圖上，直接按回車用默認的n。如果你想把外磁場箭頭繪制在圖上就輸入y，并且輸入箭頭的位置和長度
馬上當前目錄下就有了JBM.3d文件。輸入v3d JBM.3d通過v3d程序對這個文件作圖，就蹦出了圖形窗口，如下所示，可以看到在一個矩形區域內均勻分布的每個格點位置都有一個箭頭指示此位置的磁感生電流的方向和大小。
想退出的話在圖形窗口按ESC鍵，或者在文本窗口按ctrl+C。
上面的圖非常丑陋，箭頭特別擁擠，分子結構也看不清楚，完全沒法用。下一節講怎么對JBMAP的計算進行自定義。
3.2 自定義JBMAP的運行參數
JBMAP有兩個方式指定運行參數，比如要把外磁場矢量改為0 1 1（程序自動會做歸一化成為0.0000 0.7071 0.7071），既可以運行JBMAP -B 0 1 1命令，也可以先啟動JBMAP然后直接輸入B 0 1 1后按回車，在屏幕上顯示的參數中都會看到外磁場矢量已經變更成自設的值了。
JBMAP命令后面可以直接接的選項可以輸入JBMAP -h查看，其中幾個關鍵的在這里提一下：
-o：輸出的.3d文件名，默認是-o JBM
-f：原子球大小的倍數，默認為-f 0.2。想讓作出的圖里原子球大一些就把這個改大
-B：外磁場矢量，默認為-B 0 0 1
-j：感生電流類型，默認是-j TOT，即diatropic和paratropic的總和。寫-j PARA代表只考慮paratropic部分，寫-j DIA代表只考慮diatropic部分
-q：考慮的軌道。例如只想考慮17 20 21三個分子軌道對磁感生電流的貢獻，就寫-q +3 17 20 21，只去掉15 16號分子軌道的貢獻就寫-q -2 15 16。默認是考慮所有分子軌道的貢獻（只有占據軌道才有貢獻）
-m：設置SYSMOIC計算感生電流用的CTOCD方法的具體形式，參看http://sysmoic.chem.unisa.it/MANUAL/S4.SS12.html的說明。當基組是完備的時候這些方法結果等價。有些選項是僅對于某些方法才需要設，比如用common gauge-origin的時候需要用-g自定義原點坐標（默認為0 0 0）。默認用的CTOCD-DZ2基本上對于所有情況都是合適的選擇，也不用設額外參數。關于CTOCD不同形式的介紹和差異，可以看前述的SYSMOIC原文以及感生電流綜述文章WIREs Comput Mol Sci, 6, 639 (2016)的Continuous Transformations of the Origin of the Current Density一節。
JBMAP運行后可以修改的參數當中較重要的如下：
? B：外磁場矢量，默認B 0 0 1
? JTERM：感生電流類型，默認是JTERM 0，即diatropic和paratropic的總和。1和2分別代表paratropic和diatropic
? RI和RF：分別指定格點所均勻分布的矩形盒子的X、Y、Z最小笛卡爾坐標和最大笛卡爾坐標。默認是根據體系坐標自動確定以覆蓋整個體系。對上面的苯分子的例子，從屏幕上默認的參數可見，當前對應的是RI -6.0635 -6.6921 -2.0000和RF 6.0635 6.6921 2.0000
? FATT：給感生電流乘的系數。比如大部分箭頭在圖上都顯得很不顯著，就可以設大這個值以更便于觀看
? FMM：磁感生電流過濾范圍。例如FMM 0.01 1就代表忽略感生電流的模小于0.01和大于1的格點（是對乘上FATT之前而言的），免得干擾觀看。默認0 0.1
? STEP：格點間距，數值越小箭頭越密、計算越耗時。默認0.4 Bohr
? LUNA：0和1分別代表箭頭的長度和面積正比于感生電流的模，默認是后者
? METH：計算方法，默認的2對應CTOCD-DZ2
為了方便，我推薦以命令行方式運行JBMAP。交互式程序怎么以命令行方式運行，我在《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612）里有很詳細說明。下面以這種方式結合自定義參數用JBMAP再次對苯進行計算。創建一個文本文件叫CDMAP.txt，內容如下，每一行對應一條傳給JBMAP程序的命令。此例RI和RF的Z值都是1，X和Y范圍都是-6到6 Bohr，因此格點會分布在苯分子上方1 Bohr的XY平面上的正方形區域內。此文件后面的y、n、n、n對應交互式界面里程序向你提問時依次輸入的四個指令。
B 0 0 1
RI -6 -6 1
RF 6 6 1
FMM 0.01 1
FATT 3
STEP 0.4
LUNA 1
JTERM 0
y
n
n
n
運行JBMAP -f 0.5 < CDMAP.txt，然后再運行v3d JBM.3d，看到下圖，可見清楚地把苯分子上方的磁感生電流描繪了出來。
再來看只考慮pi軌道時的磁感生電流圖的繪制。把benzene.wfx文件載入Multiwfn，進入主功能0，按照《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）說的，從HOMO開始往下挨個看軌道，就能找到三個pi軌道，序號為17、20、21。原子數多的時候也可以用《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）介紹的做法讓Multiwfn自動找出來所有pi軌道序號。運行JBMAP -q +3 17 20 21 -f 0.5 < CDMAP.txt，然后再用v3d繪圖，得到下圖，這便是苯的pi電子的感生電流圖了。相對于之前的所有電子貢獻的圖，此圖少了sigma電子造成的出現在內環的逆時針（paratropic電流）特征，而且也沒有了氫原子區域的感生電流，因為pi電子不在氫上。
3.3 v3d的使用細節
下面專門細談一下v3d界面的操作。v3d圖形窗口出現后，可以用鼠標拖動體系旋轉視角。還可以在圖形窗口處于激活狀態下按h鍵，此時文本窗口會出現可以用的所有按鍵的功能說明。其中較重要的按鍵在這里說一下。
先說視角調節按鍵：
a、b、g：視角由alpha、beta、gamma值定義，按a、b、g鍵會分別增加它們1.0。如果按住shift鍵再按，也即按A、B、G，會分別增加它們10
? f和F：調節縮放因子。按一下f圖像就會放大一點，按一下F就會放大很多
? d和D：調節視角距離。按一下d就會增加0.1，按一下D就會增加1。視角距離越遠，近大遠小的透視效果越弱、越接近正交視角的效果（據我所知v3d沒法完全用正交視角）
? n和N：設置霧化距離閾值。超過視角一定距離閾值開始就會用霧化效果使得遠處物體被屏蔽。按n增加閾值0.1，按N增加1.0
? C：令圖像居中
每次按以上按鍵調節設置的時候在文本窗口都會看到當前值。如果你想把這些鍵的操作效果反過來，比如按一次f就令圖像縮小一點，那就先按一下c鍵進入反轉操作狀態（文本窗口里會看到c變成了c-1），然后再按f就行了。再次按c則還原到普通操作狀態。
按p鍵可以把當前圖像保存成當前目錄下的.ps文件。這是矢量圖像格式，在Windows下可以用裝了ghostscript插件的irfanview或photoshop等程序打開，也可以用Acrobat distiller轉成pdf再打開。在Rocky Linux、CentOS等系統里也自帶了ps文件觀看程序。
默認情況下產生的ps文件畫出來的感生電流的箭頭會有較厚的白邊，當箭頭多的時候白邊會導致后面的箭頭被遮擋，比較難看。解決辦法是按一下c鍵先進入反轉操作狀態，然后按住k鍵一會兒使得白邊量減小（文本窗口不會提示減小到了多少），然后再導出ps文件，然后就會發現白邊非常細了。
如果想保存當前的作圖設置，按S，當前目錄下就出現了SpecialValues.txt文件，里面記錄了作圖設置，可以看到f、d等數值都記錄在內。當再次啟動v3d時，程序若發現當前目錄下有這個文件，就會自動載入里面的設置。如果你希望每次產生的ps文件里箭頭的白邊都非常細，可以啟動v3d后先直接按S保存一個記錄默認狀態的SpecialValues.txt文件，然后把里面對應白邊粗細的ridfr前面的值改為0.1。這樣每次啟動v3d就會默認用這個設置了。
4 計算鍵電流強度
對鍵截面做感生電流積分的概念我在《考察分子磁感生電流的程序GIMIC 2.0的使用》（http://www.shanxitv.org/491）中已經專門講過了，這可以定量考察穿越特定化學鍵的電流強度從而便于對比，這在SYSMOIC里叫鍵電流強度。SYSMOIC算鍵電流強度遠比用GIMIC更方便，都省得自己定義積分范圍，結果還能直接可視化。下面我們對苯的一個C-C鍵計算鍵電流強度。
直接在命令行窗口中輸入BOCUST命令運行此程序，程序首先會提示
| plotted 12 atoms
| plotted 12 bonds
====> Pair C1 - C2 distance= 1.396 Ang
這說明程序發現有12個原子對兒的距離都處在默認的距離下限（0.9埃）和距離上限（1.6埃）之內。其中C1-C2是第一個原子對，現在問你是否對這個C1-C2鍵計算鍵電流強度。這里輸入y然后按回車，很快就會出現結果
CONTRIBUTION 1 -0.6000526E+00 AU -0.1690912E-07 SI
CONTRIBUTION 3 0.1784109E+00 AU 0.5027513E-08 SI
TOTAL CS ==> -0.4216417E+00 AU -0.1188161E-07 SI
(2) CS/CS0 RELBEN FICS 0.99 3.00 0
上面顯示的兩個CONTRIBUTION分別是diatropic和paratropic對鍵電流強度的貢獻，由于方向相反，顯然符號相反。二者之和取絕對值0.4216417 a.u.就是鍵電流強度大小，換算成常用的nA/T單位就乘以28.179409，即11.88 nA/T。這個值與《考察分子磁感生電流的程序GIMIC 2.0的使用》（http://www.shanxitv.org/491）在B3LYP/6-311G*級別通過GIMIC程序得到的12.086 nA/T非常接近。上面信息中的CS/CS0 RELBEN FICS 0.99代表當前值與參考值的比值是0.99。這里參考值是程序內置的對苯在某個級別計算的C-C鍵電流強度12 nA/T。
接下來程序又問你是否對C1-C6計算鍵電流強度。由于當前所有C-C鍵都是等價的，只需要計算一個就行了，因此輸入N回車，程序就結束運行了。如果你想把所有鍵都進行計算，就接著輸入11次y回車把剩下11個鍵都算完。當前體系一共6個C-C鍵、6個C-H鍵。
當前目錄下出現了BCS.3d文件。運行v3d BCS.3d，就可以看到下圖
此圖里黑色箭頭體現了C1-C2鍵電流的大小和方向。積分截面處在鍵的正中央并與鍵垂直，圖中的曲線是這個截面上的感生電流0.001 a.u.等值線，曲線內部區域是對感生電流積分的區域，紅色和藍色分別對應paratropic和diatropic部分。可見diatropic分布在環的外側，這與之前從向量場圖中看到的順時針環電流的主要分布是一致的。圖中的紅色小點對應截面上感生電流的極值點。鍵上面標的文字99對應于前述的參考比值0.99乘以100。文字出現的位置在鍵的正中央，不容易分辨，可以在圖形窗口激活狀態下按W鍵關閉文字顯示，然后再自行ps上。
如果你想一次性把所有鍵的電流強度都計算，又不想一次次麻煩地輸入y，那么可以寫一個文本文件比如叫ally.txt，里面有12行，每一行都是y。運行BOCUST < ally.txt，這樣所有鍵就都會被計算了。此時對得到的BCS.3d繪制會看到下圖
下面介紹一下BOCUST運行時可以接的選項中比較重要的：
-o、-f、-B、-m、-q等：和前面介紹的JBMAP程序的情況完全相同。即BOCUST也可以用-B指定外磁場方向、用-q指定考慮哪些軌道，等等。
-d、-D：確定哪些原子對兒被考慮時用的距離下限和上限（埃），默認-d 0.9 -D 1.6
-e：確定鍵截面上做積分的區域用的感生電流等值線的數值（a.u.），默認-e 0.001
-CS0：參考電流強度（nA/T），默認-CS0 12
-nocd：不顯示積分區域和極值點。如果想讓圖上只有箭頭而沒有曲線和小點，應該加上這個
-l [比值]：在截面上繪制感生電流的不同數值的等值線，后面跟的是相鄰等值線數值的比值
-C [粗細]：繪制截面的邊框
例如這里計算C1-C2的sigma電子的鍵電流強度（即扣除三個pi軌道貢獻），并在截面上顯示等值線、顯示邊框，就運行BOCUST -q -3 17 20 21 -l 3 -C 5，然后輸入y回車，再N回車。之后用v3d對BCS.3d作圖看到下圖。可見整個截面，以及里面的感生電流等值線，都顯示得很清楚。由于這個截面上diatropic和paratropic的大小都差不多（畢竟苯分子沒有sigma芳香性特征），所以鍵電流強度僅為0.007 a.u.。
5 總結
此文對分析磁感生電流的程序SYSMOIC做了介紹并給出了具體例子。由例子可見，通過SYSMOIC結合Gaussian繪制感生電流向量場圖很方便，不需要借助第三方程序。SYSMOIC做鍵電流強度計算也很自動化，不需要像GIMIC那樣謹慎地定義積分區域，而且鍵電流大小和方向還能直接通過箭頭直觀顯示出來。SYSMOIC的計算速度也比較快。
本文只涉及了SYSMOIC最常用的功能，還另外有一些其它功能，比如對感生電流密度做拓撲分析，計算各種感生電流密度相關的場（ACID及其各向異性、磁屏蔽密度、vorticity張量的跡及其各向異性等）并繪制成等值面和填色平面圖、計算任意格點上的感生電流張量和導數等等。