以下是分子動力學模擬的一般性步驟, 具體的步驟和過程依賴於你研究的體系和所用的軟件, 但這並不影響我們把它當作一個入門指南.
1. 評估體系 首先需要對我們要進行模擬的體系做一個簡單的評估, 有三個問題是我們必須要明確的:
2. 選擇工具 選擇合適的模擬工具, 大前提是它能夠實現你所感興趣的目標. 這需要你非常廣泛謹慎地查閱文獻, 看看別人用這些工具都做了些什么, 有沒有和你的研究體系類似的, 相關的研究. 千萬不要做到一半才發現原來你用的工具根本就不能實現你所感興趣的idea, 切記!
在選擇合適的模擬工具時, 主要考慮下面兩點:
- 軟件的選擇. 通常與軟件主流使用的力場有關, 軟件本身也具有一定的偏向性.
- 蛋白體系: GROMACS, AMBER, NAMD均可
- DNA, RNA體系: 首選AMBER
- 界面體系: DL_POLY比較強大
- 材料體系:LAMMPS是不錯的選擇
- 力場的選擇. 力場用來描述體系中最小單元間的相互作用, 是對實驗性質或量子化學計算結果擬合后生成的經驗式. 有人會嫌它粗糙, 但它確確實實為我們模擬大系統提供了可能, 只能說關注的切入點不同罷了. 常見的有三類力場: 全原子力場, 聯合力場, 粗粒化力場. 此外還有所謂第一代, 第二代, 第三代力場的說法, 這里就不一一列舉了.
再次提醒注意: 必須選擇適合於我們所關注體系和所感興趣的性質及現象的力場.
3. 初始結構 通過實驗數據或者某些工具得到體系內的每一個分子的初始結構坐標文件. 之后, 我們需要按我們的想法把這些分子按照一定的規則或是隨機的放在一起, 從而得到整個體系的初始結構, 這也是我們模擬的輸入文件.
4. 輸入參數 得到了結構輸入文件, 我們還需要力場參數輸入文件, 也就是針對我們體系的力場文件. 這通常由所選用的力場決定, 包括電荷, 鍵合參數和非鍵參數等勢能函數的輸入參數.
5. 確定盒子 體系的大小通常由你所選用的盒子大小決定. 我們必須對可行性與合理性做出評估, 從而確定體系的大小, 這依賴於具體的體系.
6. 能量最小化 由於初始構象可能會存在兩個原子靠得太近的情況(稱之為bad contact), 所以需要在正式模擬開始的第一步對體系進行能量最小化. 比較常用的能量最小化方法有兩種, 最速下降法和共軛梯度法. 最速下降法是快速移除體系內應力的好方法, 但是接近能量極小點時收斂比較慢, 而共軛梯度法在能量極小點附近收斂效率高一些. 所以一般做能量最小化時都是先利用最速下降法進行優化, 完成之后再對得到的構象利用共軛梯度法優化一次, 這樣做能有效地保證后續模擬的進行.
7. 平衡模擬 你需要設置適當的模擬參數, 並且保證這些參數的設置與力場的構造過程相一致. 舉個簡單的例子, GROMOS力場是用范德華勢雙截斷來定義范德華參數的, 如果你用GROMOS力場的話也應該用雙截斷來處理范德華相互作用. 常見的模擬思路是, 先在NVT下限制住你的溶質(劑)做限制性模擬, 這是一個升溫的過程, 當溫度達到你設定的溫度后, 接着做NPT模擬, 此過程將調整體系的壓強進而使體系密度收斂.
如何判斷體系達到平衡是比較技術性的問題. 簡單地講可以通過以下幾種方式:
- 看能量(勢能, 動能和總能)是否收斂
- 看體系的壓強, 密度等等是否收斂
- 看體系的RMSD是否達到你能接受的范圍
- 其他經驗
8. 成品模擬 經過一段時間的平衡模擬, 在確定體系已經完全弛豫之后, 就可以開始采集數據了. 運行足夠長時間的模擬以確定我們所感興趣的現象或是性質能夠被觀測到, 並且務必確保此現象出現的可重復性.
9. 數據分析 數據拿到手后, 很容易通過一些可視化軟件得到軌跡動畫, 但這並不能拿來發文章. 真正的工作才剛剛開始——分析數據. 你所感興趣的現象或性質只是表面, 隱含在它們之中的機理才是文章的主題.
小分子過程：
reduce test.pdb > test_h.pdb
antechamber -i test_h.pdb -fi pdb -o test.mol2 -fo mol2 -c bcc -s 2
parmchk -i test.mol2 -f mol2 -o test.frcmod
#use tleap
tleap -f oldff/leaprc.ff99SB
source leaprc.gaff
AAA = loadmol2 test.mol2
list
loadamberparams test.frcmod
check AAA
saveoff AAA test.lib
saveamberparm AAA aaa.prmtop aaa.inpcrd
一個簡單模擬過程：
Ssh n1~n21 #進入節點
>tleap #xleap有圖形頁面，這個沒有，服務器上只能運行tleap
>source leaprc.ff14SB ##加載力場（官網上的那個protein應該去掉）
>foo=sequence { ACE ALA NME } ##存儲在foo單元中一個丙氨酸二肽
>desc foo ##可以查看該二肽
＞source leaprc.tip3p ##加載水盒
＞solvatebox foo TIP3PBOX 10.0 ##使用TIP3PBOX水模型，水盒表面距離蛋白10埃
>saveamberparm foo prmtop inpcrd ##保存參數拓撲文件和坐標文件
>quit ##退出tleap
一、能量最小化的腳本 (01_Min.in)
Minimize
&cntrl #起始符
imin=1, #進行最小化
ntx=1, #從inpcrd文件中讀取坐標（而不是速度）
irest=0, #不重新啟動模擬（不適用於最小化）
maxcyc=2000, #能量最小化的最大周期
ncyc=1000, # 0-ncyc 循環用最速下降算法，然后在 ncyc-maxcyc循環#轉換為共軛梯度算法
ntpr=100, #每ntpc個循環打印到amber輸出文件
ntwx=0, #沒有寫入amber軌跡文件（其實也可以寫入，=1）
cut=8.0, #非鍵截止距離
/
二、升溫的腳本 (02_Heat.in)
Heat
&cntrl
imin=0, ##不進行最小化
ntx=1, #無初始速度信息
irest=0, #不重啟run，run as a new MD
nstlim=10000, #跑MD的步數
dt=0.002, #時間步長（ps）
ntf=2, #不計算震動約束鍵的力
ntc=2, #允許震動來約束所有和氫鍵相連的鍵
tempi=0.0, #初始溫度
temp0=300.0, #最終溫度
ntpr=100, #每100步打印一次輸出
ntwx=100, #每ntwx步寫一次amber軌跡文件
ntwr=100, #每100步寫入一次rst文
cut=8.0, ntb=1, #定容的周期邊界 ntp=0, #沒有壓力控制 ntt=3, #Langevin恆溫器溫度控制
gamma_ln=2.0, #Langevin恆溫器碰撞頻率
nmropt=1, #核磁共振限制和讀取重量改變
ig=-1, #為偽隨機數生成程序隨機化種子
ntxo=1, #rst文件中的坐標，速度，盒子尺寸信息的格式，=1（ASC），=2（NetCDF）
ioutfm=1, #dcd文件中的坐標和速度的格式，=0（ASC），=1（NetCDF）
ntwv=0, #控制dcd文件中的速度信息是否寫入，=0（不寫入），=1（每隔ntwx步寫入一次速度）
iwrap=1, #處理鏡像
/
&wt type='TEMP0', istep1=0, istep2=9000, value1=0.0, value2=300.0 /
&wt type='TEMP0', istep1=9001, istep2=10000, value1=300.0, value2=300.0 /
&wt type='END' /
##最后三行允許恆溫器在整個模擬過程中改變目標溫度。對於前 #9000步，溫度將從0K增加到300K。對於步驟9001到10000，溫度將保持在300K。
##因為這個模擬非常簡短，所以NTPR和NTWX的設置非常低。使用這些設置進行長時間的MD模擬將產生非常大的輸出和軌跡文件，比普通的MD要慢。對於真正的跑MD，你需要增加NTPR和NTWX，不然時間太長。
三、跑動力學的腳本 （03_Prod.in）
Production
&cntrl
imin=0,
ntx=5, #從未格式化的inpcrd坐標文件中讀取坐標和速度
irest=1, #重啟之前的MD（意味着inpcrd文件提供初始原子速度）
nstlim=30000,
dt=0.002,
ntf=2,
ntc=2,
temp0=300.0, #恆溫器的溫度在300k
ntpr=100,
ntwx=100,
cut=8.0,
ntb=2, #使用常壓的周期性邊界條件
ntp=1, #對恆壓模擬使用Berendsen氣壓計
ntt=3,
gamma_ln=2.0,
ig=-1,
/
四、跑amber動力學sander（sander是amber的通用MD引擎）
現在我們已經有了所有的元素:參數和拓撲文件prmtop，坐標文件inpcrd，以及輸入文件01min。02 _heat。03 _prod。在此，我們准備運行實際的最小化、加熱和生產MD。
1、運行最小化
$ $AMBERHOME/bin/sander -O -i 01_Min.in -o 01_Min.out \ -p prmtop -c inpcrd -r 01_Min.rst -inf 01_Min.mdinfo
-O #如果輸出文件已經存在，就覆蓋他
-I 01_Min.in #選擇輸入文件
-o 01_Min.out #寫輸出文件
-p prmtop #選擇參數和拓撲文件
-c inpcrd #選擇坐標文件inpcrd
-r 01_Min.rst #用坐標和速度寫輸出重新啟動文件
-inf 01_Min.mdinfo #用模擬狀態寫MD信息文件
#####在sander完成之后，應該會有一個輸出文件01_Min.out，一個重啟文件01_Min.rst（將使用它繼續為系統加熱），和一個MD信息文件01_Min.mdinfo
Ps：提示prmtop不存在，我又做了一遍tleap
2、跑升溫動力學 （使用初始最小化的重新啟動文件，我們將對系統進行加熱）
$ $AMBERHOME/bin/sander -O -i 02_Heat.in -o 02_Heat.out \ -p prmtop -c 01_Min.rst -r 02_Heat.rst -x 02_Heat.mdcrd \ -inf 02_Heat.mdinfo
-c 01_Min.rst #現在對於輸入坐標我們從最小化中選擇重新啟動文件
-x 02_Heat.mdcrd #輸出動力學模擬的軌跡文件
#####在heating完成之后，我們將得到三個文件：02_Heat.out 02_Heat.rst 02_Heat.mdinfo ps：文件的內容同最小化
我們打開02_Heat.out，我們能夠看到系統信息，包括時間步長，能量和溫度。舉個例子（1000時間步長）：
NSTEP = 1000 TIME(PS) = 2.000 TEMP(K) = 30.20 PRESS = 0.0
Etot = -6945.8055 EKtot = 115.0391 EPtot = -7060.8446
BOND = 0.8441 ANGLE = 1.0099 DIHED = 9.0591
1-4 NB = 2.7105 1-4 EEL = 46.1439 VDWAALS = 1440.3323
EELEC = -8560.9443 EHBOND = 0.0000 RESTRAINT = 0.0000
Ewald error estimate: 0.6018E-03
NSTEP #所在的時間步長
TIME #模擬的總時間（包括重啟）
TEMP #系統溫度
PRESS #系統壓力
Etot #系統總能量
Ektot #系統的總動力學能量
EPtot #系統的總勢能
3、跑動力學
$ $AMBERHOME/bin/sander -O -i 03_Prod.in -o 03_Prod.out \ -p prmtop -c 02_Heat.rst -r 03_Prod.rst -x 03_Prod.mdcrd\ -inf 03_Prod.info &
### !!!注意：最后的&號表示讓程序在后台運行（因為模擬時間過長），期間可以通過$ tail -f 03_prod.out 查看跑的進程，也可以通過$ cat 03_Prod.info 查看精確的系統數據和估算的完成時間。 Ctrl+c 結束tail