一种利用遗传算法优化反应过渡态的方法与流程

[0001]
本发明涉及一种利用遗传算法优化反应过渡态的方法，属于计算化学领域。


背景技术：

[0002]
在分子构象转变或者化学反应中往往都会存在过渡态，过渡态结构对应着势能面上反应路径中能量最高点，通过最小能量路径连接反应物和产物。确定过渡态有助于了解反应机理，计算能垒推算反应速率。
[0003]
但在化学反应过程中，过渡态不稳定，存在时间极短，无法通过实验的方法获得，计算化学方法是目前寻找反应过渡态最好的办法，但目前在量子化学计算上仍有一些困难。目前寻找化学反应过渡态主要有3种方法：ts、qst2以及qst3，目前ts法是常用方法。在高斯中选用opt＝ts方法来寻找过渡态，需要提供一个初猜结构，初猜结构对优化结果有很大影响，往往会存在如下几个问题：1、初猜结构错误；2、初猜结构不是本反应的过渡态；3、初猜结构的能量不是本反应路径中过渡态能量最小值。初猜结构的不准确容易导致寻找过渡态的失败，且初猜不正确往往需要重新改变初猜结构，耗时长且不一定能够获得最优过渡态。
[0004]
遗传算法是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。即“适者生存，不适者淘汰”。主要过程：从代表问题可能潜在的解集的一个种群(population)开始，种群又是由个体(individual)组成的，每个个体都是经由基因(gene)编码的，所以个体实际上是染色体chromosome)带有特征的实体。染色体是物质遗传的载体，多种基因的集合才决定了外在的表现。因此，一开始需要将所需的表现型映射到内容控制的基因型即编码工作。由于仿基因编码比较繁琐，所以进行简化采用二进制编码，初代种群产生之后，按照适者生存和优胜劣汰的原理，逐代(generation)演化产生出越来越好的近似解，在每一代，根据问题域中个体的适应度(fitness)大小选择(selection)个体，并借助于自然遗传学的遗传算子(genetic operators)进行组合交叉(crossover)和变异(mutation)，产生出代表新的解集的种群。这个过程将导致种群像自然进化一样的后生代种群比前代更加适应于环境，末代种群中的最优个体经过解码(decoding)，可以作为问题近似最优解。因此，寻找一种更简便且准确度更高的计算化学反应过渡态的方法对于化学反应过渡态的确定以及对后续化学工业的指导具有重要意义。


技术实现要素：

[0005]
针对现有的化学计算反应过渡态的过程中存在的初猜结构不易确定、随机性较强，且寻找反应过渡态的过程耗时长，计算结果往往错误或者不是理想结果，以及计算周期长的问题，本发明提供了一种利用遗传算法优化反应过渡态的方法。本发明方法通过输入反应物的分子式自动转化为结构式再利用遗传算法寻找出反应物可能的过渡态，计算出对应的反应路径，最终选出所需的过渡态及反应路径。本发明无需人为去调分子结构，计算过
程快且结果更加准确。
[0006]
具体的，本发明的技术方案为，一种利用遗传算法优化反应过渡态的方法，所述方法包括：
[0007]
(1)输入反应物的分子式，通过rdkit将分子式转化为结构式并输出对应的.mol文件；
[0008]
(2)利用multiwfn将步骤(1)得到的.mol文件转化为.gjf文件；
[0009]
(3)将.gjf文件中的目标原子的坐标设为变量，同时，设定坐标变量的范围，并随机生成多组目标原子坐标的初值；
[0010]
(4)输入目标原子坐标变量的初值，利用基组算法进行计算，计算得到.out文件，进行判断：当结果非正常结束，则删除结果并代入其他初值进行迭代计算；当结果为正常结束，但是频率计算结果中无虚频或者有超过一个虚频，则删除结果并代入其他初值进行迭代计算；当结果为正常结束，且频率计算结果中有且只有一个虚频，则将.out文件通过文件的读写方式写入并计算该.out文件中对应的过渡态所对应的反应路径，之后代入其他初值进行迭代计算；
[0011]
(5)初值迭代完成后，整理步骤(4)计算得到的符合条件的初值，并利用遗传算法自动变异和交叉组合得到新的初值，重复步骤(4)和(5)，直至将符合条件的初值计算完毕；
[0012]
(6)迭代完成后，整理所有计算得到的过渡态及其对应的反应路径，根据反应产物来判断和确定符合条件的过渡态。
[0013]
在本发明的一种实施方式中，步骤(1)中，所述输入反应物的分子式应根据rdkit的规则进行输入，其中，所述rdkit为一款开源化学信息学软件。
[0014]
在本发明的一种实施方式中，所述multiwfn是一款量子化学波函数分析程序。
[0015]
在本发明的一种实施方式中，步骤(2)中，需要根据实际情况将选用的基组以及注意事项(例如溶剂化效应)写入.gjf文件，与传统计算方法的输入一致。
[0016]
在本发明的一种实施方式中，步骤(4)中基组算法为：只要能在高斯软件上计算的算法在此步均能写入文件并调用高斯软件；包括半经验法、从头算法、密度泛函方法、mpn方法中的任一种。
[0017]
在本发明的一种实施方式中，所述坐标变量的范围优选为其余三个原子坐标范围的最大值或更大的范围，即其余三个原子的坐标的x、y、z的最大值和最小值分别作为坐标变量的x、y、z坐标的最大值和最小值，即获得目标原子的坐标范围，或大于前述坐标范围的坐标范围。
[0018]
在本发明的一种实施方式中，所述坐标变量的初值是在坐标变量的范围内随机生成的目标原子的坐标值。
[0019]
在本发明的一种实施方式中，步骤(3)中所述的多组初值可以根据实际情况选定，例如100组、200组、300组等。
[0020]
在本发明的一种实施方式中，步骤(5)中所述遗传算法自动变异、交叉组合得到新的初值，即通过设定交叉率和变异率的参数值，然后利用遗传算法对符合条件的初值进行变异交叉组合，得到新的初值。
[0021]
在本发明的一种实施方式中，所述交叉率为0.6～1，优选0.8。
[0022]
在本发明的一种实施方式中，所述变异率小于0.1，优选0.003。
[0023]
在本发明的一种实施方式中，步骤(4)中，当结果为正常结束，且有且只有一个虚频时，将此时得到的.out文件中的能量值与之前得到的所有的满足条件的.out文件中的能量值进行比对，当其与之前的能量值差值在
±
0.00002hatree范围之内时(其中，1hartree＝2625.5kj mol-1
)，则将此次获得的能量值进行删除，且无需计算此时对应的反应路径；删除的原因是由于二者能量值差别过小，对应的反应路径一致。
[0024]
在本发明的一种实施方式中，所述反应路径的计算方式为常规的反应路径的计算方式。
[0025]
在本发明的一种实施方式中，迭代过程中的迭代次数可以根据实际需要进行自行设定。
[0026]
本发明还提供了上述利用遗传算法优化反应过渡态的方法在计算化学中的应用。
[0027]
有益效果：
[0028]
(1)本发明根据反应路径直接输入反应物分子式，无需调整预猜结构，直接调用高斯软件后根据输出文件直接判断结果是否满足过渡态的结果的要求(有且只有一个虚频)，并根据得出的满足要求的输出结果直接调用高斯软件计算出对应的反应路径，最后从所有计算结果中找到所需的反应路径对应的过渡态即可。此过程无需人工时刻关注计算结果，大大节省寻找周期，并且可以对不同路径的过渡态做比对分析。
[0029]
(2)本发明方法能够避免人为输入时候的不准确性，以及在优化过程中微调结构存在误差导致寻找优化过渡态过程中花费大量的时间甚至无果；另外，本发明方法找出的过渡态相比较人为调结构找到的过渡态能量值可能会更低，结果更加符合真实情况。
附图说明
[0030]
图1为实施例2计算过程中得到满足条件的结果，其中，(a)为迭代第1～118次的结果；(b)为迭代第119～237次的结果；(c)为迭代第238～304次的结果；(d)为迭代第305～337次的结果；(e)为迭代第338～394次的结果；(f)为迭代第395～410次的结果；(g)为迭代第411～518次的结果；(h)为迭代第519～540次的结果。
[0031]
图2为实施例2最终确认的hcooh与h2o2反应的反应路径。
[0032]
图3为实施例3最终确认的hcooh与h2o2反应的反应路径。
具体实施方式
[0033]
下面结合实施例对本发明作进一步的描述，但本发明的实施方式不限于此。
[0034]
实施例1
[0035]
一种利用遗传算法优化反应过渡态的方法，所述方法包括：
[0036]
(1)输入反应物的分子式，通过rdkit将分子式转化为结构式并输出对应的.mol文件；
[0037]
(2)利用multiwfn将步骤(1)得到的.mol文件转化为.gjf文件；
[0038]
(3)将.gjf文件中的目标原子的坐标设为变量，同时，设定坐标变量的范围，并随机生成多组目标原子坐标的初值；
[0039]
(4)输入目标原子坐标变量的初值，利用基组算法进行计算，计算得到.out文件，进行判断：当结果非正常结束，则删除结果并代入其他初值进行迭代计算；当结果为正常结
束，但是频率计算结果中无虚频或者有超过一个虚频，则删除结果并代入其他初值进行迭代计算；当结果为正常结束，且频率计算结果中有且只有一个虚频，则将.out文件通过文件的读写方式写入并计算该.out文件中对应的过渡态所对应的反应路径，之后代入其他初值进行迭代计算；
[0040]
(5)初值迭代完成后，整理步骤(4)计算得到的符合条件的初值，并利用遗传算法自动变异和交叉组合得到新的初值，重复步骤(4)和(5)，直至将符合条件的初值计算完毕；
[0041]
(6)迭代完成后，整理所有计算得到的过渡态及其对应的反应路径，根据反应产物来判断和确定符合条件的过渡态。
[0042]
进一步的，步骤(2)中，需要根据实际情况将选用的基组以及注意事项(例如溶剂化效应)写入.gjf文件，与传统计算方法的输入一致。
[0043]
进一步的，步骤(4)中所述基组算法包括半经验法、从头算法、密度泛函方法、mpn方法中的任一种。
[0044]
进一步的，所述坐标变量的范围优选为其余三个原子坐标范围的最大值或更大的范围，即其余三个原子的坐标的x、y、z的最大值和最小值分别作为坐标变量的x、y、z坐标的最大值和最小值，即获得目标原子的坐标范围，或大于前述坐标范围的坐标范围。
[0045]
进一步的，步骤(5)中所述遗传算法自动变异、交叉组合得到新的初值，即通过设定交叉率和变异率的参数值，然后利用遗传算法对符合条件的初值进行变异交叉组合，得到新的初值；所述交叉率为0.6～1，优选0.8；所述变异率小于0.1，优选0.003。
[0046]
进一步的，步骤(4)中，当结果为正常结束，且有且只有一个虚频时，将此时得到的.out文件中的能量值与之前得到的所有的满足条件的.out文件中的能量值进行比对，当其与之前的能量值差值在
±
0.00002hatree范围之内时(其中，1hartree＝2625.5kj mol-1
)，则将此次获得的能量值进行删除，且无需计算此时对应的反应路径。具体的，删除的原因是由于二者能量值差别过小，对应的反应路径一致。
[0047]
实施例2
[0048]
下面以hcooh与h2o2在常温下真空中反应为例，利用本发明所述的遗传算法优化二者反应过程中的过渡态。
[0049]
(1)输入反应物的分子式o[ch+](＝o)[oh+]o，rdkit将分子式转化为结构式并输出对应的.mol文件；
[0050]
(2)利用multiwfn将步骤(1)得到的.mol文件转化为.gjf文件，并将.gjf文件中各原子的坐标转换为笛卡尔坐标，并选用计算方法及基组为b3lyp/6-31g(d)，写入.gjf文件；
[0051]
(3)将.gjf文件中的目标原子的坐标设为变量，同时，根据其他三个原子的坐标c[0.46740000,0.12730000,0.36070000]、o[0.08460000,-0.31430000,1.47030000]、o[-0.58080000,0.21880000,-0.51860000]设定目标原子h坐标变量的范围为x[-0.6,0.6],y[-0.4,0.3],z[-0.6,1.5]，并随机生成100组初值；
[0052]
(4)输入坐标变量的初值并利用基组算法进行计算，计算得到.out文件，进行判断；当结果非正常结束，则删除结果并代入其他初值进行迭代计算；当结果为正常结束，但是频率计算结果中无虚频或者有超过一个虚频，则删除结果并代入其他初值进行迭代计算；当结果为正常结束，且频率计算结果中有且只有一个虚频，则将.out文件通过文件的读写方式写入并计算该.out文件中对应的过渡态所对应的反应路径(irc反应路径)，之后代
入其他初值进行迭代计算；
[0053]
(5)步骤(4)的初值迭代完成后，整理步骤(4)计算得到的符合条件的初值，并利用遗传算法自动变异和交叉组合得到新的初值，交叉率为0.8，各个位点的变异率为0.003，重复步骤(4)和(5)。当结果为正常结束，有且只有一个虚频时，将此时得到的.out文件中的能量值与之前得到的所有的满足条件的.out文件中的能量值进行比对，当其与之前的能量值差值在
±
0.00002hatree范围之内时(其中，1hartree＝2625.5kj mol-1
)，则将此次获得的能量值进行删除，且无需计算此时对应的反应路径。
[0054]
(6)迭代完成后，整理所有计算得到的过渡态及其对应的反应路径，根据反应产物来判断和确定符合条件的过渡态。
[0055]
最终一共计算得到了符合条件的初值31个，结果见图1(a-h)所示，能量值见表1。
[0056]
表1实施例2计算得到符合条件的初值及其对应的过渡态的能量值
[0057]
[0058][0059]
由此可见，此处共进行了540次迭代，筛选出了31个可能解，共耗时3天(此次运算在普通电脑的虚拟机上运行，若到集群计算机上会极大缩短计算时间)。结合后续过渡态的反应路径的筛选，可知第237次迭代获得的结果是最终优化后的过渡态，对应的初值坐标是x:0.33630679、y：-0.12623329、z：-0.27125001，吉布斯自由能是-341.200042hartree，对应的反应路径见图2，即双氧水上的一个氢与甲酸中羰基上的氧相连，双氧水中失去氢的氧与甲酸上的碳相连。
[0060]
该对应的最优化的过渡态的信息如下表2。
[0061]
表2实施例2中计算得到的最优过渡态所对应的四个能量值
[0062]
能量能量值(hartree)电子能与零点能之和-314.170799电子能与热能之和-341.164919电子焓和热焓之和-341.163975吉布斯自由能-341.200042
[0063]
对比例1
[0064]
当选用传统的计算方法来计算hcooh与h2o2反应的过渡态时，方法如下：
[0065]
先在gaussian软件中输入hcooh与h2o2的结构式，并调整反应物的相对位置(特别是反应位点的相对位置)后保存为.gjf文件，然后调用gaussian软件进行计算。查看输出结果会出现如下情况：1、未能正常收敛；2、收敛但不是一个虚频，不是所需的过渡态；3、收敛且只有一个虚频。如果未收敛或者不止一个虚频，再手动调整反应基团的键长、键角及二面角然后再调用高斯进行计算，如果收敛且只有一个虚频可以进行irc反应路径计算，再进行确定找出的过渡态是否为连接反应物与产物的过渡态。
[0066]
经过了多次调整结构，耗时两周，最终获得了过渡态结构，对应过渡态的反应路径见图3，该过渡态对应的吉布斯自由能为-341.199807hartree，结果见表3。
[0067]
表3实施例3中计算得到的最优过渡态所对应的四个能量值
[0068][0069][0070]
从图2和图3可以看出两者的反应路径一样，均是双氧水上的一个氢被甲酸中羰基
上的氧吸引，双氧水上失去氢离子的氧与羰基上的碳相连。但通过遗传算法自动找到的此路径对应的吉布斯自由能最低值为-341.200042hartree，而多次手动调整算出的吉布斯自由能最低值为-341.199807hartree，本发明方法计算得到的过渡态的能量值较对比例1更低，更加符合实际状况，且本发明方法无需手动调控分子结构，耗时短，优势明显。
[0071]
虽然本发明已以较佳实施例公开如上，但其并非用以限定本发明，任何熟悉此技术的人，在不脱离本发明的精神和范围内，都可做各种的改动与修饰，因此本发明的保护范围应该以权利要求书所界定的为准。