量子化学計算プログラムGaussianを用いた構造最適化と配座探索の実行
CONFLEXで行う構造最適化および配座探索は、通常CONFLEX本体に搭載されている古典的な分子力場を用いて実行しています。エネルギー極小構造を高速かつ効率的に網羅する上で分子力場は有用ですが、対象とする系に必要な原子タイプやパラメーターが無い場合は、これらを補わなければ極小構造を求めることができません。また、力場では記述するのが難しい電子状態での極小構造を求めたい場合、力場を使って算出した配座について量子化学計算を行った場合と、最初から量子化学計算を行って配座探索を実行した場合とで、得られる配座が異なることも考えられます。
CONFLEX8以降では、分子力場を用いる代わりに量子化学計算プログラムGaussianを外部プログラムとして呼び出し、最適化計算や配座探索を実行することができます。ここでは、Gaussianを呼び出して配座探索計算を行う場合の設定方法と、その結果の例を示します。
【環境設定】
外部プログラムとしてGaussian（16または09）を呼び出すためには、CONFLEXがインストールされているマシンにGaussianもインストールされている必要があります。またGaussianを実行する際には、お使いのアカウントでGaussian実行用の環境設定があらかじめ為されている必要があります。
Mac, Linuxの場合
環境変数”g16root”または”g09root”をあらかじめ設定してください。尚、環境変数「GAUSS_SCRDIR」を指定していない場合は、中間ファイルがカレント・ディレクトリに出力されます。
Windowsの場合
環境変数「GAUSS_EXEDIR」に、Gaussianがインストールされているフォルダー名（C:¥G16W、等）を設定してください。尚、キーワード「EXT_GAU_SCRDIR=」に中間ファイルを出力するフォルダーを指定しない場合は、「%GAUSS_EXEDIR%¥Scratch」が出力ファイルおよび中間ファイルを出力するフォルダーとして設定されます。
また、環境設定を含めた設定ファイルを別途用意することで、Gaussianの環境設定を計算実行時のみ有効にする方法もあります。
【Propylene glycolのB3LYP/6-31G(d)レベルによる配座探索計算】
ここではPropylene glycolを利用して、B3LYP/6-31G(d)レベルで配座探索を行う方法を説明します。
Propylene glycolの入力構造
Propylene glycolの座標データ（PropyleneGlycol.mol）
PropyleneGlycol.mol 13 12 0 0 0 0 0 0 0 0 0 0 -4.0461 -0.3459 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 -3.6894 0.1585 0.8737 H 0 0 0 0 0 0 0 0 0 0 0 0 -3.6894 0.1585 -0.8737 H 0 0 0 0 0 0 0 0 0 0 0 0 -5.1161 -0.3459 0.0000 H 0 0 0 0 0 0 0 0 0 0 0 0 -3.5328 -1.7978 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 -3.8895 -2.3022 0.8737 H 0 0 0 0 0 0 0 0 0 0 0 0 -1.9928 -1.7979 -0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0 -1.6361 -1.2935 -0.8737 H 0 0 0 0 0 0 0 0 0 0 0 0 -1.6361 -1.2934 0.8736 H 0 0 0 0 0 0 0 0 0 0 0 0 -4.0095 -2.4719 -1.1676 O 0 0 0 0 0 0 0 0 0 0 0 0 -4.5428 -3.2264 -0.9068 H 0 0 0 0 0 0 0 0 0 0 0 0 -1.5162 -3.1461 0.0000 O 0 0 0 0 0 0 0 0 0 0 0 0 -0.5757 -3.1535 0.1927 H 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 1 3 1 0 0 0 0 1 4 1 0 0 0 0 1 5 1 0 0 0 0 5 6 1 0 0 0 0 5 7 1 0 0 0 0 5 10 1 0 0 0 0 7 8 1 0 0 0 0 7 9 1 0 0 0 0 7 12 1 0 0 0 0 10 11 1 0 0 0 0 12 13 1 0 0 0 0 M END
[Interfaceから実行する場合]
PropyleneGlycol.molファイルをCONFLEX Interfaceを用いて開きます。
Calculationメニューから「CONFLEX」を選択し、開いた計算設定ダイアログのをクリックします。詳細設定ダイアログが開きます。
まず、詳細設定ダイアログの「General Settings」ダイアログにある「Calculation Type:」のプルダウンメニューから「Conformation Search」を選択します。さらに、「Estimating Method:」のプルダウンメニューから「External Program」を選択します。これにより、構造最適化などの計算が外部プログラム（Gaussian）を利用して行われるようになります。
次に、「External Program」ダイアログで、Gaussianの計算設定を行います。
「Job Title:」はコメントであり、ここでは「Propylene Glycol B3LYP/6-31G(d) optimization」とします。「%Mem=」は「128」とし、単位のプルダウンメニューから「MW」を選択します。「%NProcShared=」は「2」とします。「Route Section:」は「B3LYP/6-31G(d) Opt」とし、B3LYP/6-31G(d)での構造最適化を行うように設定します。
Gaussianの計算設定は、ファイルを用いて行うこともできます。【ファイルを用いたGaussianの計算設定】をご参照ください。
次に、「Conformation Search」ダイアログで、「Search Limit:」の値を「5.0」にし、配座探索空間の範囲を決めます。
設定が終わりましたら、詳細設定ダイアログのをクリックします。
開いたダイアログに「EDIF_HARD=1.0D-3」を追加します。これは、新たに得られた配座を保存するかどうかを判定する際のエネルギー差の閾値を0.001 kcal/molとする指定です。
設定が終わりましたらをクリックします。B3LYP/6-31G(d)レベルで配座探索が行われます。
[コマンドラインから実行する場合]
計算設定は、PropyleneGlycol.iniファイルにキーワードを記述することで行います。
PropyleneGlycol.iniファイル
CONFLEX
SEL=5
EXTERNAL_PROGRAM=GAUSSIAN
EXT_GAU_LINK0=(mem=128MW)
EXT_GAU_LINK0=(nprocshared=2)
EXT_GAU_ROUTE=("B3LYP/6-31G(d) opt”)
EXT_GAU_COMMENT=("Propylene Glycol B3LYP/6-31G(d) optimization”)
EXT_GAU_CHARGE=0
EXT_GAU_SPIN=1
EXT_PREOPT=OFF
EDIF_HARD=1.0D-3
「CONFLEX」キーワードは、配座探索計算を行うこと意味します。
「SEL=5」キーワードは、配座探索のエネルギー上限値を5.0 kcal/molと設定します。また、「EXTERNAL_PROGRAM=」は使用する外部プログラム名、「EXT_GAU_LINK0=」はLink0コマンド（複数指定可）、「EXT_GAU_ROUTE=」はルートセクション、「EXT_GAU_COMMENT=」はコメント、「EXT_GAU_CHARGE=」は系全体の電荷、「EXT_GAU_SPIN=」は系のスピン多重度をそれぞれ表しています。尚、「EXT_GAU_ROUTE=」と「EXT_GAU_COMMENT=」では他と異なり「(“」と「”)」で指定しますので、ご注意ください。「EXT_PREOPT=OFF」はGaussianを呼び出す前に分子力場による最適化を行わないようにする指定、「EDIF_HARD=1.0D-3」は新たに得られた配座を保存するかどうかを判定する際のエネルギー差の閾値を0.001 kcal/molとする指定です。
Gaussianの計算設定は、ファイルを用いて行うこともできます。【ファイルを用いたGaussianの計算設定】をご参照ください。
PropyleneGlycol.molとPropyleneGlycol.iniをフォルダに格納し、下記コマンドを実行してください。計算が始まります。
C:\CONFLEX\bin\flex9a_win_x64.exe -par C:\CONFLEX\par PropyleneGlycolenter
上記は、Windowsの場合です。他のOSにおける実行コマンドについては、本書の「実行方法」を参照してください。
計算結果
PropyleneGlycol.ls1の出力内容とGaussianの入出力ファイルを示します。Gaussianの入出力ファイル名はls1ファイルの「CONF. ID」に対応しており、最安定構造はPropyleneGlycol_00006.logおよびPropyleyeGlycol_00006.fchkに出力されています。またls1ファイルの「ORIGINAL STERIC E」は、Gaussianの計算で得られた全エネルギーをkcal/molに変換したもの（最安定構造の全エネルギーは-269.56133526 a.u.なので、-269.56133526×627.5095=-169152.2987）です。尚この計算では、Gaussian 16 Rev. A.03を使用しています。
PropyleneGlycol.ls1：
====================================================================================
CONF. ORIGINAL DISTRI- NO.
NO. ID STERIC E DELTA E BUTION INIT. REOPT. NEG.
------------------------------------------------------------------------------------
1 00000006 -169152.2987 0.0000 21.2884 T* F+ 0
2 00000017 -169152.2709 0.0278 20.3129 T* F+ 0
3 00000019 -169152.1497 0.1490 16.5548 T* F+ 0
4 00000011 -169151.9408 0.3579 11.6356 T* F+ 0
5 00000001 -169151.9312 0.3675 11.4487 T* F+ 0
6 00000008 -169151.6798 0.6189 7.4899 T* F+ 0
7 00000012 -169151.4409 0.8578 5.0044 T* F+ 0
8 00000018 -169151.1376 1.1611 2.9993 T* F+ 0
9 00000007 -169150.9622 1.3365 2.2310 T* F+ 0
10 00000016 -169149.5319 2.7668 0.1995 T* F+ 0
11 00000020 -169149.3686 2.9301 0.1515 T* F+ 0
12 00000009 -169149.2519 3.0468 0.1244 T* F+ 0
13 00000022 -169149.2130 3.0857 0.1165 T* F+ 0
14 00000021 -169149.0687 3.2300 0.0913 T* F+ 0
15 00000023 -169149.0383 3.2604 0.0868 T* F+ 0
16 00000024 -169149.0232 3.2755 0.0846 T* F+ 0
17 00000013 -169148.9887 3.3101 0.0798 T* F+ 0
18 00000003 -169148.8937 3.4050 0.0680 T* F+ 0
19 00000010 -169147.6211 4.6776 0.0079 T* F+ 0
20 00000014 -169147.4593 4.8394 0.0060 T* F+ 0
21 00000004 -169147.3564 4.9423 0.0051 T* F+ 0
22 00000015 -169147.3371 4.9616 0.0049 T* F+ 0
23 00000005 -169147.3181 4.9806 0.0048 T* F+ 0
24 00000002 -169147.2083 5.0904 0.0040 T* F+ 0
------------------------------------------------------------------------------------
MINIMUM ENERGY: -169152.2987 KCAL/MOL
AVERAGE ENERGY: -169151.9983 KCAL/MOL (BASED ON THE BOLTZMANN LAW)
====================================================================================
Gaussian入出力ファイルの一覧：
PropyleneGlycol_00001.fchk PropyleneGlycol_00008.fchk PropyleneGlycol_00015.fchk PropyleneGlycol_00001.gjf PropyleneGlycol_00008.gjf PropyleneGlycol_00015.gjf PropyleneGlycol_00001.log PropyleneGlycol_00008.log PropyleneGlycol_00015.log PropyleneGlycol_00002.fchk PropyleneGlycol_00009.fchk PropyleneGlycol_00016.fchk PropyleneGlycol_00002.gjf PropyleneGlycol_00009.gjf PropyleneGlycol_00016.gjf PropyleneGlycol_00002.log PropyleneGlycol_00009.log PropyleneGlycol_00016.log PropyleneGlycol_00003.fchk PropyleneGlycol_00010.fchk PropyleneGlycol_00017.fchk ・・・・・
Propylene glycolをMMFF94sで配座探索を行う（Search Limit: 7（SEL=7））と、21個の配座異性体が得られます。これらをB3LYP/6-31G(d)レベルで再度最適化すると、MMFF94sで5番目に安定な配座が最も安定になります（下図）。
Propylene glycolのB3LYP/6-31G(d)レベルでの最安定構造
また、MMFF94sで10番目に安定な構造をB3LYP/6-31G(d)で最適化すると他の構造と同一になる（下図）ため、最終的にB3LYP/6-31G(d)で得られた配座は20個となりました。
MMFF94sで10番目に安定な構造（左）をB3LYP/6-31G(d)で再度最適化すると、異なる配座異性体（右）として収束する
先にMMFF94s力場を用いて探索し、得られた配座異性体をGaussianで再度最適化した結果と、Gaussianのみを利用して配座探索を行った結果を比較しますと、次の図に示す4つの構造が新たに見つかったことがわかりました。
B3LYP/6-31G(d)による配座探索で新たに得られた構造
【ファイルを用いたGaussianの計算設定】
Link0コマンドやルートセクションなどを個別のキーワードで設定する以外に、ファイルで設定する方法があります。ここでは、ファイルを用いたGaussianの計算設定について説明します。
＊例1
Gaussian実行用の入力ファイルの内容を下記とします。
%mem=32MW %nprocshared=2 #HF/3-21G opt test of HF/3-21G calculation 0 1 molecular specification ...
分子構造データより上の設定を、gaussian.hdrファイルとして保存します。
gaussian.hdrファイル：
%mem=32MW %nprocshared=2 #HF/3-21G opt test of HF/3-21G calculation 0 1
[Interfaceから実行する場合]
詳細設定ダイアログの「External Program」ダイアログで、「Header File:」のチェックボックスにチェックを入れます。
次に、「Select...」をクリックして、「gaussian.hdr」ファイルを選択してください。gaussian.hdrの内容を元に、Gaussianの計算設定がなされます。
[コマンドラインから実行する場合]
PropyleneGlycol.iniファイルは、下記のように変更します。
PropyleneGlycol.iniファイル
CONFLEX
SEL=5
EXTERNAL_PROGRAM=GAUSSIAN
EXT_GAU_HEADER_FILE=(gaussian.hdr)
EXT_PREOPT=OFF
EDIF_HARD=1.0D-3
PropyleneGlycol.molとPropyleneGlycol.ini、gaussian.hdrをフォルダに格納し、コマンドを実行します。
＊例2
キーワード「Gen」を指定して基底関数を設定する時のように分子構造データの後に設定を記述する場合も、別ファイルに記述して読み込ませることが可能です。Gaussian実行用の入力ファイルの内容を下記とします。
%mem=20MW %nprocshared=2 # B3LYP/gen Opt TEST OF B3LYP/6-31G(d) CALCULATION 0 1 molecular specification ... H C O 0 6-31G(d) ****
分子構造データより上の設定を、gaussian.hdrファイルとして保存します。
gaussian.hdrファイル：
%mem=32MW %nprocshared=2 #HF/3-21G opt test of HF/3-21G calculation 0 1
分子構造データより下の設定を、gaussian.ftrファイルとして保存します。
gaussian.ftrファイル：
H C O 0 6-31G(d) ****
[Interfaceから実行する場合]
詳細設定ダイアログの「External Program」ダイアログで、「Header File:」のチェックボックスにチェックを入れます。
次に、「Select...」をクリックして、「gaussian.hdr」ファイルを選択してください。
また、「Footer File:」のチェックボックスにチェックを入れます。
次に、「Select...」をクリックして、「gaussian.ftr」ファイルを選択してください。gaussian.hdrとgaussian.ftrの内容を元に、Gaussianの計算設定がなされます。
[コマンドラインから実行する場合]
PropyleneGlycol.iniファイルは、下記のように変更します。
PropyleneGlycol.iniファイル
CONFLEX SEL=5 EXTERNAL_PROGRAM=GAUSSIAN EXT_GAU_HEADER_FILE=(gaussian.hdr) EXT_GAU_FOOTER_FILE=(gaussian.ftr) EXT_PREOPT=OFF EDIF_HARD=1.0D-3
PropyleneGlycol.molとPropyleneGlycol.ini、gaussian.hdr、gaussian.ftrをフォルダに格納し、コマンドを実行します。
【ジョブスケジューラを使ったGaussianの実行（Linuxのみ）】
CONFLEXから実行されるGaussianジョブを、ジョブスケジューラー経由で実行する方法を説明します。ここでは、ジョブスケジューラーはUniva Grid Engineとします。
[コマンドラインから実行する場合]
PropyleneGlycol.iniファイルは、下記のように変更します。
PropyleneGlycol.iniファイル
CONFLEX SEL=5 EXTERNAL_PROGRAM=GAUSSIAN EXT_JOB_COM=("qsub -sync yes”) EXT_JOB_FILE=qsub.sh EXT_GAU_LINK0=(mem=128MW) EXT_GAU_LINK0=(nprocshared=2) EXT_GAU_ROUTE=("B3LYP/6-31G(d) opt”) EXT_GAU_COMMENT=("Propylene Glycol B3LYP/6-31G(d) optimization”) EXT_GAU_CHARGE=0 EXT_GAU_SPIN=1 EXT_PREOPT=OFF EDIF_HARD=1.0D-3
「EXT_JOB_COM=」は、ジョブスケジューラーの実行コマンド（qsub -sync yes）を設定しています。ジョブスケジューラーを用いる場合は必ずインタラクティブに実行するオプション（Univa Grid Engineでは「-sync yes」）を加えてください。また、「EXT_JOB_FILE=」は、ジョブスケジューラーに対する設定などが記述されたファイルを指定します。
qsub.shファイル：
#!/bin/bash #$ -S /bin/bash #$ -N g16_from_conflex #$ -cwd #$ -V #$ -j y #$ -q all.q #$ -pe sme 2 export g16root=/usr/local
これにより、Gaussianのジョブは「qsub -sync yes qsub.sh」で実行され、ジョブスケジューラー経由でのジョブ実行となります。なお、qsub.shにはGaussianの実行コマンドがありませんが、CONFLEXが自動的に追加します。qsub.shは、PropyleneGlycol.molとPropyleneGlycol.iniと同じフォルダに格納してください。
【一時的な環境設定の有効化】
Gaussianの環境設定を計算実行時のみ有効にする方法を説明します。
Windowsの場合
[Interfaceからの実行]
詳細設定ダイアログのをクリックし、開いたダイアログに「EXT_JOB_FILE=setup.txt」を追加します。
[コマンドラインから実行する場合]
PropyleneGlycol.iniファイルは、下記のように変更します。
PropyleneGlycol.iniファイル
CONFLEX
SEL=5
EXTERNAL_PROGRAM=GAUSSIAN
EXT_JOB_FILE=setup.txt
EXT_GAU_LINK0=(mem=128MW)
EXT_GAU_LINK0=(nprocshared=2)
EXT_GAU_ROUTE=("B3LYP/6-31G(d) opt”)
EXT_GAU_COMMENT=("Propylene Glycol B3LYP/6-31G(d) optimization”)
EXT_GAU_CHARGE=0
EXT_GAU_SPIN=1
EXT_PREOPT=OFF
EDIF_HARD=1.0D-3
ここで「EXT_JOB_FILE=」は設定ファイル名を表します。設定ファイルsetup.txtの内容を以下のようにすると、CONFLEXからGaussianを呼び出す毎に、このファイルの内容とGaussian実行コマンドを含むbatファイルを作成し、計算を実行します。
setup.txtファイル：
set GAUSS_EXEDIR=C:\G16W
また、をクリックし開いたダイアログやiniファイルに「EXT_GAU_SCRDIR=D:\Scratch」を加えると、中間ファイルの保存先をD:¥Scratch以下に変更できます。加えない場合は、「%GAUSS_EXEDIR%¥Scratch」ファルダ−に保存されます。
Mac, Linuxの場合
[Interfaceから実行する場合]
詳細設定ダイアログのをクリックし、開いたダイアログに「EXT_JOB_COM=(“sh")」と「EXT_JOB_FILE=setup.sh」を追加します。
[コマンドラインから実行する場合]
PropyleneGlycol.iniファイルは、下記のように変更します。
PropyleneGlycol.iniファイル
CONFLEX SEL=5 EXTERNAL_PROGRAM=GAUSSIAN EXT_JOB_COM=("sh”) EXT_JOB_FILE=setup.sh EXT_GAU_LINK0=(mem=128MW) EXT_GAU_LINK0=(nprocshared=2) EXT_GAU_ROUTE=("B3LYP/6-31G(d) opt”) EXT_GAU_COMMENT=("Propylene Glycol B3LYP/6-31G(d) optimization”) EXT_GAU_CHARGE=0 EXT_GAU_SPIN=1 EXT_PREOPT=OFF EDIF_HARD=1.0D-3
ここで「EXT_JOB_COM=」は実行コマンド、「EXT_JOB_FILE=」は設定ファイル名をそれぞれ表します。設定ファイルsetup.shの内容を以下のようにすると、CONFLEXからGaussianを呼び出す毎に、このファイルの内容とGaussian実行コマンドを含むスクリプトファイルを作成し、shコマンドで計算を実行します。
setup.shファイル：
export g16root=/usr/local export GAUSS_SCRDIR=/tmp