GROMACS其top/下面的force fields

GROMACS不同force fields的參數會放在各自的資料夾下面。以下列出GROMACS-2018.1目前支援的force fields，並以amber99sb為例子介紹下面的每個檔案的意思：

GROMACS其top/下的共用檔案

安裝完GROMACS以後，設定各種force field的參數會放在安裝目錄下的share/gromacs/top/
共用的參數會放在檔案裡面，像是：

GROMACS的檔案格式（file format）

在跑GROMACS的時候，會需要很多檔案。這些檔案的功用會用不同的副檔名來決定。以下介紹這些檔案的功能：

• 參數檔案（Parameter files）
• mdp：存放模擬所需的各種參數，grompp與convert-tpr會用到

用ggplot2畫盒鬚圖（boxplot）

下圖是用ggplot2這個套件來畫boxplot

GROMACS移除trajectory的PBC

使用GROMACS跑分子模擬後，可能因為設定PBC (Periodic Boundary Conditions)，結構某些部分跑到格子外面。在輸出結構的時候，移除PBC的方法如下：

1. printf "1\n1\n" | gmx trjconv -s md.tpr -f md.xtc \
   -pbc mol -ur compact -o md_noPBC.xtc
2. printf "1\n1\n" | gmx trjconv -s md.tpr -f md.xtc \
   -pbc nojump -center -o md_noPBC.xtc
3. printf "22\n22\n" | gmx trjconv -s md.tpr -f md.xtc \
   -n index.ndx -pbc cluster -fit translation -o md_noPBC.xtc

前面兩個適用在沒有ligand的狀況時，只輸出protein本身結構。第3個指令則是在有ligand的時候，需要同時輸出protein + ligand結構，這時候必須要先用make_ndx的指令（參考GROMACS選取某些原子成為一群），先把兩者放在一個group中，然後再輸出。

注意，在計算RMSD的時候，必須要先把結構的PBC移除，不然會得到錯誤的答案。詳細作法參考《GROMACS用RMSD估計構型變化》

_EOF_

HPCCHINA

這個部落格介紹關於超級電腦的處理器架構與網路架構，與一些處理的問題。2017三月到五月發表了五篇文章，然後就消失了，但這五篇文章給了許多有趣的訊息。以下就是這五篇文章：

• 带你从专业的眼光看老黄的新核弹－Deep Learning适用”显卡” Volta V100
  介紹Nvidia Volta架構的顯卡
• 世界十大超级计算机深度解(tiao)析(kan)
  2016年11月的前十名架構與網路分析，還有計算效能的方法簡介
• 用你听得(bu)懂的“人”话告诉你什么是量子计算机
• 两小时完成人类DNA序列比对？看INTEL，NVIDIA，IBM大乱战
  與DNA alignment有關的比較
• IBM Deep Learning “专用” S822LC 服务器深(fu)度(yan)测评
  PowerCPU + Nvidia P100的使用，有測試GROMACS 

_EOF_

讓xmgrace搭配設定檔自動出圖

《讓xmgrace預設出圖是PNG檔》介紹如何讓xmgrace直接轉圖檔，但此方法如果是在遠端機器，因為沒有X window，會出現DISPLAY error的狀況。解法則是使用gracebat

• gracebat -nxy in.xvg -param xmgrace.param -hardcopy -printfile output.png

GROMACS計算特定一點到多點的距離

跑完分子模擬，想要看一群對很多群的最小距離，可以用pairdist這個指令。指令使用看這篇《GROMACS計算所有C-alpha的距離》。這指令麻煩的地方在於需要進入指令一行一行輸入要計算不同群，下面的指令則是直接執行版：

• gmx pairdist -s md.tpr -f md.xtc -type min -ref "group Met1" -sel "group Ile50; group Leu28; group Asp122; group Glu17" -o pairdist.xvg

上面這個指令可以計算Met1到Ile50、Leu28、Asp122、Glu17這四個殘基隨時間改變的最小距離，存在pairdist.xvg這個檔案裡面。

• Group的選擇請參考《gmx make_ndx》的說明
• -ref與-sel的語法則是參考《Selection syntax and usage》

GROMACS mdrun平行處理與GPU加速

GROMACS執行分子模擬的程式叫做mdrun，要加速模擬速度除了平行處理（parallelization）還可以用GPU加速，官方的參考資料如下：

• Getting good performance from mdrun：
      給5.x版本之後與2016、2018版本用的說明文件
• Useful mdrun features：也是GROMACS-5.x版本之後，mdrun的一些例子
• Acceleration and parallelization：版本4.6之後，但不到5.x的加速方式

GROMACS將.gro檔案轉成.pdb格式

GROMACS預設的座標存成GRO格式，但一般都會用PDB格式來看。轉換指令為：

• gmx editconf -f xxx.gro -o xxx.pdb

GROMACS 2018

GROMACS 2018在2018年1月10日釋出，用來取代GROMACS 2016。Release Notes裡面提到在效能上重要的更新有：

• PME long-ranged interactions can now run on a single GPU, which means many fewer CPU cores are needed for good performance.
• Optimized SIMD support for recent CPU architectures: AMD Zen, Intel Skylake-X and Skylake Xeon-SP.

GROMACS用RMSD估計構型變化

跑完蛋白質分子模擬後，將分子的軌跡畫出來以外，最常做的就是估計其構型的變化。通常會計算backbone的root-mean-square deviation (RMSD)，在GROMACS是這樣做：

1. gmx trjconv -s md.tpr -f md.xtc \
   -pbc mol -ur compact \
   -o md_noPBC.xtc
2. gmx rms -s md.tpr -f md_noPBC.xtc -tu ns -o rmsd.xvg 

上面指令的意思是：

1. trjconv指令會移除periodic boundary conditions (PBC)的影響，不然後面計算出來的RMSD會是錯誤的。
2. rms指令則是計算RMSD，會有兩個步驟，先選fitting group（通常就是backbone atoms）、再選計算RMSD的group（通常也是backbone atoms）

修改GROMACS的.tpr檔案

GROMACS裡面的.tpr檔案（又稱run input file）是跑分子模擬最重要的檔案。這個檔案合併了下面幾項資訊：

• .gro file：分子結構資訊，每個原子的座標
• .top file：分子topology資訊，原子在分子如何互相連結的資訊
• .mdp file：模擬的參數設定
• .ndx file (optional)：原子的index file

要產生這個檔案，要用gmx grompp去把.gro、.top、.mdp檔案合併
要看這個檔案用gmx dump，要修改用gmx convert-tpr去修改，例如說：

1. gmx convert-tpr -s md.tpr -o md_ca.tpr（選C-alpha）
2. gmx convert-tpr -s md.tpr -until 20000 -o md_upto_20ns.tpr
3. gmx dump -s md.tpr > md.tpr.ascii

GROMACS無法用CUDA加速

系統是Ubuntu 16.04，在跑mdrun的時候出現下面問題

NOTE: Error occurred during GPU detection:

CUDA driver version is insufficient for CUDA runtime version 

Can not use GPU acceleration, will fall back to CPU kernels.

用/usr/local/cuda/samples/1_Utilities/deviceQuery測./deviceQuery出現

./deviceQuery Starting... 

CUDA Device Query (Runtime API) version (CUDART static linking) 

cudaGetDeviceCount returned 35 -> CUDA driver version is insufficient for CUDA runtime version Result = FAIL

原本以為是CUDA driver與CUDA runtime version（也就是CUDA toolkit version）搭配不起來（搭配的列表看這裡），測試的指令如下：

• nvcc --version    （看Cuda compiler driver、CUDA toolkit的版本）
• nvidia-smi      （看系統安裝的硬體與driver的版本）
• cat /proc/driver/nvidia/version    （同上，會顯示NVRM version）

但不管怎樣檢查，通通都是可以配合的版本。最後才想到，當初希望機器的顯示用內顯而不要用獨顯，因為這個緣故，造成程式無法去找到GPU來跑。所以下完下面的指令：

• prime-select nvidia

這樣之後，機器就可以順利用CUDA加速了

GROMACS計算所有C-alpha的距離

參考《模擬Tryptophan cage的分子動態》，做這個蛋白質分子模擬1 ns，最後會得到configuration file "md.tpr"與trajectory file "md.xtc"這兩個檔案。使用pairdist這個指令來計算此蛋白質中任意兩兩C-alpha原子之間的距離，存放在pairdist.xvg這個檔案中：

• gmx pairdist -s md.tpr -f md.xtc -o pairdist.xvg -ref "group \"C-alpha\"" -refgrouping none -sel "group \"C-alpha\"" -selgrouping none

上面指令的參數意思分別是：

• md.tp：configuration file
• md.xtc：trajectory file
• pairdist.xvg ：任意兩兩C-alpha原子之間的距離存放在這邊。參考pairdist的說明，如果reference有N群 （r1, r2 ... rN）、selection有M群 （s1, s2 ... sN），那總共就會有1 + N*M個欄位。第一個欄位是時間，後面接的就是距離r1-s1, r2-s1 ... rN-s1 再來是 r1-s2, r2-s2 ... rN-s2、...、r1-sM, r2-sM ... rN-sM
• -ref "group \"C-alpha\""：reference points是C-alpha原子
• -sel "group \"C-alpha\""：比較的對象是C-alpha原子
• -refgrouping none：把reference group拆成一個一個比較
• -selgrouping none：把selection group拆成一個一個比較

讓GROMACS支援新的residue

GROMACS裡面想要增加carbamylated lysine這個residue（在這邊叫做KCX），參考《Adding a Residue to a Force Field》，總共需要下面幾個步驟：

1. 新增 /opt/local/share/gromacs/top/FORCEFIELD/kcx.rtp 這個檔案
2. 新增 /opt/local/share/gromacs/top/FORCEFIELD/kcx.hdb 這個檔案
3. 修改 /opt/local/share/gromacs/top/residuetypes.dat，加入KCX

這邊不需要像文件上面說明的那麼多步驟是因為沒有多增加什麼特別的原子。上面三個檔案的意思分別是：

1. 新的residue其原子的電荷分佈、還有連結的共價鍵等
2. 新的residue上面如何加入氫原子
3. 讓GROMACS支援新的residue

若要增加Ligand，請參考Justin Lemkul的《Protein-Ligand Complex》的說明

_EOF_

GROMACS選取某些原子成為一群

在分析分子模擬產生的軌跡檔（trajectory file）時，有時候只會想要看某幾個特定原子在結構上面相對的位置，這時候就必須要自己選取感興趣的原子。在GROMACS-5.1.2裡面用下面的指令：

• gmx make_ndx -f md.tpr -o md.ndx

檔案md.tpr是跑模擬時會用到的輸入檔案，包含了模擬的設定、每個原子的座標軸等等（讀取這些設定參考《看GROMACS的binary files的內容》）；檔案md.ndx則是設定完成後儲存的index file。上述指令輸入後，會進入make_ndx這隻程式的交談式環境中，首先會跳出目前儲存在md.tpr裡面的不同原子組成的群：

• System：整個系統裡面所有的原子
• Protein：系統內組成蛋白質的原子
• Protein-H：組成的蛋白質的原子，除了H atom以外
• C-alpha：蛋白質主鏈上的C-alpha atom
• Backbone：應該蛋白質主鏈上的N、C-alpha、C這三個原子
• MainChain：應該是蛋白質主鏈上的所有原子
• MainChain+Cb：上面的再加上C-beta atom
• MainChain+H
• SideChain：組成蛋白質支鏈上的原子    
• SideChain-H：組成蛋白質支鏈上的原子，除了H atom以外
• Prot-Masses      
• non-Protein       
• Other                   
• Water：水分子     
• SOL：溶液分子，通常就是水分子        
• non-Water  
• Ion：系統中的離子                  
• Water_and_ions：水與離子

大約會有上面這幾群，要注意-H的意思是減掉H atom。

GROMACS用GBSA不支援鋅離子

用GROMACS-5.1.2，想要用GBSA（Generalized Born model augmented with the hydrophobic solvent accessible surface area (SA) term）這種implicit water來跑含有zinc ion的系統，結果出現下面的錯誤訊息：

Can't do GB electrostatics; the implicit_genborn_params section of the forcefield is missing parameters for 18 atomtypes or they might be negative.
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

應該是因為沒有鋅離子相關的參數

GROMACS續跑

跑分子模擬會有兩種狀況必須要續跑：一是跑到一半系統因為其他緣故關機，另外一種則是已經跑完但想要延長模擬的時間。在此用GROMACS-5.1.2為例子，說明如何續跑。

首先是最初跑模擬的指令，結構與設定都存在 md.tpr 這個檔案中，假設跑10 ns

• gmx mdrun -deffnm md -nt 8 -v 

參數 -deffnm md 意思是說，所有產生的檔案都用md開頭，像是 md.cpt與md_prev.cpt（checkpoint file）、md.xtc與md.trr （trajectory files）、md.edr（energy file）、md.gro（最終的結構）。參數 -nt 8 意思是用八個核心去跑，參數 -v 是碎念模式


第一種狀況：電腦因為其他原因把模擬的程式給關掉，假設只跑到7.7 ns

• gmx mdrun -deffnm md -cpi md.cpt -nt 8 -v -append

這時候要從7.7 ns開始續跑。檔案 md.cpt 是所謂的checkpoint file，一段時間系統會把狀態存在這個檔案裡面，有可能跳掉前存到7.5 ns。上面的指令就會去存取 checkpoint file，然後從存取的點開始繼續往後跑。參數 -append 指定續跑的資料要附加在原本的檔案後面。
¡¡ 注意：這邊一定要指定checkpoint file，不然就是直接從頭開始跑 !!


第二種狀況：已經跑完10 ns的模擬了，想要續跑到20 ns

• gmx convert-tpr -s md.tpr -until 20000 -o md_10_20.tpr
• gmx mdrun -deffnm md_10_20 -cpi md.cpt -nt 8 -v 
• gmx trjcat -f md.xtc md_10_20.xtc -o md_20ns.xtc

第一個指令是產生跑模擬的設定檔，設定要跑到20 ns（= 20000 ps）；第二個指令則是要求從前面跑出來的檢查檔 md.cpt，也就是10 ns開始跑。兩個指令就是要求續跑 10 - 20 ns。最後一個指令則是將前面兩個trajectories給連結起來，成為一個完整的 0 - 20 ns的檔案。

看GROMACS的binary files的內容

用GROMACS跑分子模擬時，許多檔案都會存成binary files（簡單說就是給機器讀的檔案、人沒辦法讀）。這些檔案包含了：

• Portable xdr run input file：以 .tpr 結尾
• Trajectory files：以 .xtc, .trr 結尾
• Energy file：以 .edr 結尾
• Checkpoint file：以 .cpt 結尾
• Hessian matrix：以 .mtx 結尾

在GROMACS-5.1.2裡面，使用 gmx dump 這個指令，可以把binary files變成text files，這樣才能夠讀的懂內容。例如說下面的指令把 md.tpr 轉成 md.tpr.ascii：

• gmx dump -s md.tpr > md.tpr.ascii

去讀 md.tpr.ascii 就可以看到裡面包含了分子模擬的設定，分子的結構等等資訊。
目前遇到的問題是可以把 md.tpr 這個binary file轉成ascii file，但是卻沒有辦法把這個檔案轉回去。這個檔案的結構是Portable XDR format，也許要去參考grompp這個指令的原始碼，看是如何轉檔案的...

_EOF_