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

　今回は酢酸の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. 作成した酢酸二量体

0. 環境構築

　以前作成した仮想環境env1にbiopythonをinstallしておきます。

conda activate env1
conda install -c bioconda biopython

1. resn・resiを変更する

　2ACA.pdbでは酢酸の残基名（residue name；resn）がUNKとなっています。またどちらの酢酸も残基番号（residue number；resi）が1であり区別されていません。この2点をテキストエディタで直接編集することもできますが、今回はPyMOLで編集します。

PyMOLで2ACA.pdbを開く
Mouse → Selection Mode → Atomsで原子を一つずつ選択できるようにする
一方の酢酸の原子を全てクリックして選択する
Command Input Areaにalter sele, resn='ACA'と入力しEnter
Command Input Areaにalter sele, resi=1と入力しEnter
何もない所をクリックして選択を解除する
もう一方の酢酸の原子を全てクリックして選択する
Command Input Areaにalter sele, resn='ACA'と入力しEnter
Command Input Areaにalter sele, resi=2と入力しEnter
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))

# psi4で構造最適化
psi4 -i ACA_opt.in -o ACA_opt.out -n 4

計算終了後のlogファイル（ACA_opt.out）から、構造最適化後の構造をコピーしACA_opt.xyzを作成します。

　続いて構造最適化後の酢酸についてRESP電荷を求めます。githubで公開されているスクリプトの例を参考に、以下のスクリプトを作成しました。

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

# 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電荷を確認

3. AmberToolsでトポロジーファイルを作成する

　AmberToolsのantechamberでACA_opt_resp.mol2を変換し、トポロジーファイルを作成します。以下を仮想環境env2で実行します。

# .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

　初めに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で実行します。

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')

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

tleap -f leap.in

正常に終了するとleap.parm7、leap.rst7、leap.pdbが出力されます（図5）。

図5. leap.pdb