Gromacs分子動力學模擬主要可以分為以下幾個步驟,不同的體系步驟可能略有不同。
在開始之前,先簡單了解一下預平衡:
分子動力學模擬的最終目的是對體系進行抽樣,然后計算體系的能量,各種化學鍵,成分分析等等。打個比方說,我們有一個蛋白質,我們想將它放入一種溶液中(可能是水,也可能不是),然后看看這個體系的能量如何變化,蛋白質的化學鍵,與水分子形成的氫鍵等等信息,那么我們需要將蛋白質放入溶液中,映射到現實中就是講溶劑放入溶劑中,然后等體系穩定后,觀察其性質。
在MD中,這一過程不向現實中一樣是自然發生的,我們需要通過模擬是體系演化到平衡狀態,這就是預平衡。一般來說預平衡會有以下辦法:
- 蛋白質結構能量最小化:PDB文件都是從晶體中獲得的,所以蛋白質放入溶液中后必然會發生變化,這就需要對其進行能量最小化,確保蛋白質的結構是穩定結構。
- 蛋白質位置限定性模擬:有時加入溶劑后,分子間相互作用力會過大,導致蛋白質體系崩潰。這時我們需要限制蛋白質中重原子的位置,維持其結構,等溶劑分子弛豫之后再放開限制進行模擬。
- NVT預平衡,NPT預平衡:一般先做NVT模擬,減小盒子內壓力,然后再做NPT模擬。
以上步驟當然不用全做,視情況而定,不過一般蛋白質能量最小化和位置限定性NPT還是要做的。
以下是分子動力學模擬的步驟,有些步驟可以省略。
- 獲取並處理PDB文件
一般PDB文件是從網站上下載,如http://www.rcsb.org/pdb/home/home.do。獲取PDB文件后有可能還要做一些處理,如末端氫原子,結晶水,等等。視情況而定。
2. 使用pdb2gmx獲得拓撲文件
命令pdb2gmx的詳細信息可以參加http://manual.gromacs.org/programs/gmx-pdb2gmx.html。具體的命令參數我會在另一篇文章中詳述。一般而言,我們使用時會是向下面這樣:
gmx pdb2gmx -ff amber99sb-ildn -f *.pdb -o *.gro -p *.top -water tip3p
-ff 選項,制定要使用的力場;
-f選項,制定輸入的PDB文件;
-o選項,制定生成的gro文件名
-p選項,制定要生成的拓撲文件名
-water選項,制定要使用的水分子模型
注意,除了生成*.gro文件和*.top文件之外,還會生成一個posre.itp,位置限定性文件(我把它理解成position-restraints的縮寫)。
如果不使用-ff選項的話,指令運行后會讓你自行選擇力場。
3. 定義盒子
定義盒子和填充溶劑可以看做一步,在這里為了詳細就分開來說。
與前面一樣,涉及到的命令及文件都在其他文章中詳述,下文不再贅述。使用editconf命令創建盒子:
gmx editconf -f *.gro -o *.gro -c -bt cubic -d 1.2
-f:指定輸入的蛋白結構
-o:指定輸出帶盒子信息的結構文件
-c:將蛋白質置於盒子的中心,這個選項是可選的,不必須。
-d:蛋白質與模擬盒子在XYZ方向上的最小距離,一般不能小於0.9nm
-bt:指定盒子類型,這里使用了立方盒子,還可以用八面體,十二面體等。
這樣我們就得到了周期型立方格子中的蛋白質分子。
editconf命令可以用於gro文件與pdb文件的相互轉換。用-f指定源文件,-o指定所需文件名即可。
4. 蛋白質真空中的能量最小化(非必須)
一般而言這一步不是必須的,不過這里還是簡述一下。如果我們只需要在真空中進行能量最小化的化,下一步就可以直接成品模擬了。
Gromacs使用grompp指令(GROMacs Pre-Preocessor)對帶有格子信息的gro文件與蛋白質的拓撲文件,還有mdp文件進行處理,從而得到用於mdrun的輸入文件*.tpr。tpr為二進制文件。具體指令如下:
gmx grompp -f *.mdp -c *.gro -p *.top -o *.tpr
-f:指定輸入參數文件。mdp文件會有專門的文章敘述
-c:指定輸入結構文件
-p:指定輸入拓撲文件
-o:指定用於mdrun的tpr文件
運行之后我們得到*.tpr文件和參數文件mdout.mdp
然后使用mdrun命令運行能量最小化:
gmx mdrun -v -deffnm *
-v:顯示模擬過程中的信息
-deffnm:我把它理解成define-file-name的縮寫。定義輸出文件名,文件后綴會自動加上。
運行后得到日志文件*.log,全精度軌跡問價*.trr,能量文件*.edr,能量最小化的結構文件*.gro。
5. 向盒子中填充溶劑
其實這只是一小步,同上,為了詳細我把它單獨列為一步。
使用solvate命令填充溶劑,以水為例:
gmx solvate -cp *.gro -cs *.gro -o *.gro -p *.top
-cp:指定需要填充水分子的體系,即前面我們用editconf得到的帶格子的結構文件
-cs:指定要使用的水模型
-p:指定體系的拓撲文件(原蛋白質的拓撲文件),這樣solvate就可以修改體系的拓撲文件。
-o:指定填充水分子后的輸出文件
運行之后我們可以得到得到-o所指定的文件,並且-p指定的top文件也會發生改變。
6. 添加離子
向盒子中添加溶劑之后,我們得到了一個帶電荷的溶液體系,因此必須進行中和。GROMACS中添加離子的指令是genion(我把它理解成generate-ion的縮寫),但是不巧的是genion需要的輸入文件為tpr文件。跟前面一樣,這需要grompp(GROMacs Pre-Processor)來產生。grompp可以處理坐標文件和拓撲(描述分子的文件)從而產生原子級別的輸入文件,即tpr文件,tpr文件包含了體系中所有原子的參數。
為了將坐標信息(gro)和拓撲信息(top)結合起來,我們需要一個mdp文件。mdp文件通常用於進行能量最小化,這里只是簡單的生成tpr文件。
gmx grompp -f *.mdp -c *.gro -p *.top -o *.tpr
-f:指定mdp文件
-c:指定結構文件(加入溶劑后的結構文件)
-p:指定拓撲文件(還是之前生成的蛋白質拓撲,當然在加入溶劑時該文件發生了變化)
-o:指定輸出文件
得到tpr文件后,就可以在其中加入離子了:
gmx genion -s *.tpr -o *.gro -p *.top - pname * -nname * -nn *
-s:將上述生成的tpr文件作為輸入
-o:生成新的結構文件
-p:再次改變top文件,反應蛋白質結構的改變
-pname:指定要添加的陽離子名稱,后面未指定數量,即為不添加
-nname:指定添加的陰離子名稱
-nn:添加的陰離子數目
7. 能量最小化
現在,我們定義了盒子(周期性邊界條件),溶劑分子,離子。整個體系已經到達電中性。在進行模擬之前,我們必須確保體系的結構正常,原子間距離不要太近,結合構型合理。這就需要對結構進行弛豫,這一過程稱之為能量最小化(EM,energy minimization),是MD中非常重要的一步。
與前面類似,依然是需要用grompp來產生tpr文件,首先要定義一個minim.mdp文件,定義好之后:
gmx grompp -f minim.mdp -c *.gro -p *.top -o *.tpr
-f:指定mdp文件
-c:指定結構文件
-p:指定拓撲文件
-o:指定輸出的文件名
得到tpr文件后就可以進行能量最小化了
gmx mdrun -v -deffnm em
mdrun的指令與前面一樣。
我們將得到以下文件:
*.log 日志文件,記錄了能量最小化過程
*.edr 二進制能量文件
*.trr 全精度的二進制軌跡文件
*.gro 能量最小化的結構
現在,我們的體系已經處於能量最小點了,可以做一些真正的模擬了!
8. NVT平衡
NVT平衡實際上是很重要的一步,但是它的核心在於mdp文件,而mdp文件我將在另一篇文章中單獨闡述,因此這里對於NVT模擬就簡化處理。
在一開始的pdb2gmx中我們生成了一個posre.itp 文件,這里終於派上用場了!它的作用是對蛋白質中的重原子(非氫原子)施加位置限制力。施加限制之后,這些原子就不能隨便移動,除非能量非常大。這樣做的目的在於平衡蛋白質周圍分子的同時而不引起蛋白質結構的變化。
定義好mdp文件后,就可以進行模擬了。
gmx grompp -f *.mdp -c *.gro -p *.top -o *.tpr
-c指定前面生成的能量最小化的結構文件,-p依然指向那個被修改了多次的蛋白質top文件,-o指定輸出文件。
gmx mdrun -deffnm *
指定輸出文件名。
9. NPT平衡
與NVT平衡類似,關鍵在於mdp文件中,因此不再贅述,命令如下
gmx grompp -f *.mdp -c *.gro -t *.cpt -p *.top -o *.tpr
gmx mdrun -deffnm *
cpt為斷點文件(check point),詳見關於文件的文章中。
10. 成品MD
現在我們的體系已經在需要的溫度和壓強下平衡(弛豫)好了,我們可以放開位置限制並進行最終的MD,以收集數據了。
同樣,先定義mdp文件,然后運行
gmx grompp -f *.mdp -c npt.gro -t npt.cpt -p *.top -o *.tpr
gmx mdrun -deffnm *
11. 分析
暫略,做到這里再補充。