分析梳理--分子動力學(xué)模擬的常規(guī)步驟二(Gromacs)

作者,Evil Genius

有人說,人生低谷期就是在救你的命,可是一直低谷,一般人也是受不了。

今天我們開始分子動力學(xué)第二步:定義單位盒子和添加溶劑。

上一步我們生成了蛋白拓?fù)浣Y(jié)構(gòu),文章在分析梳理--分子動力學(xué)模擬的常規(guī)步驟一(Gromacs)。

gmx pdb2gmx讀取.pdb或.gro文件,并讀取數(shù)據(jù)庫文件,向分子添加氫并在GROMACS(GROMOS)或可選的.pdb中生成坐標(biāo),格式以及GROMACS格式的拓?fù)?。隨后可以處理這些文件以生成運行輸入文件。那么在這個基礎(chǔ)上,我們要進(jìn)行下一步。

第一步輸出的日志需要注意一下信息

后續(xù)我們平衡電荷需要用。

我們首先來看看參數(shù)

gmx editconf 將通用結(jié)構(gòu)格式轉(zhuǎn)換為.gro、g96或.pdb在分子動力學(xué)模擬中,通常會給體系添加一個周期性的模擬盒子.gmx editconf有許多控制盒子的選項.
Options to specify input files:

 -f      [<.gro/.g96/...>]  (conf.gro)
           Structure file: gro g96 pdb brk ent esp tpr
 -n      [<.ndx>]           (index.ndx)      (Opt.)
           Index file
 -bf     [<.dat>]           (bfact.dat)      (Opt.)
           Generic data file

Options to specify output files:

 -o      [<.gro/.g96/...>]  (out.gro)        (Opt.)
           Structure file: gro g96 pdb brk ent esp
 -mead   [<.pqr>]           (mead.pqr)       (Opt.)
           Coordinate file for MEAD

Other options:

 -[no]w                     (no)
           View output .xvg, .xpm, .eps and .pdb files
 -[no]ndef                  (no)
           Choose output from default index groups
 -bt     <enum>             (triclinic)
           Box type for -box and -d: triclinic, cubic, dodecahedron,
           octahedron
 -box    <vector>           (0 0 0)
           Box vector lengths (a,b,c)
 -angles <vector>           (90 90 90)
           Angles between the box vectors (bc,ac,ab)
 -d      <real>             (0)
           Distance between the solute and the box
 -[no]c                     (no)
           Center molecule in box (implied by -box and -d)
 -center <vector>           (0 0 0)
           Shift the geometrical center to (x,y,z)
 -aligncenter <vector>      (0 0 0)
           Center of rotation for alignment
 -align  <vector>           (0 0 0)
           Align to target vector
 -translate <vector>        (0 0 0)
           Translation
 -rotate <vector>           (0 0 0)
           Rotation around the X, Y and Z axes in degrees
 -[no]princ                 (no)
           Orient molecule(s) along their principal axes
 -scale  <vector>           (1 1 1)
           Scaling factor
 -density <real>            (1000)
           Density (g/L) of the output box achieved by scaling
 -[no]pbc                   (no)
           Remove the periodicity (make molecule whole again)
 -resnr  <int>              (-1)
            Renumber residues starting from resnr
 -[no]grasp                 (no)
           Store the charge of the atom in the B-factor field and the radius
           of the atom in the occupancy field
 -rvdw   <real>             (0.12)
           Default Van der Waals radius (in nm) if one can not be found in the
           database or if no parameters are present in the topology file
 -[no]sig56                 (no)
           Use rmin/2 (minimum in the Van der Waals potential) rather than
           sigma/2
 -[no]vdwread               (no)
           Read the Van der Waals radii from the file vdwradii.dat rather than
           computing the radii based on the force field
 -[no]atom                  (no)
           Force B-factor attachment per atom
 -[no]legend                (no)
           Make B-factor legend
 -label  <string>           (A)
           Add chain label for all residues
 -[no]conect                (no)
           Add CONECT records to a .pdb file when written. Can only be done
           when a topology (tpr file) is present

KNOWN ISSUES

* For complex molecules, the periodicity removal routine may break down,
* in that case you can use gmx trjconv.
-f:[.gro/.g96/...>] (conf.gro)輸入蛋白結(jié)構(gòu)
-O:[.grol.g96/...>] (out.gro)輸出帶模擬盒子信息的結(jié)構(gòu)文件
-c:將蛋白分子置于盒子中心
-bt:盒子類型,默認(rèn)使用長方盒子,dodecahedron、cubic等
-d:溶質(zhì)與盒子之間的距離,在X,Y,Z方向上的最小距離
到盒子邊緣的距離是一個重要參數(shù),因為我們要使用周期性邊界條件,必須滿足最小映象約定,也就是說,一個蛋白質(zhì)永遠(yuǎn)不能“看到”它自身的周期性映象(不能與其自身有相互作用),否則計算的力就會含有虛假的部分.指定溶質(zhì)與盒子之間的距離為1.0nm意味著,蛋白質(zhì)分子的任意兩個周期性映象之間的距離至少是2.0nm.對于模擬中常用的任何截斷方案,這個距離都足夠了。-d選項設(shè)定分子到盒子邊緣的最小距離,以nm為單位,它決定了盒子的尺寸.理論上在絕大多數(shù)系統(tǒng)中,-d都不能小于0.9nm。
當(dāng)你對周期性的邊界條件與盒子類型有了更多了解后,強烈推薦使用菱形十二面體晶胞,因為在周期性距離相同的情況下,它的體積大約只有立方體晶胞的71%,因此可以減少需要加入的溶劑水分子的數(shù)目。
我們多擴展一下,周期性邊界(periodic boundary condition)
考慮一個極小的立方體液滴單元, 里面有一千個原子. 其中有83=512個在內(nèi)部, 可以看作完全處于液體環(huán)境中, 而剩下的接近一半, 都在這個單元的表皮上, 顯然不會有著相同的性質(zhì). 即使是到了106 個原子這樣一個量級, 仍然有接近6%有著不同的性質(zhì). 而我們所取得這個液體單元, 必須要是對宏觀普適的抽樣, 那么怎么做呢?
這個問題可以通過由Born and von Karman 在1912年提出的周期性邊界條件來克服. 第一種說法是, 就是假想空間中有無數(shù)個相同的液體元胞的展開, 在立方體邊界上的微粒, 依然可以受到鄰近的元胞的作用. 第二種說法, 臨近的元胞都是本體的鏡像. 落實到實處就是, 一個微粒運動出了盒子, 將從盒子的另一邊再穿進(jìn)來. 有點像小時后玩的貪吃蛇.
如此一來, 盒子便沒有了邊界, 也沒有了在表面的粒子, 這樣的盒子很容易地用坐標(biāo)的形式表現(xiàn)出來, 不用去儲存真實的坐標(biāo). 但是我們?nèi)杂修k法, 去得到這個元胞在真實體系中展開的情形. 折算到盒子內(nèi)部坐標(biāo)通常稱之為wrap, 而按真實展開的叫做unwrap.
有個非常關(guān)鍵的問題是, 這個非常小的無限循環(huán)的小東西, 能不能和宏觀的系統(tǒng)有著相同的性質(zhì)? 這取決于分子間作用力的范圍和所考察的范圍. 概括的說, 對于短程的強相互作用, 作用力距離不能超過盒子的一半. 這可能有一點抽象, 請看好Fig1.13, 如果粒子1的作用力距離超過一倍的盒子長, 那么會出現(xiàn)粒子自己對自己作用, 這顯然是make sense的; 如果超過一半的盒長, 那么對于某些粒子, 粒子1在正方向?qū)ζ渥饔靡淮? 從反方向同樣會再次對它作用, make no science sense. 這個就是最小鏡像原則(minimum image convention), 任意一個粒子的作用力僅被計算一次.
以立方體為例, 講一下思路. 要注意的是, 計算坐標(biāo)和計算距離是有所不同的. 這里盒子是放在原點, 邊長為L的立方體.
計算坐標(biāo)時, 我們不但想知道其在盒子中的位置, 也想知道實際上它跑到哪去了, 這個數(shù)據(jù)可以用來計算擴散等參數(shù). 這時我們需要引入一個image flag的概念. 這個值代表著原子穿出盒子的次數(shù), 正方向穿出image flag+1, 反方向-1.
當(dāng)微粒從正方向跑出盒子, 那就意味著它跑進(jìn)了下一個盒子, 坐標(biāo)值將大于L. 此時應(yīng)該減去L的值, 將其wrap進(jìn)盒子, 同時image flag+1; 如果從反方向跑了出去, 那就需要加一個L讓他wrap進(jìn)盒子, 同時image flag-1.
計算距離時分為兩種情況, 僅僅計算wrap體系中兩個微粒間作用力, 那么就有, coord1-coord2. 如果這個距離大于1/2L, 那L減去這個距離; 如果小于-1/2L, 那就需要加上L. 總之很簡單, 在紙上畫個坐標(biāo)軸取個特殊值就能搞明白.

我們來加一下盒子

gmx editconf -f conf.gro -o newbox.gro -c -bt cubic -d 1.0

然后我們添加溶劑gmx solvate ,我們還是看一下參數(shù)

Options to specify input files:

 -cp     [<.gro/.g96/...>]  (protein.gro)    (Opt.)
           Structure file: gro g96 pdb brk ent esp tpr
 -cs     [<.gro/.g96/...>]  (spc216.gro)     (Lib.)
           Structure file: gro g96 pdb brk ent esp tpr

Options to specify input/output files:

 -p      [<.top>]           (topol.top)      (Opt.)
           Topology file

Options to specify output files:

 -o      [<.gro/.g96/...>]  (out.gro)
           Structure file: gro g96 pdb brk ent esp

Other options:

 -box    <vector>           (0 0 0)
           Box size (in nm)
 -radius <real>             (0.105)
           Default van der Waals distance
 -scale  <real>             (0.57)
           Scale factor to multiply Van der Waals radii from the database in
           share/gromacs/top/vdwradii.dat. The default value of 0.57 yields
           density close to 1000 g/l for proteins in water.
 -shell  <real>             (0)
           Thickness of optional water layer around solute
 -maxsol <int>              (0)
           Maximum number of solvent molecules to add if they fit in the box.
           If zero (default) this is ignored
 -[no]vel                   (no)
           Keep velocities from input solute and solvent
gmx solvate可以做以下兩件事之一:
1.創(chuàng)建一個充滿溶劑的盒子.可以通過指定-cs和-box選項來完成.對具有盒子信息但不含原子的結(jié)構(gòu)文件則可以通過指定-cs和-cp來實現(xiàn).
2.將溶質(zhì)分子,如蛋白質(zhì)進(jìn)行溶劑化,使其處于溶劑分子的包圍之中.-cp和-cs分別用于指定溶質(zhì)和溶劑.不設(shè)定-box時,會使用溶質(zhì)坐標(biāo)文件(-cp)中的盒子信息.
-cp[<gro/.g96/...>] (protein.gro)輸入文件,使用的蛋白質(zhì)構(gòu)型文件來自前面editconf步驟中的輸出文件
-cs [<gro/.g96/...>] (spc216.gro)輸入庫文件,來自標(biāo)準(zhǔn)安裝的GROMACS的劑構(gòu)型文件
-p[<.top>] (topol.top)指定所要修改的拓?fù)湮募?/h5>
-o[<.gro/.g96/...>] (out.gro)輸出加過劑的結(jié)構(gòu)文件
使用的spc216.gro是通用的已平衡的三位點溶劑模型.也可以使用spc216.gro作為SPC,SPC/E或TIP3P水模型的溶劑構(gòu)型,因為它們都是三位點的水模型。
注意topol.top中[molecules]的變化(SOL那一項)
solvate記錄了增加的水分子數(shù)目,并將其寫入拓?fù)湮募幸燥@示它所做的更改,注意,如果你使用其他的(非水)溶劑,solvate不會在拓?fù)湮募袑懭脒@些信息!它自動記錄更新水分子的功能是直接寫在源代碼中的.
定義好了模擬盒子后,使用solvate模塊添加溶劑,輸出包含溶劑的結(jié)構(gòu)文件solv.gro.
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

看一下添加的效果

生活很好,有你更好。

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容