4. Gaussianの入力データ

4.1. 入力データの構成

入力データの並びは、以下のように5つのセクションに分かれています。

  • LinkOコマンド
    「%」で始まる行で.chkや.rwfなどのスクラッチファイル、使用メモリ量やCPU数を指定します。

  • ルートセクション
    「#」で始まる行で、理論の選択、計算内容の設定を行います。複数行に分けることができるため、セクションの終わりを示す空白行が必須です。

  • タイトルセクション
    コメントを付けます。 このセクションも終わりを示す空白行が必須です。

  • 分子指定セクション
    分子の電荷、多重度の後にZ行列やCartesian座標を用いて分子構造を記述します。 このセクションも終わりを示す空白行が必須です。

  • 付加情報
    計算内容によっては、付加情報が必要になります。その場合には、ここに記述します。このセクションは内容によって空白行がいらない場合があります。 下記の例は、水分子の最適化構造をHF/6-31G*レベルで求める際のGaussianの入力データです。

入力データ例(h2o.dat)

1   %CPU=0-1
2   %Mem=48Mb
3   %Chk=h2o
4   # HF/6-31G* Opt=Z-Matrix Test
5  
6   h2o
7  
8   0 1
9   O
10  H,1,r1
11  H,1,r1,2,a1
12  
13  r1=0.958
14  a1=104.5

入力データ中で、タイトル行のコメント、ファイル名以外は大文字・小文字の区別がありません。 上の例では、キーワードを見やすいように大文字を使っている部分があります。 また、漢字フォントやエスケープなどの表示できない文字は使えません。全角の空白が紛れこんでいても気が付かない場合があります。 この例をもとに、各セクションの説明をしましょう。

4.1.1. Link0コマンド

Link0コマンドは、「%」で始まる行で、一項目ごとに一行ずつ記述します。 また、各行の順番は任意で、省略した場合はデフォルト値が使われます。
1行目の%CPUでは、使用するコア番号を指定します。 f_nodeの28コアを使用する場合は %CPU=0-27 と指定してください。 f_node以外を使用する場合は、バッチジョブまたはインタラクティブジョブ投入により動的に割り当てられるコア番号を指定してください。 コア番号は numactl -s コマンドのphyscpubindで確認できます。 例えば、q_nodeの7コアを確保し numactl -s コマンドで以下のように表示された場合は、 %CPU=7-13 と指定します。 (28から55まではハイパースレッディングによる論理コアであり、指定すると計算効率が低下するため指定しないことを推奨します)

$ numactl -s | grep physcpubind
physcpubind: 7 8 9 10 11 12 13 35 36 37 38 39 40 41

2行目の%Memでは、使用するメモリを指定します。この例では、48Mbにしてあります。 単位を省略するとGaussianword(=8Bytes)になります。出力リストにはGaussianword単位で、使用メモリ量が出力されます。 単位の指定には、他にMb、Mw、Gb、などが使えます。デフォルトは、100MW(=800MB)です。
3行目の%Chkでは、チェックポイントファイルのファイル名を指定します。 この例の場合、実行ディレクトリにh2o.chkというファイルが作成され、計算結果の要約がバイナリ形式で保存されます。 チェックポイントファイルは、他のグラフィックソフトで構造や分子軌道を描いたり、さらに進んだ計算を引き続き行うときに使うことができます。 省略した場合は、ファイルが保存されません。
使用するメモリの指定は重要で、場合によってはメモリ不足で計算が実行できないときがあります。 メモリの指定とともに、後述する、ハードディスク上の.rwfファイルのサイズを制限するキーワード、MaxDiskの指定によって使用するアルゴリズムが変わる場合があります。 ただし、使用メモリ・ディスク使用量をいたずらに大きくすると、必ずしも計算を速くするわけではありません。 計算する系の大きさ、使用する基底関数、そして計算手法によってメモリ及びハードディスクサイズを変えなければなりません。 限られた計算機資源を有効に使うには、無闇な計算はせず、経験を積んでいただく必要があります。

4.1.2. ルートセクション

4行目の「#」で始まる行からルートセクションが始まります。 この例では一行ですが、数行に渡って記述することができるので、セクションの終わりを示す空白行が必須です。
5行目の空白行がルートセクションの終わりを示します。
このセクションで、計算に用いる手法のキーワードを記します。 HF/6-31Gとあるのは、基底関数として6-31Gを用い、Hartree-Fockによる計算を行うことを指定しています。 Optは構造最適化、Testは結果の要約をアーカイブとして出力しないことをそれぞれ意味するキーワードです。

Gaussianでは、構造最適化をする際のデフォルトがredundantintemalcoordinate(RIC)という新しい手法になっています。 RICではデカルト座標を用い、すべての原子間の結合情報を判定する手法です。 この手法の利点は、Z行列などで明確に指定できない結合や結合角についても考慮できるため、環構造を持つ化合物などの構造最適化の収束がよいということです。 この例でOpt=Z-Matrixを単にOptとすると、RICでの計算を行うことになります。 ただし、通常の分子ではZ行列の方が収束がよいことも多く(計算時間が短い)、慣れてきたら座標の入力方法を使い分けて下さい。

4.1.3. タイトルセクション

6行目がタイトル行です。この行はコメント用に自由に使えます。 この行も複数行になっても構いません。このセクションも、7行目のように空白行で終わりを示します。

4.1.4. 分子指定セクション

8行目から分子定義のセクションが始まります。 この行の「01」という数値は、全電荷とスピン多重度をそれぞれ示しており、この場合は中性分子で一重項です。 これに続く行は分子の構造データであり、CartesianやZ行列などを用いて構造を記述します。 この例では、Z行列の形式で記述して、r1とa1の2つのパラメータを使って2つの0-H結合を等価に扱っています(C2v対称)。 Z行列の作り方は、後述します。 12行目の空白行でこのセクションの終わりを示しています。
この例では、Z行列でパラメータを使っているので、付加情報としてr1とa1を定義するセクションが必要になります。 15行目の空白行で、このセクションの終わりを示します。

4.1.5. 付加情報

計算内容によっては、この外にも付加情報を必要とする時があります。 その場合は、必要な数だけ付加情報セクションを繰り返します。 この時、空白行で終わりを示す必要があるかどうかは、キーワードによります。
テスト用入力データは/apps/t3/sles12sp2/isv/gaussian/g16/B01/g16/tests/com/というディレクトリに1152個の入力例があります。 grepコマンドなどを使って関連する入力例を探すと便利です。

4.2. 理論の選択

Gaussianに含まれる代表的な理論のキーワードを下表に示します。 これらはほんの一部で、他にもたくさんありますし、マニュアルに記載されていても現在のRevisionに含まれていないものがあります。 将来のRevisionアップで新たに計算できる手法が加わることがあります。 また、Gaussianから、分子をいくつかのレイヤーに分けて異なるレベルの理論を混在させることができる、ONIOMモデルが使えるようになりました。

代表的手法のキーワード

キーワード 手法
HF Hartree-Fock(HF)Self-Consistent Field(SCF)計算
CASSCF active電子及び軌道で構成される多配置参照SCF計算
GVB Generalized Valence Bond
MP2 2次の摂動論
MP3 3次の摂動論
MP4(SDTQ) 4次の摂動論
MP5 5次の摂動論
CIS 1電子励起の配置間相互作用
CISD 1電子及び2電子励起の配置間相互作用
QCISD Quadratic CISD
CCSD 1電子及び2電子励起のcoupled cluster法
Slater LSD交換汎関数による密度汎関数理論(DFT)計算
Xalpha Xalpha交換汎関数によるDFT計算
HFB Becke1988交換汎関数によるDFT計算
VWN Vosko、Wilk、Nusair相関汎関数によるDFT計算
VWN5 Vosko、Wilk、Nus&ir相関汎関数VによるDFT計算
LYP Lee、Yang、Parr相関汎関数によるDFT計算
PL Perdew1981相関汎関数によるDFT計算
P86 Perdew1986相関汎関数によるDFT計算
PW91 Perdew/Wang1991相関汎関数によるDFT計算
Becke3 Becke3-parameter汎関数によるDFT計算

HF以外は、ほとんどが電子相関の効果を取り入れるための手法でそれぞれに特徴があり、必要な計算機資源の大きさも大きく異なります。 また、高いレベルの手法にはそれに応じた精度の高い基底関数を選択する必要があります。 従って、いきなり高度な理論を使って計算するのは避けて、対象によって使い分ける必要があります。 最も良く使われるのは、HF、MP2あるいはB3LYPで構造最適化を行い、得られた構造を元にさらに基底関数を大きくしたり、レベルの高いMP4、CISD、QCISD、CCSDなどによる一点計算をすることです。 電子相関の効果を取り入れた手法では、内殻電子の励起は考えないfrozen-core(FC)近似がデフォルトであることに注意しましょう。 内殻電子の寄与も取り入れたいときは、例えばMP2=fullのようにします。
最近、盛んに使われるようになった密度汎関数理論(DFT)は、従来の電子相関を取り入れた手法と比べ、計算時間、使用メモリ量、ディスク使用量ともに少なくてすみ、しかもMP2と同程度かそれ以上の計算精度を与えることが多いことが知られています。 ただし、現在、まだ従来の分子軌道法や実測値との比較が盛んに行われている段階ですので、使用に当たっては注意が必要です。 いくつかの問題点も指摘されています。 できるだけ自分の研究に近く、しかも新しい論文を探していくつか読んでみることをおすすめします。 最終的にどの程度の理論や基底関数を選ぶかは、計算対象の性質に加え、使用できるメモリ使用量、ディスクの空き容量を考慮に入れなければなりません。 何も考えずにMP2などの計算をすると数十GBのディスク容量が必要になることがあります。 TSUBAME3では、バッチジョブによる計算用のスクラッチ領域としてSSDによる約1.9TBの高速テンポラリディレクトリが準備されています。 計算ノードの/scr/$JOB_ID. $Qnameに自動生成されますが、batch終了時に自動的に削除されるため、注意が必要です。
また、自分のジョブがどの程度ディスクを使うかという感覚が必要です。 あらかじめMaxDiskというキーワードで、使用するディスクの大きさを指定しておくと安全です。 デフォルトの単位は、Gaussianword(=8Bytes)なので、誤解を防ぐためにも単位を付けましょう。 単位にはMb、Mw、Gbなどが使えます。

(例) # MP2=(Full、Direct)/Aug-cc-pVQZ Scf=Direct MaxDisk=2Gb

4.3. 基底関数の選択

Gaussianに内蔵されている基底関数を下表に示します。 一般に良く用いられる基底は、3-21G、6-31G、6-311G、D95V、D95、LANL2B、LANL2DZなどです。 また、必要に応じ分極関数として、Li以降の原子にd軌道を(一般に*で表す)、さらに水素原子にもp軌道を(**で表す)加えることも良く行われます。 ただし、3-21G*は、Na以降の原子にのみ34軌道を加えた基底です(論文には、3-21G(*)と記されていることがあります)。
6-31G*あるいば6-31G**が、標準的基底関数として定評があり、最もよく使われています。 陰イオンの場合には、Li以降の原子にdiffuse関数として拡がったs軌道及びp軌道を(+で表す)、さらに水素原子には拡がった8軌道を(++で表す)加えることも良く行われます。 Gaussianでは使える元素の種類も増えるとともに、最近、電子相関用にDunningが新たに開発したcc-pVDZ、cc-pVTZ、cc-pVQZ、cc-pV5Z、cc-pV6Z(※2)という新しい基底関数が含まれています。 これらの基底の評判は良く、論文でもよく使われるようになってきています。 ただし、cc-pVTZ以上は大きな基底関数ですので、適用する分子でどの程度の数の基底となるか見積もっておくことが必要です。 どのような基底を選択するかという判断の元にもなりますから、基底関数の総数は計算の前に必ず調べておくことが必要です。
(※2)Dz、Tz、QzのD、T、Qはそれぞれ、Double、Triple、Quadrupleの略で、2倍、3倍、4倍を意味します。 5倍、6倍は一般的ではないので数字が使われていますが、そこまで拡張する時代が来るとは、誰も思っていなかったのでしょう。

Gaussian内蔵基底関数

基底関数 使用可能元素 内容
STO-3G H-Xe Popleの最小基底
3-21G H-Xe Popleのvalence double-zeta
6-21G H-Ci Popleのvalence double-zeta
4-31G 4-31G Popleのvalence double-zeta
6-31G H-Kr Popleのvalence double-zeta
6-311G H-Kr Popleのvalence tripe-zeta
D95V H-Ne Dunning-Huzinagaのvalence double-zeta(NaとMgを除く)
D95 H-Cl Dunning-Huzinagaのdouble-zeta
SHC FrとRaを除く D95V(H-Ne)+Effective-CorePotential(ECP)(Na以下の元素)
CEP-4G H-Rn 最小基底ECP
CEP-31G H-Rn valence double-zeta ECP
CEP-121G H-Rn valencetriple-zeta ECP
LANL2MB H-Ba、La-Bi STO-3G(H-Ne)+最小基底ECP(Na-Bi)
LANL2DZ H、Li-Ba、La-Bi
cc-pVDZ H-Ar 電子相関用valence double-zeta(Li、Be、Na、Mgを除く)
cc-pVDZ H-Ar 電子相関用valence double-zeta(Li、Be、Na、Mgを除く)
cc-pVTZ H-Ar 電子相関用valencetriple-zeta(Li、Be、Na、Mgを除く)
cc-pVQZ H-Ar 電子相関用valellce quadruple-zeta(Li、Be、Na、Mgを除く)
cc-pV5Z H-Ar 電子相関用valencequintllple-zeta(He、Li、Be、Na、Mgを除く)
cc-pV6Z H-N 電子相関用valence sextuple-zeta(He、Li、Beを除く)
SV、SVP、TZV H-Kr Ahlrichseの基底関数
Midix H,C,N,O,F,P,S,Cl TruhlarのMIDI!基底関数

参考のため、上表に代表的基底関数について、原子1個あたりに必要な基底関数の数を示しておきます。 これを元に、例としてベンゼンの基底関数を実際に数えてみましょう。

3-21Gでは 9 x 6(C) + 2 x 6(H) = 66個
6-31G*では 15 x 6(C) + 2 x 6(H) = 102個 6-31G**では 15 x 6(C) + 5 x 6(H) = 120個

原子1個あたりに必要な基底関数の数

STO-3G 3-21G 3-21G* 6-31G 6-31G* 6-31+G* 6-31G**
H 1 2 2 2 2 2 5
Li-Ne 5 9 9 9 15 19 15
Na-Ar 9 13 19 13 19 23 19

4.4. Z 行列

構造の記述には、対称性の導入などが比較的簡単で便利なため、古くからZ行列表示がよく使われています。 Z行列では、結合長、結合角及び二面角を用いて原子の位置を指定します。以下にいくつか例を示します。

H2O

O
H 1 r1
H 1 r1 2 a1
r1=0.947
a1=105.5

HOOH

O
O 1 r1
H 1 r2 2 a1
H 2 r2 1 a1 3 t1
r1=1.1393
r2=0.949
a1=102.2
t1=115.2

H2CO

O
C 1 r1
H 2 r2 1 a1
H 2 r2 1 a1 4 180.

r1=1.184
r2=1.092
a1=122.15

下記の2つは、ダミー原子を使った例です。 ダミー原子はXで表しますが仮想的なもので、分子軌道の計算自体に直接影響することはありません。 このように直線分子の場合や対称性を明確にする場合には、しばしばダミー原子を使う必要があります。

HCN

X
C 1 1.
N 2 r1 90.
H 2 r2 90. 3 180.

r1=1.133
r2=31.059

NH3

X
N 1 1.
H 2 r1 1 a1
H 2 r1 1 a1 3 120.
H 2 r1 1 a1 3 -120.

r1=1.0025
a1=111.68

この場合、Z行列中に直接書き込んだ1.、120.、-120.という数値はOpt=Z-Matrixで構造最適化する場合は変数として扱われません。 r1、a1のみが変数として扱われます。

X
N 1 1.
K 2 rl 1 a1
H 2 r1 1 a1 3 t1
H 2 r1 1 a1 3 -t1
r1=1.OO25
a1=111.68
t1=120.

4.5. その他

4.5.1. キーワード使用例

構造最適化

# hf/6-31g* opt=z-matrix test
# mp2=(full、direct)/6-31g* opt=zmatrix test maxdisk=2Gb

遷移状態の構造最適化

# hf/6-31g* opt=(ts、z-matrix、noeigentest) test
# hf/6-31g* opt=(ts、z-matrix、noeigentest、calcfc) test

振動解析

# hf/6-31g* freq pop=(reg、npa)

電子相関の一点計算

# mp4/6-31g** test
# cisd/6-31g**test*

分子軌道の数値を出力する場合、pop=regかpop=fullを加えます。 pop=regの場合は占有軌道5つと空軌道を5つ、pop=fullではすべての軌道を出力します。 デフォルトではHFの軌道なので、使用した手法の軌道や電子密度を知りたい場合は、densit=currentも同時に指定して置くことが必要です。 既存のチェックポイントファイルを利用する場合は、さらに、'geom=checkguess=read'の二つを加えます。
また、一つの入力ファイルに連続するジョブを含めることができます。 例えば、構造最適化に続いて振動解析を連続で行う場合、OptとFreqを同時に指定せずに、-Link1-で二つ目のジョブを繋ぐことができます。-Link1-を使えばいくらでもジョブを継続できますが、出力リストが大きくなりすぎることがあります。

4.5.2. 便利なGaussianユーティリティ

チェックポイントファイルから構造データを取り出すユーティリティnewzmatが用意されています。

h2o.chkからmopac形式の構造データ(h2o.inp)として取り出す方法

(例) newzmat -ichk -omopac h2o

cache形式の構造データ(h2o.cac)として取り出す方法

(例) newzmat -ichk -ocache h2o

protein data bank(PDB)形式(h2o.pdb)の構造データとして取り出す方法

(例) newzmat -ichk -opdb h2o

チェックポイントファイルをASCII形式に変換するユーティリティfomchkが用意されています。

(例)  formchk h2o

上記のように実行するとh2o.fchkというASCII形式のファイルが作られます。 このデータを使ってGaussViewなどで分子軌道を描くことが出来ます。
さらに詳しい使い方や他のユーティリティ、ワードに関してはGaussianのユーザーマニュアルを参照して下さい。