分子模擬及分子動力學簡介

1、 什么是分子模擬 分子模擬是在分子模型的基礎上用計算機做實驗，“計算機實驗” 通過模擬微觀粒子的運動來計算宏觀性質(zhì)溫度壓力黏度傳遞性質(zhì)表面張力分子間的作用模型牛頓力學量子力學統(tǒng)計力學等 分子模擬具有理論和實驗的雙重性質(zhì) 分子模擬不能完全取代實驗理論實驗模擬理論的正確性模擬參數(shù)的正確性模擬方法的選擇理論的更新與自然界相比的準確程度尺度（米）時間（秒）10-910-710-510-310-1510-910-610-31電子模擬（量化計算，dft）分子模擬（分子動力學，蒙特卡洛）顆粒方法流體力學原子結(jié)構薛定諤方程模擬電子云能量性質(zhì)，化學鍵等信息 量子化學計算一般處理幾個到幾十個原子常見軟件：gaus

2、sian，nwchem等 密度泛函（dft）可以算到上百個原子常見軟件：vasp 分子水平的模擬 以分子的運動為主要模擬對象 采用經(jīng)驗性的分子間作用函數(shù)模擬微粒之間的作用 一般情況下不考慮電子轉(zhuǎn)移效應，因而不能準確模擬化學成鍵作用 發(fā)展最早 1950s，alder，勞倫斯利物默實驗室，分子動力學模擬32個原子 1950s，metropolis，洛斯阿洛莫斯實驗室，蒙特卡洛模擬32個原子 分子級別的模擬應用的領域很廣 廣泛應用于化學，物理，生物，化工，材料，機械，治藥等領域 簡單易學 蒙特卡洛是一種優(yōu)化方法通過蒙特卡洛算法來尋求能量最優(yōu)點隨機方法 通過系綜平均來求取宏觀性質(zhì) 模擬的是平衡狀態(tài)，不

3、涉及時間效應（kmc除外） 優(yōu)點是可以跨越時間因素，缺點是得不到有關時間信息的性質(zhì) 可以模擬平衡狀態(tài)，也可以模擬中間狀態(tài) 可以獲得有關時間的信息 受時間的限制，無法模擬緩慢過程分子體系（幾百幾億）求解牛頓運動方程宏觀性質(zhì)( )du rm rdr= -& & 同時考慮原子核的運動（牛頓力學）和電子的運動（量子力學） 能同時準確模擬物理作用和化學鍵作用 目前來說cpmd可以處理的體系還很?。◣资畟€原子） 將分子基團（幾個或者幾十上百個原子）當成單個的微粒來處理 微粒之間的作用也是通過類似于分子動力學的未能函數(shù)來描述 可以模擬更長的時間跨度電子原子核原子量子級別模擬分子級別模擬cg級別模擬分子動力學

4、對勢能函數(shù)的依賴性：所有從分子動力學計算出來得到的宏觀性質(zhì)最終都取決于勢能模型( )du rm rdr= -& &分子動力學的核心：牛頓運動方程勢能（位能）模型：( )( )u rf r=rur例：甲烷，某些惰性氣體質(zhì)點處理ur方阱模型ur階梯模型鍵的振動鍵角扭矩分子內(nèi)部各原子（基團）之間的范德華力、靜電力一般要計算1-4（相隔超過兩個鍵的原子或基團對）15432qqq分子之間的范德華力分子之間的靜電力ccchhhhhhhh 10根鍵長作用 18個鍵角作用 8個扭矩作用 27個范德華力作用 27個靜電作用20()1exp() eu lda ll=- morse類鍵長模型能量阱深參數(shù)鍵長平衡鍵長

5、參數(shù) 胡克類鍵長模型20()()2ku lll=-鍵長平衡鍵長參數(shù) 胡克類鍵角模型20()()2ku lqq=-鍵角平衡鍵角參數(shù)0( ) 1cos()2nnnvunwwg=+-扭矩障礙參數(shù)扭動360度所經(jīng)過的能量最低點的次數(shù)126126( )4u rrrsse驏=-桫lennard-jones模型ur 混合規(guī)則abab1()2a babsx ss=+a babez e e=通常都取1 一般情況下只考慮點電荷之間的作用力 不考慮極化作用所帶來的長程項的作用114abnni jijijqqurpe=邋分子a的點電荷總數(shù)分子b的點電荷總數(shù) 用有限的微觀分子體系模擬實際宏觀體系的必要手段3214342

6、121計算周期邊界條件下兩個微粒之間的作用宏觀體系微觀體系主盒影子影子影子 ewald方法 particle mesh ewald 方法（pme) 計算量很大bac( , , )( , )( , )( , )u a b cu a bu a cu b c+ 處理方法(2)(3)( )( , ,)totalijijikjkuurur r r=+邋邋 通過量子化學模擬回歸得到點電荷范德華力鍵長、鍵角、扭矩力 實驗數(shù)據(jù)回歸鍵長鍵角范德華力初始化能量優(yōu)化平衡數(shù)據(jù)產(chǎn)出避免局部分子重疊，并不是動力學模擬 根據(jù)所有分子的當前坐標計算個分子的受力（位能函數(shù)） 根據(jù)受力更新分子的坐標 在此過程中收集用來計算宏觀性

7、質(zhì)的有關信息讀入模型參數(shù)，模擬控制參數(shù)去除某些可能存在的原子重疊去除某些嚴重扭曲的鍵長、鍵角、扭矩等方法最速下降法牛頓拉夫森方法其他一般優(yōu)化幾千到幾萬步 verlet法簡單易行，但是有精度損失leap-frog法，verlet法的變種 速度verlet法和verlet法相比，可以同時獲得位置、速度、加速度，而且沒有精度損失beeman算法，速度計算精度更高，但是計算量大預測-校正法算法選擇精度和運算量的折衷體系原子體系剛性分子體系非剛性分子，剛性鍵非剛性分子涉及到的分子運動平動平動、轉(zhuǎn)動平動、轉(zhuǎn)動、扭矩平動、轉(zhuǎn)動、扭矩、振動建議步長10 fs5 fs2 fs-1 fs 不同性質(zhì)的作用力采取不同

8、的步長 最常見的方法：r-respa 成功用于多種體系，簡單流體，有機分子體系，蛋白質(zhì)等 和普通速度verlet法相比，對于復雜分子體系（高分子，生物分子等），r-respa可以大大加快運算速度，步長可以達到1-2 fs 為保證原子各自運動時分子的整體性而添加的分子內(nèi)部的約束條件算法 最常見的是針對氫原子，因為氫原子很高的振動頻率 verlet算法：shake 速度verlet算法：rattle nve，最簡單的分子動力學e=k+u恒定體積，而且和外界沒有能量或者物質(zhì)交換，只有動能和勢能之間的轉(zhuǎn)換通過監(jiān)視各能量項的變化，可以檢查程序是否基本正確21/2(3)niicbimknnk t=-p p

9、體系總動能和體系溫度的關系原子總數(shù)體系總的約束條件數(shù)自由度模擬開始時，體系中各微粒的初始速度也是可以根據(jù)這個關系計算微粒 i 的動量微粒 i 的質(zhì)量bpvn k tw=+壓力與分子間的作用勢能的關系理想氣體分子間的作用造成的剩余壓力1111( )33inniiiiiiwu= -籽=邋r rrrr frrr f維利系數(shù)w的計算擴散系數(shù)二元徑向分布函數(shù)2( )()nnijij ivg rnd=-邋rrrr1( )(0)6iidtt=-rrrrxtx通過斜率求取d實際位置，不是pbc后的位置 速度調(diào)節(jié) 簡單scale v=v*(t0/t(t)1/2 berendsen法 增加調(diào)節(jié)參數(shù) 重新隨機設置速

10、度每隔一定步數(shù)，按照控制溫度隨機抽取一組高斯速度分布熱量交換nose-hoover法將交換熱源也當成體系的一部分進行積分產(chǎn)生嚴格的nvt系綜控制體積活塞類似于溫度控制通過調(diào)節(jié)體系的體積來調(diào)節(jié)壓力壓力的控制比溫度更難，因為壓力和能量的微分量（力）成正比，每一步的波動更大 氣體的模擬 小分子體系，不需要復雜的勢能模型 幾百到幾千個分子，分子分布稀疏，大部分是短程作用 一般用一臺微機就可以處理，計算時間幾分鐘幾小時 簡單的液體，不涉及太多的界面性質(zhì) 小分子體系，勢能模型不是很復雜 幾百個分子，可能涉及到靜電作用，可能需要長程校正 用微機也可以處理，計算時間一般幾小時幾天 必要性 體系越來越大 模擬時

11、間越來越長 解決辦法 制造更快的處理器 并行計算機 例子：50000原子的生物體系，1ns模擬 單個處理器：12天 16個并行處理器：1天或者 message passing interface 90年代初制定和完善的一套并行語法 支持fortran, c, c+ 簡單易學 并行效率需要1小時需要1/2小時完美的并行效率 處理器的速度遠遠超過數(shù)據(jù)傳輸?shù)乃俣龋罅康臅r間花在處理器之間的信息傳遞上了 cpu的速度幾乎是幾何級數(shù)增長 內(nèi)存的速度是代數(shù)級數(shù)增長 加快數(shù)據(jù)傳輸，盡量減少花在數(shù)據(jù)傳輸上的時間 數(shù)據(jù)傳輸硬件上的進步 算法上做文章 每一個處理器負責處理一部分原子 每一步計算每一個處理器都要接受

12、其它處理器負責處理的原子的相關信息 信息傳輸量大，使用的處理器越多并行效率效率越低 一般適合處理5-10萬左右微粒的體系n=20n=15n=610n=1620n=11151201481291961613155271741131018實際情況 按照體系的實際物理位置按區(qū)域劃分每個處理器的處理范圍 每一步計算每一個處理器只需要和相鄰的處理器交換信息 數(shù)據(jù)傳輸量小，并行效率高，適合處理大型體系（超過10萬微粒） 算法比較復雜（邊界的處理）n=201201481291961613155271741131018 主要針對與生物和化學軟材料體系 優(yōu)點 程序設計水平高，計算效率高，號稱可以有效并行到上千個處

13、理器 兼容多種輸入和輸出文件格式，有很好的分析輔助軟件vmd 有很好的維護服務 不需安裝 免費 缺點 萬一需要自己安裝的話比較麻煩/research/namd/ 主要針對生物體系，也適當兼容一般化學分子 優(yōu)點 有很好的內(nèi)置勢能模型 自定義新模型和新分子很方便 有很完善的維護網(wǎng)站 缺點 計算效率不高（收斂到16個處理器），運算速度慢 $400 主要針對生物體系，也包含部分化學體系 優(yōu)點 勢能模型更新很快 自定義新模型比較方便 維護服務很好 缺點 運算速度慢，計算效率低 $600http:/www.charmm

14、.org/ 一般性分子動力學軟件，對生物體系略有偏重 優(yōu)點 支持多種模型 免費 缺點 仍在開發(fā)中，某些方面還不完善/tinker/ 一般性分子模擬軟件 優(yōu)點 兼容當前大多數(shù)的勢能模型 編程水平高，計算效率高（比namd差，強于其他所有類似軟件） 可以模擬軟材料和固體物理系統(tǒng) 免費 缺點 維護差/sjplimp/lammps.html 一般性分子模擬軟件 優(yōu)點 界面友好 計算效率高（有兩個版本供選擇，適合于不同大小的體系） 維護服務很好 缺點 兼容性不好 100英鎊http:/www.cse.clrc.ac

15、.uk/msi/software/dl_poly/ 主要針對生物體系，也適當照顧一般化學體系 優(yōu)點 算法好，計算效率高 界面友好 維護服務好 免費軟件 缺點 兼容性不好/納米顆粒一般是人工合成物質(zhì)納米管，納米球，納米棒等等納米顆粒在生物醫(yī)學方面的應用(pnas, v99, p6645, 2002)藥物輸送，基因輸送 納米顆粒的生物安全性 (science, v300,p243, 2003)c60在魚類腦部的累計 (environmental health perspectives, v112, p1058, 2004)c60導致老鼠的肺部病癥 (tox

16、icol. sci. v77, p117, p126, 2004)用分子動力學研究c60和基因分子的作用md引擎 namd靜電力-pme恒溫 (300 k)衡壓 (1 bar)時間步長 = 2 fs離子-中和體系模擬時間 4-20 nsdna模型 amber 99 雙鏈, 12-20 堿基對 單鏈, 12-20 堿基c60, sp2 碳原子用lj 模型模擬溶劑 4600到6000水分子 tip3p模型水，離子水，離子水，離子水，離子三個方向施加周期邊界條件1 ns 的動畫起始結(jié)構 c60和dna能相互吸引，形成穩(wěn)定的雜化結(jié)構 二者的最初作用一般發(fā)生在1-2納秒之內(nèi)c60和b-dna之間的距離和

17、時間的函數(shù)c60和a-dna之間的距離與時間的函數(shù) c60和dna分子之間典型的結(jié)合能大約在15-40 kcal/mol, 取決于二者的相對作用位置和結(jié)合點 c60和dna分子之間的結(jié)合能遠大于熱波動kt (0.6 kcal/mol) 同時也遠大于在同等模擬環(huán)境下c60和c60之間的結(jié)合能 (7.5 kcal/mol)c60吸附于dna側(cè)面c60吸附于dna尾部 c60和b-dna形成穩(wěn)定的結(jié)合體，一般的結(jié)合能在20-40kcal/mol 比較穩(wěn)定的結(jié)合部位在dna的尾部或者側(cè)面（minor grooves）c60進入dna尾部 c60和a型dna的側(cè)面和尾部結(jié)合，而且，c60可以進入a型dn

18、a的尾部如果沒有c60吸附，尾部的堿基對形成氫鍵的原子對之間的距離 穿入dna尾部的c60破壞了dna尾部的堿基對的氫鍵c60進入dna尾部后，尾部的堿基對形成氫鍵的原子對之間的距離變形的角度 吸附在a型dna側(cè)面的dna使得dna堿基對之間的相對結(jié)構變形。本應平行的堿基對變形形成了一個大約20度的角 c60和單鏈dna之間的結(jié)合導致單鏈dna的嚴重變形 c60被單鏈dna的某些部分包圍，疏水的堿基平面和c60緊密接觸 c60和單鏈dna的結(jié)合能大約在20到30 kcal/mol c60和dna之間結(jié)合的驅(qū)動力是疏水作用dna親水部位的水密度分布c60和dna結(jié)合部位的水密度分布c60進入損傷部位去掉的dna片斷 c60很快進入dna的損傷部位，形成穩(wěn)定的結(jié)合體 c60和受損的dna的結(jié)合能大約在30 kcal/mol dna和c60能形成穩(wěn)定的結(jié)合體dna和c60的典型的結(jié)合能大約在20到40 kcal/mol，遠高于相同環(huán)境下c60和c60的結(jié)合能c60和dna的結(jié)合點在dna的疏水部位c60和dna的結(jié)合影響dna的結(jié)構c60的吸附使得dna的堿基對之間的角度變形c60進入a-dna的尾部，切斷了尾部堿基對的氫鍵c60能使單鏈dna嚴重變形c60進入dna的損傷部位，形成穩(wěn)定的結(jié)合體 allen & tildesley, computer simulation