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

happyknighthawk

各位老师好：
我在做Gaussian偶合常数计算的时候，发现手动提取偶合常数计算值非常耗时，尤其当涉及到要提取多个原子对的计算结果、以及当分子存在多构象（需要做构象平均时），手动提取更是耗时，于是萌生了想要写一个脚本来提取计算结果的想法。
断断续续自学python有一小阵了，但是当要自己写脚本来解决实际问题的时候，还是有点发懵。我的初步想法是：把偶合常数计算结果对应的所有行【如Total nuclear spin-spin coupling J (Hz): 下面的行】从文件中提取出来，经过一系列操作，最后用numpy生成一个二维数组。当给出一对原子编号（a,b）时，就可以通过array[a][b]索引出对应的偶合常数。
目前，“提取偶合常数计算结果对应的行”问题不大，但是“转换二维数组”捣鼓了好一阵也搞不定。
现在，想请教各位老师：（1）首先，我是不是把这个问题弄的太复杂了？（2）其次，能不能请各位老师针对这个需求给我一些关键的思路？谢谢！
（PS：附件中附了一个我之前用B972/pcJ-1做偶合常数计算的输出文件，由于关键词里面写了readatoms，仅指定计算了个别原子，所以计算结果里面大部分是0。但形式跟计算所有原子时并没有差别）
非常非常感谢！

让你变成回忆

可以参考一下Multiwfn中读取Overlap矩阵的子程序，具体可以参考Multiwfn源代码中，util.f90文件中的“subroutine readmatgau”子程序。

chrinide

用pandas的concat函数做数据拼接：把原子编号作为Dataframes的index和colomn names，做横向拼接，你就会得到一个64*64的Dataframe，就搞定了

liyuanhe211

1. 
   
2. def lower_triangle_to_full_matrix(input_lists):
   
3.     # convert a lower triangle lists (diagonal included) to a full matrix
   
4.     dim = len(input_lists[-1])
   
5.     ret = [["" for x in range(dim)] for y in range(dim)]
   
6.     for x in range(dim):
   
7.         for y in range(dim):
   
8.             if x<=y:
   
9.                 ret[x][y] = input_lists[y][x]
   
10.             else:
    
11.                 ret[x][y] = input_lists[x][y]
    
12.     return ret
    
13. 
    
14. 
    
15. 
    
16. 
    
17. def read_lower_triangle_matrix(column_count,input_lines):
    
18.     import numpy as np
    
19. 
    
20.     # Read in Gaussian output lower triangle matrix and generate a 2D-python-list
    
21.     # 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)
    
22.     # without text to number conversion
    
23. 
    
24.     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
    
25.     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)
    
26.     lines_to_merge = [[x]+list(range(column_count-5,0,-5))[:int(x/5)] for x in range(column_count)]
    
27.     lines_to_merge = [list(np.cumsum(x)) for x in lines_to_merge]
    
28.     lower_triangle = [sum([input_lines[line_count] for line_count in x],[]) for x in lines_to_merge]
    
29.     return lower_triangle_to_full_matrix(lower_triangle)

复制代码

很早之前写的，仅供参考

happyknighthawk

让你变成回忆 发表于 2019-8-26 08:23
可以参考一下Multiwfn中读取Overlap矩阵的子程序，具体可以参考Multiwfn源代码中，util.f90文件中的“subro ...

非常感谢您的指点，我记得Multiwfn是用Fortran写的，目前的确还没有相关的编程基础，等python再学熟悉一点，后续我会再学Fortran的。依然感谢您！

happyknighthawk

chrinide 发表于 2019-8-26 10:56
用pandas的concat函数做数据拼接：把原子编号作为Dataframes的index和colomn names，做横向拼接，你就会得 ...

非常感谢您，我马上去了解一下pandas的concat函数，太感谢了

happyknighthawk

liyuanhe211 发表于 2019-8-26 11:41
很早之前写的，仅供参考

非常感谢李老师，我好好研究一下您提供的代码，实在太感谢您了！

卡开发发

提供个土方法：你可以按照数据的顺序去生成index=[atom_i,atom_j]（i>=j），然后数据data直接用正则表达式抓，然后给定i、j就data(index.index([i,j]))。

1. import re
   
2. import sys
   
3. import ast
   
4. #作者：卡开发发
   
5. #日期：2019-08-26
   
6. #使用方法：read.py *.log [i,j]
   
7. #*.log为Gaussian输出文件，[i,j]分别是行和列指标
   
8. #（python指标从0开始，所以比Gaussian对应的指标减去1）
   
9. #比如LZ附件的"Total nuclear spin-spin coupling J (Hz)"的[48,30]为-0.245622D+00
   
10. #read.py sscc-test.log. [47,29]结果输出为：
    
11. #-0.245622D+00
    
12. #其余根据需要可进一步进行修改
    
13. 
    
14. #get_link用于得到特定link的输出
    
15. def get_link(num,text):
    
16.         #通过输出文件开头的Initial command:...字段获得Gaussian的目录
    
17.         gdir=re.findall('Initial\scommand:\n\s(.*)l1\.exe',text)[0]
    
18.         #选择Enter xxx.exe和Leave Link xxx之间的内容
    
19.         head='\(Enter\s'+gdir+'l'+str(num)+'\.exe\)\n'
    
20.         tail='\n\sLeave\sLink\s.?'+str(num)+'\sat\s'
    
21.         link=re.findall(head+'(.*)'+tail,text,re.S)[0]
    
22.         return link
    
23. 
    
24. #get_info用于从link1002输出的部分读取结果
    
25. def get_info(text):
    
26.         #选择density=...到最后的内容
    
27.         info_=re.findall('density=\s*\d*.\d*\n(.*)',text,re.S)[0]
    
28.         #选择...(Hz):为标题和数据的分隔
    
29.         sep='\s\D*\s\(Hz\):'
    
30.         titles=re.findall(sep,info_,re.S)
    
31.         infos_=re.split(sep,info_,re.S)[1:]
    
32.         #数据通过分割后以字典形式保存
    
33.         infos=[]
    
34.         for i in range(10):
    
35.                 #数据分割后会被按照顺序排列为一维数组
    
36.                 #[[0,0],[1,0],[1,1],[2,1]...[5,4],[6,0]...[6,4],[7,0]...]
    
37.                 data=re.findall('\s(.\d\.\d*D.\d\d)',infos_[i],re.S)
    
38.                 info={'title':titles[i],'data':data}
    
39.                 infos.append(info)
    
40.         return infos
    
41. 
    
42. #gen_index用于生成与数据一维数组对应的序号
    
43. def gen_index(natoms):
    
44.         #先产生一个下三角阵，满足i>=j
    
45.         #[0,0]
    
46.         #[1,0][1,1]
    
47.         #[2,0][2,1][2,2]
    
48.         #...
    
49.         #[9,0][9,1]... ... [9,9]
    
50.         #... ...
    
51.         mat=[]
    
52.         for i in range(natoms):
    
53.                 for j in range(natoms):
    
54.                         if i>=j:
    
55.                                 mat.append([i,j])
    
56.         #按j指标范围的顺序重新排列成一维
    
57.         #[0,0],[1,0],[1,1],[2,1]...[5,4],[6,0]...[6,4]... #k=0
    
58.         #[5,5],[6,5],[6,6],[7,5]...[10,0]...[10,9][11,0]...[11,9] #k=1
    
59.         #...
    
60.         index=[]
    
61.         for k in range(natoms//5+1):
    
62.                 for num in mat:
    
63.                         if k*5<=num[1]<(k+1)*5:
    
64.                                 index.append(num)
    
65.         return index
    
66. 
    
67. #打开文件，运行参数1用来存放输出文件名
    
68. f=open(sys.argv[1],'r')
    
69. text=f.read()
    
70. #读取link101输出内容来获取原子数目，通过原子数目生成对应的序号
    
71. link101=get_link(101,text)
    
72. natoms=int(re.findall('NAtoms=\s*(\d*)\sNQM=',link101)[0])
    
73. index=gen_index(natoms)
    
74. #读取link1002输出内容来获取我们需要的数据
    
75. link1002=get_link(1002,text)
    
76. infos=get_info(link1002)
    
77. #运行参数2用来存放我们指定的序号，以数组形式
    
78. index_=ast.literal_eval(sys.argv[2])
    
79. #输出序号对应的值
    
80. #此处infos[-1]为"Total nuclear spin-spin coupling J (Hz)"，也能换成别的
    
81. print(infos[-1]['data'][index.index(index_)])

复制代码

happyknighthawk

卡开发发 发表于 2019-8-26 17:07
提供个土方法：你可以按照数据的顺序去生成index=[atom_i,atom_j]（i>=j），然后数据data直接用正则表达式 ...

太感谢卡开发发老师了，之前实在没敢奢求有大神能帮我写出直接能用的脚本，感谢您的无私付出，我明天好好学习下这个脚本！
这两天我一直在捣鼓这个事情，也希望后面能把自己手头这个写完，哪怕代码丑一点
再次感谢您！

卡开发发

happyknighthawk 发表于 2019-8-27 01:37
太感谢卡开发发老师了，之前实在没敢奢求有大神能帮我写出直接能用的脚本，感谢您的无私付出，我明天好好 ...

没关系啊，我也是初学。大家一起讨论，程序一起研究，这样知识学到手才会有点带感