分子动力学模拟的结果是确定的么？
注：这本来是我一门课结题报告，仓促完成，写的也非常跳跃，放在这里作为一个记录留存的必要。
1.牛顿力学的迷思
1814年，大约在牛顿去世后的一个世纪，拉普拉斯有一段讲话，他说到：
「如果一瞬间我们获得了一种智慧，这种智慧能理解所有的力（而正是在这些力作用下大自然得以表现的生动活泼），能理解构成它的物体的各自的状态——一种足够庞大的智慧，我们能够将这些数据提交给她来分析，她将以同样的公式接受宇宙中最伟大的天体和最轻的原子的运动；对她来说，没有什么是不确定的，未来，就像过去一样，将呈现在它的眼前。」[1]
这当然是充满著雄心壮志的一番言论，充满著决定论的语调，一切都是确定的，我们在牛顿力学的基础上，只需要一个巨大无比算力的机器，给任何输入的数据，就再也没有我们无法确定的结果了。大约又过去了一百年，汤姆森说到：
「似乎对我来说，如果问「我们是否理解物理中的特定现象」，这个问题等同于说「我们能否为它建立一个力学模型？」」[2]
这几乎就等同于现在大家对于力学的认识了，力学某种意义上主要的工作就是建模，我们考虑主要的要素，考虑边界条件，采用大量的假设，譬如在材料力学中有梁的假设，杆的假设，弹性力学中我们只考虑线弹性小变形，有没有假设的力学模型么？假设总是存在于力学模型建立的第一步。然而事实上，就算在阿波罗飞船早已经在大约半个世纪前载著人类登月的今天，我们依然无法预测的事情太多，譬如明天的股票行情，明天的天气。虽然有大量的经济学的、气象地理学的模型，但是这些模型往往都是不确定的，这无疑给拉普拉斯那充满雄心壮志的猜想带来了巨大的打击，我们真的无法制造一个模拟未来一切发生的容器么？
2.分子动力学模拟的非线性
2.1 分子动力学模拟的发展简史
分子动力学模拟的理论基础依然是古老的牛顿力学，自从牛顿1687年发表《自然哲学的数学原理》，经典力学一直紧紧地拥抱著牛顿三大定律。在1930年，蒙特卡罗提出了蒙特卡罗模拟，最简单也是最能理解的一个蒙特卡罗模拟便是在一个正方形上随机生成点，通过统计在一个给定圆半径的1/4圆内和圆外的概率分布，就可以用来计算圆周率。如图 1所示，当点的个数n=3000时，计算得到的圆周率为3.1133。在蒙特卡罗模拟的基础上，在1957年，Alder和Wainwright使用一台IBM 704晶体管计算机模拟了硬球的弹性碰撞[3]。时至今天，我们已经能够模拟大概1000万原子级别尺度的分子动力学模拟[4]，时间尺度可以达到微妙级[5]。
图 1 用来计算圆周率的蒙特卡罗模拟
2.2 分子动力学模拟的计算原理
人们对于大自然的认识程度，或者说物理学发展进步的衡量尺度可以是人们对于物理现象的观测尺度。这也是为什么理论物理学家在研究夸克等极小的粒子，而天体物理学家专注于观察星系的演变和宇宙的边缘。这样的尺度不仅仅存在于时间，也存在于空间当中。分子动力学模拟的计算元素是原子，与第一性原理考察电子场，求解薛定谔方程不同，在分子动力学模拟的世界里，原子就是组成这个世界的最小元素，不可再分。通过第一性原理计算或者实验，我们获得分子动力学模拟的力场，在基于牛顿经典力学，我们便可求解固定时间步原子与分子的物理运动，通过获得粒子的位移矢量和动量矢量，在统计力学的方法上求得所需要计算的物理量。简单来说，分子动力学是基于牛顿力学在热力学系统下研究原子与分子的物理运动的模拟方法。
2.3 真实系统、模型和模拟
「真实系统」更像是一个哲学的概念，但是就目前我们的观测手段而言，我们可以认为我们所拍到的氢原子电子结构图[6]是最逼近真实系统的存在。在分子动力学模拟中常见的力场是范德华势，譬如在研究纳米线在碳纳米管中的扭转问题，我们便可以先研究纳米线（以NbSe3为例）中Nb原子与Se原子之间的范德华力场，这就是一个力学模型。而整个扭转问题通过分子动力学的实现便是模拟。
严格来说，分子动力学模拟不是对「真实系统」的模拟，而是对于模型的模拟，而是基于力学模型的模拟，其模拟的对象并不是真实系统，而是力学模型。某种意义上，模型往往比真实系统简单，而模拟可以得到系统的各个不同的状态，因而比真实系统更复杂。用Haile, J. M. 的话来说——「We simulate molecular models of such substances」[8].
2.4 一个撞击的例子
分子动力学的经典模型便是硬球模型，如图 2所示，是硬球在连续边界和32个硬球组成的离散边界发生弹性碰撞的轨迹图，每次弹性碰撞径向速度改变完全相反的方向，切向速度不变。在原有的轨迹图的初始条件施加微小的扰动，连续边界下的轨迹是稳定的，而离散边界是不稳定的。经过400次撞击后，连续边界的轨迹是准周期的，而离散边界的轨迹是混沌的。
图 2 连续边界和32硬球离散边界下的硬球撞击：（A1）连续边界下的原轨迹图；
（A2）连续边界下受扰动的原轨迹图；（A3）400次撞击后连续边界的轨迹图；
（B1）32硬球离散边界下的原轨迹图；（B2）32硬球离散边界下受扰动的原轨迹图；（B3）400次撞击后离散边界的轨迹图；
2.5 另一个撞击的例子
如图 3所示，为同样演算法下，当不同球撞击时间相同或者差值在计算误差范围内，不同的撞击搜寻顺序将会导致不同的碰撞结果。结合前一个撞击粒子，我们这里提出一个问题：每一次分子动力学模拟的轨迹都不相同，统计出的物理量还是可靠的么?
图 3 同样演算法下出现的不同结果
2.6 物理量的统计
图 4 储氦单壁碳纳米管静水压作用下的分子动力学模拟[9]：
NVT系综用于将SWCNT和氦原子的温度分别设定为0.001K和100K。
模拟中的起始和终止压力分别为10MPa和100MPa。
图 5 双壁碳纳米管外管热泳驱动力的统计[10]
如图4所示，是作者之前所做的工作中的一个分子动力学模拟结构示意图，在对储氦单壁碳纳米管在静水压左右下进行模拟，在模拟过程中需要控制温度和压强，自然需要统计物理量。图5，是作者硕士课题的一部分工作，对双壁碳纳米管外管热驱动力进行统计，可以发现热驱动力会发生很大的抖动。在这里需要提到分子动力学模拟很重要的时间概念，分子动力学模拟所统计的物理量是在平衡态下一段时间内的平均，计算公式如图6所示。
图 6 分子动力学模拟中物理量的统计
3.可靠性与真实性的结果
综上所述，虽然我们发现在分子动力学模拟中存在许多非线性的现象，但是很多时候并没有影响最后物理量的统计，往往我们认为当系统演化至平衡态后，统计得到的物理量是可靠的。这其中的关联不得不引人深思，因此在这里放上埃格伯特·范坎彭的一段话作为结尾：
「经验告诉我们，尽管我们对大多数微观变数一无所知，但仍有可能发现宏观行为的规律，并用一般规律来表述。因此，这些微观变数的精确值似乎并不重要，因此我们不妨对它们进行平均。物理学家的任务是解释这个奇迹是如何发生的……在我看来，这是统计力学的关键问题。」[11]
参考文献：
- Binder K, Ceperley D M, Hansen J P, et al. Monte Carlo methods in statistical physics[M]. Springer Science & Business Media, 2012.
- Baumg?rtner A, Binder K, Hansen J P, et al. Applications of the Monte Carlo method in statistical physics[M]. Springer Science & Business Media, 2013.
- Alder, B. J.; Wainwright, T. E. (1959). "Studies in Molecular Dynamics. I. General Method". J. Chem. Phys. 31(2): 459.
- Lammps, see http://www.cs.sandia.gov/~sjplimp/lammps.html
- Reddy T, Shorthouse D, Parton DL, et al. Nothing to Sneeze At: A Dynamic and Integrative Computational Model of an Influenza A Virion. Structure. 2015;23(3):584-597.
- Stodolna A S, Rouzée A, Lépine F, et al. Hydrogen atoms under magnification: direct observation of the nodal structure of stark states[J]. Physical review letters, 2013, 110(21): 213001.
- Ying, et al. why the nanowires are twisted in single-walled carbon nanotubes?. In process.
- Haile J M, Johnston I, Mallinckrodt A J, et al. Molecular dynamics simulation: elementary methods[J]. Computers in Physics, 1993, 7(6): 625-625.
- Ying P, Zhao Y, Tan H. Study on collapse controlling of single-wall carbon nanotubes by helium storage[J]. Computational Materials Science, 2019, 164: 133-138.
- 应鹏华. 储氦单壁碳纳米管坍塌行为控制研究[D]. 2018.
- McCammon J A, Harvey S C. Dynamics of proteins and nucleic acids[M]. Cambridge University Press, 1988.
推荐阅读：