计算化学公社

 找回密码 Forget password
 注册 Register
Views: 12879|回复 Reply: 14

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

[复制链接 Copy URL]

422

帖子

7

威望

4772

eV
积分
5334

Level 6 (一方通行)

BSJ Institute

发表于 Post on 2021-4-30 11:46:40 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 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文件位置 两个扫描坐标。例如:
  1. python3 anaScan2D.py test.log 'R(1,2)' 'R(1,3)'
复制代码
不一会,就会输出当前已经成功扫过的每个点的两个扫描坐标数值以及能量。这样,即使扫描没有正确结束,也能利用已有的数据对势能面作图:
  1. 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.
  2. R1= 1.3943 , R2= 1.3914 , E= 0.161499
  3. R1= 1.3956 , R2= 1.3912 , E= 0.140213
  4. R1= 1.3968 , R2= 1.3908 , E= 0.136637
  5. R1= 1.3979 , R2= 1.3906 , E= 0.143130
复制代码
2. 接下来会询问是否需要提取其中部分点的坐标,写y,按照提示输入要提取出来的两个坐标的范围,新生成的gjf文件的关键字,使用的核心数和内存大小(GB),电荷和自旋多重度,就会自动生成一系列.gjf文件。文件尾部的信息需要提前写在当前目录的end.txt中,程序会读取其内容放在各gjf文件末尾。
  1. Now are you going to generate gjf files for a range of points? y/n
  2. y
  3. Range for coord1: e.g. 1.0-2.0
  4. 1.0-1.5
  5. Range for coord2: e.g. 1.0-2.0
  6. 1.0-1.5
  7. Your kwds:
  8. opt=modredundant freq b3lyp em=gd3bj 6-31G(d)
  9. (以下略去)
复制代码
如果我在end.txt中写B 23 12 F(换行)B 22 23 F,就会生成一系列类似下面的文件,名称为test_XXX_YYY.gjf,XXX和YYY分别为该点对应的两个扫描坐标数值。
  1. %mem=24gb
  2. %nprocs=16
  3. # opt=modredundant freq b3lyp em=gd3bj 6-31G(d)

  4. TC

  5. 0 1
  6. (略去)
  7. B 23 12 F
  8. B 22 23 F
复制代码

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

评分 Rate

参与人数
Participants 7
威望 +1 eV +21 收起 理由
Reason
找雷 + 5 精品内容,感谢老哥提高我工作效率!
ginlpein + 3 十分感谢耐心解答!
ggdh + 5 谢谢分享
LEVELF + 1
Butadiene + 2
zsu007 + 5 赞!
sobereva + 1

查看全部评分 View all ratings

15

帖子

0

威望

203

eV
积分
218

Level 3 能力者

发表于 Post on 2022-3-24 18:24:47 | 显示全部楼层 Show all
本帖最后由 ginlpein 于 2022-3-24 18:27 编辑

您好,请问您这个脚本该怎么能够适用于对二面角扫描结果的提取?
计算使用的是PM7半经验计算对两个二面角进行二维扫描。按照您写的内容运行,扫描坐标不论以何种形式输入,似乎都没能提取出各点能量,直接跳到询问是否生成gjf那步了,请问是我哪里操作错误吗?或者针对二面角需要哪里作出对应修改?
是否生成gjf文件那边选择y最后也是什么都没有生成出来。
  1. D:\Relax_scan_tool\Python>python3 anaScan2D.py scan.log 'D(37,6,14,4)' 'D(38,39,81,143)'
  2. 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.
  3. Now are you going to generate gjf files for a range of points? y/n
  4. N

  5. D:\Relax_scan_tool\Python>python3 anaScan2D.py scan.log 'D61' 'D346'
  6. 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.
  7. Now are you going to generate gjf files for a range of points? y/n
  8. N
复制代码


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

422

帖子

7

威望

4772

eV
积分
5334

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-3-24 19:45:35 | 显示全部楼层 Show all
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步),传上来给我调试。

15

帖子

0

威望

203

eV
积分
218

Level 3 能力者

发表于 Post on 2022-3-25 15:56:45 | 显示全部楼层 Show all
本帖最后由 ginlpein 于 2022-3-26 23:26 编辑

您好,后面我尝试做了一个3*3的双二面角扫描,使用您的脚本依旧无法实现解析。
  1. D:\Relax_scan_tool\Python>python3 anaScan2D.py scan_test.log 'D(37,6,14,4)' 'D(38,39,81,143)'
  2. 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.
  3. Now are you going to generate gjf files for a range of points? y/n
  4. n
复制代码

而在log文件中搜索扫描坐标,确实可以搜索到结果,扫描坐标出现在如下位置及情况:
  1. ! D61   D(37,6,14,4)          -53.1497         Scan                            !
  2. ! D62   D(37,6,14,9)          127.3039         estimate D2E/DX2                !
  3. ! D63   D(41,6,14,4)          129.5928         estimate D2E/DX2                !
  4. ! D64   D(41,6,14,9)          -49.9536         estimate D2E/DX2                !
  5. ! D65   D(14,6,37,38)        -173.96           estimate D2E/DX2                !
  6. ! D66   D(14,6,37,373)          1.2749         estimate D2E/DX2                !
  7. ! D67   D(41,6,37,38)           3.3306         estimate D2E/DX2                !
  8. ! D68   D(41,6,37,373)        178.5655         estimate D2E/DX2                !
  9. ! D69   D(14,6,41,40)         175.002          estimate D2E/DX2                !
  10. ! D70   D(14,6,41,385)         -6.0362         estimate D2E/DX2                !
  11. ! D71   D(37,6,41,40)          -2.3562         estimate D2E/DX2                !
  12. ! D72   D(37,6,41,385)        176.6056         estimate D2E/DX2                !
  13. ! D73   D(67,7,16,5)          -53.1497         estimate D2E/DX2                !
复制代码
  1. ! D61   D(37,6,14,4)          -53.1497         -DE/DX =    0.0003              !
  2. ! D62   D(37,6,14,9)          125.2525         -DE/DX =    0.0001              !
  3. ! D63   D(41,6,14,4)          132.6627         -DE/DX =    0.0002              !
  4. ! D64   D(41,6,14,9)          -48.935          -DE/DX =    0.0                 !
  5. ! D65   D(14,6,37,38)        -168.2441         -DE/DX =   -0.0001              !
  6. ! D66   D(14,6,37,373)         11.9407         -DE/DX =   -0.0001              !
  7. ! D67   D(41,6,37,38)           5.8314         -DE/DX =    0.0                 !
  8. ! D68   D(41,6,37,373)       -173.9838         -DE/DX =    0.0                 !
复制代码

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

scan_test.7z

8.95 MB, 下载次数 Times of downloads: 1

422

帖子

7

威望

4772

eV
积分
5334

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-3-25 16:27:39 | 显示全部楼层 Show all
ginlpein 发表于 2022-3-25 15:56
您好,后面我尝试做了一个3*3的双二面角扫描,使用您的脚本依旧无法实现解析。

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

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

15

帖子

0

威望

203

eV
积分
218

Level 3 能力者

发表于 Post on 2022-3-25 20:57:39 | 显示全部楼层 Show all
本帖最后由 ginlpein 于 2022-3-25 20:58 编辑
Accelerator 发表于 2022-3-25 16:27
这是由于二面角为负值,原始代码中用-DE/-DX中的负号作为分隔符对行进行分割,遇到负数就没能正确识别。
...

感谢您的耐心回复!
不知道是不是因为我改的不对,我把44行与48行的'-'改成'-DE' 后对附件中的scan_test进行测试依旧无法生效,不知道问题出在哪?
  1. import os, sys
  2. print('This script is to analyze a 2D scan which is not correctly completed. \
  3.         It will extract all the completed points so that you could still draw a 2D-PES figure.')

  4. _, path, r1, r2 = sys.argv
  5. r1s = []
  6. r2s = []
  7. es = []
  8. coords = []
  9. lineNum = 0
  10. from common import readLog
  11. elementTable = ['pass','H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si','P','S','Cl','Ar','K','Ca',\
  12. 'Sc','Ti','V','Cr','Mn','Fe','Co','Ni','Cu','Zm','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr','Nb',\
  13. 'Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I','Xe','Cs','Ba',\
  14. 'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf',\
  15. 'Ta','W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn']

  16. import re
  17. with open(path) as f:
  18.         etemp = 0
  19.         isStation = False
  20.         flagGeo = 0
  21.         thisCoord = []
  22.         previousCoord = []
  23.         while True:
  24.                 line=f.readline()
  25.                 lineNum += 1
  26.                 if line == "":
  27.                         break
  28.                 if "SCF Done" in line:
  29.                         etemp = float(line.split('=')[1].split("A")[0].strip())
  30.                         continue
  31.                 if "point found" in line:
  32.                         isStation = True
  33.                         continue
  34.                 if 'Standard orientation' in line:
  35.                         flagGeo = 1
  36.                         continue
  37.                 if flagGeo == 1 and (re.match("\s*[0-9]+\s*[0-9]+\s*[0-9]+\s*\-*[0-9]+\.",line) != None):
  38.                         lsplitted = line.split()
  39.                         element = elementTable[int(lsplitted[1])]
  40.                         thisCoord.append([element] + lsplitted[-3:])
  41.                 if (r1 in line) and isStation:
  42.                         r1s.append(float(line.split(r1)[1].split('-DE')[0].strip()))
  43.                         es.append(etemp)
  44.                         continue
  45.                 if (r2 in line) and isStation:
  46.                         r2s.append(float(line.split(r2)[1].split('-DE')[0].strip()))
  47.                         continue
  48.                 if "Rotational" in line:
  49.                         flagGeo = 0
  50.                         if isStation:
  51.                                 coords.append(previousCoord)
  52.                         previousCoord = thisCoord[:]
  53.                         thisCoord = []
  54.                         isStation = False

  55. i=0

  56. for e in es:
  57.         print('R1= {:.4f} , R2= {:.4f} , E= {:.6f} '.format(r1s[i],r2s[i],e))
  58.         i+=1

  59. def withinRange(value,rangeList):
  60.         if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
  61.                 return True

  62. from common import writeGjf
  63. #def writeGjf(path,coords,kwd,core,mem,charge,mult,end,chk='')
  64. print('Now are you going to generate gjf files for a range of points? y/n')
  65. if(input() != 'y'):
  66.         exit()
  67. print('Range for coord1: e.g. 1.0-2.0')
  68. rangeR1 = input().split('-')
  69. print('Range for coord2: e.g. 1.0-2.0')
  70. rangeR2 = input().split('-')
  71. print('Your kwds:')
  72. kwd = input()
  73. print('Your cores:')
  74. core = input()
  75. print('Your memory:')
  76. mem = input()
  77. print('Your charge:')
  78. charge = input()
  79. print('Your multiplicity:')
  80. mult = input()
  81. print('Now reading end.txt for the end part of your gjf file.')
  82. end = ''
  83. try:
  84.         with open('end.txt') as f:
  85.                 end = f.read()
  86. except:
  87.         print('No end.txt found.')
  88. itemNum = 0
  89. for i in coords:
  90.         r1 = r1s[itemNum]
  91.         r2 = r2s[itemNum]
  92.         if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
  93.                 writeGjf(path.replace('.log',f'_{r1}_{r2}.gjf'),i,kwd,core,mem,charge,mult,end)
  94.         itemNum += 1
复制代码

  1. D:\Relax_scan_tool\Python>python3 anaScan2D_2.py scan_test.log 'D(37,6,14,4)' 'D(38,39,81,143)'
  2. 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.
  3. Now are you going to generate gjf files for a range of points? y/n
  4. N
复制代码

422

帖子

7

威望

4772

eV
积分
5334

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-3-25 21:23:40 | 显示全部楼层 Show all
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

15

帖子

0

威望

203

eV
积分
218

Level 3 能力者

发表于 Post on 2022-3-25 22:38:25 | 显示全部楼层 Show all
本帖最后由 ginlpein 于 2022-3-25 23:11 编辑
Accelerator 发表于 2022-3-25 21:23
破案了,是由于你是在Windows上运行的,运行命令不应该在扫描坐标前后写引号,即应该:D:\Relax_scan_too ...


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

15

帖子

0

威望

203

eV
积分
218

Level 3 能力者

发表于 Post on 2022-3-26 23:02:59 | 显示全部楼层 Show all
本帖最后由 ginlpein 于 2022-3-26 23:28 编辑
Accelerator 发表于 2022-3-25 21:23
破案了,是由于你是在Windows上运行的,运行命令不应该在扫描坐标前后写引号,即应该:D:\Relax_scan_too ...

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

具体命令行报错如下:
  1. Now are you going to generate gjf files for a range of points? y/n
  2. y
  3. Range for coord1: e.g. 1.0-2.0
  4. -20.0--40.0
  5. Range for coord2: e.g. 1.0-2.0
  6. 110.0-130.0
  7. Your kwds:
  8. # opt=modredundant nosymm geom=connectivity pm7
  9. Your cores:
  10. 28
  11. Your memory:
  12. 75GB
  13. Your charge:
  14. 0
  15. Your multiplicity:
  16. 0
  17. Now reading end.txt for the end part of your gjf file.
  18. No end.txt found.
  19. Traceback (most recent call last):
  20.   File "anaScan2D_3.py", line 98, in <module>
  21.     if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
  22.   File "anaScan2D_3.py", line 65, in withinRange
  23.     if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
  24. IndexError: list index out of range
复制代码

  1. Now are you going to generate gjf files for a range of points? y/n
  2. y
  3. Range for coord1: e.g. 1.0-2.0
  4. (-20.0)-(-40.0)
  5. Range for coord2: e.g. 1.0-2.0
  6. 110.0-130.0
  7. Your kwds:
  8. # opt=modredundant nosymm geom=connectivity pm7
  9. Your cores:
  10. 28
  11. Your memory:
  12. 75GB
  13. Your charge:
  14. 0
  15. Your multiplicity:
  16. 0
  17. Now reading end.txt for the end part of your gjf file.
  18. No end.txt found.
  19. Traceback (most recent call last):
  20.   File "anaScan2D_3.py", line 98, in <module>
  21.     if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
  22.   File "anaScan2D_3.py", line 65, in withinRange
  23.     if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
  24. IndexError: list index out of range
复制代码


我以为是74与76行分割符号的问题,尝试将“-”改为其他符号作为分割
  1. rangeR1 = input().split('-')
复制代码
比如","与"~"
  1. rangeR1 = input().split(',')
复制代码
  1. rangeR1 = input().split('~')
复制代码
结果只是命令行不在报错,但是并未生成任何gjf文件,依旧无法输出结构。具体如下:
  1. Now are you going to generate gjf files for a range of points? y/n
  2. y
  3. Range for coord1: e.g. 1.0,2.0
  4. -20.0,-40.0
  5. Range for coord2: e.g. 1.0,2.0
  6. 110.0,130.0
  7. Your kwds:
  8. # opt=modredundant nosymm geom=connectivity pm7
  9. Your cores:
  10. 28
  11. Your memory:
  12. 75GB
  13. Your charge:
  14. 0
  15. Your multiplicity:
  16. 0
  17. Now reading end.txt for the end part of your gjf file.
  18. No end.txt found.
复制代码
似乎脚本在输出结构时不支持负数输入。
由于我python水平为零,再不知道能如何改动,希望您不吝赐教。

422

帖子

7

威望

4772

eV
积分
5334

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-3-28 09:19:43 | 显示全部楼层 Show all
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('~')

15

帖子

0

威望

203

eV
积分
218

Level 3 能力者

发表于 Post on 2022-3-28 15:53:36 | 显示全部楼层 Show all
Accelerator 发表于 2022-3-28 09:19
将73到76行的分隔符改成~就好了:
print('Range for coord1: e.g. 1.0~2.0')
rangeR1 = input().split( ...

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

  1. D:\calculate\FOR_GWB\close_to_3rd\inner>python3 anaScan2D_3.py scan.log D(56,8,15,2) D(55,54,84,119)
  2. 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.
  3. R1= -46.9472 , R2= 83.1563 , E= 0.298256
  4. R1= -36.9473 , R2= 83.1563 , E= 0.299025
  5. R1= -26.9474 , R2= 83.1563 , E= 0.301488
  6. R1= -16.9474 , R2= 83.1563 , E= 0.305661
  7. R1= -6.9474 , R2= 83.1563 , E= 0.311482
  8. R1= 3.0526 , R2= 83.1563 , E= 0.311322
  9. R1= 13.0526 , R2= 83.1563 , E= 0.317850
  10. R1= 23.0526 , R2= 83.1563 , E= 0.325014
  11. R1= 33.0527 , R2= 83.1563 , E= 0.333035
  12. R1= 43.0527 , R2= 83.1563 , E= 0.341850
  13. R1= 53.0527 , R2= 83.1563 , E= 0.350731
  14. R1= 63.0527 , R2= 83.1563 , E= 0.359668
  15. R1= 73.0527 , R2= 83.1563 , E= 0.369227
  16. R1= 73.0527 , R2= 93.1563 , E= 0.369834
  17. R1= 63.0527 , R2= 93.1563 , E= 0.360400
  18. R1= 53.0527 , R2= 93.1563 , E= 0.351517
  19. R1= 43.0526 , R2= 93.1563 , E= 0.343413
  20. R1= 33.0527 , R2= 93.1563 , E= 0.333904
  21. R1= 23.0527 , R2= 93.1563 , E= 0.325531
  22. R1= 13.0526 , R2= 93.1563 , E= 0.317831
  23. R1= 3.0526 , R2= 93.1563 , E= 0.311608
  24. R1= -6.9474 , R2= 93.1563 , E= 0.306729
  25. R1= -16.9474 , R2= 93.1563 , E= 0.303177
  26. R1= -26.9474 , R2= 93.1563 , E= 0.300998
  27. R1= -36.9473 , R2= 93.1563 , E= 0.300176
  28. R1= -46.9471 , R2= 93.1563 , E= 0.300717
  29. R1= -46.9471 , R2= 103.1563 , E= 0.301174
  30. R1= -36.9473 , R2= 103.1562 , E= 0.300646
  31. R1= -26.9474 , R2= 103.1562 , E= 0.301486
  32. R1= -16.9474 , R2= 103.1563 , E= 0.303695
  33. R1= -6.9474 , R2= 103.1563 , E= 0.307286
  34. R1= 3.0526 , R2= 103.1563 , E= 0.312163
  35. R1= 13.0526 , R2= 103.1563 , E= 0.318459
  36. R1= 23.0527 , R2= 103.1563 , E= 0.326116
  37. R1= 33.0527 , R2= 103.1563 , E= 0.334860
  38. R1= 43.0527 , R2= 103.1563 , E= 0.344302
  39. R1= 53.0527 , R2= 103.1563 , E= 0.353789
  40. R1= 73.0527 , R2= 103.1563 , E= 0.370614
  41. R1= 73.0527 , R2= 113.1563 , E= 0.371791
  42. R1= 63.0527 , R2= 113.1563 , E= 0.362528
  43. R1= 53.0527 , R2= 113.1563 , E= 0.353673
  44. R1= 43.0527 , R2= 113.1563 , E= 0.345912
  45. R1= 33.0526 , R2= 113.1563 , E= 0.338379
  46. R1= 23.0526 , R2= 113.1563 , E= 0.331041
  47. R1= 13.0526 , R2= 113.1563 , E= 0.319474
  48. R1= 3.0526 , R2= 113.1563 , E= 0.313089
  49. R1= -6.9474 , R2= 113.1563 , E= 0.308112
  50. R1= -16.9474 , R2= 113.1563 , E= 0.304521
  51. R1= -26.9473 , R2= 113.1562 , E= 0.302332
  52. R1= -36.9473 , R2= 113.1562 , E= 0.301517
  53. R1= -46.9471 , R2= 113.1562 , E= 0.302052
  54. R1= -46.9472 , R2= 123.1562 , E= 0.299854
  55. R1= -36.9473 , R2= 123.1562 , E= 0.300146
  56. R1= -26.9474 , R2= 123.1563 , E= 0.302254
  57. R1= -16.9474 , R2= 123.1563 , E= 0.306162
  58. R1= -6.9474 , R2= 123.1563 , E= 0.311791
  59. R1= 3.0526 , R2= 123.1563 , E= 0.316973
  60. R1= 13.0526 , R2= 123.1563 , E= 0.320516
  61. R1= 23.0527 , R2= 123.1563 , E= 0.328606
  62. R1= 33.0527 , R2= 123.1563 , E= 0.335376
  63. R1= 43.0527 , R2= 123.1563 , E= 0.344265
  64. R1= 53.0527 , R2= 123.1563 , E= 0.353539
  65. R1= 63.0527 , R2= 123.1563 , E= 0.363006
  66. R1= 73.0526 , R2= 123.1563 , E= 0.372874
  67. R1= 73.0526 , R2= 133.1563 , E= 0.374837
  68. R1= 63.0526 , R2= 133.1563 , E= 0.364563
  69. R1= 53.0527 , R2= 133.1563 , E= 0.354669
  70. R1= 43.0527 , R2= 133.1563 , E= 0.345287
  71. R1= 33.0527 , R2= 133.1563 , E= 0.336386
  72. R1= 23.0527 , R2= 133.1563 , E= 0.328478
  73. R1= 13.0526 , R2= 133.1563 , E= 0.321955
  74. R1= 3.0527 , R2= 133.1563 , E= 0.316932
  75. R1= -6.9474 , R2= 133.1562 , E= 0.313386
  76. R1= -16.9474 , R2= 133.1562 , E= 0.307243
  77. R1= -26.9473 , R2= 133.1562 , E= 0.305265
  78. R1= -36.9473 , R2= 133.1561 , E= 0.304705
  79. R1= -46.9472 , R2= 133.1561 , E= 0.301513
  80. R1= -46.9472 , R2= 143.1561 , E= 0.304642
  81. R1= -36.9473 , R2= 143.1562 , E= 0.305020
  82. R1= -26.9473 , R2= 143.1562 , E= 0.307327
  83. R1= -16.9474 , R2= 143.1562 , E= 0.311434
  84. R1= -6.9474 , R2= 143.1562 , E= 0.317234
  85. R1= 3.0526 , R2= 143.1562 , E= 0.318421
  86. R1= 13.0527 , R2= 143.1563 , E= 0.324911
  87. R1= 23.0527 , R2= 143.1563 , E= 0.333027
  88. R1= 33.0527 , R2= 143.1563 , E= 0.338979
  89. R1= 43.0527 , R2= 143.1563 , E= 0.347943
  90. R1= 53.0526 , R2= 143.1563 , E= 0.357507
  91. R1= 63.0526 , R2= 143.1563 , E= 0.367679
  92. R1= 73.0526 , R2= 143.1563 , E= 0.378483
  93. R1= 73.0526 , R2= 153.1562 , E= 0.383749
  94. R1= 63.0526 , R2= 153.1562 , E= 0.372724
  95. R1= 53.0526 , R2= 153.1562 , E= 0.362320
  96. R1= 43.0527 , R2= 153.1562 , E= 0.352556
  97. R1= 33.0527 , R2= 153.1562 , E= 0.343399
  98. R1= 23.0527 , R2= 153.1562 , E= 0.335176
  99. R1= 13.0527 , R2= 153.1562 , E= 0.328373
  100. R1= 3.0527 , R2= 153.1562 , E= 0.323190
  101. R1= -6.9473 , R2= 153.1562 , E= 0.319705
  102. R1= -16.9473 , R2= 153.1562 , E= 0.317776
  103. R1= -26.9473 , R2= 153.1562 , E= 0.313159
  104. R1= -36.9472 , R2= 153.1562 , E= 0.312718
  105. R1= -46.9472 , R2= 153.1561 , E= 0.313721
  106. R1= -46.9472 , R2= 163.1561 , E= 0.320420
  107. R1= -36.9472 , R2= 163.1561 , E= 0.319487
  108. R1= -26.9473 , R2= 163.1561 , E= 0.320033
  109. R1= -16.9473 , R2= 163.1561 , E= 0.322125
  110. R1= -6.9473 , R2= 163.1562 , E= 0.325612
  111. R1= 3.0527 , R2= 163.1562 , E= 0.330571
  112. R1= 13.0527 , R2= 163.1562 , E= 0.337085
  113. R1= 23.0527 , R2= 163.1562 , E= 0.341124
  114. R1= 33.0527 , R2= 163.1562 , E= 0.349786
  115. R1= 43.0527 , R2= 163.1562 , E= 0.359296
  116. R1= 53.0526 , R2= 163.1562 , E= 0.369245
  117. R1= 63.0526 , R2= 163.1562 , E= 0.379761
  118. R1= 73.0526 , R2= 163.1562 , E= 0.390817
  119. R1= 73.0526 , R2= 173.1562 , E= 0.399401
  120. R1= 63.0526 , R2= 173.1562 , E= 0.388648
  121. R1= 53.0527 , R2= 173.1562 , E= 0.378242
  122. R1= 43.0526 , R2= 173.1562 , E= 0.368121
  123. R1= 33.0526 , R2= 173.1562 , E= 0.357865
  124. R1= 23.0527 , R2= 173.1562 , E= 0.348798
  125. R1= 13.0527 , R2= 173.1562 , E= 0.341358
  126. R1= 3.0527 , R2= 173.1562 , E= 0.335582
  127. R1= -6.9473 , R2= 173.1562 , E= 0.331550
  128. R1= -16.9473 , R2= 173.1562 , E= 0.329229
  129. R1= -26.9473 , R2= 173.1562 , E= 0.328304
  130. R1= -36.9473 , R2= 173.1561 , E= 0.327559
  131. R1= -46.9472 , R2= 173.1561 , E= 0.327274
  132. Now are you going to generate gjf files for a range of points? y/n
  133. y
  134. Range for coord1: e.g. 1.0~2.0
  135. -20.0~-40.0
  136. Range for coord2: e.g. 1.0~2.0
  137. 110.0~130.0
  138. Your kwds:
  139. # opt=modredundant nosymm geom=connectivity pm7
  140. Your cores:
  141. 28
  142. Your memory:
  143. 75GB
  144. Your charge:
  145. 0
  146. Your multiplicity:
  147. 0
  148. Now reading end.txt for the end part of your gjf file.
  149. No end.txt found.
复制代码


这是按您提示更改分隔符之后使用的脚本
  1. import os, sys
  2. print('This script is to analyze a 2D scan which is not correctly completed. \
  3.         It will extract all the completed points so that you could still draw a 2D-PES figure.')

  4. _, path, r1, r2 = sys.argv
  5. r1s = []
  6. r2s = []
  7. es = []
  8. coords = []
  9. lineNum = 0
  10. from common import readLog
  11. elementTable = ['pass','H','He','Li','Be','B','C','N','O','F','Ne','Na','Mg','Al','Si','P','S','Cl','Ar','K','Ca',\
  12. 'Sc','Ti','V','Cr','Mn','Fe','Co','Ni','Cu','Zm','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y','Zr','Nb',\
  13. 'Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I','Xe','Cs','Ba',\
  14. 'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb','Lu','Hf',\
  15. 'Ta','W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn']

  16. import re
  17. with open(path) as f:
  18.         etemp = 0
  19.         isStation = False
  20.         flagGeo = 0
  21.         thisCoord = []
  22.         previousCoord = []
  23.         while True:
  24.                 line=f.readline()
  25.                 lineNum += 1
  26.                 if line == "":
  27.                         break
  28.                 if "SCF Done" in line:
  29.                         etemp = float(line.split('=')[1].split("A")[0].strip())
  30.                         continue
  31.                 if "point found" in line:
  32.                         isStation = True
  33.                         continue
  34.                 if 'Input Orientation' in line:
  35.                         flagGeo = 1
  36.                         continue
  37.                 if flagGeo == 1 and (re.match("\s*[0-9]+\s*[0-9]+\s*[0-9]+\s*\-*[0-9]+\.",line) != None):
  38.                         lsplitted = line.split()
  39.                         element = elementTable[int(lsplitted[1])]
  40.                         thisCoord.append([element] + lsplitted[-3:])
  41.                 if (r1 in line) and isStation:
  42.                         r1s.append(float(line.split(r1)[1].split('-DE')[0].strip()))
  43.                         es.append(etemp)
  44.                         continue
  45.                 if (r2 in line) and isStation:
  46.                         r2s.append(float(line.split(r2)[1].split('-DE')[0].strip()))
  47.                         continue
  48.                 if "Rotational" in line:
  49.                         flagGeo = 0
  50.                         if isStation:
  51.                                 coords.append(previousCoord)
  52.                         previousCoord = thisCoord[:]
  53.                         thisCoord = []
  54.                         isStation = False

  55. i=0

  56. for e in es:
  57.         print('R1= {:.4f} , R2= {:.4f} , E= {:.6f} '.format(r1s[i],r2s[i],e))
  58.         i+=1

  59. def withinRange(value,rangeList):
  60.         if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
  61.                 return True

  62. from common import writeGjf
  63. #def writeGjf(path,coords,kwd,core,mem,charge,mult,end,chk='')
  64. print('Now are you going to generate gjf files for a range of points? y/n')
  65. if(input() != 'y'):
  66.         exit()
  67. print('Range for coord1: e.g. 1.0~2.0')
  68. rangeR1 = input().split('~')
  69. print('Range for coord2: e.g. 1.0~2.0')
  70. rangeR2 = input().split('~')
  71. print('Your kwds:')
  72. kwd = input()
  73. print('Your cores:')
  74. core = input()
  75. print('Your memory:')
  76. mem = input()
  77. print('Your charge:')
  78. charge = input()
  79. print('Your multiplicity:')
  80. mult = input()
  81. print('Now reading end.txt for the end part of your gjf file.')
  82. end = ''
  83. try:
  84.         with open('end.txt') as f:
  85.                 end = f.read()
  86. except:
  87.         print('No end.txt found.')
  88. itemNum = 0
  89. for i in coords:
  90.         r1 = r1s[itemNum]
  91.         r2 = r2s[itemNum]
  92.         if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
  93.                 writeGjf(path.replace('.log',f'_{r1}_{r2}.gjf'),i,kwd,core,mem,charge,mult,end)
  94.         itemNum += 1
复制代码

422

帖子

7

威望

4772

eV
积分
5334

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-3-28 17:18:29 | 显示全部楼层 Show all
ginlpein 发表于 2022-3-28 15:53
您好,修改分割符以后脚本不再报错,但是执行输出gjf文件的时候还是有些问题。
当两个扫描坐标范围都在 ...

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

15

帖子

0

威望

203

eV
积分
218

Level 3 能力者

发表于 Post on 2022-3-28 20:05:26 | 显示全部楼层 Show all
Accelerator 发表于 2022-3-28 17:18
范围要从小到大输入,是-40.0~-20.0。
另外自旋多重度应该是1不是0

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

15

帖子

0

威望

203

eV
积分
218

Level 3 能力者

发表于 Post on 2022-3-28 21:04:16 | 显示全部楼层 Show all
本帖最后由 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,但是没有效果,各点具体的原子坐标还是没能提取出来,只有一个空壳。

  1. %mem=1gb
  2. %nprocs=1
  3. # 1

  4. TC

  5. 0 1
复制代码
scan_test.png

422

帖子

7

威望

4772

eV
积分
5334

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-3-29 09:16:38 | 显示全部楼层 Show all
ginlpein 发表于 2022-3-28 21:04
分隔符修改以后,能够成功定位到各点了,但对大体系输出,最后的gjf文件里面没有原子坐标 ...

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

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 04:14 , Processed in 0.240342 second(s), 25 queries .

快速回复 返回顶部 返回列表 Return to list