↑のチュートリアルを行って修飾アミノ酸である「4-Hydoroxyl-Proline」のトポロジーの残基ファイルを作成する。ただし、チュートリアル通りに行かないところもあったので自分なりの方法を書きます。
なお、通常のProlineとの差は下記画像参照。
実行環境
- Ambertools22
- Ubuntu 20.04LTS
ステップ1 mol2ファイルの用意
ACE-PR4-ACEのα-helix conformationとβ-sheet conformationの2つのmol2ファイルを用意する（チュートリアルのページから落とせます）。
2つのconformationを使うことで、電荷をもっともらしいものにする。
αとβの構造作成にはpymolなどを使えば良いと思います。
ステップ2 Gaussian入力ファイルの作成
以下のantechamberコマンドで、alphaとbetaの2つのGaussianの入力ファイルを生成する。
antechamber -fi mol2 -fo gzmat -i pro4a.mol2 -o pro4a.gau
antechamber -fi mol2 -fo gzmat -i pro4b.mol2 -o pro4b.gau
今回はαとβのconformationは変えたくないためCα周りの二面角（いわゆるΦとψ）を固定して部分最適化を行います。pro4a.gauとpro4b.gauを修正しましょう。
チュートリアルでは
- キーワードセクションのoptをpoptに変更(つまり部分最適化を行う)
- 最小化中のt19とt21を "F "タグでフリーズします。
と書いてあります。
つまり、ルートセクションを
#HF/6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) popt
とし、
t19= -60.0944 F
b20= 1.1998
a20= 120.3382
t20= 144.1795
b21= 1.3494
a21= 116.3387
t21= -39.9570 F
とします。これで二面角が固定されます。
個人的には
opt=ModRedundunt
として
19 17 7 5 F
21 19 17 7 F
としても良いと思います。
ステップ3 Gaussianの実行
gaussianを実行し、alphaとbetaの2つの出力ファイルを得る。
計算環境がない方はWebMOを使うのも手かと。
ステップ4 ESPファイルを取り出しRESP電荷を計算
gaussianの出力ファイル（alphaとbeta）からESPファイルを抽出し、catコマンドでマージする。
$ espgen -i pro4a.out -o pro4a.esp
$ espgen -i pro4b.out -o pro4b.esp
$ cat pro4a.esp pro4b.esp >pro4.esp
続いてrespgenを実行するため、mol2ファイルからacファイルを生成します。これはrespgenがacファイルの入力しか受け付けないからです。
一応-at amberとamberの原子タイプを指定していますが、ここでは多分不要です。
antechamber -fi mol2 -i pro4a.mol2 -fo ac -o pro4.ac -at amber
respgenを実行し、2段階のrespフィットの入力ファイルを準備します。
$ respgen -i pro4.ac -o pro4-step1.respin -f resp1 -n 2
$ respgen -i pro4.ac -o pro4-step2.respin -f resp2 -n 2
respコマンドでRESP電荷を算出し、antechamberでRESP電荷を付与したacファイルを生成します。
$ resp -O -i pro4-step1.respin -o pro4.respout1 -e pro4.esp -t qout_stage1
$ resp -O -i pro4-step2.respin -o pro4.respout2 -e pro4.esp -q qout_stage1 -t qout_stage2
$ antechamber -i pro4.ac -fi ac -at amber -o pro4_resp.ac -fo ac -c rc -cf qout_stage2
※もとのチュートリアルではrespgenのインプットファイルをmol2にしていたが、うまくいかなかったためacファイルに切り替えた。また、２つの構造を扱うためrespgenのオプション-n 2を追加した。また、crgファイルを作らず、qout_stage2の電荷をそのままacファイルに反映させている。
なお、1つのコンフォメーションしか使わない場合は、charge methodタグをrespに設定したgaussian outputファイルから直接acファイルを生成することも可能です。
ステップ5 prepファイルの作成
pro4_resp.acとmainchain.pro4を読み込んで、prepgenコマンドでprepファイルを作成します。
mainchain.pro4は以下↓
HEAD_NAME N1
TAIL_NAME C7
MAIN_CHAIN C6
OMIT_NAME C1
OMIT_NAME C2
OMIT_NAME O1
OMIT_NAME H1
OMIT_NAME H2
OMIT_NAME H3
OMIT_NAME N2
OMIT_NAME H10
OMIT_NAME C8
OMIT_NAME H11
OMIT_NAME H12
OMIT_NAME H13
PRE_HEAD_TYPE C
POST_TAIL_TYPE N
CHARGE 0.0
mainchainファイルでは、「HEAD_NAME」は直前の残基と接続する原子、「TAIL_NAME」は直後の残基と接続する原子。「OMIT_NAME」はキャップ原子である。MAINCHAINは主鎖原子（HEAD_NAMEとTAIL_NAMEは自動的に主鎖原子）であり、PRE_HEAD_TYPEとPOST_TAIL_TYPEはそれぞれヘッドとテールに接続されている他の残基の原子の原子タイプである。CHARGE は残基の正味の電荷である。
prepgenを実行してprepファイルを生成する。
$ prepgen -i pro4_resp.ac -o pro4.prepi -f prepi -m mainchain.pro4 -rn PR4 -rf pro4.res
PRE_HEAD_TYPE is C
POST_TAIL_TYPE is N
Net charge of truncated molecule is 0.00
HEAD_ATOM 7 N1
TAIL_ATOM 19 C7
MAIN_CHAIN 1 7 N1
MAIN_CHAIN 2 17 C6
MAIN_CHAIN 3 19 C7
OMIT_ATOM 1 2 C1
OMIT_ATOM 2 5 C2
OMIT_ATOM 3 6 O1
OMIT_ATOM 4 1 H1
OMIT_ATOM 5 3 H2
OMIT_ATOM 6 4 H3
OMIT_ATOM 7 21 N2
OMIT_ATOM 8 22 H10
OMIT_ATOM 9 23 C8
OMIT_ATOM 10 24 H11
OMIT_ATOM 11 25 H12
OMIT_ATOM 12 26 H13
Number of mainchain atoms (including head and tail atom): 3
Number of omited atoms: 12
ステップ6 トポロジーファイル、座標ファイルの作成
pro4.prepiを読み込めばleapを使用してトポロジーファイルを生成することができます。ここでは、ALA-PRO4-ALAのトポロジーファイルをleapで生成してみます。
tleap.in
source leaprc.protein.ff14SB
loadamberprep pro4.prepi
aa =sequence {ALA PR4 ALA}
saveamberparm aa ala-pr4-ala.prmtop ala-pr4-ala.inpcrd
tleapを実行しましょう。
$ tleap -f tleap.in
うまくいきましか？ala-pr4-ala.prmtopとala-pr4-ala.inpcrdが生成されていればこれで完了です。
なお、"prepgen"ではなく"residuegen"コマンドを使った方法はmainchainの指定が間違っていてダメでした。誰か解決法わかる方がいたら教えてください。
参考
「Amberによる生体高分子シミュレーション入門」