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

研究者报道了一种名为Deep Potential-Range Correction (DPRc) 的新模型，可应用于溶液反应的QM/MM模拟中，有望为生命科学领域的自由能应用带来帮助。

基于分子力学（MM）的分子动力学模拟已在生命科学领域广泛应用，但从头算（ab initio）量子力学（QM）方法的应用仍然受到限制。基于ab initio 的QM/MM方法能大幅提升酶催化反应、药物发现或精准医学等领域的计算精度，但受到计算成本的制约。

近年来，DP已成功应用于多种体系，但应用于溶液体系时仍存在较大困难。近日，DeepModeling社区的研究者在Journal of Chemical Theory and Computation 上报道了 Deep Potential-Range Correction (DPRc) 模型^[1]，这种模型将 DP 作为校正项与 QM/MM 模型耦合，并用非磷酸转移酶在水溶液中的反应体系（图1）进行了测试。

图1 非磷酸转移酶的反应

基于DPRc校正的QM/MM能量如下式所示：

其中R是原子坐标，P是QM哈密顿量的单粒子密度矩阵，E_QM(R; P) 和E_MM(R)分别是QM能量和MM能量，E_QM/MM(R; P)包括QM区域和MM区域间的静电力和分子间作用，E_ML(R)是用DP训练的能量。将半经验的模型校正为从头算模型时，MM/MM 作用将保持不变；此外，为了保证能量守恒，QM/MM 作用应当随距离变化而过渡到0。因此DPRc模型将原始DP模型^[2]中平滑函数改为:

此外，DPRc模型避免了MM原子的单体能量E_i⁽⁰⁾：

其中单体能量可由下式得到：

利用改进后的DeePMD-kit^[3]和DP-GEN^[4]软件，研究者以MNDO/d和DFTB2两种半经验方法为例，在不同的截断半径和采样方法下，根据伞形抽样的模拟结果作出非磷酸转移酶体系在水溶液中的自由能曲线，与参考曲线进行比较（图2）。为了测试模型的迁移性，研究者还进一步研究了酶的不同变体的自由能曲线（图3），并用一维数据训练的模型，绘制了二维的自由能曲线（图4）。

图2 不同参数下的自由能曲线

图3 不同变体的自由能曲线。

其中，(b)和(e)没有包含在训练集中

图4 二维自由能曲线

训练数据来自一维反应坐标

结果表明，DPRc 模型在 QM/MM 作用的截断半径为6埃时可以达到较高的精度，高温数据有利于增加模型精度。此外，模型可以同时模拟4种不同变体的反应，并能良好地迁移到S1P变体及二维反应坐标中。该工作有望为药物发现和酶设计等自由能应用带来帮助。

DeePMD-kit教程1：

为不同区域设置不同的截断半径

以图1的反应为例，假设模型的type_map是C、H、HW、O、OW、P，其中HW和OW表示MM水分子中的原子。在DPRc模型中，我们需要为QM区域内的作用设置的无限大的截断半径，QM/MM作用的截断半径设置为6埃，并抹除MM区域内的相互作用，描述符可以表示为QM/QM和QM/MM的混合，关键参数如下所示：

"descriptor" :{              "type":             "hybrid",              "list" : [                  {                      "type":     "se_e2_a",                      "sel":              [6, 11, 0, 6, 0, 1],                      "rcut_smth":        40.,                      "rcut":             40.,                    "neuron":           [12, 25, 50],                    "exclude_types":    [[2,2], [2,4], [4,4], [0,2], [0,4], [1,2], [1,4], [3,2], [3,4], [5,2], [5,4]],                    "axis_neuron":      12,                      "set_davg_zero":    true,                      "_comment": " QM/QM interaction"                  },                  {                      "type":     "se_e2_a",                      "sel":              [6, 11, 100, 6, 50, 1],                      "rcut_smth":        0.50,                      "rcut":             6.00,                      "neuron":           [12, 25, 50],                      "exclude_types":    [[2,2], [2,4], [4,4], [0,0], [0,1], [0,3], [0,5], [1,1], [1,3], [1,5], [3,3], [3,5], [5,5]],                      "axis_neuron":      12,                      "set_davg_zero":    true,                      "_comment": " QM/MM interaction"                  }              ]          }

其中exclude_types表示排除的原子类型间的相互作用。因为体系中QM区域较小，仅有24个原子，因此在非周期性的数据中将rcut设置为一个比较大的数，即可确保截断半径实际上是无限的。

DeePMD-kit教程2：

指定原子的单体能量

DPRc需要指定MM区域原子的单体能量为0。当descriptor中的set_davg_zero设置为true时，fitting_net中的atom_ener参数可以设置为指定的单体能量（null表示不指定能量）：

"fitting_net": { "neuron": [240, 240, 240], "resnet_dt": true, "atom_ener": [null, null, 0., null, 0., null]}

具体参数说明，可以参见DeePMD-kit文档。点击“阅读原文”，阅读论文全文。

参考文献

1. Zeng et al. J. Chem. Theory Comput., 2021.

2. Zhang et al. NIPS’18, 4441-4451.

3. Wang et al. Comput. Phys. Commun., 2018, 228, 178-184.

4. Zhang et al. Comput. Phys. Commun., 2020, 253, 107206.

- End -

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

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

推荐阅读