Gaussian Input File データベース

[スポンサーリンク]

一般的な話題

[スポンサーリンク]

お疲れ様です！Maitotoxinです！

みなさん、計算化学好きですか？特にGaussian、好きですか？僕は好きです！自分の見えない電子の世界を、なんか見えるようにしてくれる量子化学計算ってすごいですよね〜。実際、蛍光分子を扱う以上、量子化学計算による軌道計算は欠かせないものになっています。

でも同時に、めちゃめちゃとっつきにくいものだと思います。弊ラボでも、計算設備は整っていますが、あまり自分でやろうとする方は少ないです。慣習化された汎関数や基底関数での計算なんて、数分でinputファイルを作ってソフトウェアに投げるだけなのに….。難しいんでしょうね。

このように知ってしまえば簡単なものも、知りかたがまず分からないのが専門性だと思います。僕自身、最初はunixコマンドやサーバーの繋げ方すらわからない・量子化学計算で何がわかるのかも知らない状態からGaussianの使い方を独学で学んだので、大変地獄でした。以下のサイトが役に立ったのでシェアします。

Computational Chemistry (このページが一番最強です！！)
計算化学：汎関数って何？
有機化学者のための計算化学
HPCマニュアル
Gaussian公式ホームページ

またこの書籍も役に立つと思います。

新版すぐできる量子化学計算ビギナーズマニュアル (KS化学専門書)

武次徹也

¥3,520(as of 05/15 16:33)

Amazon Rakuten Ichiba Yahoo

Amazon product information

なみぃさんの記事ゼロから始める！量子化学計算～遷移状態を求める～を読んで、最初の頃の気持ちを思い出したのと、これから学ぶ人の助けになればいいな！と思いGaussian input fileの例をここにまとめたいと思います。

まず前提。Gaussianを用いた量子計算の流れ

Gaussian計算の流れを、ざっと確認します。語弊はありますが、量子化学を知らなくても量子化学計算は始められます！ (注: もちろん勉強したほうがいいですよ？)

Gaussianの使用を人任せにする方に多いのが、ソフトウェアの使用に数学的知識・計算が必要だと思って避けているシチュエーションなのだと周りを見て感じました。

もちろん中で回っている計算を理解した方がいいですが（この本おすすめです。）、とりあえず走らせてみるのがユーザーとしての学習にはオススメします。

Gaussianの使用は非常にシンプルで

分子の構造情報(.xyzとかデカルト座標、もしくはz-matrix)の作成 (ソフトウェア: Gauss View, Mercury, Avogadoro, Chem3Dなど色々)
Gaussian Input file (.gjf)の作成 (ソフトウェア：text editもしくはGauss Viewなど)
GaussianへのSubmit・計算 (ソフトウェア: Gaussian)
Gaussian Output file .chk/.fchk/.log の確認・解析 (ソフトウェア: text edit, Gauss View, iqmolなど)

の流れで行います。本記事で扱うのは、そして最も大事なのは2. Gaussian Input file の作成の部分です。Gaussian Input fileとは我々にGaussianにどのような計算を回して欲しいのか？の全てをソフトウェアで読み取り可能な形式でまとめたテキストファイルになります。すなわち、綺麗なInputさえかけてしまえば、僕らがシュレディンガー方程式を手で解く必要はなく（そもそも不可能）！！量子化学計算は終わったも同然です!! (脚色してます。この文章で私をジャッジしないでください。)

では、Inputはどのように書くのか？Gauss Viewを使えば自動生成もできますが、Text editで作る方が色々便利です。ルールはHPCマニュアルや Gaussian公式ホームページに載ってはいるのですが、見本がないとイメージが湧きにくいですよね？そこでこのデータベースを作成しました。ぜひコピペして使って行ってください (注、ブログ記入型の都合で空行などの体裁が少し崩れているものもあります)。基本的には分子のxyz情報と最初2行を変えて、最後に2行の空行を加えることでそのままGaussianに投げられると思います。(もし間違っているものがありましたらご指摘いただけますと幸いです。)

例1) b3lyp/6-31g(d,p)レベルでの基底状態計算 (S0構造最適化, 振動数計算)

%nprocshared=28
%mem=100GB
%chk=TNT.chk
# opt freq b3lyp/6-31g(d,p)

****

0 1
C -4.59255 0.08304 -0.05852
C -3.07318 0.17974 -0.06585
C -2.35822 1.43873 -0.06606
N -2.97016 2.74777 -0.05954
O -2.21637 3.82816 -0.06039
O -4.14959 2.90090 -0.05332
C -0.94837 1.40993 -0.07316
C -0.22055 0.22236 -0.07994
N 1.20996 0.24774 -0.08700
O 1.85923 1.39310 -0.08718
O 1.83665 -0.76914 -0.09283
C -0.91178 -0.98707 -0.07972
C -2.31693 -1.04992 -0.07289
N -2.91158 -2.36710 -0.07344
O -2.14448 -3.43785 -0.08011
O -4.09243 -2.53018 -0.06800
H -5.18863 0.99767 -0.05309
H -4.91946 -0.45909 0.85345
H -4.92857 -0.45422 -0.97006
H -0.39416 2.33820 -0.07345
H -0.33521 -1.90221 -0.08503

最初の例ですのでインプットファイルの解説,書き方を以下に説明します。

%nprocshared=28　(使用するプロセッサ数。あなたのパソコンに応じて決めてください。)
%mem=100GB　(使用するメモリ数。あなたのパソコンに応じて決めてください。)
%chk=TNT.chk　(作成されるチェックポイントファイル＝アウトプットファイルの名前。)
# opt freq b3lyp/6-31g(d,p)　(ルートセクション=計算条件。この場合、opt=構造最適化、freq＝振動数計算をb3lyp/6-31g(d,p)レベルで行うと言う意味)

（空行を入れることでルートセクションの終わりを表す。）

****　（コメント行、この分子の名前とか。僕は適当に入れている。）

（空行を入れることでコメントセクションの終わりを表す。）

0 1　（電荷とスピン多重度。スペースで区切る。この場合、電荷0, 多重度1の計算）
C -4.59255 0.08304 -0.05852　(以下分子中原子のxyz座標情報。Z-matrixを書く方法もある)
C -3.07318 0.17974 -0.06585
C -2.35822 1.43873 -0.06606
N -2.97016 2.74777 -0.05954
O -2.21637 3.82816 -0.06039
O -4.14959 2.90090 -0.05332
C -0.94837 1.40993 -0.07316
C -0.22055 0.22236 -0.07994
N 1.20996 0.24774 -0.08700
O 1.85923 1.39310 -0.08718
O 1.83665 -0.76914 -0.09283
C -0.91178 -0.98707 -0.07972
C -2.31693 -1.04992 -0.07289
N -2.91158 -2.36710 -0.07344
O -2.14448 -3.43785 -0.08011
O -4.09243 -2.53018 -0.06800
H -5.18863 0.99767 -0.05309
H -4.91946 -0.45909 0.85345
H -4.92857 -0.45422 -0.97006
H -0.39416 2.33820 -0.07345
H -0.33521 -1.90221 -0.08503

（空行を入れることで構造情報の終わりを表す)

（追加の情報が必要な時は、追加。例9とか参照）

（空行を入れることで追加情報の終わりを表す)

（空行を入れることでinput fileの終わりを表す)

と、まぁ、全部の計算がこんな感じです。

例2) M062X/6-31g(d,p)レベルでの励起状態計算（S1構造最適化, 振動数計算）

%mem=100GB
%nprocshared=28
%chk=TNT.chk
# opt freq td=(singlets,nstates=10,root=1) b3lyp/6-31g(d,p)

****

0 1
C -4.59255 0.08304 -0.05852
C -3.07318 0.17974 -0.06585
C -2.35822 1.43873 -0.06606
N -2.97016 2.74777 -0.05954
O -2.21637 3.82816 -0.06039
O -4.14959 2.90090 -0.05332
C -0.94837 1.40993 -0.07316
C -0.22055 0.22236 -0.07994
N 1.20996 0.24774 -0.08700
O 1.85923 1.39310 -0.08718
O 1.83665 -0.76914 -0.09283
C -0.91178 -0.98707 -0.07972
C -2.31693 -1.04992 -0.07289
N -2.91158 -2.36710 -0.07344
O -2.14448 -3.43785 -0.08011
O -4.09243 -2.53018 -0.06800
H -5.18863 0.99767 -0.05309
H -4.91946 -0.45909 0.85345
H -4.92857 -0.45422 -0.97006
H -0.39416 2.33820 -0.07345
H -0.33521 -1.90221 -0.08503

例3) cam-B3LYP/6-31g(d,p)レベルでの励起状態計算（S0構造）(S0構造が入ったTNT.chkファイルが入ったディレクトリ内で行う)

%mem=100GB
%nprocshared=28
%chk=TNT.chk
# sp td=(singlets,nstates=10,root=1) cam-b3lyp/6-31g(d,p) Guess=Read Geom=AllCheckpoint

****

ここの “Guess=Read Geom=AllCheckpoint”は非常に使いやすく、前回行った計算結果・構造からそのまま別の計算を走らせるときに有用です。例えば励起状態計算後、NTO計算　を行う場合（遷移軌道の可視化）は例4のようになります。

例4) NTO計算 (S1遷移軌道の可視化)

%mem=100GB
%nprocshared=28
%chk=TNT.chk
# M062X/6-31g(d,p) Geom=AllCheck Guess=(Read,Only) Density=(Check,Transition=1) Pop=(Minimal,NTO,SaveNTO)

0 1

これで計算できますが、もともとあったchk fileを上書きしたくない場合は

%mem=100GB
%nprocshared=28

%oldchk=TNT.chk
%chk=TNT_NTO.chk
# M062X/6-31g(d,p) Geom=AllCheck Guess=(Read,Only) Density=(Check,Transition=1) Pop=(Minimal,NTO,SaveNTO)

0 1

にすることで、chk fileを古いのを上書きすることなく新たに生み出せるそうです。（ケムステスタッフからの寄稿）

例5) M062X/6-31g(d,p)レベルでの励起状態計算（T1構造最適化, 振動数計算）

%mem=100GB
%nprocshared=28
%chk=TNT.chk
# opt freq td=(50-50,nstates=10,root=1) b3lyp/6-31g(d,p)

****

0 1
C -4.59255 0.08304 -0.05852
C -3.07318 0.17974 -0.06585
C -2.35822 1.43873 -0.06606
N -2.97016 2.74777 -0.05954
O -2.21637 3.82816 -0.06039
O -4.14959 2.90090 -0.05332
C -0.94837 1.40993 -0.07316
C -0.22055 0.22236 -0.07994
N 1.20996 0.24774 -0.08700
O 1.85923 1.39310 -0.08718
O 1.83665 -0.76914 -0.09283
C -0.91178 -0.98707 -0.07972
C -2.31693 -1.04992 -0.07289
N -2.91158 -2.36710 -0.07344
O -2.14448 -3.43785 -0.08011
O -4.09243 -2.53018 -0.06800
H -5.18863 0.99767 -0.05309
H -4.91946 -0.45909 0.85345
H -4.92857 -0.45422 -0.97006
H -0.39416 2.33820 -0.07345
H -0.33521 -1.90221 -0.08503

これ、T1指定した構造最適化じゃなくて、一番エネルギーの低い励起状態に落ち着きます。まぁ、普通はT1に落ち着くと思います。

例5-2) tdを使わないT1構造最適化

%mem=100GB
%nprocshared=28
%chk=TNT.chk
# opt freq b3lyp/6-31g(d,p)

****

0 3
C -4.59255 0.08304 -0.05852
C -3.07318 0.17974 -0.06585
C -2.35822 1.43873 -0.06606
N -2.97016 2.74777 -0.05954
O -2.21637 3.82816 -0.06039
O -4.14959 2.90090 -0.05332
C -0.94837 1.40993 -0.07316
C -0.22055 0.22236 -0.07994
N 1.20996 0.24774 -0.08700
O 1.85923 1.39310 -0.08718
O 1.83665 -0.76914 -0.09283
C -0.91178 -0.98707 -0.07972
C -2.31693 -1.04992 -0.07289
N -2.91158 -2.36710 -0.07344
O -2.14448 -3.43785 -0.08011
O -4.09243 -2.53018 -0.06800
H -5.18863 0.99767 -0.05309
H -4.91946 -0.45909 0.85345
H -4.92857 -0.45422 -0.97006
H -0.39416 2.33820 -0.07345
H -0.33521 -1.90221 -0.08503

S0の構造最適化、振動数計算から太字で示したところが変わっています。

例6) 水素原子のみの構造最適化 (XRD結晶構造解析後の構造最適化に使用)

%nprocshared=28
%mem=100GB
%chk=TNT.chk
# opt freq b3lyp/6-31g(d,p)

****

0 1
C -1 -4.59255 0.08304 -0.05852
C -1 -3.07318 0.17974 -0.06585
C -1 -2.35822 1.43873 -0.06606
N -1 -2.97016 2.74777 -0.05954
O -1 -2.21637 3.82816 -0.06039
O -1 -4.14959 2.90090 -0.05332
C -1 -0.94837 1.40993 -0.07316
C -1 -0.22055 0.22236 -0.07994
N -1 1.20996 0.24774 -0.08700
O -1 1.85923 1.39310 -0.08718
O -1 1.83665 -0.76914 -0.09283
C -1 -0.91178 -0.98707 -0.07972
C -1 -2.31693 -1.04992 -0.07289
N -1 -2.91158 -2.36710 -0.07344
O -1 -2.14448 -3.43785 -0.08011
O -1 -4.09243 -2.53018 -0.06800
H 0 -5.18863 0.99767 -0.05309
H 0 -4.91946 -0.45909 0.85345
H 0 -4.92857 -0.45422 -0.97006
H 0 -0.39416 2.33820 -0.07345
H 0 -0.33521 -1.90221 -0.08503

原子の横の-1がfreeze(固定)、0が移動可能のコマンドになります。すなわち水素だけ動いていいよ！って状態での構造最適化になります。

おまけですが分子が大きいと書くのが大変なので、僕は以下のPythonコードでgjf fileより準備しています。

filename_input = “./A.gjf”
filename_output = “./output.gjf”

start_line_num = 9
output_text = “”

with open(filename_input, “r”) as f:
lines = f.readlines()

for line_num, line in enumerate(lines):
line = line.strip() # remove break line

if line_num+1 < start_line_num: # line_num starts from 0
output_text += “{}\n”.format(line) # add break line at the end
else:
line_array = [e for e in line.split(” “) if e != “”] # split input line by space and make an array
if len(line_array) > 1:
if line_array[0] == “H”:
line_array.insert(1, str(0)) # add “0” in case of “H”
else:
line_array.insert(1, str(-1)) # otherwise, add “1”

line_output = “\t”.join(line_array)
output_text += “{}\n”.format(line_output) # add break line at the end

# print(“=== output text ====”)
# print(output_text)

with open(filename_output, “w”) as f_out:
f_out.write(output_text)

f = open(‘output.gjf’, ‘a’, encoding=’utf-8′, newline=’\n’)
f.write(‘\n’)
f.write(‘\n’)
f.close()

例8) ブタンの回転、ポテンシャル面の計算 (ケムステスタッフからの寄稿)

%chk=butane_scan2
%mem=100GB
%nproc=36
#p opt=modredundant b3lyp/6-31g(d) pop=full

memo

0 1
C -1.81226530 -0.33960474 -0.18951040
H -2.75438476 -0.79975943 -0.40297839
H -1.10701678 -1.08952983 0.10227189
C -1.29892308 0.38635154 -1.44691537
H -2.00326075 0.26246181 -2.24281724
H -0.35589216 -0.02730756 -1.73756778
C -1.12912084 1.88597588 -1.14049402
H -0.42506352 2.00984515 -0.34434097
H -0.77209312 2.39031163 -2.01403661
H -2.07222653 2.29971720 -0.85020142
C -1.98469104 0.67604887 0.95518044
H -2.68960723 1.42622055 0.66322903
H -2.34178564 0.17172469 1.82870237
H -1.04248246 1.13588249 1.16894671

D 11 1 4 7 S 36 10.000000

例９) 原子ごとに基底関数を指定する場合

%nprocshared=8
%mem=8GB
%chk=tellurophene.chk
# opt freq M062X/genecp

****

0 1
C 1.30655 0.82344 0.00492
C 2.77391 0.99154 0.00168
N 3.31017 2.23958 -0.00402
C 4.65665 2.40580 -0.00701
N 5.47368 1.32308 -0.00430
C 4.94379 0.07449 0.00139
N 3.59724 -0.09075 0.00438
C 0.69855 -0.36588 0.01053
C -0.68897 -0.37931 0.01309
C -1.31990 0.79801 0.00976
C -2.79024 0.93767 0.01194
N -3.59244 -0.16036 0.01764
C -4.94194 -0.02122 0.01962
N -5.49593 1.21687 0.01592
C -4.70003 2.31522 0.01024
N -3.35058 2.17509 0.00826
Te -0.02269 2.46360 0.00106
H 5.07835 3.40146 -0.01154
H 5.59443 -0.78913 0.00356
H 1.25457 -1.29508 0.01308
H -1.22688 -1.31910 0.01765
H -5.57573 -0.89727 0.02416
H -5.14094 3.30252 0.00729

Te 0
lanl2dz
****
H C N 0
6-31G**
****

Te 0
lanl2dz

genもしくはgenecpがキーワードです。重原子を含む計算で使う方が多いのではないでしょうか？最後の段落で、原子ごとの基底関数、そして有効内殻ポテンシャル（ECP）を指定しています。（何言っているか分からない場合はこちらのページ参照）。ECPを指定しない場合でも原子と基底関数によっては問題ありませんが、、、どういうシチュエーションですかね？あまり思い浮かばないです。

こちらのブログが大変参考になりました。