為了更好地對gromacs進行理解和學習，在這里將借gromacs官方實例的練習進行一個記錄。

蛋白配體復合物在gromacs中的步驟可以總結(jié)成：

生成拓撲文件 ————>? 定義盒子與溶劑?————>? 加入離子（目的是讓體系的凈電荷為0）?————>? 能量最小化??————>? NVT, NPT平衡??————>?MD??————>? 分析。

其步驟和之前基礎(chǔ)的溶菌酶在水體系的平衡相差不大，之后將用官方實例作一次練習。

在本次實例中，將使用T4溶菌酶與JZ4配體進行動力學模擬

步驟一 ： 準備蛋白拓撲

首先去RCSB下載T4溶菌酶的pdb文件（PDB: 3HTB）

在pymol中，T4溶菌酶的結(jié)構(gòu)顯示為這樣

接下來，需要將對T4溶菌酶去除結(jié)晶水

grep -v HOH 3HTB.pdb > 3HTB_clean.pdb

在準備好T4溶菌酶后，我們現(xiàn)在面臨的問題是，JZ4 配體在 GROMACS 提供的任何力場中都不是可識別的實體，因此如果您嘗試將此文件傳遞給 pdb2gmx，它將給出致命錯誤。僅當力場的 .rtp（剩余拓撲）文件中存在構(gòu)建塊條目時，才能自動組裝拓撲。由于情況并非如此，我們將分兩步準備系統(tǒng)拓撲：

1. 使用pdb2gmx準備蛋白質(zhì)拓撲

2. 使用外部工具準備配體拓撲

由于我們將分別準備這兩種拓撲結(jié)構(gòu)，因此必須將蛋白質(zhì)和 JZ4 配體分別保存到不同的坐標文件中。像這樣保存 JZ4 坐標。

grep JZ4 3HTB_clean.pdb > jz4.pdb

得到的JZ4長這樣

然后只需刪除 3HTB_clean.pdb 中的 JZ4 即可。至此，準備蛋白質(zhì)拓撲結(jié)構(gòu)就變得輕而易舉了。本教程將使用的力場是 CHARMM36，可從 MacKerell 實驗室網(wǎng)站獲取。在那里下載最新的 CHARMM36 力場壓縮包和 "cgenff_charmm2gmx.py "轉(zhuǎn)換腳本，我們稍后將使用該腳本。下載與您安裝的 Python 版本（Python 2.x 或 3.x）相對應的轉(zhuǎn)換腳本版本。

在下載完力場文件后在工作目錄中解壓力場壓縮包：

現(xiàn)在您的工作目錄中應該有一個 "charmm36-jul2022.ff "子目錄。用 pdb2gmx 生成一個?T4 溶菌酶的拓撲結(jié)構(gòu)：

gmx pdb2gmx -f 3HTB_clean.pdb -o 3HTB_processed.gro -ter

系統(tǒng)將提示你進行 3 項選擇。

1. 力場

2. 水模型

3. 終點類型

本教程選擇 CHARMM36 force field（選項 1），它列在 "從當前目錄 "列表的第一位。

選擇默認水模型（CHARMM-modified TIP3P），然后為末端選擇 "NH3+"和 "COO-"。由于 N 端殘基為蛋氨酸 (MET)，pdb2gmx 會選擇與碳水化合物不兼容的末端類型，因此這種交互式選擇是必要的。您必須選擇特定于蛋白質(zhì)的末端，否則會出現(xiàn)原子名稱不匹配的致命錯誤。

這里我出現(xiàn)了個PO4的錯誤，原因是沒有刪除PO4這個小配體，在實際操作中要將其他配體刪干凈.

步驟二：準備配體拓撲

我們現(xiàn)在必須處理配體問題。但是，對于力場無法自動識別的某些種類，我們該如何為其提供參數(shù)呢？正確處理配體是分子模擬中最具挑戰(zhàn)性的任務之一。力場作者花費數(shù)年時間推導出自洽的力場，而在此框架中引入一個新種類并非易事。任何新種類的力場參數(shù)都必須以與原始力場一致的方式推導和驗證。

對于 OPLS、AMBER 和 CHARMM 力場，這種推導通常采用各種量子力學計算的形式。這些力場的主要文獻介紹了所需的程序。對于 GROMOS 力場，參數(shù)化方法并不明確，主要依賴于凝聚相行為的經(jīng)驗擬合。也就是說，為每種原子類型計算一些初始電荷和倫納德-瓊斯參數(shù)，評估其準確性并加以改進。雖然最終結(jié)果非常令人滿意，即流體與現(xiàn)實世界中的流體非常相似，但推導過程可能會非常費力。

因此，自動化工具非常受歡迎。對于每種力場，都有一些方法或軟件程序聲稱可以提供與各種力場兼容的參數(shù)。但并非所有的方法或軟件都同樣準確。請參考下表中的幾個例子：

為 JZ4 添加氫原子

在本教程中，我們將使用 CGenFF 服務器生成 JZ4 拓撲。在訪問網(wǎng)站后注冊賬戶并激活。CGenFF 需要 Sybyl .mol2 文件作為輸入，以收集基本的原子類型信息和鍵合連接性。CHARMM 也是一種全原子力場，即明確表示所有 H 原子。晶體結(jié)構(gòu)通常不會指定 H 原子坐標，因此必須將其內(nèi)置。要生成 .mol2 文件并添加 H 原子，使用 Avogadro 程序。在 Avogadro 中打開 jz4.pdb，從 "Build（構(gòu)建）"菜單中選擇 "Add Hydrogens（添加氫原子）"。Avogadro 將在 JZ4 配體上構(gòu)建所有氫原子。保存名為 "jz4.mol2 "的 .mol2 文件（文件 -> 另存為...，從下拉菜單中選擇 Sybyl Mol2）。

以上是官方的表述，在這里，我使用openbabel對配體進行加氫。

obabel -ipdb jz4.pdb? -omol2 -Ojz4.mol2 -h

使用 CGenFF 生成 JZ4 拓撲

現(xiàn)在，jz4.mol2 文件已可用于生成拓撲結(jié)構(gòu)。訪問 CGenFF，登錄賬戶，然后點擊頁面頂部的 "Upload molecule"。上傳 jz4.mol2，CGenFF 服務器將以 CHARMM "stream"文件（擴展名為 .str）的形式快速返回拓撲圖。將其內(nèi)容從網(wǎng)站保存到名為 "jz4.str "的文件中。

CHARMM stream文件包含所有拓撲信息?- atom類型、?charges和?bonded連通性。它還包含額外的結(jié)合參數(shù)，這些參數(shù)是通過類比力場未涵蓋的內(nèi)部相互作用而生成的。CGenFF 還為每個參數(shù)提供了懲罰分數(shù)，即評估所分配參數(shù)的可靠程度。低于 10 的值可直接使用。10 - 50 之間的值意味著需要對拓撲結(jié)構(gòu)進行一些驗證，而大于 50 的懲罰值通常需要手動重新參數(shù)化。這種懲罰性評分是 CGenFF 服務器最重要的功能之一。許多其他服務器生成的拓撲結(jié)構(gòu)都是 "黑盒子"，用戶只能默默信任。將整個研究項目寄托在未經(jīng)驗證的自動程序上是非常危險的。糟糕的配體拓撲會導致大量時間的浪費和不可靠的結(jié)果。請務必驗證新參數(shù)化的拓撲結(jié)構(gòu)！至少要根據(jù)力場中的現(xiàn)有分子檢查配體的電荷量和原子類型。

檢查 jz4.str 的內(nèi)容，查看電荷和新二面參數(shù)的懲罰。所有懲罰值都非常低，表明該拓撲結(jié)構(gòu)質(zhì)量非常好，可以直接用于我們的模擬。

CHARMM 格式的 JZ4 拓撲很好，但如果我們要在 GROMACS 中運行，它就沒用了。請使用我們之前從 MacKerell 網(wǎng)站下載的 cgenff_charmm2gmx.py 腳本。腳本名稱中會包含版本 _py2 或 _py3，但為簡單起見，此處省略了這些信息，但您需要使用下載的 Python 版本腳本的實際名稱。使用以下命令執(zhí)行轉(zhuǎn)換。

python cgenff_charmm2gmx_py3_nx2.py JZ4 jz4.mol2 jz4.str charmm36-jul2022.ff

在運行這個程序的時候，在我目前的python3.9的環(huán)境中會出現(xiàn)一些錯誤，需要對腳本進行針對性的修改才能跑通。要注意環(huán)境的問題。

跑完顯示：

該配體引入了不屬于現(xiàn)有力場的新鍵合參數(shù)，這些參數(shù)被寫入一個名為 "jz4.prm "的文件，其格式為 GROMACS .itp 文件。我們稍后將處理該文件，但必須注意到它的存在。現(xiàn)在配體拓撲結(jié)構(gòu)被寫入 "jz4.itp"，其中包含配體[?moleculetype?] 定義。

通過 pdb2gmx，我們得到了一個名為 "3HTB_processed.gro "的文件，其中包含經(jīng)過處理的、符合力場的蛋白質(zhì)結(jié)構(gòu)。我們還有來自 cgenff_charmm2gmx.py 的 "jz4_ini.pdb "文件，其中包含所有必要的 H 原子，并與 JZ4 拓撲中的原子名稱相匹配。使用 editconf 將此 .pdb 文件轉(zhuǎn)換為 .gro 格式：

gmx editconf -f jz4_ini.pdb -o jz4.gro

將 3HTB_processed.gro 復制到一個新的待處理文件中，例如將其命名為 "complex.gro"，因為將 JZ4 添加到蛋白質(zhì)中將生成蛋白質(zhì)配體復合物。接下來，只需復制 jz4.gro 中的坐標部分并將其粘貼到 complex.gro 中，位置在最后一行蛋白質(zhì)原子的下方和方框矢量之前，就像這樣：

由于我們已在 .gro 文件中添加了 22 個原子，因此請增加 complex.gro 的第二行以反映此更改?，F(xiàn)在坐標文件中應該有 2636 個原子。

構(gòu)建拓撲

在系統(tǒng)拓撲中加入 JZ4 配體的參數(shù)非常簡單。只需在包含位置限制文件后，在 topol.top 中插入一行 #include "jz4.itp"即可。位置限制文件的包含表示 "蛋白質(zhì) "分子類型部分的結(jié)束。

配體引入了新的二面性參數(shù)，這些參數(shù)由 cgenff_charmm2gmx.py 腳本寫入 "jz4.prm"。在 topol.top 頂部插入 #include 語句，添加這些參數(shù)：

這個 #include 語句的位置非常重要 ‘-’? 必須出現(xiàn)在任何 [ moleculetype ] 條目之前，因為在構(gòu)建任何分子之前必須先定義所有參數(shù)。它還必須出現(xiàn)在父力場的 #include 語句之后，因為在引入使用原子類型的鍵合參數(shù)之前，必須先了解所有原子類型。

最后要進行的調(diào)整是在 [ molecular ] 指令中。考慮到complex.gro中有一個新分子，我們必須將其添加到此處，如下所示：

現(xiàn)在，拓撲結(jié)構(gòu)和坐標文件與系統(tǒng)內(nèi)容一致。

下一節(jié)：添加溶劑