【求助】想用python写一个提取Gaussian偶合常数计算结果的脚本 - 计算化学公社

def lower_triangle_to_full_matrix(input_lists):
# convert a lower triangle lists (diagonal included) to a full matrix
dim = len(input_lists[-1])
ret = [["" for x in range(dim)] for y in range(dim)]
for x in range(dim):
for y in range(dim):
if x<=y:
ret[x][y] = input_lists[y][x]
else:
ret[x][y] = input_lists[x][y]
return ret
def read_lower_triangle_matrix(column_count,input_lines):
import numpy as np
# Read in Gaussian output lower triangle matrix and generate a 2D-python-list
# example input at the end of input (it should not include the title, e.g. spin-spin coupling J (Hz), but include table headers, e.g. 1 2 3 4 5)
# without text to number conversion
discard_lines = [0]+[sum(list(range(column_count+1,0,-5))[:(x+1)]) for x in range(int(column_count/5))] # Table headers to delete
input_lines = [line.split()[1:] for count,line in enumerate(input_lines.splitlines()) if count not in discard_lines] # delete table headers (vertical and horizontal)
lines_to_merge = [[x]+list(range(column_count-5,0,-5))[:int(x/5)] for x in range(column_count)]
lines_to_merge = [list(np.cumsum(x)) for x in lines_to_merge]
lower_triangle = [sum([input_lines[line_count] for line_count in x],[]) for x in lines_to_merge]
return lower_triangle_to_full_matrix(lower_triangle)

复制代码

import re
import sys
import ast
#作者：卡开发发
#日期：2019-08-26
#使用方法：read.py *.log [i,j]
#*.log为Gaussian输出文件，[i,j]分别是行和列指标
#（python指标从0开始，所以比Gaussian对应的指标减去1）
#比如LZ附件的"Total nuclear spin-spin coupling J (Hz)"的[48,30]为-0.245622D+00
#read.py sscc-test.log. [47,29]结果输出为：
#-0.245622D+00
#其余根据需要可进一步进行修改
#get_link用于得到特定link的输出
def get_link(num,text):
#通过输出文件开头的Initial command:...字段获得Gaussian的目录
gdir=re.findall('Initial\scommand:\n\s(.*)l1\.exe',text)[0]
#选择Enter xxx.exe和Leave Link xxx之间的内容
head='\(Enter\s'+gdir+'l'+str(num)+'\.exe\)\n'
tail='\n\sLeave\sLink\s.?'+str(num)+'\sat\s'
link=re.findall(head+'(.*)'+tail,text,re.S)[0]
return link
#get_info用于从link1002输出的部分读取结果
def get_info(text):
#选择density=...到最后的内容
info_=re.findall('density=\s*\d*.\d*\n(.*)',text,re.S)[0]
#选择...(Hz):为标题和数据的分隔
sep='\s\D*\s\(Hz\):'
titles=re.findall(sep,info_,re.S)
infos_=re.split(sep,info_,re.S)[1:]
#数据通过分割后以字典形式保存
infos=[]
for i in range(10):
#数据分割后会被按照顺序排列为一维数组
#[[0,0],[1,0],[1,1],[2,1]...[5,4],[6,0]...[6,4],[7,0]...]
data=re.findall('\s(.\d\.\d*D.\d\d)',infos_[i],re.S)
info={'title':titles[i],'data':data}
infos.append(info)
return infos
#gen_index用于生成与数据一维数组对应的序号
def gen_index(natoms):
#先产生一个下三角阵，满足i>=j
#[0,0]
#[1,0][1,1]
#[2,0][2,1][2,2]
#...
#[9,0][9,1]... ... [9,9]
#... ...
mat=[]
for i in range(natoms):
for j in range(natoms):
if i>=j:
mat.append([i,j])
#按j指标范围的顺序重新排列成一维
#[0,0],[1,0],[1,1],[2,1]...[5,4],[6,0]...[6,4]... #k=0
#[5,5],[6,5],[6,6],[7,5]...[10,0]...[10,9][11,0]...[11,9] #k=1
#...
index=[]
for k in range(natoms//5+1):
for num in mat:
if k*5<=num[1]<(k+1)*5:
index.append(num)
return index
#打开文件，运行参数1用来存放输出文件名
f=open(sys.argv[1],'r')
text=f.read()
#读取link101输出内容来获取原子数目，通过原子数目生成对应的序号
link101=get_link(101,text)
natoms=int(re.findall('NAtoms=\s*(\d*)\sNQM=',link101)[0])
index=gen_index(natoms)
#读取link1002输出内容来获取我们需要的数据
link1002=get_link(1002,text)
infos=get_info(link1002)
#运行参数2用来存放我们指定的序号，以数组形式
index_=ast.literal_eval(sys.argv[2])
#输出序号对应的值
#此处infos[-1]为"Total nuclear spin-spin coupling J (Hz)"，也能换成别的
print(infos[-1]['data'][index.index(index_)])

复制代码