# Gromacs
###### tags: `Gromacs`
[toc]
## Gromacs 加到環境變數
```=python
#將gromacs 指令加到環境變數
source /pkg/GROMACS/setgromacs_2018
```
## packmol
### packmol使用說明
1. 要先建立一個pdb檔
2.  要先有一個inp檔
3. 網址處打cmd目前資料夾
4.  執行inp檔
### packmol 各分子同向
```python=
# x,y,z固定角度180+-20
constrain_rotation x 180. 20.
constrain_rotation y 180. 20.
constrain_rotation z 180. 20.
```
## Gromcas 檔案
### .trr與.gro檔
.trr是有每一步的速度與位置，.gro檔是只有末位置末速度。
 .trr檔案大
### .mdp 檔(mdp)
* [介紹網址](https://manual.gromacs.org/current/user-guide/mdp-options.html#mdp-value-integrator-steep)
* 
* [leap frog](https://manual.gromacs.org/documentation/2019/reference-manual/algorithms/molecular-dynamics.html#the-leap-frog-integrator)
* [neiborsearching,cutoff_scheme](https://computecanada.github.io/molmodsim-md-theory-lesson-novice/02-Fast_Methods_to_Evaluate_Forces/index.html)
(選擇在距離多少以後不算位能)
## Gromacs 小技巧
### 查看兩個檔案有甚麼不同
```python=
#查看檔案不同處
vimdiff constrait_em.gro(檔名) em_constraint_whole.gro(檔名)
```
## Gromacs 執行指令
### 基本指令(台灣山一號)
``` python=
#將gromacs 指令加到環境變數
source /pkg/GROMACS/setgromacs_2018
#將.pdb檔換成.gro檔 且條box的大小(unit:nm)
gmx_mpi editconf -f 1g.pdb -box 1.4 -o input.gro
#將全部檔案打包成.tpr檔(且能接受5個warning)
gmx_mpi grompp -f em.mdp -c input.gro -p H2O.top -o em.tpr -maxwarn 5
#開始計算且輸出檔名為 em(subg_gmx要做一些設定)
qsub -N em sub_gmx.sh
```
### 邊界連在一起 (-f(.trr軌跡檔)-o也可以做成.gro檔)
``` python=
#將em.gro的分子在視覺化軟體裡連在一起(用在角度作圖)(input file也可以是.trr每一步)
gmx_mpi trjconv -f em.gro -s em.tpr -pbc whole -o em_whole.gro
```
### Energy (可查看各種資訊)
``` python=
##每步能量的output(output 檔名.xvg)
gmx_mpi energy -f npt.edr -o
```
### RDF
* **RDF 解釋** :粒子徑向分佈函數 g(r)是指系統中所有粒子在空間的分佈情況，以一顆粒子為座標中心，觀察其它粒子在空間中隨著r的分佈情形。
* (藍色為找到原子，深紅為參考的原子)
* 
* 
* [相關網址](https://baike.baidu.hk/item/%E5%BE%91%E5%90%91%E5%88%86%E4%BD%88%E5%87%BD%E6%95%B8/12723225) [最後一張圖](https://www.eng.buffalo.edu/~kofke/ce530/Lectures/Lecture15/sld006.htm) [wiki](https://en.wikipedia.org/wiki/Radial_distribution_function) [wiki2](https://en.wikibooks.org/wiki/Molecular_Simulation/Radial_Distribution_Functions)
```python=
##創建 編號檔 進去可選分類 (a O, a H1 and H2)
gmx_mpi make_ndx -f npt.gro -o
## 顯現出rdf分布(system(1個水分子的質心)對system(一個水分子的質心))
gmx_mpi rdf -f npt.trr -s npt.tpr -selrpos whole_mol_com -seltype whole_mol_com -o
##顯現出rdf分布，可調整要甚麼原子對甚麼原子的rdf
gmx_mpi rdf -f npt.trr -s npt.tpr -selrpos atom -seltype atom -n index.ndx -o
gmx_mpi rdf -f heating.trr -s heating.tpr -selrpos whole_mol_com -seltype whole_mol_com
```
### 複製很多個box
```python=
## 把H2.gro 變大成5*5*5，所以有125個H2
gmx_mpi genconf -f H2.gro -nbox 5 5 5 -o super_H2.gro
```
* sub_gmx.sh檔
work:是超級電腦計算用的資料夾 > /work1/(自己主機的名字)
## Gromacs 作圖
### gnuplot
```python=
#先gnuplot進入作圖模式
gnuplot
#plot "檔名" using (X軸):(Y軸) 圖線的設定
plot "energy.xvg" using 1:2 w1
```
### python
#### load anaconda
```python=
#導入anaconda 模組
module load anaconda3/5.1.10
#確認python版本
python --version
#檢查pip
which pip
用pip安裝ase
pip install ase --user
```
#### 對角度作圖
1. 第一步打ipython可以在linux裡直接打python code
2. 建一個 plot_angdist.py檔
```python=
from ase.io import read,write
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import font_manager
# plot setting
plt.rcParams['axes.linewidth'] = 2
plt.rcParams['font.size'] = 12
# read angles from ase
atoms = read('em_constraint_whole.gro')
#H0H的角度(x=0~99)
ang_dist = [ atoms.get_angle(3*x+1,3*x,3*x+2) for x in range(100) ]
# plot histogram using matplotlib
fig,ax=plt.subplots()
ax.hist(ang_dist,density=True)
#ax.set_xlabel('x value',fontsize=20)
#ax.set_ylabel('y value',fontsize=20)
fig.tight_layout()
plt.show()
```
3.開始作圖
```python=
# &是讓圖不消失
python plotangdist.py &
```
#### 鍵長作圖 (一次好幾張圖)
``` python=
# read angles from ase
atoms = read('npt-whole.gro')
bndOH_dist = [ atoms.get_distance(3*x+1,3*x) for x in range(100) ]
bndHH_dist = [ atoms.get_distance(3*x+1,3*x+2) for x in range(100) ]
ang_dist = [ atoms.get_angle(3*x+1,3*x,3*x+2) for x in range(100) ]
# plot histogram using matplotlib
fig,ax=plt.subplots(nrows=3)#三張圖
ax[0].hist(bndOH_dist,density=True)
ax[1].hist(bndHH_dist,density=True)
ax[2].hist(ang_dist,density=True)ls
```
## Gromacs 上課筆記
### 介紹
* 原子跟原子尺度，force field用數學函數表示力，關注原子間的作用力。分子堆在一起排列成固態結構。
* step 1. FF的數學函數的parameters，描述力在不同材料、環境...。
* step 2. 找到文獻中有材料的FF parameter，對於未知材料的FF parameters要開發。
### 所需檔案
* .top file (topology): 定義鍵長、鍵角、二面角...等分子特徵。
* .itp file (FF parameters): force filed裡的參數
* Ex: OPLS forcefiled 描述有機分子的力場。分成bond,non-bonded
* bond: 由分子的震動行為，可運動的自由度，鍵，鍵角，二面角。
* non-bonded: ex:靜電作用力，凡德瓦力。
* 凡德瓦力: 當2原子靠近時，兩原子的電子雲高度重合會產生排斥力。
* OPLS立場各項的計算方式
* bond angle: 從OPLS database
* dihedral: 用QM dihedral scanning，轉分子的二面角記錄不同二面角時的能量
* atomic charge: reproduce QM electrostatic potential (E=-delta* V)
* ==.top file==
* 可以包含其他檔案(.itp)
* [ molecular type ]
* name (定義分子名稱) EX: other
* nrexcl (跟幾個鍵的原子會產生作用力) ex: 3
* [ atoms ]
* nr (定義分子中不同原子的編號) ex: 1,2
* type (分子中不同原子的種類(苯環中不同的C) ex: N,CN)
* resnr (group) ex: 1,1
* residue (殘留物)ex: LCZ,LCZ
* [ bond ]
* funct (查gromacs manual) 不同的鍵有不同數學函數 ex: 1
* [ pair ] [ angle ] [ dihedral ]
* [ molecular ] 要與 [ molecular type ]中的名子互相對應: 描述系統中有幾個分子要做模擬。
* ==.atp file== (要定義幾種atom type是一個難關)
* 區分幾種不同atom type，定義分子中甚麼原子可以定義成一樣的(不同位置的C，有些可以定義成一樣的但有些不能定義成一樣的)
* ==.rtp file== (residue topology)
* 把所有分子堆疊在一起形成一個大分子(像積木一樣)
* 從可以旋轉的鍵去定義不同的分子(積木)
* [ bonded type ]
* bonds (ex:1(harmonic)) angles (ex:1(harmonic)) dihedrals (ex:3) impropers (ex:2)
* [ ME1 ] (residue name)
* [ atoms ] (define atom)
* [ bonds ] (define bond) (angles程式會幫你定義)
* ==.pdb file==
* X,Y,Z element symbol (與rtp檔的 atom相互對應)
Sign in with Wallet
Connect another wallet