Metropolis蒙特卡羅方法、動力學蒙特卡羅方法、分子動力學方法這三種模擬方法有何特點與差異？

01-17

分子動力學是一個確定性方法，對於整個動力學過程給出模擬；Monte Carlo是一個基於隨機數的統計學方法；Metropolis Method是馬氏鏈蒙特卡洛方法(MCMC)。

這個問題很大，詳細展開足夠在一本分子模擬的教材里寫滿幾章了，了解它們最好的方法應該是去啃書，比如分子模擬的經典教科書："Understanding Molecular Simulation, From Algorithms to Application" by Frenkel Smit。

這裡先討論MD和Metropolis MC, 暫不討論Kinetic Monte Carlo (KMC)。這是因為傳統的MD和MC更具可比性（為簡單起見，這裡MC就是指狹義上的使用Metropolis演算法的蒙特卡羅，而這裡談到的MD僅限於平衡態的MD模擬，暫不討論非平衡態的NEMD），兩種方法產生的初衷都是為了計算統計力學裡的系綜平均：假設系綜所對應的相空間的概率分布為 $P(p,q)$ （p為動量，q為坐標），那麼對於任意平衡態熱力學量M, 它的系綜平均 $langle M angle$ 可以表示為

$langle M angle equiv int dp dq P(p,q) M(p,q)$ (1)

說得更直白一點，傳統的MD和MC基本上是在干同一件事：算積分。只不過最早的MD是在微正則系綜（NVE系綜）里算積分，而最早的MC是在正則系綜（NVT系綜）里算積分。

先說MC。眾所周知，正則系綜里相空間的概率分布就是玻爾茲曼分布

$P(p,q) = frac{e^{-eta mathcal{H}(p,q)}}{int dp dp e^{-eta mathcal{H}(p,q)} } = frac{e^{-eta mathcal{H}(p,q)}}{Z}$ (2)

其中 $mathcal{H}$ 是體系的哈密頓量， $Zequiv int dp dp e^{-eta mathcal{H}(p,q)}$ 是配分函數， $eta equiv 1/k_B T$ . 以一維線性諧振子 $mathcal{H}(p,q) = frac{1}{2}kq^2 + frac{p^2}{2m}$ 為例，這個相空間概率分布就是下圖所示，其實就是個高斯分布。

知道了概率分布，那麼正則系綜中任意平衡態熱力學量M的系綜平均 $langle M angle$ 就是下面這個積分

$langle M angle equiv int dp dq P(p,q) M(p,q) = frac{int dp dp M(p,q) e^{-eta mathcal{H}(p,q)}}{Z}$ (3)

對於（3）式中的積分，絕大多數情況下是沒有解析解的，因此只能通過數值求解。最簡單粗暴的用蒙特卡洛方法求解（3）式積分的方法就是在整個積分區域上均勻地隨機打點，求得每一點處的 $M(p,q)e^{-eta mathcal{H}(p,q)}$ 值，然後對所有點求和再除以配分函數即可。例如對於一維線性諧振子，我們可以如下在相空間上均勻地打點

這種演算法，的確非常簡單，但效率很低。因為位於相空間中 $e^{-eta mathcal{H}}$ 的值非常小的區域，系統處於這個區域的概率非常低，也就是說這些區域大體上來說對系綜平均 $langle M angle$ 的貢獻是非常小的，然而在採樣過程中，由於我們是隨機均勻打點，這樣 $e^{-eta mathcal{H}}$ 的值非常小的點和 $e^{-eta mathcal{H}}$ 的值非常大的點的權重是一樣的。換句話說，我們浪費了不少時間在那些權重很小的點上。那麼是否有這樣演算法，能夠使產生的點的分布符合玻爾茲曼分布，也就是說，在 $e^{-eta mathcal{H}}$ 的值非常小的區域少打點，在 $e^{-eta mathcal{H}}$ 的值非常大的區域多打點（如下圖所示），然後把每一點處的M值相加，就能得到我們需要的系綜平均呢？Metropolis演算法就是為了實現這樣的目標產生的。關於Metropolis演算法的細節，我不在這裡展開，很多書籍里都有描述。

再說MD。MD是基於牛頓方程的，而從牛頓方程是可以很自然導出孤立體系在保守力作用下是能量守恆的，因此最初的MD是在微正則系綜(NVE)里計算系綜平均。根據統計力學的等概率假設，對於總能量為E的微正則系綜，其相空間概率密度分布為

$P(p,q) = frac{1}{Omega (E)}$ (4)

其中 $Omega (E)$ 是體系處於能量為E時的態密度(density of state)

$Omega (E) = int dp dq delta (E - mathcal{H}(p,q))$ (5)

那麼熱力學量M在微正則系綜里的系綜平均 $langle M angle$ 就是

$langle M angle equiv int dpdq P(p,q) M(p,q) = frac{1}{Omega (E)} int_V dp dq M(p,q)$ （6）

其中V是積分的區域，也就是當能量為E時系統所佔據的相空間。拿一維線性諧振子作為例子，當能量為E時其在相空間中區域對應於下圖中的圓環，那麼這個系綜平均就是根據（6）式在對應於能量為E的圓環區域上進行積分。MD所做的，就是讓粒子跑遍這個圓環上的每一個點（歷經各態，ergodicity)，把每一點上M的值加起來求平均，得到的就是系綜平均(time average equals ensemble average, this is what we call ergodicity)。

綜上，MD和MC都可看作是為計算系綜平均的採樣方法。採用Metropolis演算法的MC不局限於在正則系綜里的採樣，可以應用於其他系綜，比如NPT和巨正則系綜。至於微正則系綜，理論上也可以使用Metropolis演算法，只是由於各個態之間能量相等，概率也相等，這時使用Metropolis演算法不能起到加速採樣的作用，等於白用。同樣地，MD的採樣也不局限於微正則系綜，當把一個孤立體系和一個恆溫器(thermostat)或者恆壓器(barostat)耦合起來，就可以進行正則系綜或者NPT系綜里的模擬。當然，如果耦合的恆溫器或者恆壓器是以extended Lagrangian形式實現的話（比如Nose-Hoover Thermostat), 廣義上來講MD模擬還是可看作在一個大的微正則系綜中進行。

MD和MC的比較

演算法的難易程度。一般說來，MD的演算法比MC要複雜一些。這是因為，對於MC而言，每產生一幀新的構型只要計算能量然後根據acceptance ratio移動粒子就足夠了。但對於MD，每產生一幀新的構型所需要的不只是能量，還需要每個粒子的速度和受力來對牛頓運動方程進行數值積分，這就加大了計算量。如果還要考慮為了提高計算效率的Verlet neighbor list等設置，MD的演算法會更複雜一些。
所能計算性質的多少。很明顯，MD能比MC計算更多的性質，MC只能計算一些與動力學無關的平衡態性質，比如能量、比熱容等等，MD既可以計算與動力學無關的性質，也可以計算與動力學相關的性質，比如擴散係數、熱導率、黏度、態與態之間的躍遷速率，光譜等等。這是因為，使用Metropolis演算法的MC，為了快速得到平衡態的分布，把整個系統簡化為一條馬爾可夫鏈，第N個構型只和第N－1個構型有關聯，與第N－2或者更早的構型沒有什麼關聯。在實際情況中不是這樣的，系統是有記憶的，在 $t_0$ 時刻處於相空間中 $(p_0,q_0)$ 點的系統與在 $t_1$ 時刻處於 $(p_1,q_1)$ 的系統是有關聯的，它們之間的關聯可以用條件概率 $P(p_1,q_1;t_1|p_0,q_0;t_0)$ 來描述。很多動力學性質都和這樣的條件概率有關係，而MD產生的數據是可以計算這些概率分布的。所以世界上沒有免費的午餐，MC雖然演算法簡單，但大體上只能計算靜態性質，如果要計算動力學性質還是要上MD。
可使用的軟體。可使用的MD的軟體遠多於MC的軟體。我個人覺得很大程度上是由於MD能夠計算的性質遠比MC多，並且數據可視化後給人更多直觀的感受。比如說，你可以把MD模擬蛋白質摺疊的軌跡文件用軟體可視化之後做成動畫，動畫就比較直觀地演示了蛋白質結構隨時間的變化。可是如果你用MC做同樣一個體系的模擬，你把軌跡可視化之後看到的可能就是一堆原子像無頭蒼蠅一樣亂動，最後穩定在某個構型附近。一般來講，MD軟體的開發也不難，一旦把積分牛頓方程的engine寫好，往裡面添加新的相互作用勢並不難，這樣每添加一點新的功能，MD軟體就可以適用於更多的體系。相比起MD軟體的百花齊放，分子模擬中學術界通用的MC軟體並不多，很多課題組都是自己寫MC代碼，但適用範圍僅限於自己研究的體系。
增強抽樣(enhanced sampling)。這個問題 @JimKarrey 在他的回答中已經提到了。這是分子模擬中的一個大坑，不管MD還是MC都會面臨在系綜中採樣不足的問題。這方面知乎上討論已經很多，可以參考下面的鏈接，在此不做詳細描述。

計算化學領域中有哪些技術可以被稱為是當前的黑科技？ - 知乎用戶的回答

兔子走遍世界——物理學視角下的採樣方法 - 生命的設計原則 - 知乎專欄

最後說一下KMC，我在這方面不是很懂，期盼有高人來解答。

Kinetic Monte Carlo(KMC)和主要模擬平衡態性質的MD和Metropolis MC相差就比較遠了，它主要模擬系統隨時間的演化而非平衡態性質。KMC的基本思想是把整個系統劃分成若干個態(state)，並獲得態與態之間單位時間的躍遷幾率（這個可以通過分子模擬或者過渡態理論來計算），然後通過具體的演算法決定每一步往哪一個態躍遷，就可以模擬系統在各個態之間隨時間演化。我沒有做過KMC的研究，並不太了解KMC的演算法，僅僅粗略看過一兩篇用KMC模擬擴散的文章。以下的例子是我粗淺的理解，可能不太準確。

我們想用KMC模擬粒子在二維晶格上的擴散，並求出擴散係數。假定在單位時間步長內，粒子只能跳到它最相鄰的格點。如圖所示，粒子某時刻處於A格點，那麼粒子下一時刻只有可能跳到B,C,D,E格點，單位時間內粒子從A跳到這些格點的躍遷幾率分別為 $r_{AB},r_{AC},r_{AD},r_{AE}$ 。結合這些躍遷幾率，KMC的演算法就可以通過生成隨機數來決定粒子下一步所在的格點。在模擬足夠長時間後我們可以算出粒子的均方位移(mean square displacement, MSD)，均方位移和擴散係數成正比關係，因而可以求出擴散係數。

當粒子只能在一維情況下運動（即只能從A到B或者從A到D), 並且 $r_{AB} = r_{AD}$ ，這種情況下KMC其實就是在模擬一維隨機遊走(random walk)。

KMC的具體應用應該還很多，比如晶體生長什麼的。

MD和MC了解不深，上面也有比較詳細的回答了，我就不丟人了。KMC這塊我稍微了解一些，自己也寫了幾千行的KMC代碼，說說自己的理解吧。

1、特點：

Kinetic Monte Carlo (KMC)是專門用於模擬體系長時間演化的一種MC方法，它的特點是每一步都對應著某個原子從一個組態（能量極小態）遷移到另一個組態，例如缺陷擴散、團簇分解等反應。兩個組態之間的躍遷勢必要跨越一個能壘 $E_b$ ，而跨越這個能壘的能量是由熱擾動提供的。因此發生兩次躍遷的平均時間間隔正比 $exp(E_b/k_bT)$ 。兩次躍遷之間，原子只是在能量極小位置附近振動，這類振動是KMC不關心的，是被忽略掉的。

很顯然，KMC適用於模擬「由組態躍遷主導的體系的長時間演化」，例如缺陷在晶體內部和表面的擴散、聚集、分解等。而類似於蛋白質摺疊這種並沒有明顯組態變化的物理體系，就不適合用KMC來進行模擬。

2、演算法簡介：

核心思想和MC還是一樣的，利用偽隨機數來模擬熱擾動，只不過在KMC中僅僅考慮引起組態躍遷的熱擾動而已。

一般情況下，KMC模擬中每個可能的躍遷能壘都要提前輸入（有的演算法會調用分子靜力學等方法實時計算能壘）。在確定所有躍遷途徑的能壘後，給出每個躍遷的反應速率，例如： $v^i=v_0exp(-E_b^i/k_bT)$

然後根據速率隨機選擇躍遷途徑，選中第i個躍遷途徑的概率為：

$p(i)=v^i/sumlimits_{j}v^j$

執行該躍遷反應，更新反應帶來的躍遷途徑及其能壘的變化，並按下式增加模擬時間：

$Delta t=1/sumlimits_{j}v^j$

以上為一個完整的Monte Carlo步，程序會重複以上過程，直至達到所需時間或者步數為止。

3、用途：

前面已經說了，KMC適用於模擬「由組態躍遷主導的體系的長時間演化」。最常用的領域有晶體缺陷演化和晶體生長。KMC模擬的時間尺度主要與溫度以及體系中最低的躍遷能壘有關，一般來說，KMC模擬是能達到宏觀時間尺度（秒~年）的。

以ising模型的模擬為例，傳統的演算法會進行一次次反轉自旋的嘗試，根據metropolis判據選擇是否接受該嘗試。然而在低溫下，這種嘗試大部分是被拒絕掉的。假設此時拒絕的概率為999/1000，那麼metropolis演算法要進行1000次嘗試才能得到一次有效反轉，而KMC則相反，它不管三七二十一先反轉了，然後根據反轉概率（躍遷速率），給出此次反轉所需時間期望值，即1000次嘗試對應的時間。這樣能有效濾掉無用嘗試，大大增加模擬的時間尺度。

當然有利也有弊，由於體系中除去組態躍遷之外的所有信息都被忽略掉了，因此要求體系的演化軌跡必須是由這些已知的組態躍遷主導，而無法考慮未知的組態躍遷。舉例來說，假設我們用KMC 模擬單個空位V在晶體內的擴散係數，此時考慮的組態躍遷只有空位遷移一個。然而實際情況下，晶體中可能有位錯、晶界等缺陷，這些缺陷能與空位相互作用影響其擴散係數，如果沒有提前將該機制納入考慮的話，模擬中是體現不出該影響的。因此，正確的設計模擬所關心的物理機制是KMC模擬的核心。

總的來說，KMC一般不是用來做機理研究的，而是在已知組態躍遷機制的情況下，模擬多種躍遷機制對體系長時間演化軌跡的影響。它常與DFT、MD、MS等較為「精細」的方法結合，形成一個從原子到微米尺度，從飛秒到年尺度的多尺度模擬方法。例如NM上的這篇文章：http://www.nature.com/nmat/journal/v4/n1/abs/nmat1286.html

謝邀。

做模擬的目的一般來講，一個是理解、觀察過程的微觀機理，再一個是從微觀狀態統計某些過程的宏觀量。

首先說分子動力學（MD）模擬，他的基本框架就是經典力學的運動方程積分的確定論的方法，但體系本身的演化是一個自然的馬爾可夫鏈（多體過程就不要想著有確定的解析解了）。MD在恆溫恆壓控制的演算法方面也會引入大量隨機過程。常見的是Langevin dynamics （恆溫）和Langevin piston（恆壓），做MD模擬的看到這兩個名詞簡直再熟悉不過。

做MD的模擬的目的是什麼？首先是在微觀層面上觀察到分子運動，進而可以理解某些過程的機理，這一點是MD的巨大優勢。其次，但也是最重要的，MD需要實現各態遍歷性（Ergodicity）從而能夠給出一些宏觀統計量的值。每一個宏觀統計量，像自由能，熵，焓，擴散係數，結合常數等等，都對應著某些微觀量的系綜統計平均。要得到一個準確的系綜統計平均，就要實現各態遍歷性。如何保證一個MD模擬可以實現各態遍歷？這個問題直到二十世紀末才直到得到回答。一個MD方法應該具有辛性（Symplecticity，不要問我這個詞什麼意思，它就是音譯）才能遍歷各態，到這一步，對MD的基礎理論的研究已經屬於數學上的辛幾何了，有志青年可以研究研究，個人感覺經典力學發展到這一步，一點都不比量子力學簡單，數學形式上可以還要更加艱深。

解決了各態遍歷的可實現性，現在的問題就是如何實現各態遍歷性，或一定條件下的（偽）各態遍歷。於是現在就有了多種加速採樣方法。包括像Umbrella Sampling, Metadynamics, Adaptive Biasing Force之類的方法。個人已被ABF玩死。這些演算法的確是「加速」了，但是想要得到一個好的統計結果，需要下花很多工夫，操作起來比直接的MD模擬麻煩多了。

MD是很好，結果準確（相對蒙卡而言），過程直觀。缺點就是，計算資源消耗極大，做計算一定要找個土豪錢多的組，多買或者多租計算節點，要不然就會像我這樣被坑慘了。

再說蒙卡方法，這是一類基於概率論的方法。我沒有真正做過蒙卡的模擬，就憑我的了解說一點。講錯了輕拍。

大家都熟悉的一個蒙卡方法的例子是這樣的：一張正方形白紙上畫了一個圓，現在要求計算這個圓的面積。正方形的面積是容易計算的，很好辦。圓就難多了。採用蒙卡的思想的話，可以往紙上隨機地扔圖釘。當圖釘數量很多時，可以統計在圓內和圓外的圖釘的比例，從而大致估計圓的面積佔據整張正方形白紙面積的比例，於是就能估算得圓的面積了。這裡，圓的面積就是我們要求得的統計量，而扔圖釘則是我們產生的隨機過程。大夥還可以用類似的方法寫個程序計算圓周率的值。

蒙卡方法應用於分子模擬中，則是這樣的：隨機產生體系的狀態，在每個狀態（即每一次採樣）中計算統計量的值。當採樣數足夠多的時候，計算這些統計量的Boltzmann系統平均，然後就得到了相應的宏觀物理量的估計值。因此可以看到，蒙卡方法應用於分子模擬中，是直接從各態遍歷性下手的。產生一堆狀態，並不太管體系在這些狀態里演化的每一步的細節。因此，蒙卡很難給出一個直觀的體系在狀態間演化的過程，而是需要你根據蒙卡的統計結果間接推測演化的路徑和機理。

由於蒙卡是隨機方法，所以每次估計出來的值都不會一樣。但對於一個已經斂的採樣過程，再增加大量的採樣數目，最後的數值也只會在一個較小的合理區間內波動。

一般的蒙卡方法應用於簡單的分子體系是很好辦的。比如一個純水的盒子這樣的。隨機產生水分子的分布就可以了。但是這種簡單的蒙卡方法的局限性也是非常明顯的，即它會產生大量的不合理的體系狀態，這些狀態在統計上的貢獻極低，從而浪費了大量的計算資源。我們知道，一個物理量F的系統平均是這樣計算的：

$<F> = frac{int{F(p,q)ullet e^{-frac {E(p,q)} {k_{B}T}} }} {int e^{-frac {E(p,q)} {k_{B}T}}}$

能量越高的地方，權重越低。要知道，一個體系的自然演化路徑，總是天然地要找能量最低路徑。只有這些能量最低路徑上的採樣才是最有用，最有效的。於是就有了Metropolis蒙卡方法。

Metropolis蒙卡方法則是先得到體系的一個合理狀態，然後儘可能地沿著一個合理的路徑進行採樣，從而大大提高採樣的效率。

下圖是我從一部經典的分子動力學教材中截取的插圖（來自《Understanding Molecular Simulation, From Algorithms to Application》，作者Daan Frenkel, Berend Smit）：

想要測量尼羅河的平均深度，我們不必在整個地圖上隨機采點測量（左）。站在河裡面，沿著合理的路徑（河道）來採樣才是最高效的。

簡單粗暴的蒙卡方法對於生物蛋白質這樣的體系來說是非常致命的。簡單粗暴地隨機產生構象的話，很容易就破壞了蛋白質的二級結構，而你希望採樣到的正確的構象則往往無法恢復。這時候，如果考慮蛋白質本身的行為性質，循著它本身構象變化的路徑採樣，那樣才會得到比較正確的結果。（這一段完全是我的推測，期待用蒙卡做過蛋白質模擬的人來細說。）比較MD和MC，除了理論框架的確定論和概率論的差別以外，一個重要的區別是，MD天然地產生一個馬爾可夫鏈，而MC的馬爾可夫性是需要額外的方法來補充的。

本科畢設就是編程實現蒙特卡洛模擬稀薄流體。

分子動力學與蒙特卡洛模擬放大的區別就是100％與95%的差別。

後來又了解到飛行力學中的蒙特卡洛方法計算飛行器軌跡，不禁覺得概率的思想很重要，保證覆蓋性的同時大大降低了計算量，加快了計算時間。

ps知乎上各種高人太多了，這麼偏的問題都有人回答...

分子動力學就是牛頓運動方程的方法，取的時間步驟越短，越精確。

蒙特卡洛是另一種方法，先產生大量大量的樣本。

metropolis是一種從樣本中選擇結果的方法

動力學蒙卡是另一種選擇結果的方法。（你說得動力學蒙卡是指基於格林函數的蒙卡嗎？）

取得樣本越大，越精確。

以上幾個回答基本上把這幾種模擬方法都解釋的差不多了，懂的自然懂，不懂的刷知乎也沒啥意義了，趕緊找幾本演算法或simulation的書看看吧～上面都有詳細的介紹。

我想補充的是kinetic monte carlo（KMC）的精確翻譯應該是運動學蒙特卡洛，因為如果是動力學dynamic 的話，意味著相鄰兩個步態一定是動量能量守恆的。而kmc顯然無法做到這一點。所以發明這個演算法的人用的是kenetic 而非dynamic 。我導師還提議過動理學蒙特卡洛這個名字，不過感覺有些中二～

MC和MD都是抽樣方法，一個是隨機的模擬，一個是確定性的模擬；

MD能研究動力學過程，非平衡體系；MC有難度，即使是動力學，也是偽動力學。

想了解MD的基本原理，這有個網易雲課程，可以看看：分子動力學模擬入門引導網路課程

只用過密度泛函方法_(:з」∠)_