Gaussian 16 に実装されているONIOMを使ってタンパク質PHBH内の水酸化反応の遷移エネルギーを求める方法のチュートリアル

環境

UNIX環境（macOSまたはLinux）
AmberTools 22がインストールされている（一部AmberToolsのプログラムを使用する）
AMBER/GROMACSなどでMDを動かす方法を習得している
Gaussian 16が使える環境にある
GaussView 6がインストールされていることが強く推奨される

目的

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

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

TAO packageのインストール

こちらの使用条件は

UNIX環境（Linux, Macなど）
GaussView
Perl（普通のUNIX環境ならデフォルトで使えるはず）
VMD, PyMOLなどのタンパク質viewerソフト

です。Windowsの方はストアからUbuntu 22.04をインストールするなりして環境を整えてください（適当）。

インストールは単純ですが、一例を書いてみます。私はよくホームディレクトリのappsディレクトリにプログラムをインストールすることが多いので、そこに入れてみます。

# appsディレクトリを作成して移動
mkdir ~/apps ; cd ~/apps
# taopackageをGitHubからダウンロード
git clone https://github.com/YoshitakaMo/taopackage.git
cd taopackage
### taopackageのインストール場所を変更したい場合、install.shを編集して以下の部分を変更する
# 
# mkdir bin
# USERPATH=~/apps/taopackage/bin
# このUSERPATHの部分を書き換える
# 
###
# taopackageのインストール
./install.sh

終わったらご自身のお使いの.zshrcの最後の方に以下の行を加えてPATHを設定しておきます。

# taopackage
export PATH="${HOME}/apps/taopackage/bin:$PATH"

ONIOM法についての重要な参考文献

ONIOM法をやって論文を書くならまず目を通しておきたい論文。金属原子、特に遷移金属が含まれる場合にはDFTの計算レベルや基底関数の選択について他の量子化学計算上のテクニックや知識が必要となるが、ONIOMの流れ自体はこのチュートリアルと通じるものがあると分かってもらえると思う……。この他にも色々ONIOM法を使って研究した論文はあるので、検索して目を通してみよう。

TAO Packageの元論文

Matrix Metalloproteinase 2 Inhibition: Combined Quantum Mechanics and Molecular Mechanics Studies of the Inhibition Mechanism of (4-Phenoxyphenylsulfonyl)methylthiirane and Its Oxirane Analogue

諸熊先生のONIOMについてのReview記事

The ONIOM method: its foundation and applications to metalloenzymes and photobiology

ONIOM法を開発した故諸熊奎治先生のJACS論文

Gaussianの環境設定およびジョブの実行方法

MacまたはLinux OSを使っている場合の設定です。まずホームディレクトリ以下の~/.bash_profileに以下のような記述を加えます(bashを使用中の場合)。

~/.bash_profile

#
# Gauusian 16
#
export g16root=/Applications
. $g16root/g16/bsd/g16.profile
export GAUSS_SCRDIR=/Users/YoshitakaM/Scratch

ここはお使いのシェルに応じて編集するファイルを変えます（.tcshrc,.zshrcなど）。Linuxの場合はこれで問題ありませんが、Macユーザーの場合は.bashrcを編集するとなぜかHomebrewのアップデートが動作しなくなるという謎現象が発生するので、必ず~/.bash_profileにしておきます。

Macの場合は/Applications以下にGaussian 16本体の/Applications/g16とGaussView 6/Applications/gvを置いておきます。ディレクトリ名を変えると認識してくれない気がします。中のgviewはDockにでも置いておきましょう。Linuxの場合、上の設定が通っていればGaussian 16の実行コマンドであるg16やGaussviewを開いてくれるgvが動作するはずです。Macの人はシェルスクリプトを書いて動作させましょう。

Gaussian 16の計算を動作させるシェルスクリプトの例はこんな感じ。

run_g16.sh

#!/bin/sh

# Gaussian 16 Rev C01以降では環境変数GACSS_CDEFとGAUSS_MDEFで
# それぞれ使用するCPUコアとメモリの数を指定することができる。
# この場合インプットのgjfファイル内にそれらの値を記述する必要がなくなる。
# お使いのPCのスペックをもとに値を書き込む。
# GAUSS_CDEFは32コアならば"0-31"。
# GAUSS_MDEFは理論上の最大値を使用しないようにする。（例：RAMを32GB積んだマシンなら、設定値は25GBとかそのへん）
export GAUSS_CDEF="0-31"
export GAUSS_MDEF="25GB"

input=foobar

g16 < ${JOB}.gjf > ${JOB}.log ||  { echo "g16 optimization failed." ; exit 1 ; }

これでGaussian 16のインプットファイルfoobar.gjfをg16に計算させ、結果をfoobar.logファイルに書き出すというシェルスクリプト。Gaussian 16は基本的に標準出力に結果を書き出そうとします。

ONIOM法を使った遷移エネルギー決定のおおまかな流れ

だいたいこんな感じ

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

特に遷移状態での構造決定が難所です。これについては遷移状態の構造の予想がついていれば、それに近い点から遷移構造への構造最適化ができるのですが、予想がついていない場合は様々な可能性を試さなければならないのが計算化学の大変なところです。実験と違ってあらゆる可能性を考えなければ論文の査読で落とされかねないのですよね。

構造最適化によって構造が収束したあと、最低でも1つ確かめることがあります。それは標準振動解析の値であり、反応前・反応後の構造では構造最適化計算終了後の振動解析の固有値がすべて正の値になっていることが必要です（Gaussian 16の計算ファイルの場合、logファイル内に記述されている）。また、遷移状態の構造では、1つだけ負の固有値となっており（振動数が虚数解となり、これが遷移状態の山のてっぺんに相当する）、他の固有値は正となることが必要です。これを確かめることが最低限です。もちろん、その遷移状態が現実に即しているかどうかはまた別に検討する必要があります。詳しくは実践で見てみましょう。

実践ONIOM法（普通の構造最適化編）

Step 1. MDで10 nsほど動かす

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

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

中の解説をしますと、まず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ディレクトリでコマンドtleap -f leap.in を実行し、AmberTools 22のLEaPを動作させれば、AmberTools用のトポロジーファイルと初期座標ファイルが生成されます。後は普通のMDプロトコルに従ってAMBERやGROMACSを使って10nsほどNPTシミュレーションを動かしておきます。

Step 2. 5 Åより遠方の水分子を除外する

次に、MDが終わった後のタンパク質について、タンパク質＋リガンドをセンタリングした状態でPDBファイル形式に出力します。上記のダウンロード素材集ではこれを1iuw_10ns.pdbとして用意しました。まずはこれをPyMOLかなにかで見てください。タンパク質の周りに水が存在し溶媒ボックスいっぱいに詰まっていることがわかります。

計算で全部の水分子を入れてやるのは大変なので、タンパク質から5 Aより外の水分子はここで除外しておきます。
色んなやり方がありますが、PyMOLのコマンドを使って実行する方法を紹介します。

まずはPyMOL上で、タンパク質＋リガンドの残基をsequence一覧のところからマウスを使って範囲選択します。

次に、select byres (resn WAT within 5 of sele) or seleをコンソールに打ち込むと、現在選択されていた範囲(先のコマンドのsele部分に相当)に加えて、その周囲5Åの水分子(resn WAT within 5 of sele部分に相当)をまとめて選択することができます。byresをつけているのは原子単位ではなく残基単位で選択するためです。

画面を見て選択範囲が正しくタンパク質＋リガンド＋5Å内の水分子となっていることを確認したら、この選択範囲にある原子を別ファイル~/Desktop/init.pdb（デスクトップ上のinit.pdb）に出力します。コマンドではsave ~/Desktop/init.pdb, seleです。

この~/Desktop/init.pdbを改めてpymolで開いてみると、たしかに5Åより外の水分子が削除されていることがわかります。

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

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

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

AmberToolsのtleapコマンドを用いて上記init.pdbからinit.parm7, init.rst7, そしてref.pdbを作り出す。
コアファイルリストcorelist.txtを作成する。
init.parm7とinit.rst7、そしてpdb2oniom.pyスクリプトを用いてGaussian 16のONIOM計算実行用ファイルのoniominput1.gjfファイルを作り出す

の操作をやっていきます。

3-1. AmberToolsの`tleap`コマンドの実行

以下のleap_step3_1.inを使ってinit.pdbからinit.parm7, init.rst7, そしてref.pdbを作り出すことを行います。これらのファイルは後にpdb2oniom.pyスクリプトでGaussian 16のONIOM計算実行用ファイルを作成するために使います。この作業はAmberToolsでtleapを使ってMDシミュレーションのトポロジーファイル.parm7と初期座標ファイル.rst7を作った経験のある方ならば慣れているものでしょう。

leap_step3_1.in

#Use Amber ff14SB force field
source leaprc.protein.ff14SB
source leaprc.water.tip3p
source leaprc.gaff2
loadamberparams frcmod.ionsjc_tip3p
loadamberprep ../leap.wat/anion_pohb.prep
loadamberprep ../leap.wat/peroxide2.prep

#read a PDB file
mol = loadPDB init.pdb
center mol
# 以下のboxサイズはタンパク質の大きさに合わせて適切なボックスサイズを指定する
# わからなかったら150 150 150くらいでたぶん良い。
set mol box {150 150 150}

charge mol

saveAmberParm mol init.parm7 init.rst7

savePDB mol ref.pdb
quit

このファイルを作成します。ここで、loadamberparams, loadamberprepのところにはMDシミュレーションの操作のときと同様に追加のパラメータファイル（prepやmol2形式のファイル）へのファイルパスを設定します。つまり、../leap.wat/anion_pohb.prepや../leap.wat/peroxide2.prepの部分は適切なファイルパスに変更してください。そしてmolの部分にも先程作成したinit.pdbへのファイルパスを設定します。
今回は、init.pdbに存在するタンパク質＋リガンド分子・水分子についてのトポロジーファイルと初期座標ファイルが欲しいので、系の電荷を中和するためのNa⁺, Cl^-イオンを設置したり水分子を追加で撒くためsolvateBoxのコマンドは必要ありません。

これを作成したら、ターミナルから（AMBER/AmberToolsの実行コマンドへのパスが環境変数PATHに追加されている状態で）tleap -f leap_step3_1.inを実行します。その後、エラー無く0バイトでないinit.parm7, init.rst7とref.pdbファイルが作成されれば成功です。このあたりはMDシミュレーションのチュートリアルも参考にしてください。

3-2. AMBER/AmberToolsでの構造最適化計算（minimization）を行う

今作られたファイルをinit.parm7, init.rst7を用いて、AMBERを用いた簡単なエネルギー最小化の構造最適化計算を行います。このあたりは分子動力学シミュレーションの流れと全く同じです。有償ソフトウェアのAMBERを持っている場合はsanderまたはsander.MPIコマンドの代わりにpmemd（pmemd.MPI）で代用することができ、こちらのほうが計算が早いです（が、minimization計算ならばどちらを使っても1〜2分以内には終わるはずです）

この作業はチュートリアルファイルのstep3/minimizeディレクトリで行います。サンプルのrun.shファイルを動かして、構造最適化を実行します（要AmberTools22）。

$ cd step3/minimize
$ ./run.sh

この計算でエラーが無ければ、構造最適化されたあとの座標ファイルmin2.rst7ファイルが得られます。これとinit.parm7を3-4のpdb2oniom.pyのインプットファイルに使います。

3-3. コアファイルリスト`corelist.txt`を作成

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

corelist.txt

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

のようにテキスト形式で記述します。書き方のフォーマットは上のように、3文字の残基名と残基番号をともに記述します。ただし、この残基番号は仕様上先程生成されたref.pdbの番号に揃える必要があることに注意してください（生物学的な残基番号と必ずしも一致しない）。
コア領域に設定する基準ですが、まずこの反応系の場合ですと、FAPとPHAは化学反応の中心点なので含めます。次にTYR201, PRO293、TYR385は、文献情報から水素原子の移動反応に大きく関わるとのことなので、含めます。ARG214ですが、これは反応の中心場となるPHAのCOO^-と密接に塩橋を形成しています。このような電子雲を広げられそうな残基は、領域に含めた方がより良い計算精度につながります。 SER212, TYR222, Gly295も入れればよかったかな。できる限り広いほうが良いのですが、QM領域の原子数が増えると原子数の2~3乗に比例してONIOMの計算速度が低下していきますので注意。

このコアリストファイルを用いて次に進みます。

3-4. pdb2oniom.pyの実行

https://github.com/BILAB/pdb2oniom にて私が作成して公開しているpdb2oniom.pyを利用して、init.parm7とその座標に対応する構造最適化済のmin2.rst7（parm7に対応していればrst7ファイルは何でもよい。init.rst7も使用可能）, corelist.txtを入力としてGaussian 16でのONIOMインプットファイルを生成します。
まず適当な方法でここに置いてあるpdb2oniom.pyスクリプトをダウンロードし、どこかへ置いておきます。今回は~/apps/以下に置くことにします。ターミナルコマンドで以下のように打ちます。

mkdir -p ~/apps
curl -o ~/apps/pdb2oniom.py https://raw.githubusercontent.com/BILAB/pdb2oniom/main/pdb2oniom.py
chmod +x ~/apps/pdb2oniom.py

また、お使いのPython環境に合わせてライブラリをインストールしておきます。

python3 -m pip install absl-py numpy scipy ParmEd

そして、これを用いて（init.parm7とinit.rst7があるディレクトリ上で）以下のコマンドを実行します。

python3 ~/apps/pdb2oniom.py -p init.parm7 -r minimize/min2.rst7 --resid corelist.txt -o oniominput1.gjf

これにて、エラーがなければスムーズにONIOMインプットファイルoniominput1.gjfとoniominput1.gjf.onbが生成されます。

$ python3 ~/apps/pdb2oniom.py -p init.parm7 -r minimize/min2.rst7 --resid corelist.txt -o oniominput1.gjf

I1226 01:27:50.008807 8521628288 pdb2oniom.py:341] pdb2oniom.py ver.0.2; Amber parm and rst7 files to Gaussian ONIOM input file.
I1226 01:27:50.410809 8521628288 pdb2oniom.py:345] Total charge of low layer is calculated -7
I1226 01:27:50.410898 8521628288 pdb2oniom.py:346] Core residues list file corelist.txt provided.
I1226 01:27:50.410926 8521628288 pdb2oniom.py:347] All residues within 999.0 Å from core region are free to move (0) during geometry optimization.
I1226 01:27:50.461310 8521628288 pdb2oniom.py:137] 12983 atoms are set to be movable during geometry optimization.
I1226 01:27:50.923890 8521628288 pdb2oniom.py:359] Opening file oniominput1.gjf for output ...
I1226 01:27:50.924405 8521628288 pdb2oniom.py:362] Successfully wrote oniominput1.gjf file.
I1226 01:27:50.924434 8521628288 pdb2oniom.py:363] Opening file oniominput1.gjf.onb for output ...
I1226 01:27:50.924645 8521628288 pdb2oniom.py:367] Successfully wrote oniominput1.gjf.onb file.
I1226 01:27:50.924666 8521628288 pdb2oniom.py:368] pdb2oniom.py ends.

Step 4. ONIOMインプットファイルの意味

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

oniominput1.gjf

%chk=oniominput1.chk（チェックポイントファイル名の指定）
%mem=60GB （計算で使用するメモリ量の指定）
%nprocshared=16（計算で使用する並列プロセッサ数の指定）

#p opt=(quadmac,calcfc,maxstep=10) oniom(pm6:amber=softonly)
scf=qc geom=connectivity nosymm iop(2/15=3) Amber=(FirstEquiv,print)
（Gaussianで用いる計算オプションの指定はすべてここ。以下で詳細を述べる。オプションを設定し終わるまで空白行を使わない。）
　　　　　　　（←1行空白行を入れる）
ONIOM inputfile generated by pdb2oniom.py from init.parm7 and min2.rst7.
Please use GaussView 6 to modify the qm region and add the link atom info.（コメント欄。自由記述可。空白行は使わない。）
　　　　　　　（←1行空白行）
-7 1 0 1 0 1　　（系の全電荷とスピン多重度の指定。以下で詳細を述べる。）
N-N3-0.159200     0      7.222    11.448    28.600 L　（以降原子のデータ）
H-H-0.198400      0      7.638    10.543    28.434 L
H-H-0.198400      0      7.584    11.817    29.468 L
H-H-0.198400      0      7.285    11.995    27.753 L
C-CX-0.022100     0      5.774    11.252    28.832 L
H-HP-0.111600     0      5.519    11.608    29.830 L
C-C2C-0.086500    0      5.045    12.084    27.694 L
H-HC-0.012500     0      3.968    12.180    27.836 L
H-HC-0.012500     0      5.480    13.081    27.761 L
（以下省略）

まずはこれについて述べます。Gaussian 16のインプットファイルの詳細な記法はここで述べるよりは本家を見たほうが絶対いいのでそちらも参照。重要なところだけ述べます。

%mem……例えばスパコンで計算させる場合、使える限りのメモリを使った方がお得。ただし、限界量よりは必ず少し小さめの値を入れておく。仕様を調べて値を入力しよう。現在の計算環境だと64GBメモリを積んでるマシンもあるので、それより小さい値、60GBとかにしておく。

%nprocshared……これも最近のスパコンなどでは16コアとか24コアとか32コアとか普通にあるので、計算させるマシン環境に応じて%nprocshared=16とか24とかにする。確かGaussian 16はノード間MPIはあまりできなかったような気がする（？）ので、MDみたいにノードをまたいでは指定しない。Lindaが使えるならノード間MPIもできるがYayoiでは対応していない。

#p opt=(quadmac,calcfc,maxstep=10) oniom(pm6:amber=softonly) scf=qc geom=connectivity nosymm iop(2/15=3) Amber=(FirstEquiv,print) ここのオプションが最大の肝。

まずpは「詳細な出力を表示させる」オプション。だいたいいつも入れておく。

opt=(quadmac,calcfc,maxstep=10)……構造最適化を実行するoptオプションとそれに付随するもの。quadmacはQM/MM(ONIOM)をつけるときに推奨されるもの、calcfc、maxstep=10はhttps://www.hpc.co.jp/chem/software/gaussian/help/keywords/opt/ この辺を読んでね

oniom(pm6:amber=softonly)……これはQM領域の計算レベルがpm6で、MM領域の計算レベルをAMBER力場で計算させ、全体エネルギーをONIOM式で外挿するという指定。amber=softonlyは、gjfのファイルの下の方に分子力場のパラメータを明示的にすべて記述することで、それらを利用するようにする設定。ここをhardfirstにした場合は、Gaussian 16に入っている随分昔のAMBER力場のパラメータparm96を利用することになるが、古いのであまり推奨できない。QMの計算レベルだが、今回の系のように金属原子などが存在しない、かつ大員環が存在しない場合には、適切な計算レベルであるb3lyp/6-31g(d)をかける前にpm6などの低い計算レベルを適用しておき、動作チェックをしておいたほうがよい（→Step 8.を参照）。この辺はQMの知識次第で応用可能。

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

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

Amber=(FirstEquiv,print)……AMBER力場における追加の設定。FirstEquivは後ほど言及するが、重複する力のパラメータがインプットgjfファイル内に定義されている場合に、最初のものを優先して読み込み、あとのものを無視する設定。printはデバッグ用のオプションで、インプットgjfファイルを通じて定義された力場パラメータ情報をlogファイルの最初に出力する（性質上なくても問題ない）。

次に系全体の電荷の設定ですが、出力された時点では-7 1 0 1 0 1となっています。このうち、左端の数値（-7）だけはインプットのparm7ファイルから自動的に計算されます。残りの値は各自手動で適切に変更する必要があります。

これは左から(real系の総電荷) (real系の総スピン多重度) (model系の高い計算レベルでの総電荷)
(model系の高い計算レベルでの総スピン多重度) (model系の低い計算レベルでの総電荷)
(model系の低い計算レベルでの総スピン多重度)を指定します。ここでreal系ってのはここではタンパク質全体、model系はQMに設定した領域、計算レベルの高い低いはQM, MMレベルを表しています。（公式ページの"Per-layer Charge and Spin Multiplicity" タブ参照）

タンパク質全体の電荷はMDのtleapコマンドで系全体が0となったとき、Na+がCl-よりも7つ多かったことから、タンパク質の総電荷は-7、スピン多重度は（特に変な原子がないので）1です。

（2023年4月24日）インプットファイル中の ImpTrsで始まる行は削除し、oniom(XXXX:amber=softfirst)にしておいてください。

Step 5. QM領域に設定する原子の指定

上述のONIOMインプットファイルの意味をだいたい理解したら、まずはQM領域として計算させる原子をONIOMインプットファイル内で指定します。これはGaussViewを用いて設定したほうが早い上にミスも少ないので、こちらで設定することにします。

今回の例では、補酵素FADの過酸化状態FADHOOHのペルオキシドOH⁺が求電子剤となって基質p-OHBの芳香族に攻撃する、という化学反応なので、まずこれらに関わりそうな原子をQM領域に設定します。上のcorelist.txtでやっていた操作とは意味が異なることに注意しましょう。

GaussViewを用いている場合は前回の最後にできたファイルoniominput.gjfをそれで開きます。次に、Edit -> Atom Listを開くと、すべての原子の接続情報や原子種（Atom type）などのパラメータが一覧になっています。今回のファイルの場合、6246行目からFADHOOHのパラメータ、次いでp-OHB、水分子(HOH)と続いています（この行は20種類のアミノ酸以外は空欄で表示されます）。

このLayerのところに書かれているLow, Highを切り替えることで、High layer（QM領域）とLow layer（MM領域）を設定していきます。
ここでONIOMにおいて、QM領域に設定するときのTIPSを紹介します。以下のことを守って領域を設定しましょう。

反応に最も関わる部分をQM領域にまず設定する。ここではFADHOOHのOH⁺部分とp-OHB。
それに伴い、電子雲が広がっていきそうな塩橋部分、水素結合ネットワーク部分は含める。（SER212, ARG214, TYR201, TYR222, TYR385）
過去の論文を読んで、遷移状態に重要となることがわかっている残基は含める（PRO293、THR294）
芳香環は可能な限りすべて含める（FADHOOHとか例えばポルフィリン環全部とか）
各アミノ酸の側鎖については、側鎖の末端だけでなくCH₂基まで含めたほうがよい（Link atomの兼ね合い）。
プロリン残基をQM領域に含まなければならない場合、プロリン全体をQM領域に含め、さらにC末端側の隣のアミノ酸のNH-CHまでは含める（ペプチド結合の兼ね合い）

これがすべて設定し終えたら、Ctrl+G（Command+G）を押してインプットファイルの保存を行います。(input2.gjf)

こうしてできたinput2.gjfをテキストエディタで開いてずっと下っていくと、例えば3350行めあたりにはこんな感じになっています。

さて一行ずつ解説していきますと、

C-CX--0.263700……これは記号-で意味が句切られており、それぞれ左から原子C（炭素）、Amber力場形式のatom type（原子種）：CX（これはCα原子に割り当てられている事が多い）、最後に点電荷-0.263700ということを示しています。
他の例N-N2--0.529500で見ますと、原子は窒素、Amber atom typeはN2、点電荷は-0.529500ということです。これらはpdb2oniom.pyを使って自動生成されましたので、この辺は初心者はいじらないでください。

次の空白文字のあとに来る数字は0または-1となっています。0ならばその原子は構造最適化中で動くことが許可されている、-1ならば完全に固定されている（＝構造最適化中で動かない）ことになっています。これはpdb2oniom.pyのオプション-near 7 -resid corelist.txtとすると、corelistに書かれた残基から周囲7Åの残基が0に、それより外の原子が-1に設定されます。ただデフォルトでは999Åに設定されており、全原子が0に設定されて構造最適化対象になっています。このフラグの一括変更は後に紹介するsetmvflgコマンドでも行うことができ、後で変更する場合はそちらを使うことを推奨します。

その後の3つ組の数字はcartesianの$x,y,z$座標です。いずれの場合も、1字以上の任意の空白で区切ります。

さて、最後に書かれているLまたはHは、それぞれMM領域またはQM領域を表しています（Low/High layer）。先ほどHigh layerに設定した原子はここでHになっているはずです。ここで、エディタの3345行め（molecule specificationセクションにおける3335行め）にはL H-H1 3338となっています。これはLink atomを表していて、この原子自体はLow layer（MM領域）なのですが、ちょうどこの原子の隣にはHigh layerに設定した原子が存在しており、そちら側から見た時の計算処理をHigh layerのatom type:H1として計算し、molecule specificationセクションにおける3338行め（インプットファイル全体の3338行め、ではない）の原子に接続するということを表しています。ただ、このLink atomのatom typeの設定は後ほど手動で修正することになりますので注意。

上の図のように、Low layer（描画上はほっそい線）の3335番目の原子がHigh Layer（ボール状）の3338番目の原子に接続しているのがわかると思います。これらはそれぞれARG214のCγ、Cδ炭素原子です。3338の原子はMM領域からみたときにLink atom部分に設定されたatom typeに応じて力場の計算が行われることになります。注意点として、ONIOMにおいてはHigh layer（QM領域）に設定されている原子においてもlow layerレベルの力場パラメータが計算上設定されていることに留意してください。

今のチュートリアルではLow/High layerの原子の設定をGaussViewの上で行うやり方を紹介していましたが、実際にはテキストエディタでこのgjfファイルを開き、molecule specificationセクションの右端に書かれたLまたはHを書き換えれば同じことができますので、仕組みに慣れてくるとそちらでやった方が早いかもしれません。ここでLink atomの記述文法はGaussViewで修正したgjfファイルを開いて保存し直すと適切に書き換えてくだるので、その機能も利用すると便利です。

こうして無事QMとMM領域を設定できれば次に進みます。

Step 6. QM領域・MM領域に対応する電荷・スピン多重度情報を入れる

次にやるべきは、各領域(layer)の総電荷情報とスピン多重度を正確に入力することです。画面のように、座標が始まる前の行が-7 1 0 1 0 1となっていますが、ここは左から「Real系(Low-layer)の総電荷」「Real系(Low-layer)のスピン多重度」「Model系(High-layer)の総電荷」「Model系(High-layer)のスピン多重度」「Model系(Low-layer)の総電荷」「Model系(Low-layer)のスピン多重度」というふうになっています。ここでReal/Model系というのはタンパク質の全体領域/QM領域だけだと思ってください（2層式ONIOM計算の場合）。High/Low layerについてはここではQM計算レベル/MM計算レベルということです。スピン多重度というのは、不対電子の数から計算される数です。例えばFe³⁺で高スピン状態であれば、不対電子は5つあり、$2 \times (5/2) + 1 = 6$です。不対電子がない場合は$2\times (0/2)+1=1$です。

QM領域を構成する電荷とスピン多重度を丁寧に考えて記述します。今回のケースはQM領域の総電荷は-1となり、スピン多重度情報について、今回は不対電子がないので全部1です。つまり-7 1 -1 1 -1 1ということになります。もしアニオンや遷移金属原子、金属触媒がQM領域にある場合にはスピン多重度にはよく注意しましょう。

Step 7. Link atomのatom type定義の手動修正

GaussView用いてQM領域を再設定して保存し終わったら、それをテキストエディタで開いてLink atom部分のatom typeを手動で修正していきます。これはテキストエディタ上でL HをサーチすることでLink atom定義部分に飛ぶことができます。Link atomの設定は大変重要で、計算誤差を小さくするために、可能な限りsp3炭素–sp3炭素の間でQM/MMの境界を設定するようにし、Link atomは水素原子のatom typeに設定してください。例えばGLN残基ですと

 N-N--0.415700    0   55.516145     58.120294     39.140041    L
 H-H-0.271900     0   56.191534     57.533215     38.668463    L
 C-CX--0.003100   0   55.805461     58.495724     40.527475    L
 H-H1-0.085000    0   54.871223     58.542282     41.090271    L
 C-C2C--0.003600  0   56.704958     57.429029     41.175117    L H-HC 3176
 H-HC-0.017100    0   57.729416     57.600316     40.854173    L
 H-HC-0.017100    0   56.652163     57.546842     42.258385    L
 C-C2C--0.006082  0   56.334841     55.985201     40.791690    H
 H-HC-0.016965    0   55.254055     55.866739     40.925788    H
 H-HC--0.006049   0   56.557869     55.838379     39.730222    H
 C-C-0.684925     0   57.001127     54.910356     41.651831    H
 O-O--0.592044    0   56.871659     54.858651     42.867956    H
 N-N--0.907966    0   57.734731     53.981313     40.957641    H
 H-H-0.343885     0   57.703691     54.044107     39.934833    H
 H-H-0.382599     0   57.696986     53.025417     41.319102    H
 C-C-0.597300     0   56.464328     59.882616     40.575887    L
 O-O--0.567900    0   56.040210     60.753752     41.327674    L

こんな感じで、Link atomがGLN残基のCβ炭素とCγ炭素の間でQM(H)とMM(L)領域を分けています。そしてLink atomのatom typeをHCに設定しています。

これを全部のLink atom定義部分にて修正します。終わったものをinput3.gjfとして保存します。

Step 8. テスト実行

テストランを通じて、ONIOM計算を完全に開始することができるかどうかを検証します。ONIOM計算で最も設定がややこしいのは力場の定義部分ですが、これが正しく定義されているかどうかはGaussian 16を使って実際に計算してみるまでわかりにくい上にGaussian 16上でSCF Doneが来るまでわからないということもあるので、はじめに使うQM領域の計算レベルはpm6など半経験的分子軌道法で計算量が低いものにしておいた方がよく、その後DFTのb3lyp/6-31g(d)などに設定してあげる方が良いです。

Gaussian 16を使ってinput3.gjfについて計算を始めます。

$ g16 < input3.gjf > input3.log

これでgjfファイルに問題がなければ、計算が落ちずに続いているときにinput3.logのファイル内に、

SCF Done:  E(RPM6) = -0.394892611893     A.U. after   23 cycles
            NFock= 22  Conv=0.86D-08     -V/T= 0.9988
 KE=-3.174653381021D+02 PE=-1.190928014153D+04 EE= 6.001125295994D+03
 Leave Link  502 at Wed Apr 20 20:16:35 2016, MaxMem=  5368709120 cpu:        18.4
 (Enter /home/apps/g16/l601.exe)
 Copying SCF densities to generalized density rwf, IOpCl= 0 IROHF=0.

こんな感じでSCF Doneという行が現れているとひとまず成功です。この確認はターミナルから

# grepコマンドで文字検索を行う
$ grep "SCF Done" input3.log

という風にやってもOKです。これが見えましたら、この計算を止めてinput3.gjfのroute sectionを書き換え、pm6の計算レベルをb3lyp/6-31g(d)など適切なものにしてあげましょう。これをinput4.gjfとして以下の項に進むことにします。

余談ですが、pdb2oniom.pyを使うことでRoute sectionに出力されているAmber=(FirstEquiv,print)のFirstEquivを外して計算しようとすると、すぐにエラー落ちすることになります。

 Define LH2       23
 Define LNA       24
 Define LCD       25
 Define LH5       26
 Define LNC       27
 Define LNB       28
 Define LCC       29
 Include all MM classes
 Multiple matching entries:  320  322  321  323
 ImpTrs NB    NA    CR    H5     1.1 180.0 2.0
 ImpTrs NA    NB    CR    H5     1.1 180.0 2.0
 Multiple matching entries:  571  574  573  575
 ImpTrs C     CA    CA    HA     1.1 180.0 2.0
 ImpTrs CA    C     CA    HA     1.1 180.0 2.0
 Multiple matching entries:  575  579  578  580
 ImpTrs C     CA    CA    HA     1.1 180.0 2.0
 ImpTrs CA    C     CA    HA     1.1 180.0 2.0
 Multiple matching entries: 1112 1114 1113 1115
 ImpTrs NB    NA    CR    H5     1.1 180.0 2.0
 ImpTrs NA    NB    CR    H5     1.1 180.0 2.0
 Multiple matching entries: 1559 1562 1561 1563
 ImpTrs C     CA    CA    HA     1.1 180.0 2.0
 ImpTrs CA    C     CA    HA     1.1 180.0 2.0
 Multiple matching entries: 1563 1567 1566 1568

... (中略)

 Multiple matching entries: 6172 6176 6175 6177
 ImpTrs C     CA    CA    HA     1.1 180.0 2.0
 ImpTrs CA    C     CA    HA     1.1 180.0 2.0
 Inconsistent parameter file.
 Error termination via Lnk1e in g16/l101.exe at Sat Dec 24 11:42:10 2022.
 Job cpu time:       0 days  0 hours  1 minutes 44.7 seconds.
 Elapsed time:       0 days  0 hours  0 minutes  6.6 seconds.
 File lengths (MBytes):  RWF=      6 Int=      0 D2E=      0 Chk=      1 Scr=      1

これはinput3.gjfの最後の方に追記した力場パラメータのImpTrsで始まるところに、ある4つの原子についての二面角の異なるパラメータ定義が重複しているため、この問題を解消しなければならないということを表しています。大抵の場合、これらについて片方を削除すれば解決します。見た目には多くの重複定義が存在しているようになっていますが、実際に削除する必要があるImpTrs定義はわずかですので、手動で重複定義を削っていってください。（余談終わり）

Step 9. ONIOM計算の本番実行

以降は計算レベルをpm6から適切なものに変更したinput4.gjfを用いた本番の構造最適化計算を行います。このONIOM計算はとても時間がかかるのでどこかしらのスパコン上で計算をさせるのが一般的です。したがって各スパコンの上でqsub, pjsubコマンドでGaussian 16ジョブを投入することになります。ジョブの投入の仕方などは各スパコンの計算マニュアルを参照してください。

いずれにしろGaussian 16の計算は

$ g16 < input4.gjf > input4.log

のような形で計算を始めることになるはずです。

ONIOM構造最適化計算のモニタリング

上記で作成したinput4.gjfを用いたONIOM計算がスパコン上などで実行されている状態であると仮定します（計算レベルをpm6から適切なものに変更するのを忘れないでください）。すると、input4.logにはどんどん計算ログが書き込まれていくことになると思いますが、このファイルをずっと眺めるのは初心者にはとてもつらいものです。何が書かれているかまったくわからないことでしょう。

エネルギー確認

ONIOMの構造最適化ジョブを走らせているlogファイルに対して、TAO packageに含まれるoniomjobコマンドはその進行をチェックするのに使えます。以下の例はinput4.logの、構造最適化ジョブがまだ終了していないlogファイルです。次のようにコマンドを打つことで、logファイル中のエネルギーの概要を得ることができます。

$ oniomlog -o -i input4.log

このコマンドのアウトプットは次の通り:

ONIOMLOG 1.0 : Summary and monitor tool for Gaussian ONIOM log files.


Gaussian out file is input4.log

Optimization -o was used
Final ONIOM energy report for a two-layer ONIOM calculation (Hartree):
  This corresponds to the last complete step of an unfinished geometry optimization job.

High Level Model:           -2877.924883
Low  Level Model:               0.360940
Low  Level Real :             -64.827383
Low  Level Real-Model :       -65.188323

ONIOM Energy :              -2943.113207


Dipole moment (Debye) (X, Y, Z) is
(         6.4002,        -10.6791,         67.9132).
Total Dipole moment (Debye) is         69.0450.

Attention: this Gaussian calculation did not terminate normally!

OPT flag -o was turned on. Energy of each step will be printed.
Energies along the optimization path (kcal/mol):
Step Number          ONIOM              High Model          Low (Real-Model)
   1             15536.857674            63.573356         15473.284318
   2               293.603348            51.399887           242.203460
   3               244.103870            36.025721           208.078149
   4               143.893481            15.944185           127.949295
   5                88.800006            51.580556            37.219450
   6                66.819528            15.737673            51.081855
   7                52.156415             7.494022            44.662393
   8                84.734517            38.111923            46.622594
   9               120.220961            85.625889            34.595072
  10               117.215455            82.823643            34.391812
  11                22.975480            10.386616            12.588865
  12                13.769774             4.390699             9.379075
  13                21.004404             7.145781            13.858623
  14                13.643621             3.344198            10.299423
  15                14.034374             3.955874            10.078499
  16                12.190951             2.790695             9.400256
  17                10.530386             1.054049             9.476337
  18                18.214029            10.131328             8.082701
  19                10.069752             0.997369             9.072382
  20                 9.326937             0.840606             8.486331
  21                 8.666753             0.725519             7.941234
  22                10.566923             0.869868             9.697054
  23                 8.509402             0.524309             7.985093
  24                 7.593707             0.095390             7.498317
  25                 6.813307            -0.110568             6.923876
  26                 6.805409            -0.192464             6.997874
  27                 6.364374            -0.236690             6.601064
  28                 6.168418            -0.333938             6.502355
  29                 4.365153            -0.463076             4.828229
  30                 3.560172            -0.722323             4.282495
  31                 1.762191            -0.997131             2.759322
  32                 0.236510             0.039953             0.196557
  33                 0.000000             0.000000             0.000000
There are 33 steps of optimization in this job.

最後のステップのONIOMエネルギーは単位hartreeで示されます。フラグ-oが指定された時、構造最適化の経路はこのoniomlogの最後に示されています。
ここにはこれまで7ステップの最適化が行われていることを表しています（ステップ1は初期構造）。
ONIOMエネルギー、QM領域のエネルギー（highモデル）とMM計算、MMReal-MMModel（Low model）からの寄与が最後のステップをリファレンスとしてkcal/molで表示されている（そのため、一度もまだSCF Doneが行われていない場合にはErrorが返されます）。フラグ-haが付いている場合，全エネルギーがkcal/molの相対エネルギーではなく単位hartreeでの合計エネルギーで表示されます。これらのエネルギー値から、計算が最小値に向かってステップが進んでいることが簡単にわかります。

上の例のように、ONIOMの構造最適化ステップが進むにつれて前ステップとのエネルギー差が小さくなっていくことがとても望ましい状態です。これは構造最適化によって構造エネルギーの極小点に少しずつ収束していっている状態です。しかし、場合によってはこのoniomlogで表示されるエネルギー情報がこのようになっておらず、

差分エネルギー表示がバグってしまっている場合
差分エネルギーが0に近づいているが、ステップごとにプラスマイナス方向に振動している場合

があります。①はこんなoniomlog -o -i中でこんな感じにStep 2以降で途中からバグって表示されるケースです（Step 1でバグ表示が現れる場合は、まだ構造最適化ステップが1つも進んでいない状態であるため、無視して良いです）。

Dipole moment (Debye) (X, Y, Z) is
(        35.1489,         10.1764,         65.4345).
Total Dipole moment (Debye) is         74.9712.

Attention: this Gaussian calculation did not terminate normally!

OPT flag -o was turned on. Energy of each step will be printed.
Energies along the optimization path (kcal/mol):
Step Number          ONIOM              High Model          Low (Real-Model)
OPT flag -o was turned on. Energy of each step will be printed.
Energies along the optimization path (kcal/mol):
Step Number          ONIOM              High Model          Low (Real-Model)
   1             15536.857674            63.573356         15473.284318
   2               293.603348            51.399887           242.203460
   3               244.103870            36.025721           208.078149
   4               143.893481            15.944185           127.949295
   5                88.800006            51.580556            37.219450
   6                66.819528            15.737673            51.081855
   7          37035728.380000         32811.683381      33819343.202888
Argument "******************" isn't numeric in subtraction (-) at /users/moriwaki/apps/taopackage/ESPT/Oniom.pm line 401, <LOGFILE> line 96723.
Argument "******************" isn't numeric in subtraction (-) at /users/moriwaki/apps/taopackage/ESPT/Oniom.pm line 401, <LOGFILE> line 96723.
Argument " ********************" isn't numeric in subtraction (-) at /users/moriwaki/apps/taopackage/ESPT/Oniom.pm line 443, <LOGFILE> line 96724.

これはエネルギー計算が発散してしまったときによくあり、直前までに出てきた構造をよく見ると系の中のある2つの原子が近づきすぎていたり、すでに原子座標が爆発して表示不可になっている場合があります。初期構造(原子の座標)やLink atomの設定、構造最適化する原子グループ設定があまり良くないことが原因として考えられます。まずはLink atom, 原子グループの設定を変えてみて、次に改善されなければ初期構造（1iuw_10ns.pdbやinit.rst7に含まれる座標）を新しく作り直して入れ直すと問題が解消されることがあります。②のケースでは、opt=(maxstep=10)の部分をmaxstep=5とかmaxstep=2とかにすることで改善することがあります。ただ、初期構造を作り直すことでも解決するかもしれません。

新規インプットファイルの生成

以前のONIOMジョブで得られたlogファイルに保存されている最適化構造（途中の構造を含む）の座標をもとに、新たにGaussianのONIOMジョブを作りたい時には、これもoniomlogコマンドで同様に行うことができます。

$ oniomlog -oi -t input4.gjf -i input4.log -fo newinput.gjf -fn 8

フラグ-oiによってプログラムoniomlogに新たなGaussian ONIOMインプットファイルを新規に生成させることができます。このプログラムではinput4.gjfをテンプレートとして(フラグ-tで指定する)、input4.log(フラグ-i)の構造最適化中に現れた8番目の構造（フラグ-fn 8）を抽出します。そして新たなGaussian ONIOMインプットファイルnewinput.gjf (フラグ-fo)を生成します。フラグ-fnの後の数字は指定不要で、デフォルトは0, 0は最後の構造を表します。-fn <int>で指定する数値は、oniomlog -o -i input4.log中で表示される構造最適化ステップの番号を見ながら選びます。

適切な理論レベル(level of theory)を見つけることと適切な状態に収束された構造最適化を得るためには、たくさんの計算回数を必要とすることがほとんどなので、これを活用しましょう。

また、時々ONIOM構造最適化のジョブが終わるまでにやけに時間がかかることがあります。ある場合においては、構造最適化が収束しないことがあります。例として、（先の説明にも出ましたが）たくさんの最適化ステップの末に2つのよく似た構造の間で最適化が循環してしまい、設定された最適化ステップ数を超えてしまう前に総エネルギーがとても小さな範囲の間で往復してしまう場合です（振動状態）。oniomlogはこの場合に特に有効で、GaussianのONIOMジョブがエラーメッセージにより終了してしまうずっと前に、構造最適化計算が振動状態になっているかどうかを調べることができます。これが起きた場合、oniomlogは最適化に沿った構造をチェックするのに使えるうえ、必要な変更を加えた新たなONIOMインプットファイルを生成して別の最適化ジョブに備えるのにも使えます。（例えば異なる理論レベルを使うことや，構造最適化における最大ステップ数の設定変更など）。

計算機のコストと時間を節約するため、ぜひoniomlogを使った積極的なモニタリングを推奨します。

最適化フラグの再設定

例えば最適化中に可動するコア領域を、中心から999Åではなく7Åのすべての残基を動かすというように変更したいという場合など、ONIOMジョブにおける構造最適化フラグを再設定したい場合には、setmvflgを使うと簡単です。サンプルは以下の通り。

$ setmvflg -b -i input4.gjf -onb oniominput1.gjf.onb -resid corelist.txt -near 7 -o input4_newflg.gjf

setmvflgのためにはONIOMインプットファイルであるinput4.gjf, pdb2oniomプログラムによってPDBファイルからONIOMインプットファイルを生成するときに作られたONBファイルであるoniominput1.gjf.onbおよびコア残基リストファイルcorelist.txtを必要とします。このコマンドでinput4.gjfにおける最適化フラグを最設定します。
input4_newflg.gjfでは、構造最適化中に動ける領域がコア領域から7Åの残基すべてとなるように再設定されたものになっています。

$ setmvflg -b -i input4.gjf -onb oniominput1.gjf.onb -resid corelist.txt -near 7 -o input4_newflg.gjf

SETMVFLG 1.0 : Resets optimization flag for a Gaussian ONIOM input file (0/-1)

Use coordinates from file oniominput1.gjf.onb to reset optimization flag.

Reading ONIOM input file input4.gjf...
There are 12983 atoms in file input4.gjf.
Done!

Reading ONIOM ONB file oniominput1.gjf.onb...
There are 12983 atoms in file oniominput1.gjf.onb.
There are 2608 residues in file oniominput1.gjf.onb.
Done!

Output file name is input4_newflg.gjf

In file input4.gjf, there are 12983 atoms, 2608 residues subject to geometry optimization.

After reset, there are 2846 atoms, 267 residues subject to geometry optimization.

Open input4_newflg.gjf to write new ONIOM input file...

Successfully wrote input4_newflg.gjf file.

SETMVFLG ends.

このアウトプットでは、input4.gjfでは2608残基（12983原子）が動いていたのに対し、再設定後のinput4_newflg.gjf中では267残基（2846原子）が動くよう再設定されたことがログメッセージで示されています。ここで、setmvflgコマンドは-resid corelist.txtから見て-near 7Å以内のフラグを0に、それより外を-1に設定しようとしますが、フラグ-bを用いた場合、onbファイルに書かれた最初期の座標を基準として使用することになり、このフラグをつけない場合、入力とするGaussianインプットファイル中の座標の方を基準として最適化フラグの再設定を行います。どちらの場合においても、新しく生成されたインプットファイルは与えられたインプットファイルと同じ座標を持つことになります。

部分構造確認

（この機能は私はあまり使わないのですが）このプログラムoniomlogは、構造最適化経路に沿って得られる構造の一部を迅速に抽出することもできます。

$ oniomlog -o -i input4.log -s input4_partqm.xyz

フラグ-oによって構造最適化経路に沿って変化するQM領域だけを抜き出し、それをinput4_partqm.xyzという名前でファイル形式xyz形式で保存します。
もしフラグ-oを指定しない場合は、最後の構造のみが保存されることになります。

$ oniomlog -g -l 1 -o -i input4.log -s input4_partqmmov.xyz

フラグ-g -l 1 -oをつけることで，構造最適化経路に沿ってQM領域とすべての可動系（gjfファイル内で座標xyzの前に0を指定した原子群）の部分を抽出することができ、input4_partqmmov.xyzに保存することができます。

input4_partqm.xyzとinput4_partqmmov.xyzは比較的小さなファイルで、VMDなどの可視化ソフトで手軽に見ることができます。
フラグ-l, -gとその組み合わせは様々な方法で，ONIOM計算から構造を抽出するのに使うことができる。詳しくはoniomlog -hで表示されるヘルプを参照。

ここまでの構造最適化計算がうまく進まないとき

以下の点を試します。

QM領域に入れている原子は多すぎやしませんか？入れすぎるとどんどん計算は重くなりますが、かと言って適切にQM領域を入れていないと計算結果が大きく変わることもあります。
計算レベルを高いものにしていませんか？b3lyp/6-31G(d)レベルは割とサクサク進みますが、b3lyp/6-31+g(d,p)になると精度が良くなる代わりに、QM領域に入れている原子数が多いとほとんど進みにくくなります。計算機の性能が高ければ進むこともあります。b3lyp/6-31G(d)である程度計算を進めておき、途中から計算レベルを上げるという技もあります。
異常終了したときに、logファイルの最後の方に現れるエラーメッセージが"FormBX had a problem."だった場合は、単純にそこから計算を再開させれば動くことが多いです。Gaussianの計算の再開を読んでください。
家庭用の計算機上で何もエラーメッセージが無いのに計算が異常終了している場合は、計算機の性能が足りていない場合があります。もっとよい計算機上で計算させてみると問題なく動くことがあります。

実践ONIOM法（遷移状態（Transition State）の探索）

上記の方法を行えば、とりあえずReactant State（反応物の状態）と、それからProduct State（反応後の生成物の状態）についてそれぞれ構造最適化をONIOM法で行うことができるはずです。なぜならば、構造最適化は開始構造からみてエネルギーの値とエネルギー微分（勾配ベクトル）を求め、エネルギーの最小化（実際には極小点へと限りなく近づく操作）を行っていっているからです。これを指定している命令がgjfファイル内の最初の方に書くopt=(~~)のあたりです。

さて、これに比べると遷移状態の構造探索は結構難しいものとなります。教科書などでの説明は、（そもそも紙が平面な都合もあって、）遷移状態は反応物と生成物の間に存在する一つの山のように描かれていることが多いのですが、実際にはこれは数学的には鞍点(saddle point)であり、ある1つのベクトル（ここでは反応物と生成物の間に存在する反応座標）に対してのみ上に凸な状態で、他の無数のベクトルに対しては下に凸なエネルギー地形となっています。
ここで概念的ですが反応物の状態を化学式で[A-B C]とし、遷移状態を[A…B…C]、生成物を[A B-C]というような書き方で表すとします。-は共有結合みたいなものだと考え、空白は十分距離が離れていると考えます。ここで量子化学計算の手法で遷移状態の構造を求めたい場合、以下のようにして行われることが多いです。

系の中で、大きく移動するはずの（または、移動しても良い）原子グループ（上の概念図でBに相当する）を頭のなかで決める
A点からB点までの距離を固定し、Bの位置をC点に近づけるような方向で大きくしながら何点かとり（または逆にB点とC点の間を固定し、この値を徐々に小さくしていく形で）、それぞれで構造最適化を行う。
前段階の操作でエネルギーの山となる点が見えてきたら、その近くの構造をインプットとして、opt=(TS)などをつけた遷移状態探索を開始する。
最後にIRC(lntrinsic Reaction Coordinate; 固有反応座標)の計算をつけて、反応物・生成物の状態に近づくことができれば終了。

という流れです。最近だとさらに便利なTS探索法があるようなのですが、ONIOMの上でそれを行う方法がまだ私には理解できていなかったので紹介できませんでした。

それではPHBHを例にとって実際にやってみましょう。とはいえ、今回のPHBHの系の場合、すでに過去の論文などで反応座標や遷移状態がすでに判明しています。まずは練習ということでこれを再現するためのチュートリアルをここでは紹介することになりますが、実際にはあらゆる可能性の遷移状態を考え、それぞれを比較検討することでもっともらしい反応座標を探索することになることを覚えておきましょう。

反応座標のscan

ここでは上述の流れのうち、1番目と2番めに相当する操作を行います。今回、Bに相当する「移動するはずの」原子グループはペルオキシドのOH+となります。そしてAとCはそれぞれFADと基質p-OHBに相当します。これを画像で示したものが次の図です。

GaussViewの機能を使うと簡単に原子団の座標をいじることができるので、それを使ってO6264とO6265の距離を図のように1.4600前後に設定してみます。ちなみにReactant
stateの構造最適化では、この2点間の距離はだいたい1.46前後となっています。

さてこの設定した距離を固定した状態で構造最適化を行うには、インプットファイルに次のように設定します。

Tools > Redundant Coordinatesを選んで……

図のように、bondについてのRedundant Coordinatesの設定を行います。ここではscan coordinatesを選び、1.46 Åから大きくなる方向に0.5Åずつ距離を変更しながら10回繰り返す操作、というふうに設定しています。ここは各自の好きなようにやりたいだけ設定してください。

このscan coordinatesという設定は、実際にGaussianを計算させたときには「まず初期値で与えられた2原子間の距離を固定して構造最適化する」→「続いて、その最適化された構造を初期座標として、その2点間の距離を設定された値だけ変更し、改めて構造最適化する」→「この繰り返しを設定された値分だけ回す」という流れになります。

これにてOKボタンをおし、この設定でインプットファイルを保存してgjfファイルの中身を見ますと、このようになっているはずです。

testscan1.gjf

%chk=testscan1.chk
%mem=80GB
%nprocshared=32
#p opt=(modredundant,maxstep=10,quadmacro)
oniom(b3lyp/6-31+g(d,p):amber=hardfirst) nosymm scf=qc geom=connectivity
iop(2/15=3)

Here is comment section.

-7 1 -1 1 -1 1
H-HC-0.076010   -1   47.64772000   55.50134400   77.54871900 L
C-CT--0.190264  -1   47.34974700   54.75746000   78.28462700 L
H-HC-0.076011   -1   47.90441500   53.83286400   78.13675700 L
H-HC-0.076010   -1   47.55837800   55.15213100   79.27407200 L
C-C-0.512403    -1   45.86841500   54.49005100   78.15935600 L
……
（座標情報セクション）
O-OW--0.834000  -1   13.77095200   48.53563200   41.44696900 L
H-HW-0.417000   -1   13.97077900   48.89671900   42.31293000 L
H-HW-0.417000   -1   13.76948500   47.57428000   41.59171700 L

1 2 1.0
2 3 1.0 4 1.0 5 1.0
3
4
……
（結合情報セクション）
12990 12991 1.0 12992 1.0
12991
12992

B 6264 6265 S 10 0.500000

VDW NH 1.8240  0.1700
VDW P5 2.1000  0.2000
VDW CD 1.9080  0.0860
VDW ND 1.8240  0.1700

上記のように、opt=(modredundant)とB 6264 6265 S 10 0.500000という情報が追加されており、これが実際にscan coordinateを動かす上で効いてくるインプット情報です。optなどと一緒に併記されているnosymmやiop=(2/15=3)も外さずにきちんとつけてあることを確認しましょう。なければ、ここで付け足してください。

そして、これを確認したら、実際にGaussian 16を使ってエネルギー構造最適化をしかけます。ただし、これまでのエネルギー構造最適化とは異なり、1つの点で構造最適化が終わっても続けて10ステップ分最適化を続けるので、計算終了までの時間はその分長くなります。また、途中でエラー落ちする確率も大幅に上がります。1点について構造最適化計算が収束する（下の例でConverged?がすべてYESになる）と、logファイルに以下のようなメッセージが追加されます。

         Item               Value     Threshold  Converged?
 Maximum Force            0.000250     0.000450     YES
 RMS     Force            0.000020     0.000300     YES
 Maximum Displacement     0.000862     0.001800     YES
 RMS     Displacement     0.000142     0.001200     YES
 Maximum MM Force         0.000010     0.000450     YES
 RMS     MM Force         0.000000     0.000300     YES
 Predicted change in Energy=-5.575057D-08
 Optimization completed.
    -- Stationary point found.

このメッセージの後、Gaussianプログラムは次の1点の構造に向けて自動的に距離を再設定し構造最適化を再開します。つまり、上の例だと11回分-- Stationary point found.のメッセージが現れることになります。grepコマンドでログファイルの中から"Stationary"という文字列を検索してみればよいと思います。
エラー落ちしてしまった場合は、あまりに構造がおかしくなっていないか確かめながら、改めて計算を仕掛けなおしてください。覚えておいてほしいのは、設定した2原子の距離の固定が適切に働きながら、生成物状態から反応物状態の方向へと向かっているというのさえ見えていればだいたいOKということです。この操作により、適切に計算が行われていれば、横軸にその固定した距離、縦軸にエネルギーをとってプロットしてやれば、どこかで極大点が見えるはずです。これが遷移状態の候補となります。上記の操作を通じて、まずはこれを発見しましょう。

TS（遷移状態）の構造探索

先の項の操作を行うことで、答えを言ってしまえばだいたい1.86 Åの距離が離れたあたりでエネルギーの極大点が現れているはずです。この状態の全体の構造は図のようになっています。

この構造を起点にして、opt=TSを使った遷移状態の構造最適化に入ります。ちなみにTSの計算はopt=calcfcを必要とします。
ここでさらに追加オプションでおすすめなのがopt=noeigentestです。こちらは公式によれば、計算量が非常に大きいときに推奨されるとしています。今回のような大きい系ではつけておいたほうが最初の方は計算が進みやすくなるので、収束の最後の方でopt=noeigentestを解除して再計算すると良いかもしれません。

EigenTest requests and NoEigenTest suppresses testing the curvature in Berny optimizations. The test is on by default only for transition states in internal (Z-matrix) or Cartesian coordinates, for which it is recommended. Occasionally, transition state optimizations converge even if the test is not passed, but NoEigenTest is only recommended for those with large computing budgets.

したがって、TS計算のためのインプットファイルのルートセクションは以下のような感じになります。

#p opt=(calcfc,ts,noeigentest,maxstep=10,quadmacro)
oniom(b3lyp/6-31+g(d,p):amber=hardfirst) nosymm scf=qc geom=connectivity
iop(2/15=3)

これをつけて、後はTSの構造最適化がうまくConvergedすることを祈りながら待ちます。

TSの構造最適化後の確認

TSの構造最適化計算が収束するとこんな感じになりました

最初1.86Åくらいからはじめましたが、1.92Å前後まで伸びていました。しかしながら、これが真なる遷移状態構造だと決めつけるのはまだ早いです。少なくとも次の2点を確認する必要があります。

得られたTS構造に対し、振動数計算freqをつけた計算を行って虚振動が1つだけとなっていることを確認する
得られたTS構造に対し、IRC計算を行い、望ましい原系(Reactants)、生成系(Products)へと座標が変化していっているかどうかを確認する。

IRC計算はルートセクションでopt=(~~)を外し、代わりにIRC=(calcfc,forward,maxpoints=10,StepSize=10)のような感じで設定します。ここで、forwardをreverseにすることでTS状態からの順方向、逆方向の構造遷移を観察することができます（どちらかがReactants, もう一方がProductsに対応します）。

3つの状態のエネルギー値の比較

これまでの計算がすべてうまく行っていれば、Reactant State, Transition State, Product Stateの3つを求めることができたと思います。TSからのIRC計算がうまく行っていれば、これらが反応座標上でうまくつながっており、もっともらしい反応経路が描けたと思います。

さて、こうして3つの状態のログファイルが求まったとき、次のようにして簡単に3つのエネルギー状態を比較することができます。まずgtable.inというファイルを作成し（ファイル名はなんでも良い）、その中にこれまでの3つのログファイルへのパスを記します。

<Label>HOGEHOGE</Label>

<ONIOM>
reactant/reactant1_1.log
ts/ts1_1.log
product/product1_1.log
</ONIOM>

そしてTAO Packageに含まれているコマンドでgaussiantable -i gtable.in -o gtable.outと入力すると、即座にエネルギーをまとめた表を次のように返してくれます。

=================================================
Output for block 1: HOGEHOGE
Absolute values in Hartee:
            ONIOM Total               Model High                   Model Low                    Real Low              Real-Model Low

reactant1_1.log
           -2600.483170             -2532.265295                    8.970611                  -59.247264                  -68.217875

ts1_1.log
           -2600.472865             -2532.263599                    0.271774                  -67.937492                  -68.209267

product1_1.log
           -2600.546542             -2532.357580                   18.174134                  -50.014827                  -68.188962

-------------------------------------------------
Relative values in kcal/mol:

            ONIOM Total               Model High                   Model Low                    Real Low              Real-Model Low

reactant1_1.log
               0.000000                 0.000000                    0.000000                    0.000000                    0.000000

ts1_1.log
               6.466701                 1.064807                -5458.602939                -5453.201045                    5.401893

product1_1.log
             -39.766064               -57.909523                 5775.298089                 5793.441547                   18.143458

=================================================

上の例の場合、Reactantに比べてTSエネルギーが相対的に$6.46 \textrm{kcal/mol}$ほど高く、Productが相対的に$39.76 \textrm{kcal/mol}$ほど低いということがわかりました。