編集記録
[24.04.14] 仮想環境env1の構築について修正

　Amberのtutrialに沿って、4-hydroxy-Proline (PR4) のトポロジーファイルを作成します。その際、Gaussianの代替として無料で利用可能なPsi4でRESP電荷を求めてみます。

　タイトルに関連して、@tacomaさん、@Ag_smithさんの記事がとても参考になりました。Gaussianが利用可能であればそちらを参考にしてください。

Step0 : 環境構築

　環境構築にはminicondaを使いました。ただ、PCによっては失敗するようです。env1はStep2,3で、env2はStep4,5で使います。[24.04.14] env1の構築でrdkitのインストールに失敗することが多かったため、仮想環境内のpipを使ってみました（方法2）。

env1（方法1）

# rdkit 
conda install -c conda-forge rdkit python=3.10 -y 
# psi4 1.7
conda install -c psi4 psi4=1.7 -y
# resp
conda install -c psi4 resp -y
# openbabel
conda install -c conda-forge openbabel -y

env1（方法2）

# python 
conda install python -y 
# psi4
conda install -c conda-forge psi4 -y
# resp
conda install -c conda-forge resp -y
# openbabel
conda install -c conda-forge openbabel -y
# rdkit (pip)
pip install rdkit

env2

# Ambertools
conda install -c conda-forge ambertools -y

Step1 : 2つのコンフォメーションを作成

　ACE-PR4-NMEの.xyzファイルを異なるコンフォメーションで2つ作成します。PyMOLを使うと簡単に作成できます。

　pro4a.xyzは二面角N-CA-C‐Nが-60°のコンフォメーションです。

pro4a.xyz

27
pro4a
H          0.58800        1.62700        1.97400
C          0.63100        2.24000        1.08300
H          1.62100        2.13300        0.65300
H          0.48400        3.27700        1.34700
C         -0.42500        1.87100        0.06400
O         -1.02300        2.71600       -0.54700
N         -0.66900        0.54900       -0.15300
C         -1.81000        0.15300       -0.97500
H         -1.55200        0.14400       -2.03000
H         -2.62500        0.85000       -0.84400
C         -2.13600       -1.25000       -0.47500
H         -2.62400       -1.85100       -1.23600
O         -2.90200       -1.22000        0.70000
C         -0.76700       -1.79600       -0.09600
H         -0.23900       -2.11200       -0.98900
H         -0.83200       -2.63300        0.58400
C         -0.06500       -0.58400        0.54300
H         -0.28300       -0.55400        1.60200
C          1.45400       -0.68900        0.43100
O          2.10400       -1.08600        1.35800
N          1.98700       -0.39800       -0.77400
H          1.41200        0.05800       -1.44500
C          3.41500       -0.45400       -1.00200
H          3.93900        0.35000       -0.49500
H          3.80700       -1.39400       -0.64200
H          3.60000       -0.37900       -2.06500
H         -3.73500       -0.80000        0.53600

　pro4b.xyzは二面角N-CA-C‐Nが120°のコンフォメーションです。

pro4b.xyz

27
pro4b
H          0.23300        2.93400        0.49600
C         -0.82500        2.78400        0.67800
H         -1.36000        3.69700        0.46200
H         -0.96300        2.54700        1.72900
C         -1.40600        1.69400       -0.20300
O         -2.43700        1.87500       -0.78700
N         -0.68900        0.53600       -0.33900
C         -1.34400       -0.61900       -0.96800
H         -0.90100       -0.80800       -1.93500
H         -2.39300       -0.40300       -1.09300
C         -1.13100       -1.80700       -0.00400
H         -0.53900       -2.57400       -0.48400
O         -2.33400       -2.35700        0.47300
C         -0.37900       -1.21500        1.19100
H          0.32300       -1.90900        1.63600
H         -1.11100       -0.93800        1.94100
C          0.28000        0.05500        0.64100
H          0.45100        0.78700        1.41500
C          1.59300       -0.29700       -0.06100
O          1.64400       -1.15500       -0.90200
N          2.68200        0.38500        0.34400
H          2.57200        1.14300        0.97600
C          3.97800        0.17900       -0.27000
H          4.19700       -0.87800       -0.31600
H          4.01200        0.58100       -1.27600
H          4.73000        0.66900        0.33300
H         -2.76100       -2.83900       -0.22100

　原子の順番が入れ替わらないように注意してください。 Step3でRESP電荷の値がおかしくなります。PyMOLであれば一方のコンフォメーションを作成した後、BuilderのFix/Sculptを使ってもう一方を作成するとうまくいきます。

Step2 : 部分構造最適化

　作成した2つのコンフォメーションをHF/6-31G*で構造最適化します。その際、コンフォメーションを変えないために2つの二面角 (N-CA-C‐N、C-N-CA-C) を固定します。

　.inファイルを下記のpro4a.inのように作成し、部分構造最適化します。

pro4a.in

memory 4 GB

molecule Mol {
0 1
H          0.58800        1.62700        1.97400
C          0.63100        2.24000        1.08300
H          1.62100        2.13300        0.65300
H          0.48400        3.27700        1.34700
C         -0.42500        1.87100        0.06400
O         -1.02300        2.71600       -0.54700
N         -0.66900        0.54900       -0.15300
C         -1.81000        0.15300       -0.97500
H         -1.55200        0.14400       -2.03000
H         -2.62500        0.85000       -0.84400
C         -2.13600       -1.25000       -0.47500
H         -2.62400       -1.85100       -1.23600
O         -2.90200       -1.22000        0.70000
C         -0.76700       -1.79600       -0.09600
H         -0.23900       -2.11200       -0.98900
H         -0.83200       -2.63300        0.58400
C         -0.06500       -0.58400        0.54300
H         -0.28300       -0.55400        1.60200
C          1.45400       -0.68900        0.43100
O          2.10400       -1.08600        1.35800
N          1.98700       -0.39800       -0.77400
H          1.41200        0.05800       -1.44500
C          3.41500       -0.45400       -1.00200
H          3.93900        0.35000       -0.49500
H          3.80700       -1.39400       -0.64200
H          3.60000       -0.37900       -2.06500
H         -3.73500       -0.80000        0.53600
}

set {
   basis 6-31G*
   opt_type min
   geom_maxiter 300
   G_CONVERGENCE GAU
}

set optking {
  frozen_dihedral = ("
    21 19 17 7
    19 17 7 5
  ")
}

E = optimize('hf')
print(str(E))

psi4 -i pro4a.in -o pro4a.out -n 2

　pro4b.xyzの構造最適化ではpro4a.inの座標部分だけを変更します。終了後は.outファイルから最適化後の構造を抜き出し、新たな.xyzファイル (pro4a_opt.xyz, pro4b_opt.xyz) として保存します。PyMOLなどで最適化後の構造がおかしな構造になっていないかの確認をしておきましょう。

　memory 4 GB と -n 2は使用するPCに合わせて変更してください。他のアミノ酸などで構造最適化が収束しなかった場合 (もしくはエラーが出た場合) は、オプションにintrafrag_step_limit 0.1やopt_coordinates bothを追加するとうまくいくことがあります。

Step3 : RESP電荷の計算

　psi4のRESP電荷計算では、"resp"というplug-inを使用します。

元論文：https://onlinelibrary.wiley.com/doi/10.1002/qua.26035
GitHub：https://github.com/cdsgroup/resp

　大部分はrespのexamplesを参考に、.mol2ファイルへのRESP電荷の付与はMaking it rainを参考に以下のRESP.pyを作成しました。

RESP.py

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

args = sys.argv

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

Mol1, xyz1 = xyz2Mol(args[1], int(args[3]), int(args[4]))
Mol2, xyz2 = xyz2Mol(args[2], int(args[3]), int(args[4]))

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

mol1 = psi4.geometry(xyz1)
mol1.update_geometry()
mol1.set_name('conformer1')

mol2 = psi4.geometry(xyz2)
mol2.update_geometry()
mol2.set_name('conformer2')

molecules = [mol1, mol2]

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,
           'RESTRAINT'         : True,
           'IHFREE'            : False,
           'WEIGHT'            : [1, 1],
           }

# Call for first stage fit
charges1 = resp.resp(molecules, options)
print("\nStage One:")
print('RESP Charges')
print(charges1[1])

# Change the value of the RESP parameter A
options = {}
options['RESP_A'] = 0.001
resp.set_stage2_constraint(molecules[0], charges1[1], options)


# Add constraint for atoms fixed in second stage fit
options['grid'] = []
options['esp'] = []
for mol in range(len(molecules)):
    options['grid'].append('%i_%s_grid.dat' %(mol+1, molecules[mol].name()))
    options['esp'].append('%i_%s_grid_esp.dat' %(mol+1, molecules[mol].name()))

# Call for second stage fit
charges2 = resp.resp(molecules, options)

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

# Save coords to xyz file 
psi4out_xyz = name + '_resp.xyz'
mol1.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"
outfile_pdb = name + ".pdb"
print("Finished. Saved compound with partial charges as mol2 file: %s" % outfile_mol2)
obConversion.WriteFile(mol, outfile_mol2)

　構造最適化の時と同様に、set_num_threads(nthread=2) と set_memory(4GB)は使用するPCに合わせて変更してください。

　2つのコンフォメーションを考慮してRESP電荷を計算します。

# 引数は conf.1、conf.2、電荷、スピン多重度 です。
python RESP.py pro4b_opt.xyz pro4a_opt.xyz 0 1

　pro4b_opt_resp.mol2が出力されます。座標は先に指定したコンフォメーション (conf.1) のものを使用しています。

Step4 : .prepiファイルの作成

　pro4b_opt_resp.mol2をAmberToolsのantechamberで.acファイルに変換します。

antechamber -fi mol2 \
            -i pro4b_opt_resp.mol2 \
            -fo ac \
            -o pro4.ac \
            -at amber \
            -pf y

　-atについて、ligandではgaffやgaff2ですが、アミノ酸ではff14SBやff19SBなどを適用するためにamberです。-pf yで中間ファイルを削除しています。

　続いて、prepgenでACE-PR4-NMEのアミノ酸部分を指定するために.mcファイルを作成します。AtomNameはpro4.acと一致させます。

pro4.mc

HEAD_NAME N
TAIL_NAME C6 
MAIN_CHAIN C5
OMIT_NAME C
OMIT_NAME C1
OMIT_NAME O
OMIT_NAME H
OMIT_NAME H1
OMIT_NAME H2
OMIT_NAME N1
OMIT_NAME H9
OMIT_NAME C7
OMIT_NAME H10
OMIT_NAME H11
OMIT_NAME H12
PRE_HEAD_TYPE C
POST_TAIL_TYPE N
CHARGE 0.0

　最後に、AmberToolsのprepgenで.prepiファイルを作成します。

prepgen -i pro4.ac -o pro4.prepi -f prepi -m pro4.mc -rn PR4

Step5 : tleapで確認

　作成した.prepiファイルをleapで読み込めるかの確認をします。leap.inを作成し、tleapを実行します。

leap.in

source leaprc.protein.ff14SB
loadamberprep pro4.prepi
aa =sequence {ALA PR4 ALA}
savePDB aa ala-pr4-ala.pdb
saveamberparm aa ala-pr4-ala.prmtop ala-pr4-ala.prmcrd
quit

tleap -f leap.in

　問題無ければala-pr4-ala.pdb、ala-pr4-ala.prmtop、ala-pr4-ala.prmcrdが出力されます。

おまけ

　RESP電荷の値をPsi4で求めた場合とGaussianで求めた場合とで比較してみます。Psi4で求めた場合、Gaussianで求めた場合よりも電荷の絶対値が小さくなるようです。