[发明专利]一种基于约束正则模式的复杂函数极小值搜索方法

说明书

技术领域

本发明属于计算化学和物理技术领域，具体涉及一种快速的复杂函数极小值搜索方法，可用于遍历复杂分子、周期性晶体体系的势能面。

背景技术

结构预测和反应路径搜索作为当代化学和物理计算模拟研究的核心任务，对于理解和预测材料的热力学性质和动力学性质具有不可替代的作用。虽然分子动力学模拟作为一种搜索势能面和模拟化学反应过程的常规工具已经得到广泛的应用，但是这一方法的预测能力在处理高维度复杂势能面体系、或具有高反应活化能的化学过程体系时会大幅下降。比如分子动力学方法，不能很好处理蛋白质折叠或碳纳米管的生长过程。

对于这些具有较高反应活化能的化学过程，一种比较常用的方法是先搜索化学过程所对应的过渡态，再找到经过这一过渡态的反应通道。这一方法最明显的缺点就是它高度依赖于科研人员通过化学直觉所预先猜测的可能的反应通道，而这种预先猜测在高维体系中是几乎不可能做出的。另一类常用的方法是通过在预先定义的反应坐标上增加限制条件（比如添加偏置势函数）来增加高活化能过程的采样概率，比如metadynamics方法（Phys. Rev. Lett. 2003，90, 238302）和伞状抽样方法（umbrella sampling）（J. Chem. Phys. 1998，109, 7737）。由于需要根据体系预先定义限制条件，这些方法通常只能用于搜索单步化学反应，不能用于遍历势能面上可能的能量极小点。

与前面方法相反，全局最优化方法，比如BH方法（Basin-Hopping 方法，J. Phys. Chem. A 1997，101, 5111）和遗传算法（Genetic Algorithms in Search, Optimization and Machine Learning， Addison-Wesley Reading, MA, 1989）是以得到函数的最小值为最终目标，也就是势能面上的最稳定结构。这些方法普遍采用了激进的结构形变，来达到跳出势能面上能量极小值区域的目的，所以这些方法完全忽略了势能面上不同极小值间互变的反应通道信息，而且对于复杂势能面的处理效率也大大降低。

发明内容

本发明的目的是为提供一种简便普适的势能面函数搜索方法，克服以往方法中对于高维度势能面的复杂体系预测能力大大降低的缺点。

本发明的另一个目的是提供一种方法，可以同时用于能量极小值的结构预测和反应通道的预测。

为实现上述目的，本发明基于偏置势函数驱动分子动力学思想，通过构造随机的运动模式，并沿着优化得到的约束正则模式，使构型从势能面上一个能量极小值向另一个能量极小值转化，最后利用Metropolis 蒙特卡洛方法确定新找到的能量极小值构型是否被接受。其中描述势能面的能量函数的求解方法包括经验力场方法，量子力学从头算方法，第一性原理密度泛函方法。

本发明方法示意图如图1所示，具体步骤为：

利用输入的原子坐标q_i, i=1,2,…n，已知的势能面能量函数E=f(q_i), i=1,2,…n，能量函数对坐标的一阶导数F=g(q_i), i=1,2,…n，求解能量所有极小值对应的原子坐标，以及连接这些极小值的反应通道，具体步骤如下：

  步骤一、输入一个能量极小值坐标A，并随机产生一个模式N0；

 步骤二、在随机产生的模式N0基础上，通过约束双子方法优化得到一个约束的正则模式N；

步骤三、在原有的势函数基础上，新添加该N方向的一个偏置高斯势函数，则总势能面V_tot可以表示为真实势能面V_real和偏置势函数V_G的和，即                                                ，沿着该N方向外推当前坐标，到由偏置势函数和原有势函数叠加产生的能量极小值附近；

步骤四、能量极小值优化当前坐标；

步骤五、判定是否达到最多的偏置高斯势函数个数（NG）的上限，或达到预设的能量阈值E_r （满足E_max< E_r）；达到则转到步骤六，否则重复步骤二、步骤三、步骤四；

步骤六、去掉所有添加的偏置高斯势函数，能量极小值优化得到新的坐标体系B；