python如何输出计算结果_【小工具】利用Python脚本从Gaussian计算结果中提取信息...

1.前言

高斯(Gaussian)是一个功能强大的量子化学综合软件包，所有从事计算化学相关领域的科研工作者应该都使用或者了解过这个软件。它的输出文件(.log文件)是一个文本文件，可以利用文本工具打开，其中包含了这个计算任务的所有计算过程以及计算结果，而我们关心的可能只是其中的一小部分信息。我们可以借助简单的python脚本将我们需要的信息提取出来。本文将介绍如何使用python脚本提取优化任务中的最后一个结构，以及这个结构的电子能，自由能校正等信息。

2.工具

GaussView 6

Gaussian 16

Python 3.7

3.创建计算任务

我们先利用GaussView构建一个高斯计算任务，笔者用乙苯的优化作为例子来说明。手搭的乙苯结构

保存计算任务，用文本编辑工具比如notepad++来设定为这个计算任务分配的计算资源。我这里设置的是8核，6GB，用的方法/基组是B3LYP/6-31G(d)。用notepad++设置计算资源

接着执行这个计算任务，计算的结果会存在.log文件中。这个log文件就是我们后面需要从中读取信息的文件。

4.从高斯的计算结果文件中读取信息

在读取计算结果之前，我们首先要确认计算任务是否正常结束。用notepad++打开log文件，拉到文件最后看到正常结束的计算任务log文件末尾

任务正常结束。我们打开GaussView可以查看优化得到的结构和这个结构的电子能数据和热力学矫正数据，这个任务总共优化了17步，我们用脚本读取的数据就是第17步的这个结构的坐标和这个结构对应的能量数据，即下图右边红框数据。用GaussView查看计算结果

log文件中的信息非常多，因此为了从log中得到这些信息，我们首先需要知道这些信息在log文件的哪些地方。

首先是结构的坐标，它们都存在log文件中的Standard orientation这个标签下。因此如果我们要获得最终优化好的结构的三维坐标，我们只需要找到log文件中的最后一个Standard orientation标签，并把这个标签下面的内容读取出来。Stardard orientation标签下的内容

实现逻辑非常简单，首先用脚本找到所有的Standard orientation行的行号，在这个文件中有18个，我们选择最后最后一个行号，我们假设是i，根据Standard orientation标签的结构特点，坐标的起始行数就是i+5，如果我们知道这个分子有n个原子，那么坐标的结尾行号就是i+5+n。而log文件中有一些行中在NAtoms后面写明了这个分子有多少个原子。从NAtoms中得知这个分子有多少原子

因此，很容易的我们就能写出读取坐标的代码

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])

通过以上代码我们将三维坐标信息存在了coord变量里，元素信息存在了atom_type变量里，然后重新写了一个Coord字符串保存可以写进文件里的分子坐标。成功将坐标信息提出并存入字符串变量Coord

接下来我们用类似的方法读取我们需要的能量信息。其实操作和上面提取坐标的操作大同小异，我们只需要知道这些能量信息在哪些行，然后定位到这些行然后读取这些信息即可。由于每一个优化的结构都会给出一个电子能，因此对于电子能我们只取最后一个结构的电子能，而其他热力学校正的数据只会出现一次，因此直接读就完事了。记录能量信息的行

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中保存了最后一个结构的电子能，TCH中保存了焓校正，TCG中保存了自由能校正。

5.结束语

这个脚本是笔者一年多前刚进入现在所在的课题组的时候写的，是一个入门级的python脚本，当然，目前这个版本的脚本比一年前的版本要高明不少因此也简洁不少。只要会python的文件操作并且了解log文件的构成就能出一堆类似的脚本。以上展示的只是读取这些信息的关键步骤，这些读出来的信息后续如何使用读者们可以自己发挥。

如果这个脚本对你有帮助，请麻烦你给这篇文章点个赞。如果有任何问题可以在评论区留言，也可以私信我。