文章來源：“分子動力學(xué)”公眾號

鏈接：https://mp.weixin.qq.com/s/nQ0MdLmIXHPB8VLSe7CTZQ

一、進行模擬的基本步驟

以下是進行分子動力學(xué)模擬的最基本步驟，來自Gromacs維基。具體的步驟可能會因為模擬目標的不同而存在差異，這里只做基本的概述。

十分清楚要什么樣的模擬系統(tǒng)性質(zhì)和現(xiàn)象，明白模擬的目標。

選擇適當?shù)墓ぞ咭詫崿F(xiàn)模擬目標。這需要熟知其他研究人員對相似模擬系統(tǒng)的模擬方法等資料，要好好讀別人的論文，主要包括：

• 模擬使用的軟件，當然還要考慮使用的力場，有時因為力場關(guān)系，需要更改使用的軟件。

• 力場描述了系統(tǒng)中原子的性質(zhì)和相互作用，選擇一個適合自己研究目標的力場至關(guān)重要，還是需要多看論文。根據(jù)sensenbobo血淚史和別人的經(jīng)驗，模擬核酸，可以采用amber力場，水分子可以使用TIP3P；模擬蛋白質(zhì)，可以采用OPLS-AA，水分子使用TIP4P（重申一下，我本人沒有跑過測試）。

• 獲取模擬分子的坐標文件（PDB），或者自己生成一個。自己生成很好啊，強烈鼓勵，但是要十分明白有一定的危險性，要明白每一步模擬系統(tǒng)的現(xiàn)象是為什么。

• 根據(jù)模擬軟件，生成模擬系統(tǒng)初步坐標文件。

• 獲取或者生成系統(tǒng)的拓撲文件，比如Gromacs的pdb2gmx命令，或者使用網(wǎng)絡(luò)服務(wù)器PRODRG（google一下這幾個字母），amber的leap命令和antechamber，NAMD的psfgen或者VMD。請參考相關(guān)手冊。

• 給模擬系統(tǒng)添加一個盒子（盒子的翻譯方法我十分不喜歡，箱子顯得太大，邊界顯得太難理解，但是盒子聽起來想到蛋糕，我們是嚴肅的科學(xué)家。娃哈哈。）在盒子中添加水分子，添加離子如NA和CL之類（添油加醋，生活就是這樣進行的）。添加任何東西到你的系統(tǒng)中，就要相應(yīng)的修改一下拓撲文件，以保持與系統(tǒng)對應(yīng)，不然會出錯哦。

• 進行能量最優(yōu)化。這樣做是消去系統(tǒng)中原子之間距離太近造成的高能量。必要的時候，分子可以現(xiàn)在真空中能量最優(yōu)話，然后再做溶劑環(huán)境中能量最優(yōu)話。

• 選擇恰當?shù)哪M參數(shù)模擬平衡模擬系統(tǒng)。需要慎重考慮力場的作用，必要的時候保持系統(tǒng)體積溫度不變（NVT系綜）進行位置限制性平衡系統(tǒng)，然后在保持系統(tǒng)溫度壓強不變（NPT系綜）修正系統(tǒng)的密度。然后再選擇正確的系綜（NPT，NVT，NET，BT？？）進行數(shù)據(jù)收集。（真煩！）

• 跑足夠長時間的模擬，使系統(tǒng)達到理想的平衡。使用模擬軟件的工具或者其他軟件觀察系統(tǒng)的性質(zhì)，如溫度，能量瞪瞪。使用適當?shù)膮?shù)進行數(shù)據(jù)采集模擬，注意模擬銜接過程中各個系統(tǒng)性質(zhì)保持不變，如果原子速度之類。這里還要考慮力場的作用，同時要考慮怎么測量描述模擬目標。

• 模擬足夠長的時間使模擬目標明顯出現(xiàn)（可能也出現(xiàn)跟你期望相反的結(jié)果，出現(xiàn)率應(yīng)該百分之五十，加上個人感情因素，可能會達到百分之八十。保重?。?。

• 分析模擬結(jié)果，解釋得到結(jié)果。要好好利用視圖工具，同時要十分注意漂亮的三維圖說明不了什么問題。圖表，比較，統(tǒng)計，一個都不能少。就這樣吧。

二、Gromacs使用－新手入門篇

其實gromacs的思路比較簡單，大體就是以下的流程,在此先對em進行一下總結(jié)：

2.1 convert the pdb-file to a gromacs structure file(.gro / .pdb) and a gromacs topology file(.top)

這一步需要gromacs的預(yù)處理程序pdb2gmx 標準的輸入命令應(yīng)該是：

pdb2gmx -ignh -f speptide.pdb -p speptide.top -o speptide.gro

這個命令之后有一個force field的選擇 系統(tǒng)給出了11種力場，具體以后實際計算中如何選擇力場還是應(yīng)該多參閱文獻和gromacs說明書的。這步之后可以在屏幕看到：

--------- PLEASE NOTE ------------ You have succesfully generated a topology from: speptide.pdb. The G43a1 force field and the spc water model are used. Note that the default mechanism for selecting a force fields has changed, starting from GROMACS version 3.2.0 --------- ETON ESAELP ------------

至此，完成了從pdb到gro和top的轉(zhuǎn)變。另外，我也看過有的作這一步命令時候并不進行top文件的轉(zhuǎn)換，因為在genbox（對盒子加水溶劑化時還要對top文件進行修改）

2.2 solvate the peptide in water

這一步需要用到editconf和genbox兩個程序，命令的標準書寫應(yīng)該是：

editconf -bt cubic -f speptide.gro -o speptide_box.gro -c -d 0.5

可以在屏幕上看到：

Read 191 atoms Volume: 8.17315 nm^3, corresponds to roughly 3600 electrons No velocities found system size : 1.774 3.372 1.367 (nm) diameter : 3.517 (nm) center : 2.650 1.453 2.417 (nm) box vectors : 1.774 3.372 1.366 (nm) box angles : 90.00 90.00 90.00 (degrees) box volume : 8.17 (nm^3) shift : -0.391 0.806 -0.158 (nm) new center : 2.259 2.259 2.259 (nm) new box vectors : 4.517 4.517 4.517 (nm) new box angles : 90.00 90.00 90.00 (degrees) new box volume : 92.18 (nm^3)

因為程序的默認，所以這一步命令可以寫作：

editconf –f speptide –o –d 0.5

這里－f之后的sp文件是默認gro文件或者pdb文件，－o后面生成的output file是放入盒子的體系的gro文件 按照這樣的格式，輸出的文件被默認命名為－out.gro

體系放入盒子之后在盒子中注入水

genbox -cp speptide_box.gro -cs -p speptide.top -o speptide_water.gro

也可簡寫為 ：

genbox –cp out –cs –p sprptide –o b4em Screen-- Output configuration contains 9035 atoms in 2967 residues Volume : 92.1813 (nm^3) Density : 994.449 (g/l) Number of SOL molecules: 2948

之后看加水之后的top文件

用命令

tail speptide.top

可以看到：

[ system ] ; Name Protein in water [ molecules ] ; Compound #mols Protein 1 SOL 2948

其中SOL后面的數(shù)字2948就是genbox加入盒子的水分子的數(shù)目。

Editconf程序的另一個用途是講gro文件轉(zhuǎn)化回pdb這時可以講speptide_water.gro轉(zhuǎn)化回pdb觀察

editconf –f speptide_water.gro –o speptide.pdb

拖回本機 用spbdv或者vmd觀察 。

下一步是生成索引文件，需使用 make _ ndx 來生成一個。

make_ndx –f b4em

2.3 perform an enery minimization of the peptide in solvent

現(xiàn)在模擬系統(tǒng)已經(jīng)差不多準備好了。在我們開始動力學(xué)之前，我們必須進行能量最小化，以減輕系統(tǒng)中可能存在的任何不良接觸(原子重疊導(dǎo)致明顯的排斥，從而在模擬中引起數(shù)值問題)。

grompp –v –f em.mdp –c speptide_water.gro –p speptide.top –o speptide_em.tpr

或者簡單地使用這個命令:

grompp –v –f em –c b4em –p speptide –o em

在這個命令之后，壞的接觸已經(jīng)被移除。所以我們現(xiàn)在可以盡量減少能量消耗。

mdrun –v –s speptide_em.tpr –o speptide_em.trr –c after_em.gro –g emlog.log

或者

mdrun –v –s em –o em –c after_em –g emlog

屏幕可見：

Steepest Descents: Tolerance (Fmax) = 2.00000e+03 Number of steps = 100 Step= 0, Dmax= 1.0e-02 nm, Epot= -8.02562e+04 Fmax= 2.70468e+04, atom= 6792 Step= 1, Dmax= 1.0e-02 nm, Epot= -8.69641e+04 Fmax= 1.24531e+04, atom= 6657 Step= 2, Dmax= 1.2e-02 nm, Epot= -9.48067e+04 Fmax= 5.54072e+03, atom= 9009 Step= 3, Dmax= 1.4e-02 nm, Epot= -1.02743e+05 Fmax= 2.86957e+03, atom= 8766 Step= 4, Dmax= 1.7e-02 nm, Epot= -1.08941e+05 Fmax= 1.53075e+04, atom= 150 Step= 5, Dmax= 2.1e-02 nm, Epot= -1.10083e+05 Fmax= 1.31653e+04, atom= 150 Step= 6, Dmax= 2.5e-02 nm, Epot= -1.10564e+05 Fmax= 2.16454e+04, atom= 150 Step= 7, Dmax= 3.0e-02 nm, Epot= -1.11528e+05 Fmax= 1.82398e+04, atom= 150 Step= 9, Dmax= 1.8e-02 nm, Epot= -1.12767e+05 Fmax= 3.87548e+03, atom= 120 Step= 11, Dmax= 1.1e-02 nm, Epot= -1.13598e+05 Fmax= 1.23013e+04, atom= 120 Step= 12, Dmax= 1.3e-02 nm, Epot= -1.14464e+05 Fmax= 5.37656e+03, atom= 120 Step= 13, Dmax= 1.5e-02 nm, Epot= -1.14787e+05 Fmax= 1.74494e+04, atom= 120 Step= 14, Dmax= 1.9e-02 nm, Epot= -1.15874e+05 Fmax= 8.35820e+03, atom= 120 Step= 16, Dmax= 1.1e-02 nm, Epot= -1.16420e+05 Fmax= 6.65421e+03, atom= 120 Step= 17, Dmax= 1.3e-02 nm, Epot= -1.16807e+05 Fmax= 1.08106e+04, atom= 120 Step= 18, Dmax= 1.6e-02 nm, Epot= -1.17248e+05 Fmax= 1.12834e+04, atom= 120 Step= 19, Dmax= 1.9e-02 nm, Epot= -1.17425e+05 Fmax= 1.45145e+04, atom= 120 Step= 20, Dmax= 2.3e-02 nm, Epot= -1.17590e+05 Fmax= 1.67089e+04, atom= 120 Step= 22, Dmax= 1.4e-02 nm, Epot= -1.18845e+05 Fmax= 2.65542e+03, atom= 121 Step= 23, Dmax= 1.7e-02 nm, Epot= -1.19198e+05 Fmax= 2.39072e+04, atom= 120 Step= 24, Dmax= 2.0e-02 nm, Epot= -1.20588e+05 Fmax= 6.17310e+03, atom= 120 Step= 26, Dmax= 1.2e-02 nm, Epot= -1.20901e+05 Fmax= 1.04224e+04, atom= 120 Step= 27, Dmax= 1.4e-02 nm, Epot= -1.21229e+05 Fmax= 8.89909e+03, atom= 120 Step= 28, Dmax= 1.7e-02 nm, Epot= -1.21346e+05 Fmax= 1.51125e+04, atom= 120 Step= 29, Dmax= 2.1e-02 nm, Epot= -1.21686e+05 Fmax= 1.29390e+04, atom= 120 Step= 31, Dmax= 1.2e-02 nm, Epot= -1.22216e+05 Fmax= 3.21419e+03, atom= 120 Step= 33, Dmax= 7.5e-03 nm, Epot= -1.22523e+05 Fmax= 6.52631e+03, atom= 120 Step= 34, Dmax= 8.9e-03 nm, Epot= -1.22796e+05 Fmax= 6.31916e+03, atom= 120 Step= 35, Dmax= 1.1e-02 nm, Epot= -1.22998e+05 Fmax= 8.04922e+03, atom= 120 Step= 36, Dmax= 1.3e-02 nm, Epot= -1.23167e+05 Fmax= 9.93619e+03, atom= 120 Step= 37, Dmax= 1.5e-02 nm, Epot= -1.23317e+05 Fmax= 1.10029e+04, atom= 120 Step= 39, Dmax= 9.3e-03 nm, Epot= -1.23857e+05 Fmax= 1.88770e+03, atom= 810 writing lowest energy coordinates. Steepest Descents converged to Fmax < 2000 in 40 steps Potential Energy = -1.2385705e+05 Maximum force = 1.8876959e+03 on atom 81 Norm of force = 1.3663420e+04 NOTICE!!!!!!!!! If the potential enery after minimization is lower than -1.1e+05kJ/mol, it is acceptable and the structure can be used for MD caculations

三、gromac的一些筆記

用了一個星期，基本操作差不多了。自頂向下大概是這樣：1．核心步驟當然是mdrun了，當然后續(xù)的數(shù)據(jù)分析才是對理論水平的要求，先不討論。這里需要一個tpr輸入文件——由grompp生成的計算體系和計算控制詳細信息的二進制文件。計算生成trr軌跡文件、xtc壓縮的軌道文件、edr能量等狀態(tài)參量記錄文件、log計算記錄等。后續(xù)處理靠這些輸出文件和 gmx提供的命令。2.grompp預(yù)處理文件，將三種文件整合到一個tpr文件作為mdrun的輸入。分別是：

• mdp文件：計算控制文件，里面要給md模擬時的計算參數(shù)加以描述，控制時間、精度、溫度什么的。

• gro文件：分子坐標文件，記錄要模擬的體系中全部（要計算的）原子的坐標、名稱等。

• top文件：力場參數(shù)文件，記錄gro文件中的原子之間相互作用立場參數(shù)和原子分子類別信息的文件。

3．mdp文件一般可以看手冊根據(jù)需要自己填寫，或者看模版。4．gro、top文件，用editconf加周期邊界元胞，genbox加溶劑。那么最初的gro、top文件如何得到呢？

• 如果是從網(wǎng)站上直接當下來的pdb文件，用pdb2gmx轉(zhuǎn)化。（有時需要先去掉里面的水，有時需要加-ignh不考慮多余的氫，有時需要加一些離子（genion）來平衡蛋白質(zhì)上的電荷）

要是還要加入自己的分子，比如在MS什么的軟件中先畫了一個，存成pdb文件，這時不能直接轉(zhuǎn)，因為殘基gmx不認識，立場庫里面沒有描述，需要自己寫top文件。不過這很累，還要自己計算，又不會，怎么辦呢，網(wǎng)上推薦用ProDrg軟件（只找到網(wǎng)絡(luò)版的網(wǎng)址google一下就有），把pdb文件 copy進去，run一下就出來gro，和itp文件了（還有其他的一坨，用不到）。把自己gro文件里的內(nèi)容copy到轉(zhuǎn)出來的gro里，總原子數(shù)和序號改一下，再把itp文件include到老top文件中，[ molecules ] 標簽下面添上加入的分子名和數(shù)，就可以grompp了。