使用Gaussian時的幾個實用腳本和命令

使用Gaussian時的幾個實用腳本和命令

First release: 2014-Nov-5 Last update: 2022-Jun-6

本文提供幾個平時Gaussian計算時常用的簡單小腳本和命令，對于實際研究很有用處，也希望讀者能舉一反三。另外也非常建議讀者看看《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612），里面對腳本編寫有深入淺出的介紹和不少例子，仔細看過后就能輕松理解下文介紹的各種腳本并隨意改編。

1 依次執行多個指令

比如要依次執行g09 < 1.gjf > 1.out、g09 < 2.gjf > 2.out、g09 < 3.gjf > 3.out，可以只輸入一條命令，每條命令間用分號隔開：
g09 < 1.gjf > 1.out;g09 < 2.gjf > 2.out;g09 < 3.gjf > 3.out

也可以寫一個文本文件比如t.sh，里面寫上
g09 < 1.gjf > 1.out
g09 < 2.gjf > 2.out
g09 < 3.gjf > 3.out

然后用chmod +x *給它加上可執行權限，再輸入./t.sh運行即可（后同）。

如果不讓指令依次執行，而是同時執行，把每行命令后面都加上&即可。

2 一次性執行當前目錄下所有Gaussian輸入文件

把以下內容存到比如runall.sh里，然后執行之即可。會把每個.gjf文件用g09執行，生成同文件名但后綴為.out的輸出文件。腳本會提示正在運行哪個文件，運行完之后還會輸出用時多少。
#!/bin/bash
icc=0
nfile=`ls ./*.gjf|wc -l`
for inf in *.gjf
do
((icc++))
echo Running ${inf} ... $$icc of $nfile$
time g09 < ${inf} > ${inf//gjf/out}
echo ${inf} has finished
echo
done

PS：如果是windows環境，需要寫DOS批處理腳本，實現方式見《從高斯windows下的批量執行談DOS批處理文件》（http://www.shanxitv.org/6）。

3 一次性把當前目錄下所有chk文件轉換為fchk文件

把以下內容存到比如chk2fch.sh里，然后執行之即可。會依次調用formchk把當前目錄下每個chk轉換為同文件名的.fchk文件。
#!/bin/bash
for inf in *.chk
do
formchk ${inf}
done

4 一次性執行當前文件夾下所有子目錄里的Gaussian輸入文件

此腳本會依次進入當前目錄下的每個下一級目錄，并把其中的.gjf都依次執行，產生的.out文件的文件名和所在位置都和輸入文件相同。把以下內容存到比如folder_runall.sh里，然后執行它之即可。
#!/bin/bash
shopt -s nullglob
for i in `ls -F |grep /`
do
cd $i
for inf in *.gjf
do
echo Running ${inf} ...
time g09 < ${inf} > ${inf//gjf/out}
echo ${inf} has finished
echo
done
cd ..
done

上面的腳本只能進入下一級的目錄，而再下一級的目錄（或者更深的目錄）里的.gjf文件，以及當前目錄下的.gjf文件都不會執行。如果想把這些.gjf也都執行，應當把以下內容存到比如allfolder_runall.sh里，然后執行它。為清楚起見，每次進入新的目錄時屏幕上都會提示。
#!/bin/bash
shopt -s nullglob
for i in `ls -R |grep :|tr : " "`
do
echo
cd $i
echo "****** Entered" $i folder;echo
for inf in *.gjf
do
echo Running ${inf} ...
time g09 < ${inf} > ${inf//gjf/out}
echo ${inf} has finished;echo
done
cd - > /dev/null
done

5 字符替換

下面這個方法可以批量修改計算的級別和任務類型。

例如將當前目錄下包括任意級子目錄下的.gjf中的M062X替換為B3LYP，執行：sed -i "s/M062X/B3LYP/g" `grep M062X -rl *|grep .gjf`

如果要替換的字符有*、/這樣的符號，需要前面加上\避免被sed轉義。比如6-31G**需要寫成6-31G\*\*，M062X/cc-pVTZ需要寫成M062X\/cc-pVTZ。另外，如果有括號或空格出現，那么grep后面應該用雙引號括住。
例如把MP2/6-311+G(2d,p)都替換成M062X/6-31G*：sed -i "s/MP2\/6-311+G(2d,p)/M062X\/6-31G\*/g" `grep "MP2/6-311+G(2d,p)" -rl *|grep .gjf`
例如把opt freq都替換成NMR關鍵詞：sed -i "s/opt freq/NMR/g" `grep "opt freq" -rl *|grep .gjf`

如果只想替換當前目錄下的.gjf，則把-rl改成-l即可。如果想處理所有文件而不僅限于.gjf，則把|grep .gjf部分刪掉即可。

6 批量刪除文件

用下面的命令可以刪除當前目錄和任意級子目錄下的所有.out文件
find ./ -name "*.out"|xargs rm -f

7 顯示幾何優化收斂情況

這個很簡單。比如C4H8.out是幾何優化輸出文件，執行
grep Converged C4H8.out -A4
就會把當前的優化收斂情況輸出出來。

如果執行
grep Converged C4H8.out -c
就會輸出匹配的次數，也就是相當于顯示優化到了第幾步了。

也可以執行grep -E "out of|Converged" C4H8.out -A4，這樣步數和收斂情況都會輸出。

8 做Counterpoise任務時只用一半的BSSE校正能

筆者在《談談BSSE校正與Gaussian對它的處理》（http://www.shanxitv.org/46）專門談過BSSE校正問題。在J. Chem. Theory Comput., 10, 49 (2014)文中，作者建議對于<= aug-cc-pVTZ檔次基組的計算時如果用Counterpoise方式考慮BSSE，應當只用一半BSSE校正能。但是直接用Gaussian的Counterpoise關鍵詞只會產生不考慮和考慮完整的BSSE校正能的情況。如果把以下腳本放到當前目錄，而且當前目錄里有一批使用了Counterpoise關鍵詞的.log輸出文件，程序就會計算出只使用了一半BSSE校正能的結果。

for filename in `ls -v *.log`
do
echo $filename":"
Edimer=`grep "Counterpoise corrected energy =" $filename | awk -F = '{print $2}'`
EBSSE=`grep "BSSE energy =" $filename | awk -F = '{print $2}'`
Esum=`grep "sum of fragments =" $filename | awk -F = '{print $2}'`
echo "($Edimer-$EBSSE*0.5-$Esum)*627.51" |bc | awk '{printf "%6.2f", $0}'
echo " kcal/mol"
done

結果顯示在屏幕上，如下所示
1.log:
-3.03 kcal/mol
2.log:
-31.46 kcal/mol
3.log:
-20.91 kcal/mol
...略

9 將當前目錄下所有Gaussian輸出文件(out)轉換為輸入文件(gjf)

參見下文，需要利用Multiwfn程序。
一鍵把所有gjf文件轉成xyz文件、把所有Gaussian輸出文件轉成gjf文件的腳本
http://www.shanxitv.org/530（http://bbs.keinsci.com/thread-16161-1-1.html）

10 顯示當前目錄下所有out文件最后一次SCF Done的能量

比如當前目錄下有一大批分子做幾何優化的輸出信息，我們想把最后一次含有SCF Done的行連同輸出文件名顯示出來，由此得到每個體系最終結構下的電子能量，可以用以下腳本

#!/bin/bash
for inf in *.out
do
tac $inf | grep -m 1 "SCF Done" | tr '\n' ' '
echo $inf
done

輸出信息示例：
SCF Done: E(RB3LYP) = -619.115825500     A.U. after    1 cycles C16H16.out
SCF Done: E(RB3LYP) = -696.562886104     A.U. after    1 cycles C18H18.out
SCF Done: E(RB3LYP) = -773.962700591     A.U. after    1 cycles C20H20.out
...略