AMBER分子動(dòng)力學(xué)簡(jiǎn)例

1、AMBER分子動(dòng)力學(xué)簡(jiǎn)例AMBER分子動(dòng)力學(xué)簡(jiǎn)例（一）概述以下是使用AMBER包的簡(jiǎn)單教程，希望對(duì)開始學(xué)習(xí)分子動(dòng)力學(xué)的同學(xué)有用處。申明一下，以下教程原版來自網(wǎng)上，是最最基本的教程，同時(shí)也非常實(shí)用，有非常好的借鑒意義。AMBER分子動(dòng)力學(xué)程序包是加州圣弗蘭西斯科大學(xué)（University of California San Francisco,UCSF)的Peter A Kollman和其同事編的，程序很全，現(xiàn)在已經(jīng)發(fā)展到版本9.0。AMBER功能涵蓋種類非常多的生物分子，包括蛋白、核算以及藥物小分子。軟件詳細(xì)情況請(qǐng)瀏覽.以下是AMBER軟件包中四個(gè)

2、主要的大程序：Leap：用于準(zhǔn)備分子系統(tǒng)坐標(biāo)和參數(shù)文件，有兩個(gè)程序：xleap：Xwindows版本的leap，帶GUI圖形界面。tleap：文本界面的Leap。Antechamber：用于生成少見小分子力學(xué)參數(shù)文件的。有的時(shí)候一些小分子Leap程序不認(rèn)識(shí)，需要加載其力學(xué)參數(shù)，這些力學(xué)參數(shù)文件就要antechamber生成。Sander：MD數(shù)據(jù)產(chǎn)生程序，即MD模擬程序，被稱做AMBER的大腦程序。Ptraj：MD模擬軌跡分析程序。學(xué)習(xí)項(xiàng)目本教程研究的題目是腦下垂體荷爾蒙之一的oxytocin，需要X光衍射晶體結(jié)構(gòu)文件1NPO.PDB。該文件包含了該荷爾蒙和其運(yùn)載蛋白的復(fù)合物，可以從蛋白數(shù)據(jù)庫

3、下載。PDB文件是不包含氫原子的，Leap程序會(huì)自動(dòng)的加上PDB文件缺少的東西。當(dāng)?shù)谝淮问褂肞DB文件的時(shí)候，要十分留意文件包含的信息，所以PDB文件缺少的殘疾、側(cè)鏈或者添加的變異都在這個(gè)地方記錄?？梢杂梦谋鹃喿x程序閱讀PDB文件頭部的信息，即以REMARKS開頭的信息文本行。PDB一個(gè)重要的信息是SSBOND記錄，該記錄說明結(jié)構(gòu)中二硫鍵的位置，這樣的信息在使用Leap程序建立分子結(jié)構(gòu)的時(shí)候需要。在本教程中，我們將比較oxytocin在真空和溶液中分子動(dòng)力學(xué)的差異，如果沒有二硫鍵，將影響整個(gè)結(jié)果。整個(gè)過程在Linux系統(tǒng)下完成，如對(duì)Linux不熟悉，請(qǐng)翻閱有關(guān)Linux書籍。建立項(xiàng)目目錄，并進(jìn)

4、入該目錄：mkdirmyprojcd  myproj下載1NPO.PDB文件到該目錄下，使用文本閱讀器閱讀該P(yáng)DB文件。從在文件的開頭部分可以得到該文件是一個(gè)二聚物的PDB文件，刪除文件中A、C、和D鏈對(duì)應(yīng)的文本行，剩下的就是oxytocin的晶體機(jī)構(gòu)了，保存為oxyt.pdb文件。使用Swiss PDB viewer分析以下oxyt.pdb文件，可以得知該pdb文件缺少了一個(gè)側(cè)鏈。如果使用Deep View，它會(huì)自動(dòng)給文件加上側(cè)鏈。重新保存pdb文件為oxyt.pdb文件。xleap和tleap程序功能是一樣的，都是準(zhǔn)備分子結(jié)構(gòu)的坐標(biāo)文件和拓?fù)湮募?。xleap啟動(dòng)一個(gè)X界面，比較慢

5、。tleap純文本，要快一些，我們選擇tleap：tleap-s-fleaprc.ff03其中l(wèi)eaprc.ff03是AMBER的2003力場(chǎng)文件。力場(chǎng)是一個(gè)很重要的文件，定義了分子、原子和殘疾等的信息，請(qǐng)查閱有關(guān)資料。oxy=loadpdb oxyt.pdb該命令讀入oxyt.pdb文件到oxy變量中。bondoxy.1.SG  oxy.6.SG連接oxy變量的第一和第六個(gè)殘疾的SG原子，即是連接二硫鍵?？梢允褂胏heck命令檢查oxy是否完好，沒有特殊情況的話，最后應(yīng)該是“OK”，表示分子系統(tǒng)狀態(tài)是好的，可以保存。check  oxy接下來就是保存分子系統(tǒng)了：savea

6、mberparmoxyoxy_vac.topoxy_vac.crd其中oxy_vac.top和oxy_vac.crd分別是分子系統(tǒng)真空狀態(tài)下的拓?fù)浜妥鴺?biāo)文件。接下來要給分子系統(tǒng)添加水環(huán)境，solvateoctoxyTIP3PBOX9.0 該命令讓Leap程序使用TIP3PBOX水模型,水環(huán)境距離分子為9納米。也可以用solvatebox命令，但是solvatebox添加的水環(huán)境是一個(gè)立方體，不是八面體。一般要求水表面到蛋白距離要達(dá)到8.5納米，避免MD過程中蛋白沖出水環(huán)境，但是水環(huán)境越大計(jì)算時(shí)間將越長(zhǎng)。如果MD過程內(nèi)含PME，那么水箱長(zhǎng)度要大于2倍截矩（cutoff），截矩在MD配置

7、文件中定義。加溶劑之后，要計(jì)算系統(tǒng)是否為電中性，使用命令：charge oxy如果現(xiàn)實(shí)不是電中性，要使用addions添加反性電荷，如Cl或者Na等。最后保存溶劑環(huán)境的系統(tǒng)拓?fù)浣Y(jié)構(gòu)和坐標(biāo)文件，并退出Leap程序：saveamberparmoxyoxy.topoxy.crdquit也可以把命令集中到一個(gè)，讓tleap讀取。如將上面的命令集中到以下文件中"oxy.leaprc"-sourceleaprc.ff03oxy=loadpdb oxyt.pdbbond  oxy.1.SGoxy.6.SGcheck  oxysaveamberparm oxy oxy_

8、vac.top oxy_vac.crdsolvateoct oxy saveamberparm oxy oxy.top  oxy.crdquit將以上兩橫線之間的命令保存為"oxy.leaprc"，然后使用以下命令一步搞定：tleap-s-foxy.leaprc到此，分子系統(tǒng)準(zhǔn)備已經(jīng)完成，目錄中會(huì)多了一個(gè)leap.log文件，這是Leap程序的記錄文件，它包含了leap程序所做的一切，包括自動(dòng)的和手動(dòng)命令?，F(xiàn)在已經(jīng)有了可以進(jìn)行分子動(dòng)力學(xué)的基本文件了，再編寫sander的控制文件，就可以進(jìn)行模擬了。這些將在后文中介紹。   &

9、#160;AMBER分子動(dòng)力學(xué)簡(jiǎn)例（二）分子動(dòng)力學(xué)（1）真空模式真空模式分子動(dòng)力學(xué)模擬將使用NVT系宗分兩步進(jìn)行，即系統(tǒng)能量最優(yōu)化和分子動(dòng)力學(xué)過程。1、系統(tǒng)能量最優(yōu)化。我們將使用淬火能量最優(yōu)化解除系統(tǒng)內(nèi)原子之間的不正常相互作用，這些原子之間的高能量相互作用如果不消除，可能影響后續(xù)的分子動(dòng)力學(xué)過程。因?yàn)閯?dòng)力學(xué)過程是能量梯度變化的，太高的能量壁壘可能讓MD局限在某一個(gè)能量局部最小化位置中。Sander程序是分子動(dòng)力學(xué)模擬程序，它的主要功能是能量最優(yōu)化，動(dòng)力學(xué)模擬和NMR優(yōu)化計(jì)算。我們必須給Sander一個(gè)運(yùn)行的配置文件，使其按照我們的要求進(jìn)行計(jì)算。能量最優(yōu)化的配置文件如下："min_va

10、c.in"oxytocin:  initial minimization prior to MD&cntrlimin = 1,maxcyc =  500,ncyc = 250,ntb = 0,igb = 0,cut =  12/-Sander程序的基本命令參數(shù)如下：sander  O i in o out p prmtop c inpcrd r restrt -ref refc x mdcrd v mdvel e  mden inf mdinfo所以我們的命令可以如下：sander O i  min_vac.in o

11、 min_vac.out p oxy_vac.top c oxy_vac.crd r oxy_vacmin.rst  &可以使用more命令或者tail命令查看min_vac.out的輸出內(nèi)容，那是能量最優(yōu)化的記錄。2、分子動(dòng)力學(xué)模擬。Sander程序的配置文件為：md_vac.inoxytocin  MD in-vacuo, 12 angstrom cut off, 250 ps&cntrlimin = 0, ntb =  0,igb = 0, ntpr = 100, ntwx = 500,ntt = 3, gamma_ln = 1.0,temp

12、i =  300.0, temp0 = 300.0,nstlim = 125000, dt = 0.002,cut =  12.0/-命令為：sander  O i md_vac.in o md_vac.out p oxy_vac.top c oxy_vacmin.rst r oxy_vacmd.rst  x oxy_vacmd.mdcrd ref oxy_vacmin.rst inf   &MD的計(jì)算過程一般比較久，真空相對(duì)與溶劑中要快以下，250ps的模擬大概在一個(gè)主頻為2。0GHz的Linux單機(jī)上運(yùn)行5分鐘。

13、以上配置文件中用到很多參數(shù)，這些參數(shù)這是sander程序參數(shù)的一小部分，以下將相應(yīng)解釋。如果要了解sander的其他參數(shù)，請(qǐng)閱讀AMBER用戶指南。&ctrl和"/"：sander的參數(shù)一般要求出現(xiàn)在這兩個(gè)標(biāo)識(shí)符號(hào)之間，參數(shù)以及這兩個(gè)標(biāo)識(shí)符被稱做控制模塊。cutoff：以納米為單位的截矩。即超出截矩范圍的非鍵連接相互作用將不計(jì)。ntr：原子位置能量抑制位，1表示抑制，0表示不抑制。imin：能量最優(yōu)化標(biāo)志位。1表示sander將進(jìn)行能量最優(yōu)化，0表示讓sander進(jìn)行分子動(dòng)力學(xué)模擬。macyc：能量最優(yōu)化次數(shù)。ncyc：便是經(jīng)過多少次能量?jī)?yōu)化以后，能量?jī)?yōu)化從淬火過程

14、變?yōu)樘荻茸兓^程。ntmin：能量?jī)?yōu)化方法標(biāo)志位。0表示前10個(gè)能量最優(yōu)化為淬火過程，然后進(jìn)行梯度能量?jī)?yōu)化；1表示ncyc次淬火過程，然后進(jìn)行梯度能量?jī)?yōu)化，為默認(rèn)值；2表示只進(jìn)行淬火過程。dx0：表示啟動(dòng)模擬步長(zhǎng)。dxm：最大優(yōu)化步數(shù)。drms：梯度能量?jī)?yōu)化標(biāo)準(zhǔn)，默認(rèn)值為1.0E-4  kcal/mol.A。更多參數(shù)將在后文中解釋。   AMBER分子動(dòng)力學(xué)簡(jiǎn)例（三）分子動(dòng)力學(xué)（2）  水環(huán)境中的分子動(dòng)力學(xué)模擬溶劑環(huán)境中的分子動(dòng)力學(xué)模擬分為以下四步進(jìn)行：1、溶劑環(huán)境能量最優(yōu)化。這一步保持溶質(zhì)（蛋白）不變，去除溶劑中能量不正常的范德華相

15、    互作用。2、整系統(tǒng)能量最優(yōu)化。去除整個(gè)系統(tǒng)中能量不正常的相互作用。3、有限制的分子動(dòng)力學(xué)。保持蛋白質(zhì)不動(dòng)，溶解溶劑的不同層，同時(shí)逐漸將系統(tǒng)溫度從0K提升到300K。4、整系統(tǒng)分子動(dòng)力學(xué)模擬。在一個(gè)大氣壓，300K的環(huán)境下整個(gè)系統(tǒng)分子動(dòng)力學(xué)模擬。可以得到成果的分子動(dòng)力學(xué)模擬。#1、溶劑環(huán)境能量最優(yōu)化。該步驟的配置文件min1.in如下：-oxytocin:  initial minimisation solvent + ions&cntrlimin = 1,maxcyc =  1000,ncyc = 500,ntb = 1,ntr

16、 = 1,cut = 10/Hold the  protein fixed500.0RES 1  9ENDEND該過程保持肽鏈不動(dòng)，其中500.0單位是kcal/mol，表示作用在肽鏈上使其不動(dòng)的力?！癛ES  1 9”表示肽鏈殘基數(shù)目，因?yàn)槲覀儗W(xué)習(xí)使用的oxytocin有9個(gè)殘基。模擬命令如下：sander O i min1.in  o min1.out p oxy.top c oxy.crd r oxy_min1.rst ref oxy.crd  &2、整系統(tǒng)能量最優(yōu)化。配置文件min2.in如下：-oxytocin: 

17、initial minimisation whole system&cntrlimin = 1,maxcyc =  2500,ncyc = 1000,ntb = 1,ntr = 0,cut =  10/-命令如下：sander  O i min2.in o min2.out p oxy.top c oxy_min1.rst r oxy_min2.rst  &3、有限制的分子動(dòng)力學(xué)。第 一步分子動(dòng)力學(xué)保持蛋白分子位置不變，但是不是完全固定每個(gè)原子，同時(shí)緩解蛋白分子周圍的水分子，是溶劑環(huán)境能量?jī)?yōu)化。在這個(gè)步驟中，我們將主要目的是對(duì) 特定的原子

18、使用作用力使其能量?jī)?yōu)化。我們要優(yōu)化溶劑環(huán)境，至少需要10ps，我們將使用20ps用來優(yōu)化我們上兩步制作的分子系統(tǒng)的周期性邊界的溶劑環(huán) 境。命令配置文件md1.in如下：-oxytocin:  20ps MD with res on protein&cntrlimin = 0,irest = 0,ntx =  1,ntb = 1,cut = 10,ntr = 1,ntc = 2,ntf = 2,tempi =  0.0,temp0 = 300.0,ntt = 3,gamma_ln = 1.0,nstlim = 10000, dt =  0.002,

19、ntpr = 100, ntwx = 500, ntwr = 1000/Keep protein fixed with  weak restraints10.0RES 1  9ENDEND-上述參數(shù)解釋如下：ntb  = 1：表示分子動(dòng)力學(xué)過程保持體積固定。imin = 0：表示模擬過程為分子動(dòng)力學(xué)，不是能量最優(yōu)化。nstlim =  #：表示計(jì)算的步數(shù)。dt = 0.002：表示步長(zhǎng)，單位為ps，0.002表示2fs。temp0 =  300：表示最后系統(tǒng)到達(dá)并保持的溫度，單位為K。tempi = 100：系統(tǒng)開始時(shí)的溫度。 gam

20、ma_ln =  1：表示當(dāng)ntt3時(shí)的碰撞頻率，單位為ps-1（請(qǐng)參考AMBER手冊(cè)）ntt = 3：溫度轉(zhuǎn)變控制，3表示使用蘭格氏動(dòng)力學(xué)。tautp =  0.1：熱浴時(shí)間常數(shù)，缺省為1.0。小的時(shí)間常數(shù)可以得到較好的耦聯(lián)。vlimit =  20.0：保持分子動(dòng)力學(xué)穩(wěn)定性速度極限。20.0為缺省值，當(dāng)動(dòng)力學(xué)模擬中原子速度大于極限值時(shí)，程序?qū)⑵渌俣冉档偷綐O限值以下。comp = 44.6：  溶劑可壓縮單位。ntc = 2：Shake算法使用標(biāo)志位。1表示不實(shí)用使用，2表示氫鍵將被計(jì)算，3表示所有鍵都將被計(jì)算在內(nèi)。tol =  #.#：坐標(biāo)

21、位置重新設(shè)置的幾何位置相對(duì)容忍度。我 們將使用一個(gè)較小的作用力，10kcal/mol。在分子動(dòng)力學(xué)中，當(dāng)ntr1時(shí)，作用力只需要5-10kcal/mol（我們需要引用一個(gè)坐標(biāo)文件做 分子動(dòng)力學(xué)過程的比較，我們需要使用"-ref"參數(shù)）。太大的作用力同時(shí)使用Shake算法和2fs步長(zhǎng)將使整個(gè)系統(tǒng)變得不穩(wěn)定，因?yàn)榇蟮淖饔昧κ瓜到y(tǒng) 中的原子產(chǎn)生大頻率的振動(dòng)，模擬過程并步需要。運(yùn)行命令如下：sander  O i md1.in o md1.out p oxy.top c oxy_min2.rst r oxy_md1.rst x  oxy_md1.mdcrd re

22、f oxy_min2.rst inf  &進(jìn)行MD的運(yùn)行時(shí)間一般較長(zhǎng)，可以使用程序的并行版本提交集群計(jì)算。在主頻位2.0GHz的P4單機(jī)上，大概需要一個(gè)鐘頭?？梢噪S時(shí)查看md1.in文件的程序輸出。4、整系統(tǒng)分子動(dòng)力學(xué)模擬。這一步中，我們將進(jìn)行整個(gè)系統(tǒng)的分子動(dòng)力學(xué)模擬，而不對(duì)某些特定原子位置進(jìn)行限制。因?yàn)橹R(shí)一個(gè)小例子，我們將只進(jìn)行250ps的MD計(jì)算。配置文件md2.in如下：-oxytocin:  250ps MD&cntrlimin = 0, irest = 1, ntx = 7,ntb = 2, pres0 = 1.0, 

23、ntp = 1,taup = 2.0,cut = 10, ntr = 0,ntc = 2, ntf = 2,tempi =  300.0, temp0 = 300.0,ntt = 3, gamma_ln = 1.0,nstlim = 125000, dt =  0.002,ntpr = 100, ntwx = 500, ntwr =  1000/-上面一些參數(shù)解釋如下:ntb=2：表示分子動(dòng)力學(xué)過程的壓力常數(shù)。ntp1：表示系統(tǒng)動(dòng)力學(xué)過程各向同性。taup  = 2.0：壓力緩解時(shí)間，單位為ps。pres1：引用1個(gè)單位的壓強(qiáng)。使用以下命令進(jìn)行MD：sa

24、nder O i md2.in  o md2.out p oxy.top c oxy_md1.rst r oxy_md2.rst x oxy_md2.mdcrd ref  oxy_md1.rst inf   &模擬時(shí)間較長(zhǎng)，在P4單機(jī)2.0GHz的上需要7.5個(gè)小時(shí)。到 此，模擬全部完成，接下來要對(duì)得到的數(shù)據(jù)進(jìn)行分析。主要數(shù)據(jù)文件.out文件，包含系統(tǒng)能量、溫度，壓力等等；.mdcrd文件，是分子動(dòng)力學(xué)軌跡文件， 可以求系統(tǒng)蛋白的RMSD，回轉(zhuǎn)半徑等等。數(shù)據(jù)分析根據(jù)不同的研究目的不同而不同，我們將在后文中進(jìn)行一些簡(jiǎn)單的分析。 &#

25、160; 使用pdbviewer軟件 load structure的時(shí)候，假設(shè)側(cè)鏈缺失會(huì) 有提示錯(cuò)誤信息。直接使用xleap產(chǎn)生該結(jié)構(gòu)的top  crd后，再使用ambpdb重新產(chǎn)生structure即可加上側(cè)鏈mber中生成小分子模板（轉(zhuǎn)）第一步：生成小分子模板蛋白質(zhì)中各氨基酸殘基的力參數(shù)是預(yù)先存在的，但是很多模擬過程會(huì)涉及配體分子，這些有機(jī)小分子有很高的多樣性，他們的力參數(shù)和靜電信息不可能預(yù)存在庫文件中，需要根據(jù)需要自己計(jì)算生成模板。amber中的antechamber 程序就是生成小分子模板的。生成模板要進(jìn)行量子化學(xué)計(jì)算，這一步可以由antechamber中附帶的mo

26、pac完成，也可以由gaussian完成，這里介紹用gaussian計(jì)算的過程。建議在計(jì)算前用sybyl軟件將小分子預(yù)先優(yōu)化，不要用gaussian優(yōu)化，大基組從頭計(jì)算進(jìn)行幾何優(yōu)化花費(fèi)的計(jì)算時(shí)間太長(zhǎng)。gaussian計(jì)算的輸入文件可以用antechamber程序直接生成，生成后去掉其中關(guān)于幾何優(yōu)化的參數(shù)即可將小分子優(yōu)化后的結(jié)構(gòu)存儲(chǔ)為mol2各式，上傳到工作目錄，用antechamber程序生成gaussian輸入文件，命令如下：antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat這樣可以生成49.in文件，下載到windows環(huán)境，運(yùn)行g(shù)auss

27、ian計(jì)算這個(gè)文件，如果發(fā)現(xiàn)計(jì)算時(shí)間過長(zhǎng)或者內(nèi)存不足計(jì)算中斷，可以修改文件選擇小一些的基組。獲得輸出文件49.out之后將它上傳到工作目錄，再用antechamber生成模板，命令如下：antechamber -i 49.out -fi gout -o 49mod.mol2 -fo mol2 -c resp運(yùn)行之后就會(huì)生成一個(gè)新的mol2文件，如果用看圖軟件打開這個(gè)文件會(huì)發(fā)現(xiàn)，原子的顏色很怪異，這是因?yàn)閙ol2的原子名稱不是標(biāo)準(zhǔn)的原子名稱，看圖軟件無法識(shí)別。下面一步是檢查參數(shù)，因?yàn)榭赡軙?huì)有一些特殊的參數(shù)在gaff中不存在需要程序注入，命令如下：parmchk -i 49mod.mol2 -f

28、mol2 -o 49mod這樣那些特殊的參數(shù)就存在49mod這個(gè)文件中了第二步：處理蛋白質(zhì)文件amber自帶的leap程序是處理蛋白質(zhì)文件的，他可以讀入PDB各式的蛋白質(zhì)文件，根據(jù)已有的力場(chǎng)模板為蛋白質(zhì)賦予鍵參數(shù)和靜電參數(shù)。PDB格式的文件有時(shí)會(huì)帶有氫原子和孤對(duì)電子的信息，但是在這種格式下氫原子和孤對(duì)電子的命名不是標(biāo)準(zhǔn)命名，力場(chǎng)模板無法識(shí)別這種不標(biāo)準(zhǔn)的命名，因此需要將兩者的信息刪除ATOM 12 1H ARG A 82 12.412 8.891 34.128 1.00 0.00 H在PDB各式里面，氫原子的信息會(huì)在第13或者14列出現(xiàn)H字符，可以應(yīng)用grep命令刪除在13或者14列出現(xiàn)H的行命

29、令如下：grep -v '.H' 1t4j.pdb > xgrep -v '.H' x > 1t4j_noh.pdb除了刪除氫和孤對(duì)電子，還應(yīng)該把文件中的結(jié)晶水、乙酸等分子刪除，這些分子的信息常常集中在文件的尾部，可以直接刪除。處理過之后的蛋白質(zhì)文件，只包括各氨基酸殘基和小分子配體的重原子信息，模擬需要的氫原子和水分子將在leap中添加接下來需要進(jìn)一步整理蛋白質(zhì)文件，主要關(guān)注氨基酸的不同存在形態(tài)和小分子的原子名稱。半胱氨酸有兩種存在形態(tài)，有的是兩個(gè)半胱氨酸通過二硫鍵相連，有的則是自由的，兩種半胱氨酸在力場(chǎng)文件中是用不同的unit來表示的，這相當(dāng)于是兩

30、個(gè)完全不同的氨基酸，需要手動(dòng)更改蛋白質(zhì)文件中半胱氨酸的名字，橋連的要用CYX，自由的用CYS,可以通過查看晶體的PDB文件來查看以二硫鍵存在的半胱氨酸殘基。組氨酸有若干種質(zhì)子態(tài)，和半胱氨酸一樣，也需要查閱文獻(xiàn)確定這些質(zhì)子態(tài)，并更改殘基名稱最后需要修改的是配體分子的原子名，這是工作量最大的事情，仔細(xì)觀察可以發(fā)現(xiàn)，一般PDB文件中配體的各個(gè)原子名稱，和我們上面通過antechamber 生成的49mod.mol2中原子名稱并不一致，這會(huì)造成leap無法識(shí)別這些原子，無法為這些原子賦予力參數(shù)和靜電參數(shù)，因此需要手動(dòng)更改蛋白質(zhì)文件中配體分子的原子名稱。進(jìn)行這一步可以同時(shí)用看圖軟件打開49mod.mol

31、2和蛋白質(zhì)文件，隱藏蛋白質(zhì)文件中除配體分子以外的所有結(jié)構(gòu)，旋轉(zhuǎn)兩個(gè)文件，讓他們姿態(tài)相近，以方便觀察，并且在各自均標(biāo)識(shí)原子名稱，然后用文字編輯軟件打開蛋白質(zhì)文件，核對(duì)看圖軟件中兩個(gè)分子對(duì)應(yīng)的原子名稱，手動(dòng)逐一修改。修改原子名稱需要極大的耐心和細(xì)心，如果發(fā)生錯(cuò)誤下一步的操作會(huì)無法繼續(xù)。我現(xiàn)在想到的也只有這個(gè)笨辦法，如果誰還有別的好法子，歡迎告訴我?，F(xiàn)在文件的準(zhǔn)備工作都已經(jīng)完成，該開始正式的模擬了第三步：生成拓?fù)湮募妥鴺?biāo)文件用amber進(jìn)行分子動(dòng)力學(xué)模擬需要坐標(biāo)和拓?fù)湮募?，坐?biāo)文件記錄了各個(gè)質(zhì)點(diǎn)所座落的坐標(biāo)，拓?fù)湮募涗浟苏麄€(gè)體系各質(zhì)點(diǎn)之間的鏈接狀況、力參數(shù)電荷等信息。這兩個(gè)文件是由leap 程序

32、生成的amber中有兩個(gè)leap程序，一個(gè)是純文字界面的tleap，一個(gè)是帶有圖形界面的Xleap。但是amber的圖形界面做得很差，用遠(yuǎn)程登錄使用圖形界面不僅麻煩而且慢，所以我比較偏愛使用tleap，兩個(gè)leap的命令是完全一樣的，其實(shí)用哪一個(gè)都無所謂。啟動(dòng)tleap，在shell里輸入命令tleap，leap就啟動(dòng)了，shell里會(huì)顯示-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/prep to search path.-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap

33、/lib to search path.-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/parm to search path.-I: Adding /usr/local/amber8/antechamber-1.23/dat/leap/cmd to search path.Welcome to LEaP!(no leaprc in search path)>這個(gè)>是leap的提示符下面要調(diào)入庫文件。amber是模擬生物分子的好手，主要就是依靠專門為蛋白質(zhì)多糖核酸量身訂做的amber力場(chǎng)，力場(chǎng)的所有參數(shù)都存儲(chǔ)在庫文件里，

34、所以打開leap第一件事便是調(diào)入庫文件。amber提供了很多種庫，這里我們只用到兩個(gè)庫，gaff和02庫，輸入命令：>source leaprc.gaff>source leaprc.ff02之后兩個(gè)庫就調(diào)入進(jìn)來了這時(shí)可以用list命令看看庫里都有什么：這里面羅列的就是庫里面的unit，包括20種氨基酸、糖以及核酸還有一些常見離子的參數(shù)下面一步是調(diào)入配體分子的模板，首先補(bǔ)全參數(shù)，輸入命令：>loadamberparams 49mod然后讀入模板文件,輸入命令：>MOL = loadmol2 49mod.mol2其中MOL是unit的名字，要保證這個(gè)名字和pdb文件中配體

35、的殘基名完全一致，否則系統(tǒng)仍然無法識(shí)別pdb文件中的小分子現(xiàn)在再輸入list命令會(huì)發(fā)現(xiàn)庫里面多了一個(gè)unit:那個(gè)就是配體分子的模板下面讀入pdb文件，輸入命令：>comp = loadpdb 1t4j_noh.pdb如果輸入這個(gè)命令之后，屏幕上出現(xiàn)了大量的創(chuàng)建unit或者atom的信息，如下所示，則說明上面一步的pdb文件處理沒有做好，還需要重新處理，通常這種情況都發(fā)生在配體分子上面，有時(shí)則是因?yàn)榈鞍踪|(zhì)中存在特殊殘基。Creating new UNIT for residue: FRJ sequence: 1Created a new atom named: O36 within re

36、sidue: .RCreated a new atom named: S33 within residue: .RCreated a new atom named: O35 within residue: .RCreated a new atom named: N34 within residue: .R如果屏幕僅僅顯示添加原子，這說明輸入的PDB文件中缺失了部分重原子，leap根據(jù)模板自動(dòng)補(bǔ)全了這些缺失的原子，這種情況不會(huì)影響進(jìn)一步的計(jì)算Added missing heavy atom: .R.AAdded missing heavy atom: .R.AAdded missing heav

37、y atom: .R.AAdded missing heavy atom: .R.A根據(jù)體系的具體情況，還可能要將成對(duì)的CYX殘基中的二硫鍵相連，有時(shí)候還會(huì)鏈接其他原子，比如將糖基鏈接在氨基酸殘基上，用bond命令可以完成，命令用法如下：>bond comp.35.SG comp.179.SG其中comp是剛才讀入的分子名稱，35和179是殘基序號(hào)，SG是CYX殘基模板中硫原子的名稱，用comp.35.SG這樣的語法就可以定位一個(gè)原子如果希望進(jìn)行考慮溶劑效應(yīng)的動(dòng)力學(xué)模擬，可能還需要為體系加上水，加水有很多種命令，這里只列舉一個(gè)：>solvatebox comp TIP3PBOX 1

38、0.0solvatebox命令是說要加上一個(gè)方形的周期水箱，comp指要加水的分子，TIP3PBOX是選擇的水模板名稱，10.0是水箱子的半徑有的體系總電荷不為0，為了模擬穩(wěn)定，需要加入抗衡離子，系統(tǒng)會(huì)自動(dòng)計(jì)算體系的靜電場(chǎng)分布，在合適的位置加上離子，命令如下：>addions comp Na+ 0意思是用鈉離子把體系總電荷補(bǔ)平，還可以選擇其他庫里面有的離子。完成這一步之后就可以輸出拓?fù)湮募妥鴺?biāo)文件了，但是為了方便起見，在運(yùn)行最后一步之前要先把leap里加工好的分子單獨(dú)存成一個(gè)庫文件，以后可以直接調(diào)用這個(gè)庫文件，免得重復(fù)上面的操作：>saveoff comp 1taj.off這樣就

39、生成了一個(gè)off文件在那里，下面輸出拓?fù)湮募妥鴺?biāo)文件>saveamberparm comp 1t4j.prmtop 1t4j.inpcrdChecking Unit.Building topology.Building atom parameters.Building bond parameters.Building angle parameters.Building proper torsion parameters.Building improper torsion parameters.total 1 improper torsion appliedBuilding H-Bond

40、 parameters.Not Marking per-residue atom chain types.Marking per-residue atom chain types.(Residues lacking connect0/connect1 -these don't have chain types marked:res total affectedCMET 1)(no restraints)>quit現(xiàn)在準(zhǔn)備好拓?fù)湮募妥鴺?biāo)文件，接下來就可以開始能量?jī)?yōu)化和動(dòng)力學(xué)模擬了。如果愿意的話還可以用ambpdb這個(gè)命令生成一個(gè)pdb文件，直觀地看一看生成了什么東西，命令如下：s

41、nowyowlslocalhost actualamber$ ambpdb -p 1t4j.prmtop <1t4j.inpcrd> kankan.pdb| New format PARM file being parsed.| Version = 1.000 Date = 09/08/06 Time = 16:36:09snowyowlslocalhost actualamber$可以下載之后用看圖軟件欣賞，如果加了溶劑盒子的話，看的時(shí)候可要小心，會(huì)恨嚇人的樣子。第四步：能量?jī)?yōu)化用amber進(jìn)行分子動(dòng)力學(xué)模擬需要坐標(biāo)和拓?fù)湮募?，這在上一步已經(jīng)完成了，分別生成了1t4j.prmto

42、p 和1t4j.inpcrd兩個(gè)文件，下面一步就是動(dòng)力學(xué)模擬之前的能量?jī)?yōu)化了。由于我們進(jìn)行的起始結(jié)構(gòu)來自晶體結(jié)構(gòu)或者同源模建，所以在分子內(nèi)部存在著一定的張力，能量?jī)?yōu)化就是在真正的動(dòng)力學(xué)之前，釋放這些張力，如果沒有這個(gè)步驟，在動(dòng)力學(xué)模擬開始之后，整個(gè)分子可能會(huì)因此散架。能量?jī)?yōu)化由sander模塊完成，運(yùn)行sander至少需要三個(gè)輸入文件，分別是分子的拓?fù)湮募鴺?biāo)文件以及sander的控制文件?，F(xiàn)在分子的拓?fù)湮募妥鴺?biāo)文件已經(jīng)完成，需要建立輸入文件，min_1.inInitial minimisation of our structures&cntrl   imin=

43、1, maxcyc=4000, ncyc=2000,   cut=10, ntb=1,ntr=1,   restraint_wt=0.5   restraintmask=':1-283'/文件首行是說明，說明這項(xiàng)任務(wù)的基本情況；   &cntrl與/之間的部分是模擬的參數(shù)其中imin=1表示任務(wù)是能量?jī)?yōu)化，maxcyc=4000表示能量?jī)?yōu)化共進(jìn)行4000步，ncyc=2000表示在整個(gè)能量?jī)?yōu)化的4000步中，前 2000步采用最陡下降法，在2000步之后轉(zhuǎn)換為共軛梯度法，如果模擬的時(shí)候不希望

44、進(jìn)行方法的轉(zhuǎn)換，可以再加入另一個(gè)關(guān)鍵詞NTMIN，如果NTMIN =0則全程使用共軛梯度法，NTMIN=2則全程使用最陡下降法，此外還有=3和=4的選項(xiàng)，分別是xmin法和lmod法，具體情況可以看手冊(cè)。第二行的cut=10表示非鍵相互作用的截?cái)嘀担瑔挝皇前＃?ntb=1表示使用周期邊界條件，這個(gè)選項(xiàng)要和前面生成的拓?fù)湮募鴺?biāo)文件相匹配，如果前面加溶劑時(shí)候用的是盒子水，就設(shè)置ntb=1，如果加的是層水，那就應(yīng)該選擇ntb=0；ntr=1表示在能量?jī)?yōu)化的過程中要約束一些原子，約束的是哪些原子呢？后面有。第三行和第四行都是關(guān)于約束原子的信息，restraint_wt=0.5限定了約束的力常數(shù)，在這里約束原子就是把原子用一根彈簧拉在固定的位置上，一旦原子偏離固定的位置，系統(tǒng)就會(huì)給他施加一個(gè)回復(fù)力，偏離的越遠(yuǎn)，回復(fù)力越大，回復(fù)力就是由這個(gè)力常數(shù)決定的，單位是Kcal/(mol*A)。 restraintmask=':1-283'表示約束的是從1到283號(hào)殘基，在這個(gè)分子中，1-283號(hào)殘基是蛋白質(zhì)上的氨基酸殘基，從284號(hào)開始，就都是水了，所以這個(gè)命令的意思就是，約束整個(gè)蛋白質(zhì)，自由地優(yōu)化溶劑分子。因?yàn)槿軇┓肿邮乔懊娴膖leap自動(dòng)給加上的，所以一定要額外多關(guān)注一些。下面運(yùn)行sander來執(zhí)行這個(gè)優(yōu)化：snowyowlslocalhost actualamb