KST48：更强大的MECP程序

人心向背

Accelerator 发表于 2022-5-18 16:41
我没用过Gaussian的QM/MM，需要研究。最近不一定有时间，如果要新增功能的话应该要下半年了。

好的，谢谢大佬，期待大佬的新增功能。

Accelerator

今日更新：增加了一个结构插值功能（LST）
与过渡态搜索一样，有些交叉点非常难找，特别是如果同时涉及到许多化学键的改变时；而这在作者所研究的电子转移与化学反应的耦合中非常常见。如果只涉及一两根化学键改变，通过键长扫描可以比较方便地处理，而一旦化学键太多，此法就行不通了。于是有这样一种思路：在两个态的极小值点之间进行插值，对于插入的中间结构，计算两个态的能量，取能量差最小的作为交叉点的初始结构。于是加入了此功能，通过在输入文件中设置LST1和LST2两部分，程序会对两个分子首先进行对齐，然后插值并生成一系列gjf文件。
事实证明这个功能并不好用，简单场合还行，如果几何结构变动很大，插值出来的点结构就很不合理。虽然原理上也可以借助QST2或CI-NEB的思想进一步处理，但考虑到交叉点计算的单步耗时往往很长，这些算法并不实用。因此目前该功能仅停留在实验阶段。

Freeman

Accelerator 发表于 2022-6-16 15:35
今日更新：增加了一个结构插值功能（LST）
与过渡态搜索一样，有些交叉点非常难找，特别是如果同时涉及到 ...

如果势能面是左边的情况，插值是没问题的，但如果是右边的情况，插值就没用了。所以找过渡态的技巧大概不能直接套用在mecp上。

CopperL

请问这个报错是为什么呢？是因为体系太大了吗？我用小分子体系比如乙烯、亚甲基卡宾等同样的输入只修改原子坐标测试是没问题的。（报错如下所示，输入文件见附件）
Now Entering BFGS Step 0
Traceback (most recent call last):
  File "kst48.py", line 1076, in <module>
    main()
  File "kst48.py", line 1021, in main
    conv_step, GEOM = runOpt(GEOM, flag='hybrid', state1 = state1, state2 = state2, mult1 = mult1, mult2 = mult2)
  File "kst48.py", line 893, in runOpt
    X1, G1, B1, E1, E2 = BFGS(X0, G0, B0, n_step, state1 = state1, state2 = state2, mult1 = mult1, mult2 = mult2)
  File "kst48.py", line 800, in BFGS
    XNew = propagationBFGS(X0, G0, B0)
  File "kst48.py", line 687, in propagationBFGS
    dk = -numpy.linalg.solve(Bk, Gk)
  File "<__array_function__ internals>", line 6, in solve
  File "/usr/local/lib64/python3.6/site-packages/numpy/linalg/linalg.py", line 394, in solve
    r = gufunc(a, b, signature=signature, extobj=extobj)
ValueError: solve1: Input operand 1 has a mismatch in its core dimension 0, withgufunc signature (m,m),(m)->(m) (size 252 is different from 257)

Accelerator

CopperL 发表于 2022-9-24 09:56
请问这个报错是为什么呢？是因为体系太大了吗？我用小分子体系比如乙烯、亚甲基卡宾等同样的输入只修改原子 ...

我用PBE0-D3BJ/STO-3G测试这个例子，在我的机器上表现正常。是否可以将你运行目录下的全部文件打包发上来？

CopperL

Accelerator 发表于 2022-9-25 16:18
我用PBE0-D3BJ/STO-3G测试这个例子，在我的机器上表现正常。是否可以将你运行目录下的全部文件打包发上来 ...

感谢回复。
文件打包后的大小超出上传上限，我把运行目录下的文件上传了某度网盘，链接是

链接：https://pan.baidu.com/s/1AV1lcWUaXpeTJUOFBNmZWw?pwd=pvk5
提取码：pvk5

麻烦您下载查看，谢谢。

Accelerator

CopperL 发表于 2022-9-26 09:00
感谢回复。
文件打包后的大小超出上传上限，我把运行目录下的文件上传了某度网盘，链接是

是由于输出文件0_A.log中几何结构部分的格式与我编写和使用时遇到的不同，导致没有正确读取。只要将kst48.py的314行改为如下即可：

1.             elif 'Distance matrix' in l or 'Rotational constants' in l:

复制代码

Github上的代码已更新。感谢报告这一bug。

ABQTrap

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

老师好，请问KST48可以只做几何约束优化，不用MECP功能吗？因为BAGEL没有几何约束优化的功能。谢谢！

Accelerator

ABQTrap 发表于 2023-3-29 10:46
“最后一部分设置几何约束。可以固定键长和键角，也可以在3N-7维空间内进行一维或二维扫描。” ...

找交叉点的梯度是两个态的梯度进行四则运算得到的，所以原则上可以使得程序只采用一个态的梯度。但是这样需要你自己修改代码，比较折腾。
http://bbs.keinsci.com/thread-27086-1-1.html提到的方法是否可以做限制性优化，我没有试过。

ABQTrap

Accelerator 发表于 2023-3-30 22:11
找交叉点的梯度是两个态的梯度进行四则运算得到的，所以原则上可以使得程序只采用一个态的梯度。但是这样 ...

好的，谢谢老师！

柚子

老师，您好！在此我有两个问题想要请教。①请问出现这样的报错是什么原因呢？不知道是不是我输入文件的原因，但运行例子也是这样的情况。
  File "kst48.py", line 54
    header_1 = f'%chk=a.chk\n%nprocs={NProcs} \n%mem={Mem} \n# {Method} {Td1} nosymm\n\n Title Card \n\n{Charge} {Mult1}'
                                                                                                                        ^
SyntaxError: invalid syntax


②请问这个程序默认搜索的是S1和T1的MECP吗？我本科毕设是TADF方向的，想要研究机理，需要这个结构点的信息。


期待老师的回复。

Accelerator

柚子 发表于 2023-4-4 10:50
老师，您好！在此我有两个问题想要请教。①请问出现这样的报错是什么原因呢？不知道是不是我输入文件的原因 ...

把输入文件发上来。
找哪几个态的交叉点在输入文件里设置，没有默认一说。如果想用Gaussian找S1和T1的交叉点，需要在td1或td2里写跟TD-DFT有关的关键字。

ABQTrap

老师您好，我是刚接触KST48和BAGEL。

我尝试了GAMESS中exam40的MECP例子（CH2），用KST48做MECP，方法是normal，BAGEL用的是CASPT2。我想请问KST48输出的两个态的能量信息默认是CASSCF吗？因为我感觉输出的能量好像不是CASPT2能量。比如几何优化最后一步，KST48在terminal上的输出能量是E1 = -38.89377848，E2 = -38.90046691，BAGEL中A.log的能量 CASPT2 energy:  -38.972166738288，B.log的能量为-38.971255607506。我附上了输出信息，麻烦老师看一下。因为我不是很懂BAGEL，可能是我搞错了。谢谢！

Accelerator

ABQTrap 发表于 2023-4-10 16:34
老师您好，我是刚接触KST48和BAGEL。

我尝试了GAMESS中exam40的MECP例子（CH2），用KST48做MECP，方法是 ...

是的，由于在编写时主要使用CASSCF进行测试，目前读取的是CASSCF能量。本周我会对代码进行更新，以支持CASPT2。

ABQTrap

Accelerator 发表于 2023-4-10 16:50
是的，由于在编写时主要使用CASSCF进行测试，目前读取的是CASSCF能量。本周我会对代码进行更新，以支持CA ...

好的，谢谢老师，太牛了