摘要

分子模拟（Molecular Simulation）是一种利用计算机以原子水平的分子模型来模拟分子结构与行为，进而计算分子体系的各种物理、化学性质的方法。通过该方法，研究人员可以在微观尺度下预测实验无法检测的关键结构与性质、探究物理化学过程背后的机理以及有效降低实验次数和成本。因此，分子模拟技术目前在智能制药、先进材料、生物科技、石油化工等领域都有非常广泛的应用，是科学计算领域中一个非常重要的方向。作为系列文章的第一篇，本文将从源头出发，详细阐述分子模拟技术的发展历程和现状，同时揭示其底层机理。

分子模拟技术在智能制药、先进材料、生物科技、石油化工等领域都有非常广泛的应用。2017年的一份统计报告指出，世界上最常用的15款高性能计算软件中有超过一半都属于分子模拟领域[1]。图1中展示的是两个典型的分子模拟研究案例，分别为生物医药领域的新冠药物发现研究和先进材料领域的碳纳米管材料研发工作。而进入二十一世纪以来，计算能力的大幅提升使得分子模拟的速度和准确度都有着显著的增加。可以预见的是，在未来的工业界，计算化学与分子模拟将会发挥更重要的作用。

图1：分子模拟研究案例：(a) 新冠治疗药物筛选 (b) 碳纳米管材料性能研究
来源：https://developer.nvidia.com/blog/molecular-dynamics-gpu-applications-help-fight-covid-19-via-ngc/

从本质上来说，分子模拟研究的目标是通过一个体系的化学组成和其所处的环境来预测其物理/化学以及动力学性质。如图2所示，一个完整的分子模拟工作流包括分子建模、构建力场、分子模拟、结果分析等几大部分[2]。根据模拟过程的时间空间尺度以及精确度需求，可以大致将分子模拟分成第一性原理模拟和经验力场模型模拟两大类。图3展示了计算机模拟仿真研究的各种类别[3]。图中的横轴为空间尺度，范围覆盖皮米(picometer,  10^-12米)到米；纵轴则为时间尺度，覆盖飞秒(femtosecond, 10^-15秒)至小时的范围。从图中的左下角到右上角，随着体系尺寸和模拟时间的增加，计算过程中所需要的精确度随之下降，而模拟方法从基于量子力学的第一性原理模拟转变为基于实验数据的经验力场模型模拟方法。当体系的时空尺寸进一步增加到宏观尺度以后，一般需要将流体粒子看作连续的介质来进行模拟计算，因此这部分的模拟研究一般采用基于纳维-斯托克斯(N-S)方程的计算流体力学模拟方法(Computational Fluid Dynamics, CFD)。CFD也是计算模拟仿真领域中很重要的一部分，关于CFD与GPU的结合在往期智算芯闻中有过详细介绍，感兴趣的读者可以点击文末链接前往。本文将重点介绍微观尺度下的分子模拟方法。

图2：分子模拟流程中的关键步骤
来源：Gartner and Jayaraman, “Modeling and Simulations of Polymers: A Roadmap”, Macromolecules 2019, 52, 755−786

图3：计算机模拟方法分类
来源：Keith etc. “Combining Machine Learning and Computational Chemistry for Predictive Insights into Chemical Systems”, Chem. Rev. 2021, 121, 9816−9872

作为一种具有大数量吞吐量、大规模并行计算特点的科学计算应用，分子模拟程序很适合通过GPU来进行加速。目前大部分主流的分子模拟应用软件中都包含GPU加速功能，在大规模体系的模拟中可以获得巨大的收益。当然，不同软件的GPU加速功能的实现方式和完备程度也不尽相同，其中包括仅支持将部分计算过程转移到GPU上计算的LAMMPS[4]等,支持全流程GPU计算的GROMACS[5]， NAMD[6]等, 以及完全基于GPU重新设计数据结构和算法的新型模拟引擎，如OpenMM[7], HOOMD-blue[8]等。本系列将分上下两篇，上篇详细梳理分子模拟技术的起源发展和技术细节，下篇则重点介绍分子模拟中的GPU加速应用现状及未来发展方向。

分子模拟技术的发展历程

01 第一性原理分子模拟

前文中提到，根据底层计算原理来区分，分子模拟技术可以被归为两大类别：第一性原理模拟和经验力场模型模拟。本小节首先介绍第一性原理分子模拟的方法的起源与发展。第一性原理这个词翻译自拉丁语 ab initio, 最初是由亚里士多德(Aristotle)提出的一个哲学理论，其寓意为从事物的本质出发来分析和解决问题。而在计算化学领域中，第一性原理计算又称从头算，指从最基本的物理学定律出发，不外加假设与经验拟合的推导与计算。也就是说，该方法可以在不需要任何实验参数的情况下，仅通过原子结构和物理常量来精准计算化学体系的各种性质。如何做到这一点呢？那就得归功于二十世纪最伟大的物理理论——量子力学。量子力学相信大家都听说过，可以说是如雷贯耳，但究竟什么是量子力学？著名物理学家理查德·费曼(Richard Feynman)曾经说过，“没人真正理解量子力学。” 一语成谶，时至今日，关于量子力学的原理也没有一个确定的答案，甚至目前还在研究量子力学底层原理的物理学家都已经是凤毛麟角了。但神奇的是，这一套没有人真正理解的理论，却可以通过简洁的公式精准预测微观粒子的各种行为，成为现代科技的基石。接下来我会花简短的篇幅来介绍一下这个颠覆所有人认知的理论的发展历程。

谈到量子力学的起源，那就不得不提及那个改变物理学发展轨迹的经典实验：双缝干涉实验。双缝干涉实验是由托马斯·杨(Thomas Young)在1801年提出的一个验证光的波动性的实验。实验本身非常简单：将一束光照向两个狭缝，然后在后面的探测屏上观察图样。在当时，主流的观点是光是由确定的粒子组成的。假如光是粒子，则探测屏上观查到的图样应该是两条亮纹。但在实际实验中，人们观察到的却是一种明暗相间的干涉条纹，并且条纹亮度可以用波动理论完美解释。这个实验结果促使光的波动说被广泛接受，并在此后的100多年里成为了主流。然而，在二十世纪初，一位叫爱因斯坦(Einstein)的德国人使用光的量子假说完美解释了光电效应，有力地证明了光也是一种粒子，他也因此获得了1921年的诺贝尔奖。这让人们又产生了疑惑，光到底是粒子还是波呢？因此有人改进了双缝实验来寻找答案。如图4所示，人们通过在狭缝后装置探测器来检测光子通过哪一条狭缝，或者其是否同时通过两条狭缝，以确定光的物理本质。此时令人震惊的实验结果发生了：在进行探测之后，探测屏上原本的干涉图样会完全消失，替代显示出的是两个单缝图样的简单总和。也就是说，探测这一行为竟然改变了光的物理实在。这一发现完全打破了经典物理学对世界本质的理解，在当时的物理界掀起了千层浪，甚至引发了唯心/唯物论的哲学讨论。更为神奇的是，埃尔文·薛定谔(Erwin Schrödinger)于1926年提出的波动方程可以完美预测量子的这些奇特行为。然而，如何解释这个波函数方程仍然是未解之谜。图5展示了1927年索尔维(Solvay)会议的合照，该会议是量子力学发展史上的一次重要事件，几乎所有当时世界上最优秀的科学家都参加了会议。为了解释薛定谔方程，物理学家们提出过很多假说，其中最有名的当属由玻尔(Bohr)、玻恩(Born)和海森堡(Heisenberg)等人提出的哥本哈根诠释。该理论认为量子体系波函数的真正含义是概率，一个事件发生的概率应该是波函数的绝对值平方。而测量的动作则会造成波函数塌缩，将原本的量子态按照机率塌缩成一个确定的状态。由此可以得到的推论是：当我们不观测宇宙时，宇宙就是一团概率波的叠加，而当我们观测宇宙时，宇宙才会出现一个确定的状态。听起来很荒谬？没错。然而这套理论可以解释和预测大部分(不是所有)量子现象，目前仍然是量子世界的主流诠释。可以想象，当时这个颠覆性的理论必然有着许多反对者，其中就包括伟大的爱因斯坦和量子波动方程之父薛定谔。爱因斯坦坚信世界是确定存在的，曾说过“上帝不掷色子”来反驳概率诠释。（玻尔的回应也很有趣，“爱因斯坦，别试图告诉上帝他应该做什么。”） 而薛定谔同样提出了著名的思想实验，“薛定谔的猫”，试图证明概率诠释的荒谬。决定论和概率论之间的论战是一段非常有趣的历史，因篇幅有限这里就不展开了，有兴趣的读者们可以自行百度。

图4：双缝干涉实验。左图为开启探测器的实验结果，右图为关闭探测器的实验结果。
来源：My State Of Mind by viewzone's editor Gary Vey

图5：1927年索尔维会议合影
来源：Conferencia de Solvay 1927 (1ª parte) - The Big Bang Physics

言归正传，第一性原理模拟的本质就是通过求解薛定谔方程来得到体系的能量和力场，从而预测粒子的运动轨迹。薛定谔方程的形式非常简洁：Ĥψ=Eψ，其中 ψ 为体系的波函数，Ĥ为表征波函数总能量的哈密顿算符，E 为体系的能量。那为什么这个方程这么难求解呢？因为波函数 Ψ(r₁, r₂, …, r_N) 是由体系内所有粒子之间的距离同时决定的，其维度为3N。事实上，只有最简单的双粒子体系（氢原子、氦离子）才有准确的解析解。而只要体系中的粒子总数达到三个或三个以上，我们就必须使用数值方法来求解。薛定谔方程的精确求解所需要的计算量是非常惊人的，目前最精确的完全组态相互作用（Full Configuration Interaction, FCL）方法[9]的时间复杂度为 O(N!) , 其中 N 为体系中的粒子（包括原子核和电子）总数。而这只是计算单个构象的能量所需的计算量，在进行分子动力学（molecular dynamics）模拟的过程中，每一步模拟都需要计算体系的能量和作用力。一步模拟的步长一般为1飞秒（femtosecond, 10^-15秒），通常一个有意义的化学过程发生的时间尺度至少需要达到纳秒（nanosecond, 10^-9秒）级别，而一些生物学过程，如蛋白质折叠，甚至需要几秒才能完成。因此，在当前的计算能力下，想要准确地对一个大体系的化学过程进行完全精确的模拟是不可能的。正如保罗·狄拉克（Paul Dirac）所说，“我们已经完全了解了物理和化学世界的基础定律，剩下的难点就是如何去求解这个定律推导出来的复杂方程。” 一直以来，量子化学领域的科学家们都在致力于通过近似和优化的方法来缩短求解薛定谔方程所需的时间，也取得了很多重要的成果。表1中总结了一些比较有名的方法以及其时间复杂度。其中最常用的哈特里-福克（Hartree-Fock, 简称HF）方法[10]通过自洽场求解的算法将时间标度降低到了 O(N⁴), 与FCL方法相比有了很大的提升。当然，需要指出的是，该方法做了很多假设，并且完全忽略了电子相关能，因此其精度要低于FCL方法。

表1中介绍的方法虽然各有不同，但其本质都是基于体系波函数来求解薛定谔方程，因此同属于波函数方法。而1965年由华裔科学家沈吕九（Sham）和其导师科恩（Kohn, 1988诺贝尔奖获得者）提出的密度泛函（Density Functional Theory, DFT）方法[11]则完全打破了使用波函数作为基础变量的常规思维。DFT方法的核心思想是使用电子密度取代波函数作为研究的基本量。因为波函数 Ψ(r₁, r₂, …, r_N) 是一个 3N 维度的函数，而电子密度泛函 ρ(r) 则只包含3个变量。很明显，用 ρ 取代 Ψ 可以极大地简化计算过程。这种变量替代处理的正确性可以在数学上得到严谨的证明，即体系的能量 E(ρ) 是一个仅仅关于电子密度的泛函。通过DFT自洽场的算法，求解过程的时间复杂度可降低到 O(N³) 。目前的一些最新算法可以达到 O(N^1.5) ，甚至在某些特殊情况下可以做到线性标度 O(N) 。因为其计算效率高、算法直接、易编程等优点，DFT方法成为了当今最常用的第一性原理计算方法。

表1：基于波函数的第一性原理求解方法小结

02 经验力场模型分子模拟

上一部分我们介绍了基于量子力学的第一性原理模拟方法。由于第一性原理计算的成本非常昂贵，通常情况下该方法只能模拟上至几千个原子的体系。但生物、材料研究中涉及的体系通常都远远大于这个限制, 动辄达到数万乃至数百万级原子规模。因此，对于这些大规模的体系，我们一般使用经验力场模型方法来进行模拟研究。该方法与第一性原理模拟方法的区别在于：其通过原子间势函数（interatomic potential）来估算体系中分子和原子的相互作用力，而不是直接求解薛定谔方程来精确计算体系力场。原子间势函数的种类有很多，但几乎都是通过对体系的实验数据和量子力学计算结果进行拟合得到的。势函数 E(r) 作为一个从原子空间坐标到体系总能量的映射，一般而言由几种能量组成：化学键伸缩势、化学键角弯曲势、二面角扭曲势以及非键相互作用能如范德华势、库仑势等。传统的势函数通常只适用于体系的物理变化过程模拟，但近年来出现的一些新型势函数也可以对一些特定的化学反应过程进行拟合。表2中总结了一些常见的势函数及其应用范围[12-17]。由于引入了外部数据，经验力场模型模拟方法不再属于第一性原理的范畴，其精度也有明显下降。然而，该方法带来的计算速度和效率提升是非常显著的，很适合进行大体系下的模拟研究。另外，这种力场计算方法的局部特性使其具有较好的拓展性，比较容易实现大规模并行计算，也成为了通用GPU在科学计算领域中的一个重要应用方向。

简单来说，基于经验力场模型的分子动力学模拟过程可以被分解为以下几个步骤：Ⅰ. 根据体系中每个原子的空间坐标，通过原子间势函数计算出体系的能量与每个原子所受的作用力；Ⅱ. 根据原子的初始位置以及上一步中计算出的作用力，通过积分的方式计算出经过一个较短的时间间隔(步长，通常为1飞秒)以后的原子坐标以及速度；Ⅲ.得到新的原子坐标数据以后，重复Ⅰ和Ⅱ的计算过程，完成下一步模拟。通常对一次完整的物理化学过程进行仿真需要通过许多步模拟。模拟过程中首先需要让体系进入平衡状态，之后才能收集可靠的数据进行处理分析。

表2：常见的势函数

03 机器学习力场模型模拟方法

前文中我们介绍了两种经典的分子模拟方法。总的来说，第一性原理模拟的精度很高，但是计算速度极慢；相反的，经验力场模型分子模拟的特点则是速度快、模拟效率高，但是精度较低。因此，这两种方法中间存在着一个巨大的鸿沟，即快与准不可兼得。近年来，随着机器学习和人工智能技术的飞速发展，在分子模拟领域也出现了很多关于机器学习力场模型（Machine Learning Potentials,简称MLP）的研究[18]。MLP可以直接学习从输入（体系化学组成，原子坐标等）到输出（体系性质）的关系，从而避免复杂的量子力学计算过程。因此，MLP方法具备同时达到高精度和高效率的潜力。如图6所示，许多学者认为MLP方法可以填补第一性原理与经验力场模拟方法之间的空白[19]。MLP研究中比较有代表性的工作包括鄂维南院士团队设计的深势（Deep Potential，简称DP）模型[20]。DP模型通过巧妙的方法构建了体系的描述符（descriptor），使其在保持简单的形式的同时使体系能量满足平移、旋转和置换对称性，从而能更准确地描述体系的特征。同时他们也改进了神经网络结构中的损失函数（loss function）,在其中加入体系能量、体系力场等重要数值，从而使得训练出的模型具备更强的预测能力。经过对水分子第一性原理计算数据的学习，DP模型成功在超级计算机上完成了超一亿水分子体系的分子动力学模拟[21]，显著拓展了分子模拟的体系尺寸边界（见图7），也因此获得了2020年度的戈登贝尔奖（Gorden Bell Prize）。当然，最理想化的MLP模型需要具备学习底层量子力学机理的能力，即模型需要具备可解释性。但是这一点目前在大部分机器学习模型中都还难以实现。另外，目前MLP方法还存在着模型局限性高、可靠训练数据集匮乏、长程力处理困难等问题。尽管如此，作为一种新兴方法，机器学习分子力场方法还是给分子模拟领域打开了一个全新的视角，即AI+高性能计算+传统科学计算的融合，这也将是科学计算领域未来的一个重要发展方向。

图6：机器学习方法的精度与速度介于第一性原理与经验力场方法之间
来源：Unke etc. “Machine Learning Force Fields”, Chem. Rev. 2021, 121, 16, 10142–10186

图7：深势模型成功拓展分子模拟体系边界
来源：深势科技官网

小结

分子模拟作为一种具有悠久历史的科学计算方法，长期以来在先进材料、生物科技以及药物发现等领域发挥着重要的作用。自诞生以来，分子模拟方法就一直面临着精度和效率不可兼得的困局。尽管计算机技术的飞速发展带来了算力的大幅提升，但是目前人们在面对传统第一性原理计算中的“维度灾难”问题时还是显得束手无策。而经验力场方法虽然计算效率较高，可以模拟大规模分子系统的运动，但是该方法依赖实验数据的特性决定了其难以达到量子力学精度。想要实现精度和效率的统一，我们必须从底层算法以及软硬件层面同时进行改进和优化，其中GPU提供的加速能力也是至关重要的一个环节。虽然目前的技术距离这个目标尚有差距，但假如我们真的能在量子力学精度下实现对介观乃至宏观体系的模拟，那许多困扰人类已久的未解之谜都将找到答案。

预告

GPU的大规模并行计算能力在各种分子模拟方法中都存在着很大的发挥空间。事实上，无论是经典的分子模拟应用，还是在“AI+Science”理念下诞生的新型模拟软件，想要达到理想的性能都离不开GPU的加速。本系列文章的下篇将通过案例介绍GPU助力分子模拟的实现方式和亟待解决的难题，同时对未来方向做出展望。

参考资料

1.NVIDIA Expands HPC Footprint at SC17, as Study Shows GPU Acceleration Key for Scientific Computing. Available from: https://blogs.nvidia.com/blog/2017/11/13/hpc-apps-top500/.

2.Gartner, T.E. and A. Jayaraman, Modeling and Simulations of Polymers: A Roadmap. Macromolecules, 2019. 52(3): p. 755-786.

3.Keith, J.A., et al., Combining Machine Learning and Computational Chemistry for Predictive Insights Into Chemical Systems. Chem Rev, 2021. 121(16): p. 9816-9872.

4.LAMMPS Molecular Dynamics Simulator. Available from: https://www.lammps.org/.

5.Gromacs. Available from: https://www.gromacs.org/.

6.NAMD, Scalable Molecular Dynamics. Available from: https://www.ks.uiuc.edu/Research/namd/.

7.OpenMM, High performance, customizable molecular simulation.; Available from: https://openmm.org/.

8.Hoomdblue. Available from: http://glotzerlab.engin.umich.edu/hoomd-blue/.

9.Ross, I.G., Calculations of the energy levels of acetylene by the method of antisymmetric molecular orbitals, including σ-π interaction. Transactions of the Faraday Society, 1952. 48(0): p. 973-991.

10.Hartree, D.R. and W. Hartree, Self-consistent field, with exchange, for beryllium. 1935. 150(869): p. 9-33.

11.Kohn, W. and L.J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects. Physical Review, 1965. 140(4A): p. A1133-A1138.

12.Brooks, B.R., et al., CHARMM: the biomolecular simulation program. J Comput Chem, 2009. 30(10): p. 1545-614.

13.Brooks, B.R., et al., CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. 1983. 4.

14.Gunsteren, W.F.v., X. Daura, and A.E. Mark. GROMOS Force Field. 2002.

15.Jorgensen, W.L., J.D. Madura, and C.J. Swenson, Optimized intermolecular potential functions for liquid hydrocarbons. 1984. 106:22: p. Medium: X; Size: Pages: 6638-6646 2016-05-16.

16.Jorgensen, W.L., D.S. Maxwell, and J.J.J.o.t.A.C.S. Tirado-Rives, Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. 1996. 118: p. 11225-11236.

17.Weiner, S.J., et al., An all atom force field for simulations of proteins and nucleic acids. J Comput Chem, 1986. 7(2): p. 230-252.

18.Behler, J., Four Generations of High-Dimensional Neural Network Potentials. Chemical Reviews, 2021. 121(16): p. 10037-10072.

19.Unke, O.T., et al., Machine Learning Force Fields. Chemical Reviews, 2021. 121(16): p. 10142-10186.

20.Zhang, L., et al., Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics. Physical Review Letters, 2018. 120(14): p. 143001.

21.Jia, W., et al. Pushing the limit of molecular dynamics with ab initio accuracy to 100 million atoms with machine learning. in SC20: International conference for high performance computing, networking, storage and analysis. 2020. IEEE.