Psi4/AmberToolsで作成したトポロジーファイルをつかいMDシミュレーション

　以前投稿した記事ではPsi4/AmberToolsを使ったトポロジーファイルの作成方法を紹介しました。本記事では、低分子のトポロジーファイルを同方法で作成し、MDシミュレーションで動かす直前までを解説します。なお本記事の内容を応用すると、基質‐酵素複合体のMDシミュレーションもできるようになります。

 https://qiita.com/wai_yuito/items/872a4db3c56426251064

　今回は酢酸の2量体を扱います。初期構造はベンゼンをPyMOLで編集して作ってみました（図1）。これからこの`2ACA.pdb`を編集していきます。なお`CONECT`の行は不要なので削除しておきました。

 ```:2ACA.pdb
HETATM    1  C01 UNK     1       2.001   0.067  -0.002  0.00  0.00           C  
HETATM    2  C04 UNK     1      -2.008  -0.068   0.014  0.00  0.00           C  
HETATM    3  C07 UNK     1      -3.546   0.002  -0.001  0.00  0.00           C  
HETATM    4  C08 UNK     1       3.539  -0.006   0.002  0.00  0.00           C  
HETATM    5  H01 UNK     1      -3.920   0.018   1.022  0.00  0.00           H  
HETATM    6  H02 UNK     1      -3.945  -0.869  -0.521  0.00  0.00           H  
HETATM    7  H03 UNK     1      -3.863   0.909  -0.517  0.00  0.00           H  
HETATM    8  H04 UNK     1       3.906  -0.022  -1.024  0.00  0.00           H  
HETATM    9  H05 UNK     1       3.860  -0.911   0.517  0.00  0.00           H  
HETATM   10  H06 UNK     1       3.942   0.866   0.518  0.00  0.00           H  
HETATM   11  H07 UNK     1       1.022   1.520  -0.903  0.00  0.00           H  
HETATM   12  H08 UNK     1      -1.027  -1.527  -0.867  0.00  0.00           H  
HETATM   13  O01 UNK     1      -1.269   1.149   0.034  0.00  0.00           O  
HETATM   14  O02 UNK     1      -1.383  -1.348   0.007  0.00  0.00           O  
HETATM   15  O03 UNK     1       1.217  -1.129   0.007  0.00  0.00           O  
HETATM   16  O04 UNK     1       1.334   1.331  -0.015  0.00  0.00           O  
END
 ```

図1. 作成した酢酸二量体
![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/ccfb2263-cab4-8d01-049c-b48840bcae99.png)

# 0. 環境構築
　以前作成した仮想環境[env1](https://qiita.com/wai_yuito/items/872a4db3c56426251064)にbiopythonをinstallしておきます。
```bash
conda activate env1
conda install -c bioconda biopython
```

# 1. resn・resiを変更する
　`2ACA.pdb`では酢酸の残基名（residue name；resn）がUNKとなっています。またどちらの酢酸も残基番号（residue number；resi）が1であり区別されていません。この2点をテキストエディタで直接編集することもできますが、今回はPyMOLで編集します。
 
1. PyMOLで`2ACA.pdb`を開く
2. `Mouse → Selection Mode → Atoms`で原子を一つずつ選択できるようにする
3. 一方の酢酸の原子を全てクリックして選択する
4. Command Input Areaに`alter sele, resn='ACA'`と入力しEnter
5. Command Input Areaに`alter sele, resi=1`と入力しEnter
6. 何もない所をクリックして選択を解除する
7. もう一方の酢酸の原子を全てクリックして選択する
8. Command Input Areaに`alter sele, resn='ACA'`と入力しEnter
9. Command Input Areaに`alter sele, resi=2`と入力しEnter
10. `File → Export Molecule… → Save…`によりファイル名`2ACA_res.pdb`として保存する

 ```:2ACA_res.pdb
HETATM    1  O01 ACA     1      -1.269   1.149   0.034  0.00  0.00           O  
HETATM    2  O02 ACA     1      -1.383  -1.348   0.007  0.00  0.00           O  
HETATM    3  C04 ACA     1      -2.008  -0.068   0.014  0.00  0.00           C  
HETATM    4  C07 ACA     1      -3.546   0.002  -0.001  0.00  0.00           C  
HETATM    5  H01 ACA     1      -3.920   0.018   1.022  0.00  0.00           H  
HETATM    6  H02 ACA     1      -3.945  -0.869  -0.521  0.00  0.00           H  
HETATM    7  H03 ACA     1      -3.863   0.909  -0.517  0.00  0.00           H  
HETATM    8  H08 ACA     1      -1.027  -1.527  -0.867  0.00  0.00           H  
HETATM    9  C01 ACA     2       2.001   0.067  -0.002  0.00  0.00           C  
HETATM   10  O03 ACA     2       1.217  -1.129   0.007  0.00  0.00           O  
HETATM   11  O04 ACA     2       1.334   1.331  -0.015  0.00  0.00           O  
HETATM   12  C08 ACA     2       3.539  -0.006   0.002  0.00  0.00           C  
HETATM   13  H04 ACA     2       3.906  -0.022  -1.024  0.00  0.00           H  
HETATM   14  H05 ACA     2       3.860  -0.911   0.517  0.00  0.00           H  
HETATM   15  H06 ACA     2       3.942   0.866   0.518  0.00  0.00           H  
HETATM   16  H07 ACA     2       1.022   1.520  -0.903  0.00  0.00           H  
END
```

# 2. Psi4でRESP電荷を求める
　初めに酢酸のXYZファイルを入手します。PyMOLで`2ACA_res.pdb`を開き一方の酢酸を`extract object`で切り抜いたのち、`Export Molecule…`で`ACA.xyz`として保存します。

　次に`ACA.xyz`を構造最適化します。今回はB3LYP/6-31+g(d,p)で計算します。inputファイル（`ACA_opt.in`）を作成し、仮想環境env1で計算します。
```:ACA_opt.in
memory 4 GB

molecule Mol {
0 1
C 2.001000 0.067000 -0.002000
O 1.217000 -1.129000 0.007000
O 1.334000 1.331000 -0.015000
C 3.539000 -0.006000 0.002000
H 3.906000 -0.022000 -1.024000
H 3.860000 -0.911000 0.517000
H 3.942000 0.866000 0.518000
H 1.022000 1.520000 -0.903000
}

set {
   basis 6-31+G**
   opt_type min
   geom_maxiter 300
   g_convergence gau
}

E = optimize('b3lyp')
print(str(E))
```
```bash
# psi4で構造最適化
psi4 -i ACA_opt.in -o ACA_opt.out -n 4
```
計算終了後のlogファイル（`ACA_opt.out`）から、構造最適化後の構造をコピーし`ACA_opt.xyz`を作成します。

　続いて構造最適化後の酢酸についてRESP電荷を求めます。githubで公開されているスクリプトの例を参考に、以下のスクリプトを作成しました。
```python:ACA_resp.py
constraint_group_list =[[5, 6, 7]]
constraint_charge_list = [1, 2, 3, 8]

from openbabel import openbabel as ob
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import rdDetermineBonds
import os
import sys
import psi4
import resp

args = sys.argv

def cleanUp(psi4out_xyz):
    deleteTheseFiles = ['1_default_grid.dat','1_default_grid_esp.dat','grid.dat','timer.dat']
    deleteTheseFiles.append(psi4out_xyz)
    for fileName in deleteTheseFiles:
        if os.path.exists(fileName):
            os.remove(fileName)

def xyz2Mol(file_name, c, m):
    Mol = Chem.MolFromXYZFile(file_name)
    rdDetermineBonds.DetermineBonds(Mol, charge = c)
    xyz = Chem.MolToXYZBlock(Mol)
    xyz = '\n'.join(xyz.split('\n')[2:])
    c_m = '%s %s\n' %(c, m)
    xyz = c_m + xyz
    return Mol, xyz

obConversion = ob.OBConversion()
obConversion.SetInAndOutFormats("xyz", "mol2")

name = args[1].rstrip('.xyz')

Mol, xyz = xyz2Mol(args[1], int(args[2]), int(args[3]))

psi4.set_num_threads(nthread=4)
psi4.set_memory('4GB')

mol = psi4.geometry(xyz)
mol.update_geometry()


# https://github.com/cdsgroup/resp/blob/master/examples/example1.py

options = {'VDW_SCALE_FACTORS' : [1.4, 1.6, 1.8, 2.0],
           'VDW_POINT_DENSITY' : 1.0,
           'RESP_A'            : 0.0005,
           'RESP_B'            : 0.1,
           }

# Call for first stage fit
charges1 = resp.resp([mol], options)
print('Electrostatic Potential Charges')
print(charges1[0])
print('Restrained Electrostatic Potential Charges')
print(charges1[1])

# Change the value of the RESP parameter A
options = {}
options['RESP_A'] = 0.001

# Add constraint for atoms fixed in second stage fit
constraint_charge = []
for i in constraint_charge_list:
    constraint_charge.append([charges1[1][i-1], [i]])
options['constraint_charge'] = constraint_charge
options['constraint_group'] = constraint_group_list
options['grid'] = ['1_%s_grid.dat' %mol.name()]
options['esp'] = ['1_%s_grid_esp.dat' %mol.name()]

# Call for second stage fit
charges2 = resp.resp([mol], options)

# Get RESP charges
print("\nStage Two:")
print('RESP Charges')
print(charges2[1])


#https://colab.research.google.com/github/pablo-arantes/making-it-rain/blob/main/Partial_Charges.ipynb

### save coords to xyz file 
psi4out_xyz = name + '_RESP.xyz'
mol.save_xyz_file(psi4out_xyz,1)

### read xyz file and write as mol2 
mol = ob.OBMol()
obConversion.ReadFile(mol, psi4out_xyz)

### write as mol2 
outfile_mol2 = name + "_resp.mol2"
obConversion.WriteFile(mol, outfile_mol2)

### set new partial charges
count = 0
newChg_temp = charges2[1]
for atom in ob.OBMolAtomIter(mol):
    newChg = newChg_temp[count]
    atom.SetPartialCharge(newChg)
    count += 1

### write as mol2 
outfile_mol2 = name + "_resp.mol2"
print("Finished. Saved compound with partial charges as mol2 file: %s" % outfile_mol2)
obConversion.WriteFile(mol, outfile_mol2)

## clean up
cleanUp(psi4out_xyz)
```
スクリプトの上から2行では、酢酸において等価な原子などを指定しています（図2）。XYZファイルに対応した原子番号は`酢酸をクリックして選択 → 右クリック → label → atom identifiers → index`で表示できます。

図2. constraint_group_listとconstraint_charge_list
![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/e04d4bb6-5289-0209-5172-8cae75b2d7b9.png)


```bash
# psi4でRESP電荷を計算
python RESP.py ACA_opt.xyz 0 1
```
引数1は最適化後の構造、引数2は電荷、引数3はスピン多重度です。計算が正常に終了すると`ACA_opt_resp.mol2`が出力されます。PyMOLで割り当てられたRESP電荷を表示し、異常がないことを確認します（図3）。電荷は`酢酸をクリックして選択 → 右クリック → label → other properties → partial charge`で表示できます。

図3. ACA_opt_resp.mol2のRESP電荷を確認
![fig3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/c29014dc-d22d-b35d-af83-62d4df8e04ee.png)


# 3. AmberToolsでトポロジーファイルを作成する
　AmberToolsのantechamberで`ACA_opt_resp.mol2`を変換し、トポロジーファイルを作成します。以下を仮想環境[env2](https://qiita.com/wai_yuito/items/872a4db3c56426251064)で実行します。
```bash
# .mol2形式 → .ac形式
antechamber -fi mol2 -i pro4b_opt_resp.mol2 \
            -fo ac -o aca.ac \
            -at gaff -rm ACA -pf y

# .ac形式 → .prepi形式
antechamber -fi ac -i aca.ac \
            -fo prepi -o aca.prepi \
            -pf y

# .ac形式 → .pdb形式
antechamber -fi ac -i aca.ac \
            -fo pdb -o aca_template.pdb \
            -pf y
```
`aca.prepi`がトポロジーファイルです。`aca_template.pdb`は次の工程で`2ACA_res.pdb`のatom nameを変更するために使います。

# 4. atom nameを変更する
　作成したトポロジーファイルを`2ACA_res.pdb`に適用するには、`2ACA_res.pdb`中の酢酸のatom nameが`aca.prepi`で定めたものと一致している必要があります（図4）。

図4. 2ACA_res.pdbとaca.prepi
![fig4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/0aef8f24-aebe-9362-ead1-ea20b72024ce.png)


　初めにPyMOLやテキストエディタを使い、`2ACA_res.pdb`中の2つの酢酸をそれぞれ別のPDBファイル（`ACA_1.pdb`、`ACA_2.pdb`）として保存します。

　次にantechamberで出力した`aca_template.pdb`（正しいatom nameをもつ）に`ACA_1.pdb`、`ACA_2.pdb`のatom nameをそろえます。RDkitで共通構造を抽出できることを利用し、次のスクリプトを作成しました。これを仮想環境env1で実行します。
```python:edit_an.py
import sys
from rdkit import Chem
from rdkit.Chem import rdFMCS
from Bio.PDB import *

args = sys.argv

pdb1 = args[1]
pdb2 = args[2]
newpdb = pdb1.split('.pdb')[0] + '_editan.pdb'

mol1 = Chem.MolFromPDBFile(pdb1, removeHs=False)
mol2 = Chem.MolFromPDBFile(pdb2, removeHs=False)

mols = [mol1, mol2]

res = rdFMCS.FindMCS(mols, bondCompare=rdFMCS.BondCompare.CompareAny)
mcs = Chem.MolFromSmarts(res.smartsString)

mol1_match = mol1.GetSubstructMatch(mcs)
mol1_match_list = [x for x in mol1_match]

mol2_match = mol2.GetSubstructMatch(mcs)
mol2_match_list = [x for x in mol2_match]

pdb_parser = PDBParser()
structure = pdb_parser.get_structure('X', pdb1)
pdb1_res = structure.get_list()[0].get_list()[0].get_list()
pdb1_atoms = []

for res in pdb1_res:
    resname = res.get_resname()
    atoms = res.get_list()
    for atom in atoms:
        res_atom = [resname, atom.get_fullname(),  atom.get_coord()]
        pdb1_atoms.append(res_atom)

pdb_parser = PDBParser()
structure = pdb_parser.get_structure('X', pdb2)
pdb2_res = structure.get_list()[0].get_list()[0].get_list()
pdb2_atoms = []
ref_resn = []

for res in pdb2_res:
    resname = res.get_resname()
    ref_resn.append(resname)
    atoms = res.get_list()
    for atom in atoms:
        res_atom = [resname, atom.get_fullname(), atom.get_coord()]
        pdb2_atoms.append(res_atom)

f = open(pdb1, 'r')
f_l = f.readlines()
f.close()

i = 0

for n, m in zip(mol1_match_list, mol2_match_list):
    i = i + 1
    f_l[n] = f_l[n].replace(pdb1_atoms[n][1], pdb2_atoms[m][1])
    
with open(newpdb, 'w') as f:
    f.writelines(f_l)
    f.write('END')
```
```bash
python edit_an.py ACA_1.pdb aca_template.pdb
python edit_an.py ACA_2.pdb aca_template.pdb
```
引数1はatom nameを直したいPDB、引数2は正しいatom nameをもつPDBです。atom nameが`aca_template.pdb`とそろった`ACA_1_editan.pdb`、`ACA_2_editan.pdb`が出力されます。

　続いて`ACA_1_editan.pdb`、`ACA_2_editan.pdb`を一つのPDBファイルにまとめます。両者をPyMOLで開き、`Export Molecule…`から`2ACA_res_an.pdb`として保存します。ここまでで問題が無ければ、`2ACA_res_an.pdb`は次のようになるはずです。2つの酢酸の座標が入れ替わっていても、問題ありません。なお`CONECT`の行は削除しました。
```:2ACA_res_an.pdb
HETATM    1  C   ACA     1      -2.008  -0.068   0.014  0.00  0.00           C  
HETATM    2  O   ACA     1      -1.269   1.149   0.034  0.00  0.00           O  
HETATM    3  C1  ACA     1      -3.546   0.002  -0.001  0.00  0.00           C  
HETATM    4  O1  ACA     1      -1.383  -1.348   0.007  0.00  0.00           O  
HETATM    5  H   ACA     1      -3.920   0.018   1.022  0.00  0.00           H  
HETATM    6  H1  ACA     1      -3.945  -0.869  -0.521  0.00  0.00           H  
HETATM    7  H2  ACA     1      -3.863   0.909  -0.517  0.00  0.00           H  
HETATM    8  H3  ACA     1      -1.027  -1.527  -0.867  0.00  0.00           H  
HETATM    9  C   ACA     2       2.001   0.067  -0.002  0.00  0.00           C  
HETATM   10  O   ACA     2       1.217  -1.129   0.007  0.00  0.00           O  
HETATM   11  C1  ACA     2       3.539  -0.006   0.002  0.00  0.00           C  
HETATM   12  O1  ACA     2       1.334   1.331  -0.015  0.00  0.00           O  
HETATM   13  H   ACA     2       3.906  -0.022  -1.024  0.00  0.00           H  
HETATM   14  H1  ACA     2       3.860  -0.911   0.517  0.00  0.00           H  
HETATM   15  H2  ACA     2       3.942   0.866   0.518  0.00  0.00           H  
HETATM   16  H3  ACA     2       1.022   1.520  -0.903  0.00  0.00           H  
END
```
　最後に1つ目の最後の行と2つ目の最初の行との間に`TER`を挿入します。
```
HETATM    7  H2  ACA     1      -3.863   0.909  -0.517  0.00  0.00           H  
HETATM    8  H3  ACA     1      -1.027  -1.527  -0.867  0.00  0.00           H  
TER
HETATM    9  C   ACA     2       2.001   0.067  -0.002  0.00  0.00           C  
HETATM   10  O   ACA     2       1.217  -1.129   0.007  0.00  0.00           O  
```

# 5. LEaPに通す
　本記事の最後に、作成した`2ACA_res_an.pdb`がLEaPを通ることを確認します。`leap.in`を作成し仮想環境env2でtleapを実行します。
 ```:leap.in
loadamberprep ./aca.prepi
source leaprc.gaff
source leaprc.water.tip3p

sys = loadPdb 2ACA_res_editan.pdb

solvateBox sys TIP3PBOX 10.0

charge sys
check sys

saveAmberParm sys leap.parm7 leap.rst7
savepdb sys leap.pdb

quit
 ```
 ```bash
 tleap -f leap.in
 ```
 正常に終了すると`leap.parm7`、`leap.rst7`、`leap.pdb`が出力されます（図5）。

 図5. leap.pdb
![fig5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/15e5f457-eca2-0b5c-9dd8-cd04d1822727.png)