齊 楠，王圣軍

(陜西師范大學(xué) 物理學(xué)與信息技術(shù)學(xué)院，陜西 西安 710119 )

計(jì)算機(jī)模擬現(xiàn)已成為實(shí)驗(yàn)方法和理論研究以外的一種有效的物理研究手段，而分子動(dòng)力學(xué)方法作為一種確定性模擬方法，其在計(jì)算機(jī)模擬中顯得尤為重要.在許多計(jì)算物理教材[1-3]中，對于模擬由大量粒子組成的系統(tǒng)，分子動(dòng)力學(xué)方法都是尤為重要的內(nèi)容.

分子動(dòng)力學(xué)方法是指利用每個(gè)分子的運(yùn)動(dòng)方程通過統(tǒng)計(jì)方法來計(jì)算系統(tǒng)的性質(zhì).其以粒子間的相互作用為出發(fā)點(diǎn)，在計(jì)算機(jī)上求解運(yùn)動(dòng)方程的數(shù)值解.在這一求解過程中，涉及的物理量包括能量E、力F、質(zhì)量m、溫度T、位置r、時(shí)間t、速度v等.然而，在系統(tǒng)中，這些物理量的大小在國際單位制下具有過于小的數(shù)量級，這不利于在計(jì)算機(jī)上進(jìn)行運(yùn)算.于是，在分子動(dòng)力學(xué)方法中，首先需要對物理量重新標(biāo)度.這樣，對于重新標(biāo)度后的物理量的大小，其數(shù)量級在1附近，并且能夠得到形式簡單的運(yùn)動(dòng)方程.

在分子動(dòng)力學(xué)模擬中，重新標(biāo)度的方法是理解分子動(dòng)力學(xué)方法的重要環(huán)節(jié).在重新標(biāo)度的過程中，需要考慮系統(tǒng)的運(yùn)動(dòng)特征，這具有一定的物理內(nèi)涵.然而其在多數(shù)教材中顯得過于簡略，例如沒有介紹時(shí)間單位的物理內(nèi)涵.另外在模擬過程中標(biāo)度因子是一個(gè)尤為關(guān)鍵的量.在微正則系綜中，給定初始位置和初始速度以后，在系統(tǒng)趨衡過程中需要調(diào)節(jié)速度，使系統(tǒng)的溫度達(dá)到給定系統(tǒng)的溫度.對速度的重新標(biāo)度通過標(biāo)度因子實(shí)現(xiàn).然而，與之相關(guān)的書籍[1-3]往往直接給出了標(biāo)度因子的數(shù)學(xué)表達(dá)式，沒有明確地指出如何合理地對其進(jìn)行推導(dǎo).這使得分子動(dòng)力學(xué)模擬的原理難以被學(xué)生充分理解.

接下來，本文將對標(biāo)度因子的表達(dá)式予以推導(dǎo)，并且對物理量重新標(biāo)度的方法給出合理的說明.這樣會(huì)得到重新標(biāo)度的運(yùn)動(dòng)方程、物理量和標(biāo)度因子，幫助初學(xué)者掌握完整的模擬原理，以便于在面對新的情景時(shí)推導(dǎo)相應(yīng)的公式開展計(jì)算機(jī)模擬.

1 重新標(biāo)度前的公式

在分子動(dòng)力學(xué)模擬過程中，首先要設(shè)定模擬所采用的模型.這里考察最基本的情況微正則系綜，并且采用單原子粒子系統(tǒng)，分子間的相互作用勢采用Lennard-Jones勢.系統(tǒng)的哈密頓量表示為

(1)

其中分子間相互作用勢u(rij)為

(2)

這里rij是粒子i與粒子j的距離，ε和σ是兩個(gè)參數(shù).

由于統(tǒng)計(jì)物理系統(tǒng)的特性，精確地控制初態(tài)是困難的.模擬過程中人為設(shè)置的系統(tǒng)初態(tài)不會(huì)具有所要求的能量，并且系統(tǒng)的狀態(tài)不是一個(gè)平衡態(tài).為了使該系統(tǒng)的能量達(dá)到一個(gè)確定值，在系統(tǒng)趨衡的過程中調(diào)節(jié)系統(tǒng)能量，使其逐步達(dá)到所要求的狀態(tài).熱力學(xué)平衡狀態(tài)下該系統(tǒng)具有給定的溫度T，根據(jù)能量均分定理有

(3)

這里kB是玻耳茲曼常量.

(4)

(5)

從而得到標(biāo)度因子為

(6)

為了進(jìn)行計(jì)算機(jī)模擬，通常對系統(tǒng)進(jìn)行無量綱化處理.教科書中的速度標(biāo)度因子都是在無量綱化的模型中給出的.通常它的形式與式(6)的形式并不相同，并且它依賴于無量綱化的處理方式.

2 物理量的重新標(biāo)度

粒子i的牛頓運(yùn)動(dòng)方程為

(7)

(8)

于是

(9)

其中rij=ri-rj.對于計(jì)算機(jī)模擬，需要用到式(6)、(7)、(9).

這里需要對物理量進(jìn)行重新標(biāo)度，進(jìn)而把分子動(dòng)力學(xué)模擬中用到的物理量變成無量綱的量，使運(yùn)動(dòng)方程具有簡潔的形式.

另外，如果在模擬中仍然采用國際單位制，那么粒子的質(zhì)量m、位置r和基本時(shí)間步長h的數(shù)量級均很小，可能會(huì)導(dǎo)致舍入誤差的產(chǎn)生.于是，對物理量進(jìn)行重新標(biāo)度的另一個(gè)原則是：避免過于小或過于大的數(shù)量級出現(xiàn)在計(jì)算機(jī)模擬中[4].

為了得到具有簡潔形式的運(yùn)動(dòng)方程，本文并不是直接對力學(xué)中常用的基本物理量進(jìn)行重新標(biāo)度，再導(dǎo)出其它物理量.在分子動(dòng)力學(xué)模擬中，選取一些反映相互作用和運(yùn)動(dòng)特征的物理量，首先對它們進(jìn)行恰當(dāng)?shù)闹匦聵?biāo)度.進(jìn)而把這些重新標(biāo)度的物理量作為新的基本物理量[5,6]，由新的基本物理量導(dǎo)出其它物理量、標(biāo)度因子和運(yùn)動(dòng)方程.

現(xiàn)在對物理量重新標(biāo)度，并以星號(hào)*表示重新標(biāo)度后物理量的數(shù)值.從后文可以看到，一方面物理量的數(shù)值的數(shù)量級不會(huì)遠(yuǎn)遠(yuǎn)偏離1，另一方面運(yùn)動(dòng)方程經(jīng)過了無量綱處理.

根據(jù)式(2),容易得到-ε是分子間相互作用能的最小值，這個(gè)最小值在r=21/6σ處.參數(shù)ε決定了兩個(gè)粒子處于平衡位置的相互作用能，參數(shù)σ決定了兩個(gè)粒子處于平衡位置的距離.規(guī)定σ為長度單位，ε為能量單位，得到變換關(guān)系:

l=l*σ

(10)

E=E*ε

(11)

考慮粒子i與粒子j間相互作用的時(shí)間.在其質(zhì)心系中，將兩粒子所構(gòu)成的系統(tǒng)的總能量E表示為E=aε，則a為一個(gè)無量綱的參量.由該系統(tǒng)能量守恒[7]得到

(12)

(13)

上式為系統(tǒng)運(yùn)動(dòng)時(shí)長t關(guān)于兩粒子間距離rij的微分表達(dá)式.容易得到兩粒子間由最小距離r1處首次運(yùn)動(dòng)到最大距離r2處經(jīng)歷的時(shí)間T[7]為

(14)

圖1 的數(shù)值模擬結(jié)果

粒子間相互作用的時(shí)間也可以從以下的角度考慮.當(dāng)一個(gè)粒子的振動(dòng)范圍很小時(shí)，可以對系統(tǒng)中其它粒子對該粒子的作用做線性近似.考慮粒子i在單個(gè)粒子的勢場中的運(yùn)動(dòng).粒子i的位置用ri表示，考察粒子i的運(yùn)動(dòng)周期.由式(8)，粒子i在勢場中受到的作用力f(ri)為

(15)

則r0=21/6σ處是粒子i振動(dòng)的平衡位置.將f(ri)在r0處做泰勒展開，保留線性項(xiàng)，得到

(16)

其中x=ri-r0.于是粒子i的牛頓運(yùn)動(dòng)方程為

(17)

(18)

這里時(shí)間單位中包含一個(gè)常數(shù)48.它的作用是把相互作用力代入牛頓第二定律的時(shí)候，得到的運(yùn)動(dòng)方程形式更加簡單.

(19)

于是m*=48，即重新標(biāo)度后質(zhì)量m具有m*=48的數(shù)值大小.

由式(7)、(9)、(10)、(18)，得到重新標(biāo)度后的牛頓運(yùn)動(dòng)方程，其表示為

(20)

(21)

這里h*是重新標(biāo)度后的基本時(shí)間步長.

(22)

(23)

由式(11)、(23)，重新標(biāo)度后平均動(dòng)能表示為

(24)

由式(6)、(22)、(24),重新標(biāo)度后標(biāo)度因子β表示為

(25)

這里得到的標(biāo)度因子與相關(guān)書籍中的結(jié)果是一致的[1-3].當(dāng)導(dǎo)出了所有相關(guān)的公式以后，為了簡便，忽略所有無量綱物理量的*號(hào).

3 雙原子分子與量子分子動(dòng)力學(xué)方法

原則上，由大量粒子組成的系統(tǒng)都可應(yīng)用分子動(dòng)力學(xué)方法模擬.這個(gè)系統(tǒng)包括的不只是點(diǎn)粒子系統(tǒng)，也包括具有內(nèi)部結(jié)構(gòu)的粒子組成的系統(tǒng).如果按照系統(tǒng)所遵循的運(yùn)動(dòng)規(guī)律不同，分子動(dòng)力學(xué)方法可分為經(jīng)典分子動(dòng)力學(xué)方法和量子分子動(dòng)力學(xué)方法.經(jīng)典分子動(dòng)力學(xué)方法中，系統(tǒng)中每個(gè)分子各自都服從經(jīng)典運(yùn)動(dòng)定律.有時(shí)必須考慮系統(tǒng)中的量子效應(yīng)，這需要用量子力學(xué)來處理.這種考慮微觀組態(tài)變化影響的分子動(dòng)力學(xué)模擬稱為量子分子動(dòng)力學(xué)模擬，也稱為從頭計(jì)算的分子動(dòng)力學(xué)模擬或第一性原理計(jì)算[8].量子分子動(dòng)力學(xué)方法處理的是多電子的薛定諤方程:

Hψ(r,R)=Eψ(r,R)

(26)

其中，r和R分別代表所有電子和原子核的坐標(biāo)[8].基于密度泛函理論的第一性原理的計(jì)算方法由Car與Parrinello于1985年提出并用于實(shí)際計(jì)算.在過去的幾十年中，分子動(dòng)力學(xué)模擬發(fā)展很快，并在許多領(lǐng)域得到了廣泛的應(yīng)用[8,9].

根據(jù)分子動(dòng)力學(xué)方法的基本步驟，系統(tǒng)中相互作用勢函數(shù)的確定尤為關(guān)鍵.在分子動(dòng)力學(xué)方法中，通常是基于經(jīng)驗(yàn)數(shù)據(jù)和經(jīng)驗(yàn)規(guī)律來選取勢函數(shù).除了常用的Lennard-Jones勢，還有Morse勢、Kihara勢、諧振子勢、偶極子勢等.其中，Morse勢極好地描述了雙原子分子的振動(dòng)，其數(shù)學(xué)形式如下[10]:

V(r)=D(e-2αx-2e-αx)

(27)

在薛定諤方程的數(shù)值求解中，通常以埃為單位來度量所有的長度，以電子伏為單位來度量一切能量[11].在這樣一種單位制下，方程中的物理量和常數(shù)具有較好的數(shù)量級.在所謂的量子分子動(dòng)力學(xué)方法中，也可以類似地對系統(tǒng)做無量綱化處理.只需要對長度和能量分別進(jìn)行重新標(biāo)度.這種無量綱化處理方式與前文所述是一致的.

因此，確定好勢函數(shù)和運(yùn)動(dòng)方程后，總是可以按照本文的無量綱化處理方式對物理量進(jìn)行重新標(biāo)度，使運(yùn)動(dòng)方程更加簡潔.值得注意的是，其中對基本物理量的選取和重新標(biāo)度依賴于具體的物理模型，這具有一定的物理內(nèi)涵.

4 結(jié)語

本文首先給出了趨衡過程中速度標(biāo)度因子在國際單位制下的數(shù)學(xué)表達(dá)式.然后，說明了如何把物理量重新標(biāo)度得到無量綱公式，使得在模擬中計(jì)算更加方便.本文通過物理量的重新標(biāo)度，給出了時(shí)間單位的物理意義，并且給出了教材中常用的標(biāo)度因子相應(yīng)數(shù)學(xué)表達(dá)式的推導(dǎo)過程.文中的方法可以幫助初學(xué)者掌握完整的模擬原理，以便于在面對新的情景時(shí)推導(dǎo)相應(yīng)的公式開展計(jì)算機(jī)模擬.