Psi4/AmberToolsで非天然アミノ酸のトポロジーファイルを作成

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

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

https://ambermd.org/antechamber/pro4.html

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

https://qiita.com/tacoma/items/d7237f4afbb7f18ff60c

https://qiita.com/Ag_smith/items/430e9efb32a855d4c511

# Step0 : 環境構築
　環境構築にはminicondaを使いました。ただ、PCによっては失敗するようです。env1はStep2,3で、env2はStep4,5で使います。[24.04.14] env1の構築でrdkitのインストールに失敗することが多かったため、仮想環境内のpipを使ってみました（方法2）。
```bash: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
```
```bash: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
```
```bash:env2
# Ambertools
conda install -c conda-forge ambertools -y
```

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

　`pro4a.xyz`は二面角N-CA-C‐Nが-60°のコンフォメーションです。
```tex: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°のコンフォメーションです。
```tex: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`のように作成し、部分構造最適化します。
```tex: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))
```
```bash
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](https://github.com/cdsgroup/resp/tree/master/examples)を参考に、.mol2ファイルへのRESP電荷の付与は[Making it rain](https://pablo-arantes.github.io/making-it-rain/)を参考に以下の`RESP.py`を作成しました。

```python: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電荷を計算します。
```bash
# 引数は 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ファイルに変換します。
```bash
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`と一致させます。
```txt: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ファイルを作成します。
```bash
prepgen -i pro4.ac -o pro4.prepi -f prepi -m pro4.mc -rn PR4
```

# Step5 : tleapで確認
　作成した.prepiファイルをleapで読み込めるかの確認をします。`leap.in`を作成し、tleapを実行します。
```txt: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
```
```bash
tleap -f leap.in
```
　問題無ければ`ala-pr4-ala.pdb`、`ala-pr4-ala.prmtop`、`ala-pr4-ala.prmcrd`が出力されます。

# おまけ
　RESP電荷の値をPsi4で求めた場合と[Gaussianで求めた場合](https://ambermd.org/antechamber/pro4.prepi)とで比較してみます。Psi4で求めた場合、Gaussianで求めた場合よりも電荷の絶対値が小さくなるようです。
![画像1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3302733/cd4a9d07-1d8a-5a01-b557-b289c3a08309.png)