适用于更多程序的sobMECP修改版(9.23更新MECI支持)

wal

6月左右的时候尝试过把sobMECP接上MOKIT，但因为sobMECP迭代部分比较零落，没搞明白。最近遇到了需要MRSF-TDDFT搜索MECP的情况，但GAMESS不原生支持，于是再次把目光瞄向了sobMECP
GitHub：banemecp
主程序： baneMECP.f (35.79 KB, 下载次数 Times of downloads: 7)

2025-9-24 17:58 上传 Uploaded

点击下载
Click to download banemecp_static (1.25 MB, 下载次数 Times of downloads: 5)

2025-9-24 17:58 上传 Uploaded

点击下载
Click to download
配置文件与示例输入： banemecp.zip (9.4 KB, 下载次数 Times of downloads: 7)

2025-9-14 14:43 上传 Uploaded

点击下载
Click to download
一些算例 sobmecp原例.zip (846.98 KB, 下载次数 Times of downloads: 9)

2025-9-14 15:16 上传 Uploaded

点击下载
Click to download MRSF-TDDFT.zip (3.34 MB, 下载次数 Times of downloads: 12)

2025-9-14 15:18 上传 Uploaded

点击下载
Click to download

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


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

1. [main]
   
2. env='source ~/.bashrc
   
3. module load g16'
   
4. # As a suffix, forms ${ACTUAL_INP}
   
5. INPUT_SUFFIX=gjf
   
6. # As a suffix, forms ${ACTUAL_OUT}
   
7. OUTPUT_SUFFIX=log
   
8. RUN_CMD='echo ">>>>> Start g16 calculation"
   
9. g16 ${ACTUAL_INP} > ${ACTUAL_OUT}
   
10. echo "g16 done. <<<<<<"'

复制代码

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

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

复制代码

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



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

1. %Prog g16
   
2. %Geom init.xyz
   
3. %Control
   
4. maxcyc=100
   
5. tmpdir=.
   
6. keeptmp=true
   
7. restart=false
   
8. end

复制代码

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

1. %InpTmplt1
   
2. %chk=state1.chk
   
3. # B3LYP/def2svp force [guess(read)]
   
4. 
   
5. state1
   
6. 
   
7. 0 1
   
8. [geom]
   
9. 
   
10. 
    
11. end
    
12. 
    
13. %InpTmplt2
    
14. %chk=state2.chk
    
15. # B3LYP/def2svp force [guess(read)]
    
16. 
    
17. state2
    
18. 
    
19. 0 3
    
20. [geom]
    
21. 
    
22. 
    
23. end

复制代码

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

1. %GrpTmplt
   
2. state1=scf
   
3. state2=td
   
4. end
   

复制代码

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

1. %InpTmplt1
   
2. echo ">>>>> Start gamess calculation1"
   
3. cp ../gamess.sh .
   
4. bash gamess.sh [xyzfile] --iroot 1 --basename astep[nstep]_state1 >> mrsfout.tmp 2>&1
   
5. end
   
6. 
   
7. %InpTmplt2
   
8. echo ">>>>> Start gamess calculation2"
   
9. cp ../gamess.sh .
   
10. bash gamess.sh [xyzfile] --iroot 2 --basename astep[nstep]_state2 >> mrsfout.tmp 2>&1
    
11. cp opt.trj ..
    
12. 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

wal

可以使用此贴中的程序将opt.trj转换成带有收敛情况信息的Gaussian文件

wal

更新bpupd算法测试版，此部分参考了xmecp的实现，感谢kalinite老师
跑了一个Harvey算法滑铁卢的乙烯跑出来了，不过笔者刚刚接触Fortran纯生瓜蛋子，实现是否有问题还不得而知，有时候AI给我偷偷加料我确实是看不出来的如果有大佬有兴趣看代码发现了一些愚蠢的逻辑请务必指出
ethene.zip (3.7 MB, 下载次数 Times of downloads: 3)

2025-9-24 15:58 上传 Uploaded

点击下载
Click to download