Go to Qiita Advent Calendar 2024 Top
LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@wai_yuito(山田 惟人)

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

Last updated at Posted at 2024-04-13

　以前投稿した記事では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. 作成した酢酸二量体
fig1.png

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で編集します。

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

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

# 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

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
fig4.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で実行します。

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
fig5.png

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?