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

ABINIT(spin, DFT+U, GW, BSE, RF)

 必要な部分のみ記載されている。下記の部分だけでも多くの計算が可能である。
----------------------------------------------------------------------------------------------------------------------------

The lesson on spin in ABINIT  はスピンに関連したプロパティを提供する: スピン偏極した計算とスピン軌道相互作用。

□ 1. 強磁性金属: BCC Fe
 スピンを入れた計算をする場合には nsppol で 2 を指定します。

 spinat は初期の磁化の方向を指定します。繰り返し計算(SCF計算)をしている間に、最適な磁化の方向が決められます。では、なぜ初期値を設定しておくことが必要かというと、電子のスピンが揃った初期配置とバラバラなスピンの初期配置で計算すると最終的な結果が異なる場合があるためです。

□ 2. 反強磁性金属: FCC Fe
 

  • nspden ディフォルト値は nsppol の値となっているのであまり深刻になる必要はありません。
  下記の場合に注意してください。
 ・ general collinear magnetization
    nspden 2
    nsppol 2
    nspinor 1
 ・ antiferromagnetism (反強磁性)
    nspden 2
    nsppol 1
    nspinor 1
 ・ vector magnetization
    nspden 4
    nsppol 1
    nspinor 2
 (nspden=2 は電荷密度について2つの独立した成分を持つことを意味します。一方、nsppol=1は平面波について1つの独立した成分を持つことを意味します)
 

□ 3. FCC Fe をもう一度見直す
□ 4. 強磁性の例

□ 5. スピン軌道カップリング
 ノルム保存及びウルトラソフト擬ポテンシャルを用いた場合には下記をvariableで2を設定します。

 擬ポテンシャルにPAWを用いるときには下記のvariableで1を設定します。

 この場合、nspinorは自動的に2にされます。

□ 5. 磁化とスピン軌道回転

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

The lesson on DFT+U はABINITを用いてDFT-Uをどのように実行するかを示す。そして、NiOの電子状態密度分布(DOS)を計算へと導く。

□ 0. DFT+U法の短い要約

 Abinitにて, 2つの二重カウンティング補正が備えられている :

-The Full localized limit (FLL) (see A. Lichtenstein et al PRB 52, 5467 (1995))

-The Around Mean Field (AMF) (see M. T. Czyzyk and G. A. Sawatzky PRB 49, 14211 (1994))

 いくつかの系に対して、この結果はおそらく2重カウンティングの方法に依存する。しかしながら、2つの手法は一般的に似た結果を与える。


□ 1. LDAを用いたNiOの基底状態計算
 NiOでは反強磁性(antiferromagnetism)での計算になっています。そのため、nspden, nsppol, nspinor が明記されています。上記の反強磁性金属: FCC Feを見てください。


□ 2. FLL 二重カウンティングでのDFT+Uを用いたNiOの基底状態計算
 tldau_1.in が基底状態計算の入力ファイルであり、これに下記の部分を書き加えたのが tldau_2.in となります。variableの記述から分かるとおり、PAWの擬ポテンシャルを用いた入力ファイルになっています。

# LDA+U
usepawu   1
lpawu   2 -1
upawu  8.0 0.0 eV
jpawu  0.8 0.0 eV

 jpawu = upawu * 0.1 とします。

  • usepawu DFT+Uの計算方法に2種類の方法が用意されている。
   1の場合: Full Localized Limit (FLL) (または Atomic Limit) 二重カウンティングが用いられる
   2の場合:  Around Mean Field (AMF) 二重カウンティングが用いられる
      AMFは金属や弱い強相関系に対して安定に動作する
   ※ WIEN2kは強相関系やフルポテンシャル法で安定なSICを用いており、J=0とする。
  • lpawu DFT+Uを導入する軌道を記述します。-1はDFT+Uを導入しない。2はd軌道、3はf軌道に導入します。
  • upawu U
  • jpawu J


□ 3. NiO結晶に対するLDA+Uでの密度行列の初期化
□ 4. AMF 二重カウンティングでのDFT+Uを用いたNiOの基底状態計算
□ 5. LDA+Uでの状態密度の射影

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

The lesson on the use of PAW (PAW1)  はノルム保存擬ポテンシャルの代わりとしてABINITで実装されたPAW法を提供する。PAWはCPU時間でかなり有利な点を持つ。

□ 0. PAW法の要約

□ 1. ABINITでのPAWの利用

□ 2. 平面波のカットオフ(ecut)での収束
 通常の場合におけるように、収束パラメーターは平面波基底のサイズを定義しているカットオフです。

1.a ノルム保存ケースにおけるダイヤモンドでのecutの収束性を計算する
 入力ファイル tpaw1_2.in はダイヤモンド(実験値の体積)に対するecutでの収束を計算するためのデータが含まれています。8 Haから 24 Ha まで 2 Haステップでecutを増加させるために、9データセットがあります。
 tpaw1_2.files (標準的なノルム保存擬ポテンシャル)を用いて実行させてください。

abinit < tpaw1_2.files > tmp-log

 次のような値が得られます(出力ファイル tpaw1_2.out)。

    etotal1  -1.1628880677E+01
    etotal2  -1.1828052470E+01
    etotal3  -1.1921833945E+01
    etotal4  -1.1976374633E+01
    etotal5  -1.2017601960E+01
    etotal6  -1.2046855404E+01
    etotal7  -1.2062173253E+01
    etotal8  -1.2069642342E+01
    etotal9  -1.2073328672E+01

1 mHaレベルで、etotal の収束がecut=24 Haでも達していないことを確認してください。

1.b PAWケースにおけるダイヤモンドでのecutの収束性を計算する
 セクション1.a.と同様に、同じ入力ファイルを用います。
 再び、tpaw1_2.files の最後の行を、6c.pspnc から 6c.lda.atompaw に変更します。実行して、その出力ファイル(tpaw1_2.outA)を開いてみてください。
  次の値が得られます。
    etotal1  -1.1474828697E+01
    etotal2  -1.1518675625E+01
    etotal3  -1.1524581240E+01
    etotal4  -1.1525548758E+01
    etotal5  -1.1525741818E+01
    etotal6  -1.1525865084E+01
    etotal7  -1.1525926864E+01
    etotal8  -1.1525947400E+01
    etotal9  -1.1525954817E+01

 確認してみてください。
 1 mHa での etotal の収束は 12 <= ecut <= 14 Ha (etotal4 は 最後の値の1m Ha以内です)
 0.1 Ha での etotal の収束は 16 <= ecut <= 18 Ha (etotal6 は 最後の値の0.1m Ha以内です)
 
 それで、同様な入力ファイルを用いて、ダイヤモンドでのPAW計算は、ノルム保存擬ポテンシャル計算と比べてより低いカットオフを要求します。

※ PAWはウルトラソフト擬ポテンシャルと似ているために、一般的に、ノルム保存擬ポテンシャルよりも少ないカットオフで高い精度を実現できます。
※ このチュートリアルで分かると思いますが、PAWの計算は .filesでの擬ポテンシャルをPAWに指定しなおすだけです。非常に簡単に高精度で計算できるようになりました。
※ 2014年に、VASPに匹敵する精度を持つPAW擬ポテンシャルが(http://www.abinit.org/downloads/PAW2 /JTH-TABLE/index.html)にて提供されています。JTHは論文(フリーでありがたい)を見ると、ecut は 15 から 20 Haで良いことが分かります。15 Ha では Δ1が急激に小さくなって 2.187 meV、20 Ha では Δ1が 0.888 meVです。20 Haと40 Haは大きな違いがありません。通常は ecut 15 Ha を用い、He, B, C, N, O, F と Ne (気になる場合は遷移金属も含め)では ecut 20 Ha を用いれば良いでしょう。k点数は 6750/原子数 にし、 フェルミディラック分布のスメアリングとなる occopt 3 と tsmear 0.002 (Ha) を設定します。私は入力ファイルにpawecutdg を明記しなかったのですが、出力ファイルには 60 Ha が設定されていました。しかし、pawecutdgにecutと同じ値を入力しても同じ結果でした。Fe2VAlでの擬ギャップの再現性はDOS計算でのk点数の違いによるものでした。
※ ABINITがまだVASPに負けている主な要因としては、構造最適化の速度とそのデータの引き継ぎ(時間を気にしなければ良いだけですが……)、フォノンの計算が素人では非常に困難であるなどでしょう。

□ 3. ダブルグリッドFFTカットオフ(pawecutdg)での収束
※ pawcutdg は ecut で自動的に決められますので、いまのところ省略します。

□ 4. DOSへのPAWの寄与をプロットする
□ 5. PAWの部分波基底の完全性をテストする
□ 6. PAWの結果の妥当性のチェック
□ 7. ABINITにおけるPAWの追加のコメント
7.a SCFサイクルでのミキシングの手法; 全エネルギーの分解
7.b PAW球のオーバーラップ
7.c PAWでの体積のプリント
7.d 追加のPAW入力variables
7.e PAW+U
 研究している系が強相関電子を含んでいるならば、LDA+Uが便利です。usepawulpawu, upawu and jpawu で制御できます。
 (省略)
 PAW+Uが50%の頻度で難しくなるが、SCFサイクルの収束、または、大域的最小点を保証が達せられる。usedmatpudmatpawu の使用はそれを助ける。

----------------------------------------------------------------------------------------------------------------------------
The first lesson on GW はプラズモンポールモデルを用いてGW近似(コーンシャムLDAバンド構造よりもより良い)でSi(半導体)の準粒子バンドギャップの計算を取り扱う。
□ 1. GW近似を用いたGamma点でのSiのバンドギャップの計算
□ 2. 収束の研究のための準備: コーンシャム構造(KSS file)とスクリーニング(SCR file)
□ 3. 自己エネルギーを計算するための波動関数における平面波の数に基づく収束(任意)
□ 4. Σxを計算するための波動関数における平面波の数による収束性(任意)
□ 5. Σcを計算するための波動関数における平面波の数による収束性(重要)
□ 6. スクリーニングを計算するための波動関数における平面波の数による収束性(任意)
□ 7. スクリーニングを計算するためのバンドの数による収束性(重要)
□ 8. スクリーニングマトリクスの次元での収束性(重要)
□ 9. Gamma点におけるバンドギャップに対するGW補正の計算
□ 10. GWコードにおける高度な特徴

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

The second lesson on GW はプラズモンポールモデルを用いずGW近似(コーンシャムLDAバンド構造よりもより良い)でAlの準粒子バンドギャップの計算を取り扱う。
 このレッスンではGW近似の範囲内でDFTコーンシャム固有値の自己エネルギー補正を得る方法に焦点を当てる。金属の場合においては、プラズモンポールモデルを用いない。Alのバンド幅とフェルミエネルギーを計算する。

□ 1. 予備的なコーンシャムバンド構造計算
□ 2. スクリーニングファイルの計算
□ 3. フェルミエネルギーと価電子バンドの底を探す
□ 4. GW特殊関数の計算とAlのプラズモンサテライト

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

Parallelism of Many-Body Perturbation calculations (GW) は正確な電子構造の計算(多体効果を含んだ準粒子バンド計算)をスピードアップさせる。
 このレッスンでは、ABINITのGW部分をどのように並列計算するかを示すことに狙いを定める。私たちは典型的なGoWo計算の異なったステップでの並列化手法を議論していく。ABINITでのN CPUの並列化計算は次のコマンドでなされる。
  mpirun -np N abinit < abinit.files
 

1. 並列計算でのKSSの作成 (約 179.8秒)
 GWチュートリアルの最初のレッスンで、私たちはどのようにコーンシャム構造ファイル(KSS file)を生成するかを学びました。いま、私たちは基底状態部分により実行さえるk点並列を利用して同様な計算を行う。
 まず最初に、tmbt_1.fileをワーキングディレクトリ(Work_mbt)にコピーする。

cd cd $HOME/abinit-7.6.4/tests/tutoparal/Input
mk Work_mbt
cd Work_mbt
cp ../tmbt_1.files . 

 入力ファイル ~abinit/tests/tutoparal/Input/tmbt_1.in を好みのエディタで開いて、その構造を見てください。
 最初のデータセットは基底状態の電子密度を得るためにSCF計算が実行されます。二番目のデータセットは電子密度ファイルを読み込んで、多くの空の状態を含んだコーンシャムバンド構造を計算します。

# DATASET 2 : KSS ファイルを作成する
iscf2      -2      # NSCF計算
getden2    -1      # 前回計算された電子密度を読み込む
tolwfr2    1d-12   # NSCFサイクルでの収束条件
kssform2    3      # 共役勾配法アルゴリズム (大規模系に対して勧められるオプション)
nband2      160    # NSCFサイクルで(占有と空の)バンドを計算する数
nbdbuf2     10     # バッファに対するバンド数。大きなバッファはNSCF計算のステップ数を減少させるのに役立つ
nbandkss2   150    # KSSファイルに格納されるバンド数 (収束した状態のみ記述される)

 私たちはすでにGWチュートリアルの最初のレッスンでこれらのvariableに出くわしています。
 KS ハミルトニアン(kssform=1, これはディフォルトです)の直接対角化を採用する代わりに、共役勾配法(kssform=3)でNSCFサイクルを解くことは価値があります。
 直接対角化の並列計算は非効率的でメモリも多く必要とすることから、私たちは共役勾配法を選択します。
 NSCF計算はk点とスピンでの標準的な並列処理で実行されます。バンドの(nkpt x nsppol)ブロックがノード(=CPU)に配分されます。このテストは シフトしない4x4x3グリッド(ブリュアンゾーン中の48 k点、規約表現上のくさびにおける9 k点に折りたたまれる)を用います。それゆえ、理論的な最大のスピードアップは9です。
 ABINIT を N CPU コアで次のコマンドを用いて走らせてみましょう。
※ NはCPUのコア数よりも小さくて3の倍数の値を入力してください。

mpirun -np N abinit < tmbt_1.files >& tmbt_1.log &

 しかし、心に留めておいてください。プロセッサのアイドルを避けるために、CPU数は9に分割しなければならない。実行の最後で、コードは次のGW計算で必要とする tmbt_1o_KSSを生成します。
 3つのノード(=CPU)を用いて、wall clock time は約1.5分です。

 tail tmbt_1.out

-
- Proc.   0 individual time (sec): cpu=        209.0  wall=        209.0

================================================================================

 Calculation completed.
.Delivered    0 WARNINGs and   5 COMMENTs to log file.
+Overall time at end (sec) : cpu=        626.9  wall=        626.9

 参考の出力ファイルは、~tests/tutoparal/Refs, under the name tmbt_1.out で得られます。
 収束したGW結果をえるために、150 バンドは十分ではありません。あなたの計算資源に比例させてバンド数を増やしてください。

2. Adler-Wiser展開を用いた並列計算でのスクリーニングの計算 (約384.8秒)
 チュートリアルのこのパートで、私たちはAdler-Wiser法を用いてRPA分極率を計算します。その基本式はGWノートのこのセクションで議論します。
 最初に、 tmbt_2.file をワーキングディレクトリにコピーしてください。それから、先のステップで作成されたKSS ファイルへのシンボリックリンクを作ります。※ WFKもシンボリックリンクを作らないと動作しません。

 ln -s tmbt_1o_DS2_KSS tmbt_2i_KSS 
 ln -s tmbt_1o_DS2_WFK tmbt_2i_WFK 

 入力ファイル ~abinit/tests/tutoparal/Input/tmbt_2.in を開いて、その構造を議論します。
 スクリーニング計算をコントロールしているパラメータのセットが下記にまとめられています。

optdriver   3   # スクリーニングの計算
irdkss      1   # KSS ファイルを読み込む
symchi      1   # ブリュアンゾーンでの積分の高速化を図るために対称性を考慮に入れる
awtr        1   # 時間反転を利用する。Mandatory when gwpara=2 is used.
gwpara      2   # 並列化をバンドに対して行う
ecutwfn     24  # 波動関数のカットオフ
ecuteps     8   # 分極率のカットオフ
nband       50  # (24の占有バンド) RPA展開におけるバンドの数
inclvkb     2   # Correct treatment of the optical limit.

 variableの大部分は既にGWチュートリアルで議論しています。いくつか追加の説明をするにふさわしいvariableは gwparaawtr です。
 gwpara はスクリーニングを計算するために必要な並列化のアルゴリズムを選択します。二つの異なった手法が実行されます。

  • gwpara=1 → ブリュアンゾーンのk点に関する並列化
  • gwpara=2 → メモリ配分をともなったバンドに関する並列化

 各手法は利点と欠点を有する。それはvariableのドキュメントで議論されている。このチュートリアルで、私達はgwpara=2に焦点を当てる。これは最も良いMPIスケーラビリティをともなったアルゴリズムである。そして、一番肝心で、要求されるメモリの顕著な減少を許す唯一のものです。
 オプション awtr=1はシステムが時間反転対称性を示すことを明確にする。その結果、明確に計算される遷移(共鳴遷移のみ要求される)の数を半分にする。gwpara=2を用いるとき、awrt=1は必須です。

 並列計算を走らせる前に、実装の重要な技術の詳細を議論するのは価値がある。我々の目的に対して、gwpara=2がスクリーニングの部分で用いられるということで充分である。そのコードは各プロセッシングユニットが占有バンドのフルセットを所有するような波動関数を配分する。一方、空の状態はノード(CPU)間で配分される。
1. 各ノードはRPA分極率への部分的な寄与を計算する
2. 各ノードの部分的な結果が集められる。
3. マスターノードは逆誘電行列を得るために逆行列を計算し、最終的な結果をファイルに書く。
 アルゴリズムの最初と二番目のステップ両方で、プロセッサ数と対応することが期待される。ステップ3では、対照的に、シーケンシャルに実行される。それゆえ、全ての規模で、特に、大きなスケーリング行列(大きな npweps または ω点での大きな周波数 )の場合において、それは有害な効果を持つ。
 使用される最大のCPU数は分極率を計算するために用いられる空の状態の数によって決定される。実際、この行列は配分されない。それゆえに、各ノードは(npweps2 x nomega x 16 bytes)で与えられるテーブルを格納するための十分なメモリを有しなければならない。ここでnomega は計算される全ての周波数の数である。
 下記の図で示されるBarcelona Supercomputing Centerで実施されたテストではMPIアルゴリズムの一番目と二番目の部分が非常によいスケーリングを持っていることが示される。このルーチン cchi0 と cchi0q0、そこではRPA展開(ステップ1と2)の計算が512プロセッサまで線形に増加する。多くのプロセッサで観測された全てのスピードアップの低下は主に並列化されていない計算の部分によるものである。すなわち、KSS ファイル(rdkss)と 逆行列(qloop)の読み込みである。
 その点で、実行での最も重要な技術の詳細をカバーしなければならない。我々は最終的に次のコマンドを用いてN CPUでABINITを計算させる。※ こちらのNはCPUのコア数よりも小さければよいようです。3や4でも動作。

mpirun -np N abinit < tmbt_2.files >& tmbt_2.log &

 異なったプロセッサの数を用いて入力ファイル tmb_2.inを計算させる。そして、各プロセッサ数での計算時間の記録をつける。その結果、我々は実行でのスケーラビリティをテストする。上記の図で要求されたパフォーマンス分析はテストケースとして、ZnOを用いてPAWで得られている。しかし、SiO2でのような類似の振る舞いが得られる。

 出力ファイルの結果を見てみましょう。このチュートリアルは主にどのように効率的なMPI計算を走らせるかに焦点をあわせているから、SiO2に対する収束の研究を行わない。入力ファイルで用いられる多くのパラメータは既に収束しているに近い。k点サンプリングと空の状態の数が増加するのみである。GWチュートリアルの最初のレッスンで記述されている手順を追うようにして標準的な収束のテストを実行するように入力ファイルを書き換えてみてください。
 主な出力ファイルで、どのようにバンドがノード間で配分されるかを報告しているセクションが存在する。シーケンシャル計算においては、我々は下記を得る。

 screening : taking advantage of time-reversal symmetry
 Maximum band index for partially occupied states nbvw =    24
 Remaining bands to be divided among processors   nbcw =    26
 Number of bands treated by each node ~   26

 最後の行で報告された値は、計算が多くのノードで計算されたとき減少する。
 プロセッサ数での波動関数の増加に対してメモリが割り当てられる。logファイルからこの情報を引き出すために、grepを用いることができる。シーケンシャルでの計算について、我々は下記を得る。

$ grep "Memory needed" tmbt_2.log

  Memory needed for storing ug=         29.5 [Mb]
  Memory needed for storing ur=        180.2 [Mb]

 ug は軌道のフーリエ成分を格納するために用いられる内部のバッファを意味する。軌道のフーリエ成分のサイズはnpwwfnで線形的に増加する。urは実空間のFFTメッシュ上での軌道を格納している配列である。urのサイズはFFTボックス内での全k点数に線形に増加するということを記憶に留めておいてください。数は、普通、平面波 (npwwfn).の数より大きい。GWコードにおけるFFT分割数は次のコマンドを用いて主要な出力ファイルから引き出される。

$ grep setmesh tmbt_2.out  -A 1
 setmesh: FFT mesh size selected  =  27x 27x 36
          total number of points  =    26244

 GWノートのこのセクションにおいて議論されたように、高速フーリエ変換は実行の最もCPUの部分の集中的な部分の一つである。この理由から、このコードは入力 variale fftgw を提供する。 fftgw はバッファの効率化のためにFFT点の数を減少させるために用いられる。入力 variable gwmem の二つ目の数値は、代わりに、実空間の軌道の格納を支配する。そして、計算時間の増加の費用でコストの高い配列urの格納をさけるために用いられる。

 スクリーニング計算で要求される計算の努力はq点数で線形的に増える。GWノートのこのセクションで示されたように、このコードはスクリーニング関数の対称性を利用する。その結果、規約表現上のブリルアンゾーンが明確に計算されるようになる。一方、大きなq点数は収束した結果に達することが求められる。典型的な例は、金属におけるGW計算、または、Bethe-Salpeter 方程式を用いた光物性である。

 もし、十分なプロセッシングユニットが利用可能であれば、q点サンプリングによる線形のファクターは、q点の計算をvariable nqptdm and qptdm を用いて各個の実行に分割することで明確に吸収される。その結果は mrgscr を用いることで一つのバイナリファイルに集められる。

3. ヒルベルト変換法を用いた並列計算でのスクリーニングの計算 (約651.4 秒)
※ 上記の2.より下記の理由で効率的な手法で同じ計算をしているので、1.で計算して得られた同じ出力ファイルを用いている。
 GWノートで議論されたように、多くの周波数が求められるとき、Adler-Wiser展開を基にしたアルゴリズムは最善のものではない。この段落では、高密度の周波数メッシュでRPA分極率を計算するために、どのようにヒルベルト変換法を用いるかを議論する。そのコードに実装した方程式はGWノートのセクションで文書化されている。
 通常、私達は tmbt_3.file をワーキングディレクトリにコピーする。そして、KSSファイルのシンボリックリンクを作る。

 ln -s tmbt_1o_DS2_KSS tmbt_3i_KSS 
 ln -s tmbt_1o_DS2_WFK tmbt_3i_WFK 

 入力ファイルは ~abinit/tests/tutoparal/Input/tmbt_3.in である。それを開いて、我々はその構造を見る。
 アルゴリズムを支配しているもっとも重要なパラメータのスナップショットが下記で示される。

gwcalctyp   2    # GWの種類を選択する。2はワンショトでのContour-deformation techniqueによるGW計算
spmeth      1    # 1はスペクトル法、ディフォルトの0はAdler-Wiserでの計算
nomegasf  100    # スペクトル関数の点数
gwpara      2    # 2は並列化をバンドに対して行う。ディフォルトの1は並列化をk点に対して行う。
awtr        1    # 時間反転を利用する。Mandatory when gwpara=2 is used.
freqremax  40 eV # 自己エネルギーの数値積分を実行するための誘電体行列を計算するために用いられる最大の実周波数
nfreqre    20    # 実軸上の周波数メッシュ数
nfreqim     5    # 自己エネルギーの数値積分を実行するための誘電体行列を計算するために用いられる最小の実周波数

 入力ファイルはAdler-Wiser計算で用いられたものに似ている。入力 variable spmeth はスペクトル法を可能にする。nomegasf はスペクトル関数で用いられる線形メッシュでのω′点数、つまり、スペクトル関数に対する方程式でのω′数を定義する。
 省略
 次のコマンドを用いて、N CPUコアでABINITを実行させる。

mpirun -np N abinit < tmbt_3.files >& tmbt_3.log &

 プロセッサ数を変えてスケールをテストする。CPU数が伝導の状態数を分割するとき、計算での仕事の配分は良いバランスである。記憶に留めておいてください。
 スペクトル関数を格納するために要求されるメモリはログファイに記載されています。

$ grep "sf_chi0q0" tmbt_3.log
 memory required by sf_chi0q0:           1.0036 [Gb]

 多くのプロセッサが用いられたとき、どのようにこの配列のサイズが減少するかに注意してください。
 下記の図はAlder-Wiserとhilbert変換法を用いて計算されたSiO2の電子エネルギー損失関数を示します。これらの結果を再現してみてください。(EELFはtmbt_3o_EELFファイルに記載されている。収束に到達するには、より多くのk点サンプリングが要求される)

4. 並列計算でのワンショットGW補正の計算 (約280.6秒)
 この最後の段落で、我々はgwpara=2を用いて並列化してどのようにG0W0補正を計算するかを議論する。自己エネルギー行列要素を計算するために基礎方程式がGWノートのこの部分で議論される。
 計算を走らせる前に、 tmbt_4.fileをワーキングディレクトリにコピーする。それから SCRとKSSファイルに対してシンボリックリンクを作る。

ln -s tmbt_1o_DS2_KSS tmbt_4i_KSS 
ln -s tmbt_2o_SCR     tmbt_4i_SCR
ln -s tmbt_1o_DS2_WFK tmbt_4i_WFK 

※ 金属では下記を用います。
ln -s tmbt_1o_DS2_KSS tmbt_4i_KSS 
ln -s tmbt_3o_SCR     tmbt_4i_SCR
ln -s tmbt_1o_DS2_WFK tmbt_4i_WFK 

 入力ファイル ~abinit/tests/tutoparal/Input/tmbt_4.in を開いてみましょう。
 計算に最も重要なパラメーターは下記で示されています。

optdriver   4            # 4は自己エネルギー計算でのルーチン"Sigma"を実行する.
irdkss      1            # KSSファイルを読み込む
irdscr      1            # SCRファイルを読み込む
gwcalctyp   0 ppmodel 1  # ppmodelのディフォルトである1はプラズモンポールモデル。金属では除く。
#gwcalctyp  2            # 金属の場合に用いる。ヒルベルト変換法を用いた並列計算でのスクリーニングの計算を行う。
gwpara      2            # 2は並列化をバンドに対して行う。ディフォルトの1は並列化をk点に対して行う。
symsigma    1            # 1とすると自己エネルギー行列要素の対称化を可能にする。ディフォルトは0
ecutwfn    24            # 波動関数のカットオフ
ecuteps     8            # 相関部分のカットオフ
ecutsigx   20            # 交換部分のカットオフ
nband       50           # 相関部分のバンド数

 我々の目的に対して、入力ファイルではプラズモンポールモデル近似を用いた標準的なone-shot 計算を定義することで十分であると言える。
 この場合においても、我々は並列計算を実行するために gwpara=2 を用いる。しかしながら、自己エネルギー部分で用いられる軌道の分布は、スクリーニング計算に用いたものから僅かに異なっている。
 以下では、私たちは簡単に波動関数の配分するために用いられている2つのステップのプロセスを簡単に記述する。

  1. 各ノードはQP補正が計算された状態を読み込み、そして、メモリ内に格納する。 (kptgw and bdgw によって明確にされた状態のリスト).
  2. バンド(nband でバンド数が定義される)は次のパーティションの概念を用いて分割される。

 採用された部分的な分割の効力によって、相関部分の計算がCPU数でスケールすることが予測される。使用されるプロセッサの最大数は、 nband によって制限される。しかしながら、ステップ2においてバンドが分割されたとき、プロセッサのサブセットのみが占有状態を受け取っていることに注意してください。結果ととして、交換部分における理論的な最大スピードは、実行に含まれる異なったMPIノード上での占有状態の利用によって制限される。
 最良の場合のシナリオはQP補正が全ての占有状態を要求したときです。このケースにおいては、代わりに、各ノードは自己エネルギーの部分を計算することができ、そして、ほとんど線形のスケーリングが達成される。最悪の場合のシナリオでは、準粒子補正はわずかな状態(例えばバンドギャップ計算)で Ncpu >> Nvalence についてのみ要求されるときである。このケースにおいては、代わりに、Nvalence プロセッサのみが交換部分の計算に参加する。

 要約しますと: プロセッサ数が nband を割れるとき、相関部分のMPI計算は効率がよい。交換部分の最適なスケーリングは占有状態のフルセットを各ノードが占有するときのみ得られる。このステップは最悪の場合のシナリオに近いということに注意してください。自己エネルギー行列要素(csigme)の計算は64プロセッサまでよくスケールする。

 下記の二つの図はプロセッサ数の関数としてシグマ部分のスピードアップを示す。自己エネルギーは nband = 1024 (205 占有状態)を用いて5つの準粒子状態に対して計算される。CPUの数が大きとき、スケーリングはバランスを崩した占有状態の配分のために線形の振る舞いから離れる。実行でのスケールのない部分(init1, rdkss)はアムダール則に従って全体のスピードアップが制限される。

 大規模な配列が配分されたために、その実装は良いメモリのスケーラビリティを示しています。スクリーニングのサイズはノード数で増加していない。ディフォルトで、各CPUは最適な計算を得るためにメモリに全てのq点と周波数についての全てのスクリーニング行列を格納する。大規模行列の場合、単一のq点のみがメモリに格納されるout-of-core solutionを選ぶことが可能である。そして、データは外部のSCRファイルから読み込まれる(遅い、しかしメモリ要求はより小さい)。
 どのように負荷を効率的に配分するかを知るために、次のコマンドを用いて最終的に計算を走らせましょう。

mpirun -np N abinit < tmbt_4.files >& tmbt_4.log &

 各プロセッサ数での時間を記録を保持してください。それで自己エネルギー部分のスケーラビリティをテストすることができる。これらのテストの結果は収束し ていないことに注意してください。よく収束した計算はブリュアンゾーンでサンプルした6x6x6のk点メッシュとスクリーニング行列に対して10 Haが求 められます。そのQP(準粒子)結果は空の状態の数に関して非常にゆっくりと収束します。0.1 eV精度の範囲内でQPのギャップを収束させるために、スクリーニングにて1200 バンド、自己エネルギーの計算に800の状態を含まなければならなかった。
 α-SiO2のLDAバンド構造とG0W0エネルギーバンドとの比較は下の図で示されている。Gamma点での直接ギャップは、湾ショットG0Wo法を用いたとき、LDAの6.4 eVから約9.5 eVまで異常に開いている。この結果を再現することが勧められる(理論的なLDAパラメーターで実行されたこの計算を考慮している。一方、このチュートリアルでの全ての入力ファイルは実験で得られた結晶構造が用いられている)。
※ tmbt_4.inの nkptgw  5 -> nkptgw 9, bdgw  24 25 -> gdgw1 30 としても動作した。その場合の計算時間は 9798.2秒であった。
※ バンド図は Perturbative Calculationでのk点と最後のEの値をプロットすればよい。
※ PAW(JTK)でも動作する。


5. 効率的な並列計算のための基本的なルール

  1. 「上手くいかない可能性のあると、動かない」ことは覚えておいてください。それで、あなたが入力ファイルを書くとき、短くシンプルであるようにしてみてください。
  2. 一つのことをして、そして、それを上手に行う:
  3. 同様な入力ファイルにおける optdriver の異なった値の使用を避ける。各ランレベルはメモリとCPUタイムを配分するために違ったアプローチを採用し、それ故、各データセットでバランスの取れた実行をするプロセッサ数を見つけることはほとんど不可能です。
  4. 素数定理:
  5. テストしたパラメーターがMPIアルゴリズムと影響がないときだけ、並列計算での収束の研究が実行されなければならない。例えば、スクリーニングにおけるバンド数での収束の研究は gwpara=2 を用いたとき分離した入力ファイルで実行されなければならない。
  6. 過ぎたるは及ばざるがごとし(少ない方がよい):
  7. 大規模計算は可能ならばいつでも小さな実行に分割する。例えば、スクリーニング計算は、q点で分割することが出来る。自己エネルギーの計算は容易に kptgwbdgw で分割できる。
  8. 慎重に:
  9. 依存するCPU時間と要求されたメモリがどのようにテストされたパラメーターに依存するかを評価するために収束の研究を用いてください。収束したパラメーターで最終計算に着手するとき、計算資源の評価を用いることは非常に有用です。

----------------------------------------------------------------------------------------------------------------------------
The lesson on BSE はBethe-Salpeter式を用いてSiの巨視的な誘電関数の計算を取り扱う。

1. 準備段階(KSSとSCRファイルの生成)
 下記のコマンドでワークディレクトリを作成してください。

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

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

cd Work_bs

 このディレクトリの中に下記のコマンドで入力ファイルをコピーします。

cp ../tbs_1.files .

 次のコマンドで直ぐに計算を始めましょう。

abinit < tbs_1.files >& tbs_1.log &
 
 その結果、コードが実行されている間に、入力ファイルを解析してみましょう。
 入力ファイルは
~abinit/tests/tutorial/Input/tbs_1.in に配置されています。ヘッダーは計算の簡潔な記述がされています。そのため、それを注意深く読んでください。これから地道に計算について議論していきますから、もし幾つかの部分が明確にならなくても気にしないで下さい。

 入力ファイルは次のBethe-Salpeter計算のために必要になる2つのKSSファイルとscreeningファイルを生成します。最初のデータセットはk点をシフトさせない4x4x4のグリッドでの基底状態の計算を行います(ブリルアンゾーンには64 点があります。規約表現上では8 点になります)。それから、NSCF計算で2つのKSSファイルを生成するために、データセット2と3で基底状態の計算で得られた電子密度分布を用います。その計算には共役勾配法が用いられます(
kssform=3)。

 データセット2で計算されたKSSファイルはGamma点を中心とした4x4x4のk点メッシュで、100バンドを含んでいることに注意してください。一方、データセット3で作られたKSSファイルはメッシュの点の位置を( 0.11 0.21 0.31 )だけシフトさせた4x4x4のメッシュで10バンドだけ有しています。


shiftk3    0.11 0.21 0.31  # This shift breaks the symmetry of the k-mesh.

 Gamma点を中心としたk点メッシュは規約表現上の領域で8点を含んでいます。一方、シフトさせたk点メッシュは規約表現上でのブリルアンゾーンで64点となっている結晶の対称性を壊します(ブリルアンゾーンと規約表現上のブリルアンゾーンは一致します)。シフトさせたk点メッシュは明らかに非効率的です。そのため、あなたはなぜこの奇異なサンプリングを用いたのか、そして、2つの異なるKSSファイルの生成を必要としたのか不思議に思うでしょう。

 実際に、この手法はGWチュートリアルで行われた方法から非常に異なっています。しかし、そのようなことを行う正当な理由が存在します。光学スペクトルはブリルアンゾーンのサンプリングでゆっくりと収束し、そして、基底状態やGW計算で一般的に用いられる対称的なk点メッシュよりも対称性を壊したk点メッシュは
nkpt において速い収束を導くことが予測されます。

 これが奇異なシフトを用いている理由です。しかし、何故2つのKSSファイルを必要とするのでしょうか。何故、私たちは単純にスクリーニング計算するためにシフトしたk点でのKSSファイルを用いないのでしょうか?

 この理由は、シフトしたk点上で多くの空のバンドを用いて行われる
スクリーニング計算が多くのメモリを占有するためです。そのサイズは( nband * nkpt )で大きくなります。そして、k点数を減少させるために用いられる対称性は少ない。

 要約すると、対称性のあるk点サンプリングと100バンドによるKSSはスクリーニングを計算するために用いられる。一方、シフトしたk点と10バンドによるKSSはBethe-Salpeter方程式を解くために用いられるtransition spaceを構築するために用いられる。その2つのk点メッシュはちょうどシフトだけ異なっていて、それゆえ、それらはq点の同じセットを作る(スクリーニングにおけるq点のリストはKSSファイルのk点の間での全ての可能な差として定義される)。

 これは、BSの実行において、transition space はシフトしたk点メッシュで構築されるのだけれども、 私たは対称的なメッシュを用いて生成されたSCRファイルを用いることを意味します。
 このかなり技術的な点を明らかにするために必要とされた長い議論の後、私たちは tbs_1.inの最後のデータセットで実行されるスクリーニング計算を分析するために漸く前進する。
 SCRファイルは
nband=100 and ecuteps 6.0 Haを用いてデータセット4で計算される。GWチュートリアルの最初のレッスンにおいて、0.0 1eVの範囲内で収束したQPエネルギーを得るためにこれらの値が見つけられる。それで、私たちはSCRファイルが収束し、そして、Bethe-Salpeter部分における収束テストを実行するためにそれが安定に用いられることを確信しています。
 合理的な理由から、Wのstatic limitが計算されます。

nfreqre4  1     # Only the static limit of W is needed for standard BSE calculations.
nfreqim4  0

 実際、Bethe-Salpeter方程式の標準的な形式において、スクリーニングした相互作用のstatic limitはBSハミルトニアンのクーロン部分を構築するために必要とされます。 単一の周波数を用いることは私たちにスクリーニング部分でのCPU時間を節約することを許します。しかし、このSCRファイルは、ppmodel=3,4に対応しているプラズモンポールモデルを用いているBethe-Salpeter計算かGW計算のどちらかで用いられるだけであることを記憶に留めておいてください。
 この時点で、この計算は完了します。しかし、いまだ私たちは次の段に移動する前にしなければならないことが一つあります。
 私たちが言ったように、私たちはシフトしたk点メッシュ上でのKSSファイルとBS計算に必要なSCRファイルを必要とします。それらを消さないように。より意味のある名称を用いて、これらの貴重なファイルの名称を変えるのも良いアイディアです。例えば、

mv tbs_1o_DS2_KSS 444_gamma_KSS
mv tbs_1o_DS3_KSS 444_shifted_KSS
mv tbs_1o_DS3_WFK 444_shifted_WFK
mv tbs_1o_DS4_SCR 444_SCR

 今後、k点サンプリングはもはや変えることはないことは記憶に留めておいてください。BS入力で記述されるk点のリストは、KSSファイルを生成するために用いられたものと等しくなければなりません。もし私たちがtransition spaceを構築するために用いられるk点サンプリングを変えたいならば、2つの新しいKSSファイルと一つの新しいSCRファイルはscrathから生成しなければならない。

2. Tamm-Dancoff近似による吸収スペクトルの計算
 このセクションではどのようにHaydock繰り返し法を用いてTamm-Dancoff近似(TDA)による標準的な励起子の計算を実行するかを示すことを意図している。その入力ファイルは~abinit/tests/tutorial/Input/tbs_2.inです。
 このジョブを走らせる前に、私たちはtbs_1.inで生成された出力結果を用いてこの計算に接続しなければならない。
 Unixコマンドを用いて、

ln -s 444_shifted_KSS tbs_2i_KSS
ln -s 444_shifted_WFK tbs_2i_WFK
ln -s 444_SCR tbs_2i_SCR

 KSSとSCRファイルに対する2つのシンボリックリンクを作ります。これを行う理由は私たちが入力ファイルを議論した後に明らかになるでしょう。
 このジョブは最近の計算機で1-2分で足りる。それで、入力ファイルを調査する前に、それを走らせることは価値があります。
 ~abinit/tests/tutorial/Input/tbs_2.filesファイルをワーキングディレクトリにコピーします。そして、計算を走らせます。

abinit < tbs_2.files >& tbs_2.log &

 バックグラウンドにジョブをおきます。その結果、私たちは tbs_2.in を分析することができます。l
 あなたの好きなエディターで~abinit/tests/tutorial/Input/tbs_2.inを開いてください。そして、次のセクションに進みます。そこで、私たちは典型的なBS計算を支配しているもっとも重要なvariableを議論します。
 最初に、私たちはBSEルーチンを呼ぶために optdriver=99 の設定を必要とします。

optdriver  99        # BS calculation

 variable  irdkss and irdscr  はABINITの"ird" variable と似ています。そして、それら先の段落で生成されたファイルを読み込むために用いられます。

irdkss  1     # Read the KSS file produced in tbs_1
irdscr  1     # Read the SCR file produced in tbs_1

 そのコードは入力ファイルKSSファイルとSCRファイルを見つけることを期待します。それらのファイルの名称は、tbs_2.files で明記したプレフィックスに従って構築されます。これは私たちがこのコードを走らせる前に2つのシンボリックリンクを作らなければならない理由です。
 それから、私たちは励起子ハミルトニアンをどのように構築するかを明記する5つのvariableのリストを有しています。

bs_calctype       1       # L0 is constructed with KS orbitals and energies.
soenergy          0.8 eV  # Scissors operator used to correct the KS band structure.
bs_exchange_term  1       # Exchange term included.
bs_coulomb_term  11       # Coulomb term included using the full matrix W_GG'
bs_coupling       0       # Tamm-Dancoff approximation.

 その bs_calctype=1 は独立粒子の分極率がKSSファイルから読まれるコーンシャム軌道とエネルギーで構築される。自己エネルギーをシュミレートするために、KSエネルギーは、エネルギー soenergy= 0.8 eV の scissors operator で補正される。これは私たちにtransition spaceに含まれる各状態に対する面倒なGW計算を回避させることを許します。scissors operator の利用は、Siに対して合理的な近似です。しかし、初期のKSバンド構造の単純な剛体シフトの観点からGW補正がシミュレートできないようなより複雑な系でそれは失敗します。
 残りの3つのvariable はどのように励起子ハミルトニアンを構築するかを明記します。bs_exchange_term=1 はプログラムに要点である交換部分を計算することを教えます。それゆれ、この計算は局所場効果を含んでいます。variable bs_coulomb_term は異なったオプションを選択するために用いられる。そのオプションはクーロン項に対して利用できる(どうか幾らか時間をかけてvariableの記述とBS_notes に関連する記述を読み込んでください)。最後に、bs_coupling=0 は非対角で結合しているブロックを無視します(Tamm-Dancoff 近似)。このパラメーターの個々の組み合わせは、局所場効果を含んだTamm-Dancoff近似の範囲内でBethe-Salpeter計算に対応しています。
 それから、私たちはtransition spaceを構築するために用いるバンドも明記も有します。

bs_loband         2
nband             8

 この場合、ギャップ周りの全てのバンド、そのインデックスは2から8の間です、基底形に入れられます。
 巨視的な誘電関数ついての周波数メッシュは bs_freq_mesh  を通して明記されます。

bs_freq_mesh 0 6 0.02 eV  # Frequency mesh.

 この3つの実数の項は[0, 6] eV範囲で0.02 eVステップの線形メッシュを定義します。そのメッシュ中にある周波数の数はCPU時間に重大な効果を有することはありません。k点数と連動して周波数領域が記述される。結果として bs_loband and nband は正確な収束の研究の対象となる。
 それから、私たちは巨視的な誘電関数を計算するために用いられるアルゴリズムを定義し、コントロールするパラメーターを有する。

bs_algorithm        2      # Haydock method.
bs_haydock_niter   100     # Max number of iterations for the Haydock method.
bs_haydock_tol     0.05    # Tolerance for the iterative method.
zcut               0.15 eV # Complex shift to avoid divergences in the continued fraction.

 bs_algorithm は巨視的な誘電関数を計算するために用いられるアルゴリズムを明記する。この場合、私たちは繰り返しの haydock 法を用いる。その繰り返しの最大数はbs_haydock_niter によって与えられる。その繰り返しのアルゴリズムは光学スペクトルの2つの連続的な評価の間での差がbs_haydock_tol より小さくなったとき停止する。その入力 variable zcut は連分数における発散を避けるために複素(僅かな虚軸方向)へのシフトを与える。物理的な観点から、このパラメーターは吸収ピークの実験的なブロードニングを模擬する。このテストにおいて、k点メッシュの粗さが原因で、私たちはHaydockアルゴリズムの収束を促進するためにディフォルトの0.1 eVよりも僅かに大きな値を用いなければならない。理想的に、k点数の増加に対する zcut  の値を減少させて収束の研究をしなければならない。
 k点サンプリングはvarialbeのセットによって記述される。

# Definition of the k-point grid
kptopt 1                  # Option for the automatic generation of k points,
ngkpt  4 4 4              # This mesh is too coarse for optical properties.
nshiftk 1
shiftk    0.11 0.21 0.31  # This shift breaks the symmetry of the k-mesh.
chksymbreak 0             # Mandatory for using symmetry-breaking k-meshes in the BS code.

 kptopt, ngkpt, nshiftk, and shiftk の値はKSSファイルについてグリッドを記述するために用いられた値と等しくなければならない。chksymbreak=0 は対称性が壊れた上でのチェックを避けるために用いられる。そうしないとコードはストップする。
  入力ファイルの最後のセクションは
 
ecutwfn 8.0               # Cutoff for the wavefunction.
ecuteps 2.0               # Cutoff for W and /bare v used to calculate the BS matrix elements.
inclvkb 2                 # The Commutator for the optical limit is correctly evaluated.
 
 カーネル行列要素とダイポールオペレーターの行列要素を計算するために用いられるパラメーターを明記する。私たちはGWの最初のレッスンでのこれらのvariableに既に出くわしている。それで、それらの意味は詳しいと思う。BSコードでのこれらのvariableによって振る舞われる役割のより詳細な議論はBS_notesで見つけられる。
 
 出力ファイルtbs_2.outは計算と繰り返し法が収束していないなら報告されるVARNINGsの基本的なパラメーターを報告している。それらの構造を理解するために時間を使ってください。
 次の質問に答えられますか?
 

  1. 基底関数に含まれる遷移は幾つですか?
  2. 光学的な限界を評価するために幾つの方向を用いますか?
  3. 連続的なfractionで用いられるローレンチアンブロードニングの値は何ですか?
 
 主要な出力ファイルでのこの脱線の後、私たちはその計算の出力データを解析するために話を進めます。
 5つの異なったファイルに最も重要な結果が格納されます。
 
• tbs_2o_BSR
• tbs_2o_HAYDR_SAVEhandy
• tbs_2o_RPA_NLF_MDF
• tbs_2o_GW_NLF_MDF
• tbs_2o_EXC_MDF
 
 In what follows, 私たちはフォーマットと各出力ファイルの内容の簡単な記述を提供する。
 
• tbs_2o_BSR:
 
 このバイナリファイルは共鳴ブロックの上三角形(エルミート行列である。それ故、非冗長部分が計算される。そして、ファイルにセーブされる)を格納する。そのBSRファイルはgetbsreso or irdbsresoを用いて先の計算から再実行するために用いられる。もし、収束が達成しなかったり、または、異なったzcut値でHaydock計算を実行したりするならば、この再計算の可能性はHaydock 法を再計算するのに便利である。共鳴と結合ブロックを格納するために、そのコードは2つの異なったファイルを用いることから、先に存在するTDAの上で結合を含むことが求められるならば、getbsresoirdbsresoは便利である(BSCは結合項を格納しているファイルについて用いられるプレフィックスである)。
 
• tbs_2o_HAYDR_SAVE:
 
 それはHaydock法の結果を含んでいるバイナリファイルです。:三角行列の係数と三つのベクトルが繰り返し法で用いられる。もし、収束に達しなかったならば、アルゴリズムを再計算するために、通常それが用いられる(関連する入力variable gethaydock and irdhaydockを見てください)。
 
• tbs_2o_RPA_NLF_MDF and tbs_2o_GW_NLF_MDF
 
 局所場効果を用いないRPAスペクトルがKSエネルギーとGWエネルギーでそれぞれ用いられる(簡略記憶記号: NLFは非局所場の略語です。一方、MDFは巨視的な誘電関数の略語です)。
 
• tbs_2o_EXC_MDF:
 
 励起子効果をともなう巨視的な誘電関数を報告しているフォーマット化されたファイル。このファイルは私たちの計算の最も重要な結果を含んでいることから、そのフォーマットを議論するために幾らかの時間を費やすのは価値あることです。
 
 最初に私たちは計算の基本的なパラメーターを報告しているヘッダーを持つ。
 
# Macroscopic dielectric function obtained with the BS equation.
#  RPA L0 with KS energies and KS wavefunctions     LOCAL FIELD EFFECTS INCLUDED
# RESONANT-ONLY calculation
# Coulomb term constructed with full W(G1,G2)
# Scissor operator energy =  0.8000 [eV]
# Tolerance =  0.0500
# npweps  = 27
# npwwfn  = 283
# nbands  = 8
# loband  = 2
# nkibz   = 64
# nkbz    = 64
# Lorentzian broadening =  0.1500 [eV]
 
 それから、q点のリストは入射光子の方向を与えている。
 
#  List of q-points for the optical limit:
# q =  0.938821, 0.000000, 0.000000, [Reduced coords]
# q =  0.000000, 0.938821, 0.000000, [Reduced coords]
# q =  0.000000, 0.000000, 0.938821, [Reduced coords]
# q =  0.000000, 0.813043, 0.813043, [Reduced coords]
# q =  0.813043, 0.000000, 0.813043, [Reduced coords]
# q =  0.813043, 0.813043, 0.000000, [Reduced coords]
 
 ディフォルトによって、そのコードはq空間での6つの異なった方向を考慮した巨視的な誘電関数を計算する(逆格子の3つの基本ベクトルと3つのカーテシアン軸)。入力variables gw_nqlwlgw_qlwlを用いて特定の方向を明記することが可能である。
 
 それから、異なった方向に対する周波数の関数のような巨視的で誘電関数の実数と虚数部分をもつセクションが来る。
 
# omega [eV]    RE(eps(q=1)) IM(eps(q=1) RE(eps(q=2) ) ...
0.000  1.8026E+01  0.0000E+00  1.7992E+01  0.0000E+00  1.4292E+01  0.0000E+00  1.3993E+01 0.0000E+00  1.7117E+01  0.0000E+00  1.7080E+01  0.0000E+00
  .... .... ...
 
 あなたの好むソフトを用いてデータを視覚化することができます。例えば、
 
$ gnuplot
gnuplot>  p "tbs_2o_EXC_MDF" u 1:3 w l
 
 最初のq点に対する巨視的な誘電関数の虚数部分をプロットする(吸収スペクトル)。あなたは下記で報告されるものに似たグラフィックスを得るでしょう。
 
 これらの結果は収束していないことに注意してください。私たちはこのチュートリアルの次の段落へ収束テストについての議論を延期する。
 
 スペクトルの最も重要な特徴は3.4と4.3 eV近傍にある2つの局在したピークの存在です。これらの特性とBS Kernelによって振る舞われる役割を理解するために、励起子スペクトルと局所場効果を用いずに得られたRPAの結果と比較することは便利です。
 
 gnuplotコマンドのシークエンスを用いて、
 
$ gnuplot
gnuplot>  p   "tbs_2o_EXC_MDF"     u 1:3 w l
gnuplot>  rep "tbs_2o_RPA_NLF_MDF" u 1:3 w l
gnuplot>  rep "tbs_2o_GW_NLF_MDF"  u 1:3 w l
 
 3つの異なった手法で得られた吸収スペクトルをプロットする。最後の結果は下記の図で報告されます。
 
 RPA-KSスペクトルは、よく知られるDFTのバンドギャップ問題のために、実験での光学的なしきい値を過小評価する。最も重要なことだが、最初のピークの強度は過小評価される。局所場効果が計算で正確に含まれるときでも問題は解決されない。
 
 soenergyを用いてシミュレートしたQP補正を用いたRPA-GWの結果は、RPA-KSを上回るような幾らかの顕著な改善を示さない。:そのRPA-GWスペクトルはギャップの開きのために高い周波数の方へとシフトする。しかし、2つのスペクトルの形状は非常に良く似ている。特に、最初のピークの強度はいまだ過小評価されている。
 
 対照的に、BS kernelの含有は最初のピークの大きさに加えて光学的なしきい値の双方において重要な変化を導く。この簡単な分析は私たちにSiの吸収スペクトルの最初のピークがRPAの範囲内で正確に記述されない強い励起子特性を持つことを教えてくれる。
 
  • 連続的な分数での発散を避けるために用いられるローレンチアンブロードニングの値 zcut を変える。それで、適切な値を用いて _BSR と _HAYDR_SAVEファイルからのHaydockアルゴリズムを再スタートさせる。最終的なスペクトルでのブロードニングの主要な効果は何ですか?ブロードニングに依存して収束するために繰り返し回数は幾つ必要ですか?
  • 局所場効果を用いないBS スペクトルを計算するために、bs_exchange_termbs_coulomb_term に対する適切な値を用いてください。局所場効果を用いたときと用いないときで得られた結果を比較してください。
  • そのコードが前回の計算で生成した共鳴ブロックで読み込み、そして、直接的な対角化(再計算するためにirdbsresoを用いる。しかし、共鳴ブロックを有するファイルの名称を変えることを忘れないでください)に基づく手法を用いてスペクトルを計算しするように、入力ファイルを tbs_2.in 変更してください。transition spaceでの遷移数の関数として2つのアルゴリズムによって必要となるCPU時間を比較してください。どちらが最もよいスケールですか?

 
 励起子スペクトルを収束させることは多くの異なったパラメーターの注意深い分析を必要とする:
 
bs_loband
nband
ecutwfn
ecuteps
ngkpt
nshiftk
shiftk
 
 ブリルアンゾーンでのk点x価電子帯のバンド数x遷移空間で含まれる伝導帯のバンド数で二次的に必要となるメモリがスケールすることから、精度と計算効率の良い妥協点を探すことは非常に重要です。
 
 価電子帯のバンド数と遷移空間に含まれる伝導帯のバンド数で重要な効果を持つ周波数範囲の選択を最初に行う。観測での周波数範囲に近いエネルギーでの遷移のみが寄与すると期待されるから、光学スペクトルはGW補正よりもバンド数でより早く収束することが期待される。
 
 ecutwfnはoscillator行列要素の精度に影響するのみであるから、ecutwfnは通常二次的な役割を振る舞う。私たちはKSSファイルを生成するために用いられるecut値よりも僅かに大きい値へecutwfn を設定することによって初期の基底系の打ち切りを回避することを提案する。この種の問題はプロセッサの数を増やしたり、代替的に、gwmemの適切な選定をしたりすることで通常解決されるのだけれども、メモリ問題を経験したとき、それは初期の平面波の基底系を打ち切る。
 
 既に説明したように:光学スペクトルはブリルアンゾーンサンプリングでゆっくりと収束する。k点数での収束は、それ故、最も重要であり、収束の研究での退屈な部分である。この理由から、本研究では既に見つかっている他のパラメーターに対して収束した値をもう一度用いる。

 ※ 残りの記述は収束性に関する記述であるため省略する。図を見てどの程度で収束しているかを知っておくとよい。k点数を決めるときには、zcut の値を実験値に近づけてから、k点を小さい値から変えていくと良いだろう。
----------------------------------------------------------------------------------------------------------------------------
The lesson Response-Function 1 (RF1)  はABINITでの基本的な応答関数を提供する。与えられた例はAlAs(絶縁体)の動的で誘電的な特性の研究である: ガンマ点でのフォノン、誘電率、ボルン有効電荷、LO−TO スプリッティング、ブリルアンゾーン内でのフォノン。DDBの作成が示される。


1. AlAsの基底状態での構造
 2つのタイプの原子を用いるために、次の入力ファイルvariablesを見てみましょう。
 
ntypat
typat
 
 tolvrsの値はstringentである。本実行で決定される波動関数が応答関数計算の開始点として後に用いられる。しかしながら、このサンプルファイルでのステップ数nstepは15に設定される。そして、目標となるtolvrsに届くためにはこれが十分でないことが分かる。プロダクションランにおいて, 目標であるtolvrsに届くのに十分である大きなnstepの値を選択する。本チュートリアルで、自動テストに関連して軽便であることを懸案事項としているために、大きなnstep値を許可しない。幾つかのチュートリアルサンプルでのこの小さな問題はa side note to the answer to question 1 of lesson 1において簡単に指摘される。それで、チュートリアルでの全てのサンプルは盲目的に追従しないでください。あなた自身であなたの計算をチェックしてください。
 
 あなたは固定したecut (=3Ha)とkptrlatt (8x8x8 Monkhorst-Pack grid)で定義されるk点グリッドで計算を行いましょう。双方の収束パラメーターに関連する収束テストを実寿命で行うことを暗示する。これらの選択の精度とfifth section of this tutorialの終わりに擬ポテンシャルの選択の議論を延期する。それらは受け入れられる結果を与える。よい精度ではないけれど、しかし、より重要で、その速度はチュートリアルにおいて合理的である。
 
 あなたが計算を走らせる(数秒)と、そして、次のエネルギー値を最後のechoセクションで得る。
 
etotal   -9.7626837450E+00
 
 しかしながら、この全エネルギーのより精度ある値で後に返信される。この最後のechoの前の1ダースほどの行で見つけられる。
>>>>>>>>> Etotal= -9.76268374500289E+00
 
 出力ファイルもまた双方の原子の力が消えるのを指摘する。
 
 あなたがちょうど作成した実行は基底状態の配置を定義するために考察される。その上では応答関数が計算される。この基底状態の実行での主要な出力は波動関数ファイルtrf1_1o_WFKです。それは既にあなたがtrf1_1i_WFKとして名称を変えています。実際、次の実行で入力ファイルとしてそれが用いられる。
(省略)

2. 全エネルギーの二次微分のフローズンフォノン計算
 私たちは異なった手法で原子の変位に関する全エネルギーの二次微分を計算することに目標を定める。その目的に対して、あなたは最初にrespfn_helpファイルのsections 0 and the first paragraph of section 1を最初に読み込んでください(省略)。
 
 私たちは後にrespfn_helpファイルのセクション1で紹介されている異なった入力パラメーターの記述をより詳細に説明する。
 
 いまのところ、応答関数理論の結果と直接比較を実行可能にするために、私たちはreducedの座標の最初の軸に沿って、Al原子の変位を摂動として選択する。
 
 Work_rf1に~abinit/tests/tutorespfn/Input/trf1_2.inにファイルをコピーする。これはあなたの入力ファイルです。あなたはそれを記述し、簡単にファイルに関する二つの変化を見る。xredの変化、そして、平面ファイルを読み込み、入力variable irdwfkを用いる。それから、それを実行する。基底状態の構造に関して対称性は低くなる。その結果、k点数は大きくそして、もちろん、CPU時間も増加する。
 
 この実行から、全エネルギー値とreduced座標の変化に関連する全エネルギー値の勾配を得ることが可能になる。
 
rms dE/dt=  3.5517E-03; max dE/dt=  5.0079E-03; dE/dt below (all hartree)
    1       0.005007930232      0.002526304574      0.002526304574
    2      -0.005007868293     -0.002526274885     -0.002526274885
 ...
    >>>>>>>>> Etotal= -9.76268124105780E+00
 
 最初の軸に沿ってAl原子のreduced座標の変化はさらに小さい(1/1000)。そして、私たちは有限の差分式の助けを借りてreduced座標に関連する全エネルギーの二次微分の評価を行う。
 
 私たちは全エネルギー差から最初に始める。全エネルギーはこの摂動に関して対称的である。その結果、それは線形項を持たない。先に実行された基底状態の値と本研究での摂動値の差は、それ故、(1/1000)倍の全エネルギーの二次微分(2DTE)の平方の1/2である。これらの数から、2DTEは5.00791024hartreeである。
 
 代替的に、私たちはreduced勾配から開始する。最初のreduced軸に沿ったAl原子の変位に関連するReduced勾配の値は0.005007930232(hartree)である。first orderで、この量はreduced座標の差による2DTEの生成である。2DTEの評価は、それ故、5.007930232hartreeです。他の評価との一致はさらに良い。
 
 しかしながら、
 
3.
4.
5.
6.
----------------------------------------------------------------------------------------------------------------------------

The lesson Response-Function 2 (RF2)  は前のレッスン RF1で紹介したDDBの分析法を提供する。

 レッスン RF1 によって得られた追加情報は、単純な”Sum-Over-State”近似で、 The lesson on Optic 、周波数依存の線形な光学的誘電関数と二次の非線形光感受率を得るユーティリティへの扉を開く。


----------------------------------------------------------------------------------------------------------------------------
■ GW calculation type: gwcalctyp

0-9: 1-shot
10-19: self-consistent (energy only)
20-29: self-consistent (energy and wavefunction)

0: standard Plasmon-Pole model GW calculation  (1-shot)
2: GW calculation using numerical integration (contour deformation method) (1-shot)

10: standard Plasmon-Pole model GW calculation (self-consistent (energy only))
12: GW calculation using numerical integration (contour deformation method) (self-consistent (energy only))

20: standard Plasmon-Pole model GW calculation (self-consistent)
22: GW calculation using numerical integration (contour deformation method) (self-consistent)

HSE06: 105 (1-shot), 125 ( self-consistent)
PBE0: 205 (1-shot), 225 ( self-consistent )
B3LYP: 305 (1-shot), 325 ( self-consistent )
----------------------------------------------------------------------------------------------------------------------------


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