張麗娟

（晉中學院信息技術與工程學院，山西晉中030619）

（編輯 郭繼榮）

0 引言

分子動力學模擬方法（molecular dynamicssimulation，MDS）是一種計算機模擬實驗的方法，是當前人們研究凝聚態(tài)系統(tǒng)的重要工具，原子的運動軌跡和原子運動過程中的各種微觀細節(jié)都可以通過這種方法得到，它對理論計算和實驗起到了有力的補充作用.這種方法是按體系內(nèi)每個粒子所遵從的動力學規(guī)律即牛頓運動規(guī)律來建立粒子的運動方程組，通過解方程組可以得到系統(tǒng)中每個粒子在每個時刻的坐標和動量，也就是能確定粒子的軌跡，根據(jù)經(jīng)典力學和統(tǒng)計物理理論，已知粒子的坐標和動量我們可以用統(tǒng)計計算的方法得到系統(tǒng)的各種性質(zhì)，比如：系統(tǒng)的各種熱力學性質(zhì)、粒子的均方位移、關聯(lián)函數(shù)等.所以分子動力學模擬能夠讓我們看到系統(tǒng)在一段時間內(nèi)的發(fā)展過程[2].

1 模擬的基本原理

對我們要研究的系統(tǒng)首先要選取作用勢的模型，對于一個粒子間的相互作用勢能是球?qū)ΨQ的多粒子體系，它的哈密頓量可以寫為：

其中，rij是第i個粒子和第j個粒子之間的距離.微正則系綜受到的約束有四個：能量守恒、體積不變、粒子數(shù)不變、動量恒等于零.系統(tǒng)中粒子的運動方程組可以從哈密頓量推出，如下所示：

采用有限差分法將微分方程化為有限差分方程來求解該方程組，數(shù)值求解該方程組，我們需要將此方程組的求解變成求解以下方程組：

從以上方程組可以得出，t+h時刻的坐標位置可以從前面t和t-h時刻這兩步的空間坐標位置及t時刻的作用力計算得到.若令：

則式（3）可以寫成如下形式：

需要注意的是，在第n+1步算出的速度是前一時刻，即第n步的速度，因而動能的計算比勢能的計算落后一步.

2 模擬的步驟設計

根據(jù)上面的原理我們設計的微正則系綜（粒子數(shù)恒定、體積恒定、能量恒定）的分子動力學模擬步驟如下：

（4）計算第n步的速度

在上面的算法中，動能的計算比勢能的計算要落后一步，而且這種算法不是自啟動的，只有給出初始位置｝和空間位置｝才能求出微分方程組的解，

按照改進后的算法模擬的步驟如下：

（5）重復進行第（2）步

3 平衡化的方法

大多數(shù)情況下，精確的初始條件是很難得到的，那么我們可以先給出一個合理的初始條件，接下來在模擬的時候不斷調(diào)節(jié)系統(tǒng)的能量使它達到給定的值.其方法為：首先求解系統(tǒng)的運動方程組，從中解出若干步的結(jié)果，接下來根據(jù)解出的結(jié)果計算體系的能量，看能量是否守恒而且等于我們給定的值，如若不然，那么我們就需要通過調(diào)整速度來使能量守恒，也就是用一個標度因子乘以速度，標度因子一般取為[2]：

然后回到第一步，對下一時刻的運動方程求解.反復進行上面的過程，直到系統(tǒng)達到平衡.

采用對速度標度的辦法，可以使速度發(fā)生很大變化[2]，為了消除速度標度帶來的后果，模擬過程中必需要有足夠的時間讓體系再次建立平衡，在到達趨衡階段以后，必須檢驗粒子的速度分布是否符合麥克斯韋-玻爾茲曼分布[5].

4 計算結(jié)果

計算的時候取256個氬原子作為我們的模擬系統(tǒng)，粒子間的相互作用勢為雷納德-瓊斯勢[6]：

其中，-ε是兩分子處于穩(wěn)定平衡狀態(tài)時的勢能值，也即勢能的最小值，當r＝216σ 時，U（r ）＝-ε；σ 是兩分子間作用勢為零時它們的間距，即當r＝σ時，U（r）＝0.該體系的粒子限制在一個立方體的箱子中，為了使密度守恒我們采用周期性邊界條件，邊界上采用最小像力約定.我們采用自然單位制，長度和時間的標度單位分別為 σ 和 （m σ248ε ）12，對氬原子的時間單位為3×10-12s，這樣就使得運動方程為無量綱形式.我們研究相圖上的一點 （T*， ρ*）＝ （0.722，0.831 34），立方體箱子的尺寸為6.75.初始條件定為：系統(tǒng)中每一個粒子處于面心立方格子的格點上，而速度按相應溫度下的玻爾茲曼分布抽樣取值，為了比較勢能的不同截斷對模擬結(jié)果的影響，我們?nèi)c＝2.5和rc＝3.6，速度標度因子為：

反復上面的速度調(diào)節(jié)，直到系統(tǒng)能量達到給定值，下面是計算結(jié)果：

4.1 能量演化過程

以下是256個氬原子系統(tǒng)在T*＝0.722，ρ*＝0.83134，rc＝2.5時，動能，勢能和總能量的演化過程.

圖1 動能演化圖

圖2 勢能演化圖

圖3 總能量演化圖

從上面各圖可以看出，微正則系綜的分子動力學模擬中，總能量是恒定的，從總能量的演化過程圖可以明顯地看出模擬的平衡化過程，其平衡化過程是通過對粒子速度的調(diào)節(jié)跳躍式地達到的.動能和勢能不是守恒量，但幾百步之后它們能給出持續(xù)的確定的平均值.

4.2 軌跡分析

我們在計算中初始條件令各個原子位于一個面心立方格子的格點上，模擬中隨著時間的推移系統(tǒng)內(nèi)原子的分布逐漸趨于無序，當原子分布完全處于無序狀態(tài)時，也即體系達到穩(wěn)定狀態(tài)，這時我們就可以計算體系的各個物理量從而研究體系的各種性質(zhì).圖4是體系內(nèi)原子的分布隨時間的變化圖，從圖中可以看出，初始狀態(tài)體系的原子分布處于有序狀態(tài)，模擬10步之后體系的有序狀態(tài)逐漸被打破，至到最后體系內(nèi)原子的分布完全處于無序狀態(tài).

圖4 各個時期256個氬原子的位置

4.3 徑向分布函數(shù)

徑向分布函數(shù)是一個用來研究物質(zhì)的有序性的參數(shù)，指給定某個粒子的坐標，其他粒子在空間的分布幾率，是反映物體結(jié)構(gòu)的特征物理量[2].分子動力學計算徑向分布函數(shù)的方法為[7，8]：

其中，N為分子數(shù)目，T為計算的時間步數(shù)，δr為設定的距離差，ΔN為介于r→r+dr間的分子數(shù)目.圖5是溫度為 （T*， ρ*）＝ （0.722，0.83134 ）和 （T*， ρ*）＝ （2.53，0.636 ）時氬原子體系的徑向分布函數(shù)，從 （T*， ρ*）＝（0.722，0.831 34）圖中可以看出，當r＞3.5σ時徑向分布函數(shù)的值趨近于1，當r＜0.89σ時徑向分布函數(shù)的值為零.這說明在我們所計算的氬粒子系統(tǒng)中，粒子之間的最小距離為0.89σ，而且平均看來，當粒子與粒子的距離超過3.5σ時，體系中粒子的分布是均勻的.在r~1.1σ處徑向分布函數(shù)有最大值，這就說明體系中在距離粒子為1.1σ的地方附近其他粒子出現(xiàn)的幾率是最大的，此處的局域密度應該是最大的.曲線中其他峰值雖比第一個峰值要小，但是徑向分布函數(shù)在此處有極大值，說明此處粒子分布的局域密度要大于平均密度.徑向分布函數(shù)反映出體系中原子的聚集特性，可借此了解體系的結(jié)構(gòu)[2].

圖5 不同溫度下的徑向分布函數(shù)

從兩曲線的對比可以看出，（T*， ρ*）＝ （2.53，0.636 ）時的第一峰值要比 （T*， ρ*）＝ （0.722，0.831 34）時的低.這說明：體系的配位數(shù)隨著溫度的升高而降低，有序度降低，這是符合熱力學規(guī)律的.（T*， ρ*）＝（2.53，0.636 ）時的曲線比 （T*， ρ*）＝ （0.722，0.831 34）時的曲線有一向r減小方向的平移，原因是：溫度高的原子動能比較大能夠克服它們之間的排斥作用，原子之間可以靠得更近.兩條曲線的共同特點是：隨著r的增大，體系的各個峰值逐漸減小，說明氬粒子系統(tǒng)短程有序、長程無序的特征.

5 結(jié)論

分子動力學模擬方法的應用已經(jīng)取得了許多重要的成果，可以廣泛地用于物理、化學、生物、材料、醫(yī)學等各個領域.本文通過對256個氬粒子系統(tǒng)進行微正則系綜分子動力學模擬得出了系統(tǒng)的能量演化過程、粒子的軌跡、系統(tǒng)的徑向分布函數(shù)，對了解體系特征有重要意義.

［1］陳正隆，徐為人，湯立達.分子模擬的理論與實踐［M］.北京：化學工業(yè)出版社，2007.

［2］Frenkel&Smit.分子模擬—從算法到應用［M］.汪文川，譯.北京：化學工業(yè)出版社,2002.

［3］D.W.Heemann.理論物理學中的計算機模擬方法［M］.秦克誠，譯.北京：北京大學出版社，1996.

［4］馬文淦.計算物理學［M］.北京：科學出版社，2005.

［5］汪志誠.熱力學統(tǒng)計物理［M］.北京：高等教育出版社，2003.

［6］毋志民，何煥典，羅強，等.銀、鈷和鉑原子納米團簇熔化過程的分子動力學模擬［J］.原子與分子物理學報，2004（51）：213~215.

［7］Bixon M，Jortner J.Energetic and thermodynamics sizeeffectsin molecular clusters［J］.J.Chem.Phys.,1989，91（3）：1631~1642.

［8］Schmidt M,Kusche R,et al.Negativeheat capacity for aclusterof 147 sodiumatoms［J］.Phys.Rev.Lett.,2001，86（7）：1191~1194.