动力学蒙特卡洛方法及相关讨论

《动力学蒙特卡洛方法及相关讨论》由会员分享，可在线阅读，更多相关《动力学蒙特卡洛方法及相关讨论（15页珍藏版）》请在装配图网上搜索。

1、动力学蒙特卡洛方法(KMC)及相关讨论星期二，2010-05-1101:05satchel1979动态模拟在目前的计算科学中占据着非常重要的位置。随着计算能力和第一原理算法的发展，复杂的动态参数(扩散势垒、缺陷相互作用能等)均可利用第一原理计算得出。因此，部分复杂的体系动态变化，如表面形貌演化或辐射损伤中缺陷集团的聚合-分解演变等，已可以较为精确的予以研究。KMC动力学蒙特卡洛方法(kineticMonteCarlo)原理简单，适应性强，因此在很多情况下都是研究人员的首选。此外，KMC在复杂体系或复杂过程中的算法发展也非常活跃。本文试图介绍KMC方法的基础理论和若干进展。KMC方法基本原理在原

2、子模拟领域内，分子动力学(moleculardynamics,MD)具有突出的优势。它可以非常精确的描述体系演化的轨迹。一般情况下MD的时间步长在飞秒(LO-*s)量级，因此足以追踪原子振动的具体变化。但是这一优势同时限制了MD在大时间尺度模拟上的应用。现有的计算条件足以支持MD到10ns,运用特殊的算法可以达到10火s的尺度。即便如此，很多动态过程，如表面生长或材料老化等，时间跨度均在s以上，大大超出了MD的应用范围。有什么方法可以克服这种局限呢？当体系处于稳定状态时，我们可以将其描述为处于3N维势能函数面的一个局域极小值(阱底)处。有限温度下，虽然体系内的原子不停的进行热运动，但是绝大部分

3、时间内原子都是在势能阱底附近振动。偶然情况下体系会越过不同势阱间的势垒从而完成一次演化”，这类小概率事件才是决定体系演化的重点。因此，如果我们将关注点从原子”升格到体系”，同时将原子运动轨迹”粗化为体系组态跃迁”，那么模拟的时间跨度就将从原子振动的尺度提高到组态跃迁的尺度。这是因为这种处理方法挨弃了与体系穿越势垒无关的微小振动，而只着眼于体系的组态变化。因此，虽然不能描绘原子的运动轨迹，但是作为体系演化，其组态轨迹”仍然是正确的。此外，因为组态变化的时间间隔很长，体系完成的连续两次演化是独立的，无记忆的，所以这个过程是一种典型的马尔可夫过程(Markovprocess)，即体系从组态i到组态j

4、,一J这一过程只与其跃迁速率，有关。如果精确地知道屈久我们便可以构造一个随机过程，使得体系按照正确的轨迹演化。这里正确”的意思是某条给定演化轨迹出现的几率与MD模拟结果完全一致(假设我们进行了大量的MD模拟，每次模拟中每个原子的初始动量随机给定)。这种通过构造随机过程研究体系演化的方法即为动力学蒙特卡洛方法(kineticMonteCarlo,KMC)1指数分布与KMC的时间步长在KMC模拟中，构造呈指数分布的随机数是一个相当重要的步骤。这一节中我们对此进行讨论。因为体系在势能面上无记忆的随机行走，所以任意单位时间内，它找到跃迁途径的概率不变，设为因此在化f+.)区间内，体系不发生跃迁的概率为

5、理明(&)=1一福阈血+。(&户类似的，在区间内，体系不发生跃迁的概率为凡睫(2以)=(1一”twAt+(&)=1-2虹+0(St)2以此类推，当t=时，在f+丁；区间内，体系不发生跃迁的概率为KW=(1一虹味&+0(RJ)因此，当r趋于乂时，体系不发生跃迁的概率为玲町了)=1血一*(1-电。保=：(?xp(-Artrtr)(1)这一行为类似于原子核的衰变方程。从方程(1)我们可以得到单位时间内体系跃迁概率P(匚。从方程(1)的推导过程可以看出体系的跃迁概率是一个随时间积累的物理量，因此什一)对时间积分到某一时刻丫必然等于孔她伊)，也即旭=8(1-F耻瘠#)/乩因此我们立即可以得到1女小是体系

6、处于态时所有可能的跃迁途径的速率幻一，之和，即(3)对于每个具体的跃迁途径态u,上述讨论均成立。因此，我们可以定义单位时间内体系进行ij跃迁的概萃P汀为腊（4）单位时间内体系的跃迁概率呈指数分布这一事实说明KMC的时间步长扉也应是指数分布。因此我们需要产生一个指数分布的随机数序列。这一点可以非常容易的通过一个（0,1平均分布的随机数序列f转化得到:r=1-尺闾从而&=-ln(1-r)=ln(r)(5)最后一步是因为1和T的分布相同。Sf也可以通过上述步骤从方程（4）得到。计算跃迁速率过渡态理论（TST）k扣决定了KMC模拟的精度甚至准确性。为避开通过原子轨迹来确定化的做法（这样又回到了MD的情

7、况）,一般情况下采用过渡态理论（transitionstatetheory,TST）进行计算2。在TST中，体系的跃迁速率决定于体系在鞍点处的行为，而平衡态（势阱）处的状态对其影响可以忽略不计。如果大量的相同的体系组成正则系综，则在平衡状态下体系在单位时间内越过某个垂直于iij跃迁途径的纵截面的流量即为加。简单起见，假设有大量相同的一维双组态（势阱）体系，平衡状态下鞍点所在的假想面（对应于流量最小的纵截面）为F=，则TST给出该体系从组态A迁出到B的速率为5,6(6)方程(6)中表示在组态A所属态空间里对正则系综的平均。表示只考虑体系从组态A迁出而不考虑迁入A的情况(后一种情况体系也对通过纵截

8、面的流量有贡献)。根据普遍公式设体系的哈密顿量H为职施+3),即可分解为动能和势能，同时设粒子坐标件q时体系处于组态A。则方程(6)可写为如-二辱尸Mdj囱十|心口时疗m十口上,.*应*)|-一=*(碧?)“俯顼)A上式中无限小量是为了将莎函数全部包含进去。最后一项对于&函数的系综平均可以直接通过MetropolisMonteCarlo方法计算出来：计算粒子落在勺q+3范围内的次数相对于Metropolis行走总次数的比例/s。方程最后等于ST濡(8)5：将上述讨论扩展到3维情况非常直接，这里只给出结果，详细讨论请参阅文献9(9)其中f(R)是纵截面方程，|侦/代表3维情况中粒子流动方向与截面

9、丁法向不平行对于计数的影响。简谐近似下的过渡态理论(hTST)虽然上一节已经给出了TST计算跃迁速率的方法，但是在具体工作中，允新更多地是利用简谐近似下的过渡态理论(harmonicTST,hTST)通过解析表达式给出。根据TST，跃迁速率届为3二(10)其中为在跃迁ij中体系在鞍点和态j处的自由能之差3玲=穹丁砰_旧_向=Aj-T项将上式代入方程(10),可以得到-(11)hTST认为体系在稳态附近的振动可以用谐振子表示，因此其配分函数是经典谐振子体系的配分函数。分别写出体系在态i和鞍点处的配分函数事和:根据Boltzmann公式,并将配分函数代入，则方程(11)得3IVII-?,虹=一ex

10、p(13)n噂is.I方程(13)在通常的文献上经常可以见到。声子谱可以通过Hessian矩阵对角化或者密度泛函微扰法(DFPT)求出，而就是j的势垒，可以通过NEB或者drag方法求出。因此，方程(13)保证了可以通过原子模拟(MD或者DFT方法)解析地求出辰。事实上这个方程有两点需要注意。首先虽然方程(10)中出现了普朗克常数但是在最终结果中方被抵消了。这是因为TST本质上是一个经典理论，所以充分考虑了统计效应后/r不会出现1。其次，方程(13)表明对于每一个跃迁过程，鞍点处的声子谱应该单独计算。这样会大大增加计算量，因此在绝大部分计算中均设前置因子为常数，不随跃迁过程而变化。具体数值取决

11、于体系，对于金属而言，一般取10:-Hz。KMC几种不同的实现算法点阵映射到目前为止，进行KMC模拟的所有理论基础均已具备。但是前面所进行的讨论并没有联系到具体的模型。KMC在固体物理中的应用往往利用点阵映射将原子与格点联系起来。从而将跃迁（事件）具象化为原子一格点关系的变化。比如空位（团）/吸附原子（岛）迁移等等。虽然与实际情况并不完全一致，但这样做在很多情况下可以简化建模的工作量，而且是非常合理的近似。很多情况下体系中的原子虽然对理想格点均有一定的偏离，但是并不太大（pQOlfm）,因此这种原子一点阵映射是有效的。这种做法的另一个好处是可以对跃迁进行局域化处理。每条跃迁途径只与其近邻的体系

12、环境有关，这样可以极大的减少跃迁途径的数目，从而简化计算1。需要指出的是，这种映射对于KMC模拟并不是必须的。比如化学分子反应炉或者生物分子的生长等等，这些情况下根本不存在点阵。无拒绝方式KMC的实现方法有很多种，这些算法大致可以分为拒绝（rejection）和无拒绝（rejection-free）两种范畴。每种范畴之下还有不同的实现方式。本文只选择几种最为常用的方法加以介绍。I.直接法直接法（directmethod）是最常用的一种KMC算法，其效率非常高。每一步只需要产生两个在（，L之间平均分布的随机数1和其中被用来选定跃迁途径，2确定模拟的前进时间。设体系处于态f,将每条跃迁途径J想象成

13、长度与跃迁速率辰成正比的线段。将这些线段首尾相连。如果血烦落在线段中，这个线段所代表的跃迁途径上就被选中，体系移动到态九，同时体系时间根据方程（5）前进。总结其算法如下：1. 根据方程（4）计算体系处于态时的总跃迁速率临=为柘；选择随机数r_,AIjk,W*打1左寻找途径J,满足=1；体系移动到态元，同时模拟时间前进况=一上瓦（乓；重复上述过程。需要指出的是，虽然一般步骤4中的由根据方程（5）生成，但是如果将其换为&=上并不会影响模拟结果。在文献5和6中均采用这种方式。II.第一反应法第一反应法（firstreactionmethod,FRM）在思路上比直接法更为自然。前面说过，对于处于稳态i

14、的体系而言，它可以有不同的跃迁途径J可以选择。每条途径均可以根据方程（4）给出一个指数分布的”发生时间此,也即从当前算起i第一次发生的时间。然后从如中选出最小值（最先发生的”第一反应”），体系跃迁到相应的组态Jmin，模拟时间相应地前进屏打m.。总结其算法如下：1. 设共有条反应途径，生成*个随机数口，门舟_Ilij/r根据公式勺瓦给出每条路径的预计发生时间；找出择H的最小值沁,体系移动到态mm，同时模拟时间前进重复上述过程。可以看出，这种算法的效率比直接法低下，因为每一步KMC模拟需要生成个随机数。通常情况下KMC模拟需要10；步来达到较好的统计性质，如果每一步都需要生成21个随机数，则利用

15、这种方法需要一个高质量的伪随机数发生器，这一点在比较大时尤为重要。III. 次级反应法次级反应法（nextreactionmethod,NRM）是FRM方法的一种衍生方法，其核心思想是假设体系的一次跃迁并不会导致处于新态的体系对于其他跃迁途径的舍弃（比如充满可以发生种化学反应的分子，第一种反应发生并不会造成别的反应物的变化），这样体系还可以选择者们中的次小值此扁从而跃迁到态j新d,模拟时间前进&圮e-褊心如果这次跃迁还可以满足上述假设条件，再重复上述过程。理想情况下，平均每一步KMC模拟只需要生成1个随机数。这无疑会大大提高效率以及时间跨度。但是实际上NRM的假设条件很难在体系每次跃迁之后都得

16、到满足，在固体物理的模拟中尤其如此，因此其应用范围集中于研究复杂化学环境下的反应过程。试探-接受/拒绝方式这一大类算法虽然在效率上不如直接法，但是它们所采用的试探-接受/拒绝在形式上更接近MetropolisMC方法，而且可以很方便的引入恒定步长，即快固定。因此有必要进行详细的介绍。IV. 选择直接法选择直接法在决定体系是否跃迁方面和MetropolisMC方法形式上非常相像，均是通过产生随机数和预定的阈值比较决定事件是否被采纳。具体算法如下:1.设共有条反应途径，选择反应速率最大值炉识设为。生成在Q-UI均匀分布的随机数r；2.如果jk,则体系跃迁至新态j,否则保持在态；4.模拟时间前进5.

17、重复上述过程。这种方法的长处在于每一步只需要生成一个随机数。但是缺点也很明显，对于反应速率相差太大，尤其是只有一个低势垒途径（与其他途径相比昼过大）的体系来讲，这种方法的效率会非常低下。某些情况下，这种低效率问题可以通过如下方法改进：将全部途径按照化j的大小分为几个亚组，每个亚组选定一个上限人璋但是这一步骤在整个KMC模拟过程中可能需要重复很多次，因此并不能完全解决问题。事实上低势垒在KMC中是个普遍的问题。这一点在后面还要简要提及。V. 恒定步长法与上述四种方法不同，恒定步长法（constanttimestepmethod,CTSM）中体系的前进时间是个给定的参数citedawnkaski。

18、在理想情况下，CTSM与直接法效率相同，每一步只需产生两个随机数。具体算法如下：1. 给定恒定时间步长时；将所有途径（共有丁个）设为长度恒为LM的线段，生成在一山L均匀分布的随机数门，选择途径2. J=Int（门M）+1；生成在一1一均匀分布的随机数如果2舟，则体系跃迁至新态j,否则保持在态i,模拟时间前进重复上述过程。实际模拟中，鬲需要满足（1）小于屏疏职（见”第一反应法。，以及对于峋最大的途径，接受率大致在0.5。其中第一个条件保证了所有的迁移途径发生概率都小于1,第二个条件则保证体系演化的效率不会过于低下。CTSM是非常行之有效的一类KMC算法，但是选择金时需要特别的注意以保证效率。鬲决

19、定于具体体系以及模拟温度。这在一定程度上增加了CTSM的实现及使用难度。低势垒问题前面已经指出，低势垒的途径需要特别注意。如果体系在演化过程中一直存在着势垒较其他途径低很多的一个或几个途径，会对模拟过程产生不利的影响。这个问题被称之为低势垒问题。低势垒途径对于KMC模拟最直接的影响就是大大缩短了模拟过程所涵盖的时间跨度。这一点可以从方程（5）中看出。更为深刻的影响在于，这些由低势垒的途径联系起来的组态会组成一个近似于封闭的族。体系会频繁的访问这些态，而其他的对于体系演化更为重要的高势垒途径被选择的概率非常低，这显然会降低KMC的模拟效率。例如，吸附原子在高指数金属表面扩散，其沿台阶的迁移所对应

20、的势垒要远低于与台阶分离的移动。这样,而不会选择离开台阶在平台上扩散。这显然不KMC模拟的绝大部分时间内吸附原子都在台阶处来回往复,是我们希望看到的情形。一种解决办法是人为地将这些低势垒加高以降低体系访问这些组态的几率，但是无法预测这种干扰是否会造成体系对于真实情况的严重偏离。另一种选择是利用NRM或者CTSM进行模拟，但是其效果如何尚待检测。如果考察体系的势能面，这类低势垒的途径一般处在一个”超势阱”之中。体系在这个超势阱中可以很快的达到热平衡，所需时间要短于从其中逸出的时间。如果可以明确的知道超势阱所包含的组态以及从超势阱逸出的所有途径，我们就可以按照Boltzmann分布合理的选择其中一

21、条途径，使得体系向前演化。但是如何确定哪些组态包含在超势阱之中以及体系是否已在其中达到热平衡本身就是两个难题。对于第一点，Mason提出可以利用Zobrist密钥法标定访问过于频繁的组态7;Novotny则提出通过建立及对角化一个描述体系在这些组态间演化的传递矩阵来解决第二点8。对这个问题的详细讨论已超出了本文的讨论范围，请参阅文献7以及8o实体动力学蒙特卡洛方法OKMC上述的KMC都假设任何时候原子均处于其理想点阵格子上。但是很多情况下这种点阵映射是无效的，比如间隙原子或者位错。这类结构缺陷的运动在材料的辐射损伤和老化过程中扮演着非常重要的角色。而且与单个原子或者空位的运动相比，这类缺陷的运

22、动时间跨度更长，也更为复杂，比如间隙原子团和空穴的湮没，间隙原子团的解构/融合，或者位错的攀移/交滑移等等。传统的KMC算法很难有效的处理这类问题，一方面是因为时间跨度太大，另一方面这类缺陷各自均可视为独立的实体(object),其运动更近似于系统激发，因此单个或几个原子运动的积累效果很多情况下并不能有效地反应这些实体的整体运动。实体动力学蒙特卡洛方法(ObjectkineticMonteCarlo,OKMC)就是为了处理这类问题而被提出的。OKMC在算法上与普通的KMC完全一样。需要注意的地方是在OKMC中并不存在原子点阵。所有的实体在一个真空的箱子中按照其物理实质离散化运动，比如位错环的最

23、小移动距离是其Burgers矢量大小，方向则为Burgers矢量方向；空位的移动距离为第一近邻或第二近邻的原子间距，等等。模拟过程中我们需要追踪该实体的形心，从而决定其位置、移动距离等等。此外，OKMC中对于跃迁速率k的确定也和普通的KMC有所区别。本文前面已经指出，化可以表达为虹=-AE时/&丁)的形式。普通的KMC假定MJ为常数，不同途径的化j由上E,J决定。但是在OKMC的模拟中，3E”的直接确定非常困难，因此一般的策略是对于特定的事件（包括实体自身的运动以及不同实体间的反应等），跃迁势垒保持恒定，而将前置因子视为实体规模（所包含的原子/空位数目）的函数，通过MD模拟得出，一般而言可以表

24、示为形如侦（m）=M血T）771的表达式，其中fj和X是拟合参量，m是实体规模。最后需要注意的是在OKMC的模型中，实体有空间范围，因此需要一个额外的参数Entity来表征其空间半径（假设为球形分布，否则Entity的数目多于一个）。在模拟不同实体间的反应时，需要特别考虑其形心的间距，如果小于”反应距离，即Ri+反应一定进行，否则认为两个实体互相独立。Domain利用OKMC研究了Fe-Cu合金的辐射损伤9，在模拟中考虑了间隙原子（空位）的聚合、间隙原子（空位）团的发射、间隙原子-空位湮没、空位团对杂质的捕获、表面对于空位（团）的捕获、甚至辐射轰击引起的间隙原子（空位）萌生、增殖等等事件。从中

25、可以看出，对于OKMC,一个棘手的问题是需要预先想到所有的事件。此外，OKMC所需要的所有参量基本上不可能通过原子模拟直接获得，人为的设定参数不可避免。这些参数会在多大程度上决定OKMC的准确程度无法预先得知。需要根据现有的实验数据进行修改、调试。这些困难都限制了OKMC的普及。但是如前所述，这种方法可以有效地进行大尺度的时间（天）和空间模拟（m以上），而且对于缺陷的描述更为直接和符合直观，因此在材料研究中同样占有重要的地位。KMC的若干进展等时蛙跳算法（-leapKMC）引入这类算法前，我们先简要介绍两个常用的离散分布：泊松分布(PoissonDistribution,PD)以及二项式分布(

26、BinomialDistribution,BD)泊松随机数P腿M定义为给定事件发生率届.，以及观测时间生数目，则恰好等于t的概率是一个泊松分布:t下事件发生的数目。如果用n代表给定的发(14)P(nkij,r)=n=-飞如也即如果产生一个泊松随机数序列胃仕泌,则这个序列符合泊松分布PD。需要指出，PgjG是无界的，范围是任意非负整数与其类似，二项式随机数顼Vp）定义为重复N次独立的成功率均为P的伯努利实验的成功数。如果给定成功数则倒-Vp）恰好等于7!的概率是一个二项式分布：厂I一I（15）为了和本文中的标号一致，我们将跃迁iJ的成功率表示为，将方程（15）重新写为M站，丁）=nl（JV-n）

27、!（峪t/A）（1-蛎3）与PD不同，BD中的口是有界的，为0到A之间的任意整数。可以看出，如果将这两种随机数理解为给定跃迁路径（发生率为依j）在一定的时间步长（t）内发生的次数，则可以立即运用于粒子数空间N内的KMC中，其时间范围可以得到很大提高。这就是等时蛙跳算法r-leapKMC10,11。r-leapKMC方法最早由Gillespie提出，通过PD方程（14），在给定时间步长t下决定每个跃迁途径J发生的次数nJ,然后将体系移到这些跃迁累计发生后产生的新态。因为每一步模拟体系不止发生一次跃迁，所以模拟的速度可以大大加快。我们以多种反应物在化学反应炉中的演化为例加以详细说明。设在炉内共有A

28、种分子$】S,在时刻f各自的个数为，则在粒子数空间网中M（r）（川T构成一个矢量X（f）,或称为-个组态L总共有种反应路径用。=1M）。对于给定的用，反应速率加区（坎是占据态X（f；的函数。此外，我们单独定义-个矢量与=X-X,其中X；由X通过反应坊而得，即X。X：。因此她的元素勺代表反应用所引起的*种分子的数目变化。由此建立算法如下：PD-leapKMC101.给定恒定时间步长T；2.对于每条反应途径月.按照方程(14)生成泊松随机数序列,按照模拟步数从序列中找出每种反应发生的次数3.AJX(f+时=X(f)+4.模拟时间前进；5.重复上述过程。按照J=1更新体系;Gillespie仔细考虑

29、了丁的选择条件，称为蛙跳条件(leapcondition)r=min(17)其中M外(X=Lh(X左息=1攻x)=勇伽(Xj稣(X),J=1,um=i如前所述，没有上限，因此即使F满足方程(17)，在模拟过程中也可能会出现某种分子总数为负数的情况，这显然不符合实际，也是PD-T-leapKMC的一个弱点。Tian和Burrage提出可以用二项式分布BD取代PD,因为刿-Vp有上限，所以可以有效的解决这个问题。此外，他们对于某种分子参与多种反应的情况也进行了考虑，从而提高了T-leapKMC的稳定性和普适性。其算法如下：BD-leapKMC11I(log2A)15:将所有作为树叶(不足2整数次籍

30、的叶子由0填补)，每两片叶子之和作为父节点，依次类推直至树根M一株二叉树构建完毕后，生成一个随机数门，由树根开始寻找s=ntot，若四不大于左子节点临,沿左分支向下寻找；否则设&=财卜队，沿右分支向下寻找，直至树叶j,体系按途径演化。Slepoy和Thompson等进一步提出分流-拒绝(composition-rejection,CR)方法以实现搜索用时与途径总数无关的算法16:(1)先找出mir.和*max,按照(min心)将I条途径分为。个组，=1啤项卡二：，(2)然后生成随机数门，按照上述二叉树寻找=所落入的组别G.，,(3)再生成两个随机数匕和广3,设=其中Ng.为该组中包含的途径数，

31、f=法皿心，如果f如，则选择途径L否则重复步骤(3),直至有一条途径被选中为止。可以看出，CR算法虽然搜索速度很快，但是每一步KMC需要产生至少4个随机数(F1用于确定前进时间)，因此需要高质量的随机数发生器。不过对于跃迁途径复杂的体系演化而言，CR的：)(1)效率无疑是很有吸引力的。1 A.F.Voter,itRadiationEffectsinSolids(Springer2006)p.1-24.3,107(1935).2 H.Eyring,J.Chem.Phys.(NICSerices,Vol.42,P.Kratzer,MultiscaleSimulationMethodinMolecu

32、larScienceForschungszentrum,Jlich2009)p.51-76.4E.J.Dawnkaski,D.SrivastavaandB.J.Gamson,J.Chem.Phys.102,9401(1995).5A.F.VoterandJ.D.Doll,J.Chem.Phys.80,5832(1984).6A.F.Voter,Phys.Rev.B34,6819(1986).7D.R.Mason,T.S.HudsonandA.P.Sutton,Comp.Phys.Comm.165,37(2005).8M.A.Novotny,Phys.Rev.Lett.74,1(1994);Er

33、ratum75,1424(1995).9 C.Domain,C.S.BecquartandL.Malerba,J.Nucl.Mater.335,121(2004).10 D.T.Gillespie,J.Chem.Phys.115,1716(2001).11 T.TianandK.Burrage,J.Chem.Phys.121,10356(2004).12 G.HenkelmanandH.Jonsson,J.Chem.Phys.115,9657(2001).13 G.HenkelmanandH.Jonsson,J.Chem.Phys.111,7010(1999).14 O.Trushin,A.Karim,A.KaraandT.S.Rahman,Phys.Rev.B72,115401(2005).15 M.A.GibsonandJ.Bruck,J.Phys.Chem.A104,1876(2000).16 A.Slepoy,A.P.ThompsonandS.J.Plimpton,J.Chem.Phys.128,205101(2008).资料来源量子化学网