---
# System prepended metadata

title: GROMACS 建立 FF parameters

---

# GROMACS 建立 FF parameters
[toc]
## OPLS 
* 有機，密度(density)準，二面角的呈現是DFT:平面(因為N跟H的吸引力，二面角是平面,XTB:雙峰(因為忽略N跟H的吸引力，不是平面)
* 一大氣壓動能5.多KJ/mol(3/2*N*KB*T)，無法穿越高於此的二面角
## 步驟
1. 先用gv畫結構 (.com) 檔 (最下面要打molecule(ex: benzene.wfx)

* 下面是.com檔範例
![](https://i.imgur.com/vz0BfED.png)

```
%mem = 80000MB
#P opt wb97xd/6-311g(d, p)
nosymm density=current output=wfx pop=NO

Title Card Required
```
2. 使用gaussian優化結構
```=linux
# -k是產生.sh script
subg16 -m 20GB -p 4 -k molecule(ex: benzene).com
# vi .sh檔，把chk改wfx
sbatch  /./././().sh
```
![](https://i.imgur.com/x63MGtI.png)


3. 用obabel把log 檔換成pdb檔
```
obabel -ilog file.log -O file.pdb
```

4. 建atp檔和rtp檔

`.atp範例
`
![](https://i.imgur.com/HttLfgQ.png)

`.rtp範例
`
![](https://i.imgur.com/OADh1Ng.png)


5. 打開molden，把pdb檔的residule name對應到rtp檔
`.pdb`
![](https://i.imgur.com/Upmgkue.png)


### dihedral scaning
1. 加到.pdb檔
``` =linux
# 2 1 26 27 dihedral(接不同rtp) 
# 0.0 360.0 10.0 ( 0~360角度，掃角度每10.0度，3是dihedral的function)
REMARK SCAN  2  1 26 27 0.0 360.0 10.0  3
```
2. 資料夾SCANS => 資料夾2_1_26_27(連接2個.rtp的4個原子) => 

3. 在MATERIALS_LIBRARY的資料夾用python畫dihedral  (要改compound name)
```=linux
# 改compound name
vi dihedral.py
python dihedral.py
```
3. 到SCANS資料夾裡的4個原子座標資料夾寫一個script執行
```=python
for a in {1..36}
do

    /home/link/myscripts/g16/subg16 -p 4 -m 20GB QM_$a.com
done

```
```=python
subg16 -p 12 -m 60GB QM_$a.com


```
#### 使用openbabel轉成連續的xyz檔
1. 寫一個script把所有的log檔轉成xyz檔
* 轉QM
![](https://i.imgur.com/l8ZuAa1.png)
```=linux
rm traj.xyz
touch traj.xyz
for i in {1..36}
do
    obabel -ilog QM_$i.log -O QM_$i.xyz
    cat traj.xyz QM_$i.xyz >> traj.xyz
done
```
* 轉QM和MD
```
rm traj_QM.xyz
rm traj_MD.xyz
touch traj_QM.xyz
touch traj_MD.xyz

for i in {1..36}
do
    obabel -ilog QM_$i.log -O QM_$i.xyz
    cat traj_QM.xyz QM_$i.xyz >> traj_QM.xyz
    obabel -ipdb confout_$i.pdb -O MD_$i.xyz
    cat traj_MD.xyz MD_$i.xyz >> traj_MD.xyz

done
```
2. 用molden 打開traj.xyz看正不正常

3. 檢查是否有全部算完check.py
```=python
for i in range(1, 37):
    file_name = "QM_{}.log".format(i)
    try:
        with open(file_name, "r") as f:
            lines = f.readlines()
            print(file_name)
            if 'Normal' not in lines[-1]:
                print('Not successful!!!!!!!!!!!!!!!!!!!!!!!')
            else:
                print('OK')
    except FileNotFoundError:
        print("File {} not found.".format(file_name))
```
### 建立.itp
[LipParGen](http://zarbi.chem.yale.edu/ligpargen/)
轉成.top可看harmonic常數
若.top編號沒有對齊可以轉成.gro看編號

### DDEC6 atomic charge 



1. 複製job_control.txt到QM_GAS_PHASE
```=linux
#job_control.txt的地址
/home/link/training/4_Naomi/01_FF_Parameterization/EXAMPLE/mCP/QM_GAS_PHASE/job_control.txt
```
2. job_control.txt改molecular(ex: PPD).wfx

3. load chargemol模組 且執行
```=linux
#load chargemol
module load tools/chargemol/2017

nohup Chargemol_09_26_2017_linux_serial < job_control.txt &
```

#### DDEC6好處
*　描述electrostatic potential 不錯，且可得Lennard Jones值(sigma,epsalon) form。
*　
### sigma,eps 
```= linux
#複製到QM的資料夾
/home/link/training/4_Naomi/01_FF_Parameterization/EXAMPLE/mCP/QM_GAS_PHASE/gen_VDW.py
```

1. 改 gen_VDW.py 讀.pdb,.rtp 寫成dic{atom name : atom type}
2. dic {atom type :[1,2,5,6](第幾個原子是此atom type)}

### MD,QM dihedral fitting
1. 先換python的版本
```=linux
#轉成舊版本
conda activate fitting
#解除舊版本
conda deactivate
```

2. plot MD dihedral
```=linux
#在MATERIAL plot
python dihedral.py
```
3. 把fit出來的二面腳參數改到.itp裡
4. 在6_1_10_11資料夾再新增資料夾並且把QM_* 移到:新資料夾

### 創一個叫MD_Single的資料夾 (和QM_GAS_PHASE同一層)
#### FF
1. 創FF、XTB資料夾
2. 進FF再創EM、NVT-300K資料夾
3. 進EM創single_precision資料夾
4. 進single_precision把之前的.ff資料夾、.top檔、.gro檔複製過來 (有.ff資料夾、.top檔、.gro檔、gen.sh、subgmx.sh、==em.mdp==)
5. 修改 gen.sh
![](https://i.imgur.com/nn1ucwG.png)
6. 改.gro檔 box大小
![](https://i.imgur.com/bZ8jOhF.png)
7. ```bash gen.sh```，產生em.tpr
8. ```sbatch -J em subgmx.sh```
9. 進NVT-300K資料夾，複製老師的檔案(==nvt.mdp==、gen.sh、subgmx.sh)、複製single_precision資料夾的檔案(em.gro、PPD.top、PPD.ff資料夾)
10. ```bash gen.sh```，產生nvt.tpr
11. ```sbatch -J nvt subgmx.sh```
12. ```gmx_mpi trjconv -f nvt.trr -pbc whole -s nvt.tpr -o traj-pbc.gro -center```
13. ```obabel -igro traj-pbc.gro -O traj.xyz```
14. 將traj.xyz複製到自己電腦，用VMD

#### XTB
1. 把 gaussian 優化完的xyz檔丟到MD_SINGLE/XTB，改名為input.xyz
2. 資料夾內要有：input.xyz, job.sh , md.inp
3. sbatch -J input job.sh
```
sbatch -J input job.sh
```
4. 算完把xtb.trj改名為traj.xyz
```mv xtb.trj traj.xyz```
## xtb跟ff的比較
1.把xtb跑完的輸出檔xtb.trj改檔名變成traj.xyz
2.traj-pbc.gro用obabel轉成traj.xyz
```obabel -igro traj-pbc.gro -O traj.xyz```
3.跑python
```import os
import math
from ase.io import read, write
import numpy as np
import matplotlib.pyplot as plt


# Edit the font, font size, and axes width
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 2

DIHEDRALS = ['2_1_10_12']
fig, ax = plt.subplots()
for i, dih in enumerate(DIHEDRALS):
    al = [ int(x)-1 for x in dih.split('_')  ]
    print ('computing Dih:',dih)
    # read traj from gromacs
    snapshot_morph = read('../FF/NVT-300K/traj.xyz',':')
    DIH_ff = [ traj.get_dihedral(al[0],al[1],al[2],al[3]) for traj in snapshot_morph ]
    # read traj from xtb
    snapshot_morph = read('../XTB/NVT-300K/traj.xyz',':')
    DIH_xtb = [ traj.get_dihedral(al[0],al[1],al[2],al[3]) for traj in snapshot_morph ]
    if i > 1:
        ax[i].hist(DIH_ff,bins=200,density=True, alpha=0.5,label='ff')
        ax[i].hist(DIH_xtb,bins=200,density=True, alpha=0.5,label='xtb')
        ax[i].set_title(dih)
        ax[i].set_xlabel(r'Dihedral ($^o$)',fontsize=18)
        ax[i].set_ylabel('Distribution',fontsize=18)
    else:
        ax.hist(DIH_ff,bins=200,density=True, alpha=0.5,label='ff')
        ax.hist(DIH_xtb,bins=200,density=True, alpha=0.5,label='xtb')
        ax.set_xlabel(r'Dihedral ($^o$)',fontsize=18)
        ax.set_ylabel('Distribution',fontsize=18)
        ax.set_title(dih)

ax.legend()
plt.tight_layout()
plt.show()
```
4.比較結果
![](https://i.imgur.com/tpyaxYb.png)

## Dihedral使用function 8
* table_gen.py
```python=
import numpy as np
import matplotlib.pyplot as plt

dihedrals = ['1_5_8_9', '3_4_10_11','2_3_12_14', '14_12_16_19', '3_12_14_15', '12_16_19_20']
angles = [i - 360 if i > 180 else i for i in range(0, 356, 5)]
hat_to_kjmol = 2625.5
compound = 'AA'


for index, di in enumerate(dihedrals):
    qm = []
    md = []
    for i in range(1, 73):
        with open('./' + di + '/QM_' + str(i) + '.pdb', 'r') as f:
            lines = f.readlines()
            for line in lines:
                if 'QM' in line:
                    word = line.split()
                    qm.append(float(word[3]) * hat_to_kjmol)
                    break
        with open('./' + di + '/MD_' + str(i) + '.pdb', 'r') as f:
            lines = f.readlines()
            for line in lines:
                if 'Potential' in line:
                    word = line.split()
                    md.append(float(word[3]))
                    break
    qm = np.array(qm) - qm[36]
    md = np.array(md) - md[36]
    function_x = qm - md
    g_function_x = -(np.gradient(function_x, edge_order=2))

    fig, ax = plt.subplots()
    ax.scatter(angles, function_x, label = 'f(x)', marker = 'x')
    ax.scatter(angles, g_function_x, label = "-f'(x)", marker = '^')
    ax.scatter(angles, qm, label = "qm")
    ax.scatter(angles, md, label = "md")
    ax.legend()
    ax.set_title(di)
    plt.show()

    info = {}
    for i, (a, b, c) in enumerate(zip(angles, function_x, g_function_x)):
        line = str(a) + '\t\t' + '{:.15f}'.format(b) + '    \t\t' + str(c)
        info[i + 1] = line
    with open('../FORCEFIELD/'+ compound + '.ff' + '/table_d' + str(index) + '.xvg', 'w') as f:
        first = info[37].replace('180', '-180')
        f.write(first + '\n')
        for i in range(1, 73):
            if i <= 35:
                f.write(info[i + 37] + '\n')
            else:
                f.write(info[i - 35] + '\n')
```
* subgmx.sh
```=bash
#!/bin/bash
#SBATCH -A MST111353
#SBATCH -J gromacs
#SBATCH -o ./gmx.stdout_%j
#SBATCH -e ./gmx.stderr_%j
#SBATCH --mail-type=end
#SBATCH --mail-type=fail
#SBATCH --mem=5000
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1
##SBATCH --time=4-00:00:00
#SBATCH --partition=ctest


module purge
module load app/gromacs/2019.6

#!/bin/bash

job=$SLURM_JOB_NAME

ff='./HA.ff/'

table_files="$ff""table_d6.xvg $ff""table_d7.xvg $ff""table_d8.xvg $ff""table_d9.xvg $ff""table_d10.xvg"

srun gmx_mpi mdrun -tableb $table_files -v yes -s ./$job.tpr -deffnm $job
```

## note
### 複製charge
1. CTRL+V(VISUAL BLOCK)
2. Shift + g : 跳至最下方
3. Shift +  $ : 跳行尾 
4. y ($複製)
5. p (貼上)
gmx_mpi trjconv -f heating.trr -pbc whole -s heating.tpr -o heating_1ns.gro -center -b 2300 -e 2300
### 取代指令
1. :%s/\<11\>/C/c

### 產生報錯
```=python
vi SCRIPTS/gromacs/run.py
```
![](https://hackmd.io/_uploads/H1meS_T03.png)
self.output改為 ''
gmx_d pdb2gmx -f /home/fangting/Work/MATERIALS_LIBRARY/LIBRARY/AP-core_qforce/STRUCTURE/AP-core_qforce.pdb -ff AP-core_qforce -p AP-core_qforce.top -o AP-core_qforce.gro


gmx_d pdb2gmx -f /home/fangting/Work/MATERIALS_LIBRARY/LIBRARY/1BA/STRUCTURE/1BA.pdb -ff 1BA -p 1BA.top -o 1BA.gro

gmx_mpi grompp -v  -f em.mdp -c 1BA.gro -o em.tpr -p 1BA.top -maxwarn 5
gmx_mpi grompp -v  -f nvt.mdp -c em.gro -o nvt.tpr -p AR3_10para.top -maxwarn 5