JACS：准确高效考虑核量子效应的分子动力学新方法

量子力学自诞生至今已有百年的历史，然而在化学和生物学领域当中，经典的牛顿力学图像依然被广泛使用。

这一方面当然是因为经典力学更贴近我们的生活直觉，不过另一方面，由于质量的显著不同，原子核的量子效应要比电子小得多，因此，我们往往只需要把电子用量子力学来处理，而原子核则可以只用经典力学来处理。

图1. 水分子的经典振动模式

可是在很多体系中，人们又发现原子核的量子效应不能被忽略。

一个有趣的例子是，重水尝起来是甜的。这被认为是当氘取代氕之后，重水分子与我们人体内甜味受体的结合方式与普通的水分子不再相同，我们因而感受到了甜味。

另一个例子是，假如不存在核的量子效应，水中的氢与氧会结合得更紧密，从而使得水的pH值会在8以上。当然，这并不是说我们的水就变成了碱性，而是水中氢离子和氢氧根离子的浓度都会减少一个数量级。

事实上，像氢核这样比较轻的原子核，其量子效应往往是非常显著的，显著到足以改变体系的物理和化学性质。

我们知道液态水和很多重要的生物分子中都包含大量的氢原子与氢键，当我们模拟这样的体系时，简单地把原子核用经典力学处理会产生一些偏差。这也就导致即使是从头算分子动力学（ab initio molecular dynamics, AIMD）在处理氢原子相关的过程时也会有显著的误差。

在这个时候，一个常用的工具是基于费曼路径积分的分子动力学 (比如path-integral molecular dynamics, ring-polymer molecular dynamics, centroid molecular dynamics)。可惜的是，这些基于路径积分的方法往往需要大量的时间和存储资源，尤其是当研究动力学性质的时候。

图2. 水分子中氢原子核的经典处理(a)与量子处理(b)示意图

近日，来自威斯康星大学麦迪逊分校理论化学研究所的杨阳团队开发了一种新的理论方法，可以较为准确而且高效地在分子动力学模拟中包含原子核的量子效应。这种方法的基本思想是：虽然原子核由于量子效应不应该再用点电荷来处理，但是和电子相比，每个原子核的密度分布仍然是相对比较局域的，因此原子核的位置仍然可以用它的位置期望值来表示。在这个假设下，作者基于多组分密度泛函理论，变分求解特定原子核期望值下的原子核的密度和电子密度，进而得到一个基于原子核位置期望值的能量面。这个新的能量面不同于过去传统基于电子结构计算的势能面——它包含了原子核的量子效应。经过一系列数学推导，作者得到一套与牛顿力学非常相似的分子动力学方程，区别是能量和力的计算要基于这种新的含原子核量子效应的能量面。作者把这种新的分子动力学方法称“cNEO-MD”(constrained nuclear-electronic orbital molecular dynamics)。

图3. 水分子及其氘代分子的红外振动谱

在此基础上，作者计算了一系列小分子的振动谱。图3展示了水分子及其氘代分子的红外谱，左边为AIMD的计算结果，右边为cNEO-MD的计算结果，灰色的虚线代表实验值作为对比。对于三种水分子高频区域的两个伸缩振动模式，AIMD平均高估了约100波数，而cNEO-MD的平均误差能减小一个数量级至约为17波数。另外一点值得注意的是，cNEO-MD的结果中HDO分子在2700到2800波数之间的振动峰发生了裂分，这是因为HDO分子弯曲振动的倍频和OD伸缩的基频都出现在这附近，二者发生了费米共振。实验上以及cNEO-MD的计算中都能够观察到这个裂分，然而AIMD由于对振动频率的预测不够准确，不能够观察到这个费米共振。

图4. 正丁醇分子的红外谱

除水分子之外，作者也计算了甲酸、甲醇和正丁醇等分子来测试cNEO-MD。其中正丁醇分子的红外谱如图4所示，从上往下依次为AIMD，cNEO-MD和实验结果。与实验结果相比，AIMD对于高频区O-H和C-H伸缩振动的频率有着明显的高估，而cNEO-MD能够非常准确地描述这些振动模式的频率，而且峰强和谱型也和实验谱基本一致。由此我们可以看到cNEO-MD的准确性，以及在这些氢原子参与的振动模式中核量子效应的重要性。

我们可能会有这样的疑虑，把原子核也做量子处理会不会带来难以负担的计算开销呢？毕竟单是用量子的方式处理电子已经够麻烦了。关于这个问题，作者对cNEO-MD的计算开销做了基准测试。他们发现对于论文中计算的分子，虽然所有的氢原子核都是量子处理的，cNEO-MD的计算量仅为过去AIMD的1.5-3倍。对于这样具有较高准确性的方法，这个计算效率已经很可观了，而考虑到未来其与力场开发以及机器学习结合的可能性，我们有理由期待这种分子动力学方法未来能够被应用到更多的化学和生物体系中的模拟当中。

这一成果近日发表在J. Am. Chem. Soc.上，通讯作者为威斯康星大学麦迪逊分校化学系与理论化学研究所的杨阳，文章的第一作者许熙今年年中将加入北京师范大学珠海校区先进材料研究中心。论文中计算所使用的代码在开源量子化学程序PySCF的基础上完成，源码可以通过GitHub下载  ^[1]。

原文（扫描或长按二维码，识别后直达原文页面，或点此查看原文）：

Molecular Dynamics with Constrained Nuclear Electronic Orbital Density Functional Theory: Accurate Vibrational Spectra from Efficient Incorporation of Nuclear Quantum Effects

Xi Xu, Zehua Chen, and Yang Yang*

J. Am. Chem. Soc., 2022, DOI: 10.1021/jacs.1c12932

杨阳

https://yang.chem.wisc.edu 

参考资料：

1. PySCF

https://github.com/theorychemyang/pyscf