计算化学公社

 找回密码 Forget password
 注册 Register
Views: 16191|回复 Reply: 44
打印 Print 上一主题 Last thread 下一主题 Next thread

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

  [复制链接 Copy URL]

482

帖子

10

威望

6936

eV
积分
7618

Level 6 (一方通行)

BSJ Institute

本帖最后由 Accelerator 于 2024-11-6 22:02 编辑

注:KST48长期以来不断更新,添加了许多新功能,本文已明显过时。Github主页上有全面的使用方法介绍,并对多种典型应用场合给出了案例,请读者参考Github主页上的最新用法。

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


使用方法:



程序只有一个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 16
威望 +1 eV +70 收起 理由
Reason
shuiningzhu + 5 好物!
Kmetsch + 4 好物!
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

482

帖子

10

威望

6936

eV
积分
7618

Level 6 (一方通行)

BSJ Institute

45#
 楼主 Author| 发表于 Post on 2024-11-6 22:05:50 | 只看该作者 Only view this author
本帖最后由 Accelerator 于 2024-11-6 22:06 编辑
meng0612 发表于 2024-11-6 21:33
请问老师,目前能算S0和S1的势能面交叉吗

Github主页有案例教程,里面直接给出了2个相关的例子:
S1-T1 MECP of a TADF molecule (如果你使用TD-DFT计算,可采用书写td关键字的方法)
ORCA CASSCF: S0/S1 crossing point of CH2 (如果你使用ORCA进行CASSCF计算,可采用书写tail1和tail2关键字的方法)

216

帖子

4

威望

1498

eV
积分
1794

Level 5 (御坂)

44#
发表于 Post on 2024-11-6 21:47:47 | 只看该作者 Only view this author
Accelerator 发表于 2022-1-24 19:31
程序的介绍性文章投到JOSS上后居然由于代码量太少被拒了(投稿时的版本大约800行,JOSS中送审的程序中行数 ...

JOSS的人真的不懂,我本来打算不想搞太麻烦,然后直接上JOSS,JOSS直接拒稿,然后索性又鼓捣三个月补充了很多功能和测试,现在在JCTC送审

7

帖子

0

威望

404

eV
积分
411

Level 3 能力者

43#
发表于 Post on 2024-11-6 21:33:02 | 只看该作者 Only view this author
请问老师,目前能算S0和S1的势能面交叉吗

44

帖子

0

威望

3400

eV
积分
3444

Level 5 (御坂)

42#
发表于 Post on 2024-5-4 09:25:30 | 只看该作者 Only view this author
赞!很好用的程序。不知作者能否考虑添加像sobMECP中将优化过程几何结构输出到traj.xyz的功能,这个对于判断结构变化有时还是很有用的

11

帖子

0

威望

537

eV
积分
548

Level 4 (黑子)

41#
发表于 Post on 2023-11-1 21:09:53 | 只看该作者 Only view this author
请问目前支持圆锥交叉点的搜索吗

482

帖子

10

威望

6936

eV
积分
7618

Level 6 (一方通行)

BSJ Institute

40#
 楼主 Author| 发表于 Post on 2023-10-20 12:58:32 | 只看该作者 Only view this author
搁置了许久,终于将这段时间内大家积攒的问题和诉求进行了一次更新。
本次更新内容:
新增功能
1. 原子固定:通过在输入文件内设置形如fixedAtoms = 1, 2-3的选项,可以实现原子固定,具体方法为优化过程中直接设置其梯度为0.编号从1开始,形如2-3的定义中包含原子3.
2. 对Gaussian的ONIOM计算的支持。需要设置如下选项:
  1. isONIOM = true # if you are using the ONIOM feature of Gaussian, set this true.
  2. chargeAndMultForONIOM1 = 0 1 0 1 0 1 0 1 0 1 0 1 #only useful for ONIOM calculation
  3. chargeAndMultForONIOM2 = 0 1 0 1 0 1 0 3 0 3 0 3 #only useful for ONIOM calculation
复制代码

此时先前的charge, mult1, mult2等可忽略。随后在geom部分定义分层信息,形如:
  1. C        2.02254521   -0.36259181   -0.01365118 M H 8
  2. H        2.37786899   -1.37187093   -0.01446684 M
  3. H        2.37789587    0.14134325   -0.88810836 M
  4. C        0.48255054   -0.36054374   -0.01015114 L H 1
  5. H        0.12722676    0.64873538   -0.00933548 L
  6. H        0.12322871   -0.86446478   -0.88299213 L
  7. H        0.12719988   -0.86447880    0.86430604 L
  8. N        2.51619296    0.32971093    1.18548551 H
  9. O        2.72397867   -0.30706606    2.22615545 H
  10. O        2.72401165    1.54934090    1.15437446 H
复制代码

注意原子符号后面要直接跟着X坐标。
3. 支持外链几何坐标文件
通过形如
  1. *geom
  2. @a.log
  3. *
复制代码

的方法,可以自动读取a.log中的最后一个几何坐标。支持Gaussian的.log文件、.gjf文件,以及.xyz文件(支持xyz轨迹,自动读取最后一个结构)。
4. Bug修复:修复了CASPT2读取能量不正确的问题。

3

帖子

0

威望

53

eV
积分
56

Level 2 能力者

39#
发表于 Post on 2023-9-19 19:26:03 | 只看该作者 Only view this author
wzkchem5 发表于 2023-9-18 19:39
如果自己解决了自己的问题,最好把解决方案也发一下,这样后人遇到类似问题就可以直接搜到这个帖子,不用 ...

我输入文件写错了,不是别的原因,这种也需要发一下吗

1万

帖子

0

威望

9868

eV
积分
22108

Level 6 (一方通行)

38#
发表于 Post on 2023-9-18 19:39:26 | 只看该作者 Only view this author
jwx-yao 发表于 2023-9-18 12:26
好的好的,谢谢谢谢,又试了试,解决问题了,感谢老师回答

如果自己解决了自己的问题,最好把解决方案也发一下,这样后人遇到类似问题就可以直接搜到这个帖子,不用再问一遍
Zikuan Wang
山东大学光学高等研究中心 研究员
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员
Google Scholar: https://scholar.google.com/citations?user=XW6C6eQAAAAJ
ORCID: https://orcid.org/0000-0002-4540-8734
主页:http://www.qitcs.qd.sdu.edu.cn/info/1133/1776.htm
GitHub:https://github.com/wzkchem5
本团队长期招收研究生,有意者可私信联系

3

帖子

0

威望

53

eV
积分
56

Level 2 能力者

37#
发表于 Post on 2023-9-18 19:26:24 | 只看该作者 Only view this author
wzkchem5 发表于 2023-9-18 17:15
但凡问报错问题,必须上传完整的输入输出文件
可能不是你orca输入文件写的问题,而是你的环境配置问题, ...

好的好的,谢谢谢谢,又试了试,解决问题了,感谢老师回答

1万

帖子

0

威望

9868

eV
积分
22108

Level 6 (一方通行)

36#
发表于 Post on 2023-9-18 17:15:17 | 只看该作者 Only view this author
jwx-yao 发表于 2023-9-18 09:10
请问老师可以发一个orca计算的实例吗,我想用orca算老是报错...

但凡问报错问题,必须上传完整的输入输出文件
可能不是你orca输入文件写的问题,而是你的环境配置问题,那样的话别人给orca计算的实例也没用,必须你给出你这里具体报的是什么错才行
Zikuan Wang
山东大学光学高等研究中心 研究员
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员
Google Scholar: https://scholar.google.com/citations?user=XW6C6eQAAAAJ
ORCID: https://orcid.org/0000-0002-4540-8734
主页:http://www.qitcs.qd.sdu.edu.cn/info/1133/1776.htm
GitHub:https://github.com/wzkchem5
本团队长期招收研究生,有意者可私信联系

3

帖子

0

威望

53

eV
积分
56

Level 2 能力者

35#
发表于 Post on 2023-9-18 16:10:08 | 只看该作者 Only view this author
请问老师可以发一个orca计算的实例吗,我想用orca算老是报错...

482

帖子

10

威望

6936

eV
积分
7618

Level 6 (一方通行)

BSJ Institute

34#
 楼主 Author| 发表于 Post on 2023-8-16 11:55:30 | 只看该作者 Only view this author
yxdd98 发表于 2023-8-15 20:49
您好,请问KST48有约束原子坐标的功能吗?我的体系是簇模型,需要冻结边上的原子,我试着用了orca里的const ...

目前不能固定原子位置。我尽量在9月推出一个新版本,加入原子坐标冻结的功能,以及先前楼上提到的CASPT2读取也一直没改,在这个版本内一同加上。

评分 Rate

参与人数
Participants 1
eV +2 收起 理由
Reason
yxdd98 + 2 赞!

查看全部评分 View all ratings

87

帖子

0

威望

2002

eV
积分
2089

Level 5 (御坂)

33#
发表于 Post on 2023-8-15 20:49:18 | 只看该作者 Only view this author
您好,请问KST48有约束原子坐标的功能吗?我的体系是簇模型,需要冻结边上的原子,我试着用了orca里的constraints字段,但发现计算过程中本该冻结的原子位置还是会变化。

13

帖子

0

威望

633

eV
积分
646

Level 4 (黑子)

32#
发表于 Post on 2023-5-17 20:18:24 | 只看该作者 Only view this author
jmliu 发表于 2023-5-16 19:19
请问老师能发一个bagel的输入案例吗?我按照git下来的那个文件写的,试了好多次都没法运行。一直提示我的输 ...

问题已经解决了,期待大佬的新增caspt2的优化功能。

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

GMT+8, 2026-2-22 03:14 , Processed in 0.184058 second(s), 25 queries , Gzip On.

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