# 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檔範例
```
%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
```
3. 用obabel把log 檔換成pdb檔
```
obabel -ilog file.log -O file.pdb
```
4. 建atp檔和rtp檔
`.atp範例
`
`.rtp範例
`
5. 打開molden，把pdb檔的residule name對應到rtp檔
`.pdb`
### 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
```=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
6. 改.gro檔 box大小
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.比較結果
## 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
```
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
Sign in with Wallet
Connect another wallet