此文為dddc_redsnow發(fā)表于biolover上的關(guān)于分子動力學的系列原創(chuàng)文章，相當經(jīng)典與精彩，特此將系列文章整合，一起轉(zhuǎn)載，望學習動力學的新手們共同學習，提高進步，在此特向dddc_redsnow本人表示感謝。

動力學系列之一（gromacs，重發(fā)）


在老何的鼓勵下，發(fā)一下我的gromacs上手手冊（我?guī)�人時用的，基本半天可以學會gromcas）
######################################################
# Process protein files step by step #
######################################################
pdb2gmx -f 2th_cap.pdb -o 2th_cap.gro -p 2th_cap.top -ignh -ter
nedit 2th_cap.top
editconf -f 2th_cap.gro -o 2th_cap_box.gro -d 1.5
genbox -cp 2th_cap_box.gro -cs -p 2th_cap.top -o 2th_cap_water.gro
make_ndx -f 2th_cap_water.gro -o 2th_cap.ndx
genpr -f 2th_cap_water.gro -n 2th_cap.ndx -o 2th_cap_All.itp
genpr -f 2th_cap_water.gro -n 2th_cap.ndx -o 2th_cap_M.itp
genpr -f 2th_cap_water.gro -n 2th_cap.ndx -o 2th_cap_C.itp
nedit Flavo.itp
grompp -f em.mdp -c 2th_cap_water.gro -p 2th_cap.top -o prepare.tpr
genion -s prepare.tpr -o 2th_cap_water_ion.gro -np 1 -pq 1
#####################################################
# Minimize step by step #
# 1. minimization fixing whole protein #
# 2. minimization fixing maincharin of protein #
# 3. minimization fixing Ca of protein #
# 4. minimization without fix #
#####################################################
grompp -np 4 -f em.mdp -c 2th_cap_water_ion.gro -p 2th_cap.top -o minimize_water.tpr
mpirun -np 4 mdrun -nice 0 -s minimize_water.tpr -o minimize_water.trr -c minimize_water.gro -e minimize_water.edr -g minimize_water.log &
grompp -np 4 -f em.mdp -c minimize_water.gro -p 2th_cap.top -o minimize_sidechain.tpr
mpirun -np 4 mdrun -nice 0 -s minimize_sidechain.tpr -o minimize_sidechain.trr -c minimize_sidechain.gro -e minimize_sidechain.edr -g minimize_sidechain.log &
grompp -np 4 -f em.mdp -c minimize_sidechain.gro -p 2th_cap.top -o minimize_sidechain_ex.tpr
mpirun -np 4 mdrun -nice 0 -s minimize_sidechain_ex.tpr -o minimize_sidechain_ex.trr -c minimize_sidechain_ex.gro -e minimize_sidechain_ex.edr minimize_sidechain_ex.log &
grompp -np 4 -f em.mdp -c minimize_sidechain_ex.gro -p 2th_cap.top -o minimize_all.tpr
mpirun -np 4 mdrun -nice 0 -s minimize_all.tpr -o minimize_all.trr -c minimize_all.gro -e minimize_allx.edr -g minimize_all.log&
editconf -f minimize_all.gro -o minimize_all.pdb
#####################################################
# Raise temperature step by step #
# 1. raise from 0K to 100K fixing whole protein #
# 2. raise from 100K to 200K fixing whole protein #
# 3. raise from 200K to 300K fixing whole protein #
# 4. balance fixing maincharin of protein #
# 5. balance fixing Ca of protein #
#####################################################
grompp -np 4 -f heat.mdp -c minimize_all.gro -p 2th_cap.top -o temperature100K.tpr
mpirun -np 4 mdrun -nice 0 -s temperature100K.tpr -o temperature100K.trr -c temperature100K.gro -e temperature100K.edr -g temperature100K.log &
grompp -np 4 -f heat.mdp -c temperature100K.gro -p 2th_cap.top -o temperature200K.tpr
mpirun -np 4 mdrun -nice 0 -s temperature200K.tpr -o temperature200K.trr -c temperature200K.gro -e temperature200K.edr -g temperature200K.log &
grompp -np 4 -f heat.mdp -c temperature200K.gro -p 2th_cap.top -o temperature300K.tpr
mpirun -np 4 mdrun -nice 0 -s temperature300K.tpr -o temperature300K.trr -c temperature300K.gro -e temperature300K.edr -g temperature300K.log &
g_energy -f temperature300K.edr -s temperature300K.tpr -o temperature300K.xvg
grompp -np 4 -f heat.mdp -c temperature300K.gro -p 2th_cap.top -o T300K_M.tpr
mpirun -np 4 mdrun -nice 0 -s T300K_M.tpr -o T300K_M.trr -c T300K_M.gro -e T300K_M.edr -g T300K_M.log &
grompp -np 4 -f heat.mdp -c T300K_M.gro -p 2th_cap.top -o T300K_Ca.tpr
mpirun -np 4 mdrun -nice 0 -s T300K_Ca.tpr -o T300K_Ca.trr -c T300K_Ca.gro -e T300K_Ca.edr -g T300K_Ca.log &
grompp -np 4 -f heat.mdp -c T300K_Ca.gro -p 2th_cap.top -o T300K_MD.tpr
mpirun -np 4 mdrun -nice 0 -s T300K_MD.tpr -o T300K_MD.trr -c T300K_MD.gro -e T300K_MD.edr -g T300K_MD.log &
至此，全部搞定，可以跑了。PJ的東西放出來敏感，原創(chuàng)的教學東西也可以收精華吧


動力學系列之二(amber基礎(chǔ)篇)


上次寫了一個，結(jié)果一下在因為發(fā)貼的原因丟了，加上最近忙的要命，不過還是吃完飯干活前，靜下心寫一下amber，amber遠比gromacs要復雜(不過還比不上charmm，估計第一次看charmm腳本的人一定覺得自己選錯了行當，類fortran的語言加上自創(chuàng)的格式，唉，又是一段心酸的往事...)。我把amber分成兩部分講，第一部分講基礎(chǔ)，第二部分講應用。這里先講第一部分。
1. 關(guān)于整體概括. amber的文件格式有三種最為重要：top, crd and pdb。pdb可以生成top and crd, top and crd也可以轉(zhuǎn)換成pdb，這些都是ascii碼文件，可以手動編輯(這也是DNA+Protein時候有時不得不作的事情，巨大的工作量)。top文件中是拓撲文件，相當于gromacs的top，不過直接把lib的參數(shù)搞過來，不要再調(diào)用lib了。crd文件是坐標和速度文件，類似于gro文件，pdb不用我說了吧。
2. amber中的概念: 所有的單元從原子，小分子，殘基，堿基等等一直到大分子都是unit，熟悉C++,JAVA的朋友想一下object就理解了，amber的xleap的操作就像一個外包體，利用object的獨立性進行關(guān)聯(lián)，組合或者分拆object, 實現(xiàn)不同的功能，每一個unit都是獨立單元。而標準的氨基酸，堿基就像class，你自己的蛋白就是object。每個class對應多個有自己名字的object。這些都是在xleap中實現(xiàn)的，xleap就是利用這種關(guān)系來導入，修改蛋白或者DNA/RNA。
3. amber自帶gaff力場（比gromacs那個網(wǎng)頁算的小分子正規(guī)多了），適用于大部分有機小分子，當然修改力場文件可以讓你獲得更多的支持，尤其是特殊的原子比如說鹵素，金屬原子。電荷是采用的resp擬和，而不是原子的partial charge，這個的優(yōu)點見JACS的原始文獻，不再贅述，主要是考慮極化。
4. 通過xleap搞定了大分子，小分子，水，溶劑離子，跑得時候就是sander了，當然8以后的pmemd是更好的選擇，不過只能跑pme方式的動力學。
5. 軌跡文件amber做的不好，太大，沒有壓縮。斷點續(xù)跑也不如gromacs，半個坐標的情況時有出現(xiàn)，只見原子叉叉滿天飛啊，那叫一個壯觀。分析工具首推ptraj（不是我推的，是他們，我個人喜歡carnal），加各種輔助工具可以實現(xiàn)gromacs的各種功能。
6. 優(yōu)點： 上次說過了，軌跡穩(wěn)定，重復性好，可信度高。有些朋友非要跑出個地動山搖，驚濤駭浪的軌跡，不推薦您使amber，估計等到它跑到那一天，您也畢業(yè)了。缺點：真的打算用amber嗎，打算跑10ns嗎？打算有什么特別有意思的想法嗎？恭喜您，您可以考慮買硬盤了。amber的存儲量大概是1ns->1.5GB, 如果想發(fā)一個好一點的文章，闡述一個精細的機理，沒有10條軌跡搞不定的，算算是多少硬盤....我的一個題平均120G軌跡。
洋洋灑灑寫了這么多，不過amber的真正強大之處還不在這里，它的強大在于它的應用廣泛，有著各種的“人無我有，人有我精”的套件，下次我在應用篇里面詳細講一下


動力學系列之三(amber應用篇)


上海的天氣真是一天一變啊，又是大雨滂沱了。反正出不去了，索性把amber的第二部分寫完，大家周末也有個消遣。
上回說amber的基礎(chǔ)應用談到了他有很多應用的套件，這里我把最常用的以及比較容易發(fā)好文章的幾個部分講解以下：
1。首先要說的就是用于藥物設(shè)計平衡的sander/pmemd，一般來說，用amber作一下平衡是同源模建的重要可選步驟之一，尤其適用于結(jié)構(gòu)不確定性大的loop和缺乏明確模板的地方，當然，如果您要是有macromodel，還是用那個，最近的文章發(fā)現(xiàn)申稿人很喜歡那個，可能是因為那個公司比較得到大家的認可吧
2。MMPBSA, 計算自由能的很好選擇，KUNZ他們組就喜歡拿這個說事，作為虛篩的終結(jié)器。當然了方法的想法是很好的，不過還是有一些缺點的，比如對共軛體系，可極化嚴重的分子計算結(jié)果很不準確。MMPBSA流行的時候只做一個這個舊可以發(fā)Protein Sci. 現(xiàn)在基本不可能了，不過還是很好的一個計算自由能的方法，平均半個小時一個分子
3。NMODE，有的朋友問要什么樣的分析工具，NMODE絕對是一個可選的東東，具體的原理我不講了，很簡單，矩陣二次求導做頻率，底頻高幅運動就是它分析的對象，Maccoman他們組只做這個就在JMB上發(fā)了無數(shù)的文章（人家是開山鼻祖，我們做，呵呵，估計就是BJ也很難�。�
4。Alan掃描：預測突變影響的工具，很粗糙，但是如果你計算資源不是很多的話，這個是推薦的方法。原理沒什么可講的，想用的去看manual，不到1頁
5。TMD，我最近做的東西，和什么SMD了，F(xiàn)MD了。。。。一大堆兄弟姐妹啊，不過也快做爛了，可以模擬蛋白大幅運動的過渡態(tài)
6。TI。目前公認的計算自由能最好方法，缺點是相當?shù)暮臅r間，一個師兄做n年，結(jié)果發(fā)現(xiàn)不適合他的體系，韶華逝水啊，沒有時間的朋友千萬不要碰
上面的這些方法都是最常用的，也是文章中暴光頻率最高，精通2到3個，一個好體系，綜合起來，可以沖擊JMB或者JACS。呵呵，外面雨終于小了，希望大家看了能有幫助，會一個技術(shù)不難，難的是把這個技術(shù)用的恰到好處，下次有時間聊聊怎么綜合使用各種工具分析問題


動力學系列之四（能量分析）


明天就是5.1，可是手里還是一大堆的活，實驗也是詭異的要命，時光真是快啊。前面對兩個基本軟件介紹了一下，本來想講講charmm和macromodel，后來想想，沒有什么意思，還是講一下分析方法對大家更適用，如果大家真的想聽這兩個軟件，后面我再補充。動力學分析方法多種多樣，但是能量永遠是不變的話題，很多問題的說明都是以能量為基礎(chǔ)的，這里說的能量包括各種力場能量，溶劑化能，自由能等。我主要對自由能講解一些如何進行計算。在amber中，可以通過mmpbsa來計算這個相對自由能，在免費的軟件gromacs中，能量的計算（不知道有心的朋友是否自己加和過他提供的能量項，是不是缺點兒什么啊^v^）相當?shù)拇植�，雖然提供了計算自由能的工具，但是源碼中卻是空，不知道到哪個版本才能實現(xiàn)。amber是比較好用，但是也是一個收費軟件，特別是想發(fā)文章的朋友要注意了。如果是不想買amber的朋友，一個中間路線就是用gromacs跑軌跡，用autodock算能量（自由能），雖然不如mmpbsa原理上準確，但是忽略部分熵效應也在變動不大的體系也還不錯。下面這個腳本是我自己寫的，可以實現(xiàn)這一功能，希望能對大家有個幫助： #!/bin/sh #'$1' is set for different list name #'$2' is set for the ligand name #prepare the $2.bnd in advance with command: #/~root/autodock/dist305/bin/pdbtoatm lig.pdbq| ~root/autodock/dist305/bin/atmtobnd >lig.bnd for snapshot in `cat $1` do rec=protein_mn_hoh1.${snapshot} lig=$2.${snapshot} echo ${rec} #assign the atomic solvation parameters to PDBQ-formatted version of your macromolecule, ~root/bin/addsol ${rec}.pdbq ${rec}.pdbqs ~root//bin/autotors -A $2.bnd ${lig}.pdbq ${lig}.out.pdbq > ${lig}.out.${rec}.epdb.com ~root/autodock/dist305/bin/autodock3 -p ${lig}.out.${rec}.dpf -l ${lig}.out.${rec}.epdb.log -c > $1_Binding ~root//bin/Docked.awk ${lig}.out.${rec}.epdb.log >> $1_Docked ~root//bin/Inter_score.awk ${lig}.out.${rec}.epdb.log >> $1_Inter_Energy ~root//bin/Intra_score.awk ${lig}.out.${rec}.epdb.log >> $1_Intra_Energy #rm -f *.map *.fld *.xyz *.err *.pdbqs *.glg *.log *.com *.gpf *.dpf Binding.awk #!~root/bin/gawk -f BEGIN{ name=ARGV[1] split(name,array,".") } { if($0 ~ /Estimated Free Energy of Binding/) {print array[1]" "$9" " $10} } Docked.awk #!~root/bin/gawk -f BEGIN{ name=ARGV[1] split(name,array,".") } { if($0 ~ /Final Docked Energy/) {print array[1]" "$7" " $8} } 返回小木蟲查看更多