DV-Xα分子軌道法による第一原理計算

密度汎関数法の一つであるDV-Xα分子軌道法を用いて第一原理計算を行う方法を説明します。

準備すべき事項は以下のようなものですが、 分子, 結晶 (実際には巨大な分子として扱う) の電子状態が視覚的に分かりますので, 数値計算を行いながら理論体系を学ぶことをお勧めします。

[内容]
  1. DV-Xα分子軌道法の概要
  2. DV-Xα分子軌道計算の流れ
  3. DV-Xα計算プログラムの使用法

1. DV-Xα分子軌道法の概要

DV-Xα分子軌道法は G.S.Painter, D.E.Ellis (1970), H.Adachi, M.Tsukada, C.Satoko (1978) らによって 提案・開発された第一原理計算の手法です。

分子軌道法とは?

原子における電子の波動関数は、Hでは1s1, Siでは1s22s22p63s23p2のように 複数の水素原子モデルの波動関数で表されます。 この1s, 2pなどを原子軌道といいます。 (実際には電子同士の相互作用があるため波動関数は少し変化しますが...)
それでは分子の場合はどうでしょうか。 例えば水素分子H2(2つの水素原子をA, Bとします)では
ψ(r) = CAχA(r) +CBχB(r)
のように原点をA, Bの水素原子の位置にずらした2つの1s波動関数χA(r), χB(r)の線形結合として表せると考えられます。 係数CA, CBは、 変分原理に基づいてシュレディンガ─方程式を解いて得ます。 もちろん、電子同士のクーロン斥力が働きますからもう少し丁寧に扱う必要がありますが、 その場合は2s, 2p, 3s, ....など (場合によっては平面波なども)もっとたくさんの波動関数の線形結合を作ってやります。

第一原理計算とは?

実験などで得たパラメータ(有効質量、吸収スペクトルなど)を用いずに、 基本原理に基づき物理量の計算を行うものです。 入力するのは原子番号と原子の座標のみです。

バンド計算とどこが違うの?

結晶中の原子の数は1モルで6x1023個もあります。 常識的な時間内に計算を行うことを考えると、 扱うことのできる原子の数に限りがあります。 そこである種の近似をするわけですが、 その一つが「原子配列の周期性を利用したバンド理論」、 もう一つが「結晶に含まれる原子の一部を、計算可能なサイズの巨大分子として切り出して 計算を行うクラスター法」です。 DV-Xα分子軌道法で結晶を扱うにはこのクラスター法を用います。

密度汎関数法とは?

水素原子のような一電子のシュレディンガー方程式は厳密に解くことができますが、 He以上の原子番号の原子や分子のような 多電子系のシュレディンガー方程式は厳密には解けません。 これでは先に進めないので近似を導入してやる訳です。 ある一つの電子に着目して、 その電子以外のすべての電子から受ける「古典的クーロンポテンシャル」、 パウリの排他則に基づく電子同士の相互作用である「交換エネルギー」、 電子同士が影響し合って生じる運動エネルギーの変化に相当する「相関エネルギー」 をすべてポテンシャルの項に押し込めてしまう(等価ポテンシャルで置き換える)、 という方法を採ります。
ここで「おや?」と思った人は鋭い。 そうです、波動関数を決定しないと電子の位置が分からないではないですか! それなのにどうやって他の電子からのクーロンポテンシャルや交換エネルギー、 相関エネルギーを計算するのでしょうか?
密度汎関数法では、 電子の固有エネルギーの和、 電子間のクーロンポテンシャル、交換・相関エネルギー を電子密度(位置の関数)の汎関数として表してやります。 つぎに、これらのエネルギーの総和が停留(極小)になるように、 変分原理に基づき密度関数を決定します。

具体的には以下のようにします。
まず密度関数ρin(r)を仮定します。 一電子シュレディンガー方程式に繰り込むべきポテンシャルの項 電子同士のクーロンポテンシャル φ(r)=∫[ρin(r')/|r-r'|]dr' と交換・相関ポテンシャル vxcin] を求めます。
次に、veff = v + φ + vxcをポテンシャル項として シュレディンガー方程式を解いて、波動関数の2乗 Σnii|2outを求めます。
ρout≠ρin であれば、 ρin を ρout で置き換えてはじめに戻って計算を繰り返します。
このようにして ρoutin になるまで繰り返して (セルフコンシステントな手法といいます)、 電子密度、波動関数を得ます。

Xα法とは?

局所密度近似法(Local Density Approximation; LDA)の一つです。 これは、 交換・相関エネルギーを局所的な密度関数の汎関数 Exc=-A∫ρ(r)4/3dr として表す方法です。

DV法とは?

離散変分(Discrete Variational)数値積分法のことです。 ∫χ*(r)v(r)χ(r)drの積分を計算するのに、 原子核付近あるいは電荷密度の大きい領域を中心に疑似乱数積分で求める方法です。 ψ(r)=Σckχk(r) において原子軌道関数を基底χk(r)にとります。
基底関数として他には平面波型関数、スレーター型、ガウス関数型などがあります。

2. DV-Xα分子軌道計算の流れ


3. DV-Xα計算プログラムの使用法

DV-Xα法による計算の大まかな流れは次の通りです。 計算の大部分は、研究室のUnixホストにログインして実行します。
(ログインの方法Unixの使い方研究室のUnixホストの情報 などをチェックしてください。)
MS-DOSやWindows上でも実行できますので、 自分のパソコンで計算してみたい人は阿部に聞いて下さい。

[入力ファイルの作成]

計算の過程でたくさんのファイルができますから、 新しいディレクトリを作ってそこでプログラムを実行した方がいいでしょう。

● 原子番号、原子座標などを記述したファイルを作成して、 DV-Xα分子軌道計算の入力ファイル作成の準備をする。
"f01" というファイルを作成して "makef05" を実行することにより、 DV-Xα分子軌道計算プログラムの入力ファイル "f05" を作成します。
まず、"f01" を作成してください。 NaCl の場合の "f01" の例を以下に示します。

| Z ||NEQ||   X    ||   Y    ||   Z    |
   11    1   0.00000   0.00000   0.00000
   17    2   2.82000   0.00000   0.00000
   17    2  -2.82000   0.00000   0.00000
   17    2   0.00000   2.82000   0.00000
   17    2   0.00000  -2.82000   0.00000
   17    2   0.00000   0.00000   2.82000
   17    2   0.00000   0.00000  -2.82000
---------------------------------------------
|NEQ||  CHG   ||U/D||   RD   ||   VD   |    1 
    1   1.00000
    2  -1.00000
---------------------------------------------
    0     Unit     (0:angstrom  1:atomic)
    0     Spin     (0:non-spin  1:spin  )
    1     M.P.     (0:No        1:Yes   )
    0     Sample Point ( <100000, =0 autoset )
Z
原子番号。
NEQ
等価な原子につける番号。上の例ではNaの周囲の6個のCl原子は等価ですから全て "2" としています。
X, Y, Z
原子の座標。単位はÅまたはa.u.(atomic unit, 1a.u.は水素原子のボーア半径)。
CHG
イオン結晶の場合は各原子(NEQ番号)の電荷を設定します。
Unit
座標の単位の設定。
Spin
スピンを考慮するかどうの設定。
M.P.
マーデルングポテンシャルを取り入れるかどうかの設定。上の例では "1" を選んでいますから、"f03" という設定ファイルを用意する必要があります。(後述)
Sample Point
DV数値積分のサンプル点数。通常は "autoset" で構いません。
それぞれの数値は "|" で挟まれた範囲に書く必要があります。 小数点以下の桁数は上の例より少なくて構いませんが、必ず小数点は書いてください。 (整数はダメです)
その他の例が /usr/local/dvxa/ex/ 以下にありますからコピーして使用してください。

● DV-Xα分子軌道計算の入力ファイルを作成する。
"makef05" を実行して "f05" というDV-Xα分子軌道計算の入力ファイルを作成します。

makef05
同時に何も書かれていない "f25" というファイルも作成されます。 このファイルは対称軌道を考慮して効率よく計算するために使われますが、 何も書かれていなければ対称軌道計算は行われません。 対称軌道計算を行うには後述の "symOrb" というプログラムを 実行する必要があります。
ファイルの中身については、 参考文献[1]を参照して下さい。

○ マーデルングポテンシャルのデータファイルを作成する(必要な場合のみ)
"f03" というファイルを作成します。
上述のNaClの場合、次のようにします。

    8    2    0
   1.00000  -1.00000
       10.6580000000       10.6580000000       10.6580000000
        0.0000000000        0.0000000000        0.0000000000    1
        0.5000000000        0.0000000000        0.0000000000    2
        0.0000000000        0.5000000000        0.0000000000    2
        0.0000000000        0.0000000000        0.5000000000    2
        0.0000000000        0.5000000000        0.5000000000    1
        0.5000000000        0.0000000000        0.5000000000    1
        0.5000000000        0.5000000000        0.0000000000    1
        0.5000000000        0.5000000000        0.5000000000    2
    2    2    2
             0.00000             0.00000             0.00000
             2.00000             2.00000             2.00000
1行目    8     2     0
第1カラム:Unit Cell内の原子数、第2カラム:原子の種類、 第3カラム:欠陥として取り除く原子の数。
2行目    1.00000   -1.000000
順に1番目の原子(Na+)の価数、2番目の原子(Cl-)の価数。 (8F10.5:10文字毎、小数点以下5桁、8個まで)
3行目
Unit Cellのx, y, z方向の大きさ(格子定数)。単位はa.u.(atomic unit)。 1 a.u. = 0.529177Å。
4〜11行目
第1〜3カラム:原子のx, y, z座標。単位は3行目で指定した長さ。(3F20.10, I5:20文字毎、小数点以下10桁の3個浮動小数点と5文字の整数。) 第4カラム:原子の種類。
12行目
x, y, z方向へのUnit Cellの並進拡張。(3I5: 5文字毎の整数3個)
13行目
どうしてもマーデルングポテンシャル用の格子と、実際のクラスターの原点をずらしたいとき、マーデルングポテンシャル用の格子座標の原点を入力。(3F20.10)
14行目
並進拡張した点電荷をどこまで含めるかを設定。 並進拡張したときに電荷分布が非対称になる場合に設定する。 上の例では "2 2 2" を設定しているので、±2, ±2, ±2以下の座標にのみ 点電荷を置くようになる。
詳しくは参考文献[1]を参照して下さい。

○ 対称軌道データを作成する(必要な場合のみ)
"symOrb" という Mathematica プログラムを用いて作成します。 現時点ではMathematicaがないため準備中。 詳しくは参考文献[1]を参照して下さい。

[数値計算の実行]

次に計算を始めます。 必要なファイルは "f05" と "f25" です。 マーデルングポテンシャルを考慮する場合は "f03" も必要です。

次のコマンドを実行すると計算が開始されます。 実行オプション "n" は計算サイズを表しています(後述)。 "." はカレントディレクトリ(現在のディレクトリ)を表しています。

dvscat n .
10原子ぐらいであれば数分程度で終了します。

収束状況は "f36" というファイルに記録されます。 第4カラム目の "SHOKICHI" と第5カラム目に各軌道の ポピュレーション(占有率)の初期値と終了値が書かれていますが、 この差が0.00001以下になると、 セルフコンシステントであると判断され計算が終了します。
"f06z" というファイルには今までの繰り返し計算の記録が書かれています。 1行目が計算が実行された回数を表していて、 収束した場合には "converged" と記録されます。 "dvscat" を実行すると10回の繰り返し計算が行われますが、 "converged" と記録されていない場合はさらに "dvscat" を繰り返して実行します。

100原子ぐらいのクラスター計算

まず、作成したクラスターモデルがどの程度の変数サイズを必要とするか確かめます。 "f05", "f25" があるディレクトリで
dimchk n
を実行すると、
1. NAT 2. NDAT 3.NCORE 4. NAOS
7 2 12 98

5. ISO 6. MM 7. KSIZ 8. NSIZ
98 98 4851 9604
と表示されます。これらの1〜8は次の表に示した変数サイズに対応しています。 全ての値が表の変数サイズを越えないように n, w, g, s, u, h を選択してください。

通常の計算の場合nサイズwサイズgサイズ
スピンを考慮に入れた場合sサイズuサイズhサイズ
1. 原子数500999999
2. 種類分けの数159999
3. 種類ごとの原子軌道数の和70999999
4. 全体の原子軌道数78012002400
5. ある対称ブロックでの対称軌道数の最大値3005001500
6. 対称軌道を構成する基底の延べ数120003000030000
7. 重なり積分などを1次元化したときの要素数の和53000150000700000
8. ある対称ブロックでの係数の数の最大値900002500002250000
9. マーデルングポテンシャルの数500009999999999

例えば "g" を選択したら、次のように実行します。 gサイズの場合、計算に時間がかかりますから(数時間程度)、 最後に "&" をつけてバックグラウンドで実行します。

dvscat g . &

[計算結果の表示]

状態密度

まず、"maked04" を実行して "d04" というファイルを作成します。
maked04
"d04" は以下のようになっています。目的に合わせて編集してください。
GAUS
 -15.00000  10.00000   0.20000    0
   0.5       0.5       0.5       0.5       0.5       0.5       0.5       0.5    
   0.5       0.5       0.5       0.5    
    6    0    0    0    1    2    3    0    0    0    4    5    6
   0.00000   0.00000   0.00000   1.00000   1.00000   1.00000   0.00000   0.00000
   0.00000   1.00000   1.00000   1.00000


    2  1Na  1Na  1Na  1Na  1Na  1Na  2Cl  2Cl  2Cl  2Cl  2Cl  2Cl

  1 Na 1s   1 Na 2s   1 Na 2p   1 Na 3s   1 Na 3p   1 Na 3d   2 Cl 1s   2 Cl 2s 
  2 Cl 2p   2 Cl 3s   2 Cl 3p   2 Cl 3d 
1行目:分布関数の設定 (A4:4文字)
エネルギー準位に幅をもたせるのに使う分布関数の設定。 ガウス関数("GAUS")かローレンツ関数("LORE")を選択します。
2行目:エネルギーの設定 (3F10.5, I5:10文字分, 小数点以下5桁の実数3個と5文字分の整数)
第1,2カラム:表示するエネルギー範囲の下限と上限(eV)、 第3カラム:計算を行う刻み幅(eV)、 第4カラム:"0" だと計算値でプロット、"1" だとフェルミ準位を 0eV としてプロット。
3〜4行目:各原子軌道の半値幅 (8F10.5:10文字分, 小数点以下5桁の実数8個)
分布関数の半値半幅(半値全幅の半分)(eV)を設定。 各カラムは12〜13行目(最後の2行)にの原子軌道に対応。
5行目:出力する部分状態密度の設定 (16I5:5文字分の整数)
第1カラム:出力する部分状態密度の総数、 第2カラム以降:各軌道の部分状態密度を出力ファイルの何カラム目に対応させるかを設定する。 (12〜13行目(最後の2行)にの原子軌道に順番に対応)。 "0" だと出力されない、"n" (n≧1)だと第(n+1)カラムに部分状態密度が出力される。 出力ファイルの第1カラムはエネルギーである。
6〜7行目:各軌道の部分状態密度の重み係数
各軌道の部分状態密度の重み係数を設定する。 順番に12〜13行目(最後の2行)にの原子軌道に順番に対応。 等価な原子の数 (NEQ) の逆数にすると原子当たりの状態密度が計算できる。
原子軌道を表すラベル
原子の種類、元素名、軌道が占めされている。
"d04" の作成が終わったら、次のように "dos" を実行します。 実行オプションである "計算サイズ" には "dvscat" に用いたサイズと同一のものを選びます。
dos n
計算が終了すると "dd7" というファイルができている筈ですので、 使い慣れたグラフソフトなどでプロットしてください。 "dd7" の第1カラムはエネルギー, 第2カラム以降は "d04" で指定した部分状態密度です。
"d07" というファイルはWindows用の "DVPLOT" というプログラムで表示するためのファイルです。

状態密度のプロット例 (NaCl)

波動関数、電子密度

ある平面上の波動関数、電子密度の等高線表示をしてみましょう。
まず、"wavnum" を実行して波動関数の番号を調べます。
実行オプションは "dvscat" の計算サイズと同一にします。
wavnum n
   Input Ef Mode.
1 ←0か1を入力します。
    0:計算時の固有エネルギーで表示、1:フェルミ準位を0eVとして表示

"wavnum" というファイルが出来ている筈です。 例としてフェルミ準位とその一つ上の ファイルの中から第2カラムが "0.0000 EV" となっている行を探します。 次のように実行すると見つけられます。
grep -2 "0.0000 EV" wavnum
第3カラムを見るとフェルミ準位は59番目であることが分かります。 その1つ上の準位は60番目です。
 E( 57)        -0.0362 EV        57
 E( 58)        -0.0177 EV        58
 E( 59)         0.0000 EV        59
 E( 60)         8.7985 EV        60
 E( 61)        10.7062 EV        61
次に、"makec04" を実行して質問に答えていきます。
makec04
  >Select mode
   0: input 3 atoms
   1: input 2 atoms & X axis vector
0 ←3原子の座標で平面を決定するか、
    2原子と横軸のベクトルで平面を決定するかを選択します

>input 1st atom
1 ←1つ目の原子をファイル "f01" の番号にしたがって入力します
>input 2nd atom
2 ←2つ目の原子の番号を入力。1つ目の原子の右真横に設定される。
>input 3rd atom
4 ←3つ目の原子の番号を入力。
    1つ目と2つ目の原子の上側に位置するように設定される。

>input x length
20 ←縦軸の範囲を入力。単位は a.u.
    0からこの値までをプロットします。

>input y length
20 ←横軸の範囲を入力。
>input 1st atom position (x,y)
10,10 ←1つ目の原子の座標を指定。
>input Mesh number
101 ←メッシュ数を入力。(縦横の分割数)+1
>reverse x,y ? 0: No 1: Yes
0 ←縦軸と横軸を入れ替えるかどうか答えます。
    通常は0と答えてください。

次のような "c04" というファイルが出来ている筈ですので、 6行目に表示する波動関数の個数と番号を順に設定してください。 フォーマットは 16I5 (5文字分の整数、16個まで) です。
   -3    2    3    0
  101  101
   -0.1000000000E+02   -0.1000000000E+02    0.0000000000E+00
    0.2000000000E+00    0.0000000000E+00    0.0000000000E+00
    0.0000000000E+00    0.2000000000E+00    0.0000000000E+00
    2   59   60

input atoms       :     1    2    4
x,y length        :    20.0000   20.0000
1st atom position :    10.0000   10.0000
reverse x,y       :    No
次に "contr" と "cmtxt" を実行して、 波動関数の出力ファイル "w??.txt" (??には波動関数の番号が入る) と 電荷密度の出力ファイル "chg.txt" を得ます。 "contr" の実行オプションは "dvscat" の計算サイズと同一にします。
contr n
cmtxt
  >Input number of columns (0:same as C04)
0 ←通常は0と答える
"chg.txt", "w59.txt" "w60.txt" というファイルができている筈です。 "wavplot" を実行すると波動関数や電荷密度の濃淡図が得られます。 出力ファイルはPostscript "???.ps" と PNG画像 "???.png" です。 3つ目の例のように第2,第3オプションに最小値と最大値を指定することもできます。 デフォルトでは -0.1 と 0.1 です。
wavplot w59.txt
wavplot w60.txt
wavplot chg.txt 0 0.2

HOMO (Highest Occupied Molecular Orbital) (フェルミ準位)
波動関数、電子密度のプロット例:NaCl
LUMO (Lowest Unoccupied Molecular Orbital) (フェルミ準位の1つ上の準位)
電子密度


参考文献

  1. 「はじめての電子状態計算」足立裕彦   監修, 小和田善之, 田中功, 中松博英, 水野正隆   著, 三共出版 (1998)
  2. 「量子材料科学入門 ─DV-Xα法からのアプローチ─」足立裕彦   著, 三共出版 (1991)
  3. DV−Xαホームページ
  4. 「密度汎関数法とその応用」菅野暁   監修, 里子允敏, 大西楢平   著, 講談社 (1994)