1、Chapter 7 Monte Carlo及相關(guān)方法7.1 Monte Carlo模擬7.2 分子動力學(xué)模擬7/11/2022第七章 Monte Carlo方法27.1 Monte Carlo模擬派什么用場；用計(jì)算機(jī)來模擬統(tǒng)計(jì)的隨機(jī)抽樣(Sampling)過程有非常廣的應(yīng)用，只要能把問題轉(zhuǎn)化為隨機(jī)抽樣，即可應(yīng)用應(yīng)用于事件的幾率已知，但結(jié)果很難預(yù)知的過程名字的由來Monaco（摩納哥）公國的一個城市，以賭博聞名。從大約1944年開始發(fā)展和完善7/11/2022第七章 Monte Carlo方法37.1 Monte Carlo模擬7/11/2022第七章 Monte Carlo方法4Monte C

2、arlo的要素Probability distribution functions (pdfs)研究的對象必須可以轉(zhuǎn)化為一些幾率分布函數(shù)來描述隨機(jī)數(shù)(Random Number)產(chǎn)生均勻分布的等間隔的隨機(jī)數(shù)Sampling rules （抽樣規(guī)則）如何取樣（物理上的測量如何進(jìn)行）7/11/2022第七章 Monte Carlo方法5Monte Carlo的要素, continued.如何評價(scoring)感興趣的量把每次取樣綜合起來得到感興趣的量的值（物理量的平均值）誤差估計(jì)與隨機(jī)取樣的樣本數(shù)目或是其它量的關(guān)系盡量減少變量數(shù)目以節(jié)省Monte Carlo模擬的時間并行化的考慮蒙特卡羅方法的

3、收斂性，誤差 蒙特卡羅方法作為一種計(jì)算方法，其收斂性與誤差是普遍關(guān)心的一個重要問題。收斂性誤差減小方差的各種技巧 效率收斂性 蒙特卡羅方法是由隨機(jī)變量X的簡單子樣X1，X2，XN的算術(shù)平均值： 作為所求解的近似值。由大數(shù)定律可知，如X1，X2，XN獨(dú)立同分布，且具有有限期望值（E(X)），則 即隨機(jī)變量X的簡單子樣的算術(shù)平均值 ，當(dāng)子樣數(shù)N充分大時，以概率1收斂于它的期望值E(X)。誤差 蒙特卡羅方法的近似值與真值的誤差問題，概率論的中心極限定理給出了答案。該定理指出，如果隨機(jī)變量序列X1，X2，XN獨(dú)立同分布，且具有有限非零的方差2 ，即 f(X)是X的分布密度函數(shù)。則 當(dāng)N充分大時，有如下

4、的近似式 其中稱為置信度，1-稱為置信水平。 這表明，不等式 近似地以概率 1成立，且誤差收斂速度的階為O(N-1/2) 。 通常，蒙特卡羅方法的誤差定義為上式中與置信度是一一對應(yīng)的，根據(jù)問題的要求確定出置信水平后，查標(biāo)準(zhǔn)正態(tài)分布表，就可以確定出。幾個常用的與的數(shù)值：關(guān)于蒙特卡羅方法的誤差需說明兩點(diǎn)：第一，蒙特卡羅方法的誤差為概率誤差，這與其他數(shù)值計(jì)算方法是有區(qū)別的。第二，誤差中的均方差是未知的，必須使用其估計(jì)值來代替，在計(jì)算所求量的同時，可計(jì)算出 0.50.050.003 0.67451.963減小方差的各種技巧 顯然，當(dāng)給定置信度后，誤差由和N決定。要減小，或者是增大N，或者是減小方差2。

5、在固定的情況下，要把精度提高一個數(shù)量級，試驗(yàn)次數(shù)N需增加兩個數(shù)量級。因此，單純增大N不是一個有效的辦法。另一方面，如能減小估計(jì)的均方差，比如降低一半，那誤差就減小一半，這相當(dāng)于N增大四倍的效果。因此降低方差的各種技巧，引起了人們的普遍注意。后面課程將會介紹一些降低方差的技巧。 效率 一般來說，降低方差的技巧，往往會使觀察一個子樣的時間增加。在固定時間內(nèi)，使觀察的樣本數(shù)減少。所以，一種方法的優(yōu)劣，需要由方差和觀察一個子樣的費(fèi)用（使用計(jì)算機(jī)的時間）兩者來衡量。這就是蒙特卡羅方法中效率的概念。它定義為其中c是觀察一個子樣的平均費(fèi)用。顯然 越小，方法越有效。 7/11/2022第七章 Monte Ca

6、rlo方法13一般程序產(chǎn)生隨機(jī)數(shù)轉(zhuǎn)變?yōu)殡S機(jī)變量（抽樣）利用變量來模擬產(chǎn)生結(jié)果Analyze Output是否繼續(xù)7/11/2022第七章 Monte Carlo方法141.隨機(jī)數(shù)的產(chǎn)生真正的隨機(jī)數(shù)是很難產(chǎn)生的偽隨機(jī)數(shù)(Pseudo Random Number)產(chǎn)生的方法丟骰子，擲硬幣搖獎做好的表格計(jì)算機(jī)遞推方法7/11/2022第七章 Monte Carlo方法15對隨機(jī)數(shù)發(fā)生器(Random Number Generator, RNG)的要求產(chǎn)生0,1間均勻分布的隨即序列各隨機(jī)數(shù)間是獨(dú)立的數(shù)目巨大，很少重復(fù)可重復(fù)得到有效，并且容易移植7/11/2022第七章 Monte Carlo方法16R

7、NG-乘同余法Zi = mod(aZi-1+b,m) m: 通常是計(jì)算機(jī)可以表示的最大正整數(shù)，如(231 1)a 正的奇數(shù)，如 75 = 16807b 步長，整數(shù)Zo (Starting) = Zi-1 -初值，種子7/11/2022第七章 Monte Carlo方法17RNG-乘同余法隨機(jī)數(shù)i注意，取a和m時須小心，避免循環(huán)7/11/2022第七章 Monte Carlo方法18乘同余法的改進(jìn)：混合偽隨機(jī)數(shù)序列乘同余法的缺點(diǎn)：偽隨機(jī)數(shù)之間的相關(guān)性比較強(qiáng)改進(jìn)混合偽隨機(jī)數(shù)首先產(chǎn)生一組整數(shù)隨機(jī)數(shù)列，然后從其中取出一個整數(shù)，根據(jù)這個整數(shù)再產(chǎn)生新的隨機(jī)數(shù)列結(jié)果相當(dāng)于將一個偽隨機(jī)數(shù)列按另一個偽隨機(jī)數(shù)列的

8、大小重排，因而減弱了所用的隨機(jī)數(shù)序列的相關(guān)性。缺點(diǎn)：增加計(jì)算機(jī)時7/11/2022第七章 Monte Carlo方法192.把隨機(jī)數(shù)變?yōu)樽兞侩S機(jī)數(shù)與隨機(jī)變量要一一對應(yīng)逆變換的方法產(chǎn)生一個隨機(jī)數(shù)由逆函數(shù)X=F-1()求得變量X011020300.77/11/2022第七章 Monte Carlo方法20分布函數(shù)種類均勻分布指數(shù)分布高斯分布三角形分布7/11/2022第七章 Monte Carlo方法21均勻分布(U-L)+L : 隨機(jī)數(shù)L: 變量的下限；U: 變量的上限例：在公共汽車站臺某一路車到達(dá)的時間7/11/2022第七章 Monte Carlo方法22三角形分布7/11/2022第七章

9、Monte Carlo方法23指數(shù)型分布如Boltzmann分布7/11/2022第七章 Monte Carlo方法24高斯型（正則）分布如考試成績7/11/2022第七章 Monte Carlo方法25例子：對函數(shù)的積分定域積分：產(chǎn)生隨機(jī)數(shù)， 0,1抽樣：得到隨機(jī)變量f()，根據(jù)模擬：得到感興趣的量7/11/2022第七章 Monte Carlo方法26例子：對函數(shù)的積分7/11/2022第七章 Monte Carlo方法273. 抽樣方法：理論上U可取到無窮大，實(shí)際很大的值對配分函數(shù)的貢獻(xiàn)值為0分布函數(shù)重要性實(shí)際體系配分函數(shù)的構(gòu)型積分e-U/kTU7/11/2022第七章 Monte Ca

10、rlo方法28重點(diǎn)抽樣法G(x)滿足G(1)=1, G(0)=0，將G(x)看作是區(qū)間0,1上的分布函數(shù)，則g(x)為幾率密度函數(shù)有目的的抽樣減小方差，加速收斂7/11/2022第七章 Monte Carlo方法29重點(diǎn)抽樣法，continued積分值方差如7/11/2022第七章 Monte Carlo方法30重點(diǎn)抽樣法，continuedBoltzmann分布的重點(diǎn)抽樣7/11/2022第七章 Monte Carlo方法31Boltzmann分布的重點(diǎn)抽樣，continued現(xiàn)在，分布函數(shù)變成了7/11/2022第七章 Monte Carlo方法32Boltzmann分布的Metropol

11、is抽樣方法思想：是抽樣盡可能按實(shí)際（也就是平衡態(tài)）的分布進(jìn)行做法：構(gòu)建一個Markov鏈，每改變一個構(gòu)型，即改變體系中分子的坐標(biāo)pi,qi，改變一次，即為走一步(A Step)每走一步，得到相互作用能(勢能）U，即判斷是否合理，保留合理的走法求物理量的期望值7/11/2022第七章 Monte Carlo方法33Metropolis抽樣方法產(chǎn)生隨機(jī)數(shù)0,1改變每個分子的坐標(biāo)步長計(jì)算得到體系的勢能Unew及勢能變化Unew半經(jīng)驗(yàn)的方法計(jì)算勢能，Lenard-Jones勢函數(shù)精確從頭算計(jì)算7/11/2022第七章 Monte Carlo方法34Metropolis抽樣方法, continue另外

12、再產(chǎn)生隨機(jī)數(shù)rand(0,1)判斷true，保留改變，繼續(xù)對新坐標(biāo)改變false，重新改變構(gòu)型完成足夠數(shù)量的循環(huán)次數(shù)后，得到體系物理量（熱力學(xué)函數(shù)）的平均值7/11/2022第七章 Monte Carlo方法35Metropolis抽樣方法, continue步長rmax的注意事項(xiàng)不能太大，否則很多次改變不合理不能太小，否則效率差，需走很多步才能探索完整個構(gòu)型空間7/11/2022第七章 Monte Carlo方法364. 分子的Monte Carlo模擬比原子體系復(fù)雜改變分子質(zhì)心坐標(biāo)的同時還要考慮分子的旋轉(zhuǎn)不成功的step更多剛性分子只須考慮分子的整體旋轉(zhuǎn)繞x軸選旋轉(zhuǎn)角的坐標(biāo)變換7/11/2

13、022第七章 Monte Carlo方法37分子的Monte Carlo模擬，continue任意旋轉(zhuǎn)，Euler angles7/11/2022第七章 Monte Carlo方法38分子的Monte Carlo模擬，continue任意旋轉(zhuǎn)的坐標(biāo)變換實(shí)際上對cos的做法是：7/11/2022第七章 Monte Carlo方法39分子的Monte Carlo模擬，continue使用Euler角的缺點(diǎn)上面的矩陣含有6個三角函數(shù)，計(jì)算機(jī)計(jì)算耗時解決辦法：變量替換7/11/2022第七章 Monte Carlo方法40分子的Monte Carlo模擬，continue坐標(biāo)變換矩陣即變?yōu)椋阂虼诵枰a(chǎn)

14、生四個隨機(jī)數(shù)實(shí)現(xiàn)對q0,q1,q2,q3進(jìn)行改變：產(chǎn)生-1,1間的隨機(jī)數(shù)1, 2，使S1= 12+221；產(chǎn)生-1,1間的隨機(jī)數(shù)3, 4，使S2= 32+421；變換的基失(q0,q1,q2,q3)為：7/11/2022第七章 Monte Carlo方法41柔性分子或聚合物高分子的模擬分子的內(nèi)部自由度太多，非常困難固定某些內(nèi)部自由度，如鍵長和鍵角，而放開兩面角高分子的模擬請參考相關(guān)專著，如：楊玉良, 張紅東 著：高分子科學(xué)中的 Monte Carlo 方法，復(fù)旦大學(xué)出版社，19937/11/2022第七章 Monte Carlo方法425. “Biased” Monte Carlo Metho

15、dsExplore the most important part of the phase spaceSingle solute molecule surrounded by a large number of solvent moleculesPreferential sampling接近溶質(zhì)分子的溶劑分子移動比離得遠(yuǎn)的快7/11/2022第七章 Monte Carlo方法43Preferential sampling選取一個區(qū)間(cutoff region)，在這個區(qū)間外的分子移動慢而區(qū)間內(nèi)的分子移動快選取一個幾率p0，對一個分子，如果在cutoff region內(nèi)則移動反之，產(chǎn)生一個0

16、-1間的隨機(jī)數(shù)p1p1 p0，不動Cutoff region的選取也可以不固定7/11/2022第七章 Monte Carlo方法44Preferential sampling, continue需要注意的問題移出cutoff區(qū)間內(nèi)的分子比移進(jìn)來的多Force biased Monte Carlo原子或分子按照它所受的力的方向移動Calculate forcesA: 常數(shù)fi: 力riG: 按Gaussian分布，且平均值為0，variance為：2A7/11/2022第七章 Monte Carlo方法456. 不同系綜的Monte Carlo方法最簡單的就是正則系綜即NVT系綜NPT系綜：改

17、變原子或分子坐標(biāo)的同時改變體積體積改變后，體系的勢能也要重新計(jì)算，對于簡單的相互作用勢能，下面的技巧：scaled coordinates是很有用的（Lold是體系所處立方體的邊長）:7/11/2022第七章 Monte Carlo方法46NPT系綜上式可以重新寫為：體積改變后體系的新勢能為：7/11/2022第七章 Monte Carlo方法47NPT系綜, continue則勢能變化：對于NPT系綜，判斷一個step要不要繼續(xù)是根據(jù)H：7/11/2022第七章 Monte Carlo方法48巨正則系綜(VT系綜）化學(xué)勢可寫作z：活度，常數(shù)de Broglie 波長7/11/2022第七章

18、Monte Carlo方法49巨正則系綜, continue需要用到隨機(jī)數(shù)作改變的量(取樣)：粒子的運(yùn)動，用Metropolis方法取樣粒子的湮滅(destroy)粒子的生成(create)7/11/2022第七章 Monte Carlo方法50巨正則系綜, continue一些注意事項(xiàng)粒子產(chǎn)生的幾率和粒子湮滅的幾率要一樣對于一些比較密集的體系，產(chǎn)生和湮滅成功的幾率比較小7/11/2022第七章 Monte Carlo方法517. 計(jì)算化學(xué)勢按照化學(xué)勢的定義N很大7/11/2022第七章 Monte Carlo方法52化學(xué)勢，continue對于實(shí)際的氣體或液體，氣體或液體的構(gòu)型積分：單個分子

19、（理想氣體）配分函數(shù)7/11/2022第七章 Monte Carlo方法53化學(xué)勢，continue則化學(xué)勢，令，7/11/2022第七章 Monte Carlo方法54化學(xué)勢，continue化學(xué)勢，7/11/2022第七章 Monte Carlo方法55: Configurational bias Monte Carlo Method參考：Siepmann, J. I. A Method for the Direct Calculation of Chemical Potentials for Dense Chain Systems. Mol. Phys. 1990, 70(6), 1145-1158.Siepmann, J. I.; Frenkel, D. Configur