456

布朗運動、郎之萬方程式、與布朗動力學 (Brownian Motion, Langevin Equation, and Brownian Dynamics) 文/王子瑜、曹恒光 在西元 1905 年前後，愛因斯坦(Albert Einstein)除研究狹義相對論外，也鑽研布朗運動。在愛因斯坦發表狹義相對論的 百年紀念之際，本文首先概述布朗運動並簡述愛因斯坦和史摩勒丘司基(Smoluchowski)如何以機率觀念詮釋布朗粒子的平均 行為。但現代對布朗運動的理論描述常採用較易瞭解的朗之萬(Langevin)理論，其特點是我們可以將許多流體分子碰撞布朗 粒子的效應簡化成一隨機熱擾動力，然後追蹤單一布朗粒子的運動軌跡。這種方法解決了直接以牛頓動量守恆方程追蹤系統 中所有布朗粒子與流體分子軌跡(稱為分子動力學)的困境：懸殊的時間尺度差異(timescales separation)。這種簡化布朗粒子周 遭流體貢獻的電腦模擬計算，我們稱之為布朗動力學，它已被廣泛地應用於膠體科學與生物物理領域。本文將就朗之萬方程 及布朗動力學作一簡單的介紹。

一、布朗運動、機率、與郎之萬方程式

粒子隨著時間往各個方向運動而遠離原點，將粒子的 移動距離先取平方，再取平均，我們將發現位移平方

西 元 1827 年 ， 英 國 的 植 物 學 家 勞 伯 ‧ 布 朗

的平均與所經過的時間t成正比；在同一時間條件下，

(Robert Brown)，在顯微鏡下觀察到懸浮在水中的花粉

布朗粒子遠離原點的快慢則代表布朗擴散係數D。愛

粒子，會不停地進行連續但不規則的運動。這種類似

因斯坦可說是第一位以定量理論詮釋布朗運動，稍後

生命體的運動特徵引發科學家們研究微小粒子的運動

史摩勒丘司基也發表以機率平衡方程式描述布朗運

行為。經過許多的實驗與探討，科學家發現這現象應

動。在無外力作用下，一維空間中的布朗運動可寫成，

該是微小粒子受到週遭液體分子從四面八方的連續撞

∂P ∂2 P = D 2 。 解 算 上 式 可 得 到 Einstein∂t ∂x

擊，而產生連續但不規則地隨機移動，這種移動我們 稱之為布朗運動(Brownian motion)。布朗運動具有下

Smoluchowski Equation, 〈x2〉 = 2Dt。

列的特性：(1)粒子的運動永不停止；(2)溫度的改變會 影響粒子的運動；(3)粒子的運動沒有固定的軌跡，其

根據愛因斯坦的布朗運動理論，法國物理學家皮

運動軌跡呈鋸齒狀；(4)粒子的大小影響粒子的運動速

林(Jean Perrin)進行膠體粒子(colloidal particles)的重力

度；(5)粒子的成份或密度不會影響粒子的運動。

沉降與布朗擴散的平衡實驗。密度比水重的膠體粒子 會因重力沉降至容器底端。但布朗運動(擴散)會使膠

西元1906年，愛因斯坦在發表狹義相對論後，也

體粒子往上懸浮。兩者平衡的結果會產生粒子濃度分

發表了他以『機率』的觀念探討布朗運動的定量結果。

佈。膠體粒子的化學位能(chemical potential, µ)可寫成

值得一提的是，雖然他以光電效應獲得諾貝爾物理

µ = µ0 + mgx + kBT lnc(x)，其中m和c分別代表膠體粒子

獎，他曾因布朗運動理論而被提名。根據愛因斯坦的

的質量與濃度。當系統處於熱力學平衡時，我們可得

研究分析，粒子的運動雖然不規則，但是布朗運動在

到c = c0 exp(−mgx/kBT)，其中c0代表粒子在x =0的濃

長時間下的平均移動行為會呈現常態分佈，可視為布

度。同樣的結果也可由愛因斯坦或史摩勒丘司基的布

朗粒子的擴散行為。依據布朗粒子在時間t與位置x時

朗運動理論求得。在重力作用的狀況下，布朗粒子處

的機率P(x,t)，我們可得到兩個重要結論，分別是(1)

在 位 置 x 的 機 率 P(x,t) 須 遵 守 方 程

粒子的位移平均為零(即〈x〉 = 0)。由於位移的向量特

∂P ∂  g ∂P  ，其中ζ表示粒子在流體中運動 = − P+ D  ∂t ∂x  ζ ∂x 

性，重覆將布朗粒子從原點釋出的實驗，由於布朗運 動的等向性，我們不難瞭解平均位移為零的結果。(2)

所受的摩擦阻力。平衡時，該方程的解為

物理雙月刊（廿七卷三期）2005 年 6 月 456

P=P0exp(−gx/ζD)。兩相比較，我們可得到布朗擴散係

將方程式(1)乘上x，我們得到

 d  dx   dx  2  d2x m 2 x = m  x −   dt  dt  dt   dt  

數為 D = − k BT ，該式稱為Nernst- Einstein Equation。 mζ 皮林的實驗觀察並量測到膠體濃度的分佈，實驗結果

(2)

dx = −mζ x + R(t ) x dt

證明了愛因斯坦的理論並由此求得理想氣體常數與亞 佛加厥常數。皮林因相關的實驗工作獲頒一九二六年

由於熱擾動力R(t)與粒子所處位置x並無相關性，所以

諾貝爾物理獎。

R(t ) x = 0 。同時熱力學平衡時，系統中粒子的平  dx  m   2  dt 

均動能代表溫度，即 1

2

=

1 。將上述兩項 k BT 2

事實帶入方程式(2)可解得 x2 = 2

 k BT  1 t − (1 − e−ζt )  mζ  ζ 

(3)

其中ζ。當 t >> ζ−1 ，我們可得到 x 2 = 2 k BT t 。將這 mζ 結果代入Einstein-Smoluchowski Equation，〈x2〉 = 2Dt， 我們也可推得Nernst-Einstein Equation D = k BT 。對於 mζ

上圖為 Perrin 實驗的示意圖。在一杯溶液中散佈著許多的膠 體粒子(0.29µm)，原本膠體應受到重力的影響而全部沉到容 器底部。但受到水分子的碰撞，膠體粒子進行布朗運動，而 使粒子懸浮於水中而不至於全部都沉到底部，因而呈現粒子 濃度分佈。

一半徑為a的球形粒子以速度v在黏度為η的流體中運 動，史托克斯(Stokes)透過解算動量守恆方程得到流體 力 學 阻 力 為 mζv ， 其 中 ζ =

雖然愛因斯坦是第一個以定量理論描述布朗運動

6πηa 代表摩擦係數 m

(friction coefficient) 。 將 該 結 果 代 入 Nernst-Einstein

與擴散的關係，但他與史摩勒丘司基是以機率平衡觀

Equation可得到Stokes-Einstein equation， D = k BT 。 6πηa

念來描述許多布朗粒子的平均行為，而非一顆布朗粒 子的行為。在西元1908年，郎之萬(Paul Langevin)發表 了可描述單一布朗粒子運動軌跡的方程，我們現在稱

綜而言之，兩類方法可用來描述布朗粒子在外加

之為『朗之萬』方程式(Langevin Equation)。雖然郎之

力場下的隨機運動。第一種方法是以機率平衡方程

萬對於布朗運動分析推導的方法與愛因斯坦不同，但

Fokker-Planck Equation來描述粒子在時間t、位置x、速

其軌跡的平均會與愛因斯坦直接透過機率所得到的結

度v時的機率P(x,v,t)；第二種方法則是透過Langevin

2

2

果吻合。郎之萬是依據牛頓第二定律(md x/dt = ΣF)，

Equation來描述粒子隨著時間t改變的運動軌跡。這些

考慮一個布朗粒子在運動時，同時受到流體的阻力

研究方法除了被使用在瞭解布朗運動外，也被運用到

mζ(dx/dt)與流體分子因熱運動與其碰撞的熱擾動力

其它熱擾動扮演重要角色的研究領域，例如化學反應

R(t)。 『朗之萬』方程式的推導如下，

動 力 學 (chemical dynamics) 和 生 物 奈 米 科 技

d2x dx m 2 = −ζ + R(t ) dt dt

(bio-nanotechnology)。前者例如，解算由Fokker-Planck (1)

Equation衍生而來的Kramers Equation，可求取分子越 過能量障壁的反應速率常數。後者例如，以原子力顯

物理雙月刊（廿七卷三期）2005 年 6 月 457

微鏡探討配位體與受體的生物鍵作用強度

mi

(ligand-receptor interaction)。實驗發現分開它們所需的 拉力會隨分開速度的升高而增加。這項特性也可透過

d 2 xi dx = −mi ζ i i + ∑ Fi j + R i 2 dt dt j

(4)

這裡 ri 與 mi 分別為布朗粒子 i 的位置與質量， ζ i

解算Langevin Equation而瞭解。

則代表該粒子的摩擦係數(friction coefficient)。假設摩 擦係數與粒子的位置和速度無關，同時摩擦效應具等

二、布朗動力論 (Brownian Dynamics) 『郎之萬』方程式的特性是將流體小分子與布朗

向性，則 ζ i 是純量，可由 Stokes’ Law ( ζ i = 6π ai η )決 mi

粒子的熱力學作用力表述成隨機的熱擾動作用力，而

定；其中 ai 為粒子 i 的半徑， η 為溶液的黏度。Fij 是

非平衡(相對運動)的作用力則以流體力學作用力表

粒子 i 受到系統中其他粒子 j 所施予的作用力(例如帶

述。由於『郎之萬』方程式具有簡易瞭解的特性，現

電粒子之間的庫倫作用力)，Ri 為流體分子碰撞布朗粒

在布朗運動理論的推導常採用其方式；此外，因為是

子 i 所施予的隨機熱擾動作用力。此隨機作用力的平

從牛頓第二定律推衍而來，『郎之萬』方程式中可以

均作用力為零，〈Ri(t)〉=0，且其共分散(covariance)為

直接地引入其它作用於布朗粒子的外力場。由於『郎

R i (t ) R j (t ) = 2ζ i k B T δij δ(t )I ，這裡 I 為 3 × 3 的單

之萬』方程式只代表單一布朗粒子的行為，集體行為

位張量(unit tensor)，k B 是波茲曼常數，T 為絕對溫度，

須透過許多軌跡的平均。過去，由於平均過程涉及大 量運算，相較於代表集體平均行為的機率平衡方程，

δi j 是 Kronecker delta (δij=0, 若 i≠j；δij=1, 若 i=j)。當

『郎之萬』方程式並未具有顯著優勢。現今，由於電 腦運算速度的快速發展，應用牛頓第二定律的分子動

流體的黏滯阻力很大或是僅對長時間的結構動態有興

力學和應用『郎之萬』方程的布朗動力學的『電腦模

趣時，我們可以忽略在方程式(4)左邊的慣性項，方程

擬』已被廣泛地採用來研究熱擾動力非常重要的膠體

式 簡 化 成 “ 位 置 朗 之 萬 方 程 ”(Position Langevin

與生物系統。

Equation)。

 dx i D  = i  ∑ Fi j + R i  dt k BT  j 

對一系統進行電腦模擬時，直接且嚴謹的作法是 追蹤系統內所有分子的運動軌跡。不幸的是，在大部 份的情形下，該法相當耗費計算時間。退而求其次，

(5)

方程式(5)對時間 t 積分，Ermak 等人得到簡單的計算

將所有流體分子對粒子的兩個主要貢獻，流體力學阻

方程

力與熱擾動力直接放至牛頓第二定律方程中，可大量 減少所需模擬的粒子(即流體分子)，因而加快模擬的

x i (t + ∆t ) = xi (t ) +

速度。藉由這種方式，研究者能夠直接研究布朗粒子 間的相互作用，及如何受外加力場的影響；利用電腦

Di k BT

∑ F (t )∆t + Z ( ∆t ) ij

i

(6)

j

注意熱擾動所造成的隨機位移 Zi 的大小是時間間隔∆t

動畫顯示，我們可以觀察系統中每個布朗粒子的運動

的函數。若位移隨機分佈呈高斯分佈，則其分散為

軌跡及其動態過程。這種採用『郎之萬』方程式模擬

Zi Z i ( ∆t ) = 2 Di ∆tI 。一般而言，方程式(6)是布朗

粒子的軌跡行為稱為布朗動力學。從流體分子的觀點 而 言 ， Brownian Dynamics 算 是 一 種 coarse-grained

動力學模擬的主宰方程，可得到布朗粒子隨時間變化

model。

的軌跡。下圖是簡略的布朗動力論演算法(Brownian Dynamics algorithms)的流程圖。

在含許多布朗粒子的系統中，遵循動量守恆概念 的『郎之萬』方程式可直接推展為

物理雙月刊（廿七卷三期）2005 年 6 月 458

頓第二運動定律追蹤求得，mi d2xi/dt2 = ΣjFij(|xj-xi|)。 給定系統中 N 個布朗粒

在這模擬過程中，布朗粒子所受的熱擾動力與流體阻

子的位置 xi(t=0)與∆t

力會透過所受合力項而自然地呈現。原則上，只要經 過足夠長時間的計算和取樣，可獲得粒子動態或平衡 態的結果。然而，流體分子與布朗粒子碰撞的時間尺 度遠小於布朗粒子運動的特徵時間，所以分子動力學

計算布朗粒子 i 在時間 t

模擬必須採用較小的時間間距∆t(short-time steps)來處

時，所受到的作用力 Fi

理較快的運動(fast motion)。由於絕大部份的時間都花

與隨機位移 Zi。

費在計算流體分子與布朗粒子的碰撞，為探討較慢的 布朗粒子移動的變化(evolution of slower modes)需要 非常長的模擬時間。由於系統粒子數龐大再加上截然

將作用力 Fi、隨機位移

不同的兩種時間尺度(timescales separation)，通常分子

Zi、與時間 t 時粒子 i 的

動力學的電腦運算需要花費相當長的時間，模擬複雜

位置 xi(t)代入方程式(6)

的生物分子動輒花費數個月更是司空見慣的事情，因 此昂貴的模擬計算時間是此方法的缺點。 蒙地卡羅法是在二次大戰後所發展出來的模擬方

求取時間 t+∆t 時，粒子

法，基本上可視為求算高維度積分的一種數值方法。

i 的位置 xi(t+∆t)。

應用於研究熱力學平衡系統時，等價於求算統計力學 中的核心：配分函數(partition function)。一九五三年，

三、布朗動力學(Brownian Dynamics)與分子動力

Mertropolis 等人提出重點取樣(importance sampling)的

學 (Molecular Dynamics) 和 蒙 地 卡 羅 法

方法，簡化並加速了蒙地卡羅法的記算。簡而言之，

(Monte Carlo method)的比較

蒙地卡羅法中的粒子是在龐大的相空間(phase space)

由於電腦運算能力的突飛猛進，分子模擬已被廣

上移動，並非在真正的時空中動態行為，所以所得到

泛地應用於各式各樣的研究中，例如蛋白質結構的變

的結果通常僅能求取平衡態的物理性質。相對於分子

性、聚電解質溶液的行為、或是去氧核醣核酸在多價

動力學，由於蒙地卡羅法缺乏真實時間尺度(timescales)

鹽類溶液中的結構變化特性等。分子模擬常採用布朗

的限制，很容易採用 coarse-grained model，將系統中

動力學、分子動力學與 蒙地卡羅法等計算方法。雖然

的溶劑的貢獻簡化成熱運動。同時因為熱力平衡態與

這些方法的主宰方程式不同，且各有其優缺點，但對

動態路徑無關，蒙地卡羅法不考慮流體阻力的影響。

同一系統其計算結果應相同。以下我們將簡單地介紹

以前述例子而言，模擬電解質溶液的平衡態時，只需

分子動力學與蒙地卡羅法，並與布朗動力學比較。

要考慮鹽類離子的行為；水分子的貢獻則表現在離子 的熱擾動與靜電作用的介電係數。因此，蒙地卡羅法

分子動力學模擬通常採用全原子模型(all atom

可以顯著地減少所需模擬粒子的數目，進而減少電腦

model)或將數個原子視為一個粒子的統一原子模型

運算所花費的時間。

(united atom model)。假若我們要模擬探討電解質在水 中的溶解行為，即鹽類在水中擴散的過程，我們必須

三、結語

同時計算所有的鹽離子與水分子的運動行為。知道粒 子之間的作用力，粒子動態軌跡與時間的關係可由牛

布朗運動的理論探討導致兩類研究方法的產生，

物理雙月刊（廿七卷三期）2005 年 6 月 459

分別為機率平衡方程式與『朗之萬』方程式。前者描

動的時間尺度)，則鍵的斷裂主要是因為外力已克服最

述布朗粒子在相空間(phase space)中某位置(x,v,t)的機

大鍵強度，所以斷裂拉力值接近鍵最大強度。在有限

率，著名的機率平衡方程包括 Smoluchowski Equation,

但足夠長的實驗時間下，由於熱擾動總是扮演角色，

Fokker-Planck Equation, Kramers Equation。後者則是

所以所量測的斷裂拉力值總是小於鍵最大強度。同樣

採用類似牛頓運動方程的策略，直接描述布朗粒子的

的方法可以說明以原子力顯微鏡拉力分開雙股去氧核

運動軌跡。為克服熱擾動與粒子移動這兩者時間尺度

醣核酸的實驗和以原子力顯微鏡量測平滑表面摩擦係

的巨大差異，朗之萬將流體對布朗粒子的影響簡化成

數的特性。

兩部份：隨機熱擾動和流體黏滯力。

參考文獻 由於是以力(動量)守恆為基礎，『朗之萬』方程 1. Ando, T., Meguro, T. and Yamato, I., J. Comput. Chem. 1, 3 (2002). 2. Chang, R. and Yethiraj, A., J. Chem. Phys. 116, 12 (2002). 3. Chen, J. C. and Kim, A. S., Adv. Colloid Interface Sci. 112, 159 (2004). 4. Ebeling, W., Cond. Matt. Phys. 7, 3 (2004). 5. Ermak, D. L., J. Chem. Phys. 62, 15 (1975). 6. Ermak, D. L. and Yeh, Y., Chem. Phys. Lett. 24, 243 (1974). 7. Ermak, D. L. and McCammon, J. A., J. Chem. Phys. 69, 1352 ( 1978). 8. Evans, E., Annu. Rev. Biophys. Biomol. Struct. 30, 105 (2001). 9. Lee, S. and Karplus, M., J. Chem. Phys. 81, 12 (1984). 10. Lemons, D. S. and Gythiel, A., Am. J. Phys. 65, 11 (1997). 11. Muthukumar, M. and Liu, S., J. Chem. Phys. 116, 22 (2002). 12. Pugnaloni, L. A., Ettelaie, R. and Dickinson, E., Langmuir 19, 6 (2003). 13. Rossky, P. J., Doll, J. D. and Friedman, H. L., J. Chem. Phys. 69, 15 (1978). 14. Turq, P. and Harold, F. L. and Friedman, J. Chem. Phys. 66, 7 (1977).

式很容易瞭解並配合研究系統的特性增減新的作用 項。正如同分子動力學的電腦模擬，這種特質也使得 『朗之萬』方程式演化成『布朗動力學』的電腦模擬。 『布朗動力學』可算是一種介於『分子動力學』和『蒙 地卡羅法』的電腦模擬方法。它保有分子動力學動量 守恆的原則，但也包含蒙地卡羅法中粒子受熱擾動隨 機移動的特性。所以布朗動力學可以研究系統中粒子 的動態行為，但卻抽取出周遭流體的影響，顯著地減 少模擬所需的運算時間。 目前布朗動力論模擬，已廣泛地應用在熱擾動扮 演不可或缺角色的科學和工程領域中。以前述的生物 鍵強度量測為例，生物鍵如同一位能井，當以原子力 顯微鏡嘗試拉斷生物鍵時，其過程即類似一布朗粒子 在外力的協助下逃離位能井。 先考慮兩種極端的情形：(1)若不考慮外力的協 助，生物鍵也能因熱擾動而斷裂(布朗粒子依靠熱擾動 逃出位能井)。若鍵結夠強，這種事件的機率極低(或 過程很久)。(2)若不考慮熱擾動的貢獻，外力必須大於 鍵的最大強度才能斷裂該鍵。換言之，該非布朗粒子 必須克服位能曲面上最陡處的阻力才能逃脫。在真實 的原子力顯微鏡的拉力實驗中，生物鍵同時受到周遭

作者簡介

其它原子的熱擾動與隨時間上升的外力的協助。兩者 的交互作用，使得鍵斷裂瞬間的拉力(rupture force)值

王子瑜就讀於中央大學化學工程與材料工程所博士

會隨拉力上升的速度(loading rate)而變。當拉力上升的

班。

極端地緩慢時，鍵的斷裂(布朗粒子的逃脫)主要是因

曹恒光教授現任於中央大學化學工程與材料工程所，

為長時間的熱擾動，所以斷裂拉力值很小。反過來，

研究領域為膠體與界面科學。

當拉力上升的速度非常快時(誇張地說，甚至小於熱擾

e-mail: [email protected]

物理雙月刊（廿七卷三期）2005 年 6 月 460