上回我們將拓撲文件都準備完全，并修改了top文件。接下來這一節(jié)是定義盒子添加溶劑。

此時，工作流程與其他 MD 模擬一樣。我們將定義單元格并在其中注入水。

gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0?

gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

在查看 solv.gro 時，您可能會想，為什么 editconf 沒有生成所需的十二面體單胞形狀，因為系統(tǒng)看起來是這樣的：

GROMACS 程序總是使用數值效率最高的坐標表示法，即把所有坐標重新包扎成三菱單元。mdrun 所執(zhí)行的物理計算可以通過不同的坐標包裹方式等效進行，因此首選效率最高的方式。在生成 .tpr 文件后，可以恢復所需的單元格形狀。

添加離子

現在我們有了一個包含帶電蛋白質的溶解系統(tǒng)。pdb2gmx 的輸出結果告訴我們，該蛋白質的凈電荷為 +6e（基于其氨基酸組成）。如果您在 pdb2gmx 的輸出中忽略了這一信息，請查看 topol.top 中 [ atoms ] 指令的最后一行；它應該是（部分）"qtot 6"。由于生命不是以凈電荷形式存在的，因此我們必須在系統(tǒng)中添加離子。

使用任何 .mdp 文件，使用 grompp 生成?.tpr 文件。我使用 .mdp 文件來運行能量最小化，因為它們需要的參數最少，因此最容易維護。

該mdp文件長這樣：

; LINES STARTING WITH ';' ARE COMMENTS

title ? ? = Minimization ; Title of run

; Parameters describing what to do, when to stop and what to save

integrator ? ? = steep ; Algorithm (steep = steepest descent minimization)

emtol ? ? = 1000.0? ; Stop minimization when the maximum force < 10.0 kJ/mol

emstep? ? ? ? ? = 0.01? ? ? ; Energy step size

nsteps ? ? = 50000 ? ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions

nstlist ? ? = 1 ? ? ; Frequency to update the neighbor list and long range forces

cutoff-scheme? ?= Verlet

ns_type ? ? = grid ; Method to determine neighbor list (simple, grid)

rlist ? ? = 1.0 ; Cut-off for making neighbor list (short range forces)

coulombtype ? ? = cutoff ; Treatment of long range electrostatic interactions

rcoulomb ? ? = 1.0 ; long range electrostatic cut-off

rvdw ? ? = 1.0 ; long range Van der Waals cut-off

pbc? ? ? ? ? ? ?= xyz ; Periodic Boundary Conditions

接下來運行命令：

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr

我們現在將 .tpr 文件傳遞給 genion

gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

選擇替換溶劑分子的group選擇SOL，也就是輸入15

在以前的 GROMACS 版本中，用 -pname 和 -nname 指定的離子名稱是針對特定力場的，但在 4.5 版中已經標準化。指定的原子名稱始終是大寫的元素符號以及[ moleculetype ]。殘基名稱可以附加或不附加電荷符號（+/-）。如果遇到困難，請參閱 ions.itp 了解正確的命名方法。

?[ molecular ] 現在應該如下所示：

[ molecules ]

; Compound ? ? ? ?#mols?

Protein_chain_A ? ? 1?

JZ4 ? ? ? ? ? ? ? ? 1?

SOL ? ? ? ? ? ? 10228?

CL ? ? ? ? ? ? ? ? ?6

能量最小化

現在系統(tǒng)已組裝完畢，使用 grompp進行操作，em.mdp文件如下：

gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr

確保在運行 genbox 和 genion 時更新了 topol.top 文件，否則你會收到很多討厭的錯誤消息（“坐標文件中的坐標數與拓撲不匹配”等）

現在我們可以調用?mdrun?來執(zhí)行?EM

gmx mdrun -v -deffnm em

現在，我們的系統(tǒng)處于能量最小值，可以開始真正的動力學研究了。

平衡

平衡我們的蛋白質配體復合物與平衡任何其他含有水中蛋白質的系統(tǒng)非常相似。在這種情況下，有一些特殊的考慮因素：

1. 對配體施加限制

2. 溫度耦合基團的處理

約束配體

為了限制配體，我們需要為其生成一個位置限制拓撲。首先，為 JZ4 創(chuàng)建一個只包含其非氫原子的索引組：

gmx make_ndx -f jz4.gro -o index_jz4.ndx

然后，執(zhí)行 genrestr 模塊，選擇新創(chuàng)建的索引組（即 index_jz4.ndx 文件中的第 3 組）：

gmx genrestr -f jz4.gro -n index_jz4.ndx -o posre_jz4.itp -fc 1000 1000 1000

現在，我們需要在拓撲中加入這些信息。根據我們希望使用的條件，有幾種方法可以做到這一點。如果我們只想在限制蛋白質的同時限制配體，那么請在拓撲中的指定位置添加以下幾行：

位置很重要！您必須在拓撲結構中調用 posre_jz4.itp。jz4.itp 中的參數定義了配體的 [ moleculetype ] 指令。分子類型以包含水拓撲（tip3p.itp）結束。在其他地方調用 posre_jz4.itp 會試圖將位置限制參數應用于錯誤的分子類型

如果您希望在平衡過程中進行更多控制，即獨立限制蛋白質和配體，您可以在另一個 #ifdef 塊中控制配體位置限制文件，如圖所示：

在后一種情況下，為了同時限制蛋白質和配體，我們需要在 .mdp 文件中指定 define = -DPOSRES -DPOSRES_LIG。如何處理您的系統(tǒng)取決于您。這些例子只是為了說明 GROMACS 提供的靈活性。對于標準平衡程序，同時限制蛋白質和配體可能就足夠了。您自己的需求可能會有所不同。

溫度調節(jié)

適當控制溫度耦合是一個敏感問題。將每個分子類型與自己的恒溫組耦合是個壞主意。例如，如果這樣做

tc-grps = Protein JZ4 SOL CL

您的系統(tǒng)很可能會爆炸，因為溫度耦合算法不夠穩(wěn)定，無法控制只有幾個原子的原子團（即 JZ4 和 CL）所產生的動能波動。請勿單獨耦合系統(tǒng)中的每一個種類。

典型的方法是設置 tc-grps = Protein Non-Protein，然后繼續(xù)。但是，"非蛋白質 "組也包括 JZ4。由于 JZ4 和蛋白質的物理聯系非常緊密，最好將它們視為一個整體。也就是說，為了溫度耦合的目的，將 JZ4 與蛋白質歸為一組。同樣，我們插入的少量 Cl- 離子也被視為溶劑的一部分。為此，我們需要一個特殊的索引組來合并蛋白質和 JZ4。我們可以使用 make_ndx 模塊，提供完整系統(tǒng)的任意坐標文件來實現這一目的。這里我使用的是 em.gro，即我們系統(tǒng)的輸出（最小化）結構：

gmx make_ndx -f em.gro -o index.ndx

將 "Protein "組和 "JZ4 "組合并，其中">"表示 make_ndx 提示符：

現在我們可以設置 tc-grps = Protein_JZ4 Water_and_ions，以達到我們想要的 "蛋白質非蛋白質 "效果。

使用此 .mdp 文件繼續(xù)進行 NVT 平衡。

該mdp文件如下：

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

gmx mdrun -deffnm nvt

NVT 模擬完成后，使用此 .mdp 文件繼續(xù)進行 NPT模擬：

該文件長這樣：

gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr

gmx mdrun -deffnm npt

MD模擬

完成這兩個平衡階段后，系統(tǒng)已在所需的溫度和壓力下進行了良好的平衡?，F在我們可以解除位置限制，運行生產 MD 進行數據采集。這個過程就像我們之前看到的一樣，因為我們將利用檢查點文件（在本例中現在包含保存壓力耦合信息）來 grompp。我們將運行 10-ns MD 仿真。所需要的文件長這樣

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_10.tpr

gmx mdrun -deffnm md_0_10

這一步時間要非常非常久，如果不用gpu加速的話要好幾個小時甚至幾天。所以做好心理準備。

下一節(jié)是分析部分