Gaussianで量子化学計算を行うためのインプットファイルをカテゴリーごとにまとめました。
chloroethaneやiodoethaneを題材としたインプットファイルになっており、キーワードに関する説明も記載しています。
（３次元構造はAvogadroを使用して作成しています。）
各インプットファイルのキーワードは組み合わせて使うこともできますので、この記事で紹介した組み合わせ以外も必要に応じて試してみてください。
キーワードの使い方や、キーワード同士を組み合わせた応用例など、Gaussianを使い始めたばかりの方のサポートのなれば幸いです。
・Gaussianのインプットファイルの書き方
・構造最適化やrelaxed scanなど有機化合物の構造や反応機構の計算に必要なインプットファイル
・計算がうまく進まないときに試せるキーワードの紹介
・ECPや半経験量子化学などDFT計算以外のインプット
Gaussianのインプットファイルの作成と基本的なキーワードについては以下の記事をご覧ください。
Avogadroを用いた構造作成については以下の記事をどうぞ
構造最適化関連
例1) 構造最適化+振動計算, トルエン溶媒（SMD）
%nprocshared=4 %mem=10GB %chk=opt_EtCl.chk # opt freq wb97xd/Def2SVPP scrf=(smd,solvent=toluene) EtCl_optimization 0 1 C -1.80378 0.68701 -0.20289 Cl -0.02447 0.65684 -0.28348 H -2.19182 1.35309 -1.00208 C -2.26025 1.19563 1.15574 H -2.19182 -0.34014 -0.36821 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 !ここでインプットは終わりです。
optとfreqを両方指定すると、構造最適化を行った後にその構造に対する振動計算が行われます。
振動計算は最適化された構造（局所安定構造や遷移状態）に対して同じ計算レベルで行う必要があるうえ、虚振動の数は構造が遷移状態なのか局所安定構造なのかを判定するのにも必須です。
したがって、構造最適化するのであれば振動計算も同時に行うよう、optとfreqを使うのがよいでしょう。
溶媒モデルとしてSMDモデルのトルエンを指定するためにはscrfキーワードに対してさらに追加のキーワードを入れて、
scrf=(smd,solvent=toluene)
のようにします。
他の溶媒モデルの使い方については以下の公式サイトを参照してください。
SCRF: Gaussian.com
例2) 収束基準を厳しくした構造最適化+振動計算
%nprocshared=4 %mem=10GB %chk=opt_EtCl_tight.chk # opt=tight freq wb97xd/Def2SVPP scrf=(smd,solvent=toluene) EtCl_optimization_tight 0 1 C -1.80378 0.68701 -0.20289 Cl -0.02447 0.65684 -0.28348 H -2.19182 1.35309 -1.00208 C -2.26025 1.19563 1.15574 H -2.19182 -0.34014 -0.36821 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 !ここでインプットは終わりです。
例1)と同じ計算ですが、optオプションにtightを追加することで収束基準（構造最適化の終了を判断する基準）を小さくして、より厳密な構造最適化をしています。
opt=vtight
とすると、さらに厳しい収束基準を適用することができます。
例3) 遷移状態(TS)最適化
%nprocshared=4 %mem=10GB %chk=opt_EtCl_ts.chk # opt=(calcfc,ts,maxstep=10) freq m062x/6-31G(d,p) EtCl_ts 0 1 C -1.80378 0.68701 -0.20289 Cl 0.54435 0.67342 -0.37105 H -1.93628 1.42627 -0.90576 C -2.26025 1.19563 1.15574 H -1.97576 -0.30573 -0.27096 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 !ここでインプットは終わりです。
TSの最適化では、構造最適化の前に振動計算を行い、力の定数を求める必要があります。
したがってoptキーワードに対して、
・遷移状態の最適化を指定する「ts」
・入力構造に対して構造最適化前に力の定数を計算することを指定する「calcfc」
が必須となります。
また、遷移状態周辺では緻密な構造最適化が必要になることが多いため、構造最適化の1ステップにおいて原子を動かす大きさをデフォルトより小さくするのがおすすめです。
この値を変更するためには、optキーワードに対して
・maxstep=10（デフォルト値は30）
とすると適用可能です。
必要に応じて5にするなど変更するとよいです。
例4) TS最適化（力の定数読み込み） + NoEigenTest
%nprocshared=4 %mem=10GB %chk=opt_ts_noeig.chk %oldchk=opt_EtCl_ts.chk # opt=(readfc,ts,maxstep=10,noeigentest) freq m062x/6-31G(d,p) EtCl_ts_noeig 0 1 C -1.80378 0.68701 -0.20289 Cl 0.54435 0.67342 -0.37105 H -1.93628 1.42627 -0.90576 C -2.26025 1.19563 1.15574 H -1.97576 -0.30573 -0.27096 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 !ここでインプットは終わりです。
このインプットファイルで注目したいキーワードは
・opt=readfc
・opt=noeigentest
の2点です。
これらは独立して使用できますので、どちらか一方のみを使うことも多いです。
opt=readfcは、TS最適化で必須となる力の定数を、opt=calcfcのときのように最初に計算するのではなく、chkファイルから読み込むことを指定しています。
すなわち、一度freq等の計算により振動計算を行った計算結果がある場合にのみ利用できる方法で、そのときに作成されたchkファイルを用いることとなります。
読み込むファイルの指定は%chkでもよいのですが、このようにするとこれから行う計算結果が元のchkファイルに書かれてしまいます。
管理人は元のファイルは保持しておきたい主義なので、Link0コマンドの
%oldchk=opt_EtCl_ts.chk
で古いチェックポイントファイル（opt_EtCl_ts.chk）を読み込み元としたうえで、
%chk=opt_ts_noeig.chk
で今回の計算ジョブによりできるchkファイル名を指定しています。
次に、opt=noeigentestについてです。
例3のようなインプットファイルを用いたTS最適化においては、しばしば以下のようなエラーに遭遇して構造最適化が止まってしまいます。
Gaussianの構造最適化では、各ステップごとにエネルギーの曲率を求め、虚振動の数が指定した数と一致しそうかをチェックする機能（EigenTest）があります。
opt=noeigentestはこのEigenTestを行わなくするオプションです。
上記エラーメッセージは、虚振動が1つ（TSは虚振動1つなので）であることを指定しているが、EigenTestの結果、虚振動が2つとなりそうであるというようなことを言っています。
本当に2次の鞍点（虚振動2つ）の構造になってしまうのであれば問題なのですが、単に構造最適化が不十分であり、そのまま最適化が進んでいけば解消されることもあります。
そのため、このEigenTestの機能を切ってしまえば最適化が進みます。
振動計算の結果を注視する必要はありますが、大きい系においてはC-C結合の回転に起因する小さい虚振動が残りやすい気がします。
上記エラーに遭遇するようであれば使用してみてください。
反応経路関連
例5) relaxed scan
%nprocshared=4 %mem=10GB %chk=scan_C_Cl.chk # opt=modredundant wb97xd/Def2SVPP scrf=(smd,solvent=toluene) scan_C_Cl 0 1 C -1.80378 0.68701 -0.20289 Cl -0.02447 0.65684 -0.28348 H -2.19182 1.35309 -1.00208 C -2.26025 1.19563 1.15574 H -2.19182 -0.34014 -0.36821 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 B 1 2 S 10 0.1 !ここでインプットは終わりです。
chloroethaneのC-Cl結合を0.1Å刻みで伸ばし、10個の最適化構造を計算する（10 steps）インプットファイルです。
C-Cl結合は各ステップごとに指定した長さに固定化されますが、それ以外の結合長、結合角、二面角はその状態で最適化されます。
relaxed scanを行うためには
opt=modredundant
を指定したうえで、付加情報セクションにスキャンする結合に関する情報を記載します。
今回の例では
B 1 2 S 10 0.1
となっています。
左から、順にその値の意味を説明します。
・B:
指定する座標が結合長（bond）であることを意味します。
他にA（Angle）とD（Dihedral）をよく使います。
・1 2:
分子指定セクションに記載した「1」番目と「2」番目の原子の結合長であることを意味します。
この場合ではC-Cl結合を指定したことになります。
Angleのときは値が3つ、Dihedralのときは値が4つ必要です。
・S:
指定した座標に対してrelaxed scanを行うことを指定。
Fとするとその座標の値を固定することができます。
・10:
relaxed scanのときに指定する値で、10 stepsスキャンすることを意味します。
・0.1:
同じくrelaxed scanのときに指定する値で、1 stepごとに+0.1 Å変化させることを意味します。減少させたいときは-0.1のようにマイナスをつければよいです。
ModRedundantオプションでは他にも様々なことができますので、詳細については以下の公式サイトを見てみてください。
ModRedundant: Gaussian.com
例6) 角度の固定を加えたrelaxed scan
%nprocshared=4 %mem=10GB %chk=scan_with_const.chk # opt=modredundant wb97xd/Def2SVPP scan_C_Cl_with_constraint 0 1 C -1.80378 0.68701 -0.20289 Cl -0.02447 0.65684 -0.28348 H -2.19182 1.35309 -1.00208 C -2.26025 1.19563 1.15574 H -2.19182 -0.34014 -0.36821 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 B 1 2 S 10 0.1 D 8 4 1 3 F !ここでインプットは終わりです。
基本的には例5)と変わりませんが、付加情報セクションに以下の2行が記載されています。
B 1 2 S 10 0.1 D 8 4 1 3 F
1行目は1番目の原子と2番目の原子の結合距離（B 1 2 )のスキャン（S）
2行目は8番目、4番目、1番目、3番目の原子からなる二面角(D 8 4 1 3)の固定（F）
をそれぞれ意味しており、スキャンに伴いこの二面角が変化しないように意図されています。
反応により複数の座標が変化しうる場合、片方を固定化するrelaxed scanが有効な場合があります。
例7) 固有反応座標（IRC）計算 + 構造情報をchkファイルから読み込み
%nprocshared=10 %chk=irc.chk %oldchk=opt_ts.chk # irc=(rcfc,maxcyc=50,maxpoints=30,stepsize=30) wb97xd/6-31+g(d) scrf=(smd,solvent=chloroform) geom=allcheck irc
反応の遷移状態が得られた場合、その遷移状態が確かに望みの出発物と生成物を繋ぐ遷移状態であるかを確認する必要があります。
この目的で、遷移状態の構造を起点として虚振動の方向に沿った反応経路（IRC: 固有反応座標）の計算を行います。
IRC計算では、TSの構造と、その構造における振動計算の結果が必要です。
振動計算の結果はTSの構造最適化時に通常は行うはず（虚振動が1つであることを確認する）ですので、改めて計算する必要はないことがほとんどです。
chkファイルから振動計算の結果を読み込む際にはoptのときと似ていますが、ircキーワードに対してrcfcのキーワードを追加して
irc=rcfc
とすればよいです。
また、３次元構造も、以下のキーワードを用いてchkファイルから取得することができます。
geom=allcheck
このキーワードについては、IRC計算だけでなく様々な計算で使用できます。
（後述のエネルギー計算でも用いています。）
%chkや%oldchkについては、例4)で解説しています。
計算レベル関連
例8) トリプルゼータ基底でのエネルギー計算（chkファイルからの構造読み込み）
%nprocshared=4 %mem=10GB %chk=energy_high.chk %oldchk=opt_EtCl.chk # m062x/def2tzvp scrf=(smd,solvent=toluene) geom=allcheck scrf=(smd,solvent=toluene) energy_highlevel
ある程度大きい分子を対象とした場合、構造最適化計算をダブルゼータ基底（6-31G(d)やDef2-SV(P)）で行ってから、その構造に対して一点エネルギー計算をトリプルゼータ基底（6-311G(d)やDef2-TZVP）を行うことが多いです。
このような場合は上記インプットファイルで計算を実行できます。
%oldchkで指定したchkファイルからダブルゼータ基底で構造最適化した構造を読み込むため
geom=allcheck
を入れています。
ギブズエネルギー（G）を求めたい場合はダブルゼータ基底で計算した振動計算の結果から得られる零点エネルギー（ZPE）を足すようにしましょう。
例9) 分散力補正の追加
%nprocshared=4 %mem=10GB %chk=energy_high_D3.chk %oldchk=opt_EtCl.chk # b3lyp/def2tzvp scrf=(smd,solvent=toluene) geom=allcheck scrf=(smd,solvent=toluene) EmpiricalDispersion=GD3BJ EtCl_GD3BJ
DFTの弱点である長距離相互作用の記述を改善するため、最近では分散力補正を行うのが一般的になってきています。
上記インプットファイルふはB3LYP汎関数での計算において、
EmpiricalDispersion=GD3BJ
とすることで、GD3BJ分散力モデルを適用したインプットファイルです。
一点計算だけでなく構造最適化などの計算でも同様に適用できます。
GD3BJ以外にもキーワードとして
#Petersson-Frisch分散力モデル EmpiricalDispersion=FD #Grimme D3分散力モデル EmpiricalDispersion=GD3
が使えます。
汎関数自体にすでに分散力補正項が加えられている場合もあるので、そのような場合は使用しない方がよいのかなと思います。
例10) MP2法での計算
%nprocshared=4 %mem=10GB %chk=energy_MP2.chk %oldchk=opt_EtCl.chk # MP2/def2tzvp scrf=(smd,solvent=toluene) geom=allcheck scrf=(smd,solvent=toluene) EtCl_MP2
計算手法としてDFT計算ではなく、2次のMøller–Plesset摂動法（MP2）を用いるインプットファイルです。
DFTよりも少しコストの高い計算になります。
3次の摂動法であればMP3、4次の摂動法であればMP4となり、コストと精度が上がっていきます。
（MP3はあまり使われないようです。）
基底関数の指定はDFTと同様です。
例11) 半経験的量子化学（PM7）での構造最適化
%nprocshared=4 %mem=10GB %chk=opt_PM7.chk # opt freq PM7 EtCl 0 1 C -1.80378 0.68701 -0.20289 Cl -0.02447 0.65684 -0.28348 H -2.19182 1.35309 -1.00208 C -2.26025 1.19563 1.15574 H -2.19182 -0.34014 -0.36821 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 !ここでインプットは終わりです。
半経験的手法のうち、PM7法を使った構造最適化計算です。
半経験的量子化学では基底関数の指定は必要ありません。
PM7以外の方法については以下の公式サイトを参照してください。
Semi-Empirical Methods: Gaussian.com
例12) 有効内核ポテンシャル（ECP）の追加 （ダブルゼータ基底）
%nprocshared=4 %mem=10GB %chk=opt_EtCl.chk # opt freq wb97xd/genecp EtI_ecp_double_zata 0 1 C -1.80378 0.68701 -0.20289 I 0.54435 0.67342 -0.37105 H -1.93628 1.42627 -0.90576 C -2.26025 1.19563 1.15574 H -1.97576 -0.30573 -0.27096 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 I 0 lanl2dz **** H C 0 6-31G(d) **** I 0 lanl2dz
第5周期以降の原子は、電子の数が多くその扱いも複雑になります。
このような原子についての計算コストを抑えたり、計算の安定性を向上させる目的で、有効内核ポテンシャル（ECP）が適用されます。
ECPの特徴は以下の通りです。
・内核の電子を一つずつ扱わずにポテンシャルのように表現する
・価電子は通常の基底関数と同様に扱う
上記インプットファイルではヨウ素原子に対して、ダブルゼータ基底相当のECPであるLanl2DZを適用し、他の原子についてはダブルゼータ基底の6-31G(d)を適用しています。
ECPと他の元素の基底のレベルは合わせたほうがよいです。
設定方法ですが、基底関数のキーワードとして
wb97xd/genecp
のように「genecp」キーワードを加えます。
その後付加情報セクションに以下の記述をします。
I 0 lanl2dz **** H C 0 6-31G(d) **** I 0 lanl2dz
18行目から23行目までは、基底関数の設定になります。
それぞれの原子の後の数字は「0」を指定すると全ての原子、「1」を指定するとその元素の中の1番目ものという指示になります。
ヨウ素の基底関数はLanl2DZ、水素と炭素の基底関数は6-31G(d)ということになります。
「****」は意味のある文字列なので変更しないようにしてください。
25行目、26行目までは、ECPの設定になります。
基底関数のセクションと同じ「Lanl2DZ」が記載されていてややこしいですが、こちらがECPの定義になりますので、忘れてはいけません。
Lanl2DZやSDDをECPとして用いる場合は基底関数と同じキーワードを記載するのみでECPの適用が可能です。
SDDを使うインプットファイルは次の例13で示します。
例13) トリプルゼータ基底でのECP (SDDの使用)
%nprocshared=4 %mem=10GB %chk=opt_ecp_triple.chk # opt freq wb97xd/genecp EtI_ecp_triple_zata 0 1 C -1.80378 0.68701 -0.20289 I 0.54435 0.67342 -0.37105 H -1.93628 1.42627 -0.90576 C -2.26025 1.19563 1.15574 H -1.97576 -0.30573 -0.27096 H -1.88077 2.22565 1.32472 H -1.88077 0.52978 1.95957 H -3.36977 1.21060 1.19572 I 0 SDD **** H C 0 def2tzvp **** I 0 SDD
例12で、ECPと他の元素の基底関数のレベルは合わせて方が良いと述べました。
ダブルゼータ基底に対応するLanl2DZではなく、トリプルゼータ基底を適用したい場合はSDDを用いるとよいです。
キーワードがSDDになっているくらいで例12と同様です。
水素と炭素の基底関数がトリプルゼータ基底のDef2-TZVPになっていることを確認してください。
終わりに
Gaussianのインプットファイルについて、主なものをまとめて記載しました。
これ以外のインプットファイルについても今後追加して公開するかもしれません。
以下の公式サイトから、今回紹介しなかったキーワードについても知ることができますので、興味のある方はどうぞ。
Gaussian 16 Users Reference: Gaussian.com
コメント・質問等はX (@chemmodelcomp)にてご連絡ください。
