Gaussian二维扫描输出文件的后处理脚本

Accelerator · 发表于 Post on 2021-4-30 11:46:40

本帖最后由 Accelerator 于 2022-3-29 09:17 编辑

使用Gaussian进行二维柔性扫描是非常恼人的工作，有如下几点原因：
1. 经常会在花费大量时间之后，扫到一些点时发生构型不收敛的问题而报错退出，此时用GaussView无法直观显示已经扫过的势能面的能量（只能显示折线），极为令人沮丧。
2. 即使整个二维扫描过程正常结束，用GaussView能显示势能面的曲面，有时也无法正确识别曲面上点的构型，会发现点住曲面上某个点后显示的构型与两个坐标数值对应的不符。
3. 经常只关心势能面的一部分区域，全面扫描会带来大量机时浪费。（我昨天在http://bbs.keinsci.com/thread-22864-1-1.html提到了一种思路，先用低耗时方法进行完整扫描再提取感兴趣的构型，并且大家对使用UFF进行扫描的注意事项进行了解答。）

为了解决上述问题，我写了一个Python脚本，分享出来。它可以读取Gaussian的二维扫描输出文件，使用方法如下：
1. 执行python3 anaScan2D.py log文件位置两个扫描坐标。例如：

python3 anaScan2D.py test.log 'R(1,2)' 'R(1,3)'

复制代码

不一会，就会输出当前已经成功扫过的每个点的两个扫描坐标数值以及能量。这样，即使扫描没有正确结束，也能利用已有的数据对势能面作图：

This script is to analyze a 2D scan which is not correctly completed. It will extract all the completed points so that you could still draw a 2D-PES figure.
R1= 1.3943 , R2= 1.3914 , E= 0.161499
R1= 1.3956 , R2= 1.3912 , E= 0.140213
R1= 1.3968 , R2= 1.3908 , E= 0.136637
R1= 1.3979 , R2= 1.3906 , E= 0.143130

复制代码

2. 接下来会询问是否需要提取其中部分点的坐标，写y，按照提示输入要提取出来的两个坐标的范围，新生成的gjf文件的关键字，使用的核心数和内存大小（GB），电荷和自旋多重度，就会自动生成一系列.gjf文件。文件尾部的信息需要提前写在当前目录的end.txt中，程序会读取其内容放在各gjf文件末尾。

Now are you going to generate gjf files for a range of points? y/n
y
Range for coord1: e.g. 1.0-2.0
1.0-1.5
Range for coord2: e.g. 1.0-2.0
1.0-1.5
Your kwds:
opt=modredundant freq b3lyp em=gd3bj 6-31G(d)
（以下略去）

复制代码

如果我在end.txt中写B 23 12 F（换行）B 22 23 F，就会生成一系列类似下面的文件，名称为test_XXX_YYY.gjf，XXX和YYY分别为该点对应的两个扫描坐标数值。

%mem=24gb
%nprocs=16
# opt=modredundant freq b3lyp em=gd3bj 6-31G(d)
TC
0 1
（略去）
B 23 12 F
B 22 23 F

复制代码

anaScan2D.py (2.8 KB, 下载次数 Times of downloads: 23)

ginlpein · 发表于 Post on 2022-3-24 18:24:47

本帖最后由 ginlpein 于 2022-3-24 18:27 编辑

您好，请问您这个脚本该怎么能够适用于对二面角扫描结果的提取？
计算使用的是PM7半经验计算对两个二面角进行二维扫描。按照您写的内容运行，扫描坐标不论以何种形式输入，似乎都没能提取出各点能量，直接跳到询问是否生成gjf那步了，请问是我哪里操作错误吗？或者针对二面角需要哪里作出对应修改？
是否生成gjf文件那边选择y最后也是什么都没有生成出来。

D:\Relax_scan_tool\Python>python3 anaScan2D.py scan.log 'D(37,6,14,4)' 'D(38,39,81,143)'
This script is to analyze a 2D scan which is not correctly completed. It will extract all the completed points so that you could still draw a 2D-PES figure.
Now are you going to generate gjf files for a range of points? y/n
N
D:\Relax_scan_tool\Python>python3 anaScan2D.py scan.log 'D61' 'D346'
This script is to analyze a 2D scan which is not correctly completed. It will extract all the completed points so that you could still draw a 2D-PES figure.
Now are you going to generate gjf files for a range of points? y/n
N

复制代码

我的计算体系在200~500原子的级别，由于比较结构大，导致二维柔性扫描结果文件普遍大于1G，经常超过2G。超大的结果文件导致不论计算成功完成与否，gaussview都无法正常读取结果，Result-scan中无法正确显示或输出结果。尝试cclib等程序也不太支持柔性扫描，所以一直在寻找可以用于解析计算结果的方法，希望楼主赐教。

Accelerator · 发表于 Post on 2022-3-24 19:45:35

ginlpein 发表于 2022-3-24 18:24
您好，请问您这个脚本该怎么能够适用于对二面角扫描结果的提取？
计算使用的是PM7半经验计算对两个二面角 ...

脚本的原理是对log文件每一行如果找到“D(37,6,14,4)”等字样就记录接下来出现的最近一次能量，所以我怀疑D(37,6,14,4)，D(38,39,81,143)等是否在log文件内容里有出现？
或者你生成一个体积比较小的测试文件（比如每个坐标只扫描2步），传上来给我调试。

ginlpein · 发表于 Post on 2022-3-25 15:56:45

本帖最后由 ginlpein 于 2022-3-26 23:26 编辑

您好，后面我尝试做了一个3*3的双二面角扫描，使用您的脚本依旧无法实现解析。

D:\Relax_scan_tool\Python>python3 anaScan2D.py scan_test.log 'D(37,6,14,4)' 'D(38,39,81,143)'
This script is to analyze a 2D scan which is not correctly completed. It will extract all the completed points so that you could still draw a 2D-PES figure.
Now are you going to generate gjf files for a range of points? y/n
n

复制代码

而在log文件中搜索扫描坐标，确实可以搜索到结果，扫描坐标出现在如下位置及情况：

! D61 D(37,6,14,4) -53.1497 Scan !
! D62 D(37,6,14,9) 127.3039 estimate D2E/DX2 !
! D63 D(41,6,14,4) 129.5928 estimate D2E/DX2 !
! D64 D(41,6,14,9) -49.9536 estimate D2E/DX2 !
! D65 D(14,6,37,38) -173.96 estimate D2E/DX2 !
! D66 D(14,6,37,373) 1.2749 estimate D2E/DX2 !
! D67 D(41,6,37,38) 3.3306 estimate D2E/DX2 !
! D68 D(41,6,37,373) 178.5655 estimate D2E/DX2 !
! D69 D(14,6,41,40) 175.002 estimate D2E/DX2 !
! D70 D(14,6,41,385) -6.0362 estimate D2E/DX2 !
! D71 D(37,6,41,40) -2.3562 estimate D2E/DX2 !
! D72 D(37,6,41,385) 176.6056 estimate D2E/DX2 !
! D73 D(67,7,16,5) -53.1497 estimate D2E/DX2 !

复制代码

! D61 D(37,6,14,4) -53.1497 -DE/DX = 0.0003 !
! D62 D(37,6,14,9) 125.2525 -DE/DX = 0.0001 !
! D63 D(41,6,14,4) 132.6627 -DE/DX = 0.0002 !
! D64 D(41,6,14,9) -48.935 -DE/DX = 0.0 !
! D65 D(14,6,37,38) -168.2441 -DE/DX = -0.0001 !
! D66 D(14,6,37,373) 11.9407 -DE/DX = -0.0001 !
! D67 D(41,6,37,38) 5.8314 -DE/DX = 0.0 !
! D68 D(41,6,37,373) -173.9838 -DE/DX = 0.0 !

复制代码

而扫描坐标在log中出现的次数恰好是扫描点数n的n+1次。但是确实无法读取出个点的能量结果。
附件中是我测试用做的3*3扫描的log文件。我的结果文件即使3*3的扫描规模也达到了80多M的级别，似乎不可避免的就是会很大。

Accelerator · 发表于 Post on 2022-3-25 16:27:39

ginlpein 发表于 2022-3-25 15:56
您好，后面我尝试做了一个3*3的双二面角扫描，使用您的脚本依旧无法实现解析。

而在log文件中搜索扫描坐 ...

这是由于二面角为负值，原始代码中用-DE/-DX中的负号作为分隔符对行进行分割，遇到负数就没能正确识别。
只需要将第44行改为：
r1s.append(float(line.split(r1)[1].split('-DE')[0].strip()))
48行也同上修改即可。

ginlpein · 发表于 Post on 2022-3-25 20:57:39

本帖最后由 ginlpein 于 2022-3-25 20:58 编辑

Accelerator 发表于 2022-3-25 16:27
这是由于二面角为负值，原始代码中用-DE/-DX中的负号作为分隔符对行进行分割，遇到负数就没能正确识别。
...

感谢您的耐心回复！
不知道是不是因为我改的不对，我把44行与48行的'-'改成'-DE' 后对附件中的scan_test进行测试依旧无法生效，不知道问题出在哪？

import os, sys
print('This script is to analyze a 2D scan which is not correctly completed. \
It will extract all the completed points so that you could still draw a 2D-PES figure.')
_, path, r1, r2 = sys.argv
r1s = []
r2s = []
es = []
coords = []
lineNum = 0
from common import readLog
elementTable = ['pass','H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si','P','S','Cl','Ar','K','Ca',\
'Sc','Ti','V','Cr','Mn','Fe','Co','Ni','Cu','Zm','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr','Nb',\
'Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I','Xe','Cs','Ba',\
'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf',\
'Ta','W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn']
import re
with open(path) as f:
etemp = 0
isStation = False
flagGeo = 0
thisCoord = []
previousCoord = []
while True:
line=f.readline()
lineNum += 1
if line == "":
break
if "SCF Done" in line:
etemp = float(line.split('=')[1].split("A")[0].strip())
continue
if "point found" in line:
isStation = True
continue
if 'Standard orientation' in line:
flagGeo = 1
continue
if flagGeo == 1 and (re.match("\s*[0-9]+\s*[0-9]+\s*[0-9]+\s*\-*[0-9]+\.",line) != None):
lsplitted = line.split()
element = elementTable[int(lsplitted[1])]
thisCoord.append([element] + lsplitted[-3:])
if (r1 in line) and isStation:
r1s.append(float(line.split(r1)[1].split('-DE')[0].strip()))
es.append(etemp)
continue
if (r2 in line) and isStation:
r2s.append(float(line.split(r2)[1].split('-DE')[0].strip()))
continue
if "Rotational" in line:
flagGeo = 0
if isStation:
coords.append(previousCoord)
previousCoord = thisCoord[:]
thisCoord = []
isStation = False
i=0
for e in es:
print('R1= {:.4f} , R2= {:.4f} , E= {:.6f} '.format(r1s[i],r2s[i],e))
i+=1
def withinRange(value,rangeList):
if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
return True
from common import writeGjf
#def writeGjf(path,coords,kwd,core,mem,charge,mult,end,chk='')
print('Now are you going to generate gjf files for a range of points? y/n')
if(input() != 'y'):
exit()
print('Range for coord1: e.g. 1.0-2.0')
rangeR1 = input().split('-')
print('Range for coord2: e.g. 1.0-2.0')
rangeR2 = input().split('-')
print('Your kwds:')
kwd = input()
print('Your cores:')
core = input()
print('Your memory:')
mem = input()
print('Your charge:')
charge = input()
print('Your multiplicity:')
mult = input()
print('Now reading end.txt for the end part of your gjf file.')
end = ''
try:
with open('end.txt') as f:
end = f.read()
except:
print('No end.txt found.')
itemNum = 0
for i in coords:
r1 = r1s[itemNum]
r2 = r2s[itemNum]
if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
writeGjf(path.replace('.log',f'_{r1}_{r2}.gjf'),i,kwd,core,mem,charge,mult,end)
itemNum += 1

复制代码

D:\Relax_scan_tool\Python>python3 anaScan2D_2.py scan_test.log 'D(37,6,14,4)' 'D(38,39,81,143)'
This script is to analyze a 2D scan which is not correctly completed. It will extract all the completed points so that you could still draw a 2D-PES figure.
Now are you going to generate gjf files for a range of points? y/n
N

复制代码

Accelerator · 发表于 Post on 2022-3-25 21:23:40

ginlpein 发表于 2022-3-25 20:57
感谢您的耐心回复！
不知道是不是因为我改的不对，我把44行与48行的'-'改成'-DE' 后对附件中的scan_test ...

破案了，是由于你是在Windows上运行的，运行命令不应该在扫描坐标前后写引号，即应该：D:\Relax_scan_tool\Python>python3 anaScan2D_2.py scan_test.log D(37,6,14,4) D(38,39,81,143)这样写。不然引号会被识别为坐标的一部分。
在Linux上，是要带着引号的。

另外你这个分子比较大，log里没有输出标准朝向，如果要读取几何坐标，还应该将第36行左右的Standard Orientation改成Input Orientation

ginlpein · 发表于 Post on 2022-3-25 22:38:25

本帖最后由 ginlpein 于 2022-3-25 23:11 编辑

Accelerator 发表于 2022-3-25 21:23
破案了，是由于你是在Windows上运行的，运行命令不应该在扫描坐标前后写引号，即应该：D:\Relax_scan_too ...

按照您的提示，现在脚本已经能够成功提取扫描各店计算结果！
十分感谢您的耐心解答！您的脚本真的解决了我的燃眉之急！

ginlpein · 发表于 Post on 2022-3-26 23:02:59

本帖最后由 ginlpein 于 2022-3-26 23:28 编辑

Accelerator 发表于 2022-3-25 21:23
破案了，是由于你是在Windows上运行的，运行命令不应该在扫描坐标前后写引号，即应该：D:\Relax_scan_too ...

按照您的指导，在能量输出和部分结构输出方面均能正常执行功能。
但今天遇到负角度二面角与键角结构分析时，发现脚本无法解析并输出含负数扫描范围的结构，请问该如何解决？

具体命令行报错如下：

Now are you going to generate gjf files for a range of points? y/n
y
Range for coord1: e.g. 1.0-2.0
-20.0--40.0
Range for coord2: e.g. 1.0-2.0
110.0-130.0
Your kwds:
# opt=modredundant nosymm geom=connectivity pm7
Your cores:
28
Your memory:
75GB
Your charge:
0
Your multiplicity:
0
Now reading end.txt for the end part of your gjf file.
No end.txt found.
Traceback (most recent call last):
File "anaScan2D_3.py", line 98, in <module>
if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
File "anaScan2D_3.py", line 65, in withinRange
if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
IndexError: list index out of range

复制代码

Now are you going to generate gjf files for a range of points? y/n
y
Range for coord1: e.g. 1.0-2.0
(-20.0)-(-40.0)
Range for coord2: e.g. 1.0-2.0
110.0-130.0
Your kwds:
# opt=modredundant nosymm geom=connectivity pm7
Your cores:
28
Your memory:
75GB
Your charge:
0
Your multiplicity:
0
Now reading end.txt for the end part of your gjf file.
No end.txt found.
Traceback (most recent call last):
File "anaScan2D_3.py", line 98, in <module>
if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
File "anaScan2D_3.py", line 65, in withinRange
if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
IndexError: list index out of range

复制代码

我以为是74与76行分割符号的问题，尝试将“-”改为其他符号作为分割

rangeR1 = input().split('-')

复制代码

比如","与"~"

rangeR1 = input().split(',')

复制代码

rangeR1 = input().split('~')

复制代码

结果只是命令行不在报错，但是并未生成任何gjf文件，依旧无法输出结构。具体如下：

Now are you going to generate gjf files for a range of points? y/n
y
Range for coord1: e.g. 1.0,2.0
-20.0,-40.0
Range for coord2: e.g. 1.0,2.0
110.0,130.0
Your kwds:
# opt=modredundant nosymm geom=connectivity pm7
Your cores:
28
Your memory:
75GB
Your charge:
0
Your multiplicity:
0
Now reading end.txt for the end part of your gjf file.
No end.txt found.

复制代码

似乎脚本在输出结构时不支持负数输入。
由于我python水平为零，再不知道能如何改动，希望您不吝赐教。

Accelerator · 发表于 Post on 2022-3-28 09:19:43

ginlpein 发表于 2022-3-26 23:02
按照您的指导，在能量输出和部分结构输出方面均能正常执行功能。
但今天遇到负角度二面角与键角结构分析 ...

将73到76行的分隔符改成~就好了：
print('Range for coord1: e.g. 1.0~2.0')
rangeR1 = input().split('~')
print('Range for coord2: e.g. 1.0~2.0')
rangeR2 = input().split('~')

ginlpein · 发表于 Post on 2022-3-28 15:53:36

Accelerator 发表于 2022-3-28 09:19
将73到76行的分隔符改成~就好了：
print('Range for coord1: e.g. 1.0~2.0')
rangeR1 = input().split( ...

您好，修改分割符以后脚本不再报错，但是执行输出gjf文件的时候还是有些问题。
当两个扫描坐标范围都在正数区间的时候，可以正常输出；但是两个范围任一涉及负数区间的时候脚本情况就是如下，不报错，但是也没有输出任何gjf文件。
不知道是哪的问题，还是我负数输入方式的问题？

D:\calculate\FOR_GWB\close_to_3rd\inner>python3 anaScan2D_3.py scan.log D(56,8,15,2) D(55,54,84,119)
This script is to analyze a 2D scan which is not correctly completed. It will extract all the completed points so that you could still draw a 2D-PES figure.
R1= -46.9472 , R2= 83.1563 , E= 0.298256
R1= -36.9473 , R2= 83.1563 , E= 0.299025
R1= -26.9474 , R2= 83.1563 , E= 0.301488
R1= -16.9474 , R2= 83.1563 , E= 0.305661
R1= -6.9474 , R2= 83.1563 , E= 0.311482
R1= 3.0526 , R2= 83.1563 , E= 0.311322
R1= 13.0526 , R2= 83.1563 , E= 0.317850
R1= 23.0526 , R2= 83.1563 , E= 0.325014
R1= 33.0527 , R2= 83.1563 , E= 0.333035
R1= 43.0527 , R2= 83.1563 , E= 0.341850
R1= 53.0527 , R2= 83.1563 , E= 0.350731
R1= 63.0527 , R2= 83.1563 , E= 0.359668
R1= 73.0527 , R2= 83.1563 , E= 0.369227
R1= 73.0527 , R2= 93.1563 , E= 0.369834
R1= 63.0527 , R2= 93.1563 , E= 0.360400
R1= 53.0527 , R2= 93.1563 , E= 0.351517
R1= 43.0526 , R2= 93.1563 , E= 0.343413
R1= 33.0527 , R2= 93.1563 , E= 0.333904
R1= 23.0527 , R2= 93.1563 , E= 0.325531
R1= 13.0526 , R2= 93.1563 , E= 0.317831
R1= 3.0526 , R2= 93.1563 , E= 0.311608
R1= -6.9474 , R2= 93.1563 , E= 0.306729
R1= -16.9474 , R2= 93.1563 , E= 0.303177
R1= -26.9474 , R2= 93.1563 , E= 0.300998
R1= -36.9473 , R2= 93.1563 , E= 0.300176
R1= -46.9471 , R2= 93.1563 , E= 0.300717
R1= -46.9471 , R2= 103.1563 , E= 0.301174
R1= -36.9473 , R2= 103.1562 , E= 0.300646
R1= -26.9474 , R2= 103.1562 , E= 0.301486
R1= -16.9474 , R2= 103.1563 , E= 0.303695
R1= -6.9474 , R2= 103.1563 , E= 0.307286
R1= 3.0526 , R2= 103.1563 , E= 0.312163
R1= 13.0526 , R2= 103.1563 , E= 0.318459
R1= 23.0527 , R2= 103.1563 , E= 0.326116
R1= 33.0527 , R2= 103.1563 , E= 0.334860
R1= 43.0527 , R2= 103.1563 , E= 0.344302
R1= 53.0527 , R2= 103.1563 , E= 0.353789
R1= 73.0527 , R2= 103.1563 , E= 0.370614
R1= 73.0527 , R2= 113.1563 , E= 0.371791
R1= 63.0527 , R2= 113.1563 , E= 0.362528
R1= 53.0527 , R2= 113.1563 , E= 0.353673
R1= 43.0527 , R2= 113.1563 , E= 0.345912
R1= 33.0526 , R2= 113.1563 , E= 0.338379
R1= 23.0526 , R2= 113.1563 , E= 0.331041
R1= 13.0526 , R2= 113.1563 , E= 0.319474
R1= 3.0526 , R2= 113.1563 , E= 0.313089
R1= -6.9474 , R2= 113.1563 , E= 0.308112
R1= -16.9474 , R2= 113.1563 , E= 0.304521
R1= -26.9473 , R2= 113.1562 , E= 0.302332
R1= -36.9473 , R2= 113.1562 , E= 0.301517
R1= -46.9471 , R2= 113.1562 , E= 0.302052
R1= -46.9472 , R2= 123.1562 , E= 0.299854
R1= -36.9473 , R2= 123.1562 , E= 0.300146
R1= -26.9474 , R2= 123.1563 , E= 0.302254
R1= -16.9474 , R2= 123.1563 , E= 0.306162
R1= -6.9474 , R2= 123.1563 , E= 0.311791
R1= 3.0526 , R2= 123.1563 , E= 0.316973
R1= 13.0526 , R2= 123.1563 , E= 0.320516
R1= 23.0527 , R2= 123.1563 , E= 0.328606
R1= 33.0527 , R2= 123.1563 , E= 0.335376
R1= 43.0527 , R2= 123.1563 , E= 0.344265
R1= 53.0527 , R2= 123.1563 , E= 0.353539
R1= 63.0527 , R2= 123.1563 , E= 0.363006
R1= 73.0526 , R2= 123.1563 , E= 0.372874
R1= 73.0526 , R2= 133.1563 , E= 0.374837
R1= 63.0526 , R2= 133.1563 , E= 0.364563
R1= 53.0527 , R2= 133.1563 , E= 0.354669
R1= 43.0527 , R2= 133.1563 , E= 0.345287
R1= 33.0527 , R2= 133.1563 , E= 0.336386
R1= 23.0527 , R2= 133.1563 , E= 0.328478
R1= 13.0526 , R2= 133.1563 , E= 0.321955
R1= 3.0527 , R2= 133.1563 , E= 0.316932
R1= -6.9474 , R2= 133.1562 , E= 0.313386
R1= -16.9474 , R2= 133.1562 , E= 0.307243
R1= -26.9473 , R2= 133.1562 , E= 0.305265
R1= -36.9473 , R2= 133.1561 , E= 0.304705
R1= -46.9472 , R2= 133.1561 , E= 0.301513
R1= -46.9472 , R2= 143.1561 , E= 0.304642
R1= -36.9473 , R2= 143.1562 , E= 0.305020
R1= -26.9473 , R2= 143.1562 , E= 0.307327
R1= -16.9474 , R2= 143.1562 , E= 0.311434
R1= -6.9474 , R2= 143.1562 , E= 0.317234
R1= 3.0526 , R2= 143.1562 , E= 0.318421
R1= 13.0527 , R2= 143.1563 , E= 0.324911
R1= 23.0527 , R2= 143.1563 , E= 0.333027
R1= 33.0527 , R2= 143.1563 , E= 0.338979
R1= 43.0527 , R2= 143.1563 , E= 0.347943
R1= 53.0526 , R2= 143.1563 , E= 0.357507
R1= 63.0526 , R2= 143.1563 , E= 0.367679
R1= 73.0526 , R2= 143.1563 , E= 0.378483
R1= 73.0526 , R2= 153.1562 , E= 0.383749
R1= 63.0526 , R2= 153.1562 , E= 0.372724
R1= 53.0526 , R2= 153.1562 , E= 0.362320
R1= 43.0527 , R2= 153.1562 , E= 0.352556
R1= 33.0527 , R2= 153.1562 , E= 0.343399
R1= 23.0527 , R2= 153.1562 , E= 0.335176
R1= 13.0527 , R2= 153.1562 , E= 0.328373
R1= 3.0527 , R2= 153.1562 , E= 0.323190
R1= -6.9473 , R2= 153.1562 , E= 0.319705
R1= -16.9473 , R2= 153.1562 , E= 0.317776
R1= -26.9473 , R2= 153.1562 , E= 0.313159
R1= -36.9472 , R2= 153.1562 , E= 0.312718
R1= -46.9472 , R2= 153.1561 , E= 0.313721
R1= -46.9472 , R2= 163.1561 , E= 0.320420
R1= -36.9472 , R2= 163.1561 , E= 0.319487
R1= -26.9473 , R2= 163.1561 , E= 0.320033
R1= -16.9473 , R2= 163.1561 , E= 0.322125
R1= -6.9473 , R2= 163.1562 , E= 0.325612
R1= 3.0527 , R2= 163.1562 , E= 0.330571
R1= 13.0527 , R2= 163.1562 , E= 0.337085
R1= 23.0527 , R2= 163.1562 , E= 0.341124
R1= 33.0527 , R2= 163.1562 , E= 0.349786
R1= 43.0527 , R2= 163.1562 , E= 0.359296
R1= 53.0526 , R2= 163.1562 , E= 0.369245
R1= 63.0526 , R2= 163.1562 , E= 0.379761
R1= 73.0526 , R2= 163.1562 , E= 0.390817
R1= 73.0526 , R2= 173.1562 , E= 0.399401
R1= 63.0526 , R2= 173.1562 , E= 0.388648
R1= 53.0527 , R2= 173.1562 , E= 0.378242
R1= 43.0526 , R2= 173.1562 , E= 0.368121
R1= 33.0526 , R2= 173.1562 , E= 0.357865
R1= 23.0527 , R2= 173.1562 , E= 0.348798
R1= 13.0527 , R2= 173.1562 , E= 0.341358
R1= 3.0527 , R2= 173.1562 , E= 0.335582
R1= -6.9473 , R2= 173.1562 , E= 0.331550
R1= -16.9473 , R2= 173.1562 , E= 0.329229
R1= -26.9473 , R2= 173.1562 , E= 0.328304
R1= -36.9473 , R2= 173.1561 , E= 0.327559
R1= -46.9472 , R2= 173.1561 , E= 0.327274
Now are you going to generate gjf files for a range of points? y/n
y
Range for coord1: e.g. 1.0~2.0
-20.0~-40.0
Range for coord2: e.g. 1.0~2.0
110.0~130.0
Your kwds:
# opt=modredundant nosymm geom=connectivity pm7
Your cores:
28
Your memory:
75GB
Your charge:
0
Your multiplicity:
0
Now reading end.txt for the end part of your gjf file.
No end.txt found.

复制代码

这是按您提示更改分隔符之后使用的脚本

import os, sys
print('This script is to analyze a 2D scan which is not correctly completed. \
It will extract all the completed points so that you could still draw a 2D-PES figure.')
_, path, r1, r2 = sys.argv
r1s = []
r2s = []
es = []
coords = []
lineNum = 0
from common import readLog
elementTable = ['pass','H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si','P','S','Cl','Ar','K','Ca',\
'Sc','Ti','V','Cr','Mn','Fe','Co','Ni','Cu','Zm','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr','Nb',\
'Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I','Xe','Cs','Ba',\
'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf',\
'Ta','W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn']
import re
with open(path) as f:
etemp = 0
isStation = False
flagGeo = 0
thisCoord = []
previousCoord = []
while True:
line=f.readline()
lineNum += 1
if line == "":
break
if "SCF Done" in line:
etemp = float(line.split('=')[1].split("A")[0].strip())
continue
if "point found" in line:
isStation = True
continue
if 'Input Orientation' in line:
flagGeo = 1
continue
if flagGeo == 1 and (re.match("\s*[0-9]+\s*[0-9]+\s*[0-9]+\s*\-*[0-9]+\.",line) != None):
lsplitted = line.split()
element = elementTable[int(lsplitted[1])]
thisCoord.append([element] + lsplitted[-3:])
if (r1 in line) and isStation:
r1s.append(float(line.split(r1)[1].split('-DE')[0].strip()))
es.append(etemp)
continue
if (r2 in line) and isStation:
r2s.append(float(line.split(r2)[1].split('-DE')[0].strip()))
continue
if "Rotational" in line:
flagGeo = 0
if isStation:
coords.append(previousCoord)
previousCoord = thisCoord[:]
thisCoord = []
isStation = False
i=0
for e in es:
print('R1= {:.4f} , R2= {:.4f} , E= {:.6f} '.format(r1s[i],r2s[i],e))
i+=1
def withinRange(value,rangeList):
if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
return True
from common import writeGjf
#def writeGjf(path,coords,kwd,core,mem,charge,mult,end,chk='')
print('Now are you going to generate gjf files for a range of points? y/n')
if(input() != 'y'):
exit()
print('Range for coord1: e.g. 1.0~2.0')
rangeR1 = input().split('~')
print('Range for coord2: e.g. 1.0~2.0')
rangeR2 = input().split('~')
print('Your kwds:')
kwd = input()
print('Your cores:')
core = input()
print('Your memory:')
mem = input()
print('Your charge:')
charge = input()
print('Your multiplicity:')
mult = input()
print('Now reading end.txt for the end part of your gjf file.')
end = ''
try:
with open('end.txt') as f:
end = f.read()
except:
print('No end.txt found.')
itemNum = 0
for i in coords:
r1 = r1s[itemNum]
r2 = r2s[itemNum]
if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
writeGjf(path.replace('.log',f'_{r1}_{r2}.gjf'),i,kwd,core,mem,charge,mult,end)
itemNum += 1

复制代码

Accelerator · 发表于 Post on 2022-3-28 17:18:29

ginlpein 发表于 2022-3-28 15:53
您好，修改分割符以后脚本不再报错，但是执行输出gjf文件的时候还是有些问题。
当两个扫描坐标范围都在 ...

范围要从小到大输入，是-40.0~-20.0。
另外自旋多重度应该是1不是0

ginlpein · 发表于 Post on 2022-3-28 20:05:26

Accelerator 发表于 2022-3-28 17:18
范围要从小到大输入，是-40.0~-20.0。
另外自旋多重度应该是1不是0

噢！感谢大神！
平时都分析二面角大小，陷入刻板印象了，一直没反应过来-20大于-40。。。
感谢楼主的耐心解答！

ginlpein · 发表于 Post on 2022-3-28 21:04:16

本帖最后由 ginlpein 于 2022-3-28 21:06 编辑

Accelerator 发表于 2022-3-28 17:18
范围要从小到大输入，是-40.0~-20.0。
另外自旋多重度应该是1不是0

分隔符修改以后，能够成功定位到各点了，但对大体系输出，最后的gjf文件里面没有原子坐标，只有在命令行中输入的计算条件、电荷、多重度这些信息。命令行里面的输出还是原来那样。
我尝试把common.py里面的Standard Orientation也全部改换成Input Orientation，但是没有效果，各点具体的原子坐标还是没能提取出来，只有一个空壳。

%mem=1gb
%nprocs=1
# 1
TC
0 1

复制代码

Accelerator · 发表于 Post on 2022-3-29 09:16:38

ginlpein 发表于 2022-3-28 21:04
分隔符修改以后，能够成功定位到各点了，但对大体系输出，最后的gjf文件里面没有原子坐标 ...

要替换anaScan2D.py里的Standard Orientation，我这里是可以正常提取的。
代码里从common.py导入了一个函数readLog，但实际上并没有用到，所以common.py里的修改不会影响程序运行。
我现在在1L更新了脚本。

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[辅助/分析程序] Gaussian二维扫描输出文件的后处理脚本

评分 Rate