擴散模式的分類以及擴散係數的計算

看到一篇科技論文的報導納米粒子反常受限擴散研究獲得進展, 是有關擴散的, 就順便整理一下與擴散有關的理論知識, 這些知識在利用分子動力學模擬研究擴散現象時可能會用得上.

擴散模式的不同分類

使用MD計算粒子的擴散係數時, 一般都是基於粒子進行布朗運動, 遵循簡單的擴散模型. 實際上, 根據外界限制條件的不同, 粒子的擴散有多種模式, 不同擴散模式下粒子的均方位移MSD(mean square displacement)與時間t的關係很不一樣. 根據MSD的表現, 我們可以推測出擴散的不同機理. 研究粒子, 如量子點在細胞中的運動時, 人們總結了粒子幾種可能的擴散模式. 下面是論文Confined Lateral Diffusion Of Membrane Receptors As Studied By Single Particle Tracking (Nanovid Microscopy). Effects Of Calcium-induced Differentiation In Cultured Epithelial Cells中給出的總結.

靜態模式 Stationary mode

粒子基本不運動, 擴散係數接近零.

簡單(或布朗)擴散模式 Simple diffusion mode

粒子做簡單的布朗擴散運動

$\text{MSD}(\D t)=nD \D t$, $n$ 為粒子運動空間的維數

定向擴散(或傳輸)模式 Directed diffusion mode (transport mode)

粒子在做隨機擴散的同時沿某一方向以恆定的速度漂移

$\text{MSD}(\D t)=nD \D t+v^2 \D t^2$

限制擴散模式 Restricted diffusion mode

粒子在做布朗擴散時, 其運動空間被限制在某一範圍內, $0 \le x \le L_x$, $0 \le y \le L_y$, $0 \le z \le Lz$. 粒子的運動等價於在無限深方勢阱中的布朗擴散. 細分起來又可分為兩類: 受限模式confinement model和系索模式tethering model. 但僅根據軌跡無法區分這兩種模式.

$\alg
\left< x^2 \right>(t) &={Lx^2 \over 6}-{16 L_x^2 \over \p^4} \Sum{n=1(odd)}^\infty {1 \over n^4} \exp\left[ -{1\over2} ({n\p\sx \over L_x})^2 t \right], &\s_x^2 &=2 D_x\
\left< y^2 \right>(t) &={L_y^2 \over 6}-{16 L_y^2 \over \p^4} \Sum{n=1(odd)}^\infty {1 \over n^4} \exp\left[ -{1\over2} ({n\p\sy \over L_y})^2 t \right], &\s_y^2&=2 D_y \
\left< z^2 \right>(t) &={L_z^2 \over 6}-{16 L_z^2 \over \p^4} \Sum{n=1(odd)}^\infty {1 \over n^4} \exp\left[ -{1\over2} ({n\p\s_z \over L_z})^2 t \right], &\s_z^2&=2 D_z \
6D &=2D_x+2D_y+2D_z
\ealg$

受障模式 Obstacle-impeded diffusion mode

粒子進行自由擴散時受障礙物限制, 障礙物可以是固定的, 或者有一定的移動能力. 在這種情況下, 長程擴散會減弱但仍大於零. 這種模式比較難與限制模式區分開來. 對各向同性的碰撞機率,

$\text{MSD}(\D t)=4 D\D t+A+B\ln (C\D t)$

簡單擴散模式對應的擴散係數

大多數情況下MD研究的是簡單擴散模式, 在此假定下求擴散係數的步驟是先算MSD, 然後擬合MSD線性部分的斜率, 再除以6(二維體系除以4)就是擴散係數. 但在擬合的時候, 有一個選取擬合數據範圍的問題, 也就是確定用哪段時間範圍內的MS進行擬合. 很顯然, 使用不同時間段的數據, 擬合結果會有所不同. 擬合時不能使用MSD數據的起始部分, 因為起始部分屬於擴散弛豫過程, 除非專門研究這個過程, 擬合的初始時間要大於擴散弛豫時間; MSD數據的最後部分也不能用, 因為計算MSD時, 關聯時間越大, 數據點越少, MSD的誤差越大, 所以擬合時只能取數據中間接近線性的一部分.

如果是普通擴散過程, 在開始的一小段時間內MSD是關聯時間的二次函數, 代表無障礙的定向擴散. 隨著關聯時間的增加, MSD會很快過渡到一次函數階段, 代表正常擴散. 這個一次函數區域一般來說就是計算擴散係數的最佳區域.

由於關聯時間越大, 漲落越大, 所得MSD的誤差也越大, 所以擬合的最大關聯時間也不是越大越好. 此外, MSD曲線的光滑程度與模擬時間成正比, 一般模擬時間會取最大關聯時間的10倍或更大. 也可以多次模擬求平均值.

具體怎麼決定擬合的關聯時間範圍呢? 一種方法是選取不同的最大關聯時間, 如1 ps, 5 ps, 10 ps, 50 ps, 100 ps等作出MSD隨最大關聯時間變化的曲線, 從而確定體系的特徵擴散弛豫時間. 更好的方法是計算動擴散係數(RDC, running diffusion constant), 根據它的變化來確定關聯時間的擬合範圍

$$\text{RDC}(t) = { \rmd{\text{MSD} }(t) \over 6 \rmd t} \approx { \text{MSD}(t) \over 6t}$$

動擴散係數曲線要趨近於水平線才算收斂. 對於普通擴散, 隨關聯時間的增加, 動擴散係數會從零增加到一個定值. 對於氣體, 增加過程一般是單調的; 對於稠密液體, 可能會有局部漲落. 當然，單次模擬數據會有較大漲落, 我們很難通過一次模擬結果就判斷出這個定值是多少, 所以需要多次模擬取平均.

具體做法就是, 使用 MSD/6t 對t作圖, 這樣數據成正比的線性部分應該是一條上下稍有波動的水平線, 很容易看出來. 我們可以簡單地將這段時間內擴散係數的平均值模擬的擴散係數, 也可以對此時間段內的數據進行擬合求出擴散係數. 這兩種方法得到的擴散係數值應該差不多, 但建議採用後一種方法, 雖然麻煩點, 但標準

這種方法的一種變形是做雙對數 log(MSD)-log(t) 圖, 選取其斜率儘可能接近1的一段求擴散係數. 這很容易理解

$$\ln \text{MSD}(t) = \ln t+\ln(6D)$$

缺點在於單憑肉眼很難判斷數據斜率為1的部分.

更複雜的處理方法, 可以採用隨機抽樣一致性算法(RANSAC)或者非線性擬合方法來自動確定最佳擬合值, 其大致思想是在指定的條件下, 尋找一條直線, 使得距離直線一定範圍內的點數最多. 但這些方法使用麻煩, 所以文獻中較少看到.

MD計算設置

系綜選用: 大多用NVT, 但NPT也可以, 不過要注意NPT時盒子標度導致的位置變化

MSD計算屬於相關函數計算, 需要大量的數據點, 因此, 需要長的模擬時間以及短的輸出間隔. 只有模擬時間要足夠長才能保證擴散係數收斂.

擴散係數

單位 [D] = nm2/ps = 10-6 m2/s = 10-2 cm2/s = 103 10-5 cm2/s = 103 10-9 m2/s

1 cm2/s = 100 nm2/ps

298.2 K時, 水的擴散係數 Dwat=(2.30 ± 1.5%) x 10-9 m2/s=(2.30 ± 1.5%) x 10-5 cm2/s.

其他一些分子的擴散係數, 可查閱CRC物理化學手冊.

gmx msd計算時可以設定擬合時間, 計算後會輸出相關信息以及擴散係數, 但如果你不能確定擬合時間範圍的話, 給出的數值最好不要直接使用.

 %# Diffusion constants fitted from time 50 to 100 ps  % D[         F] = 3.3547 (+/- 0.1417) (1e-5 cm^2/s)

附: 反常擴散sub-diffusive

理論上講, 當時間足夠長時, 粒子的運動可以近似看作隨機遊走, 那麼MSD與關聯時間成線性關係. 在三維情況下, $\text{MSD}(t) = 6Dt$. 對上式兩邊取對數作圖, 斜率為1. 如果計算MSD的時間不夠長, 你會發現MSD與時間呈非線性關係 $\text{MSD}=D_\a t^\a (a \lt 1)$. 這時對兩邊取對數作圖斜率不為1. 這種情況稱為sub-diffusive行為. 直觀的解釋是粒子的遊走可以近似看作是在其他粒子形成的」溶劑籠」裡活動, 粒子長時間運動的MSD可近似看作粒子在各個溶劑籠間的活動, 因此與時間呈線性關係, 而粒子短時間運動的MSD可看作是粒子還沒有走出第一個溶劑籠, 因此是呈非線性的sub-diffusive行為. 為了克服sub-diffusive行為對計算擴散係數帶來的影響, 通常是計算足夠長時間的MSD, 將MSD/t或d(MSD)/dt對時間作圖, 考察MSD斜率隨時間的變化, 當斜率達到穩定值時就可以用來計算擴散係數了.

相關文獻

Diffusion And Subdiffusion Of Interacting Particles On Comblike Structures

O. Bénichou, P. Illien, G. Oshanin, A. Sarracino, R. Voituriez; Phys. Rev. Lett. 115(22):, 2015; 10.1103/physrevlett.115.220601

Discriminating Between Anomalous Diffusion And Transient Behavior In Microheterogeneous Environments

Alexander M. Berezhkovskii, Leonardo Dagdug, Sergey M. Bezrukov; Biophysical Journal 106(2):L09-L11, 2014; 10.1016/j.bpj.2013.12.013

參考資料