XMECP：多尺度圆锥交叉搜索艺术
XMECP: Reaching State-of-Art in Multi-Scale MECP/MECI Optimization
Jiawei Xu
Released: 2023-02-19 / Updated: 2023-02-19
本文禁止转载 / Reprint of This Article Prohibited
在笔者近几年的工作中时常遇到圆锥交叉。对于分子体系的圆锥交叉优化问题，现有的程序勉强能够满足需求，但对于蛋白质或材料体系等缺少可用的程序。去年年底至今，笔者针对这一问题开发了XMECP程序，以满足对各种尺度圆锥交叉计算的需求，经过大量测试，我们认为程序已经能够满足绝大多数研究者的日常需求。日前相关文章已经发表于……，文中包含了大量程序应用范例，即使是研究不涉及圆锥交叉的工作者，看了该文的实例后也一定能有所收获。对于有这方面需求的工作者，欢迎下载XMECP使用。如果在研究工作中使用了XMECP优化交叉点或极小点、过渡态等，请在文章正文中恰当引用XMECP程序，如：Jiawei Xu, XMECP, https://www.github.com/kalinitejxu/xmecp。
本文将拣选一些有意义的实例来讲解XMECP程序的使用方法和技巧。特别要强调的是，笔者对于很多程序或方法采用一些非常简单的例子说明可行性的行为感到非常鄙视，原因是非常简单的例子完全无法体现方法的普适性，一般也没什么实际意义。很多程序或方法仅用一些只有几个原子的小分子做测试，实际一使用发现对于稍大些的体系完全没法用，要么耗时良久，要么完全跑不了，表现地也很差劲，这是非常不务实的研究风格，脱离了广大群众对程序的实际需求。因此XMECP所展现的例子覆盖了各种常见的体系，笔者有信心认为，在使用恰当的前提下，若有XMECP仍无法处理的体系，目前一定也没有公开的程序可以处理了。而且XMECP力求灵活和方便，完全不需要任何附加门槛，掌握基本的计算化学知识的人就可以非常轻松地下载即用。对于高级玩家，还可以很方便地新增对各种程序的接口。
1. 程序运行与输入文件
XMECP程序完全由Python写成，因此机子上必须有Python3及NumPy库。除此以外，对于需要使用内坐标相关功能的用户，尤其是涉及需要求Wilson’s B matrix的步骤，还需要额外安装Tensorflow库。安装方法参见相关库的官方网站。XMECP程序通过以下命令运行：
python -u xmecp.py mol.inp mol.xyz > mol.out其中mol.inp包含了计算任务的有关设定和控制参数，mol.xyz包含分子坐标，为标准的.xyz格式。选项-u则会将程序的输出实时写入输出文件mol.out，否则仅在任务结束后一次性写入输出文件，不便于监控任务情况。对于QM/MM计算，还需要根据实际情况额外提供相关参数文件，目前QM/MM计算仅支持Tcl-ChemShell接口使用AMBER参数(.prmtop)的情况，其他的情况可能需要修改源代码或新增接口，也可能不需要。
1.1 关键词输入
XMECP的关键词部分采用类似ORCA输入文件的写法，对关键词进行分块，以“#section”和“#end”作为块起始和块终止的标志。一个简单的MECP优化任务的输入文件如下：
#JobType MECP_Harvey QM
#Section Control
JobName Mol-MECP
SysName Mol
NProc 32
#End
#Section QM
QM_Charge 0
QM_Mult 1
QM_Software ORCA
QM_Path /home/jwxu/programs/orca_5_0_3/orca
#End
#Section MECP
State1 1 0 (default)
State2 3 1 (default)
#End
#Section Opt
MaxStep 0.10 (default)
MaxCycle 128 (default)
Method GDIIS (default)
#End
#Section Template
QM_Template qm.tpl (default)
#End所有关键词均不区分大小写，除#JobType必须在第一行外，各块的先后顺序也没有要求，也可以省略不需要的块。事实上所有关键词都有其默认值，例如以上注明了default的关键词一般都可以省略，具体可以在inp_proc.py中查询，也可以根据实际情况修改。例如QM_Path仅在使用ORCA作为QM软件时需要，因此可以在inp_proc.py中修改为机子上实际ORCA可执行文件的所在路径，以后便无需再写进输入文件中。下面我们简要介绍以上关键词的含义，更多选项将在后续的实例中介绍。
1.1.1 #JobType
接受两个关键词，分别指定计算类型和尺度。计算类型包括：Opt，优化极小点结构；TS：优化过渡态结构；MECP_Harvey，使用Harvey的算法优化MECP；MECP_PF，使用penalty function方法优化MECP；MECP_BPUPD，使用branching plane updating方法优化MECP。此外源代码中还有封印的MECP_Harvey_NoEd，表示使用Harvey的方法，但并不将两态能量差作为优化收敛的标准之一，高级玩家可能知道其作用，新手则不建议尝试，以防得到完全没有意义的结果。计算尺度包括QM和QMMM，对于使用Gaussian、ORCA或Q-Chem等程序做ONIOM方法的情况，此处也应选用QM。
另外有两个待发布的功能：MD，动力学计算；Freq：频率计算。
1.1.2 #Section Control
JobName和SysName分别指定任务名称和体系名称，用于中间文件和临时文件的命名。在使用Tcl-ChemShell计算时，会自动将模板文件中的${sysname}字段替换为SysName的值。NProc指定任务可使用的最大核数。
1.1.3 #Section QM
QM_Charge和QM_Mult分别指定体系的电荷和自旋多重度，对于QM/MM计算则是指QM区域的电荷和自旋多重度。QM_Software指定使用的QM软件，注意在使用ORCA时需将QM_Path指定为实际ORCA可执行文件所在的路径。
1.1.4 #Section MECP
State1和State2分别指定圆锥交叉两态的自旋多重度和态序号，其中基态的态序号为0。对于spin-flip TD-DFT计算，则是参考态的态序号为0。
1.1.5 #Section Opt
MaxStep指定优化的步长上限，单位为a.u.。MaxCycle指定优化的步数上限。Method指定优化方法，默认为GDIIS，也可使用BFGS或GEDIIS。
1.1.6 #Section Template
QM_Template指定相应量化软件的输入文件模板。以ORCA为例，以下是一个qm.tpl的基本内容，可根据实际需求修改：
! B3LYP/G D3BJ def2-SVP def2/J RIJCOSX
%basis
NewGTO Fe
S 3
1 20.5130280 0.2347140
......
G 1
1 2.2690000 1.0
end
NewECP Fe
N_core 10
lmax f
s 2
1 20.930000000 253.749588000 2
2 9.445000000 37.922845000 2
......
f 1
1 1.000000000 0.000000000 2
end
NewAuxJGTO Fe "def2/J" end
end
%scf
MaxIter 300
end
请注意：XMECP会自动根据实际情况增减关键词，因此有些关键词是不需要写到模板文件中的，例如%pal nprocs … end、uks、%tddft … end、MOread等，都不建议写入模板，以防关键词冲突或任务无意义。
2. 使用范例：分子体系
3. 使用范例：生物体系
4. 使用范例：材料体系
5. 扩展功能
Acknowledgement
* The author thanks Dr. Jian Hao, PI Minyi Zhang and PI Chunsen Li from Fujian Institute of Research on the Structure of Matter for their help and advice in this work.
* For tests vs. Q-Chem: the author thanks Prof. Haiyan Wei from Nanjing Normal University for providing software. Test jobs were completed by Qiyan Cao during her graduation thesis.
* The author also thanks Dr. Jingxiang Zou from Nanjing University for his help with MOKIT package and multi-reference calculation.