Please wait a minute...
物理化学学报  2018, Vol. 34 Issue (10): 1179-1188    DOI: 10.3866/PKU.WHXB201803161
所属专题: 材料科学的分子模拟
论文     
关于副本交换分子动力学模拟复杂化学反应的研究
辛亮1,3,孙淮1,2,3,*()
1 上海交通大学化学化工学院,上海 200240
2 上海交通大学材料基因组联合研究中心,上海 200240
3 吉林大学无机合成与制备化学国家重点实验室,长春 130012
On the Simulation of Complex Reactions Using Replica Exchange Molecular Dynamics (REMD)
Liang XIN1,3,Huai SUN1,2,3,*()
1 School of Chemistry and Chemical Engineering, Shanghai Jiao Tong University, Shanghai 200240, P. R. China
2 Materials Genome Initiative Center, Shanghai Jiao Tong University, Shanghai 200240, P. R. China
3 State Key Laboratory of Inorganic Synthesis and Preparative Chemistry, Jilin University, Changchun 130012, P. R. China
 全文: PDF(2347 KB)   HTML 输出: BibTeX | EndNote (RIS) |
摘要:

本文研究用温度副本交换分子动力学(T-REMD)和哈密顿副本交换分子动力学(H-REMD)方法模拟复杂化学反应的问题。使用具有不同活化能和反应能的简单置换反应模型,我们检验了上述两种方法用来预测反应平衡产物的效率和应用范围。T-REMD方法对具有适度活化能(约< 20 kcal·mol-1)或者反应能量(< 3 kcal·mol-1)的放热反应是有效的。由于在相空间的不完整采样,对于同时具有高活化能和反应能量的反应其模拟效率有严重障碍,并且对于吸热反应问题更为显着。另一方面,H-REMD对一系列具有不同活化能的反应能的模型表现出色,与T-REMD相比,H-REMD可以使用更少的副本获得优异的结果。

关键词: 副本交换分子动力学复杂反应温度哈密顿    
Abstract:

A complex reaction, such as combustion, polymerization, and zeolite synthesis, involves a large number of elementary reactions and chemical species. Given a set of elementary reactions, the apparent reaction rates, population of chemical species, and energy distribution as functions of time can be derived using deterministic or stochastic kinetic models. However, for many complex reactions, the corresponding elementary reactions are unknown. Molecular dynamics (MD) simulation, which is based on forces calculated by using either quantum mechanical methods or pre-parameterized reactive force fields, offers a possibility to probe the reaction mechanism from the first principles. Unfortunately, most reactions take place on timescales far above that of molecular simulation, which is considered to be a well-known rare event problem. The molecules may undergo numerous collisions and follow many pathways to find a favorable route to react. Often, the simulation trajectory can be trapped in a local minimum separated from others by high free-energy barriers; thus, crossing these barriers requires prohibitively long simulation times. Due to this timescale limitation, simulations are often conducted on very small systems or at unrealistically high temperatures, which might hinder their validity. In order to model complex reactions under conditions comparable with those of the experiments, enhanced sampling techniques are required. The replica exchange molecular dynamics (REMD) is one of the most popular enhance sampling techniques. By running multiple replicas of a simulation system using one or several controlling variables and exchanging the replicas according to the Metropolis acceptance rule, the phase space can be explored more efficiently. However, most published work on the REMD method focuses on the conformational changes of biological molecules or simple reactions that can be described by a reaction coordinate. The optimized parameters of such simulations may not be suitable for simulations of complex reactions, in which the energy changes are much more dramatic than those associated with conformational changes and the hundreds elementary reactions through numerous pathways are unknown prior to the simulations. Therefore, it is necessary to investigate how to use the REMD method efficiently for the simulation of complex reactions. In this work, we examined the REMD method using temperature (T-REMD) and Hamiltonian (H-REMD) as the controlling variable respectively. In order to quantitatively validate the simulation results against direct simulations and analytic solutions, we performed the study based on a simple replacement reaction (A + BC = AB + C) with variable energy barrier heights and reaction energies described using the ReaxFF functional forms. The aim was to optimize the simulation parameters including number, sequence, and swap frequency of the replicas. The T-REMD method was found to be efficient for modeling exothermic reactions of modest reaction energy (< 3 kcal∙mol-1) or activation energy (ca. < 20 kcal∙mol-1). The efficiency was severely impaired for reactions with high activation and reaction energies. The analysis of the simulation trajectory revealed that the problem was intrinsic and could not be solved by adjusting the simulation parameters since the phase space sampled using T-REMD was localized in the region favored by high (artificial for speed-up) temperatures, which is different from the region favored by low (experimental) temperatures. This issue was aggravated in the case of endothermic reactions. On the other hand, the H-REMD run on a series of potential surfaces having different activation energies was demonstrated to be remarkably robust. Since the energy barrier only reduces the reaction rates, while the phase space controlled by the reaction energy differences remains unchanged at a fixed temperature, excellent results were obtained with fewer replicas by using H-REMD. It is evident that H-REMD is a more suitable method for the simulation of complex reactions.

Key words: Replica exchange    Molecular dynamics    Complex reaction    Temperature    Hamilton
收稿日期: 2018-02-08 出版日期: 2018-04-13
中图分类号:  O645  
基金资助: 国家自然科学基金(21073119);国家自然科学基金(21173146);国家自然科学基金(21473112);国家重点基础研究发展规划项目(973)(2014CB239702)
通讯作者: 孙淮     E-mail: huaisun@sjtu.edu.cn
服务  
把本文推荐给朋友
加入引用管理器
E-mail Alert
RSS
作者相关文章  
辛亮
孙淮

引用本文:

辛亮,孙淮. 关于副本交换分子动力学模拟复杂化学反应的研究[J]. 物理化学学报, 2018, 34(10): 1179-1188, 10.3866/PKU.WHXB201803161

Liang XIN,Huai SUN. On the Simulation of Complex Reactions Using Replica Exchange Molecular Dynamics (REMD). Acta Phys. -Chim. Sin., 2018, 34(10): 1179-1188, 10.3866/PKU.WHXB201803161.

链接本文:

http://www.whxb.pku.edu.cn/CN/10.3866/PKU.WHXB201803161        http://www.whxb.pku.edu.cn/CN/Y2018/V34/I10/1179

Fig 1  Left) Contour plot of the potential energy surface for model reaction AB + C ⇋ A + BC in terms of bond lengths rAB and rBC, the collision angle A―B―C is fixed at 180° which corresponds to the minimum energy path. (Right) The reactive model for a simple reaction AB + C ⇋ A + BC with different activation energies and reaction energies. Ǻ = 0.1 nm
Fig 2  The canonical heat capacities calculated for models (-3, Ea) with different activation energies in 15, 20, 25, 30, 40 and 50 kcal•mol-1. 1 kcal•mol-1 = 4.1868 kJ•mol-1
Fig 3  Convergence time measured by the number of reactants for a T-REMD simulation at 1000 K using model (-3, 25). The blue line represents instantaneous number, the green line is the 25-ps block averaged data, and the red line represents the equilibrium value. The first time when the block averaged data reaches the equilibrium line is defined as the convergence time. Color online.
EAF/ps-1 tc/ns Nex
1 1.00 1, 000
2 0.50 1, 000
4 0.25 1, 000
8 0.25 2, 000
40 0.25 10, 000
Table 1  Exchange attempt frequency (EAF), convergence time (tc) and number of exchange attempts (Nex) of T-REMD simulations on model (-3, 25).
Activation energy Reaction energy
0 -1.5 -3 -5 -10
MD 15 3.42 2.97 2.46 3.43 2.21
20 2.42 2.83 3.19 2.53 2.50
25 3.95 2.47 2.75 3.31 2.60
30 2.45 2.38 3.22 2.39 4.62
40 3.64 3.29 3.28 2.56 4.99
50 3.21 2.42 3.31 3.42 3.70
T-REMD 15 0.31 0.30 1.51 0.79 0.91
20 0.65 0.39 0.51 0.51 0.69
25 0.31 0.31 0.22 0.36 0.44
30 0.19 0.23 0.30 0.45 0.33
40 0.13 0.21 0.22 0.28 0.44
50 0.13 0.12 0.10 0.10 0.22
Table 2  Maximum autocorrelation times (ps) obtained from normal MD and T-REMD simulations for different potential models.
MD REMD Theo.
T/K lnKs tc/ps T/K lnKs tc/ps lnKr
1000 1.52(0.12) 1000 1000 1.45(0.17) 250 1.41
2000 0.69(0.09) 100 2042 0.74(0.13) 250 0.69
3000 0.40(0.08) 30 2954 0.40(0.08) 250 0.33
4000 0.37(0.07) 10 4000 0.33(0.03) 250 0.26
Table 3  Logarithm equilibrium constants (lnKs) predicted using MD and T-REMD simulations and theoretical calculations (lnKr) at different temperatures (in Kelvin) for model (-3, 25).
Fig 4  Free energy contour maps obtained from BFMD (left) and T-REMD (right) for M(-3, 25) at 4000 K. The energy scale is in kcal•mol-1.
T/K T-REMD MD
1000 17.1 16.8
1937 22.8 23.1
2991 26.6 26.5
4000 30.8 30.9
Table 4  Free energy barriers of model (-3, 25) at different temperatures, obtained by T-REMD and normal MD simulation.
Fig 5  Discrepancies between predicted and theoretical equilibrium constants in natural logarithm for different potential models (∆Er, Ea) in kcal•mol-1.
Fig 6  Replica mixing ratio obtained in T-REMD simulations of different reactive models. The left panel shows the good cases, models (-3, 25) and (0, 40), the right panel shows the bad cases, models (-10, 25) and (-5, 40). The diagonal line represents the ideal case.
Fig 7  Sampling in phase space represented in reaction quotient (Q) and temperature (in K), obtained in 2-ns T-REMD simulations for models (-3, 15), (-3, 40), (-5, 15) and (-5, 40), respectively. The purple dots represent the equilibrium constants, the small dots in different colors are simulation data from cold (blue) to hot (red) temperatures. Color online.
i Ea(i)
1 15.0
2 16.1
3 17.8
4 19.9
5 21.6
6 23.7
7 26.0
8 28.5
9 31.2
10 34.0
11 36.9
12 40.0
Table 5  Optimized activation energies in models (-5, Ea(i)) used for the H-REMD simulations.
Fig 8  The energy overlaps of neighbored replicas of H-REMD (left) and the replica-diffusive ratio (right) obtained from H-REMD simulation of model (-5, 40), the diagonal line represents the ideal case.
Model T-REMD H-REMD Theo.
(-3, 25) 1.45 ± 0.17 1.44 ± 0.15 1.48
(-5, 40) 1.77 ± 0.21 2.40 ± 0.30 2.49
(-10, 50) 1.68 ± 0.20 4.90 ± 0.50 5.01
Table 6  Comparison of H-REMD and T-REMD predictions of lnK on different potential models.
Fig 9  Number of reactant molecules (AB) versus simulation time obtained from T-REMD (left) and H-REMD (right) simulations with different initial states: all AB (exothermic reaction) and all BC (endothermic reaction) comparing with theoretical values at equilibration for potential models (±5, 40). The data points are block-averaged in 25-ps bin.
1 Steinfeld, J. I. ; Francisco, J. S. ; Hase, W. L. Chemical Kinetics and Dynamics; Prentice Hall Englewood Cliffs: Upper Saddle River, NJ, USA, 1989; Vol. 3.
2 érdi, P. ; Tóth, J. Mathematical Models of Chemical Reactions: Theory and Applications of Deterministic and Stochastic Models; Manchester University Press: Oxford Road, Manchester M13 9PL, UK, 1989.
3 Kresse G. ; Hafner J. Phys. Rev. B 1993, 47, 558.
doi: 10.1103/PhysRevB.47.558
4 Van Duin A. C. T. ; Dasgupta S. ; Lorant F. ; Goddard Ⅲ W. A. J. Phys. Chem. A 2001, 105, 9396.
doi: 10.1021/jp004368u
5 Sugita Y. ; Okamoto Y. Chem. Phys. Lett. 1999, 314, 141.
doi: 10.1016/S0009-2614(99)01123-9
6 Hukushima K. ; Nemoto K. J. Phys. Soc. Jpn. 1996, 65, 1604.
doi: 10.1143/JPSJ.65.1604
7 Fukunishi H. ; Watanabe O. ; Takada S. J. Chem. Phys. 2002, 116, 9058.
doi: 10.1063/1.1472510
8 Itoh S. G. ; Okumura H. J. Comput. Chem. 2013, 34, 622.
doi: 10.1002/jcc.23167
9 Itoh S. G. ; Okumura H. ; Okamoto Y. J. Chem. Phys. 2010, 132, 134105.
doi: 10.1063/1.3372767
10 Swails J. M. ; Roitberg A. E. J. Chem. Theory Comput. 2012, 8, 4393.
doi: 10.1021/ct300512h
11 Mori T. ; Jung J. ; Sugita Y. J. Chem. Theory Comput. 2013, 9, 5629.
doi: 10.1021/ct400445k
12 Leahy C. T. ; Kells A. ; Hummer G. ; Buchete N. V. ; Rosta E. J. Chem. Phys. 2017, 147, 152725.
doi: 10.1063/1.5004774
13 Stelzl L. S. ; Hummer G. J. Chem. Theory Comput. 2017, 13, 3927.
doi: 10.1021/acs.jctc.7b00372
14 Wallace, A. F. Replica Exchange Methods in Biomineral Simulations. In Methods in Enzymology, De Yoreo, J. J., Ed. ; Academic Press: New York, NY, USA, 2013; Vol. 532, Chapter 4, pp. 71–93.
15 Zhang W. ; Chen J. J. Chem. Theory Comput. 2013, 9, 2849.
doi: 10.1021/ct400191b
16 Bergonzo C. ; Henriksen N. M. ; Roe D. R. ; Swails J. M. ; Roitberg A. E. ; Cheatham T. E. J. Chem. Theory Comput. 2014, 10, 492.
doi: 10.1021/ct400862k
17 Mori Y. ; Okamoto Y. Phys. Rev. E 2013, 87, 023301.
doi: 10.1103/PhysRevE.87.023301
18 Petraglia R. ; Nicola A. ; Wodrich M. D. ; Ceriotti M. ; Corminboeuf C. J. Comput. Chem. 2016, 37, 83.
doi: 10.1002/jcc.24025
19 Sato S. J. Chem. Phys. 1955, 23, 592.
doi: 10.1063/1.1742043
20 Plimpton S. J. Comput. Phys. 1995, 117, 1.
doi: 10.1006/jcph.1995.1039
21 Shinoda W. ; Shiga M. ; Mikami M. Phys. Rev. B 2004, 69, 134103.
doi: 10.1103/PhysRevB.69.134103
22 Denschlag R. ; Lingenheil M. ; Tavan P. Chem. Phys. Lett. 2009, 473, 193.
doi: 10.1016/j.cplett.2009.03.053
23 Bittner E. ; Nu?baumer A. ; Janke W. Phys. Rev. Lett. 2008, 101, 130603.
doi: 10.1103/PhysRevLett.101.130603
24 Plattner N. ; Doll J. D. ; Meuwly M. J. Chem. Theory Comput. 2013, 9, 4215.
doi: 10.1021/ct400355g
25 Dupuis P. ; Liu Y. ; Plattner N. ; Doll J. D. Multiscale Model. Simul. 2012, 10, 986.
doi: 10.1137/110853145
26 Wolff U. Comput. Phys. Commun. 2004, 156, 143.
doi: 10.1016/S0010-4655(03)00467-3
27 Flyvbjerg H. ; Petersen H. G. J. Chem. Phys. 1989, 91, 461.
doi: 10.1063/1.457480
28 Katzgraber H. G. ; Trebst S. ; Huse D. A. ; Troyer M. J. Stat. Mech. Theory Exp. 2006, 03, P03018.
doi: 10.1088/1742-5468/2006/03/P03018
[1] 郭俊江, 唐石云, 李瑞, 谈宁馨. 高碳烃宽温度范围燃烧机理构建及动力学模拟[J]. 物理化学学报, 2019, 35(2): 182-192.
[2] 杨化超, 薄拯, 帅骁睿, 严建华, 岑可法. 润湿特性对超级电容器储能动力学的影响机理[J]. 物理化学学报, 2019, 35(2): 200-207.
[3] 陈文琼,关永吉,张晓萍,邓友全. 分子动力学模拟研究外电场对咪唑类离子液体振动谱的影响[J]. 物理化学学报, 2018, 34(8): 912-919.
[4] 孙成珍,白博峰. 气体分子在二维石墨烯纳米孔中的选择性渗透特性[J]. 物理化学学报, 2018, 34(10): 1136-1143.
[5] 柳平英,刘春艳,刘倩,马晶. 含偶氮苯主-客体复合物的光致异构化反应对结合能与几何构象的影响[J]. 物理化学学报, 2018, 34(10): 1171-1178.
[6] 刘夫锋,范玉波,刘珍,白姝. ZAβ3和Aβ16-40亲和作用的分子机理解析[J]. 物理化学学报, 2017, 33(9): 1905-1914.
[7] 汪秀秀,赵健伟,余刚. 孔洞和孪晶界对银纳米线形变行为联合影响的分子动力学模拟[J]. 物理化学学报, 2017, 33(9): 1773-1780.
[8] 王子民,郑默,谢勇冰,李晓霞,曾鸣,曹宏斌,郭力. 基于ReaxFF力场的对硝基苯酚臭氧氧化分子动力学模拟[J]. 物理化学学报, 2017, 33(7): 1399-1410.
[9] 曹了然,张春煜,张鼎林,楚慧郢,张跃斌,李国辉. 分子动力学模拟技术在生物分子研究中的进展[J]. 物理化学学报, 2017, 33(7): 1354-1365.
[10] 陈芳,刘圆圆,王建龙,苏宁宁,李丽洁,陈红春. 混合溶剂对β-HMX结晶形貌影响的分子动力学模拟[J]. 物理化学学报, 2017, 33(6): 1140-1148.
[11] 陈贻建,周洪涛,葛际江,徐桂英. 双链阴离子表面活性剂1-烷基-癸基磺酸钠在气/液界面聚集行为:分子动力学模拟研究[J]. 物理化学学报, 2017, 33(6): 1214-1222.
[12] 周婷婷,宋华杰,黄风雷. 冲击载荷下TATB晶体滑移和各向异性的分子动力学研究[J]. 物理化学学报, 2017, 33(5): 949-959.
[13] 彭莉娟,姚倩,王静波,李泽荣,朱权,李象远. RDX及其衍生物高温热解的反应分子动力学模拟[J]. 物理化学学报, 2017, 33(4): 745-754.
[14] 刘青康,宋文平,黄其涛,张广玉,侯珍秀. 热辅助存储磁盘硅掺杂非晶碳薄膜氧化的ReaxFF反应力场分子动力学模拟[J]. 物理化学学报, 2017, 33(12): 2472-2479.
[15] 孙怡然,于飞,马杰. 纳米受限水的研究进展[J]. 物理化学学报, 2017, 33(11): 2173-2183.