gaussian和MOKIT联用脚本几何优化失败（已解决，脚本可自取）

shenzp

后来莫名发现对于NO分子，数值梯度也可以做了，不知道是不是在无意识的时候改了一下脚本(x)。不过对于较大的体系会很吃力，gaussian的数值梯度感觉也不是特别好用，（也有可能是我的脚本的问题）所以还是慎用。数值梯度的输入文件也已放入附件中。

----------------------------分割线----------------------------------
第二天：折腾了好久，百思不得其解，偶然发现是读取梯度数据时的正则表达式错了，没有读取E指数部分，已在原文中修改（这部分代码是让chatgpt写的，chatgpt写的代码果然容易出小bug不会写代码真是痛苦）。修改后经简单测试能正常运行。另外，对于没有解析梯度的方法（如NEVPT2），经过简单测试，通过输入文件写opt=(enonly,ef)、脚本算单点的方式去做，优化会失败，在两个结构间反复横跳。（一开始其实是想做NEVPT2优化的脚本的，所以中间文件名都用了NEVPT2，尴尬）不过目前觉得还是用CASSCF优化比较现实，所以暂时也懒得去查问题出在哪儿，若有老师能指点，还是不吝赐教。帖子还是放在这儿，在MOKIT官方的结构优化功能实现前，有需要的小伙伴还是可以用有解析梯度的方法做优化和数值频率脚本和测试的输入文件都放在了附件中。

----------------------------分割线----------------------------------
各位老师好，我最近想做CASSCF的几何优化，苦于molpro的预编译版在新电脑上无法运行，其它多参考程序又不会用，于是决定斗胆写一个将gaussian和邹神的MOKIT联用做几何优化的脚本。脚本是在Gaussian-MOKIT-Molpro 联用方案 gmm.py（http://bbs.keinsci.com/forum.php ... 33837&fromuid=38665）一文的基础上修改得到的，总体思路是先调用gaussian算UHF并检验波函数稳定性，之后用MOKIT调用pyscf等程序做多参考计算，然后返回能量及梯度给原始gaussian程序，这样能充分发挥gaussian的几何优化功能和stable=opt的功能，以及MOKIT方便地做多参考的功能，计算结果也能方便地用gv查看。但是我在用这个脚本调用pyscf做多参考，NO分子作测试的时候，发现automr输出的能量和梯度都能正确读取，gaussian的输出文件中也能看到这些值，但是在优化过程中梯度很大的情况下分子结构却几乎不变（图片在最后），感觉像是gaussian的optimizer出了问题，想请教一下大家问题出在哪里？


下面是涉及到的脚本及输入文件（本人小白，可能有写错或不妥的地方，请大佬们指正）
脚本gp.py:（与gmm.py的主要区别：1. gaussian输入文件多了G_BASIS，方便混合基组的输入。2. 把生成molpro的输入文件变成了生成automr的输入文件。由于automr从fch文件中读取坐标等信息，因此输入模板中随便写就可以。3. 针对automr程序输出文件的格式，对能量和梯度的读取进行了修改。4. 这个脚本仅针对有解析梯度的方法）

1. #!/path/to/python python3
   
2. import sys
   
3. import os
   
4. import json
   
5. import re
   
6. from xml.dom.minidom import parse
   
7. 
   
8. # parse settings file
   
9. with open("gp.json") as fo:
   
10.     settings = json.load(fo)
    
11. 
    
12. G_NPROC    = settings["Gaussian"]["nproc"]
    
13. G_MEMORY   = settings["Gaussian"]["memory"]
    
14. G_LEVEL    = settings["Gaussian"]["level"]
    
15. G_BASIS    = settings["Gaussian"]["basis"]
    
16. 
    
17. M_NPROC    = settings["mokit"]["nproc"]
    
18. M_MEMORY   = settings["mokit"]["memory"]
    
19. M_LEVEL    = settings["mokit"]["level"]
    
20. M_COMMANDS = settings["mokit"]["commands"]
    
21. M_TOFIND   = settings["mokit"]["tofind"]
    
22. 
    
23. if not os.path.exists("tmp"):
    
24.     os.mkdir("tmp")
    
25. 
    
26. os.chdir("tmp")
    
27. 
    
28. # Gaussian UHF calculation input file templete
    
29. MOKIT_INP = f"""%nproc={G_NPROC}
    
30. %mem={G_MEMORY}
    
31. %chk=mokit.chk
    
32. {G_LEVEL}
    
33. 
    
34. MOKIT TASK
    
35. 
    
36. """
    
37. 
    
38. ANG2BOHR = 1.8897259886
    
39. 
    
40. ELEMENTS = [
    
41.     "Bq",
    
42.     "H ",                                                                                                                                                                                     "He",
    
43.     "Li", "Be",                                                                                                                                                 "B ", "C ", "N ", "O ", "F ", "Ne",
    
44.     "Na", "Mg",                                                                                                                                                 "Al", "Si", "P ", "S ", "Cl", "Ar",
    
45.     "K ", "Ca", "Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",                                                                                     "Ga", "Ge", "As", "Se", "Br", "Kr",
    
46.     "Rb", "Sr", "Y ", "Zr", "Nb", "Mo", "Te", "Ru", "Rh", "Pd", "Ag", "Cd",                                                                                     "In", "Sn", "Sb", "Te", "I ", "Xe",
    
47.     "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",
    
48.     "Fr", "Ra", "Ac", "Th", "Pa", "U ", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og",
    
49. ]
    
50. 
    
51. (layer, InputFile, OutputFile, MsgFile, FChkFile, MatElFile) = sys.argv[1:]
    
52. 
    
53. atoms = []
    
54. coords = []
    
55. with open(InputFile, "r") as fi:
    
56.     (atoms, derivs, charge, spin) = [int(x) for x in fi.readline().split()]
    
57.     MOKIT_INP += f"{charge} {spin}\n"
    
58.     for i in range(0, atoms):
    
59.         arr = fi.readline().split()
    
60.         atom = ELEMENTS[int(arr[0])]
    
61.         coord = [f"{(float(x) / ANG2BOHR)}" for x in arr[1:4]]
    
62.         MOKIT_INP += f" {atom}"
    
63.         MOKIT_INP += " " * 4
    
64.         MOKIT_INP += (" " * 4).join(coord)
    
65.         MOKIT_INP += "\n"
    
66. 
    
67.     MOKIT_INP += "\n"
    
68. 
    
69. MOKIT_INP = f"""{G_BASIS}
    
70. 
    
71. 
    
72. """
    
73. 
    
74. with open("MOKIT.gjf", "w") as fo:
    
75.     fo.write(MOKIT_INP)
    
76. 
    
77. print(">>> Starting UHF Calculation...")
    
78. os.system("g16 MOKIT.gjf")
    
79. os.system("formchk MOKIT.chk")
    
80. os.system("mv MOKIT.fchk MOKIT.fch")
    
81. print(">>> UHF Calculation Done!")
    
82. 
    
83. print(">>> Preparing automr input file...")
    
84. NEVPT2_INP = f"""%nproc={M_NPROC}
    
85. %mem={M_MEMORY}
    
86. {M_LEVEL}
    
87. 
    
88. mokit{{ {M_COMMANDS} }}
    
89. 
    
90. 0 1
    
91. C 0.0 0.0 0.0
    
92. 
    
93. """
    
94. with open("NEVPT2.gjf", "w") as fo:
    
95.     fo.write(NEVPT2_INP)
    
96. 
    
97. # call MOKIT to do automr calculation
    
98. print(">>> Starting automr Calculation...")
    
99. os.system(f"automr NEVPT2.gjf >& NEVPT2.out")
    
100. print(">>> automr Calculation Done!")
     
101. 
     

复制代码

# extract automr output file
print(">>> Extracting automr output file ...")
# energy part
with open("NEVPT2.out","r") as f:
    lines = f.readlines()
for line in lines:
    if M_TOFIND in line:
        match = re.search(r'[-+]\d+\.\d+',line)
        if match:
           energy = float(match.group())

with open(OutputFile, "w") as fo:
    fo.writelines(f"{energy:20.12E}{0:20.12E}{0:20.12E}{0:20.12E}\n")

# gradients part
with open("NEVPT2.out","r") as fi, open(OutputFile, "a") as fo:
    pattern = re.compile(r'Cartesian gradient \(HARTREE\/BOHR\):.*')
    lines = fi.readlines()
    for i, line in enumerate(lines):
        if pattern.match(line):
            break
    lines = lines[i+1:]
    gradients = []
    for line in lines:
        matches = re.findall(r'[-+]?\d*\.\d+E[-+]\d+', line)
        if not matches:
            break
        gradients += [float(x) for x in matches]
    for i in range(0, len(gradients), 3):
        fo.write(f"{float(gradients):20.12E}{float(gradients[i+1]):20.12E}{float(gradients[i+2]):20.12E}\n")

（后面这段总是blockquote掉，干脆就不放代码框里了）

gp.json:(这一块主要区别在于1.gaussian的基组统一采用gen或genecp。2.满足automr输入文件格式。3. 针对automr中能量输出的关键词修改tofind，如SC-NEVPT2的tofind就应为"tofind": "E(SC-NEVPT2) ="。 )

1. {
   
2.     "Gaussian": {
   
3.         "nproc": "8",
   
4.         "memory": "32GB",
   
5.         "level": "#p uhf gen stable=opt scf=(noincfock,novaracc,xqc) nosymm int=(nobasistransform)",
   
6.         "basis":[
   
7.         "n o",
   
8.         "6-31g*",
   
9.         "****",
   
10.         "",
    
11.         ""
    
12.         ]
    
13.     },
    
14.     "mokit": {
    
15.         "nproc": "8",
    
16.         "memory": "32GB",
    
17.         "level":"casscf(5,5)/gen",
    
18.         "commands": "casscf_prog=pyscf,force,ist=2,readuhf="mokit.fch"",
    
19.         "tofind": "E(CASSCF) ="
    
20.     }
    
21. }
    

复制代码

输入文件：

1. %nprocshared=1
   
2. %chk=no.chk
   
3. #p opt=(nomicro) external='python -u gp.py'
   
4. 
   
5. t
   
6. 
   
7. 0 2
   
8. n 0.0 0.0 0.0
   
9. o 0.0 0.0 1.3
   
10. 
    
11. --link1--
    
12. %nprocshared=1
    
13. %chk=no.chk
    
14. #p freq=num geom=allcheck external='python -u gp.py'
    
15. 
    
16. 
    

复制代码

mizu-bai

gmm.py 这随便糊的玩意居然还有人能在上面二次开发，真的没想到
以及以前确实想过改改接 automr，出于懒一直没写

shenzp

wang7344412 发表于 2023-5-3 12:10
谢谢楼主，我还有几个问题想问一下。Q1：NEVPT2的结构优化实现了么？Q2：CASSCF 由于缺乏动态相关，结构优 ...

1. 后来发现数值梯度又行了（但是我印象中没改脚本），保险起见，附件中的脚本又更新了一下，同时加入了NEVPT2数值梯度的例子。不过目前只测试了NO，更复杂的体系应该也可以，但是会很昂贵，而且gaussian的数值梯度的变量个数是有上限的。
2. 具体差多少我也不好说，可以看看相关的benchmark，应该也和体系有关。如果体系不大，计算资源又比较充足的话，肯定用NEVPT2好（主要是体系稍大的话NEVPT2就根本算不动，用CASSCF优化也是无奈之举）。关于计算级别的选择的讨论可以参考MOKIT手册https://jeanwsr.gitlab.io/mokit-doc-mdbook/chap3_quick.html。
3. 如果用MOKIT调用的话，即使用了NEVPT2和force，算的也是CASSCF的梯度。支持NEVPT2解析梯度的程序我不了解，这方面还是建议咨询邹神
4. 只是测试的时候这样写了，让程序自己确定也是完全可以的。

wang7344412

谢谢楼主，我还有几个问题想问一下。Q1：NEVPT2的结构优化实现了么？Q2：CASSCF 由于缺乏动态相关，结构优化结果与NEVPT2相比误差大么？Q3:Pyscf 2.1.1 是不是已经支持NEVPT2的解析梯度了？还有其他软件支持吗？Q4：Mokit 不是会自动寻找合适的活性空间么，这里还需要手动确定活性空间么？