1.前言

高斯(Gaussian)是一個功能強大的量子化學(xué)綜合軟件包，所有從事計算化學(xué)相關(guān)領(lǐng)域的科研工作者應(yīng)該都使用或者了解過這個軟件。它的輸出文件(.log文件)是一個文本文件，可以利用文本工具打開，其中包含了這個計算任務(wù)的所有計算過程以及計算結(jié)果，而我們關(guān)心的可能只是其中的一小部分信息。我們可以借助簡單的python腳本將我們需要的信息提取出來。本文將介紹如何使用python腳本提取優(yōu)化任務(wù)中的最后一個結(jié)構(gòu)，以及這個結(jié)構(gòu)的電子能，自由能校正等信息。

2.工具

GaussView 6

Gaussian 16

Python 3.7

3.創(chuàng)建計算任務(wù)

我們先利用GaussView構(gòu)建一個高斯計算任務(wù)，筆者用乙苯的優(yōu)化作為例子來說明。手搭的乙苯結(jié)構(gòu)

保存計算任務(wù)，用文本編輯工具比如notepad++來設(shè)定為這個計算任務(wù)分配的計算資源。我這里設(shè)置的是8核，6GB，用的方法/基組是B3LYP/6-31G(d)。用notepad++設(shè)置計算資源

接著執(zhí)行這個計算任務(wù)，計算的結(jié)果會存在.log文件中。這個log文件就是我們后面需要從中讀取信息的文件。

4.從高斯的計算結(jié)果文件中讀取信息

在讀取計算結(jié)果之前，我們首先要確認(rèn)計算任務(wù)是否正常結(jié)束。用notepad++打開log文件，拉到文件最后看到正常結(jié)束的計算任務(wù)log文件末尾

任務(wù)正常結(jié)束。我們打開GaussView可以查看優(yōu)化得到的結(jié)構(gòu)和這個結(jié)構(gòu)的電子能數(shù)據(jù)和熱力學(xué)矯正數(shù)據(jù)，這個任務(wù)總共優(yōu)化了17步，我們用腳本讀取的數(shù)據(jù)就是第17步的這個結(jié)構(gòu)的坐標(biāo)和這個結(jié)構(gòu)對應(yīng)的能量數(shù)據(jù)，即下圖右邊紅框數(shù)據(jù)。用GaussView查看計算結(jié)果

log文件中的信息非常多，因此為了從log中得到這些信息，我們首先需要知道這些信息在log文件的哪些地方。

首先是結(jié)構(gòu)的坐標(biāo)，它們都存在log文件中的Standard orientation這個標(biāo)簽下。因此如果我們要獲得最終優(yōu)化好的結(jié)構(gòu)的三維坐標(biāo)，我們只需要找到log文件中的最后一個Standard orientation標(biāo)簽，并把這個標(biāo)簽下面的內(nèi)容讀取出來。Stardard orientation標(biāo)簽下的內(nèi)容

實現(xiàn)邏輯非常簡單，首先用腳本找到所有的Standard orientation行的行號，在這個文件中有18個，我們選擇最后最后一個行號，我們假設(shè)是i，根據(jù)Standard orientation標(biāo)簽的結(jié)構(gòu)特點，坐標(biāo)的起始行數(shù)就是i+5，如果我們知道這個分子有n個原子，那么坐標(biāo)的結(jié)尾行號就是i+5+n。而log文件中有一些行中在NAtoms后面寫明了這個分子有多少個原子。從NAtoms中得知這個分子有多少原子

因此，很容易的我們就能寫出讀取坐標(biāo)的代碼

atomicnum2elemnt = {1:'H',6:'C'}

log_file = 'The path of your .log file' # e.g. log_file = 'C:/Users/silly/Desktop/ethylbenzene.log'

with open(log_file,'r') as fr:

lines = fr.readlines()

coord_start_index_list = []

for i,line in enumerate(lines):

if 'NAtoms=' in line:

atom_num = eval(line.split()[1])

if 'Standard orientation' in line:

coord_start_index_list.append(i+5)

coord_string = lines[coord_start_index_list[-1]:coord_start_index_list[-1]+atom_num]

coord = [[eval(item.strip().split()[3]),eval(item.strip().split()[4]), eval(item.strip().split()[5])] for item in coord_string]

atom_type = [atomicnum2elemnt[eval(item.split()[1])] for item in coord_string]

Coord = ''

for tmp_coord,tmp_atom in zip(coord,atom_type):

Coord += '%2s %15f %15f %15f\n'%(tmp_atom,tmp_coord[0],tmp_coord[1],tmp_coord[2])

通過以上代碼我們將三維坐標(biāo)信息存在了coord變量里，元素信息存在了atom_type變量里，然后重新寫了一個Coord字符串保存可以寫進文件里的分子坐標(biāo)。成功將坐標(biāo)信息提出并存入字符串變量Coord

接下來我們用類似的方法讀取我們需要的能量信息。其實操作和上面提取坐標(biāo)的操作大同小異，我們只需要知道這些能量信息在哪些行，然后定位到這些行然后讀取這些信息即可。由于每一個優(yōu)化的結(jié)構(gòu)都會給出一個電子能，因此對于電子能我們只取最后一個結(jié)構(gòu)的電子能，而其他熱力學(xué)校正的數(shù)據(jù)只會出現(xiàn)一次，因此直接讀就完事了。記錄能量信息的行

log_file = 'The path of your .log file' # e.g. log_file = 'C:/Users/silly/Desktop/ethylbenzene.log'

with open(log_file,'r') as fr:

lines = fr.readlines()

EE_list = []

TCH = 0

TCG = 0

for line in lines:

if 'SCF Done' in line:

EE_list.append(eval(line.split()[4]))

elif 'Thermal correction to Enthalpy' in line:

TCH = eval(line.strip().split()[4])

elif 'Thermal correction to Gibbs Free Energy' in line:

TCG = eval(line.strip().split()[6])

EE = EE_list[-1]

EE中保存了最后一個結(jié)構(gòu)的電子能，TCH中保存了焓校正，TCG中保存了自由能校正。

5.結(jié)束語

這個腳本是筆者一年多前剛進入現(xiàn)在所在的課題組的時候?qū)懙?#xff0c;是一個入門級的python腳本，當(dāng)然，目前這個版本的腳本比一年前的版本要高明不少因此也簡潔不少。只要會python的文件操作并且了解log文件的構(gòu)成就能出一堆類似的腳本。以上展示的只是讀取這些信息的關(guān)鍵步驟，這些讀出來的信息后續(xù)如何使用讀者們可以自己發(fā)揮。

如果這個腳本對你有幫助，請麻煩你給這篇文章點個贊。如果有任何問題可以在評論區(qū)留言，也可以私信我。

總結(jié)

以上是生活随笔為你收集整理的如何写python脚本抓取数据并计算_【小工具】利用Python脚本从Gaussian计算结果中提取信息...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。