DP还能干这个？深度剖析DP在溶液反应中的应用（二）

溶液反应是生命科学领域最主要的一类反应之一。在计算模拟中，量子力学/分子力学（QM/MM）方法能够在周期性边界条件下严格处理长程静电相互作用，因此被广泛应用于各类体系中。QM/MM方法将体系分为QM部分和MM部分，分别以精准的量子力学和快速的分子力学方法计算。

2013年，提出QM/MM方法的Martin Karplus等人，因“为复杂化学系统开发多尺度模型”，获得诺贝尔化学奖。然而，从头算（ab initio）QM/MM分子动力学（MD）模拟的计算成本仍然十分高昂，难以进行大量模拟，基于半经验QM方法的QM/MM MD又缺乏准确性。

去年10月，罗格斯大学Darrin M. York研究组基于现有的深度势能模型，结合QM/MM的原理和特点，创造性地开发了Deep Potential - Range correction（DPRc）模型，并应用于水溶液中RNA的非酶磷酸酯交换反应（图1）的研究中。该项研究发表于J. Chem. Theory Comput.期刊[1]，我们已在《DP在溶液反应中的应用（一）》一文中详细介绍。DPRc模型旨在将低精度QM/MM校正到高精度QM/MM，它的特点是：只校正QM-QM相互作用和QM-MM相互作用，但不校正MM-MM的相互作用，也不校正MM区域的原子能量；在QM区域外的截断半径内，模型校正的能量平滑过渡到0。

图1 水溶液中的非酶磷酸酯交换反应，X是氧原子，但亦有五种不同的硫代变体

昨日，同班人马接续之前的研究，再次在J. Chem. Theory Comput.期刊报道了他们的最新研究成果[2]。新的研究以同一反应体系（图1）为案例，训练了DPRc模型，将DFTB2 QM/MM[6]校正至PBE0/6-31G* QM/MM。基于DPRc模型，研究比较了DFTB2 QM/MM+DPRc与PBE0/6-31G* QM/MM的反应自由能曲线，并将DPRc模型与路径积分动力学（PIMD）模拟相结合，研究量子效应对该反应自由能曲线的影响，以及¹⁸O和³²S动力学同位素效应（KIE）对反应的影响。

如图2所示，研究者首先绘制了PBE0/6-31G* QM/MM（参考值，黑线），DFTB2 QM/MM+DPRc（红线），DFTB2 QM/MM（绿线）三种方法的反应自由能曲线。在原生（Native）反应和五种硫代变体的反应中，DFTB2 QM/MM与PBE0/6-31G* QM/MM相差甚远，而DFTB2 QM/MM+DPRc可以较为完美地复现PBE0/6-31G* QM/MM的自由能曲线，证明DPRc模型相比单纯的半经验方法具有较高的精度。

图2 PBE0/6-31G* QM/MM（黑线），DFTB2 QM/MM+DPRc（红线），DFTB2 QM/MM（绿线）三种方法的自由能曲线

DPRc模型相较于自由能扰动重加权方法的精度如何？研究者用加权热力学扰动（weighted Thermodynamic Perturbation，wTP）方法[3]，试图从DFTB2模拟中复现PBE0/6-31G*曲线，如图3所示，显然，在多个反应中，DPRc方法的精度（图2红线）均显著优于wTP方法（图3绿线），这再次坚定了研究者进一步发展DPRc方法的信心。

图3 PBE0/6-31G* QM/MM（黑线）、wTP（红点）、高斯过程回归（GPR）平滑wTP（绿线）的自由能曲线

利用相同的DFTB2 QM/MM+DPRc模型，研究者分别进行了经典MD和PIMD模拟，并比较了它们的自由能曲线。如图4所示，在量子效应下，PIMD的自由能能垒比经典MD计算的能垒平均小1.48 kcal/mol（native最大，S3'最小）。也就是说，比较反应物和过渡态震动环境的差异，过渡态总体上更“松散”（即键长较长），导致了较低的能垒。

图4 经典MD和PIMD结果比较

结合PIMD，研究者分别利用TD-FEP和Bigeleisen两种计算KIE的方法，计算了X2'和X5'原子的KIE，如表1所示。两种方法得到的KIE类似，相差不超过0.003。

表1 六种变体中X2'和X5'原子的KIE

接下来，研究者借鉴KIE的概念，与自由能面（FES）相结合，创造了沿反应坐标的同位素效应（IE）的概念。沿自由能面，将采样点视为“过渡态”，FES的最小处视为反应物，用PIMD TD-FEP进行模拟，求得IE。如图5所示，研究者绘制了O2'和O5'原子沿反应坐标的同位素效应（黑线）和自由能曲线（绿线），从而研究自由能曲线上沿反应坐标的同位素效应。可以看出，当成键环境较为“松散”时，同位素效应取得最大值，反应了同位素效应和过渡态成键结构的相关性。这与图4的结果是一致的。

图5 O2'和O5'原子沿反应坐标的同位素效应（IE，黑点、黑线）、与磷原子的距离（红线）和反应自由能曲线（FES，绿线，即图4a红线）

计算效率上，仅使用1核心英特尔至强CPU时，DFTB2+DPRc比PBE0/6-31G*快100倍；在此基础上，如果使用英伟达V100仅加速DPRc部分，则DFTB2+DPRc比PBE0/6-31G*快250倍。DPRc模型在精度逼近的情况下大大加快了计算效率，与PIMD相结合后，有望成为研究溶液中反应机制、通过KIE研究过渡态结构的强大工具。

方法

初始结构：研究者用OpenBabel软件[4]生成了图1分子的初始结构，并加入了1510个TIP4P/Ew水分子，用AMBER20软件[5]进行100ps的DFTB2 QM/MM[6]（MIO参数，下同） 的0 ~ 298 K升温模拟，得到298K下的初始结构。

ab initio QM/MM模拟：沿图1所示的反应坐标ξ将反应分为91个窗口，进行25 ps/window的PBE0/6-31G* QM/MM的伞形抽样模拟，得到了初始训练集。对于每个伞形抽样模拟，用FE-ToolKit软件[7-8]得到了相应反应坐标下的自由能曲线（下同）。

训练DPRc模型：利用DP-GEN软件[9]，进行数轮迭代，生成数据。每轮分为训练、探索和标签三步，训练步用DeePMD-kit软件[10]训练并压缩[11]DPRc模型，探索步用具有DeePMD-kit接口的AMBER20软件[1,5,10,17]进行25 ps/window的DFTB2 QM/MM+DPRc下的伞形抽样模拟，标签步选取若干最大模型偏差为0.1-0.25 eV/Å的结构，计算DFTB2 QM/MM和PBE0/6-31G* QM/MM的能量差。当最大模型偏差小于0.1 eV/Å的结构足够多时，迭代停止，得到最终的DPRc模型。

DFTB2 QM/MM+DPRc PIMD模拟：DFTB2 QM/MM+DPRc PIMD通过i-PI[12]、AMBER20、PLUMED[13]、DeePMD-kit[10]等四个软件联合完成，其中，i-PI用于PIMD，PLUMED用于伞形抽样，六个AMBER20（调用DeePMD-kit）进程并行进行DFTB2 QM/MM+DPRc计算。PIMD使用PIGLET恒温器，其参数由GLE Input Generator[14]网页生成。使用SPC/Fw水模型，对每个窗口进行16次20 ps的PIMD模拟，并在相同参数下进行经典MD模拟以进行比较。

KIE计算：可用两种方法计算O5'和O2'原子的KIE。第一种方法为Bigeleisen-Mayer方程[15]，KIE从几何静止点处的振动分析计算得到：

式中假设共有(3N-6)个自由度，   和  分别是反应物中重原子和轻原子的振动频率，  和  分别是过渡态中重原子和轻原子的振动频率,  和  分别是它们的唯一虚频（虚频的sinh值被舍去），这些频率用DL-Find软件[16]和LBFGS（Limited-memory Broyden–Fletcher–Goldfarb–Shanno）算法，对自由能最小处及过渡态的结构进行优化，计算振动频率。

第二种方法由下式给出：

其中，  和  已在法一中得到；  和  是将轻原子改为重原子时，过渡态和反应物的自由能变化，对此，对过渡态和反应物分别进行四次10 ps的DFTB2 QM/MM+DPRc PIMD模拟，得到平均量子动能，并用i-PI内置的热力学自由能扰动（TD-FEP）方法计算  和  。

代码可用性

出于对开源的支持，此项研究中，具有DeePMD-kit接口的AMBER20已在https://gitlab.com/RutgersLBSR/AmberDPRc/ 开源，对DeePMD-kit和DP-GEN软件的DPRc支持也已合入主线开发分支。至此，此项研究中使用的9项软件均为开源软件，只需学习这些软件的使用方法，即可完整复现结果。

DOI: 10.1021/acs.jctc.2c00151

参考文献

[1] Zeng, J.; Giese, T. J.; Ekesan, ¸S.; York, D. M. Development of Range-Corrected Deep Learning Potentials for Fast, Accurate Quantum Mechanical/Molecular Mechanical Simulations of Chemical Reactions in Solution. J. Chem. Theory Comput. 2021, 17, 6993–7009.

[2] Giese, T. J.; Zeng, J.; Ekesan, ¸S.; York, D. M. Combined QM/MM, machine learning path integral approach to compute free energy profiles and kinetic isotope effects in RNA cleavage reactions. J. Chem. Theory Comput. 2022.

[3] Hu,W.; Li, P.;Wang, J.-N.; Xue, Y.; Mo, Y.; Zheng, J.; Pan, X.; Shao, Y.; Mei, Y. Accelerated Computation of Free Energy Profile at Ab Initio Quantum Mechanical/Molecular Mechanics Accuracy via a Semiempirical Reference Potential. 3. Gaussian Smoothing on Density-of-States. J. Chem. Theory Comput. 2020, 16, 6814–6822.

[4] O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R.

Open Babel: An open chemical toolbox. J Cheminform 2011, 3, 33.

[5] Case, D. A.; Belfon, K.; Ben-Shalom, I. Y.; Brozell, S. R.; Cerutti, D. S.; Cheatham III, T. E.; Cruzeiro, V. W. D.; Darden, T. A.; Duke, R. E.; Giambasu, G.; ; Gilson, M. K.; Gohlke, H.; Goetz, A. W.; Harris, R.; Izadi, S.; Izmailov, S. A.; Kasavajhala, K.; Kovalenko, K.; Krasny, R.; Kurtzman, T.; Lee, T.; Le-Grand, S.; Li, P.; Lin, C.; Liu, J.; Luchko, T.; Luo, R.; Man, V.; Merz, K.; Miao, Y.; Mikhailovskii, O.; Monard, G.; ; Nguyen, H.; Onufriev, A.; Pan, F.; Pantano, S.; Qi, R.; Roe, D. R.; Roitberg, A.; Sagui, C.; Schott-Verdugo, S.; Shen, J.; Simmerling, C. L.; Skrynnikov, N.; Smith, J.; Swails, J.; Walker, R. C.; Wang, J.; Wilson, R. M.; Wolf, R. M.; Wu, X.; Xiong, Y.; Xue, Y.; York, D. M.; Kollman, P. A. AMBER 20. University of California, San Francisco: San Francisco, CA, 2020.

[6] Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys. Rev. B 1998, 58, 7260–7268.

[7] Giese, T. J.; York, D. M. FE-ToolKit: The free energy analysis toolkit. https://gitlab.com/RutgersLBSR/fe-toolkit.

[8] Giese, T. J.; Ekesan, ¸S.; York, D. M. Extension of the Variational Free Energy Profile and Multistate Bennett Acceptance RatioMethods for High-Dimensional Potential ofMean Force Profile Analysis. J. Phys. Chem. A 2021, 125, 4216–4232.

[9] Zhang, Y.; Wang, H.; Chen, W.; Zeng, J.; Zhang, L.; Han, W.; E, W. DP-GEN: A concurrent learning platform for the generation of reliable deep learning based potential energy models. Comput. Phys. Commun. 2020, 253, 107206.

[10] Wang, H.; Zhang, L.; Han, J.; E, W. DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics. Comput. Phys. Commun. 2018, 228, 178–184.

[11] Lu, D.; Jiang, W.; Chen, Y.; Zhang, L.; Jia, W.; Wang, H.; Chen, M. DP Train, then DP Compress: Model Compression in Deep Potential Molecular Dynamics. 2021.

[12] Ceriotti,M.; More, J.; Manolopoulos, D. E. i-PI: A Python interface for ab initio path integral molecular dynamics simulations. Comput. Phys. Commun. 2014, 185, 1019–1026.

[13] Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R.; Parrinello, M. PLUMED: a portable plugin for free energy calculations with molecular dynamics. Comput. Phys. Commun. 2009, 180, 1961–1972.

[14] Ceriotti,M.; Bussi, G.; Parrinello,M. Colored-Noise Thermostats à la Carte. J. Chem. Theory Comput. 2010, 6, 1170–1180.

[15] Bigeleisen, J.; Mayer, M. G. Calculation of Equilibrium Constants for Isotopic Exchange Reactions . J. Chem. Phys. 1947, 15, 261–267.

[16] Kästner, J.; Carr, J.M.; Keal, T.W.; Thiel,W.;Wander, A.; Sherwood, P. DL-FIND: An opensource geomtery optimizer for atomistic simulations. J. Phys. Chem. A 2009, 113, 11856–11865.

[17] https://gitlab.com/RutgersLBSR/AmberDPRc/

- End -

（如需转载图文请与公众号后台联系）

-------------------------------

推荐阅读