布朗运动-郎之万方程 详细»
[广告]
物理雙月刊（廿七卷三期） 2005 年6 月 456 布朗運動、 郎之萬方程式、
與布朗動力學 (Brownian Motion, Langevin Equation, and Brownian Dynamics) 文/王子瑜、 曹恒光
在西元 1 905 年前後， 愛因斯坦(Albert Einstein)除研究狹義相對論外， 也鑽研布朗運動。 在愛因斯坦發表狹義相對論的百年紀念之際，
本文首先概述布朗運動並簡述愛因斯坦和史摩勒丘司基(Smoluchowski)如何以機率觀念詮釋布朗粒子的平均行為。
但現代對布朗運動的理論描述常採用較易瞭解的朗之萬(Langevin)理論， 其特點是我們可以將許多流體分子碰撞布朗粒子的效應簡化成一隨機熱擾動力，
然後追蹤單一布朗粒子的運動軌跡。 這種方法解決了直接以牛頓動量守恆方程追蹤系統中所有布朗粒子與流體分子軌跡(稱為分子動力學)的困境：
懸殊的時間尺度差異(timescales separation)。 這種簡化布朗粒子周遭流體貢獻的電腦模擬計算， 我們稱之為布朗動力學，
它已被廣泛地應用於膠體科學與生物物理領域。 本文將就朗之萬方程及布朗動力學作一簡單的介紹。 一、 布朗運動、 機率、 與郎之萬方程式 西元1 827 年，
英國的植物學家勞伯‧ 布朗 (Robert Brown)， 在顯微鏡下觀察到懸浮在水中的花粉粒子， 會不停地進行連續但不規則的運動。
這種類似生命體的運動特徵引發科學家們研究微小粒子的運動行為。 經過許多的實驗與探討， 科學家發現這現象應該是微小粒子受到週遭液體分子從四面八方的連續撞擊，
而產生連續但不規則地隨機移動， 這種移動我們稱之為布朗運動(Brownian motion)。 布朗運動具有下列的特性： (1 )粒子的運動永不停止；
(2)溫度的改變會影響粒子的運動； (3)粒子的運動沒有固定的軌跡， 其運動軌跡呈鋸齒狀； (4)粒子的大小影響粒子的運動速度；
(5)粒子的成份或密度不會影響粒子的運動。 西元1 906年， 愛因斯坦在發表狹義相對論後， 也發表了他以『機率』 的觀念探討布朗運動的定量結果。值得一提的是，
雖然他以光電效應獲得諾貝爾物理獎， 他曾因布朗運動理論而被提名。 根據愛因斯坦的研究分析， 粒子的運動雖然不規則，
但是布朗運動在長時間下的平均移動行為會呈現常態分佈， 可視為布朗粒子的擴散行為。 依據布朗粒子在時間t與位置x時的機率P(x,t)， 我們可得到兩個重要結論，
分別是(1 ) 粒子的位移平均為零(即〈x〉 = 0)。 由於位移的向量特性， 重覆將布朗粒子從原點釋出的實驗， 由於布朗運動的等向性，
我們不難瞭解平均位移為零的結果。 (2) 粒子隨著時間往各個方向運動而遠離原點， 將粒子的移動距離先取平方， 再取平均，
我們將發現位移平方的平均與所經過的時間t成正比； 在同一時間條件下，布朗粒子遠離原點的快慢則代表布朗擴散係數D。 愛因斯坦可說是第一位以定量理論詮釋布朗運動，
稍後史摩勒丘司基也發表以機率平衡方程式描述布朗運動。 在無外力作用下， 一維空間中的布朗運動可寫成， 2 2 P P D t x ∂ ∂ = ∂ ∂ 。 解 算
上 式 可 得 到 Einstein- Smoluchowski Equation, 〈x 2 〉 = 2Dt。 根據愛因斯坦的布朗運動理論，
法國物理學家皮林(Jean Perrin)進行膠體粒子(colloidal particles)的重力沉降與布朗擴散的平衡實驗。
密度比水重的膠體粒子會因重力沉降至容器底端。 但布朗運動(擴散)會使膠體粒子往上懸浮。 兩者平衡的結果會產生粒子濃度分佈。 膠體粒子的化學位能(chemical
potential, µ)可寫成 µ = µ 0 + mgx + k B T lnc(x)， 其中m和c分別代表膠體粒子的質量與濃度。 當系統處於熱力學平衡時，
我們可得到c = c 0 exp(−mgx/k B T)， 其中c 0 代表粒子在x = 0的濃度。
同樣的結果也可由愛因斯坦或史摩勒丘司基的布朗運動理論求得。 在重力作用的狀況下， 布朗粒子處在 位 置 x 的 機 率 P(x,t) 須 遵 守 方 程 P g
P P D t x x ∂ ∂ ∂ = − + ∂ ∂ ζ ∂ ， 其中ζ表示粒子在流體中運動所 受 的 摩 擦 阻 力 。 平 衡 時
， 該 方 程 的 解 為物理雙月刊（廿七卷三期） 2005 年6 月 457 P=P 0 exp(−gx/ζD)。 兩相比較， 我們可得到布朗擴散係數為 B
k T D m = − ζ ， 該式稱為Nernst- Einstein Equation。皮林的實驗觀察並量測到膠體濃度的分佈，
實驗結果證明了愛因斯坦的理論並由此求得理想氣體常數與亞佛加厥常數。 皮林因相關的實驗工作獲頒一九二六年諾貝爾物理獎。 上圖為 Perrin 實驗的示意圖。
在一杯溶液中散佈著許多的膠體粒子(0. 29µm)， 原本膠體應受到重力的影響而全部沉到容器底部。 但受到水分子的碰撞， 膠體粒子進行布朗運動，
而使粒子懸浮於水中而不至於全部都沉到底部， 因而呈現粒子濃度分佈。 雖然愛因斯坦是第一個以定量理論描述布朗運動與擴散的關係，
但他與史摩勒丘司基是以機率平衡觀念來描述許多布朗粒子的平均行為， 而非一顆布朗粒子的行為。 在西元1 908年， 郎之萬(Paul
Langevin)發表了可描述單一布朗粒子運動軌跡的方程， 我們現在稱之為『朗之萬』 方程式(Langevin Equation)。
雖然郎之萬對於布朗運動分析推導的方法與愛因斯坦不同， 但其軌跡的平均會與愛因斯坦直接透過機率所得到的結果吻合。 郎之萬是依據牛頓第二定律(md 2 x/dt 2
= ΣF)，考慮一個布朗粒子在運動時， 同時受到流體的阻力 mζ(dx/dt) 與流體分子因熱運動與其碰撞的熱擾動力 R(t)。『朗之萬』 方程式的推導如下，
2 2 ( ) d x dx m R t dt dt = −ζ + (1 ) 將方程式(1 )乘上x， 我們得到 2 2 2 ( ) d x d dx dx m
x m x dt dt dt dt dx m x R t x dt = − = − ζ +
(2) 由於熱擾動力R(t)與粒子所處位置x並無相關性， 所以 ( ) 0 R t x = 。 同時熱力學平衡時， 系統中粒子的平均動能代表溫度， 即 2 1
1 2 2 B dx m k T dt = 。 將上述兩項事實帶入方程式(2)可解得 ( ) 2 1 2 1 t B k T x t e m
−ζ = − − ζ ζ (3) 其中ζ。 當 1 >> t ζ − ， 我們可得到 2 2 B k T x t m = ζ
。 將這結果代入Einstein-Smoluchowski Equation， 〈x 2 〉 = 2Dt，我們也可推得Nernst-Einstein
Equation B k T D m = ζ 。 對於一半徑為a的球形粒子以速度v在黏度為η的流體中運動，
史托克斯(Stokes)透過解算動量守恆方程得到流體力學阻力為mζv ， 其中 6 a m πη ζ = 代表摩擦係數 (friction
coefficient) 。 將該結果代入Nernst-Einstein Equation可得到Stokes-Einstein equation， 6 B k
T D a = πη 。 綜而言之， 兩類方法可用來描述布朗粒子在外加力場下的隨機運動。 第一種方法是以機率平衡方程 Fokker-Planck
Equation來描述粒子在時間t、 位置x、 速度v時的機率P(x,v,t) ； 第二種方法則是透過Langevin
Equation來描述粒子隨著時間t改變的運動軌跡。 這些研究方法除了被使用在瞭解布朗運動外， 也被運用到其它熱擾動扮演重要角色的研究領域， 例如化學反應動 力
學 (chemical dynamics) 和 生 物 奈 米 科 技 (bio-nanotechnology)。 前者例如， 解算由Fokker-Planck
Equation衍生而來的Kramers Equation， 可求取分子越過能量障壁的反應速率常數。 後者例如， 以原子力顯物理雙月刊（廿七卷三期） 2005
年6 月 458 微鏡 探 討 配 位 體 與 受 體 的 生 物 鍵 作 用 強 度 (ligand-receptor interaction)。
實驗發現分開它們所需的拉力會隨分開速度的升高而增加。 這項特性也可透過解算Langevin Equation而瞭解。 二、 布朗動力論 (Brownian
Dynamics) 『郎之萬』 方程式的特性是將流體小分子與布朗粒子的熱力學作用力表述成隨機的熱擾動作用力， 而非平衡(相對運動)
的作用力則以流體力學作用力表述。 由於『郎之萬』 方程式具有簡易瞭解的特性， 現在布朗運動理論的推導常採用其方式； 此外， 因為是從牛頓第二定律推衍而來，
『郎之萬』 方程式中可以直接地引入其它作用於布朗粒子的外力場。 由於『郎之萬』 方程式只代表單一布朗粒子的行為， 集體行為須透過許多軌跡的平均。 過去，
由於平均過程涉及大量運算， 相較於代表集體平均行為的機率平衡方程，『郎之萬』 方程式並未具有顯著優勢。 現今， 由於電腦運算速度的快速發展，
應用牛頓第二定律的分子動力學和應用『郎之萬』 方程的布朗動力學的『電腦模擬』 已被廣泛地採用來研究熱擾動力非常重要的膠體與生物系統。 對一系統進行電腦模擬時，
直接且嚴謹的作法是追蹤系統內所有分子的運動軌跡。 不幸的是， 在大部份的情形下， 該法相當耗費計算時間。 退而求其次，將所有流體分子對粒子的兩個主要貢獻，
流體力學阻力與熱擾動力直接放至牛頓第二定律方程中， 可大量減少所需模擬的粒子(即流體分子)， 因而加快模擬的速度。 藉由這種方式，
研究者能夠直接研究布朗粒子間的相互作用， 及如何受外加力場的影響； 利用電腦動畫顯示， 我們可以觀察系統中每個布朗粒子的運動軌跡及其動態過程。
這種採用『郎之萬』 方程式模擬粒子的軌跡行為稱為布朗動力學。 從流體分子的觀點而言， Brownian Dynamics 算是一種coarse-grained
model。 在含許多布朗粒子的系統中， 遵循動量守恆概念的『郎之萬』 方程式可直接推展為 2 2 i i i i i i j i j d d m m dt
dt = − ζ + + ∑ x x F R (4) 這裡 i r 與 i m 分別為布朗粒子i 的位置與質量， i ζ
則代表該粒子的摩擦係數(friction coefficient)。 假設摩擦係數與粒子的位置和速度無關， 同時摩擦效應具等向性， 則 i ζ 是純量， 可由
Stokes’ Law ( i ζ = 6 i i a m π η )決定； 其中 i a 為粒子i 的半徑， η 為溶液的黏度。 F ij 是粒子 i
受到系統中其他粒子 j 所施予的作用力(例如帶電粒子之間的庫倫作用力)， R i 為流體分子碰撞布朗粒子 i 所施予的隨機熱擾動作用力。
此隨機作用力的平均作用力為零， 〈R i (t)〉=0， 且其共分散(covariance)為 ( ) ( ) 2 ( ) i j i B ij t t k T
t = ζ δ δ R R I ， 這裡I 為3 3 × 的單位張量(unit tensor)， B k 是波茲曼常數， T為絕對溫度， i j δ
是Kronecker delta (δ ij =0, 若 i≠j； δ ij =1 , 若i=j)。 當流體的黏滯阻力很大或是僅對長時間的結構動態有興趣時，
我們可以忽略在方程式(4)左邊的慣性項， 方程式 簡化成 “ 位置朗之 萬 方程”(Position Langevin Equation)。 i i i j i
j B d D dt k T = + ∑ x F R (5) 方程式(5)對時間 t 積分， Ermak 等人得到簡單的計算方程 ( )
( ) ( ) ( ) + ∆ = + ∆ + ∆ ∑ x x F Z i i i ij i j B D t t t t t t k T (6)
注意熱擾動所造成的隨機位移Z i 的大小是時間間隔∆t 的函數。 若位移隨機分佈呈高斯分佈， 則其分散為 ( ) 2 i i i t D t ∆ = ∆ Z Z
I 。 一般而言， 方程式(6)是布朗動力學模擬的主宰方程， 可得到布朗粒子隨時間變化的軌跡。 下圖是簡略的布朗動力論演算法(Brownian Dynamics
algorithms)的流程圖。 物理雙月刊（廿七卷三期） 2005 年6 月 459 三、 布朗動力學(Brownian Dynamics)與分子動力學
(Molecular Dynamics) 和 蒙 地 卡 羅 法 (Monte Carlo method)的比較 由於電腦運算能力的突飛猛進，
分子模擬已被廣泛地應用於各式各樣的研究中， 例如蛋白質結構的變性、 聚電解質溶液的行為、 或是去氧核醣核酸在多價鹽類溶液中的結構變化特性等。
分子模擬常採用布朗動力學、 分子動力學與 蒙地卡羅法等計算方法。 雖然這些方法的主宰方程式不同， 且各有其優缺點， 但對同一系統其計算結果應相同。
以下我們將簡單地介紹分子動力學與蒙地卡羅法， 並與布朗動力學比較。 分子動力學模擬通常採用全原子模型(all atom model)
或將數個原子視為一個粒子的統一原子模型 (united atom model)。 假若我們要模擬探討電解質在水中的溶解行為， 即鹽類在水中擴散的過程，
我們必須同時計算所有的鹽離子與水分子的運動行為。 知道粒子之間的作用力， 粒子動態軌跡與時間的關係可由牛頓第二運動定律追蹤求得， m i d 2 x i /dt
2 = Σ j F ij (| x j -x i | )。在這模擬過程中， 布朗粒子所受的熱擾動力與流體阻力會透過所受合力項而自然地呈現。 原則上，
只要經過足夠長時間的計算和取樣， 可獲得粒子動態或平衡態的結果。 然而， 流體分子與布朗粒子碰撞的時間尺度遠小於布朗粒子運動的特徵時間，
所以分子動力學模擬必須採用較小的時間間距∆t(short-time steps)來處理較快的運動(fast motion)。
由於絕大部份的時間都花費在計算流體分子與布朗粒子的碰撞， 為探討較慢的布朗粒子移動的變化(evolution of slower
modes)需要非常長的模擬時間。 由於系統粒子數龐大再加上截然不同的兩種時間尺度(timescales separation)，
通常分子動力學的電腦運算需要花費相當長的時間， 模擬複雜的生物分子動輒花費數個月更是司空見慣的事情， 因此昂貴的模擬計算時間是此方法的缺點。
蒙地卡羅法是在二次大戰後所發展出來的模擬方法， 基本上可視為求算高維度積分的一種數值方法。應用於研究熱力學平衡系統時， 等價於求算統計力學中的核心：
配分函數(partition function)。 一九五三年， Mertropolis 等人提出重點取樣(importance sampling)的方法，
簡化並加速了蒙地卡羅法的記算。 簡而言之，蒙地卡羅法中的粒子是在龐大的相空間(phase space) 上移動， 並非在真正的時空中動態行為，
所以所得到的結果通常僅能求取平衡態的物理性質。 相對於分子動力學， 由於蒙地卡羅法缺乏真實時間尺度(timescales) 的限制， 很容易採用
coarse-grained model， 將系統中的溶劑的貢獻簡化成熱運動。 同時因為熱力平衡態與動態路徑無關，
蒙地卡羅法不考慮流體阻力的影響。以前述例子而言， 模擬電解質溶液的平衡態時， 只需要考慮鹽類離子的行為； 水分子的貢獻則表現在離子的熱擾動與靜電作用的介電係數。
因此， 蒙地卡羅法可以顯著地減少所需模擬粒子的數目， 進而減少電腦運算所花費的時間。 三、 結語 布朗運動的理論探討導致兩類研究方法的產生，計算布朗粒子i
在時間 t 時， 所受到的作用力 F i 與隨機位移Z i 。 將作用力 F i 、 隨機位移 Z i 、 與時間 t 時粒子i 的位置x i
(t)代入方程式(6) 求取時間 t+∆t 時， 粒子 i 的位置x i (t+∆t)。 給定系統中 N個布朗粒子的位置x i (t=0)與∆t
物理雙月刊（廿七卷三期） 2005 年6 月 460 分別為機率平衡方程式與『朗之萬』 方程式。 前者描述布朗粒子在相空間(phase
space)中某位置(x,v,t)的機率， 著名的機率平衡方程包括Smoluchowski Equation, Fokker-Planck Equation,
Kramers Equation。 後者則是採用類似牛頓運動方程的策略， 直接描述布朗粒子的運動軌跡。 為克服熱擾動與粒子移動這兩者時間尺度的巨大差異，
朗之萬將流體對布朗粒子的影響簡化成兩部份： 隨機熱擾動和流體黏滯力。 由於是以力(動量)守恆為基礎， 『朗之萬』
方程式很容易瞭解並配合研究系統的特性增減新的作用項。 正如同分子動力學的電腦模擬， 這種特質也使得『朗之萬』 方程式演化成『布朗動力學』
的電腦模擬。『布朗動力學』 可算是一種介於『分子動力學』 和『蒙地卡羅法』 的電腦模擬方法。 它保有分子動力學動量守恆的原則，
但也包含蒙地卡羅法中粒子受熱擾動隨機移動的特性。 所以布朗動力學可以研究系統中粒子的動態行為， 但卻抽取出周遭流體的影響， 顯著地減少模擬所需的運算時間。
目前布朗動力論模擬， 已廣泛地應用在熱擾動扮演不可或缺角色的科學和工程領域中。 以前述的生物鍵強度量測為例， 生物鍵如同一位能井，
當以原子力顯微鏡嘗試拉斷生物鍵時， 其過程即類似一布朗粒子在外力的協助下逃離位能井。 先考慮兩種極端的情形： (1 ) 若不考慮外力的協助，
生物鍵也能因熱擾動而斷裂(布朗粒子依靠熱擾動逃出位能井)。 若鍵結夠強， 這種事件的機率極低(或過程很久)。 (2)若不考慮熱擾動的貢獻，
外力必須大於鍵的最大強度才能斷裂該鍵。 換言之， 該非布朗粒子必須克服位能曲面上最陡處的阻力才能逃脫。 在真實的原子力顯微鏡的拉力實驗中，
生物鍵同時受到周遭其它原子的熱擾動與隨時間上升的外力的協助。 兩者的交互作用， 使得鍵斷裂瞬間的拉力(rupture
force)值會隨拉力上升的速度(loading rate)而變。 當拉力上升的極端地緩慢時， 鍵的斷裂(布朗粒子的逃脫)主要是因為長時間的熱擾動，
所以斷裂拉力值很小。 反過來，當拉力上升的速度非常快時(誇張地說， 甚至小於熱擾動的時間尺度)， 則鍵的斷裂主要是因為外力已克服最大鍵強度，
所以斷裂拉力值接近鍵最大強度。 在有限但足夠長的實驗時間下， 由於熱擾動總是扮演角色，所以所量測的斷裂拉力值總是小於鍵最大強度。
同樣的方法可以說明以原子力顯微鏡拉力分開雙股去氧核醣核酸的實驗和以原子力顯微鏡量測平滑表面摩擦係數的特性。 參考文獻 1 . Ando, T. ,
Meguro, T. and Yamato, I. , J. Comput. Chem. 1 , 3 (2002). 2. Chang, R. and
Yethiraj , A. , J. Chem. Phys. 11 6, 1 2 (2002). 3. Chen, J. C. and Kim, A. S. ,
Adv. Colloid Interface Sci. 11 2, 1 59 (2004). 4. Ebeling, W. , Cond. Matt.
Phys. 7, 3 (2004). 5. Ermak, D. L. , J. Chem. Phys. 62, 1 5 (1 975). 6. Ermak,
D. L. and Yeh, Y. , Chem. Phys. Lett. 24, 243 (1 974). 7. Ermak, D. L. and
McCammon, J. A. , J. Chem. Phys. 69, 1 352 ( 1 978). 8. Evans, E. , Annu. Rev.
Biophys. Biomol. Struct. 30, 1 05 (2001 ). 9. Lee, S. and Karplus, M. , J. Chem.
Phys. 81 , 1 2 (1 984). 1 0. Lemons, D. S. and Gythiel, A. , Am. J. Phys. 65, 11
(1 997). 11 . Muthukumar, M. and Liu, S. , J. Chem. Phys. 11 6, 22 (2002). 1 2.
Pugnaloni, L. A. , Ettelaie, R. and Dickinson, E. , Langmuir 1 9, 6 (2003). 1 3.
Rossky, P. J. , Doll, J. D. and Friedman, H. L. , J. Chem. Phys. 69, 1 5 (1
978). 1 4. Turq, P. and Harold, F. L. and Friedman, J. Chem. Phys. 66, 7 (1
977). 作者簡介 王子瑜就讀於中央大學化學工程與材料工程所博士班。 曹恒光教授現任於中央大學化學工程與材料工程所，研究領域為膠體與界面科學。 e-mail:
hktsao@cc. ncu. edu. tw
No comments:
Post a Comment