计算化学公社

 找回密码 Forget password
 注册 Register
Views: 5750|回复 Reply: 21

[辅助/分析程序] KST48:更强大的MECP程序

  [复制链接 Copy URL]

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

发表于 Post on 2022-1-9 16:27:04 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 Accelerator 于 2022-12-6 12:28 编辑

最近的研究工作涉及到很多势能面交叉点的计算,虽然多参考态还是得做,但用DFT找MECP要方便得多。Harvey原版的MECP程序功能十分有限,虽然套上社长的接口或easyMECP之类后易用性有所提高但还是太薄弱了,于是索性扔掉Harvey的轮子,花一天时间自己写了一个纯Python的交叉点寻找程序。支持同Gaussian和ORCA联用(写代码时考虑了xtb,不知道有没有人会用xtb找大体系的交叉点,有需求的话之后再补充)。支持不同自旋态之间或不同激发态之间(结合ORCA还可以用SF-TDDFT做)的交叉点搜索。

logo_KST48.png


使用方法:



程序只有一个Python文件kst48.py,使用时在当前目录下放置一个JOBS目录用于存储输入输出文件,并准备一个KST48的输入文件。一个示例如下:

  1. #This subset is required. It controls your quantum chemistry tasks.
  2. nprocs = 8 #processors
  3. mem = 8GB # memory to be used. change this into the maxcore value if you want to use ORCA
  4. method = wb97xd def2sv # your keywords line. It will be presented in the job files. Don't write guess=mix or stable=opt; they will be added automatically.
  5. td1 =  # keywords for TD-DFT of state A (only for Gaussian; please write it in the tail part for ORCA)
  6. td2 = # keywords for TD-DFT of state B (only for Gaussian; please write it in the tail part for ORCA)
  7. mp2 = false #set true for MP2 or doubly hybrid calculation in Gaussian
  8. charge = 0
  9. mult1 = 1 # multiplicity of state A
  10. mult2 = 3 # multiplicity of state B
  11. mode = stable #normal; stable; read; inter_read; noread

  12. #This subset is optional. It controls the convergence threshols, and the details of the GDIIS algorithm. Shown here are the default values.
  13. dE_thresh = 0.000050
  14. rms_thresh = 0.0025
  15. max_dis_thresh = 0.004
  16. max_g_thresh = 0.0007
  17. rms_g_thresh = 0.0005
  18. max_steps = 100
  19. max_step_size = 0.1
  20. reduced_factor = 0.5 # the gdiis stepsize will be reduced by this factor when rms_gradient is close to converge

  21. # This subset controls which program you are using, and how to call them
  22. program = gaussian  #gaussian or orca
  23. gau_comm = 'g16'
  24. orca_comm = '/opt/orca5/orca'
  25. xtb_comm = 'xtb'

  26. #Between *geom and *, write the cartesian coordinate of your initial geometry (in angstrom)
  27. *geom
  28. C                 -0.03266394   -2.07149616    0.00000000
  29. H                  0.76176530   -1.26928217    0.00000000
  30. H                 -0.91673236   -1.36926867    0.00000000
  31. *
  32. #If you have anything to be put at the end of the input file, write it here. This part is especially useful for ORCA: you can add anything related to your keywords here.
  33. *tail1
  34. #-C -H 0
  35. #6-31G(d)
  36. #****
  37. *
  38. *tail2
  39. *

  40. #This subset controls the constraints. R 1 2 1.0 means to fix distance between atom 1 and 2 (start from 1) to be 1.0 angstrom.
  41. *constr
  42. #R 1 2 1.0
  43. #A 1 2 3 100.0 # to fix angle 1-2-3 to be 100 degree
  44. #S R 1 2 1.0 10 0.1 # to run a scan of R(1,2) starting from 1.0 angstrom, with 10 steps of 0.1 angstrom
  45. #S R 2 3 1.5 10 0.1 # you can at most set a 2D-scan
  46. *
复制代码
必须要设置的是输入文件的第一部分,包括使用计算资源、电荷、两个态的多重度、关键字等(不要在这里写stable=opt, force, guess=read等,这些由程序自动添加)。支持5种模式,在下文讲。


第二部分可选,不写就用默认值。设置构型优化的细节。程序最开始3步用BFGS算法,后续步骤切换到GDIIS算法,基于3个历史几何结构进行外推,并使用这三者的平均BFGS矩阵作为近似Hessian。当RMS梯度降低到收敛限的10倍时每一步的步长将自动乘以reduced_factor以避免因步长过大导致的震荡。最大步长默认是0.1 angstrom,这些参数都可以改。


第三部分控制使用何种程序,第四部分是笛卡尔坐标,非常简单。第五部分tail中的内容会被写入输入文件末尾,例如Gaussian的自定义基组等。这对ORCA特别有用,ORCA中以%开头的关键字都可以放在这里。A和B两个态分别用tail1和tail2控制,在设计这一点时主要是为了给SF-TDDFT用的。


最后一部分设置几何约束。可以固定键长和键角,也可以在3N-7维空间内进行一维或二维扫描。


输入文件写好后,放在kst48.py和JOBS相同目录下,假如叫做inp,运行python3 kst48.py inp即可,信息会输出在屏幕上。


关于运行模式

程序支持5种运行模式:

1. normal,类似Harvey的MECP。

2. stable,首先基于初始构型自动进行一个stable=opt,再读取得到的波函数进行后续。注意如果用ORCA处理单重态的情况,记得在method里写UKS。

3. read, 如果已经有了正确的波函数文件,就放在和kst.py相同目录下,分别叫做a.chk和b.chk(或.gbw),设置为read模式后就会自动读取,其他类似stable模式。

4. inter_read:某些十分困难的场合,RKS直接做stable=opt会告诉说是稳定的,但实际上存在能量低得多的另一个稳定的UKS波函数,它需要读取三重态波函数并且写guess=mix甚至scf=vshift=-1000之类来加强轨道混合才能得到。使用inter_read模式,程序类似stable模式,但先算态B,然后对初始构型下的态A进行稳定性分析时读取态B的收敛波函数。用户应该自己写guess=mix等其他收敛控制关键字。

5. noread:不读取任何轨道,包括优化过程中也是每一步都重新生成轨道初猜。


示例输出信息(上述inp文件对应的亚甲基卡宾,单重态和三重态的MECP):

  1. ****KST48 PROGRAM: a powerful tool for MECP locating****
  2. ****By Yumiao Ma, BSJ Institute, 2022/01/09****
复制代码
收敛完成后,会把最终结构对应的输入输出文件复制回当前目录。


下载和引用:

使用KST48必须进行引用,包括其Github页面(https://github.com/RimoAccelerator/KST48)以及作者使用KST48的第一篇文章。
1. Y. Ma, A. A. Hussein, ChemistrySelect 2022, 7, e202202354.
2. Yumiao Ma.  KST48: A Powerful Tool for MECP location. https://github.com/RimoAccelerator/KST48, accessed on xxxx.xx.xx


关于该程序目前只对少量分子进行了测试,欢迎大家试用,发现任何bug欢迎留言。(后更新:目前已经对大量体系进行了测试,非常好用。)

      **** 2022年1月24日更新:近日来连续进行了多次更新,原本上传在附件的版本很容易造成混乱,于是删除了帖子里的附件,请用户统一到Github上下载最新的发布版本。在目前的Github中,除了程序本体和输入文件模板外,还增加了一个能够与KST48配套,在常见架构的集群队列系统中调用ORCA的脚本。
      **** 2022年2月6日更新:已增加对BAGEL的支持,可搜索任意两个能用BAGEL计算的态(包含不同自旋多重度)之间的交叉点。使用方法见Github页面。





评分 Rate

参与人数
Participants 14
威望 +1 eV +61 收起 理由
Reason
zsu007 + 5 赞!
WaterMirror + 5 GJ!
ljc050512 + 5 谢谢
hebrewsnabla + 3 GJ!
Satoru + 5 好物!
Novice + 5 好物!
hdhxx123 + 5 好物!
sobereva + 1
ChrisZheng + 5 GJ!
ggdh + 5 tql!
让你变成回忆 + 5 GJ!
zjxitcc + 3
ene + 5
wzkchem5 + 5

查看全部评分 View all ratings

6672

帖子

0

威望

4002

eV
积分
10674

Level 6 (一方通行)

发表于 Post on 2022-1-9 17:36:03 | 显示全部楼层 Show all
把level shift设成负值来得到其他SCF解,是一个挺有意思的方法。。。
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-1-9 17:51:17 | 显示全部楼层 Show all
本帖最后由 Accelerator 于 2022-1-9 17:53 编辑
wzkchem5 发表于 2022-1-9 17:36
把level shift设成负值来得到其他SCF解,是一个挺有意思的方法。。。

跑AIMD的时候需要一个能一下子就收敛到正确波函数的方法(而且不能搞guess=read,一些情况下发现前一步结构的稳定波函数传给后一步肯定会得到错误的结果),因此想到的。
但实际上效果并不好,vshift写得很负时虽然能得到UKS波函数,但经常也并不对,还是得搭配stable=opt使用。

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-1-9 20:56:03 | 显示全部楼层 Show all
发现下午上传的版本有个低级错误 之前在代码里把测试用的CH2原子列表写死了 导致读取其他分子会出错。将开头部分list_element 改成[]就好了。

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-1-15 00:05:28 | 显示全部楼层 Show all
最近程序经过了若干修改,已经投稿到了JOSS上。目前测试发现在不施加限制的情况下其收敛速度很高,放一张它和Harvey程序收敛过程的对比:
Figure2.png

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-1-24 19:31:54 | 显示全部楼层 Show all
程序的介绍性文章投到JOSS上后居然由于代码量太少被拒了(投稿时的版本大约800行,JOSS中送审的程序中行数较少的大约有1000行,所以可能1000行是投稿门槛)。因此KST48就不打算发表介绍性文章了,将来等到使用KST48的研究工作发表后将该文章作为标准引文好了。
最近在用多参考态找交叉点,BAGEL运算速度高但优化算法很差,而且不支持找不同自旋态的交叉点;而OpenMolcas优化算法不错但效率太低,于是决定在下一个版本加入对BAGEL的支持。

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-1-24 19:31:56 | 显示全部楼层 Show all
程序的介绍性文章投到JOSS上后居然由于代码量太少被拒了(投稿时的版本大约800行,JOSS中送审的程序中行数较少的大约有1000行,所以可能1000行是投稿门槛)。因此KST48就不打算发表介绍性文章了,将来等到使用KST48的研究工作发表后将该文章作为标准引文好了。
最近在用多参考态找交叉点,BAGEL运算速度高但优化算法很差,而且不支持找不同自旋态的交叉点;而OpenMolcas优化算法不错但效率太低,于是决定在下一个版本加入对BAGEL的支持。

1061

帖子

16

威望

5780

eV
积分
7161

Level 6 (一方通行)

計算化学の社畜

发表于 Post on 2022-1-24 20:34:44 | 显示全部楼层 Show all
乍一看我还以为是AKB48……
Stand on the shoulders of giants

6672

帖子

0

威望

4002

eV
积分
10674

Level 6 (一方通行)

发表于 Post on 2022-1-24 20:36:22 | 显示全部楼层 Show all
Accelerator 发表于 2022-1-24 12:31
程序的介绍性文章投到JOSS上后居然由于代码量太少被拒了(投稿时的版本大约800行,JOSS中送审的程序中行数 ...

这还不简单,多加几行注释
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

604

帖子

1

威望

4887

eV
积分
5511

Level 6 (一方通行)

发表于 Post on 2022-1-24 20:59:49 | 显示全部楼层 Show all
冰释之川 发表于 2022-1-24 20:34
乍一看我还以为是AKB48……

可能是 kousaten48

评分 Rate

参与人数
Participants 1
eV +3 收起 理由
Reason
thanhtam + 3 我很赞同

查看全部评分 View all ratings

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-1-24 22:31:45 | 显示全部楼层 Show all

88

帖子

1

威望

1935

eV
积分
2043

Level 5 (御坂)

发表于 Post on 2022-2-16 10:41:37 | 显示全部楼层 Show all
厉害

8

帖子

0

威望

445

eV
积分
453

Level 3 能力者

发表于 Post on 2022-5-18 15:23:05 | 显示全部楼层 Show all
请问能够支持qm/mm模型吗,或者以后会有这个功能吗,想找到qm/mm模型中qm的MECP

8

帖子

0

威望

445

eV
积分
453

Level 3 能力者

发表于 Post on 2022-5-18 15:25:19 | 显示全部楼层 Show all
人心向背 发表于 2022-5-18 15:23
请问能够支持qm/mm模型吗,或者以后会有这个功能吗,想找到qm/mm模型中qm的MECP

这是我们搭建的qm/mm模型

1918803_2.gjf

173.3 KB, 下载次数 Times of downloads: 6

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

 楼主 Author| 发表于 Post on 2022-5-18 16:41:12 | 显示全部楼层 Show all
人心向背 发表于 2022-5-18 15:23
请问能够支持qm/mm模型吗,或者以后会有这个功能吗,想找到qm/mm模型中qm的MECP

我没用过Gaussian的QM/MM,需要研究。最近不一定有时间,如果要新增功能的话应该要下半年了。

本版积分规则 Credits rule

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

GMT+8, 2023-2-2 23:47 , Processed in 0.708416 second(s), 25 queries .

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