最小能量路径搜索软件MEPsearcher

chenxin199261

简介
MEPsearcher基于string计算方法，用于搜索二维全局离散势能面上的能量最小路径Mminimun Energy Path(MEP),通常MEP可以完整揭示反应机理。过渡态、活化能等信息可以从MEP中获取。能量最低路径的搜索可以利用于蛋白质折叠、化学反应、相变等过程的机理研究。
下载地址:  https://github.com/chenxin199261/MEPSearcher
具体信息: http://chen-xin.tk/blog/index.php/archives/44/

string方法简介
在势能面 V(x) 上，假设路径 γ 连接极小值点a、b，γ 满足如下关系。

(∇V)⊥(γ)=0

寻找MEP，首先定义N个珠子(beads、images),{φi,i=0,1,2,…,n}


,分别对N个beads进行优化。

φ˙i=−∇V(φi)

优化过程同时将beads均匀分布。采用“Parametrization by equal arc length”方法
原理很简单，简单一句话：

将一些珠子放在势能面山坡上.松开珠子,珠子都会进入能量最低点.
如果们将珠子用线穿起来,那么珠子稳定状态下就是在能量最低路径上.

技术细节
MEPsearcher采用python开发，核心数值计算部分调用ALGlib数学库。
当前版本只支持最速下降优化方法(DS),未来版本考虑加入共轭梯度，准牛顿等方法。
插值算法采用三次样条插值(cubic spline)。ALGlib提供很多插值方法，作者并未尝试，希望大家对着手册换换其它的插值算法。

使用方法
1. 准备文件你只需要准备两个输入文件，势能面文件PES.data，初猜生成文件guess.data。

• 势能面文件写成三列X Y Energy即可。可参考test文件夹中PES.data格式
• 准备初始猜测需要对势能面有一定的直观感觉，给出起始点、终止点以及初猜珠子(beads)的个数。
  

2.参数设定需要进入MEPsearcher.py中对参数进行设定，只有三个参数需要测定，分别是最大迭代次数，收敛限以及优化步长。

• 最大迭代次数，Max_iter默认1000步，可以根据具体需要调大或者调小。
• 收敛限，用与上次优化结果的距离差加和表示。Tol默认1e-5,beads过多，可以适当调小、调大这个值。
• 优化步长值，h，是最速下降优化方法一个必需参数。这个值的设定麻烦些。过小导致优化缓慢，过大导致直接优化到势能面错误的区域。需要使用者反复调试。希望在下一个版本中，能启用其它不需要步长的优化方法。
  

3.运行赋予MEPsearch.py执行权限

1. chmod +x MEPsearcher.py

复制代码

直接运行

1. ./MEPsearcher.py

复制代码

计算过程每隔50步迭代，会打印在屏幕上。
运行结束后，产生out.data文件。
这个文件存储了MEP的最终结果。分三列表示，分别为:
X Y Energy

测试结果

smutao

高维的版本有没有 ？
比如要找电子密度分布的Maximum Electron Density Path我猜Sob的Multiwfn里应该有这个功能

sobereva

smutao 发表于 2016-2-19 23:40
高维的版本有没有 ？
比如要找电子密度分布的Maximum Electron Density Path我猜Sob的Multiwfn里应该有这 ...

那个是Multiwfn的产生键径。
具体来说是先用牛顿法根据初猜搜索到鞍点，然后由此跟踪密度下降最快路径。这类似于高斯等程序产生IRC的算法。LZ的则类似于第一性原理程序支持较多的NEB方法，已知极小点，一下子就得到IRC，原理上可以用于任意维。

chenxin199261

sobereva 发表于 2016-2-20 01:31
那个是Multiwfn的产生键径。
具体来说是先用牛顿法根据初猜搜索到鞍点，然后由此跟踪密度下降最快路径。 ...

我的不用找到极值点也可以得到那个路径，比我我初猜的那样，给出一条线。这个是string方法的一个好处。
NEB好像必须锁定两个极小值点。

chenxin199261

smutao 发表于 2016-2-19 23:40
高维的版本有没有 ？
比如要找电子密度分布的Maximum Electron Density Path我猜Sob的Multiwfn里应该有这 ...

高维的暂时没有，很容易改的。可以用一个二维数组存X,Y,Z，而不是我那样分别对X，Y操作。
我没做高维的，因为高维的路径无法可视化，意义不大。