地域密着型エリア広告配信リクルートの地域密着型広告ドコイク?アドネットワーク Ads by ドコイク?[無料でホームページを作成] [通報・削除依頼]
[無料でホームページを作成] [通報・削除依頼]

ABINIT(The Lesson 1-4)

ABINIT: the tutorials (for v.7.6 tutorial) 作成 2014/07/05 
(全てを記述することは無理なので私の経験と要約した意訳のみ

----------------------------------------------------------------------------------------------------------------------------

The lesson 1 はH2分子を取り扱う: 全エネルギー、電子のエネルギー、電子密度、結合距離、エネルギーを得る。

1.1. 全エネルギーといくつかの関連する値の計算

 ターミナルやコンソールを開いて、下記のコマンドを入力してください。

cd $HOME/abinit-7.8.1/tests/tutorial/Input

 abinit/tests/tutorial/Inputのディレクトリがターミナルやコンソールで様々な操作可能なディレクトリ(カレントディレクトリ = cd)になります。このディレクトリで、チュートリアルで用いる入力ファイルを得ることができます。ワーキングディレクトリが必要になるので、このディレクトリのサブディレクトリ( ~abinit/tests/tutorial/Input/Work)を下記のコマンドで作ります。

mkdir Work

 下記のコマンドでWorkをカレントディレクトリにします。

cd Work

 チュートリアルのほとんどの操作をこのワーキングディレクトリで行います。Workにtbase1_xを下記のコマンドでコピーしてください。最後のドットは忘れないようにしてください。これはカレントディレクトリを意味します。

cp ../tbase1_x.files .

 tbase1_x.filesは、他の名称のファイルを作成するといったAbinitコードが要求する情報を与えます。このファイルの1行目と2行目を次のように書き換えます。

tbase1_1.in
tbase1_1.out

 このような操作は他のチュートリアルでも行います。最後の方には計算で用いる元素の擬ポテンシャルの場所(アドレス)を記述します。tbase1_x.filesを閉じて、下記のコマンドでカレントディレクトリ(Work)にコピーします。

cp ../tbase1_1.in .

 この操作で、Workにコピーされたファイルの名称は tbase1_1.in になっています。それで、このWorkディレクトリで下記のコマンドを入力します。

abinit < tbase1_x.files >& log

 少し待つと、それが実行されます。下記のコマンドでWorkにあるコンテンツをみることができます。

ls

 次のようないくつかのファイルがあることが分かります。

abinit log tbase1_1.in tbase1_1.out tbase1_x.files

tbase1_xo_DDB tbase1_xo_DEN tbase1_xo_EIG tbase1_xo_WFK

 log や tbase1_1.out などの出力ファイルが作られているのが分かります。それらが正しいかチェックします。 ~abinit/tests/tutorial/Refs にあるファイルとtbase1_1.outの差分を下記のコマンドで作りましょう。

diff tbase1_1.out ../../Refs/tbase1_1.out | more

 おそらく、空白を無視することが必要になるでしょう。その場合には、diff を diff -b に換えます。コンパイラなどの違いで、いくつかの差分が得られます。多くとも1.0d-10オーダーの数値的な違いを得るでしょう。差分として表示されたものには指定したファイルのアドレス(パス)の違いもあるでしょう。この他に計算時間の違い。
 では次に tbase1_1.in を見てみましょう。約50行程で、ほとんどがコメントなので、本当に入力に必要な記述は多くありません。直ぐに全てを理解しようとしてはダメです。チュートリアルを経験して、中身の記述を理解してください。入力ファイルの記述に必要な variablesの詳細は、アルファベット順に並べられて index of all variables に記載されています。基本的な計算に必要なものからGWなどの計算のために必要なvariablesまで種類分けして記述されています。
 そこには余りにも多くの記述があって気が滅入ってしまうでしょうが、実際には、いくつかのvariablesにディフォルトの値が設定してあるため、入力ファイルに僅かなvariablesを記述するだけで計算が可能です。
 ディフォルトが存在しないか、ディフォルトを避けるべきvariablesをtbase1_1.inで見つけることができます。結晶の構造に関するものが殆どですので当たり前ですね。そこにあるvariablesは非常に重要です。少し時間を使って、そのvariablesを見てみましょう。

  • acell, 格子定数を指定します。
  • ntypat, 原子の種類の数を記述します。
  • znucl, 原子番号を記述します。ntypatで記述した分だけ記述します。
  • natom, 単位胞にある原子の数を記述します。
  • typat, 下記のxcartで座標を指定する原子の種類を記述します。原子番号ではなくて、znuclで記述した順番(左から1, 2, 3となる)で原子の種類を指定します。natomで記述した行だけ記述します。
  • xcart, 原子位置をデカルト座標で入力します。natomで記述した行だけ記述します。
  • ecut, 計算の精度に関連します。この値が大きいほど多くの平面波を用います。
  • nkpt, k点数を記述します。
  • nstep, 繰り返し計算をさせる最大の回数です。ディフォルトでは繰り返しの計算が30回で終了します。
  • toldfe, 繰り返し計算をしますが、前回と今回の全エネルギー差がこの範囲内となっていれば計算を終了させます。
  • diemac. これは巨視的な誘電率を diemacで指定することで、繰り返し計算のスピードアップをさせるものです。物質の種類によって値が異なります。金属ではディフォルト値の106、半導体では12.0、絶縁体では2.0-4.0(20万円もするCrystal Studio v.13ではこの中間の3.0 になるようにしています)、分子では1.5-3.0程度(Crystal Studio v.13では2.0となる)を指定します。

数値単位はディフォルトでBohrとHaです。学会発表や論文で使われるSI単位にしたい場合は、数値の後に、空白を入れた後、angstrom と eV をそれぞれ記述します。A*Bとある場合は、BがA個並んでいるという意味になります。

 では次に、出力ファイルを見てみましょう。logにあるWARNINGまたはCOMMENTを見てください。これらのいくつかはユーザーが経験したものです。いまのところそれらは無視してよいでしょう。あなたが自作した入力ファイルで計算が走らないときなどに見てみてください。問題を解決するための手がかりになることが多くあります。
 tbase1_1.outを見てください。次のような情報についての記述があります。
tbase1_1.out のヘッダー
要求されたメモリーの情報(入力したパラメーターでは計算できないことが分かる)
前処理した入力データの表示

chkinp : Checking input parameters for consistency.
という表示があるところまで見てください。この記述までストップしていなければ、読み込まれた入力ファイルに矛盾はないでしょう。この段階で、多くのディフォルト値が提供され、前処理が終了します。

 前処理された入力データをチェックすることはとても大切なことです。tbase1_1.outでは、tbase1_1.inのより詳細な記述になっているはず。もしあなたの想定したような入力データになっていない場合には、tbase1_1.inを見直してみてください。

- "nband": その値は 2 です.
 これは各擬ポテンシャルに記述されているvalenceの値(=価電子の数)をカウントして計算(おそらく合計を1/2倍)して、それにいくつか加えたものです。加える数は系の大きさに依存しますが1以上です。SCF計算はこれで十分です。擬ポテンシャルの読み方が分からない場合には、一度計算を走らせて出力ファイルを調べ、電子状態密度分布(DOS)の計算で、nbandの値を1.2倍程度にしてみましょう。
 余談ですが、似たコードであるPWscfでは、系が金属であれば各擬ポテンシャルに属している価電子の合計の20%を追加しています(最低4が加えられる)。
 mband は nbandの誤りですがあまり気にしなくてよいです。これは昔からあるバグ(=仕様)です。きちんと計算できます。

- "ngfft": その値は30 30 30 です.
 これは3D(三次元)のFFTグリッドの点数です。ecutと単位胞の次元から導かれます。
 
 メモリ評価セクションで指摘される最大の平面波の数 mpw: その値は 752 です。
 これは正確には正しくありません。Abinitコードがk点(0 0 0)で有効である時間反転対称性を利用して1/2倍に減らしているからです。フルセットの平面波の数は1503です。例えば、tbase1_1.out の後の方を見てください。
 Abinitコードは、 通常のディフォルトである istwfk=1 (時間反転対称性を考慮しない)の代わりに、istwfk=2 ( k点(0 0 0)に対する時間反転対称性を考慮する)によって時間反転対称性を示しています。

- "nsym": その値は 16 です.
 系の空間群の番号です。3x3の行列(symrel)で対称操作が定義されます。International Tables Vol. A などを見て対応しているか調べてみるも空間群を理解するのに良いでしょう。tnons で見るように、この例の場合には、並進を伴った対称操作はありません。このAbinitコードは自動的な対称性の解析(nsymが入力ファイルに記述されていないか、ディフォルトの0の場合)を行いました。対称性を探すときには、格子定数や同じ原子の種類での原子位置の差が1.0e-8以内にある場合は同じとみなされます。
 nsymの最大値は384で、これは歪みのない2x2x2のスーパーセルでの対称操作の最大値に対応しています。さらに大きな値を使いたい場合は、abinit.F90の msym での初期設定を変えて、再度コンパイルしてください。

- "xangst" and "xred" はプリミティブセルでの原子位置を記述するするために "xcart" の代替手法として用いられる。

 出力ファイルを読み込んだ後、あなたは次のような疑問に答えられますか?
Q1. toldfe が満たされるために、何回のSCFサイクル(繰り返し計算)が必要となるのか?
Q2. toldfe よりもより収束したエネルギーになっているか?
Q3. 各原子での力(force)の値はHa/Bohr でいくつですか?
Q4. 2つの電子状態の固有値の差はいくらですか?
Q5. e/Bohr^3で電子密度が最大の場所はどこですか? そして、その値はいくつですか?

※ 力(force)というのは物質に働く力のことです。不安定な位置にあるとこの値が大きくなります。

1.2. 原子間距離の計算(手法1)
 2つの水素原子間の最適な距離計算するには3つの手法が存在する。
異なった原子間距離で全エネルギーを計算し、それらをフィッティングして、全エネルギーの最小値である原子間距離を決定する
異なった原子間距離で力を計算し、それらをフィッティングして、力の最小値である原子間距離を決定する
エネルギーを最小化したり、または力が0を探す自動的なアルゴリズムを用いる
原子間距離に対するエネルギーと力の計算を始めましょう。この練習ではどのように多重データセットを使うのかを学びます。

 tbase1_1.inでの原子間距離は1.4 Bohr でした。1.0から2.0 Bohrまで 0.05 Bohrステップで原子間距離を変えて計算しましょう。この計算は21回必要になります。
 もしあなたがUNIXの達人であれば、入力ファイルのxcartを自動的に変えて、21回の計算結果をまとめ、便利な形式でプロットするスクリプトを書くことができれば、このような計算は簡単でしょう。
 あなたはUNIXの達人ですか? もしそうでないなら、より簡単な手法がABINITでは存在します。それは多重データセットです。それに関する詳細な説明はセクション 3.3, 3.4, 3.5 and 3.6 で見つけることができます。
 "irdwfk" and "getwfk" を見てください。特に、getwfk -1 (前回のデータセットで得られた波動関数のファイル _WFKを入力ファイルにするという指示)です。
  ~abinit/tests/tutorial/Input/tbase1_2.in ファイルを用いましょう。tbase1_1.inの時のようにワークディレクトリにコピーしてください。そして、tbase1_x.filesをワークディレクトリにコピーしてxの部分を2に変えてください。tbase1_1.filesのときのように計算してみましょう。
 出力ファイルの後の方に記載されている全エネルギーを調べてみてください。

    etotal1  -1.0368223891E+00
    etotal2  -1.0538645433E+00
    etotal3  -1.0674504851E+00
    etotal4  -1.0781904896E+00
    etotal5  -1.0865814785E+00
    etotal6  -1.0930286804E+00
    etotal7  -1.0978628207E+00
    etotal8  -1.1013539124E+00
    etotal9  -1.1037224213E+00
    etotal10 -1.1051483730E+00
    etotal11 -1.1057788247E+00
    etotal12 -1.1057340254E+00
    etotal13 -1.1051125108E+00
    etotal14 -1.1039953253E+00
    etotal15 -1.1024495225E+00
    etotal16 -1.1005310615E+00
    etotal17 -1.0982871941E+00
    etotal18 -1.0957584182E+00
    etotal19 -1.0929800578E+00
    etotal20 -1.0899835224E+00
    etotal21 -1.0867972868E+00

原文には記載されいませんが、PWscfを用いているプロの間では下記のコマンド(grep)がよく用いられます。

grep etotal tbase1_1.out

これは他の第一原理計算コードを用いた場合に役に立つことが多いです。grep と diff は非常に便利です。

全エネルギーの最小値はデータセット11と12の間にあることが分かります。
     xcart11 -7.5000000000E-01  0.0000000000E+00  0.0000000000E+00
              7.5000000000E-01  0.0000000000E+00  0.0000000000E+00
     xcart12 -7.7500000000E-01  0.0000000000E+00  0.0000000000E+00
              7.7500000000E-01  0.0000000000E+00  0.0000000000E+00

対応する水素の原子間距離は1.5 と 1.55 Bohrです。力の値も1.5 と 1.55 Bohrで消えています。

     fcart11 -5.4945071285E-03  0.0000000000E+00  0.0000000000E+00
              5.4945071285E-03  0.0000000000E+00  0.0000000000E+00
     fcart12  6.9603067838E-03  0.0000000000E+00  0.0000000000E+00
             -6.9603067838E-03  0.0000000000E+00  0.0000000000E+00

これらの2つの値から、線形挿間を用いて、最適値 1.522 Bohrを得ます。
 前回のデータセットから波動関数を読み込んだとき、SCFサイクル(繰り返し計算の回数)は6から5に下がっていることにも注意してください。

1.3. 原子間距離の計算(手法2)
 もう一つの方法は、最小値の自動計算に基づいたものです。それを行うための異なったアルゴリズムがいくつか存在します。"ionmov"が2, 3や7を見てください。今回のケースにおいて、最適化するために1自由度だけ用いる。その最もよい選択はionmov 3 です(自由度が4を越える場合には ionmov 2がより良い)。

 最適化のための時間ステップの最大値を定義しなければなりません。構造最適化の回数 "ntime" を10に設定してください。その値は十分なものです。収束条件 "tolmxf"(あなたが許す絶対値での力の最大値{残留応力})に対しては、5.0d-4 Ha/Bohrが合理的な値です。ディフォルトは5.0d-5 Ha/Bohr ( = 2.5d-3 eV/Angstrom)です。"ntime" になる前に残留力がこの値より小さくなったら構造最適化の計算を止めます(SCFは続きます)。
 各トライアルでの原子間距離で生成された力が十分に収束しているかを確認するために、SCFサイクルに対する収束条件を変えることも価値があります。実際に、toldfe の値、すなわち、1.0d-6 は全エネルギー計算には十分でしょう。しかし、当然ながら他の物性を得るための正確な計算ではありません。そこで、"toldfe" を "toldff" に変えます。 "toldff" の値は"tolmxf"の値よりも1/10にします。~abinit/tests/tutorial/Input/tbase1_3.in はその例になります。一方、~abinit/tests/tutorial/Refs/tbase1_3.out は出力ファイルの例です。計算を走らせると、下記の結果が得られるます。

    etotal   -1.1058360644E+00
     fcart    1.8270533893E-04  0.0000000000E+00  0.0000000000E+00
             -1.8270533893E-04  0.0000000000E+00  0.0000000000E+00
      ...
     xcart   -7.6091015760E-01  0.0000000000E+00  0.0000000000E+00
              7.6091015760E-01  0.0000000000E+00  0.0000000000E+00

 これらのデータに従って、最適な原子間距離は約1.522 Bohです。上記の tbase1_2.out(1.2. 原子間距離の計算(手法1))で見積もった値とよい一致です。
 もしあなたに時間があるならば、収束条件を変えて再計算させて原子間距離がどのように収束していくかをみてください。

1.4. 電子密度の計算
 先に述べたプログラムの実行で全ての幾何学的配置での電子密度が既に得られています。ここでは、この量を出力してみましょう。
 最適化した原子間距離 1.522 Bohr から始めましょう。幾何学的配置を固定して実行します。 "prtden"を 1にセットします(ディフォルトは1ですので記述しなくても電子密度の出力が作成されます)。 "prtden" の記述を正しく理解するために、セクション4にあるようなファイルのより詳細な記述を読むことも価値があります。
 ~abinit/tests/tutorial/Input/tbase1_4.in が電子密度を出力するための入力ファイルになります。
 電子密度は tbase1_xo_DENに出力されます。それを記述するのは……ダメでした! このファイルはASCIIコードを用いて書かれたものではなく、フォーマットされたものではありません。この_DENファイルは、iscf -2 (この指定は、DOSやバンド図の計算のときによく用いられる)にして再計算するときに用いられます。
 読み取るには、"cut3d"と呼ばれるユーティリティーを用いてください。オプション 9を指定して、 XCrysDenのデータを得るようにしてください。

1.5. 原子化エネルギーの計算
 原子化エネルギーは分子をその構成原子になるまで離すときに必要になるエネルギーである。離れた原子は電荷的に中性になる。今回のケースでは、まず始めに分離した原子の全エネルギーを計算する。原子化エネルギーはH2の全エネルギーとHの全エネルギーの2倍との差となる。
 分離した原子の計算にはいくつかの機微が存在する。
多くの場合、分離した原子の基底状態はスピン偏極している。 "nsppol" and "spinat" を見てください。
最も高い占有準位は同じスピンの最も低い非占有準位と縮退している。その場合、通常金属に対して適切なテクニック(lesson_base4)が用いられるべきである。
基底状態の電子密度の対称性は球ではない。その結果、Abinitコードによる対称性の自動決定は機能しなくなる。この場合には、nsym 1 にセットする
 水素については、幸運にも基底状態は1s軌道で球である。そして、最も高い占有準位は同じスピンの最も低い非占有準位は縮退はしているが、異なったスピンである。各スピンの占有状態を occopt (2にセットする) と occ を用いて手動で定義する。
 最終的に、数値誤差が打ち消すように、同じボックス、同じカットオフ、(これはあまり重要ではないけれども)分子の場合と似たボックスでの配置であると考慮して、上記の差を計算することが重要である。
 ~abinit/tests/tutorial/Input/tbase1_5.in は入力ファイルの例である。一方、~abinit/tests/tutorial/Refs/tbase1_5.out は出力ファイルの例です。
 出力ファイルを見て、スピン偏極に関連するとても小さな違いに注意してください。スピンUPとDOWNの双方の電子の固有値はが下記のように与えられます。

 Eigenvalues (hartree) for nkpt=   1  k points, SPIN UP:
 kpt#   1, nband=  1, wtk=  1.00000, kpt=  0.0000  0.0000  0.0000 (reduced coord)
  -0.26414
 Eigenvalues (hartree) for nkpt=   1  k points, SPIN DOWN:
 kpt#   1, nband=  1, wtk=  1.00000, kpt=  0.0000  0.0000  0.0000 (reduced coord)
  -0.11112

 各FFTグリッドでのスピン偏極もまた分析されています。

 Spin up density      [el/Bohr^3]
,     Maximum=    1.4053E-01  at reduced coord.    0.0000    0.0000    0.0000
,Next maximum=    1.2019E-01  at reduced coord.    0.0000    0.0000    0.9667
,     Minimum=    3.4544E-06  at reduced coord.    0.4667    0.4333    0.4333
,Next minimum=    3.4544E-06  at reduced coord.    0.5333    0.4333    0.4333
 Spin down density    [el/Bohr^3]
,     Maximum=    0.0000E+00  at reduced coord.    0.9667    0.9667    0.9667
,Next maximum=    0.0000E+00  at reduced coord.    0.9333    0.9667    0.9667
,     Minimum=    0.0000E+00  at reduced coord.    0.0000    0.0000    0.0000
,Next minimum=    0.0000E+00  at reduced coord.    0.0333    0.0000    0.0000
 Magnetization (spin up - spin down) [el/Bohr^3]
,     Maximum=    1.4053E-01  at reduced coord.    0.0000    0.0000    0.0000
,Next maximum=    1.2019E-01  at reduced coord.    0.0000    0.0000    0.9667
,     Minimum=    3.4544E-06  at reduced coord.    0.4667    0.4333    0.4333
,Next minimum=    3.4544E-06  at reduced coord.    0.5333    0.4333    0.4333
 Relative magnetization (=zeta, between -1 and 1)
,     Maximum=    1.0000E+00  at reduced coord.    0.9667    0.9667    0.9667
,Next maximum=    1.0000E+00  at reduced coord.    0.9333    0.9667    0.9667
,     Minimum=    1.0000E+00  at reduced coord.    0.0000    0.0000    0.0000
,Next minimum=    1.0000E+00  at reduced coord.    0.0333    0.0000    0.0000

 zetaはスピン密度差と電荷密度の比です。それは+1から-1の値になります。今回の水素の場合、スピンDOWNの密度は存在しません。それで、zetaは+1です。
 全エネルギーは、

    etotal   -4.7010531489E-01

一方、H2分子の全エネルギーは、

    etotal   -1.1058360644E+00

 原子化エネルギーは0.1656 Ha です。我々の得た結果と実験値を比べてみましょう。

今回の計算での理論値

  • bond length: 1.522 Bohr
  • atomisation energy at that bond length: 0.1656 Ha = 4.506 eV

実験値

  • bond length: 1.401 Bohr
  • atomisation energy: 4.747 eV

結合長は酷いです(約10%離れている)、そして、原子化エネルギーは僅かに低い(5%離れている)。
 何が間違っていたの??
 あなたは我々が議論しなかった入力パラメーターが正しいか確認しました? それは、

  • ecut 平面波の運動エネルギーに関するカットオフ
  • acell セルのサイズ
  • ixc いままで指摘していないけど、交換相関汎関数に何を用いたかを明記する
  • 擬ポテンシャル

 私たちは カットオフエネルギーとして10 Ha、10x10x10 bohr^3のスーパーセル、Teterがパラメーター化した局所密度近似(ixc 1)、Goedecker-Hutter-Teter tableからの擬ポテンシャルを用いました。
 次のレッスンで、これらのパラメータの選択の仕方についての取扱い方を説明します。

----------------------------------------------------------------------------------------------------------------------------

The lesson 2 はH2分子を再び取り扱う: 収束の研究と、LDA 対 GGA

2.1 と 2.2 ecut での収束
 結合と原子化エネルギーを何度も計算しているために、全ての関連する操作を行える単一の入力ファイルを作るのは有用です。2つのデータセット( ~abinit/tests/tutorial/Input/tbase1_3.in と ~abinit/tests/tutorial/Input/tbase1_5.in を結合させる)を用いてトライしてみてください。その入力ファイルが、 ~abinit/tests/tutorial/Input/tbase2_1.in になります。一方、~abinit/tests/tutorial/Refs/tbase2_1.out は出力ファイルです。
 計算を走らせると、次の値が得られます。

    etotal1  -1.1058360644E+00
    etotal2  -4.7010531489E-01

そして、

    xcart1  -7.6091015760E-01  0.0000000000E+00  0.0000000000E+00
             7.6091015760E-01  0.0000000000E+00  0.0000000000E+00

1実行で得た結果なのであるけれども、レッスン1で決定した値に類似しています。残留力が5.0d-4よりも小さいこともチェックしてください。収束の議論は abinit_help ファイルの セクション 7でなされています。それを読んでください。ところで、abinit_help ファイルは読んでくださいね。
 多くの収束パラメーターは既に確認されています。ecutとacellを強制的に変えます。それは次の理由からです。

  • SCFサイクルと幾何学的配置の収束は toldfe, toldff 及び tolmxf (これは他の物理的な物性に対するケースとはならない)で制御できます。
  • 大きなボックスでの孤立系で行われるようなk点の収束性の研究はない。1を越えるようなk点で得られる付加的な情報はない。
  • ボックスカットの値は自動的に2より大きい値がABINITによって選ばれる。前処理で決められたngfftの値を見てください。
  • 幾何学的配置の決定のために、ionmov=3 を用います。

 ecutに関連する収束のチェックのために、tbase2_1.in でecutの値を変えて実行を行うか、~abinit/tests/tutorial/Input/tbase2_2.inで提案されている2重ループを動作させるかの選択肢があります。ecutの値は10 Ha と 35 Haの間で 5 Haのステップで選択されます。もし、あなたが2重ループを作りたいのであれば、abinit_help ファイルにある double-loop section を再度読むことは有益です。
 計算をするために必要な計算時間とメモリが大きく増大するのが分かるでしょう。出力は次のようになります。

    etotal11 -1.1058360644E+00
    etotal12 -4.7010531489E-01
    etotal21 -1.1218716100E+00
    etotal22 -4.7529731401E-01
    etotal31 -1.1291943792E+00
    etotal32 -4.7773586424E-01
    etotal41 -1.1326879404E+00
    etotal42 -4.7899908214E-01
    etotal51 -1.1346739190E+00
    etotal52 -4.7972721653E-01
    etotal61 -1.1359660026E+00
    etotal62 -4.8022016459E-01

     xcart11 -7.6091015760E-01  0.0000000000E+00  0.0000000000E+00
              7.6091015760E-01  0.0000000000E+00  0.0000000000E+00
     xcart12  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart21 -7.5104912643E-01  0.0000000000E+00  0.0000000000E+00
              7.5104912643E-01  0.0000000000E+00  0.0000000000E+00
     xcart22  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart31 -7.3977108831E-01  0.0000000000E+00  0.0000000000E+00
              7.3977108831E-01  0.0000000000E+00  0.0000000000E+00
     xcart32  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart41 -7.3304273322E-01  0.0000000000E+00  0.0000000000E+00
              7.3304273322E-01  0.0000000000E+00  0.0000000000E+00
     xcart42  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart51 -7.3001570260E-01  0.0000000000E+00  0.0000000000E+00
              7.3001570260E-01  0.0000000000E+00  0.0000000000E+00
     xcart52  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart61 -7.2955902118E-01  0.0000000000E+00  0.0000000000E+00
              7.2955902118E-01  0.0000000000E+00  0.0000000000E+00
     xcart62  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00

 対応する原子化エネルギーと原子間距離は、

ecut    atomisation   interatomic distance
(Ha)    energy (Ha)      (Bohr)

10       .1656          1.522
15       .1713          1.502
20       .1737          1.480
25       .1747          1.466
30       .1753          1.460
35       .1756          1.459

 結合長と原子化エネルギーで、0.2%の相対精度を得るためには、カットオフが30 Haとなることが分かります。この最後の実行での値を覚えておいてください。
 うーん、30 Ha は大きなカットオフです。水素で用いる擬ポテンシャルは非常に硬い……。

2.3 acell での収束
 ecutで行ったのと同様なテクニックを acell の収束に用いていみましょう。8 8 8 から 18 18 18 まで 2 2 2 のステップで探索します。 この研究では、ecutは 10に保持します。実際、ecutの収束とacellの収束の間の影響はほとんどありません。これはどちらかといえば一般的です。 ~abinit/tests/tutorial/Input/tbase2_3.in が入力ファイルの例です。出力データ(~abinit/tests/tutorial/Refs/tbase2_3.out)は次のようになります。

    etotal11 -1.1188128741E+00
    etotal12 -4.8074164402E-01
    etotal21 -1.1058360644E+00
    etotal22 -4.7010531489E-01
    etotal31 -1.1039109441E+00
    etotal32 -4.6767804802E-01
    etotal41 -1.1039012761E+00
    etotal42 -4.6743724199E-01
    etotal51 -1.1041439319E+00
    etotal52 -4.6735895176E-01
    etotal61 -1.1042058195E+00
    etotal62 -4.6736729718E-01

     xcart11 -7.8427127284E-01  0.0000000000E+00  0.0000000000E+00
              7.8427127284E-01  0.0000000000E+00  0.0000000000E+00
     xcart12  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart21 -7.6091015760E-01  0.0000000000E+00  0.0000000000E+00
              7.6091015760E-01  0.0000000000E+00  0.0000000000E+00
     xcart22  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart31 -7.5472588642E-01  0.0000000000E+00  0.0000000000E+00
              7.5472588642E-01  0.0000000000E+00  0.0000000000E+00
     xcart32  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart41 -7.5491961064E-01  0.0000000000E+00  0.0000000000E+00
              7.5491961064E-01  0.0000000000E+00  0.0000000000E+00
     xcart42  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart51 -7.5427730675E-01  0.0000000000E+00  0.0000000000E+00
              7.5427730675E-01  0.0000000000E+00  0.0000000000E+00
     xcart52  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
     xcart61 -7.5415348434E-01  0.0000000000E+00  0.0000000000E+00
              7.5415348434E-01  0.0000000000E+00  0.0000000000E+00
     xcart62  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00

対応する原子化エネルギーと原子間距離は、
 

acell
(Bohr)
atomisation
energy (Ha)
interatomic distance
(Bohr)
8
.1574
1.568
10
.1656
1.522
12
.1686
1.509
14
.1691
1.510
16
.1694
1.508
18
.1695
1.508

原子間距離で0.2%の収束に到達するためには、 acell が 12 12 12 が必要です。原子化エネルギーはそのレベルで収束するために acell 14 14 14 が必要になります。12 12 12で、差は 0.0009 Ha = 0.024 eVです。それは実用的には十分に小さいです。最終的な実行として、 acell 12 12 12 を用います。
 ほとんどの固体において、単位胞のサイズはそれよりも小さいくなります。このスーパーセルにおいて多くの真空層が取り扱われています。それで、この擬ポテンシャルを用いたH2研究は非常に簡単ではないことが明らかです。

□ 2.4 L(S)DA での最後の計算
 今から ecut と acell の正しい値を用います。ええっと、acell 12 12 12 と ecut 30で計算出来るように、tbase2_3.in を書き換えてください。そのファイルの書き変えを最小化するために、udtset 1 2 (それは1ループに減らせます)のある2重ループ機構をまだ使うことができます。~abinit/tests/tutorial/Input/tbase2_4.in が入力ファイルの例です。~abinit/tests/tutorial/Refs/tbase2_4.outが出力ファイルの例です。
 単一のペア(ecut, acell)で計算を実行してから、CPU時間は前回の一連の計算を経た最適な値の決定に費やした時間よりも多くはありません。しかしながら、必要なメモリは増加する。
 出力データは、

    etotal11 -1.1329372051E+00
    etotal12 -4.7765320721E-01

     xcart11 -7.2662135225E-01  0.0000000000E+00  0.0000000000E+00
              7.2662135225E-01  0.0000000000E+00  0.0000000000E+00
     xcart12  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
  • 対応する原子化エネルギーは 0.1776 Ha = 4.833 eV
  • 原子間距離は 1.4532 Bohr.
  • これはL(S)DAを用いた我々の最終データです.

ixc = 1を用いました。他の局所密度近似の表式 ixc = 2, 3 ... 7 が可能です。量子モンテカルロ計算を基にして、均一電子ガスのXCエネルギーからそれらは全て出発しているために、値 1, 2, 3と7はほぼ同じ結果を与えます。他の可能な ixc = 4, 5, 6はより古い電子密度汎関数で、これらのデータを頼りにしていません。

2.4 GGAの利用
 Phys. Rev. Lett. 77, 3865 (1996) で提案されたPerdew-Burke-Ernzerhof 関数を用います。
 理論上は、GGAに対しする擬ポテンシャルを用いることが必要です。しかしながら、水素のような特殊なケース、そして非常に小さなコアを持つような(コアに1s 軌道だけを持つような)一般的な擬ポテンシャルでは、LDAとGGAから生じた擬ポテンシャルはよく似ています。
 それで、擬ポテンシャルを換えません。これは、ecutの収束テスト(ecutは60%程度の頻度で計算で用いる擬ポテンシャルの特徴となる)を再度しなくてもよいために、私たちに多くの時間を省かせてくれます。
 擬ポテンシャルに関係無く、acellの収束テストは再び行いません。それは、真空ではLDAまたはGGAにおいて似た取扱いがなされるためです。
 それで、tbase2_4.inで ixc を 11 にセットすることで、GGAの範囲内で最終的な値は容易に得られます。~abinit/tests/tutorial/Input/tbase2_5.in に例がありますので見てください。

    etotal11 -1.1621428502E+00
    etotal12 -4.9869631917E-01

     xcart11 -7.1204535408E-01  0.0000000000E+00  0.0000000000E+00
              7.1204535408E-01  0.0000000000E+00  0.0000000000E+00
     xcart12  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
  • 対応する原子化エネルギーは 0.1647 Ha = 4.482 eV
  • 原子間距離は 1.4241 Bohr.
  • GGAでの最終的なデータである.
もう一度, ここに実験データを記載してみよう。
  • 結合長: 1.401 Bohr
  • 原子化エネルギー: 4.747 eV

 GGAにおいて、実験での結合長の2%の範囲内にある。しかし、実験での原子化エネルギーの5%である。LDAにおいて、実験での結合長の4%の範囲内にある。しかし、実験での原子化エネルギーの2%である。
 LDAとGGAの典型的な精度は計算する系の種類によって変化することを忘れないでください。

※ とはいえ、最も高精度であるWIEN2kはGGAでのPBEを用いている。PAW擬ポテンシャルは、PBE用に作られているので、この事実は認識しておいて、PBEを用いると良い。
※ PWscf は PAWでもLDAなどが豊富である。ファンデルワールス力が無視できな系では、LDAでの計算が最低でも2012年ごろまで行われてきた。

----------------------------------------------------------------------------------------------------------------------------

The lesson 3 はSi結晶(絶縁体)を扱う: k点グリッドの定義、カットオフエネルギーのスメアリング、バンド構造の計算、そして、再び、収束の研究。

3.1. 固定したk点数でのSiの全エネルギーの計算
 ~abinit/tests/tutorial/Input/tbase3_x.files と ~abinit/tests/tutorial/Input/tbase3_1.in を用います。これまでは分子や孤立原子を計算してきました。Siなどの結晶では、tbase3_1.in の中身を見ると分かるように、結晶の計算で必要になる新しいvariableの記述があります。

  • rprim 格子のベクトルを記述するものです。cif2cellを使えば自動で記述されるのであまり意識する必要はありません。
  例を一緒に見てみましょう。
  rprim a1 a2 a3
            b1 b2 b3
            c1 c2 c3
  と並びます。acell A B C と記述されていれば 格子ベクトル a, b, c は下記になります。
  a = a1 * A + a2 * B + a3 * C
  b = b1 * A + b2 * B + b3 * C
  c = c1 *  C+ c2 * B + c3 * C
  これが難しい場合は、angdeg を用いてください。angdeg α β γ と記述するので分かり易いですね。
  • xred (xcart の代わりに用いられる) 原子位置を a, b, c 軸の長さを1として記述します。cif2cellを使えば自動で記述されるのであまり意識する必要はありません。
  • kptopt, ngkpt, nshiftk, shiftk, kptrlatt, (簡単ではないです ... あせらないでじっくりやってください)
  • diemac (孤立した分子と比較して, 他の値が用いられます。一方、 diemix は抑制される).

~abinit/tests/tutorial/Input/tbase3_1.in を走らせると出力ファイルから全エネルギーが得られます。

   etotal   -8.8662238960E+00


3.2. k点数に関する収束の研究の開始
 ここではブリュアンゾーンのサンプリングに関係する収束の研究を行います。分解能を増加せるように異なったグリッドを用います。次のような一連のグリッドです。

ngkpt1  2 2 2
ngkpt2  4 4 4
ngkpt3  6 6 6
ngkpt4  8 8 8

しかしながら、規約表現上のブリュアンゾーンにおけるk点数は非常に大きくなります。それは下記のようになります。

nkpt1  2
nkpt2 10
nkpt3 28
nkpt4 60

 ABINITは自動的にこのk点数をグリッドの定義と対称性から計算します。
 入力ファイルでnkptを定義するにも関わらず、グリッドから計算された値と入力ファイルを比較します。問題が検出されたとき、ABINITの振る舞いがどのようになるかを見るためにこの機会を利用しましょう。では、ngkept1 4 4 4 を提案して、nkp1 2 としてみましょう。入力ファイルの例は ~abinit/tests/tutorial/Input/tbase3_2.in です。logファイルの終わりから数十行のところにメッセージがあります。

--- !BUG
message: |
    The argument nkpt=     2, does not match
      the number of k points generated by kptopt, kptrlatt, shiftk,
      and the eventual symmetries, that is, nkpt=    10.
      However, note that it might be due to the user,
      if nkpt is explicitely defined in the input file.
      In this case, please check your input file.
src_file: getkgrid.F90
src_line: 415
...

  Action : contact ABINIT group.

 これは典型的なABINITのエラーメッセージです。それはABINITがストップした原因であるどのような問題を示しています。入力ファイルでのエラー、つまり、nkpt の誤った値をどのように対処すべきかを提案します。予測値は nkpt 10 が誤った入力ファイルのエラーを注意する前に記述されています。それから、問題が生じたソースファイルがその行数も加えて指摘されています。

 明確なk点グリッドに関するnkptの計算が容易な仕事でなく、明確で効率的なグリッドの重要な選択がさらに一層難しいために、ユーザーへのいくつかのヘルプがABINITから提供されます。ABINITは自動的に異なったk点グリッドを分析し、ベストなグリッドを提案することができます。abinit_helpに記載されています。prtkpt と kptrlenで記述される積分精度の関連情報を見てください。k点セットのリストの生成は ~abinit/tests/v2 にある異なったテストケースで行われます。

 新しい材料の研究を始めるときに、k点グリッドのリスト、そして、最低でも3つの効率的なグリッド、k点での収束の研究を最初に調べることを強くアドヴァイスします。取り扱うk点数にCPU時間が比例することを忘れないでください。10 kでは2k の5倍の時間が必要になります。

3.3. k点数に関する収束の実際の動作
 k点グリッドを理科するために、Monkhost and Pack の論文(Monkhorst and Pack paper, Phys. Rev. B 13, 5188 (1976) )を読んでください。ええっと、たぶん直ぐには……その一方で、先に述べた収束の研究をしてみましょう。

 入力ファイル ~abinit/tests/tutorial/Input/tbase3_3.in が例になります。一方、~abinit/tests/tutorial/Refs/tbase3_3.out が参考になる出力ファイルです。この出力ファイルで k点(nkp)とそれらの重み(wtk)のリストを生成するために ngkpt と shiftk が特に用いられています。kptとwtkの情報を読んでください。
出力ファイルから、単位胞あたりの全エネルギーの評価が下記のようにあります。

    etotal1  -8.8662238960E+00
    etotal2  -8.8724909739E+00
    etotal3  -8.8726017432E+00
    etotal4  -8.8726056405E+00

データセット3とデータセット4の間の差は非常に小さいです。データセット2で0.0001 Haの精度が得られているのだけれども、acellとecutを固定しての全エネルギーの収束値は -8.8726 Ha です。

※ ここまで読んできて、難しいなと思った人も多いのではないでしょうか? 全エネルギーを原子数で割って、1 meV/atom 内であればそこでk点を増やす必要はありません。入力ファイルの誤りもABINITが指摘してくれます。格子定数をn倍したときには、k点数を1/n倍にすればよいので、簡単な系で計算して、大きな系にしてみるも良いでしょう。

3.4. 格子定数の決定
 "optcell" はセルの形状と体積の自動最適化を支配する。
 セル体積の自動最適化は、次のvariablesを用いる。

optcell 1
ionmov 3
ntime 10
dilatmx 1.05
ecutsm 0.5

dilatmx and ecutsm について指示されていることを読み取ってください。

※ なんということだ! これでは初学者は参ってしまう。

  • dilatmx は格子定数の許される最大のスケールを与える。ディフォルトは1.0です。dilmatxを大きくするとCPUとメモリの無駄を導く。acellが10%内で最適化されると考えた場合には、単純に dilatmx 1.1 を用いてください。
  • ecutsm は optcell が 0でないときに重要です。ディフォルトは1.d0です。ecutsmが0でないとき、歪み矛盾しないで、ecut または acell の関数としての全エネルギー曲線が滑らかになります。推奨値は0.5です(単位はHa)。


 全てのk点についてはテストしませんが、ngkpt  2 2 2  と  4 4 4 の2つのみ行います。
 ~abinit/tests/tutorial/Input/tbase3_4.in は入力ファイルの例です。一方、~abinit/tests/tutorial/Refs/tbase3_4.out は出力ファイルの例です。
 次のような格子定数の漸近的な変化を得ます。

     acell1   1.0233363682E+01  1.0233363682E+01  1.0233363682E+01 Bohr
     acell2   1.0216447241E+01  1.0216447241E+01  1.0216447241E+01 Bohr

 次のような非常に小さな残留応力も得られます。

    strten1   1.8591719160E-07  1.8591719160E-07  1.8591719160E-07
              0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    strten2  -2.8279720007E-08 -2.8279720007E-08 -2.8279720007E-08
              0.0000000000E+00  0.0000000000E+00  0.0000000000E+00

 歪みテンソルは Ha/Bohr^3単位で得られます。そして、成分の順序は下記のようになります。

                        11  22  33
                        23  13  12

 acell1 と acell2 間の相対的な差は0.1%だけです。
 それで、14si.pspnc 擬ポテンシャルを用いたSiに対しするLDA値は 10.216 Bohr(実際には 10.21644)で、5.406 Angstrom でした。実験値は 25度で5.431 Angstrom です。

3.5. バンド構造の計算
 acell を理論値の 3*10.216 に固定します。k点(4x4x4 FCCグリッド、これは8x8x8 Monkhorst-packグリッドと等しいです)も固定します。
 8バンドで計算します(4価電子帯と4伝導帯)。

 バンド構造はブリュアンゾーンの異なったラインにそって多くの異なったk点に対するコーンシャム方程式を解くことによって計算することができます。
 コーンシャム方程式に入力するポテンシャルはこの前にSCF計算から算出されます。
 L-Gamma-X-(U)-Gammaの経路でそれぞれ10, 12, 17 で分割した場合でのバンド構造を計算します。
 この経路の各点はk空間上で下記のようになります。

L     (1/2 0 0)
Gamma (0 0 0)
X     (0 1/2 1/2)
Gamma (1 1 1)

Note:

  1. 最後のGamma点は最初のGamma点より逆格子空間の他のセルにあるものです。この選択はX-U-Gammaの経路を容易に構築させます。
  2. 還元した座標を用いてk点を明確にします。標準的な教科書にあるプリミティブ 2(原子単位セル)の入力設定と一致します。コンベンショナルな8原子の座標でL, Gamma, または X点が与えられることが分かります。上記の経路では、(1/2 1/2 1/2)-(0 0 0)-(1 0 0)-(1 1 1)、しかし、2-原子-セルではそのようなコンベンショナルな座標は使うことができない。

それで、入力ファイルにセットアップします。最初のデータセットとして通常のSCF計算、次のデータセットとして下記のvariablesを用います。

  • iscf -2, SCFをしない計算(バンド構造やDOSを計算する時に指定される) ;
  • getden -1, 1 データセット1からの電子密度の出力を用いる;
  • nband 8 ;
  • kptopt -3, ブリュアンゾーンで3つの領域を定義する;
  • ndivk 10 12 17 (3つの領域をそれぞれ10、12, 17 分割する)
  • kptbounds
                    0.5  0.0  0.0 # L point
                    0.0  0.0  0.0 # Gamma point
                    0.0  0.5  0.5 # X point
                    1.0  1.0  1.0 # Gamma point in another cell.
    
  • enunit 1, 固有値の出力をeV単位にする
  • SCF計算をしないときに許される収束条件は tolwfr です。tolwfr 1.0d-10 としてください。 これはtoldfe を抑制します。

 ~abinit/tests/tutorial/Input/tbase3_5.in が入力ファイルの例です。一方、 ~abinit/tests/tutorial/Refs/tbase3_5.out は出力ファイルです。
 2番目のデータセットが始まってバンド構造を見つけられます。

 Eigenvalues (   eV  ) for nkpt=  40  k points:
 Eigenvalues (   eV  ) for nkpt=  40  k points:
 kpt#   1, nband=  8, wtk=  1.00000, kpt=  0.5000  0.0000  0.0000 (reduced coord)
  -3.78815  -1.15872   4.69668   4.69668   7.38795   9.23867   9.23867  13.45707
 kpt#   2, nband=  8, wtk=  1.00000, kpt=  0.4500  0.0000  0.0000 (reduced coord)
  -3.92759  -0.95774   4.71292   4.71292   7.40692   9.25561   9.25561  13.48927
 kpt#   3, nband=  8, wtk=  1.00000, kpt=  0.4000  0.0000  0.0000 (reduced coord)
  -4.25432  -0.44393   4.76726   4.76726   7.46846   9.31193   9.31193  13.57737
 kpt#   4, nband=  8, wtk=  1.00000, kpt=  0.3500  0.0000  0.0000 (reduced coord)
  -4.64019   0.24941   4.85732   4.85732   7.56855   9.38323   9.38323  13.64601
 ....

 これらのデータを表示するためにグラフィカルツールが必要になります。別のファイル_EIGで、k点と固有値のリストを得ることができます。
 グラフィカルツールが無くても、L, Gamma, X と Gammaでの値を見ることができます。

 kpt#   1, nband=  8, wtk=  1.00000, kpt=  0.5000  0.0000  0.0000 (reduced coord)
  -3.78815  -1.15872   4.69668   4.69668   7.38795   9.23867   9.23867  13.45707

 kpt#  11, nband=  8, wtk=  1.00000, kpt=  0.0000  0.0000  0.0000 (reduced coord)
  -6.17005   5.91814   5.91814   5.91814   8.44836   8.44836   8.44836   9.17755

 kpt#  23, nband=  8, wtk=  1.00000, kpt=  0.0000  0.5000  0.5000 (reduced coord)
  -1.96393  -1.96393   3.00569   3.00569   6.51173   6.51173  15.95524  15.95524

 kpt#  40, nband=  8, wtk=  1.00000, kpt=  1.0000  1.0000  1.0000 (reduced coord)
  -6.17005   5.91814   5.91814   5.91814   8.44836   8.44836   8.44836   9.17755
 最後の Gamma は厳密には最初のGammaと等しいです。価電子帯のトップがGamma点であることを確認してください。価電子帯の幅は12.09 eVで、最も低い非占有状態はX点で0.594 eVです。この非占有状態はgamma点にある価電子帯のトップよりも高い。Siは間接遷移型のバンドギャップ(これは正しい)でそのバンドギャップが0.549 eV(これは定量的には全く誤っている。実験値は25度で1.17 eVである)であるとして記述されている。伝導帯の最小値はX点に対して僅かにズレている(kpt #21を見る)。このバンドギャップの過小評価は(有名なDFTバンドギャップ問題)としてよく知られている。正しいバンドギャップを得るために、コーンシャム密度汎関数法を越えるためにGW近似を用いる必要がある。これは the first lesson on the GW approximation に記述されている。
 実験データとバンド構造の説明は下記を見ると良い。
M.L. Cohen and J.R. Chelikowski
Electronic structure and optical properties of semiconductors
Springer-Verlag New-York (1988).
----------------------------------------------------------------------------------------------------------------------------

The lesson 4 はAl結晶(金属)とその表面を取り扱う: 占有数、フェルミディラック関数でのスメアリング、表面エネルギー、そして、再び、収束の研究。

※ occupt と tsmear を導入して、ブロードニングをかけた場合の調査となります。金属では収束をよくするためにこのブロードニングが用いられます(金属では占有準位と非占有準位が近くにあり、SCFサイクル毎に入れ替わってなかなか収束しません。これを導入してフェルミ端近傍を問題ない範囲でブロードニングして滑らかにすることで収束を早くすることができます)。下記を読むことは有用ですが、JTHのPAW擬ポテンシャルを用いる場合は、その論文に記載されているフェルミディラック分布のスメアリングとなる occopt 3 と tsmear 0.002 (Ha) を設定すれば良いでしょう。この設定で計算を走らせながらレッスン4を読んでいくのも良いでしょう。

□ 4.1. 固定したスメアリングとk点数でのAlの全エネルギーと格子定数の計算
□ 4.2. k点数に関する収束の研究
□ 4.3. k点数とブロードニングファクターに関する収束の研究

※ これ以降は界面を行う方向けになる。
□ 4.4. Al(100)の表面エネルギーの決定: 単位胞の方位の変化
□ 4.5. 表面エネルギーの決定: (3Al + 1真空層)スラブ計算
□ 4.6. 表面エネルギーの決定: 真空層の数の増加
□ 4.7. 表面エネルギーの決定: Al層の数の増加

----------------------------------------------------------------------------------------------------------------------------
 


アクセス数
ページビュー数