分子动力学模拟是研究复杂体系宏观现象与其微观结构关系的一个重要工具。
分子动力学模拟与理论研究、实验研究共同成为科学研究的重要手段。不同于理论研究中试图构建少量关键变量间关系以理解体系性质,分子动力学模拟直接考虑大量微观自由度的演化,对当前十分关心的复杂分子系统(如复杂流体、生物分子系统等)的研究尤为重要。另一方面,实验研究通常受到条件限制,只能测量特定条件下各种因素共同影响下的系统的少量特征,而分子动力学模拟在条件设定、分离与微观测量等方面十分便利,是实验研究的重要补充。
组成宏观体系的原子受到原子间相互作用力而发生运动,其运动规律符合牛顿力学,即每个原子的运动加速度正比于其受到的合力。统计物理学表明宏观体系表现出的性质与行为是其组成粒子运动统计平均的结果。随着现代计算机技术的快速发展,可以从分子间相互作用模型出发,模拟大量原子、分子的微观牛顿运动轨道,统计得到整个体系的宏观行为。分子动力学模拟通常考虑很多(几千到几十万个)原子的体系,根据这些原子的位置和速度,计算每个原子所受到其他原子(及外界环境)对其的作用力,从而得到这些原子经过一个很小时间步长(如几飞秒,1飞秒=10-15秒)后的速度改变及新的位置。重复这个计算很多步(如几百万步),可以得到体系在一段时间内(如几十纳秒,1纳秒=10-9秒)的微观轨道。通常实验测量的体系宏观性质(如压强、比热、扩散等)可以由这个微观轨道基于统计物理原理做平均计算得出。
分子动力学模拟进一步应用的最主要限制是,描述原子分子间相互作用的模型精确程度与模拟能达到的时间空间尺度相互矛盾,较高精度的模型通常需要更大的计算量,从而能达到的时空尺度较小,而要求模拟更大系统到更长时间尺度则一般不得不使用简化模型。在超级计算机中,基于全原子模型的分子动力学模拟一般能达到纳米微秒尺度,相比复杂体系中关键现象(如生物大分子功能性构象改变)的时空尺度常常不足。而基于更加近似简化的粗粒化模型的分子模拟,其精确程度又常常不足。为此发展高效准确的分子模型、高效的计算模拟算法,乃至高效的计算机硬件,一般性扩展分子动力学模拟的应用能力,是这个领域的一个重要方面。而针对具体问题,构建较简化但足够正确的微观模型,通过分子动力学模拟以理解宏观实验现象的微观机制,则是复杂分子体系研究的一个主要内容。
分子动力学模拟主要包括以下步骤。
①分子建模。分子动力学模拟首先要建立研究的分子系统的组成粒子及其间相互作用模型。考虑系统所有的组成原子及原子间的成键、非成键相互作用,称为全原子模型。如将一些原子团简化为单个粒子,考虑这些粒子受其他粒子相互作用而运动的,则是粗粒化模型。不同层面的粗粒化模型可以使用简化系统的微观描述与运动,从而大大简化对其进行的分子动力学模拟。随着粗粒化程度的增加,系统在被简化的同时,更多的微观细节将被忽略。全原子模型及各种粗粒化模型的相互作用一般都包含许多经验性参数,用以近似描述粒子间相互作用势。更准确的模型是基于体系全部原子核与电子(或原子实与价电子)的第一性原理模型。全原子模型可以看成是这种第一性原理模型的粗粒化结果。第一性原理模型基于量子力学原理计算体系的能量,不含经验性参数,称为量子化学计算,被广泛用于小系统的化学结构计算中。基于第一性原理的分子动力学模拟则是用此量子力学模型计算系统的能量及其对每个原子实位置的梯度,从而得到其所受的力。而原子实位置、速度在此受力下依据牛顿力学做分子动力学演化。若相互作用模型考虑电子则用量子力学计算,原子位置演化同通常的分子动力学模拟。
②环境耦合——热源与压强源。基于牛顿运动方程的分子动力学模拟的是孤立体系,系统总粒子数
、系统体积
及系统的总能量
固定,对应统计力学的微正则系综。在许多实际情况中,受关注的系综主要包括与恒温热源耦合,从而温度
恒定的正则系综(NVT系综),同时与热源及压强源耦合的NPT系综,与由外界恒定化学势的粒子源耦合的巨正则系综等。这些体系的分子动力学模拟,必须引入粒子与热源、压强源等外界环境的耦合,即粒子的运动不仅受系统内其他粒子相互作用影响,也受环境施加的作用力控制。因为系统在微观层面与环境的耦合具体形式通常未知,分子动力学模拟在计入这些环境耦合时是通过建立模型方式进行的,以保证这些耦合在宏观热力学性质等方面合理。
③初始构型。分子动力学模拟从系统的一个微观起始构型(系统所有的构成粒子的初始位置)出发。对复杂分子系统,其构成的所有原子应该有合理的初始位置,使得系统的能量较低,这是进行分子模拟的基础。这主要要求原子间没有大的位置重叠,成键原子的键长、键角等实验数据或量子化学计算相符因素等。这样的初始位置对应的分子间相互作用能量较低。
在确定起始构型后要赋予构成分子的各个原子速度,这一速度是根据玻耳兹曼分布随机生成的,由于速度的分布符合玻耳兹曼统计,因此在这个阶段体系的温度是恒定的。另外,在随机生成各个原子的运动速度后必须进行调整,使得体系总体在各个方向上的动量之和为零,即保证体系没有平动位移。
④粒子受力计算。分子间相互作用势对原子位置的梯度给出原子所受的合力。粒子间相互作用势决定体系的时间演化轨迹及其宏观表现。对不同体系及不同的研究问题,分子建模给出不同的分子间相互作用势。对非成键原子,大多使用二体相互作用,即两个原子(粒子)间的相互作用仅与它们自身的相对位置有关,与其他粒子无关。这样整个系统的势能被写成两两原子间相互作用之和。这种两体势又根据其具体相互作用形式,分别有刚球势,伦纳德-琼斯势、莫尔斯势等双体势模型。对于金属、碳硅材料,则常用多体势模型。也可直接考虑电子采用第一性原理计算相互作用。相对于二体势模型,多体势的计算量通常较大。直接使用第一性原理计算相互作用则计算量非常庞大,只能用在很小的系统中。
⑤位置速度时间积分。分子动力学模拟用有限时间差分方法计算粒子位置速度在受力下的时间演化。它选取一个小的时间步长,忽略在此时间步内粒子受力的改变,假定粒子做恒加速运动。数学上对粒子运动轨道做泰勒展开,忽略高阶项。为保证此差分方法长时间的能量稳定性,一般使用隐藏式差分方法、蛙跳法或Verlet积分方法,计算粒子位置速度经过该时间步长后的更新。分子动力学的时间步长选取很重要:太长时间步长会造成差分结果不准确,此时粒子间可能过于靠近,受力过大,数据溢出;太短的时间步长会过多计算粒子间受力,浪费计算资源,降低分子模拟过程所能达到的时间尺度。一般选取的时间步长根据体系最短的特征时间尺度定,为其几十分之一即可。一般可以通过粗粒化系统、增加粒子质量、软化相互作用势,或限制屏蔽无关的快过程等方式使得体系的最短时间尺度较长,从而以较大的时间步长,使得分子动力学模拟能达到更长的时间尺度。
⑥平衡趋近。分子动力学模拟从系统的初始构象出发,经每个时间步计算体系所有粒子所受合力,积分更新粒子的位置与速度。重复这个步骤可得到体系的微观时间演化轨道。由于体系的初始构象一般有一定的任意性,需要一段时间的演化以达到平衡态,消去初始构象的影响。在此模拟过程中,通过分析体系的温度、压强、能量,以及一些重要序参量随时间的演化,判断是否体系已经达到平衡。
⑦抽样。当模拟的体系达到平衡时,可以保存体系的微观构象以便分析。为保证足够的统计精度,以及研究系统的动力学行为,通常运行模拟较长时间,以固定时间间隔,抽取足够多的构象样本加以保存。抽样间隔时间通常是每几十到几百个时间步长抽样一次。更小的抽样时间间隔通常不必要,因为这样抽取样本间的关联较大,并不能有效增加统计精度。而太长的抽样间隔则浪费了模拟得到的样本。
⑧结果分析。体系的平衡性质及动力学行为,可从这个时间序列样本中使用统计物理原理计算分析得出。系统的一些关键量,如势能、动能、压强或密度,以及系统的关键参量,诸如简单液态中的径向分布函数(RDF)、高分子的首末端距离、相变序参量等平衡性质可以直接通过样本平均计算得到;而系统的动力学特征,如扩散系数、时间关联函数等也可由样本相应算出。
对一些常用体系,已经有多个常用软件包可以直接使用。如用于生物大分子的全原子模拟的主要有Gromacs、Amber、Namd、Charmm等。这些软件包多是源代码公开的,包括有对这些常用体系的全原子力场。其中Gromacs速度最快,Namd在大规模并行计算环境下有优势,而Amber对模拟DNA、RNA体系是首选。另外,一个常用的分子模拟软件包是Lammps,它的主要特点是程序结构良好、容易修改、使用灵活,特别适用于材料体系及通常的粗粒化模型的模拟。