计算化学公社

标题: Gaussian二维扫描输出文件的后处理脚本 [打印本页]
作者
Author: Accelerator    时间: 2021-4-30 11:46
标题: Gaussian二维扫描输出文件的后处理脚本
本帖最后由 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)

  5. TC

  7. 0 1
  8. （略去）
  9. B 23 12 F
  10. B 22 23 F
复制代码

(, 下载次数 Times of downloads: 76)


作者
Author: ginlpein    时间: 2022-3-24 18:24
本帖最后由 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

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


我的计算体系在200~500原子的级别，由于比较结构大，导致二维柔性扫描结果文件普遍大于1G，经常超过2G。超大的结果文件导致不论计算成功完成与否，gaussview都无法正常读取结果，Result-scan中无法正确显示或输出结果。尝试cclib等程序也不太支持柔性扫描，所以一直在寻找可以用于解析计算结果的方法，希望楼主赐教。
作者
Author: Accelerator    时间: 2022-3-24 19:45
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步），传上来给我调试。
作者
Author: ginlpein    时间: 2022-3-25 15:56
本帖最后由 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的级别，似乎不可避免的就是会很大。
作者
Author: Accelerator    时间: 2022-3-25 16:27
ginlpein 发表于 2022-3-25 15:56
您好，后面我尝试做了一个3*3的双二面角扫描，使用您的脚本依旧无法实现解析。

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

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

作者
Author: ginlpein    时间: 2022-3-25 20:57
本帖最后由 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.')

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

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

  58. i=0

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

  64. def withinRange(value,rangeList):
  65.         if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
  66.                 return True

  68. from common import writeGjf
  69. #def writeGjf(path,coords,kwd,core,mem,charge,mult,end,chk='')
  70. print('Now are you going to generate gjf files for a range of points? y/n')
  71. if(input() != 'y'):
  72.         exit()
  73. print('Range for coord1: e.g. 1.0-2.0')
  74. rangeR1 = input().split('-')
  75. print('Range for coord2: e.g. 1.0-2.0')
  76. rangeR2 = input().split('-')
  77. print('Your kwds:')
  78. kwd = input()
  79. print('Your cores:')
  80. core = input()
  81. print('Your memory:')
  82. mem = input()
  83. print('Your charge:')
  84. charge = input()
  85. print('Your multiplicity:')
  86. mult = input()
  87. print('Now reading end.txt for the end part of your gjf file.')
  88. end = ''
  89. try:
  90.         with open('end.txt') as f:
  91.                 end = f.read()
  92. except:
  93.         print('No end.txt found.')
  94. itemNum = 0
  95. for i in coords:
  96.         r1 = r1s[itemNum]
  97.         r2 = r2s[itemNum]
  98.         if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
  99.                 writeGjf(path.replace('.log',f'_{r1}_{r2}.gjf'),i,kwd,core,mem,charge,mult,end)
  100.         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
复制代码
作者
Author: Accelerator    时间: 2022-3-25 21:23
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
作者
Author: ginlpein    时间: 2022-3-25 22:38
本帖最后由 ginlpein 于 2022-3-25 23:11 编辑
Accelerator 发表于 2022-3-25 21:23
破案了，是由于你是在Windows上运行的，运行命令不应该在扫描坐标前后写引号，即应该：D:\Relax_scan_too ...


按照您的提示，现在脚本已经能够成功提取扫描各店计算结果！
十分感谢您的耐心解答！您的脚本真的解决了我的燃眉之急！
作者
Author: ginlpein    时间: 2022-3-26 23:02
本帖最后由 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水平为零，再不知道能如何改动，希望您不吝赐教。
作者
Author: Accelerator    时间: 2022-3-28 09:19
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('~')
作者
Author: ginlpein    时间: 2022-3-28 15:53
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.')

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

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

  58. i=0

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

  64. def withinRange(value,rangeList):
  65.         if float(value) < float(rangeList[1]) and float(value) > float(rangeList[0]):
  66.                 return True

  68. from common import writeGjf
  69. #def writeGjf(path,coords,kwd,core,mem,charge,mult,end,chk='')
  70. print('Now are you going to generate gjf files for a range of points? y/n')
  71. if(input() != 'y'):
  72.         exit()
  73. print('Range for coord1: e.g. 1.0~2.0')
  74. rangeR1 = input().split('~')
  75. print('Range for coord2: e.g. 1.0~2.0')
  76. rangeR2 = input().split('~')
  77. print('Your kwds:')
  78. kwd = input()
  79. print('Your cores:')
  80. core = input()
  81. print('Your memory:')
  82. mem = input()
  83. print('Your charge:')
  84. charge = input()
  85. print('Your multiplicity:')
  86. mult = input()
  87. print('Now reading end.txt for the end part of your gjf file.')
  88. end = ''
  89. try:
  90.         with open('end.txt') as f:
  91.                 end = f.read()
  92. except:
  93.         print('No end.txt found.')
  94. itemNum = 0
  95. for i in coords:
  96.         r1 = r1s[itemNum]
  97.         r2 = r2s[itemNum]
  98.         if withinRange(r1,rangeR1) and withinRange(r2,rangeR2):
  99.                 writeGjf(path.replace('.log',f'_{r1}_{r2}.gjf'),i,kwd,core,mem,charge,mult,end)
  100.         itemNum += 1
复制代码
作者
Author: Accelerator    时间: 2022-3-28 17:18
ginlpein 发表于 2022-3-28 15:53
您好，修改分割符以后脚本不再报错，但是执行输出gjf文件的时候还是有些问题。
当两个扫描坐标范围都在 ...

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

噢！感谢大神！
平时都分析二面角大小，陷入刻板印象了，一直没反应过来-20大于-40。。。
感谢楼主的耐心解答！
作者
Author: ginlpein    时间: 2022-3-28 21:04
本帖最后由 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

  5. TC

  7. 0 1
复制代码
作者
Author: Accelerator    时间: 2022-3-29 09:16
ginlpein 发表于 2022-3-28 21:04
分隔符修改以后，能够成功定位到各点了，但对大体系输出，最后的gjf文件里面没有原子坐标 ...

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



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3