计算化学公社

标题: 适用于更多程序的sobMECP修改版(9.23更新MECI支持) [打印本页]

作者
Author: wal 时间: 2025-9-14 15:18
标题: 适用于更多程序的sobMECP修改版(9.23更新MECI支持)
本帖最后由 wal 于 2025-9-24 17:58 编辑

6月左右的时候尝试过把sobMECP接上MOKIT，但因为sobMECP迭代部分比较零落，没搞明白。最近遇到了需要MRSF-TDDFT搜索MECP的情况，但GAMESS不原生支持，于是再次把目光瞄向了sobMECP
GitHub：banemecp
主程序： (, 下载次数 Times of downloads: 7) (, 下载次数 Times of downloads: 5)
配置文件与示例输入： (, 下载次数 Times of downloads: 7)
一些算例 (, 下载次数 Times of downloads: 9) (, 下载次数 Times of downloads: 12)

本重构版保留了原程序Harvey算法，其他部分按照我的构思让AI重新设计了一下。现在MECP.f的行为逻辑为：接受astep1.xyz、astep1.grad1、astep1.grad2(gard文件格式与xyz同，注释行写入能量值)文件，返回新xyz文件。同时，允许通过命令行参数改变收敛限、置信半径等。搜索MECP时，将由迭代器来调度MECP.x进行迭代。修改可能会降低运行效率，不过我也没啥量化程序编程经验，能跑起来就不错了XD

迭代器需要两类文件输入。第一个是配置文件，需要在[main]部分指定量子化学运行参数：

[main]
env='source ~/.bashrc
module load g16'
# As a suffix, forms ${ACTUAL_INP}
INPUT_SUFFIX=gjf
# As a suffix, forms ${ACTUAL_OUT}
OUTPUT_SUFFIX=log
RUN_CMD='echo ">>>>> Start g16 calculation"
g16 ${ACTUAL_INP} > ${ACTUAL_OUT}
echo "g16 done. <<<<<<"'

复制代码

关于梯度与能量的提取，为避免每个新程序都要动源代码，本重构版同样使用配置文件管理：

[scf]
E=SCF Done:\s+E$[^)]*$\s*=\s*([-\d\.]+)
GARD.Loacte=Forces (Hartrees/Bohr)
GARD.NLineSkip=2
GARD.TargetColumns=3,4,5
# GARD.EndBy leave blank means eny by blank line.
GARD.EndBy=----
GARD.unit=hartree_per_bohr
GRAD.Type=force

复制代码

例如，上方为适用于Gaussian16普通DFT计算的提取配置。能量会使用E=给出的正则表达式进行提取；梯度部分采用GRAD.locate作为梯度输出定位符，NLineSkip指示需要跳过的行数，TargetColumns指示目标列的位置，EndBy指示终止提取的定位符。如果定位符会在输出文件中出现多次，可以通过GARD.LoacteCount和E.LoacteCount告诉程序在定位符/正则表达式匹配第几次时开始提取。如此一来，本程序可以自然适配绝大多数量子化学程序，为任何有梯度的方法提供MECP接口。更多配置规则可以参考配置文件里的rule.md。

第二类是输入文件，告诉程序你要怎么计算两个梯度：

%Prog g16
%Geom init.xyz
%Control
maxcyc=100
tmpdir=.
keeptmp=true
restart=false
end

复制代码

文件头部分，%Prog指定了conf文件名；%Geom指定了xyz文件名(这两个不需要end结尾)；%Control部分则指定一些运行参数，例如最大迭代圈数、临时文件路径、重启、收敛限等。

%InpTmplt1
%chk=state1.chk
# B3LYP/def2svp force [guess(read)]
state1
0 1
[geom]
end
%InpTmplt2
%chk=state2.chk
# B3LYP/def2svp force [guess(read)]
state2
0 3
[geom]
end

复制代码

%InpTmplt1和%InpTmplt2告诉程序第一个和第二个输入文件怎么写。其中，[]方括号内为占位符，一些可用占位符：[geom]表示几何坐标；[nstep]表示步数；[xyzfile]表示xyz文件名；还有一些参考模板目录的rule.md。不在占位符列表但以方括号圈出的内容视为读取初猜关键词，将在第一轮忽视，第二轮开始启用。

%GrpTmplt
state1=scf
state2=td
end

复制代码

最后是%GrpTmplt，指定能量与梯度的提取配置。如果你两个计算其中一个是TD-DFT，那么你可以写个td的梯度提取配置，用state=td指定使用对应方案得到正确的能量与梯度。

%InpTmplt1
echo ">>>>> Start gamess calculation1"
cp ../gamess.sh .
bash gamess.sh [xyzfile] --iroot 1 --basename astep[nstep]_state1 >> mrsfout.tmp 2>&1
end
%InpTmplt2
echo ">>>>> Start gamess calculation2"
cp ../gamess.sh .
bash gamess.sh [xyzfile] --iroot 2 --basename astep[nstep]_state2 >> mrsfout.tmp 2>&1
cp opt.trj ..
end

复制代码

对于复杂计算，如利用MOKIT的fch2inp程序传Gaussian的RODFT轨道给GAMESS做MRSF-TDDFT，程序提供了external模式来执行脚本。external模式在配置文件不提供RUN_CMD时触发(此时程序仍然需要OUT_SUFFIX来了解用什么后缀的文件提取能量与梯度)，此时会将%InpTmplt中的内容复制到一个bash脚本内执行。你需要让脚本产生命名格式为astep[nstep]_state1.${OUT_SUFFIX}这样的输出文件，这样程序就可以找到对应文件来提取能量与梯度。

总结起来，本程序搜索MECP时，需要你在当前目录下有这些文件：
1) 主程序banemecp_static和Frotran优化器baneMECP.f;
2) 几何结构文件init.xyz
3) 配置文件xxx.conf
4) 输入文件xxx.inp
5) (external模式可能需要)你自己的梯度计算脚本
运行方法：./banemecp_static xxx.inp > xxx.out

本程序已经复现过sobMECP的几个简单例子，也成功在GAMESS的MRSF-TDDFT级别算了几个MECP和CI，确认我的修改没把算法部分搞坏。
不过，搜索MECI的时候仍然遇到了一些问题，比如在MRSF-TDDFT级别搜乙烯S0/S1交叉跑很久不收敛。咨询了一下另一个MECP程序的开发者，得知这个可能是Harvey算法处理圆锥交叉本身的缺陷。所以本程序目前看来虽然原则上已经支持MECI的搜索，但实际上仍然可能遇到各种问题。如果将来能补上BPUPD或者惩罚函数算法，才应该能比较舒服地算CI。
9.23 更新后已支持bpupd算法搜索MECI，上述缺点或可忽视了。BPUPD算法在%Control字段内通过algorithm=bpupd指定，该flag的默认值为harvey

作者
Author: wal 时间: 2025-9-16 16:09
可以使用此贴中的程序将opt.trj转换成带有收敛情况信息的Gaussian文件

作者
Author: wal 时间: 2025-9-23 23:57
本帖最后由 wal 于 2025-9-24 15:58 编辑

更新bpupd算法测试版，此部分参考了xmecp的实现，感谢kalinite老师

跑了一个Harvey算法滑铁卢的乙烯跑出来了，不过笔者刚刚接触Fortran纯生瓜蛋子，实现是否有问题还不得而知，有时候AI给我偷偷加料我确实是看不出来的

如果有大佬有兴趣看代码发现了一些愚蠢的逻辑请务必指出
(, 下载次数 Times of downloads: 3)

欢迎光临计算化学公社 (http://bbs.keinsci.com/)