0. ONIOM法の概要

前回の記事で少し触れましたが、ONIOM法の理論の説明や実際にどういう流れでやるか、ということについてはこの論文やこのReviewなどを参考にすると良いです。え？わからない？まあとりあえず実際にやってみたあとでそれから内容を理解してもええんやで。ということでここでは飛ばします。
また、このブログ記事は本家TAOPackageのページの補足となっています。

ONIOMを用いた反応遷移状態の考察の大まかな流れですが、こんな感じ
1. 結晶構造を起点に少しだけ（～10 ns）MDをして、水や基質の配置をほぐす。
2. 蛋白質・基質から5 Aより外の水分子を除外（計算量を減らすため）
3. 基質と反応に重要な残基をコア領域として選択、またそこから7 A以内の残基の構造最適化を許可、それより外の残基を動かないよう設定。
4. 構造最適化開始
5. 最適化した構造について、より高い計算レベルで系のエネルギーのシングルポイント計算をする。
6. 1~5を反応前、反応後、遷移状態について行い、各エネルギー値を比較する。

1. 今回使用する蛋白質とリガンド分子

今回、蛋白質はPDB: 1iuwのp-hydroxybenzoate hydroxylase、リガンドに4-hydroxybenzoate (p-OHB)を用います。
この蛋白質は水酸基-OHをp-OHBの3位炭素に対して付加させる反応です。
この反応は、水酸基OH+が移動する芳香族求電子置換反応の中でも珍しい例で、補酵素FADの過酸化状態FADHOOHのペルオキシドOH+が求電子剤となって芳香族に攻撃するパティーンです。
詳しい反応機構はこのReviewのScheme 2でも見てください。

今回扱うのは上記のReviewのk6ステップ（水酸化ステップ）で、ここをONIOM法でエネルギーを求める方法を紹介します。

2. 環境

以下のソフトウェアをインストールしていることが条件です。具体的には

• Linux マシンまたはMac マシン
• AmberTools 16をインストール
• GROMACS 2016.3などでMDを動かす方法を習得している
• VMD 1.9.2（MacOS X OpenGL (32-bit Intel x86)でよい）
• Gaussian 09が使える環境にある
• GaussViewを使えることが強く推奨される
• 前回の記事を読んでTaopackageをインストールしてある

という感じですかね。
ONIOM法自体にはGaussian 09さえあれば良いのですが、その前準備としてGROMACSなどのMDソフトウェア、そしてtaopackageとGaussian 09のAmber力場に対して相性が良いAmberTools 16を入れておきます。

3. MDで10 nsほど動かす

まず蛋白質構造を結晶構造から用意します。PDBの1iuwをダウンロードしてきて、補酵素FADHOOHと基質PHAの電荷を設定して、solvateでちょうどいい感じに水分子とイオンを撒いて、AmberTools 16のLEaPに通すための準備をして……っていう感じです。

が、ここではそれが主題ではないので、それらの前処理が終わったファイルを置いておくのでダウンロードして持って行ってください。

中の解説をしますと、まずPDBファイル1iuw_ligands.pdbには、補酵素FADHOOH（残基FAP A 395）と基質p-OHB（残基PHA A 396）が書かれています。FADHOOHは結晶構造中にはペルオキシド状態になっておりませんが、この後の設定で適切な位置に過酸化水素OOHを生やすことができます。
PHAの方は、水酸化反応に必須な条件を満たすため、フェノール基が脱水素した状態になっています。

leap.watディレクトリにFADHOOH, PHAに対応した残基パラメータファイルperoxide2.prep, anion_pohb.prepが入っており、MDシミュレーションに必要な電荷の情報はこれに含まれております。

以上を確認したら、leap.watディレクトリでコマンドを実行し、AmberTools 16のLEaPを動作させれば、AmberTools用のトポロジーファイルと初期座標ファイルが生成されます。後は普通のMDプロトコルに従ってGROMACSなりなんなりを使って10nsほどNPTシミュレーションを動かしておきます。

4. 5 Aより遠方の水分子を除外する

次に、MDが終わった後の蛋白質について、蛋白質＋リガンドをセンタリングした状態でPDBファイル形式に出力します。上記のダウンロード素材集ではこれを1iuw_10ns.pdbとして用意しました。まずはこれをPyMOLかなにかで見てください。蛋白質の周りに水が存在し溶媒ボックスいっぱいに詰まっていることがわかります。
計算で全部の水分子を入れてやるのは大変なので、蛋白質から5 Aより外の水分子はここで除外しておきます。
色んなやり方がありますが、VMDのコマンドを使ってターミナル上でパパッと自動化してやる方法を紹介します。

ターミナル上でVMDにパスが通っていれば、vmdと打つことで対話モードのVMDが起動するはずです。

ここで、次のような流れでやれば、蛋白質から5 Aより外の水分子を蛋白質分子と一緒に出力してくれるはずです。
ここで、XXX.pdb, 1, 999, chain Xの部分には適切な値をいれてやる必要があります。例えばloadするPDBファイルが1iuw_10ns.pdbとして、このファイルには残基1から398までが蛋白質＋基質データとなっています。
しかしながら、先ほどの1iuw_10ns.pdbにはchain IDが割り振られていなかったり、水分子が大量にあるために残基ナンバーが重複していたりするので、これでは処理上不都合なことがたくさんあります。
そこでこの前処理として、蛋白質にchain Aを付加し、残りは残基ナンバーが重複しないように異なるchian IDと残基ナンバーを与えます。この処理は../scripts/chains.plで行います。ただし最後にONIOMにインプットさせるときにはchain IDがあるとエラーになるっぽいので、これを外しておきます。

処理としては

gmx editconf -f 1iuw_10ns.pdb -o temp1.pdb -resnr 0
 
perl ../scripts/onioms/chains.pl temp1.pdb > temp2.pdb
 
vmd <<EOF
 
mol load pdb temp2.pdb
 
waitfor all
 
set all [atomselect top "(resid 0 to 397 and chain A) or (same residue as resname WAT and within 5 of (resid 0 to 397 and chain A))"]
 
\$all writepdb temp3.pdb
 
quit
 
EOF
 
perl -nle 'substr($_,21,1) = " " if ($_ =~ /^ATOM/ or /^HETATM/);print;' temp3.pdb > temp4.pdb
 
gmx editconf -f temp4.pdb -o init.pdb -resnr 0
 
#rm temp1.pdb temp2.pdb temp3.pdb temp4.pdb

これをstep3ディレクトリ上でコピペして実行すればうまく目的のinit.pdbが生成されるはずです。temp1~temp4.pdbを見ながら処理を確認してみましょう。

このinit.pdbを用いていよいよONIOMインプットファイルを生成します。

5. ONIOMインプットファイルの作成

ここでは最初に述べたとおり、基質と反応に重要な残基を「コア領域」として選択、またそこから7 A以内の残基の構造最適化を許可、それより外の残基を動かないよう設定することを行います。

ここではtaopackageに含まれるコマンドpdb2oniomを用いることで、pdbファイルからONIOMインプットファイルを自動生成してくれるのですが、このコマンドを用いる上で2点前準備が必要となります。それは
1. コアファイルリストcorelist.txtの作成
2. /path/to/taopackage/ESPT/prepfilesにAMBER用のprepファイルのコピーを入れておく
ことです。

まずコア領域に設定する残基を宣言するファイルcorelist.txtファイルを作ります。このファイルの中身は

[TYR] "201"
 
[ARG] "214"
 
[PRO] "293"
 
[TYR] "385"
 
[FAP] "396"
 
[PHA] "397"

のようにテキスト形式で記述します。書き方のフォーマットは上のように、3文字の残基名とシーケンスIDをともに記述します。これはインプットとなるinit.pdbに含まれるものと同じようにします。
コア領域に設定する基準ですが、まずFAPとPHAは化学反応の中心点なので含めます。次にTYR201, PRO293、TYR385は、文献情報から水素原子の移動反応に大きく関わるとのことなので、含めます。ARG214ですが、これは反応の中心場となるPHAのCOO-と密接に塩橋を形成しています。このような電子雲を広げられそうな残基は、領域に含めた方がより良い計算精度につながります。Gly295も入れればよかったかな。できる限り広いほうが良いのですが、広すぎるとONIOMの計算速度の低下を招きますので注意。

続いて、AMBER用prepファイルのコピーを、インストールされたtaopackageディレクトリの下に存在するESPT/prepfilesに置くことです。taopackageでインストールされたコマンドは、自動的にここにpathが通っています。
上記3. MDで10 nsほど動かすで用いた残基パラメータファイルperoxide2.prep, anion_pohb.prepをここに入れます（デフォルトでchcl3.prep, zn.prepなどが置かれているディレクトリです）。
ここでこれらのprepファイルには1つ重要な追記をしなければなりません。それは、prepファイルの原子の設定部分の最終列に、原子番号を追記するということです。具体的には

    0    0    2
 
 
 
This is a remark line
 
PHA.res
 
PHA    INT  0
 
CORRECT     OMIT DU   BEG
 
   0.00000                                                      
 
   1  DUMM  DU    M    0  -1  -2     0.000      .0        .0      .00000
 
   2  DUMM  DU    M    1   0  -1     1.449      .0        .0      .00000
 
   3  DUMM  DU    M    2   1   0     1.523   111.21       .0      .00000
 
   4  O1'   o     M    3   2   1     1.540   111.208  -180.000 -0.862070
 
   5  C1'   c     M    4   3   2     1.269    63.213     0.000  1.017970
 
   6  O2'   o     E    5   4   3     1.269   126.426   180.000 -0.862060
 
   7  C1    cc    M    5   4   3     1.534   116.787    -0.000 -0.457170
 
   8  C2    cc    M    7   5   4     1.409   122.212     0.000  0.034310
 
   9  H1    ha    E    8   7   5     1.090   116.518     0.000  0.031060
 
（以下省略）

を

    0    0    2
 
 
 
This is a remark line
 
PHA.res
 
PHA    INT  0
 
CORRECT     OMIT DU   BEG
 
   0.00000                                                      
 
   1  DUMM  DU    M    0  -1  -2     0.000      .0        .0      .00000
 
   2  DUMM  DU    M    1   0  -1     1.449      .0        .0      .00000
 
   3  DUMM  DU    M    2   1   0     1.523   111.21       .0      .00000
 
   4  O1'   o     M    3   2   1     1.540   111.208  -180.000 -0.862070 8
 
   5  C1'   c     M    4   3   2     1.269    63.213     0.000  1.017970 6
 
   6  O2'   o     E    5   4   3     1.269   126.426   180.000 -0.862060 8
 
   7  C1    cc    M    5   4   3     1.534   116.787    -0.000 -0.457170 6
 
   8  C2    cc    M    7   5   4     1.409   122.212     0.000  0.034310 6
 
   9  H1    ha    E    8   7   5     1.090   116.518     0.000  0.031060 1
 
（以下省略）

というように、部分電荷情報の右にスペースをはさんで原子番号を添えてあげることです。この操作がないと、pdb2oniomで原子をうまく割り当ててくれることができないようです。

ESPT/prepfilesにtleapで用いた時のprepファイルを置き、上のように原子番号の追記をしてやると、準備完了です。
次のコマンドを打ってみます。詳細はpdb2oniom -hで出てくるヘルプを見てください。ここでは、インプットPDBファイルをinit.pdbとし、コア残基リストcorelist.txtの周辺から7 Aの残基を可動領域にし、その他の原子を固定します。これで、Gaussian 09インプットファイルoniominput.gjfが出力されます。また同時にoniominput.gjf.onbも出力されますが、これは大事なファイルなので間違って消さないようにしましょう。

このinit.pdbにはchain IDが含まれていてはいけないことに注意してください。

6. ONIOMインプットファイルの修正と実行

さて、これでoniominput.gjfで得られたと思いますので、それを早速見てみます。中身はこんな感じになっているはずです。

%chk=f_w34h310nsp1_1.gjf.chk　（チェックポイントファイル名の指定）
 
%mem=3700MB　（計算で使用するメモリ量の指定）
 
%nprocshared=4　（計算で使用する並列プロセッサ数の指定）
 
 
 
#P ONIOM(B3LYP/6-31G(d):AMBER=HardFirst) geom=connectivity nosymm iop(2/15=3) test　（Gaussianで用いる計算オプションの指定はすべてここ。以下で詳細を述べる。オプションを設定し終わるまで空白行を使わない。）
 
　　　　　　　（←1行空白行を入れる）
 
ONIOM inputfile generated by pdb2oniom from PDB file init.pdb. No connectivity generated.
 
Please use GaussView read f_w34h310nsp1_1.gjf, and generate connectivity information.（コメント欄。自由記述可。空白行は使わない。）
 
　　　　　　　（←1行空白行）
 
0 1 0 1 0 1　　（系の全電荷とスピン多重度の指定。以下で詳細を述べる。）
 
 H-HC-0.076010   -1       10.0670000      11.8690000      28.7070000   L　（以降原子データ）
 
 C-CT--0.190264  -1        9.3690000      12.6950000      28.5680000   L
 
 H-HC-0.076011   -1        9.8460000      13.6340000      28.8470000   L
 
 H-HC-0.076010   -1        9.3060000      12.7520000      27.4820000   L
 
 C-C-0.512403    -1        7.9330000      12.4750000      29.1960000   L
 
 O-O--0.550170   -1        7.4170000      13.1080000      30.1720000   L
 
（以下省略）

まずはこれについて述べます。
Gaussian 09のインプットファイルの詳細な記法はここで述べるよりは本家を見たほうが絶対いいのでそちらも参照。
重要なところだけ述べます。
%mem……例えばスパコンで計算させる場合、使える限りのメモリを使った方がお得。ただし、限界量よりは必ず少し小さめの値を入れておく。仕様を調べて値を入力しよう。例では3700MBとなっているが、現在の計算環境だと64GBメモリを積んでるマシンもあるので、それより小さい値、%mem=55GBとか60GBとかにしておく。

%nprocshared……これも最近のスパコンなどでは16コアとか24コアとか32コアとか普通にあるので、%nprocshared=16とか24とかにする。確かGaussian 09はノード間MPIはあまりできなかったような気がする（？）ので、MDみたいにノードをまたいでは指定しない。（間違ってたらごめんなさい）

#P ONIOM(B3LYP/6-31G(d):AMBER=HardFirst) geom=connectivity nosymm iop(2/15=3) test
ここのオプションが最大の肝。全部小文字で良いはず。まずpは「詳細な出力を表示させる」オプション。だいたいいつも入れておく。

oniom(b3lyp/6-31g(d):amber=hardfirst)……これはQM領域の計算レベルとB3LYP/6-31g(d)で、MM領域の計算レベルをAMBER力場で計算させ、全体エネルギーをONIOM式で外挿するという指定。amber=hardfirstのhardfirstは通常必ず指定する。QMの計算レベルだが、今回の系のように金属原子などが存在しない、かつ大員環が存在しないなど、pm6などの低い計算レベルが適用できそうな場合は先にこっちを使ってた方が収束が早いかも。この辺はQMの知識次第で応用可能。最初からb3lyp/6-31g(d)だと計算が遅いことがある。

geom=connectivity……MM領域が存在するONIOMのような場合には必ず指定する。原子間の結合を明確に規定する必要があるが、まだこの時点のインプットファイルではその記述が存在していない。これをきちんと書くためにはGaussViewなどのサポートがほぼ必須。後にさらに詳しく解説する。

nosymmとiop(2/15=3)……Gaussianでは計算量を減らすための工夫として、計算対象の分子を最初に解析し、分子対称性を検出してからできる限り楽をして計算を終了させようとする。分子対称性は量子計算では重要だからね、しょうがないね。さらにこのとき自動的に分子を回転させ、Gaussianが計算しやすいように座標を回転させる。しかし、taopackageにとってはこの分子の自動的な回転は好ましくないため、nosymmをつけることでそれを無効化する。またiop(2/15=3)はその保険であり、nosymmを忘れていてもこれが存在するとnosymm状態になる。
taopackageの解析ツールは固定された原子たちを検出して走るプログラムをいくつか用いているため、基本的には固定しておくとよい。

test……いらない。外してよい（はず）。

次に系全体の電荷の設定ですが、出力された時点では0 1 0 1 0 1となっています。
これは左から(real系の総電荷) (real系の総スピン多重度) (model系の高い計算レベルでの総電荷) (model系の高い計算レベルでの総スピン多重度) (model系の低い計算レベルでの総電荷) (model系の低い計算レベルでの総スピン多重度)を指定します。ここでreal系ってのはここでは蛋白質全体、model系はQMに設定した領域、計算レベルの高い低いはQM, MMレベルを表しています。
（公式ページ参照→http://www.gaussian.com/g_tech/g_ur/k_oniom.htm）
蛋白質全体の電荷はMDのtleapコマンドで系全体が0となったとき、Na+がCl-よりも7つ多かったことから、蛋白質の総電荷は-7、スピン多重度は（特に変な原子がないので）1です。

さて残すところは
・QM領域の設定、それに付随するlink atomの設定
・connectivityの設定
なのですが、これらの設定はGaussviewがあれば簡単に作れます。しかしないと簡単には作れません。

……でも長くなってしまったのでここでいったん切って次回以降に回します。
せっかちなひとは本家ページを読んで自習してください。