• 詳談Multiwfn的命令行方式運行和批量運行的方法

    詳談Multiwfn的命令行方式運行和批量運行的方法

    文/Sobereva@北京科音  2021-Aug-30


    0 前言

    強大、免費的波函數分析程序Multiwfn(http://www.shanxitv.org/multiwfn)已經被計算化學研究者使用得非常普遍。雖然Multiwfn通常是以交互式方式使用,一次分析一個體系,但在Multiwfn手冊5.2、5.3節也已經明確寫明了怎么通過命令行方式使用、怎么通過腳本一次算一批體系,我的大量Multiwfn相關博文里也都用了腳本非常方便地實現自動和批量分析,但還是總有用戶問一些相關的初級問題。因此,筆者決定在本文專門十分詳細且深入淺出地講解一下Multiwfn的命令行方式運行和批量運行的方法,并給出許多類型各異的例子,令讀者充分了解到靠腳本運行Multiwfn的極度便利。本文的很多內容對于通過腳本運行其它計算化學程序也是共通、適用的,某種意義上也算是個計算化學腳本編寫入門教程。

    如果讀者還沒看過這三篇,請看一下了解Multiwfn的相關常識:《Multiwfn入門tips》(http://www.shanxitv.org/167)、《Multiwfn FAQ》(http://www.shanxitv.org/452)、《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/379)。筆者假定讀者有最基本的Windows和Linux常識。如果本文有些術語、命令你不清楚,請自行Google。

    需要強調的是,做計算化學的人,哪怕不會Fortran、C、Python等編程語言,至少也該會寫shell腳本。諸如用Multiwfn以相同的流程分析三、四個體系,那沒必要用腳本;處理10個體系,不用腳本也能接受;但如果要處理幾十、上百個體系,不利用腳本是絕對不可能的,手動一個個操作不僅極度耗時,無異于原始人,還容易在疲憊時中途不慎弄錯數據。如果要用Multiwfn對大批體系算一堆描述符(見《Multiwfn可以計算的分子描述符一覽》http://www.shanxitv.org/601)用于機器學習等目的,不會寫腳本完全不可能實現。本文涉及的腳本大多是以Bash shell(主流Linux系統默認的shell)的語法寫的,既可以在Linux下用,也可以在Windows里通過cmder模擬的Bash環境下使用(也可以在VMware等虛擬機里虛擬的linux系統下或WSL、Cygwin等環境下使用)。如果你不了解cmder的話,看《量子化學程序ORCA的安裝方法》(http://www.shanxitv.org/451)里面的cmder介紹,十分方便。Bash shell腳本編程很容易學,如果沒有相關知識的話建議看看入門教程https://linuxconfig.org/bash-scripting-tutorial-for-beginnershttps://en.wikibooks.org/wiki/Bash_Shell_Scripting,頂多兩個晚上就足夠學明白。另外還有個Linux下常用命令和Bash腳本語法的速查表可以看看:https://github.com/skywind3000/awesome-cheatsheets/blob/master/languages/bash.sh。本文Bash腳本里涉及到的Linux下的各種命令的用法都可以在這里查到清晰的介紹:https://man.linuxde.net

    在Windows下也可以寫batch腳本(批處理腳本,以.bat為擴展名)實現批量執行,但batch腳本靈活程度遠不及Bash shell腳本,所以在本文中只是很簡略涉及。本文使用的Windows下的文本終端(命令行環境)是cmd,而Powershell完全不在本文中考慮,筆者感覺設計得實在太爛。

    下文第1節先介紹怎么通過命令行方式運行Multiwfn,在第2節會給出將命令行方式運行與腳本相結合的一些實際例子,第3節介紹批量運行Multiwfn的方法,第4節給出通過腳本實現Multiwfn批量計算的一些實際例子。讀者務必在看懂本文的例子后舉一反三,靈活利用腳本將大批量或過程復雜的分析盡可能變得自動化,避免當原始人。


    1 Multiwfn的命令行方式的運行方法

    通常運行Multiwfn時都是先啟動Multiwfn,然后根據屏幕上的提示手動一行行輸入指令。如果你經常要對不同體系做同一種分析,而且每次分析時要輸入的命令都是相同的話,那么為了避免每次分析時都重新敲一遍鍵盤,十分建議通過命令行方式運行,這樣只需要一行命令就可以實現分析目的。老有些特別初級的用戶嫌Multiwfn不是完全圖形界面的程序,實際上正是Multiwfn的文本+圖形界面混合的設計,才令用戶可以方便地以純命令行方式運行Multiwfn,并進而嵌入到腳本里實現自動化。

    順便普及個常識:在Windows下進入cmd命令行窗口是按Win+R鍵,然后輸入cmd,之后可以通過cd命令進入特定目錄。筆者很建議按照http://bbs.keinsci.com/thread-22940-1-1.html的做法給Windows的右鍵添加“在此處打開命令行窗口”選項,這樣進入cmd后直接就處在當前目錄下。


    1.1 命令行窗口下Multiwfn的啟動

    先說一下命令行窗口下Multiwfn的啟動。在Linux環境下,如果你嚴格按照Multiwfn手冊2.1.2節的說明進行了安裝,就可以在任何地方直接輸入Multiwfn命令啟動。如果你是在Windows的cmd命令行窗口下,且當前所在目錄就是含有Multiwfn.exe和配置文件settings.ini的目錄,也是直接輸入Multiwfn命令就可以啟動;如果當前是在其它目錄下,需要進入Windows的“高級系統設置”- “高級”- “環境變量”,將Multiwfn.exe所在目錄加入到Path變量里,并且新建一個Multiwfnpath變量,將變量值設為settings.ini的所在目錄。之后在cmd里,無論當前處在哪個目錄,也都可以通過輸入Multiwfn命令啟動了。

    Multiwfn以命令行方式啟動時可以結合輸入文件名以及選項/參數,這在Multiwfn手冊2.2節里說了。格式是:Multiwfn [輸入文件路徑] [選項/參數]。選項/參數沒寫的話就會用默認值,可以指定的有:
    -nt [線程數]:并行運算的線程數。也可以通過settings.ini里nthreads參數設默認值
    -uf [序號]:用戶自定義函數(user-defined function)的序號。也可以通過settings.ini里iuserfunc參數設默認值
    -isilent:開啟silent模式,此時所有原本會自動蹦出圖形窗口的場合就不會自動蹦出了,這樣在批量運行Multiwfn時就不用手動一次次關閉圖形窗口了。也可以通過settings.ini里isilent參數設默認值
    -set [路徑]:settings.ini文件的路徑。沒設的話會先在當前目錄下尋找settings.ini,如果找不到則會在Multiwfnpath環境變量設定的目錄下找settings.ini

    例如可以這么運行:Multiwfn /nico/COCl2.fch -nt 36 -set /sob/tmp/settings.ini -silent。輸入文件如果沒寫的話,程序會在啟動后提示你輸入。

    1.2 以命令行方式用Multiwfn做分析

    在cmd窗口或Linux的終端中,可以通過重定向方式將某個文本文件里記錄的指令傳遞給Multiwfn,文本文件里每一行的內容就是在Multiwfn里輸入的一次指令,如果指令是直接敲回車就在相應位置寫個空行。例如,要計算拉普拉斯鍵級,并且將所有原子間的鍵級矩陣導出成文本文件,在Multiwfn中載入輸入文件后需要輸入的命令是
    9  //鍵級分析
    8  //拉普拉斯鍵級
    y  //把鍵級矩陣導出為當前目錄下的bndmat.txt

    因此,如果你要對比如/sob/Rio.molden文件計算拉普拉斯鍵級并得到鍵級矩陣文件,就用文本編輯器創建一個文本文件,例如叫LBO.txt,內容是以下三行
    9
    8
    y
    然后運行Multiwfn /sob/Rio.molden < LBO.txt,運行完畢后當前目錄下就有了我們要的bndmat.txt。

    Multiwfn運行時輸出的信息都顯示在屏幕上,如果你想把其中一部分保存下來,按照Multiwfn手冊5.4節說的,直接在文本窗口選中一個區域復制,然后粘貼到文本編輯器或Excel等程序里即可。也可以把輸出在屏幕上的所有信息都重定向到一個文本文件里,比如保存到out.txt就輸入Multiwfn /sob/Rio.molden < LBO.txt > out.txt。如果想讓輸出信息既顯示在屏幕上也保存到文本文件里,就用|tee,即Multiwfn /sob/Rio.molden < LBO.txt |tee out.txt。

    在Bash shell下還有一種實現上面例子的做法,是創建一個文本文件,比如叫foo.sh并放在當前目錄下(習俗上Bash腳本文件用.sh后綴),內容寫
    Multiwfn /sob/Rio.molden << EOF > out.txt
    9
    8
    y
    EOF
    然后給腳本用chmod +x ./foo.sh加上可執行權限,之后輸入./foo.sh運行之即可。這個腳本從第二行開始的內容會被作為命令依次傳遞給Multiwfn,直到EOF那行為止(不是必須叫EOF,腳本里兩處EOF改用其它名字也可以。用EOF,即end of file,是常見習俗)。

    在Bash shell下還有更簡單的通過命令行方式運行Multiwfn的做法,就是通過管道直接把要執行的命令傳遞給Multiwfn,這樣就免得創建記錄執行指令的文本文件了。例如上例可以通過直接輸入以下命令達到相同目的
    echo -e "9\n8\ny" | Multiwfn /sob/Rio.molden > out.txt
    這里9\n8\ny是被傳遞給Multiwfn的指令,每個\n代表按一次回車(最后一次的不用明確寫),因此傳過去的指令就是9[回車]8[回車]y[回車]。

    如果你既不想讓Multiwfn計算過程輸出的信息顯示在屏幕上,也不想存到某個文件里,可以把上面的> out.txt改為> /dev/null(對于Linux)或改為> NUL(對于cmd),這即是把輸出信息定向到了垃圾桶設備里。

    大家肯定發現了,在按照上面的做法在運行結束后,屏幕上會看到諸如這樣的信息
    forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
    Image              PC                Routine            Line        Source

    Multiwfn.exe       00007FF6BB5B938F  Unknown               Unknown  Unknown
    Multiwfn.exe       00007FF6BB56C980  Unknown               Unknown  Unknown
    ...略
    這是因為你沒有優雅地(gracefully)退出Multiwfn,而相當于運行完特定指令后直接強制把Multiwfn關了,這屬于意外退出,因此最后會輸出一堆報錯。這些信息完全不用管它!因為我們想達到的目的已經達到了。但是如果你嫌它們礙眼,不想輸出這堆破信息,那也可以加入額外的指令來優雅地退出,也就是在Multiwfn中先返回主菜單,然后輸入q。具體來說,上面那個例子改成
    echo -e "9\n8\ny\n0\nq" | Multiwfn /sob/Rio.molden > out.txt
    這里傳過去的指令中的0是從鍵級分析功能退回到主菜單的選項,q是用來優雅地退出。還一個做法是不改傳入的指令,但是利用2>把上面的報錯信息重定向到垃圾桶里從而避免顯示出來:
    echo -e "9\n8\ny" | Multiwfn /sob/Rio.molden > out.txt 2> /dev/null
    小常識:>是用來重定向標準輸出信息的,包括Multiwfn程序本身輸出的所有信息。2>是用來重定向標準錯誤信息的,像上面那些由操作系統輸出的報錯信息就屬于這類。如果標準輸出和錯誤輸出信息想全都重定向到某個文件里,最簡單的做法是用&>,但cmd不支持這個,cmd里需要寫成> out.txt 2>&1(在Linux里也可以這么寫)。

    Multiwfn手冊4.4節里有通過主功能4繪制平面圖的大量例子,顯然你也可以用上面的做法通過一行命令就讓Multiwfn的主功能4作出圖來,在主功能4后處理菜單里直接給了保存圖像文件的選項。但別忘了要么把settings.ini里的isilent設為1,要么在運行指令后面加上-silent,否則平面數據算完后會自動在屏幕上蹦出圖像,到時候還得手動關閉。

    通過以上介紹,可以看出通過命令行方式使用Multiwfn相當容易。想通過命令行方式實現某種分析時,你只需隨便找個簡單小體系在Multiwfn的交互式界面里操作一遍,把要輸入的命令都記下來,就自然知道通過命令行方式運行怎么實現了。


    2 通過腳本運行Multiwfn的一些實例

    下面提供幾個通過腳本運行Multiwfn的實例,體現出利用腳本做分析的極度便利性,其中利用了上面介紹的命令行方式運行Multiwfn的知識。

    2.1 一鍵實現晶體密度預測

    在《使用Multiwfn預測晶體密度、蒸發焓、沸點、溶解自由能等性質》(http://www.shanxitv.org/337)的1.1節筆者介紹了怎么利用Multiwfn基于分子表面靜電勢描述符預測中性分子的晶體密度,讀者請先把那篇文章里的計算流程看一遍,原理不再累述。這里我們把計算流程寫成一個Bash腳本preddens.sh,用戶只需運行比如./preddens.sh RDX.fchk命令就能用這個腳本調用Multiwfn對RDX.fchk做分子表面靜電勢分析,并自動提取數據、按照預測公式給出預測的密度值,省得自己提取數據和手算了。此腳本內容如下,也可以在http://www.shanxitv.org/attach/612/preddens.sh下載。

    #!/bin/bash
    echo "Running:" $1
    Multiwfn $1 << EOF > out.txt
    12
    0
    -1
    -1
    q
    EOF

    t1=$(grep "Estimated density" out.txt | awk -F : '{print $2}' | awk '{print $1}')
    t2=$(grep "Product of sigma^2_tot and nu" out.txt | awk -F \( '{print $2}')
    echo "M/Vm:" $t1 "g/cm^3"
    echo "nu*sig2tot:" $t2 "(kcal/mol)^2"
    echo "$t1 $t2" | awk '{ print "Predicted density: " a*$1+b*$2+g " g/cm^3"}' a=0.9183 b=0.0028 g=0.0443
    rm -f out.txt

    使用./preddens.sh RDX.fchk命令進行計算會看到以下輸出
    Running: RDX.fchk
    M/Vm: 1.8045 g/cm^3
    nu*sig2tot: 26.08649 (kcal/mol)^2
    Predicted density: 1.77441 g/cm^3
    腳本首先告訴用戶當前算的是哪個文件,算完后會從Multiwfn輸出文件里把M/Vm和νσ2tot兩項提取并顯示出來,最后利用這兩個值和擬合參數按照預測公式給出預測的密度。

    這里來詳細逐行解讀一下腳本內容:
    ? #!/bin/bash代表將這個腳本當做Bash shell腳本來執行,使用/bin/bash作為語法的解釋器
    ? echo "Running:" $1里的echo命令用來將后面的內容輸出。變量前頭要帶著$來引用,$1變量對應的是這個腳本接收到的第1個參數,也即上例的RDX.fchk。如果腳本還接上了更多的參數,分別對應$2、$3...變量。$#變量是傳入的參數總個數
    ? Multiwfn $1 << EOF > out.txt這行和下面幾行內容代表用Multiwfn執行腳本后面接的那個文件,信息輸出到當前目錄下的out.txt里。往Multiwfn里傳入的命令12代表進入主功能12定量分子表面分析功能,后面的0是開始做分析,接下來兩個-1對應返回主菜單,最后q是優雅地退出Multiwfn。不清楚這些指令的含義的話自己在Multiwfn交互式界面里照著敲一遍,看屏幕上的提示自然就明白了
    ? 之后t1=$(grep "Estimated density" out.txt | awk -F : '{print $2}' | awk '{print $1}')這行略微復雜,這里一點點解釋。t1用來記錄從Multiwfn的輸出文件里提取的M/Vm這一項,t1=$( ... )代表把$()里的命令的運行結果賦值給t1變量,grep "Estimated density" out.txt代表把Multiwfn輸出文件里含有Estimated density的行提出來,即Estimated density according to mass and volume (M/V):    1.8045 g/cm^3這一行。我們要把這里的1.8045提取出來,因此把grep命令返回的這一行通過管道方式(|符號)傳遞給后面的awk命令。awk是個特別常用的基于命令的文本處理工具,基本用法請自行Google。awk -F : '{print $2}'代表以:分隔符,把傳過來的信息的第2項(被分隔符分隔的第i部分對應$i變量)用print命令輸出出來,即    1.8045 g/cm^3。這部分信息之后又被傳遞給awk '{print $1}',awk默認情況下會以一個或多個空格作為分隔符,因此被輸出的變量$1對應的就是1.8045。
    ? t2=$(grep "Product of sigma^2_tot and nu" out.txt | awk -F \( '{print $2}')這行命令將νσ2tot值賦予t2變量。它先把輸出文件里下面這行取出來
    Product of sigma^2_tot and nu:   0.00006625 a.u.^2 (   26.08649 (kcal/mol)^2)
    然后交由awk輸出以(為分隔符的第2部分,即26.08649。注意由于(符號有特殊含義,所以用-F指定分隔符的時候前面加了個\避免它被轉義。
    ? echo "M/Vm:" $t1 "g/cm^3"這行把相關字符串和t1變量值一起輸出出來給用戶看
    ? 記錄M/Vm和νσ2tot的t1和t變量已經有了,接下來就是將它們與擬合參數進行運算,通過下面這行來實現:
    echo "$t1 $t2" | awk '{ print "Predicted density: " a*$1+b*$2+g " g/cm^3"}' a=0.9183 b=0.0028 g=0.0443
    這里echo "$t1 $t2"將M/Vm和νσ2tot的值輸出為同一行再傳遞給awk,它們就對應awk里$1和$2兩個變量了。awk命令后頭a=0.9183 b=0.0028 g=0.0443是把alpha、beta、gamma三個擬合參數分別命名為a、b、g變量,之后就可以在'{  }'擴住的awk的指令里使用了。awk的指令print "Predicted density: " a*$1+b*$2+g " g/cm^3"用來輸出雙引號擴住的字符串以及預測公式計算的結果。
    實際上不定義a、b、g變量而直接把其數值寫到awk指令里也完全可以,即下面這樣
    echo "$t1 $t2" | awk '{ print "Predicted density: " 0.9183*$1+0.0028*$2+0.0443 " g/cm^3"}'
    ? 最后out.txt就沒用了,因此用rm -f刪除之。用-f時系統不會讓你再確認是否刪除。

    實際上,上面這個腳本倒數第二行也可以寫為
    echo "Predicted density:" `echo "0.9183*$t1+0.0028*$t2+0.0443" | bc -l` "g/cm^3"
    這里` `擴住的echo "0.9183*$t1+0.0028*$t2+0.0443" | bc被當做一整條命令來執行,其結果返回在此處。用$()代替` `來擴住也完全可以。bc是Linux下的一個基于文本的計算器,echo把要做的運算表達式通過|傳遞給了bc,bc就返回了計算結果。bc結合-l選項非常重要,否則考慮的小數位數會可能會很有限導致明顯的誤差。用bc雖然很方便,但是注意在cmder里面沒法用,而前面通過awk實現運算的話既可以在Linux下使用也可以在Windows中通過cmder提供的Bash環境里正常使用。


    2.2 繪制分子表面靜電勢的腳本

    在《使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖》(http://www.shanxitv.org/443)中筆者給出了非常便利的腳本用于將Multiwfn與VMD程序結合繪制效果很好的分子表面靜電勢圖,如今在文獻中已經被廣泛使用,不了解的話先看一下。這一節說一下文中用的ESPiso.bat這個batch腳本的原理。在Multiwfn目錄下的examples\drawESP\ESPiso.bat圖標上點右鍵選“編輯”,可以看到內容如下,是依次執行一系列命令
    Multiwfn 1.fch -ESPrhoiso 0.001 < ESPiso.txt
    move /Y density.cub density1.cub
    move /Y totesp.cub ESP1.cub
    Multiwfn 2.fch -ESPrhoiso 0.001 < ESPiso.txt
    move /Y density.cub density2.cub
    move /Y totesp.cub ESP2.cub
    Multiwfn 3.fch -ESPrhoiso 0.001 < ESPiso.txt
    move /Y density.cub density3.cub
    move /Y totesp.cub ESP3.cub
    Multiwfn 4.fch -ESPrhoiso 0.001 < ESPiso.txt
    move /Y density.cub density4.cub
    move /Y totesp.cub ESP4.cub

    move /Y *.cub "D:\study\VMD193"

    其中Multiwfn 1.fch -ESPrhoiso 0.001 < ESPiso.txt就是讓Multiwfn按照當前目錄下的ESPiso.txt里的指令計算1.fch。我們只需要計算過程產生的cub文件,所以就沒有把輸出信息重定向到某個文件里而是直接顯示在屏幕上。-ESPrhoiso 0.001是特殊的參數,意義在《巨幅降低Multiwfn結合VMD繪制分子表面靜電勢圖耗時的一個關鍵技巧》(http://www.shanxitv.org/602)里做了解釋。Multiwfn整個計算過程中會在當前目錄下產生記錄電子密度格點數據的density.cub和記錄靜電勢格點數據的totesp.cub,腳本中通過move /Y命令將它們分別改名為density1.cub和ESP1.cub,/Y選項要求若當前目錄下有同名文件時總是自動覆蓋而不詢問用戶。

    以類似方式,ESPiso.bat腳本依次對當前目錄下的2.fch、3.fch、4.fch也都這么計算,而如果有的fch文件并不存在,則Multiwfn的運行會自動以報錯方式瞬間結束,也等于跳過這個文件而接著做后面的事。在整個ESPiso.bat的最后,move /Y命令把當前目錄下所有.cub文件都挪到了VMD目錄下,當前假設是D:\study\VMD193。之后,就可以用上述博文里的ESPiso.vmd腳本里定義的iso或iso2命令基于VMD目錄下的電子密度和靜電勢的cub文件繪制分子表面靜電勢著色圖了。

    再來看記錄了Multiwfn要執行的指令的ESPiso.txt里的內容,如下所示,每一行的含義都在//后面進行了注釋(后同)
    5  //進入主功能5(計算格點數據)
    1  //被計算的是電子密度
    2  //中等質量格點
    2  //把格點數據導出為cub文件,對于電子密度默認是density.cub
    0  //返回主菜單
    5  //再次進入主功能5
    12  //被計算的是靜電勢
    1  //低質量格點(比電子密度格點數據的格點質量低一個檔,這是因為計算靜電勢耗時較高,而且分子表面靜電勢變化較為平滑,因此低一個檔也不至于投影上去的靜電勢對應的顏色明顯不夠準確)
    2  //把格點數據導出為cub文件,對于靜電勢默認是totesp.cub

    在前述的《使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖》中也提到,在examples\drawESP\目錄下還有ESPiso_eV.bat和ESPiso_eV.txt,如果用它倆代替分別ESPiso.bat和ESPiso.txt,則靜電勢格點數據就會用eV而非默認的a.u.為單位。這是如何實現的?打開記錄運行指令的ESPiso_eV.txt文件可看到以下內容

    [...略,和ESPiso.txt相同]
    0  //返回主菜單
    r  //重新載入文件(r意為reload)
    totesp.cub  //載入前面剛算完的totesp.cub
    13  //主功能13,格點數據處理功能
    11  //對內存中的格點數據進行數學運算
    5  //乘以一個數
    27.2114  //Hartree到eV的轉換系數
    0  //將內存中的格點數據導出為cub文件
    totesp.cub  //還是叫totesp.cub,因此會把之前的覆蓋


    2.3 對一批激發態做空穴-電子分析

    《使用Multiwfn做空穴-電子分析全面考察電子激發特征》(http://www.shanxitv.org/434)中介紹的空穴-電子分析是如今使用相當普遍的考察電子激發特征的方法,沒看過此文的話建議先看看。在文中給出了一個腳本,可以一次性把所有激發態的基于空穴-電子分析框架定義的指數全都得到,免得自己一次一次手動操作分別算各個激發態、依次手動拷出來數據了。這一節就來解讀一下這個腳本,腳本的內容如下

    #!/bin/bash
    cat << EOF > calcall.txt
    18
    1
    D-pi-A.out
    EOF
    for ((i=1;i<=3;i=i+1))
    do
    cat << EOF >> calcall.txt
    $i
    1
    2
    0
    0
    1
    EOF
    done

    ./Multiwfn D-pi-A.fchk < calcall.txt |tee out.txt
    rm ./calcall.txt ./result.txt -f

    grep "Sr index" ./out.txt |nl >> result.txt;echo >> result.txt
    grep "D index" ./out.txt |nl >> result.txt;echo >> result.txt
    grep "RMSD of hole in" ./out.txt |nl >> result.txt;echo >> result.txt
    grep "RMSD of electron in" ./out.txt |nl >> result.txt;echo >> result.txt
    grep "H index" ./out.txt |nl >> result.txt;echo >> result.txt
    grep "t index" ./out.txt |nl >> result.txt
    echo
    echo "Finished!"

    這個腳本分為兩部分,第一部分是創建一個向Multiwfn傳遞依次對所有激發態做分析的指令文件calcall.txt,第二部分是調用Multiwfn運行并從輸出文件out.txt里提取數據。

    先看第一部分。下面這段內容是創建calcall.txt文件,命令是先進入主功能18(電子激發分析),再進入空穴-電子分析功能(子功能1),然后輸入的D-pi-A.out是Gaussian的TDDFT任務的輸出文件名字,Multiwfn會從里面讀取組態系數信息。
    cat << EOF > calcall.txt
    18
    1
    D-pi-A.out
    EOF

    下面這部分通過for ; do ... done語句進行循環,i變量從1開始每次增加1直到3,因此默認情況相當于依次分析最低3個激發態。
    for ((i=1;i<=3;i=i+1))
    do
    cat << EOF >> calcall.txt  //把這次循環產生的內容追加到之前的calcall.txt的后面
    $i  //被分析的激發態序號,對應i變量當前的值
    1   //可視化和分析空穴、電子、躍遷密度...
    2   //中等質量格點
    0   //現在定量指標已經被輸出了。從后處理菜單返回
    0   //從空穴-電子分析功能返回主功能18
    1   //再次進入空穴-電子分析,之后會再讓你輸入被分析的激發態序號
    EOF
    done

    然后腳本用Multiwfn基于上面的過程產生的指令文件calcall.txt進行計算,信息輸出到out.txt里。之后的rm ./calcall.txt ./result.txt -f把當前已經無用的指令文件刪掉,并且如果上次執行腳本后若產生了result.txt也給刪掉。

    接下來就是從out.txt里提取數據了。Multiwfn對每個激發態計算出一大堆衡量電子激發特征的定量指標,比如輸出D指數的那行內容是這樣的:
    D_x:   0.005  D_y:   0.819  D_z:   0.001    D index:   0.819 Angstrom
    腳本中下面這條語句用來把所有含有D index的行都提取出來,通過管道傳遞給nl命令來給每一行都加上行號(也即相當于在行首標記了相應激發態的序號),并追加到記錄總結果的文件result.txt末尾
    grep "D index" ./out.txt |nl >> result.txt;echo >> result.txt
    上面這行后面的echo >> result.txt用來往result.txt里添加一個空行。;符號使得多行內容可以寫在一行里顯得簡潔。

    上面這個腳本運行后得到的result.txt內容如下,看起來非常清晰
         1  Sr index (integral of Sr function):   0.51896 a.u.
         2  Sr index (integral of Sr function):   0.64906 a.u.
         3  Sr index (integral of Sr function):   0.54538 a.u.

         1  D_x:   0.522  D_y:   0.001  D_z:   0.001    D index:   0.522 Angstrom
         2  D_x:   3.481  D_y:   0.007  D_z:   0.025    D index:   3.481 Angstrom
         3  D_x:   0.574  D_y:   0.001  D_z:   0.001    D index:   0.574 Angstrom

         1  RMSD of hole in X/Y/Z:       1.443   1.160   0.437   Norm:   1.902 Angstrom
         2  RMSD of hole in X/Y/Z:       3.055   0.826   0.740   Norm:   3.251 Angstrom
         3  RMSD of hole in X/Y/Z:       0.984   1.032   0.416   Norm:   1.486 Angstrom
    略...

    空穴-電子分析功能還可以導出比如空穴、電子、躍遷密度、密度差等函數的格點數據。稍微舉一反三改寫腳本和傳給Multiwfn的命令,就可以實現一次性把所有激發態的格點數據都分別導出了。


    2.4 一鍵實現RESP電荷計算的腳本

    在《計算RESP原子電荷的超級懶人腳本(一行命令就算出結果)》(http://www.shanxitv.org/476)里我給了一個僅需一行命令就能計算出做分子動力學很常用的RESP電荷的Bash腳本,是Multiwfn目錄下的examples\RESP\RESP.sh,這里來解讀一下,我就直接把關鍵性的注釋插入到下面給出的原腳本里了。讀者請先看一下那篇博文了解一下腳本的用法,這個腳本會自動調用Gaussian先對結構用較便宜級別進行優化,然后在更好的基組下算單點任務產生波函數文件,最后交由Multiwfn計算RESP電荷。

    以井號開頭的是注釋,是給自己或者用戶看的,包括用處、版本、用法、作者信息
    # A script to calculate RESP charge based on Gaussian and Multiwfn
    # Written by Tian Lu (sobereva@sina.com), 2019-Apr-17
    # You should properly modify content to use proper calculation level and solvent
    # Example:
    # Calculating neutral singlet molecule: RESP.sh H2O.xyz
    # Calculating anionic singlet molecule: RESP.sh ani.pdb -1 1

    #!/bin/bash

    下面把Gaussian的優化和單點任務的關鍵詞賦予到相應變量里,用戶可根據實際情況修改溶劑名、計算級別。pop=MK結合IOp要求Gaussian把擬合點上的靜電勢數值輸出出來
    keyword_opt="# B3LYP/def2SVP em=GD3BJ scrf(solvent=water) opt=loose"
    keyword_SP="# B3LYP/def2TZVP em=GD3BJ scrf(solvent=water) pop=MK IOp(6/33=2,6/42=6)"
    Gaussian=g09  //將執行Gaussian用的命令指定為一個變量

    export inname=$1  //將腳本后面傳入的第一個參數(含有結構信息的文件路徑)設為inname變量
    filename=${inname%.*}  //把inname記錄的文件名的后綴去掉設為filename變量
    suffix=${inname##*.}  //把inname記錄的文件名的擴展名設為suffix變量
    if [ $2 ];then  //如果運行腳本時指定了第2個參數(凈電荷),輸出出來并且設為chg變量
     echo "Net charge = $2"
     chg=$2
    else  //腳本沒有跟著第2個參數的情況,默認電荷為0
     echo "Net charge was not defined. Default to 0"
     chg=0
    fi
    if [ $3 ];then  //如果運行腳本時指定了第3個參數(自旋多重度),輸出出來并且設為multi變量
     echo "Spin multiplicity = $3"
     multi=$3
    else  //腳本沒有跟著第3個參數的情況,默認自旋多重度為1
     echo "Spin multiplicity was not defined. Default to 1"
     multi=1
    fi

    下面的命令讓Multiwfn載入輸入文件,通過主功能100的子功能2把載入的體系坐標導出為非常簡單的xyz文件格式(tmp.xyz)
    Multiwfn $1 > /dev/null << EOF
    100
    2
    2
    tmp.xyz
    0
    q
    EOF

    下面的內容用來創建Gaussian優化任務的輸入文件gau.gjf的坐標之前的部分。可見關鍵詞、凈電荷、自旋多重度都是直接用的前面定義的相應變量。并且計算后保留gau.chk。
    cat << EOF > gau.gjf
    %chk=gau.chk
    $keyword_opt
        //空行
    test
        //空行
    $chg $multi
    EOF

    下面把tmp.xyz的第2行之后的部分(坐標部分)輸出并追加到gau.gjf最后面。NR是awk里的特殊變量,對應被處理的文件當前的行號
    awk '{if (NR>2) print }' tmp.xyz >> gau.gjf

    最后給gau.gjf末尾加兩個空行,這是命令行方式運行Gaussian所要求的
    cat << EOF >> gau.gjf
        //空行
        //空行
    EOF
    rm tmp.xyz  //現在tmp.xyz就沒用了,刪之

    下面運行Gaussian,做的是優化任務。期間大家可以監控這里產生的gau.out判斷優化是否順利,或者判斷如果腳本中途運行失敗是什么原因
    echo Running optimization task via Gaussian...
    $Gaussian < gau.gjf > gau.out

    接下來產生單點任務的輸入文件。標題行、凈電荷、自旋多重度、坐標直接從優化產生的gau.chk中讀取,并且讀取里面收斂的波函數當初猜,有助于減少本次SCF達到收斂的迭代次數。
    cat << EOF > gau.gjf
    %chk=gau.chk
    $keyword_SP geom=allcheck guess=read
        //空行
        //空行
    EOF

    下面運行Gaussian做單點任務
    echo Running single point task via Gaussian...
    $Gaussian < gau.gjf > gau.out

    下面語句的意思是如果Gaussian輸出文件gau.out里能搜到含有Normal termination的行,判斷結果就為真,就顯示Done!,否則就提示出錯并且退出腳本。exit后面跟的數字含義見https://tldp.org/LDP/abs/html/exitcodes.html
    if (grep -Fq "Normal termination" gau.out) ; then
     echo Done!
    else
     echo The task has failed! Exit the script...
     exit 1
    fi

    調用formchk程序將gau.chk轉換為同目錄下的gau.fchk,期間不在屏幕上輸出信息
    echo Running formchk...
    formchk gau.chk> /dev/null

    調用Multiwfn以gau.fchk為輸入文件,通過主功能7(布居分析)里的子功能18計算RESP電荷
    echo Running Multiwfn...
    Multiwfn gau.fchk > /dev/null << EOF
    7
    18
    8  //計算時從Gaussian輸出文件里讀取擬合點的靜電勢信息而不直接計算靜電勢
    1  //開始標準RESP電荷計算
    gau.out  //從gau.out里讀取靜電勢
    y  //導出計算完的RESP電荷以及原子坐標成為gau.chg
    0
    0
    q  //優雅地退出
    EOF

    rm gau.gjf gau.fchk gau.chk gau.out  //刪除之前留下的已無用的文件
    chgname=${1//$suffix/chg}   //把輸入文件的擴展名替換成chg并設為chgname變量
    mv gau.chg $chgname  //把gau.chg改名,使得最終得到的chg文件和輸入文件有相同的名字

    最后輸出任務成功的信息,告訴用戶產生的chg文件的名字。括號前面的\是為了避免此符號被轉義。
    echo Finished! The optimized atomic coordinates with RESP charges \(the last column\) have been exported to $chgname in current folder

    由此例可見,通過寫腳本,可以讓一連串的計算流程完全自動化,變得格外省事。


    3 Multiwfn的批量運行的方法

    使用Multiwfn處理一批文件,一種做法是像前面2.2節那樣,把要運行的一批命令手動寫入批處理文件或Bash腳本里,另一種做法是用循環語句,自動令Multiwfn對滿足條件的一批文件依次進行計算。例如,讓Multiwfn按照genELFcub.txt文件里的指令對當前目錄下所有.wfn文件依次進行計算,并且假定任務會產生ELF.cub文件,要把輸入文件的名字作為前綴加上去,通過下面的腳本就可以實現,分為兩種情況

    ? Bash腳本:建立一個.sh腳本文件,寫入以下內容,然后執行腳本
    #!/bin/bash
    for inf in *.wfn
    do
    echo Running $inf ...
    Multiwfn $inf < genELFcub.txt > /dev/null
    mv ELF.cub ${inf//.wfn}_ELF.cub
    done
    此腳本中inf變量循環當前目錄下每個.wfn文件,交給Multiwfn默默執行,屏幕上提示正在運行哪個文件。${inf//.wfn}_ELF.cub代表去除當前輸入文件的.wfn擴展名后再加上_ELF.cub。比如如果正在算的是Kanan.wfn,則產生的ELF格點文件將被自動改名為Kanan_ELF.cub。

    ? batch腳本(可直接運行于Windows):建立一個文本文件,擴展名改為bat,寫入以下內容并保存,然后雙擊執行,或者在cmd里輸入.bat文件名執行之。里面涉及的語法解釋參看《使用Gaussian時的幾個實用腳本和命令》(http://www.shanxitv.org/258)。
    for /f %%i in ('dir *.wfn /b') do (
    Multiwfn %%i < genELFcub.txt > NUL
    rename ELF.cub %%~ni_ELF.cub
    )


    4 批量運行Multiwfn的實例

    下面看幾個批量運行Multiwfn的一些實例。

    4.1 批量格式轉換

    Multiwfn支持非常豐富的輸入文件格式,而且還能把載入的信息導出成許多文件格式、產生諸多量子化學程序和一些第一性原理程序的輸入文件。利用這個特點,結合自己寫的腳本,就可以輕易實現大批量的文件格式之間自由轉換。這在《一鍵把所有gjf文件轉成xyz文件、把所有Gaussian輸出文件轉成gjf文件的腳本》(http://www.shanxitv.org/530)里給了現成的腳本作為例子,請先閱讀一下此文了解腳本用法。這里我們就看一下其中的將當前目錄下所有Gaussian輸出文件(.out)里最后一幀結構轉化為Gaussian輸入文件(.gjf)的腳本,這對應Multiwfn目錄下的examples\scripts\outgjf.sh。腳本內容和注釋如下

    #!/bin/bash
    icc=0  //初始化累加變量
    nfile=`ls *.out|wc -l`  //ls *.out獲得out文件列表,交由wc -l返回行數,故nfile是out文件的總數
    for inf in *.out  //循環當前目錄下所有out文件
    do
    ((icc++))  //令icc加1。(( ))里可以做這種整數運算。當前的icc是正在被處理的out文件的序號
    echo Converting $inf to ${inf//out/gjf} ... \($icc of $nfile\)  //其中${inf//out/gjf}是把inf變量對于的文件名的擴展名從out改為gjf后的名字。$icc of $nfile告訴用戶正在處理第幾個、總共有多少個
    Multiwfn $inf << EOF > /dev/null
    100  //主功能100
    2   //導出文件和產生輸入文件
    10  //產生Gaussian輸入文件
    ${inf//out/gjf}  //新產生的Gaussian輸入文件名
    0  //返回主菜單
    q  //優雅地退出
    EOF
    done

    運行這個腳本后,可以看到諸如這樣的運行過程提示
    ...略
    Converting B2H6.out to B2H6.gjf ... (5 of 151)
    Converting Benzaldehyde.out to Benzaldehyde.gjf ... (6 of 151)
    Converting Benzene.out to Benzene.gjf ... (7 of 151)
    ...略

    顯然,如果想讓Multiwfn批量在其它兩種文件格式間轉換,只需要修改被考慮的輸入文件對象(即把*.out改成相應的),修改傳遞給Multiwfn的指令(即進入主功能100的子功能2后選擇的選項),并把${inf//out/gjf}里的out和gjf改成相應的即可。


    4.2 提取所有體系的HOMO-LUMO gap并計算平均值

    如《使用Multiwfn觀看分子軌道》(http://www.shanxitv.org/269)所述,Multiwfn在載入fch/molden/mwfn/gms這樣記錄了占據軌道和空軌道的波函數文件后,當進入主功能0時,不僅會蹦出圖形窗口用來觀看體系結構和軌道等值面,在文本窗口還輸出了許多信息,包括HOMO-LUMO gap,如下所示。
     Note: Orbital     5 is HOMO, energy:   -0.291987 a.u.   -7.945379 eV
           Orbital     6 is LUMO, energy:    0.065333 a.u.    1.777798 eV
           HOMO-LUMO gap:    0.357320 a.u.    9.723178 eV    938.144245 kJ/mol

    這一節提供一個Bash腳本,從當前目錄下所有fch文件里提取HOMO-LUMO gap并連同文件名一起輸出到當前目錄下的gap.txt,在最后還給出HOMO-LUMO gap的平均值和最小值。此腳本連同注釋內容如下,原腳本可以在http://www.shanxitv.org/attach/612/getgap.sh下載。

    #!/bin/bash
    rm -f gap.txt  //刪除之前可能殘留的gap.txt
    avggap=0  //用于記錄HOMO-LUMO gap的平均值。先對所有體系累加,最后除以總數
    nfile=0  //記錄當前正處理的文件序號
    for inf in *.fch
    do
    ((nfile++))
    echo Processing $inf "..."
    Multiwfn $inf -silent << EOF > out.txt  //注意這里用了-silent避免蹦出主功能0的圖形界面
    0
    q  //優雅地退出
    EOF
    gapthis=$(grep "HOMO-LUMO gap" out.txt |awk '{print $5}')  //從Multiwfn輸出信息中找含有HOMO-LUMO gap的行,然后提取以空格分隔時的第5個值,即eV為單位的gap
    echo $inf: $gapthis eV >> gap.txt  //向gap.txt里輸出當前體系文件名和gap值
    avggap=$(echo $avggap+$gapthis | bc -l)  //把當前體系的gap累加到avggap變量上
    rm -f out.txt
    if [[ $nfile -eq 1 ]] ; then  //處理的是第1個文件的情況,此時nfile等于1。[[ ]]里可以做條件判斷
      gapmin=$gapthis  //暫把第1個文件的HOMO-LUMO gap當做最小值
      namemin=$inf  //記錄HOMO-LUMO gap最小的體系的文件名
    else
      if [[ $(echo "$gapthis < $gapmin" | bc) -eq 1 ]] ; then  //當前HOMO-LUMO gap小于目前最小值的情況
        gapmin=$gapthis
        namemin=$inf
      fi
    fi
    done

    echo | awk '{printf ("\n%s %10.6f %s\n","Average HOMO-LUMO gap:",t1/t2," eV")}' t1=$avggap t2=$nfile >> gap.txt  //將avggap變量除以已處理的文件總數nfile,得到平均gap,追加到gap.txt
    echo "Done! Result has been outputted to gap.txt in current folder"
    echo "Totally" $nfile "files were processed"
    echo $namemin has minimal gap of $gapmin eV >> gap.txt

    運行這個腳本之后,得到的gap.txt里就是諸如下面的內容
    HCN.fch: 18.809098 eV
    naphthalene.fch: 4.829532 eV
    NH3BF3.fch: 9.803340 eV
    pyrazine.fch: 5.411161 eV
    waterdimer.fch: 11.407130 eV

    Average HOMO-LUMO gap:  10.052052  eV
    naphthalene.fch has minimal gap of 4.829532 eV

    可見,用腳本進行批處理和統計是相當方便的。注意由于此腳本用了bc命令,所以在cmder下面沒法用。

    在上面的腳本中,不能直接用對gapthis和gapmin用邏輯運算符來判斷,因為Bash腳本里只能對整數變量直接做大小的判斷。上面腳本中是這樣寫的:
    if [[ $(echo "$gapthis < $gapmin" | bc) -eq 1 ]] ; then
    這里先用echo "$gapthis < $gapmin" | bc把兩個浮點數比大小的語句傳遞給bc這個數學運算工具,如果判斷為真,返回1,否則返回0。$(echo "$gapthis < $gapmin" | bc)代表bc返回的數值。

    腳本中下面的這行命令使用了格式化輸出,使得輸出信息格式嚴格滿足我們的要求
    echo | awk '{printf ("\n%s %10.6f %s\n","Average HOMO-LUMO gap:",t1/t2," eV")}' t1=$avggap t2=$nfile >> gap.txt
    在awk的指令中,可以通過printf以特定格式進行輸出。"\n%s %10.6f %s\n"中的\n代表換行,%s代表輸出字符串,%10.6f代表輸出的浮點數占10個位置,小數點后保留6位。后面跟著的"Average HOMO-LUMO gap:",t1/t2," eV"這三項就是套用前述的格式要輸出的內容了。awk程序本身是用于處理特定文本文件或者通過管道傳遞過來的文本的,當前我們只是借用awk來做相除運算并格式化輸出變量,所以前面用echo |傳遞一個空信息來走個形式。


    4.3 對AIMD模擬里每一幀計算原子自旋布居、鍵級和自旋密度cub文件并作圖

    4.3.1 前言&準備

    我們既可以對一批不同的體系通過腳本調用Multiwfn做波函數分析,也可以對同一個體系的特定過程(IRC、勢能面掃描、動力學模擬等)的每一幀做波函數分析,將所有幀的分析結果匯總到一起繪圖,就可以非常清晰直觀地考察化學過程中電子結構的變化了。通過腳本實現對掃描和IRC的每一幀批量用Multiwfn分析然后繪圖的過程在《制作動畫分析電子結構特征》(http://www.shanxitv.org/86)和《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》(http://www.shanxitv.org/200)都有詳細說明。

    在本節,筆者專門示例一下怎么靠腳本對ORCA跑的從頭算動力學(AIMD)里每一幀進行波函數分析然后作圖。讀者請閱讀《使用ORCA做從頭算動力學(AIMD)的簡單例子》(http://www.shanxitv.org/576)了解ORCA跑動力學的初步知識。用于本節的例子是下面這個O2與乙炔反應的動力學過程,是北京科音高級量子化學培訓班(http://www.keinsci.com)里的一部分,至于具體模擬細節和本文無關,就不多說了。

    ORCA跑AIMD有個便利之處是可以要求每隔特定步數就把當前幀的波函數導出為gbw文件,可以之后被Multiwfn所利用。在當前模擬的輸入文件中,寫上dump gbw stride 10 filename "wfn",這樣每10步就會在當前目錄下產生一個以wfn為前綴的gbw文件,即wfn000010.gbw、wfn000020.gbw等等,序號是當前的步數,乘以步長就是模擬時間。
    注:如果你是用Gaussian跑分子動力學,沒法直接通過關鍵詞要求程序每隔一定步數就產生一次記錄波函數的chk文件供Multiwfn分析,但可以先用《GauMD2xyz:將Gaussian的從頭算動力學軌跡轉換為xyz軌跡的程序》(http://www.shanxitv.org/126)把動力學輸出文件轉成xyz格式的軌跡,再用《將多幀xyz文件轉化成量子化學輸入文件的工具:xyz2QC》(http://bbs.keinsci.com/thread-12468-1-1.html)中的做法把xyz的每一幀拆分成單點任務gjf文件并要求保留chk,這樣按照《使用Gaussian時的幾個實用腳本和命令》(http://www.shanxitv.org/258)批量運行后就有了每一幀的chk文件。

    如果把Multiwfn的settings.ini里的orca_2mklpath參數設成ORCA目錄下的orca_2mkl可執行文件的路徑,那么Multiwfn直接就可以載入gbw文件并做波函數分析。但假設你當前機子里沒有和產生gbw文件相同版本的ORCA,那Multiwfn就沒法直接打開gbw,而只能在產生gbw文件的機子上用下面的腳本先把當前目錄下的gbw文件都轉化為同名的molden輸入文件(擴展名為.molden.input),之后Multiwfn就可以載入并分析了。如果讀者仔細看了本文前面的內容,這個腳本的原理自然可以輕松理解,所以不再講解。

    #!/bin/bash
    for inf in *.gbw
    do
    echo converting $inf ...
    name=${inf//.gbw}
    orca_2mkl $name -molden
    echo
    done

    現在,當前目錄下我們有了動力學過程中不同時刻的波函數文件,可以開始批量分析了。

    4.3.2 獲得原子的自旋布居隨時間的變化

    自旋密度和自旋布居在《談談自旋密度、自旋布居以及在Multiwfn中的繪制和計算》(http://www.shanxitv.org/353)中有詳細的說明,不了解的話先仔細看看。這里先考察一下動力學過程中幾個關鍵原子的自旋布居隨時間的變化,從而定量了解單電子的分布在反應過程中是如何變化的。

    運行下面的腳本,就可以實現把所有波函數文件里4號氫原子(上面圖中直接會碰上O2的那個氫原子)的原子自旋布居一次性計算并輸出出來
    #!/bin/bash
    for input in `ls *.molden.input`  //循環當前目錄下所有molden輸入文件
    do
    echo Processing $input  //提示用戶正在處理哪個文件
    Multiwfn $input << EOF > out.txt
    7  //布居分析與原子電荷計算
    5  //Mulliken方法
    1  //輸出Mulliken布居分析結果
    n  //不把原子電荷輸出為chg文件
    0  //返回主功能7
    0  //返回主菜單
    q  //優雅地退出
    EOF
    grep "Population of atoms:" -A 100 out.txt|grep "4(H )"|cut -c 41-51 >> 4.txt
    done
    rm -f out.txt

    計算Mulliken電荷時,Multiwfn輸出的原子的電子布居信息是這樣的:
     Population of atoms:
         Atom      Alpha pop.   Beta pop.    Spin pop.     Atomic charge
         1(C )      3.16257      2.79109      0.37147         0.04634
         2(H )      0.60772      0.30555      0.30217         0.08674
         3(C )      2.89559      2.89522      0.00037         0.20919
         4(H )      1.00000      0.00000      1.00000         0.00000
         5(O )      4.23044      3.90425      0.32620        -0.13469
         6(O )      4.10368      4.10389     -0.00021        -0.20757
    上面腳本倒數第二行grep "Population of atoms:" -A 100代表把Multiwfn輸出的out.txt中Population of atoms:這一行及下面足夠多的行(這里是100行,對當前體系足夠多)都提取出來,然后傳遞給grep "4(H )"命令再把其中4(H )這一行提取,最后交給cut -c 41-51命令提取第41~51列的內容,即Spin pop.下面的數值,追加到4.txt末尾。之所先經歷grep "Population of atoms:" -A 100這么一步,是因為當前的out.txt文件里包含4(H )的行并不止上述一處,在Population of atoms:這一行更前面還有其它地方出現過4(H )。

    上面的腳本運行后,產生的4.txt內容如下
      0.00097 
      0.00006 
     -0.00223 
     -0.00601 
    ...略

    類似地,把上面腳本里4(H )改成其它原子、把4.txt改成其它文件名,可以同理把O5和O6的自旋布居也提取出來。最后把H4、O5、O6的自旋布居放在一起通過Origin等程序作圖,恰當設置后即得到下圖

    由上圖的原子自旋布居的變化可見,將近30 fs的時候,H4基本上就變成一個氫原子了,其上面單電子數幾乎精確為1。從上圖還可以看到,起初O2分子上每個O都帶一個alpha單電子,在O2碰上乙炔后,單電子分布變化劇烈,到最后實際上形成了一個HCO自由基和一個CO閉殼層分子,顯然最后不帶單電子(自旋布居為0)的O6就是CO上的氧。

    4.3.3 考察鍵級隨模擬時間的變化

    鍵級是定量考察化學鍵特征的非常方便、常用的指標,在《Multiwfn支持的分析化學鍵的方法一覽》(http://www.shanxitv.org/471)中有詳細介紹。通過下面的腳本,可以考察C1-C3、C3-O6、O5-O6的Mayer鍵級是怎么隨模擬進行而發生變化的。

    #!/bin/bash
    for input in `ls *.molden.input`
    do
    echo Processing $input
    Multiwfn $input << EOF > out.txt
    9   //鍵級分析
    1   //Mayer鍵級
    n   //不導出鍵級矩陣
    0   //返回出菜單
    q   //優雅地退出
    EOF
    grep "1(C )    3(C )" out.txt|cut -c 69-78 >> 1_3.txt
    grep "3(C )    6(O )" out.txt|cut -c 69-78 >> 3_6.txt
    grep "5(O )    6(O )" out.txt|cut -c 69-78 >> 5_6.txt
    done

    Multiwfn計算Mayer鍵級輸出的信息格式是這樣的:
     #    1:    1(C )    2(H ) Alpha:   0.406768 Beta:  0.396838 Total:  0.803606
     #    2:    1(C )    5(O ) Alpha:   0.946191 Beta:  1.328004 Total:  2.274195
     #    3:    2(H )    5(O ) Alpha:   0.027651 Beta:  0.029457 Total:  0.057109
    ...略
    因此靠grep "1(C )    3(C )" out.txt|cut -c 69-78就可以把含有C1-C3鍵級的那一行的第69~78列,即考慮alpha和beta電子總效果的Total:后面那個鍵級值提取出來。

    上面的腳本運行完了之后,1_3.txt、3_6.txt和5_6.txt就分別記錄了C1-C3、C3-O6、O5-O6對各個波函數文件算的鍵級了。將數據一起導入到Origin里并繪圖,就可以得到下圖。可見此圖把乙炔和氧氣反應動力學過程中化學鍵的形成和斷裂的情況一覽無余地清晰展現了出來,顯然非常有意義。

    4.3.4 考察自旋密度等值面圖隨模擬時間的變化

    這一節我們利用Bash腳本、Multiwfn和VMD,把氧氣與乙炔反應過程中自旋密度等值面的變化動畫制作出來,結果如下所示,綠色和藍色分別對應正值和負值等值面。可見此動畫可以把單電子分布在反應過程中的變化非常生動直觀地展現出來,頗有意義。

    為了做上面的動畫,我們先運行下面的腳本,它會一次性對當前目錄下的molden輸入文件進行計算,最終在當前目錄下的cub子目錄下產生一批名為"wfn[序號].cub"的自旋密度的格點數據文件。

    #!/bin/bash
    mkdir cub  //新建cub目錄
    for input in `ls *.molden.input`
    do
    echo Processing $input
    Multiwfn $input << EOF > out.txt
    5  //計算格點數據
    5  //自旋密度
    -10  //修改延展距離
    3  //對當前情況3 Bohr的延展距離已足夠避免等值面被截斷
    2  //中等質量格點
    2  //導出cub文件成為spindensity.cub
    0  //返回主菜單
    q  //優雅地退出程序
    EOF
    mv spindensity.cub cub/${input//molden.input/cub}  //將產生的自旋密度的cub文件的名字改成和輸入文件同名
    done

    將上面腳本產生的cub目錄放到C:\目錄下。啟動VMD,將我寫的這個tcl語言的腳本http://www.shanxitv.org/attach/612/showiso.tcl里的內容粘貼到VMD命令行窗口中運行,VMD就會依次載入各個cub文件,把渲染出的自旋密度為0.02 a.u.的等值面圖保存為VMD目錄下的“[序號].bmp”文件。這個tcl腳本內容不在本文范疇以內,因此不做解釋了,VMD的tcl腳本的編寫我在北京科音分子動力學與GROMACS培訓班(http://www.keinsci.com/workshop/KGMX_content.html)里講得特別系統詳細。

    最后,用免費的視頻處理程序ffmpeg把當前目錄下所有bmp文件合并成gif動畫,分為兩步
    先基于某幀的bmp文件創建調色板文件:ffmpeg -i 00027.bmp -vf palettegen palette.png
    將所有bmp文件合并為video.gif,每秒5幀:ffmpeg -r 5 -i %05d.bmp -i palette.png -lavfi paletteuse video.gif
    這樣得到的video.gif動畫文件即前面所展示的。


    5 總結&其它

    本文非常詳細地介紹了如何通過命令行形式方便地運行Multiwfn,以及如何利用Bash腳本或Windows下的批處理文件調用Multiwfn自動做大批量分析處理。文中不僅深入淺出地詳細介紹了所有相關細節,還給出了類型多樣的實例,由淺入深,涉及到常見腳本編寫的方方面面。只要大家仔細把本文的說明看懂、把例子徹底搞懂,最好有時間時再看看前面筆者給出的Bash shell腳本編寫的入門教程網頁系統地學一下,并且遇到問題時善用google(never baidu),就可以沒有壓力、得心應手地針對自己遇到的實際問題寫出對應的腳本,令自己從繁復的重復性的操作中徹底解放,不再需要總是在Multiwfn的交互式界面中一條條手動敲入命令以及手動提取/處理/統計數據。筆者也非常鼓勵大家把經常通過Multiwfn做的操作步驟較多的分析包裝成腳本,從而能夠以最少限度的操作輕易完成整套分析過程,并且筆者也鼓勵將這樣的腳本在Multiwfn論壇http://bbs.keinsci.com/wfn(中文)或http://www.shanxitv.org/wfnbbs(英文)上向其他用戶分享。

    由于如本文所介紹的,Multiwfn能夠以命令行方式運行,還使得它可以在一些第三方程序中進行調用。大家可以看Multiwfn主頁http://www.shanxitv.org/multiwfn的Resources頁面中的Related tools部分,里面有不少Multiwfn用戶寫的Multiwfn輔助程序或者借助Multiwfn計算的程序都是這么調用的Multiwfn。

    雖然本文的腳本說的都是怎么以命令行方式調用Multiwfn運行以及自動化處理,即便你是其它程序的用戶并且對命令行和腳本編寫缺乏了解的話,應該從本文中也明顯有所獲益。比如如果你是GROMACS分子動力學程序的用戶,仔細看過本文之后應該也已經很明白怎么通過寫腳本從而自動或批量地對一個或多個體系進行建模、模擬和后處理分析了。

    筆者還有很多Multiwfn相關博文里都給出了腳本,使得通過Multiwfn進行分析或繪圖的過程盡可能簡單和自動化,這些博文列于下面,強烈推薦大家看看。只要大家把本文全都讀懂,自然也就能輕松理解這些博文里的腳本運行機制了。
    《制作動畫分析電子結構特征》(http://www.shanxitv.org/86
    《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》(http://www.shanxitv.org/200
    《使用Multiwfn一鍵批量產生各類光譜圖(含演示視頻)》(http://www.shanxitv.org/479
    《RESP2原子電荷的思想以及在Multiwfn中的計算》(http://www.shanxitv.org/531
    《計算適用于OPLS-AA力場做模擬的1.2*CM5原子電荷的懶人腳本》(http://www.shanxitv.org/585
    《使用Multiwfn做自然躍遷軌道(NTO)分析》(http://www.shanxitv.org/377
    《使用Multiwfn繪制躍遷密度矩陣和電荷轉移矩陣考察電子激發特征》(http://www.shanxitv.org/436
    《使用Multiwfn+VMD快速地繪制高質量AIM拓撲分析圖》(http://www.shanxitv.org/445
    《使用Multiwfn和VMD繪制平均局部離子化能(ALIE)著色的分子表面圖》(http://www.shanxitv.org/514
    《將Gaussian等程序收斂的波函數作為ORCA的初猜波函數的方法》(http://www.shanxitv.org/517

    最后,留給大家兩個編寫腳本的小練習,感興趣的讀者可以嘗試做一下,以檢驗有沒有充分掌握本文的內容,需要用到的腳本編寫知識沒有超越前文所介紹的
    (1)首先閱讀《使用Multiwfn超級方便地計算出概念密度泛函理論中定義的各種量》(http://www.shanxitv.org/484)了解怎么通過Multiwfn的主功能22對某個體系一次性計算出化學勢、親電指數等概念密度泛函理論涉及的所有量。假定當前目錄下有一堆中性閉殼層小分子的.xyz文件,請編寫一個腳本,調用Multiwfn的主功能22對各個分子都計算概念密度泛函理論涉及的所有量,并且把所有小分子的電負性數值連同文件名一起輸出到eleneg.txt中,也把所有小分子的親電指數數值連同文件名一起輸出到elephi.txt中,并且在這兩個文件末尾給出所有計算過的分子中最大值、最小值以及對應的文件名。
    (2)閱讀Multiwfn手冊4.4節了解怎么繪制各種平面圖。假設當前目錄下有某條IRC軌跡各個點的fch文件,名為[序號].fch,如00001.fch、00002.fch...。寫一個腳本,調用Multiwfn對這些文件繪制4,9,10三個原子構成的平面的ELF填色圖,保存成和輸入文件同名的png文件。最后,再調用ffmpeg之類的程序將它們合并為gif動畫。

    久久精品国产99久久香蕉