| 5 | 1/1 | 返回列表 |
| 查看: 11001 | 回復(fù): 16 | ||||||||||||||||
| 當(dāng)前只顯示滿足指定條件的回帖,點(diǎn)擊這里查看本話題的所有回帖 | ||||||||||||||||
御劍江湖榮譽(yù)版主 (著名寫手)
天池冶海
|
[交流]
【轉(zhuǎn)帖】用Amber+Gaussian做小分子力場(chǎng) 已有14人參與
|
|||||||||||||||
|
用Amber+Gaussian做小分子力場(chǎng),不是我的研究領(lǐng)域,不過閱讀此篇文章,發(fā)現(xiàn)不錯(cuò),轉(zhuǎn)載給此研究領(lǐng)域的蟲子們。。。 第一步:生成小分子模板 蛋白質(zhì)中各氨基酸殘基的力參數(shù)是預(yù)先存在的,但是很多模擬過程會(huì)涉及配體分子,這些有機(jī)小分子有很高的多樣性,他們的力參數(shù)和靜電信息不可能預(yù)存在庫文件中,需要根據(jù)需要自己計(jì)算生成模板。amber中的antechamber 程序就是生成小分子模板的。生成模板要進(jìn)行量子化學(xué)計(jì)算,這一步可以由antechamber中附帶的mopac完成,也可以由gaussian完成,這里介紹用gaussian計(jì)算的過程。 建議在計(jì)算前用sybyl軟件將小分子預(yù)先優(yōu)化,不要用gaussian優(yōu)化,大基組從頭計(jì)算進(jìn)行幾何優(yōu)化花費(fèi)的計(jì)算時(shí)間太長。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ù)aussian計(jì)算這個(gè)文件,如果發(fā)現(xiàn)計(jì)算時(shí)間過長或者內(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 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的行命令如下: grep -v ‘^………….H’ 1t4j.pdb > x grep -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)于是兩個(gè)完全不同的氨基酸,需要手動(dòng)更改蛋白質(zhì)文件中半胱氨酸的名字,橋連的要用CYX,自由的用CYS 組氨酸有若干種質(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.mol2和蛋白質(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è)笨辦法,如果誰還有別的好法子,歡迎告訴我。 現(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 程序生成的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/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ǔ)在庫文件里,所以打開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文件中配體的殘基名完全一致,否則系統(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: 1 Created a new atom named: O36 within residue: .R Created a new atom named: S33 within residue: .R Created a new atom named: O35 within residue: .R Created a new atom named: N34 within residue: .R 如果屏幕僅僅顯示添加原子,這說明輸入的PDB文件中缺失了部分重原子,leap根據(jù)模板自動(dòng)補(bǔ)全了這些缺失的原子,這種情況 不會(huì)影響進(jìn)一步的計(jì)算 Added missing heavy atom: .R.A Added missing heavy atom: .R.A Added missing heavy atom: .R.A Added 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 10.0 solvatebox命令是說要加上一個(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 這樣就生成了一個(gè)off文件在那里,下面輸出拓?fù)湮募妥鴺?biāo)文件 >saveamberparm comp 1t4j.prmtop 1t4j.inpcrd Checking Unit. Building topology. Building atom parameters. Building bond parameters. Building angle parameters. Building proper torsion parameters. Building improper torsion parameters. total 1 improper torsion applied Building H-Bond 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 affected CMET 1 ) (no restraints) >quit 現(xiàn)在準(zhǔn)備好拓?fù)湮募妥鴺?biāo)文件,接下來就可以開始能量優(yōu)化和動(dòng)力學(xué)模擬了。如果愿意的話還可以用ambpdb這個(gè)命令生成一 個(gè)pdb文件,直觀地看一看生成了什么東西,命令如下: [snowyowls@localhost 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:09 [snowyowls@localhost actualamber]$ 可以下載之后用看圖軟件欣賞,如果加了溶劑盒子的話,看的時(shí)候可要小心,會(huì)恨嚇人的樣子。 第四步:能量優(yōu)化 用amber進(jìn)行分子動(dòng)力學(xué)模擬需要坐標(biāo)和拓?fù)湮募,這在上一步已經(jīng)完成了,分別生成了1t4j.prmtop 和1t4j.inpcrd兩個(gè)文件 ,下面一步就是動(dòng)力學(xué)模擬之前的能量優(yōu)化了。 由于我們進(jìn)行的起始結(jié)構(gòu)來自晶體結(jié)構(gòu)或者同源模建,所以在分子內(nèi)部存在著一定的張力,能量優(yōu)化就是在真正的動(dòng)力學(xué)之前, 釋放這些張力,如果沒有這個(gè)步驟,在動(dòng)力學(xué)模擬開始之后,整個(gè)分子可能會(huì)因此散架。 能量優(yōu)化由sander模塊完成,運(yùn)行sander至少需要三個(gè)輸入文件,分別是分子的拓?fù)湮募,坐?biāo)文件以及sander的控制文件。 現(xiàn)在分子的拓?fù)湮募妥鴺?biāo)文件已經(jīng)完成,需要建立輸入文件,min_1.in Initial minimisation of our structures &cntrl imin=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ù)是能量優(yōu)化,maxcyc=4000表示能量優(yōu)化共進(jìn)行4000步,ncyc=2000表示在整個(gè)能量優(yōu)化的4000步中, 前 2000步采用最陡下降法,在2000步之后轉(zhuǎn)換為共軛梯度法,如果模擬的時(shí)候不希望進(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表示在能量優(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)化: [snowyowls@localhost actualamber]$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.i npcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out 命令中,-O表示覆蓋所有同名文件,-i min_1.in表示sander的控制文件是min_1.in,-p 1t4j.prmtop表示分子的拓?fù)湮募?br /> ,-c 1t4j.inpcrd表示坐標(biāo)文件,-ref 1t4j.inpcrd是參考坐標(biāo)文件,只有在控制文件中出現(xiàn)關(guān)鍵詞ntr=1的時(shí)候才需要給 sander指定-ref文件,這是約束原子的參考坐標(biāo),- ref 1t4j.inpcrd就是說以1t4j.inpcrd中的坐標(biāo)為準(zhǔn)進(jìn)行約束原子的優(yōu)化。 以上這四個(gè)是輸入文件。-r 1t4j_min1.rst表示經(jīng)過模擬之后新的原子坐標(biāo)會(huì)輸出到1t4j_min1.rst文件中,-o 1t4j_min1.out 則表示優(yōu)化過程中的相關(guān)信息都會(huì)寫入到1t4j_min1.out文件中。 運(yùn)行起這個(gè)命令之后,等拿到 1t4j_min1.rst文件就意味著對(duì)溶劑的優(yōu)化已經(jīng)差不多了,顯然下面還需要對(duì)蛋白質(zhì)本身進(jìn)行優(yōu)化, 這個(gè)優(yōu)化還要分兩步進(jìn)行,控制文件分別是min_2.in 和min_3.in min_2.in Initial minimisation of our structures &cntrl imin=1, maxcyc=5000, ncyc=2500, cut=10, ntb=1,ntr=1, restraint_wt=0.5 restraintmask=’:1-283@CA,N,C’ / 在這里發(fā)生變化的是約束原子的范圍, ‘:1-283@CA,N,C’表示1-283號(hào)殘基中名叫CA,N和C的原子,這些原子實(shí)際上是蛋白質(zhì)主鏈上 的原子,也就是說這一次的優(yōu)化是約束了蛋白質(zhì)主鏈上的原子之后,對(duì)溶劑和側(cè)鏈原子進(jìn)行自由優(yōu)化。 min_3.in Initial minimisation of our structures &cntrl imin=1, maxcyc=10000, ncyc=5000, cut=10, ntb=1, / 在這里不再進(jìn)行約束原子的優(yōu)化了,對(duì)整個(gè)分子進(jìn)行全原子優(yōu)化。 三步優(yōu)化的命令如下: [snowyowls@localhost actualamber]$ sander -O -i min_1.in -p 1t4j.prmtop -c 1t4j.inpcrd -ref 1t4j.inpcrd -r 1t4j_min1.rst -o 1t4j_min1.out [snowyowls@localhost actualamber]$ sander -O -i min_2.in -p 1t4j.prmtop -c 1t4j_min1.rst -ref 1t4j_min1.rst -r 1t4j_min2.rst -o 1t4j_min2.out [snowyowls@localhost actualamber]$ sander -O -i min_3.in -p 1t4j.prmtop -c 1t4j_min2.rst -r 1t4j_heat0.rst -o 1t4j_min3.out 接下來就是真正的分子動(dòng)力學(xué)模擬了…… |
物理前沿 | xuexijisuan | 分子模擬要術(shù) | amber |
分子模擬 | 精華帖轉(zhuǎn)載 | 分子對(duì)接 | Amber |
計(jì)算化學(xué) | 搜集 | 申請(qǐng) | 精華 |
分子動(dòng)力學(xué)模擬 |

木蟲 (著名寫手)

銅蟲 (正式寫手)
|
最后需要修改的是配體分子的原子名,這是工作量最大的事情,仔細(xì)觀察可以發(fā)現(xiàn),一般PDB文件中配體的各個(gè)原子名稱,和我們上面通過antechamber 生成的49mod.mol2中原子名稱并不一致,這會(huì)造成leap無法識(shí)別這些原子,無法為這些原子賦予力參數(shù)和靜電參數(shù), 因此需要手動(dòng)更改蛋白質(zhì)文件中配體分子的原子名稱。 進(jìn)行這一步可以同時(shí)用看圖軟件打開49mod.mol2和蛋白質(zhì)文件,隱藏蛋白質(zhì)文件中除配體分子以外的所有結(jié)構(gòu),旋轉(zhuǎn)兩個(gè)文件,讓他們姿態(tài)相近,以方便觀察,并且在各自均標(biāo)識(shí)原子名稱,然后用文字編輯軟件打開蛋白質(zhì)文件,核對(duì)看圖軟件中兩個(gè)分子對(duì)應(yīng)的原子名稱,手動(dòng)逐一修改。 這個(gè)是用antechamber做小分子力場(chǎng)最為困難的地方,相當(dāng)?shù)牟蝗诵曰? |

| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 本科新能源科學(xué)與工程,一志愿華理能動(dòng)285求調(diào)劑 +3 | AZMK 2026-03-28 | 6/300 |
|
|---|---|---|---|---|
|
[考研] 070305高分子化學(xué)與物理 304分求調(diào)劑 +4 | c297914 2026-03-28 | 4/200 |
|
|
[考研] 求調(diào)劑 +3 | QiMing7 2026-03-25 | 4/200 |
|
|
[考研] 085701環(huán)境工程,267求調(diào)劑 +16 | minht 2026-03-26 | 16/800 |
|
|
[考研] 0856,材料與化工321分求調(diào)劑 +12 | 大饞小子 2026-03-27 | 13/650 |
|
|
[考研] 291求調(diào)劑 +6 | HanBeiNingZC 2026-03-24 | 6/300 |
|
|
[考研] 265求調(diào)劑11408 +3 | 劉小鹿lu 2026-03-27 | 3/150 |
|
|
[考研] 求調(diào)劑 +4 | 零八# 2026-03-27 | 4/200 |
|
|
[考研] 312求調(diào)劑 +9 | 上岸吧ZJY 2026-03-22 | 13/650 |
|
|
[考研] 一志愿211 初試270分 求調(diào)劑 +6 | 谷雨上岸 2026-03-23 | 7/350 |
|
|
[考研]
|
平樂樂樂 2026-03-26 | 4/200 |
|
|
[考研] 材料考研求調(diào)劑 +3 | Dendel 2026-03-23 | 6/300 |
|
|
[考研] 資源與環(huán)境 調(diào)劑申請(qǐng)(333分) +9 | holy J 2026-03-21 | 9/450 |
|
|
[考研] 263求調(diào)劑 +6 | yqdszhdap- 2026-03-22 | 10/500 |
|
|
[有機(jī)交流]
20+3
|
FENGSHUJEI 2026-03-23 | 5/250 |
|
|
[考研] 300求調(diào)劑,材料科學(xué)英一數(shù)二 +5 | leaflight 2026-03-24 | 5/250 |
|
|
[考研] 284求調(diào)劑 +3 | yanzhixue111 2026-03-23 | 6/300 |
|
|
[考研] 269求調(diào)劑 +4 | 我想讀研11 2026-03-23 | 4/200 |
|
|
[考研] 一志愿重慶大學(xué)085700資源與環(huán)境,總分308求調(diào)劑 +7 | 墨墨漠 2026-03-23 | 8/400 |
|
|
[考研] 求調(diào)劑 +5 | Zhangbod 2026-03-21 | 7/350 |
|