齊 楠,王圣軍
(陜西師范大學 物理學與信息技術(shù)學院,陜西 西安 710119 )
計算機模擬現(xiàn)已成為實驗方法和理論研究以外的一種有效的物理研究手段,而分子動力學方法作為一種確定性模擬方法,其在計算機模擬中顯得尤為重要.在許多計算物理教材[1-3]中,對于模擬由大量粒子組成的系統(tǒng),分子動力學方法都是尤為重要的內(nèi)容.
分子動力學方法是指利用每個分子的運動方程通過統(tǒng)計方法來計算系統(tǒng)的性質(zhì).其以粒子間的相互作用為出發(fā)點,在計算機上求解運動方程的數(shù)值解.在這一求解過程中,涉及的物理量包括能量E、力F、質(zhì)量m、溫度T、位置r、時間t、速度v等.然而,在系統(tǒng)中,這些物理量的大小在國際單位制下具有過于小的數(shù)量級,這不利于在計算機上進行運算.于是,在分子動力學方法中,首先需要對物理量重新標度.這樣,對于重新標度后的物理量的大小,其數(shù)量級在1附近,并且能夠得到形式簡單的運動方程.
在分子動力學模擬中,重新標度的方法是理解分子動力學方法的重要環(huán)節(jié).在重新標度的過程中,需要考慮系統(tǒng)的運動特征,這具有一定的物理內(nèi)涵.然而其在多數(shù)教材中顯得過于簡略,例如沒有介紹時間單位的物理內(nèi)涵.另外在模擬過程中標度因子是一個尤為關(guān)鍵的量.在微正則系綜中,給定初始位置和初始速度以后,在系統(tǒng)趨衡過程中需要調(diào)節(jié)速度,使系統(tǒng)的溫度達到給定系統(tǒng)的溫度.對速度的重新標度通過標度因子實現(xiàn).然而,與之相關(guān)的書籍[1-3]往往直接給出了標度因子的數(shù)學表達式,沒有明確地指出如何合理地對其進行推導(dǎo).這使得分子動力學模擬的原理難以被學生充分理解.
接下來,本文將對標度因子的表達式予以推導(dǎo),并且對物理量重新標度的方法給出合理的說明.這樣會得到重新標度的運動方程、物理量和標度因子,幫助初學者掌握完整的模擬原理,以便于在面對新的情景時推導(dǎo)相應(yīng)的公式開展計算機模擬.
在分子動力學模擬過程中,首先要設(shè)定模擬所采用的模型.這里考察最基本的情況微正則系綜,并且采用單原子粒子系統(tǒng),分子間的相互作用勢采用Lennard-Jones勢.系統(tǒng)的哈密頓量表示為
(1)
其中分子間相互作用勢u(rij)為
(2)
這里rij是粒子i與粒子j的距離,ε和σ是兩個參數(shù).
由于統(tǒng)計物理系統(tǒng)的特性,精確地控制初態(tài)是困難的.模擬過程中人為設(shè)置的系統(tǒng)初態(tài)不會具有所要求的能量,并且系統(tǒng)的狀態(tài)不是一個平衡態(tài).為了使該系統(tǒng)的能量達到一個確定值,在系統(tǒng)趨衡的過程中調(diào)節(jié)系統(tǒng)能量,使其逐步達到所要求的狀態(tài).熱力學平衡狀態(tài)下該系統(tǒng)具有給定的溫度T,根據(jù)能量均分定理有
(3)
這里kB是玻耳茲曼常量.
(4)
(5)
從而得到標度因子為
(6)
為了進行計算機模擬,通常對系統(tǒng)進行無量綱化處理.教科書中的速度標度因子都是在無量綱化的模型中給出的.通常它的形式與式(6)的形式并不相同,并且它依賴于無量綱化的處理方式.
粒子i的牛頓運動方程為
(7)
(8)
于是
(9)
其中rij=ri-rj.對于計算機模擬,需要用到式(6)、(7)、(9).
這里需要對物理量進行重新標度,進而把分子動力學模擬中用到的物理量變成無量綱的量,使運動方程具有簡潔的形式.
另外,如果在模擬中仍然采用國際單位制,那么粒子的質(zhì)量m、位置r和基本時間步長h的數(shù)量級均很小,可能會導(dǎo)致舍入誤差的產(chǎn)生.于是,對物理量進行重新標度的另一個原則是:避免過于小或過于大的數(shù)量級出現(xiàn)在計算機模擬中[4].
為了得到具有簡潔形式的運動方程,本文并不是直接對力學中常用的基本物理量進行重新標度,再導(dǎo)出其它物理量.在分子動力學模擬中,選取一些反映相互作用和運動特征的物理量,首先對它們進行恰當?shù)闹匦聵硕?進而把這些重新標度的物理量作為新的基本物理量[5,6],由新的基本物理量導(dǎo)出其它物理量、標度因子和運動方程.
現(xiàn)在對物理量重新標度,并以星號*表示重新標度后物理量的數(shù)值.從后文可以看到,一方面物理量的數(shù)值的數(shù)量級不會遠遠偏離1,另一方面運動方程經(jīng)過了無量綱處理.
根據(jù)式(2),容易得到-ε是分子間相互作用能的最小值,這個最小值在r=21/6σ處.參數(shù)ε決定了兩個粒子處于平衡位置的相互作用能,參數(shù)σ決定了兩個粒子處于平衡位置的距離.規(guī)定σ為長度單位,ε為能量單位,得到變換關(guān)系:
l=l*σ
(10)
E=E*ε
(11)
考慮粒子i與粒子j間相互作用的時間.在其質(zhì)心系中,將兩粒子所構(gòu)成的系統(tǒng)的總能量E表示為E=aε,則a為一個無量綱的參量.由該系統(tǒng)能量守恒[7]得到
(12)
(13)
上式為系統(tǒng)運動時長t關(guān)于兩粒子間距離rij的微分表達式.容易得到兩粒子間由最小距離r1處首次運動到最大距離r2處經(jīng)歷的時間T[7]為
(14)
圖1 的數(shù)值模擬結(jié)果
粒子間相互作用的時間也可以從以下的角度考慮.當一個粒子的振動范圍很小時,可以對系統(tǒng)中其它粒子對該粒子的作用做線性近似.考慮粒子i在單個粒子的勢場中的運動.粒子i的位置用ri表示,考察粒子i的運動周期.由式(8),粒子i在勢場中受到的作用力f(ri)為
(15)
則r0=21/6σ處是粒子i振動的平衡位置.將f(ri)在r0處做泰勒展開,保留線性項,得到
(16)
其中x=ri-r0.于是粒子i的牛頓運動方程為
(17)
(18)
這里時間單位中包含一個常數(shù)48.它的作用是把相互作用力代入牛頓第二定律的時候,得到的運動方程形式更加簡單.
(19)
于是m*=48,即重新標度后質(zhì)量m具有m*=48的數(shù)值大小.
由式(7)、(9)、(10)、(18),得到重新標度后的牛頓運動方程,其表示為
(20)
(21)
這里h*是重新標度后的基本時間步長.
(22)
(23)
由式(11)、(23),重新標度后平均動能表示為
(24)
由式(6)、(22)、(24),重新標度后標度因子β表示為
(25)
這里得到的標度因子與相關(guān)書籍中的結(jié)果是一致的[1-3].當導(dǎo)出了所有相關(guān)的公式以后,為了簡便,忽略所有無量綱物理量的*號.
原則上,由大量粒子組成的系統(tǒng)都可應(yīng)用分子動力學方法模擬.這個系統(tǒng)包括的不只是點粒子系統(tǒng),也包括具有內(nèi)部結(jié)構(gòu)的粒子組成的系統(tǒng).如果按照系統(tǒng)所遵循的運動規(guī)律不同,分子動力學方法可分為經(jīng)典分子動力學方法和量子分子動力學方法.經(jīng)典分子動力學方法中,系統(tǒng)中每個分子各自都服從經(jīng)典運動定律.有時必須考慮系統(tǒng)中的量子效應(yīng),這需要用量子力學來處理.這種考慮微觀組態(tài)變化影響的分子動力學模擬稱為量子分子動力學模擬,也稱為從頭計算的分子動力學模擬或第一性原理計算[8].量子分子動力學方法處理的是多電子的薛定諤方程:
Hψ(r,R)=Eψ(r,R)
(26)
其中,r和R分別代表所有電子和原子核的坐標[8].基于密度泛函理論的第一性原理的計算方法由Car與Parrinello于1985年提出并用于實際計算.在過去的幾十年中,分子動力學模擬發(fā)展很快,并在許多領(lǐng)域得到了廣泛的應(yīng)用[8,9].
根據(jù)分子動力學方法的基本步驟,系統(tǒng)中相互作用勢函數(shù)的確定尤為關(guān)鍵.在分子動力學方法中,通常是基于經(jīng)驗數(shù)據(jù)和經(jīng)驗規(guī)律來選取勢函數(shù).除了常用的Lennard-Jones勢,還有Morse勢、Kihara勢、諧振子勢、偶極子勢等.其中,Morse勢極好地描述了雙原子分子的振動,其數(shù)學形式如下[10]:
V(r)=D(e-2αx-2e-αx)
(27)
在薛定諤方程的數(shù)值求解中,通常以埃為單位來度量所有的長度,以電子伏為單位來度量一切能量[11].在這樣一種單位制下,方程中的物理量和常數(shù)具有較好的數(shù)量級.在所謂的量子分子動力學方法中,也可以類似地對系統(tǒng)做無量綱化處理.只需要對長度和能量分別進行重新標度.這種無量綱化處理方式與前文所述是一致的.
因此,確定好勢函數(shù)和運動方程后,總是可以按照本文的無量綱化處理方式對物理量進行重新標度,使運動方程更加簡潔.值得注意的是,其中對基本物理量的選取和重新標度依賴于具體的物理模型,這具有一定的物理內(nèi)涵.
本文首先給出了趨衡過程中速度標度因子在國際單位制下的數(shù)學表達式.然后,說明了如何把物理量重新標度得到無量綱公式,使得在模擬中計算更加方便.本文通過物理量的重新標度,給出了時間單位的物理意義,并且給出了教材中常用的標度因子相應(yīng)數(shù)學表達式的推導(dǎo)過程.文中的方法可以幫助初學者掌握完整的模擬原理,以便于在面對新的情景時推導(dǎo)相應(yīng)的公式開展計算機模擬.