程序只有一个Python文件kst48.py,使用时在当前目录下放置一个JOBS目录用于存储输入输出文件,并准备一个KST48的输入文件。一个示例如下:
- #This subset is required. It controls your quantum chemistry tasks.
- nprocs = 8 #processors
- mem = 8GB # memory to be used. change this into the maxcore value if you want to use ORCA
- 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.
- td1 = # keywords for TD-DFT of state A (only for Gaussian; please write it in the tail part for ORCA)
- td2 = # keywords for TD-DFT of state B (only for Gaussian; please write it in the tail part for ORCA)
- mp2 = false #set true for MP2 or doubly hybrid calculation in Gaussian
- charge = 0
- mult1 = 1 # multiplicity of state A
- mult2 = 3 # multiplicity of state B
- mode = stable #normal; stable; read; inter_read; noread
- #This subset is optional. It controls the convergence threshols, and the details of the GDIIS algorithm. Shown here are the default values.
- dE_thresh = 0.000050
- rms_thresh = 0.0025
- max_dis_thresh = 0.004
- max_g_thresh = 0.0007
- rms_g_thresh = 0.0005
- max_steps = 100
- max_step_size = 0.1
- reduced_factor = 0.5 # the gdiis stepsize will be reduced by this factor when rms_gradient is close to converge
- # This subset controls which program you are using, and how to call them
- program = gaussian #gaussian or orca
- gau_comm = 'g16'
- orca_comm = '/opt/orca5/orca'
- xtb_comm = 'xtb'
- #Between *geom and *, write the cartesian coordinate of your initial geometry (in angstrom)
- *geom
- C -0.03266394 -2.07149616 0.00000000
- H 0.76176530 -1.26928217 0.00000000
- H -0.91673236 -1.36926867 0.00000000
- *
- #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.
- *tail1
- #-C -H 0
- #6-31G(d)
- #****
- *
- *tail2
- *
- #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.
- *constr
- #R 1 2 1.0
- #A 1 2 3 100.0 # to fix angle 1-2-3 to be 100 degree
- #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
- #S R 2 3 1.5 10 0.1 # you can at most set a 2D-scan
- *
复制代码必须要设置的是输入文件的第一部分,包括使用计算资源、电荷、两个态的多重度、关键字等(不要在这里写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):
- ****KST48 PROGRAM: a powerful tool for MECP locating****
- ****By Yumiao Ma, BSJ Institute, 2022/01/09****
复制代码收敛完成后,会把最终结构对应的输入输出文件复制回当前目录。
下载和引用:
1. Y. Ma, A. A. Hussein, ChemistrySelect 2022, 7, e202202354.
关于该程序目前只对少量分子进行了测试,欢迎大家试用,发现任何bug欢迎留言。(后更新:目前已经对大量体系进行了测试,非常好用。)
**** 2022年1月24日更新:近日来连续进行了多次更新,原本上传在附件的版本很容易造成混乱,于是删除了帖子里的附件,请用户统一到Github上下载最新的发布版本。在目前的Github中,除了程序本体和输入文件模板外,还增加了一个能够与KST48配套,在常见架构的集群队列系统中调用ORCA的脚本。
**** 2022年2月6日更新:已增加对BAGEL的支持,可搜索任意两个能用BAGEL计算的态(包含不同自旋多重度)之间的交叉点。使用方法见Github页面。