AkaiKKR

 ここでは全電子計算法であるAkaiKKRを用いて状態密度分布(DOS)やバンド図を描く方法などを解説する。※ AkaiKKR = Machikaneyama である。
------------------------------------------------------------------------------
■ Akai-KKRの使い方(学習方法)
 使用方法は『計算機マテリアルデザイン入門』、大阪大学出版会 によく書かれている。これで自習することも可能。CMDのビギナーズコースでは質問をしない限り書籍と同じ内容を学習することなるのでがっかりしないで欲しい。下記の動画や記述を見て十分であれば、赤井先生や小倉先生、佐藤先生に質問するつもりでなければ、Akai-KKRのためにCMDに参加する価値は大きく減少するだろう。
・Youtubeの動画: https://www.youtube.com/playlist?list=PL9r3BfacS1mG_kipTNqiz_1VwWaEkzHpF
・「ねがてぃぶろぐ(http://gomisai.blog75.fc2.com/blog-entry-600.html)」は素晴らしい。
------------------------------------------------------------------------------
■ Akai-KKRとSPR-KKRの違い(27/Sep/2016)
◇ オススメのコードは?
・最初に扱うのはAkai-KKRの方が良いです。なぜならば、コンパイルし易いのと、計算での収束性が良いこと、バンド分散の計算も安定していることが挙げられます。
・マニアックな計算をしたいという方はSPR-KKRです。フルポテンシャル法が使えるので、水素原子が入るような隙間が沢山ある結晶構造の計算に適しています。
(Akai-KKRもそうですが、マフィンティンポテンシャルを使う場合は、水素原子が入りそうな隙間には原子番号0の原子を指定して計算の精度を高めることができます。GaAsやZrNiSnなどが代表例)
※ 注意点として、SPR-KKRのはマニュアルを読んだだけでは分からないことが多くあります。セミナーの内容を調べたり、開発者に相談したりする苦労も必要になります。

◇ Akai-KKRとSPR-KKRが出来ること
・CPA計算が可能(単位胞で元素置換した系の計算が可能)
・SCFとDOSの計算が可能
・全エネルギーやスピンモーメントを計算可能
 得られたデータから相の安定性や形成エネルギーを算出可能

◇ Akai-KKRが得意
・収束が安定している
・バンド分散の計算が安定している

◇ SPR-KKRが可能なこと
・非等価な原子位置(Akai-KKRでのnatmに対応)を十数以上入力可能
(Fe2VAlのnatmは4なのでAkai-KKRでも十分余裕がある)
・WIEN2kのように空間群から入力ファイルを作れる
・フルポテンシャル計算が可能
・LDA+U, GGA+Uが計算可能
・数多くの物性が計算可能(XAS, XPS(core-level and valence band), ARPES, XES, MXO, AES)
 Residual Resistivity of Alloys, Spin- and Orbital Momentsなど
 http://ebert.cup.uni-muenchen.de/index.php?option=com_content&view=article&id=8&catid=4&Itemid=7
 や
 http://ebert.cup.uni-muenchen.de/index.php?option=com_remository&Itemid=20&func=startdown&id=51&lang=en
 のマニュアルを見ると良い。
・クラスターや界面の計算が可能
------------------------------------------------------------------------------
■ k点数はどう考えればよい?
 「ねがてぃぶろぐ」さんの[AkaiKK使い方]方が私よりも精力的に調べてくださっているので見てください。
「ねがてぃぶろぐ」さんは本当にありがたいです。Killedの対処法などディープなところなどもいつか調べて頂けたら嬉しいです。
 bzqltyについて: 全エネルギーにおいて1 meV/atom程度で十分と聞くので、bzqlty 4だと単金属の計算では1桁足りないかもしれませんね。形成エネルギーは実験値を0.1 eV程度内で予測すると言われるので、0.01 eVまで計算できれば十分かなと思っています。格子定数が長くなるとk点の分割数は少なくて済みます。例えばFe2VAlは5.67 Åなので、Feが2.87 Åとすると、Fe2VAlでは4*5.67/2.87=7.9からbzqlty 8と同精度の計算になると思います(結晶構造による対称性が考慮されて、special k point{Monkhorst-Pack method}となるようにしていればもっとk点数は少なくて済みます)。とは言え、単金属の計算では精度が不安になります。Feをbzqlty 4で計算して、マーナハンの状態方程式にあてはめてみました。もし計算誤差が大きくて曲線上に乗らなければ困ってしまいますが、問題はなさそうですね。sdftyp pbeで計算しているので格子定数は長めになり易いと思います。実験値との差は(2.881/2.866-1)*100=0.52%で3%内に収まっていますね。
Condition (reltyp nrl) Lattice constant [bohr] Total energy [Ry] Spin moment [μB/atom]
go: bzqlty 4 5.4400 -2522.8160225 2.36864
go: bzqlty 4 > dos: bzqlty 8 5.4400 -2518.1759027 0.37377
go: bzqlty 8 5.4400 -2522.8159776 2.32490
go: bzqlty 16 5.4400 -2522.8160602 2.34245
bzqlty 4と16のスピンモーメントの差: (1-2.36864/2.34245)*100 = -1.11%
十数%も違うわけではないので、この差なら問題ないかなと思います。
------------------------------------------------------------------------------
■ MTポテンシャルを用いた弱点が下記のHP[T1]に詳細に記されている。
・ 少し歪ませたらどうなるか? という問題や異方性の高い物質に弱い。しかし、不純物計算がスーパーセル法を用いずにできるという利点がある。
・ CPAはすべての欠陥の配置を平均化したポテンシャルを採用するため、欠陥の濃度が非常に希薄な系は正しく記述できない。また、力の計算が困難であるため、欠陥が生じたことによる局所的な原子位置の緩和の効果は取り入れられない[T3]。
[T1] http://pmt.sakura.ne.jp/wiki/images/Cmdtxt1.pdf
[T2] http://ann.phys.sci.osaka-u.ac.jp/jp/re_sub1.html
[T3] http://www.murata.co.jp/zaidan/annual/pdf/k01/2004/a31130.pdf
------------------------------------------------------------------------------
AkaiKKRの基本的な使用方法
※ case は彼方が付ける計算する系の簡単な名称(私の場合、計算する系の略称や組成とすることが多い)

■ AkaiKKRの解凍とインストール (Intel fortran case)
1. Akai Groupにメールを送り、AkaiKKRのデータを頂く。
2. tar zxvf cpa2002v*.tar.gz と tarminal 上で入力して解凍する。
  x と y はヴァージョンで異なる。
3. cd cpa2002v*
4. touch source/*.f
5. make
6. ifort source/gpd.f -o gpd
http://act.jst.go.jp/content/h13/scom_net/S01/PageMain.html
※ touch ファイルの履歴を更新するという意味。
※ 古いオブジェクトファイル(.o のあるファイル)が残っている場合の対処方法のもうひとつの方法として、makefileの最後に下記を追加(rmの前はTab)し、make clean と入力した後に、make を入力する。
clean:
    rm -f -r $(objs)
※ ifort compiler options
-openmp: OpenMP* 宣言子に基づいてコンパイラーがマルチスレッド・コードを生成するようにします (-fopenmp と同じです)。
-openmp-stubs: シーケンシャル・モードで OpenMP* プログラムをコンパイルします。OpenMP* 宣言子は無視され、OpenMP* スタブ・ライブラリーがリンクされます。(シーケンシャル)
-pc80: 内部 FPU 精度を 64 ビットの仮数に設定します。(デフォルト)
-mp1: 浮動小数点の精度を上げます (速度に与える影響は -mp よりも少ないです)。
-mieee-fp: -mp と同じです。-mno-ieee-fp で無効にできます。
-i-dynamic: インテルが提供するライブラリーを動的にリンクします。-shared-intel を使用してください。
-mcmodel=medium: コードを最初の 2GB までに制限するようコンパイラーに指示します。データは制限されません。
-O2: 処理速度を最大限に最適化します。(デフォルト)
-[no-]ip  ファイル内の単一ファイル IPO を有効 (デフォルト)/無効にします。

■ AkaiKKRの解凍とインストール (gfortran case)
1. Akai Groupにメールを送り、AkaiKKRのデータを頂く。
2. tar zxvf cpa2002v*.tar.gz と tarminal 上で入力して解凍する。
  x と y はヴァージョンで異なる。
3. cd cpa2002v*
4. touch source/*.f
5. gedit makefile
  omp = -fopenmp
  nomp = -openmp-stubs 
  fort = gfortran
  flag = -mpc80 -mieee-fp -m64 -O2 -march=native -funroll-loops
6. make
7. gfortran source/gpd.f -o gpd
http://act.jst.go.jp/content/h13/scom_net/S01/PageMain.html
※ touch ファイルの履歴を更新するという意味。
※ 古いオブジェクトファイル(.o のあるファイル)が残っている場合の対処方法のもうひとつの方法として、makefileの最後に下記を追加(rmの前はTab)し、make clean と入力した後に、make を入力する。
※ 上手くいかない場合は nomp =  も試してみるとよい。-dy オプションはエラーが出ることもある。
clean:
    rm -f -r $(objs)
※ gfortran compiler options
-mpc: Set 80387 floating-point precision (-mpc32, -mpc64, -mpc80)
-mieee-fp: 浮動小数点比較に IEEE 規格を使う
-Bdynamic, -dy, -call_shared: Link against shared libraries
-m64: Generate 64bit x86-64 code
-mavx: Support MMX, SSE, SSE2, SSE3, SSSE3, SSE4.1, SSE4.2 and AVX built-in functions and code generation
※ intel と似た機能のオプションは、gfortran ではmが最初に付いたものになっていることがある。

■ intel : gfortran options
-xHost : -march=native
-ipo : -flto
-no-prec-div : -ffast-math
-mp (=-mieee-fp) : -mieee-fp
-assume noprotect_parens : -fno-protect-parens
(default) : -fstack-arrays
-fast (=-xHOST -O3 -ipo -no-prec-div -static) : -march=native -O3 -flto -ffast-math  -static
-funroll-loops: -funroll-loops
-[no-]unroll-aggressive : -funroll-all-loops
[1] http://thatcadguy.blogspot.jp/2011/03/gnu-fortran-gfortran-vs-intel-fortran.html 
[2] http://www.hpc-sol.co.jp/support/intel_compiler_recipes/20121212_optimization.html 
[3] http://www.softek.co.jp/SPG/Pgi/TIPS/option_for_PDF.html 
[4] http://www.cqpub.co.jp/interface/column/freesoft/2002/200208/08-4.htm
[5] http://www.isus.jp/article/compileroptimization/intel-compilers-101/ 

■ SCF 計算
□ 基本的な系の計算
※ 詳細は http://kkr.phys.sci.osaka-u.ac.jp/pdf/machikaneyama2000_memo.pdf を読むとよい。
1. inにある入力ファイルの例やこのHPの下の方に書かれている入力ファイルを参考にして、計算したい系の入力ファイルを作成する。
2. 書き換えるところは以下の通り。※1
  ・ go の横にあるポテンシャルファイルの名称を data/case にする。
  ・ brvtyp で 計算する結晶構造に対応した記号を入力する。
  ・ a や c/a などは構造を表現できる最低限でよい。aはbohr単位。
    0を入力すると、MJW(qvolum.f)での全ての原子の体積(a.u.)にconcの値を掛けて総和し、それに(1/3)乗した結果を用いる(chklat.fで計算される)。
  ・ ewidth は1.0〜1.5程度。上手く計算出来なければ値を増減してみる。
  ・ reltyp は 非相対論 nrl や スカラー相対論 (すなわちスピン軌道相互作用を落とした相対論的計算)sra で良い。迷ったらsraにする。MCDの場合 srals に(z方向を量子化軸{スピンをz軸方向に揃える}とした相対論計算{スピン軌道相互作用})。
  ・ bzqlty は k点に相当するもので 4 程度でよい。
  ・ ntyp が 原子の種類
  ・ type の行は ntyp で記述した行だけ記述する(不純物やドープでは下記のように記述量が増える)。任意の名称で良い。
  ・ mxl がマフィンティン球内で考慮する l の最大値。
  ・ anclr が原子番号
  ・ conc が 濃度または占有率
  ・ natom が 座標を記述する行数
  ・ atomicx で a, b, c軸に対する原子座標を記述する。通常はa軸の長さを基本単位とする。各ベクトルの長さに対応させたければ、0.5a 0.5b 0.5c などと記述する。
  ・ atomtyp は type で記述した名称と同じにする。
3. 書き換えたらSAVEする。
4. export OMP_NUM_THREADS=CPU数
5. specx < in/case
6. err= が -6.0 程度になるまで specx < in/case を何度も入力して繰り返す。
7. total energy= が表示されるのでメモっておこう。
※1 入力ファイルはbrvtyp の記述が優先され、c/aなどはAkaiKKRコード側で最適な値に変えてくれる。
※2 errは簡単な概念として err ≒ log√[Σ(Vin-Vout)^2] と考えると良い。
※最初にcのある記述はコメントとして無視される。
※ sdftyp では mjwasa というのもある。ASAはMT球よりも球を大きくしたもので、ASAでの球は計算する単位胞の体積と一致するように設定される。MT球ではポテンシャルがフラットな部分が現れるが、ASAではそれが緩和され、よい結果を与える場合がある。

□ 不純物を入れた場合(concを0にする)の計算
1. 基本的な系の計算で作成した入力ファイルで、ncmp を不純物を入れる個数+1の値にする。
2. ncmp の数だけ直ぐ次の行に anclr (不純物の原子番号)と conc(濃度) を書き入れる。
3. Fe に Ir を不純物(0%)として入れた場合の例を下記に示す。
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Fe          2        1       0.0     2      26   100.0
                                                           77        0.0
c------------------------------------------------------------

□ 他の元素をドープをした場合(concを0にしない)の計算
1. 基本的な系の計算で作成した入力ファイルで、ncmp を不純物を入れる個数+1の値にする。
2. ncmp の数だけ直ぐ次の行に anclr (不純物の原子番号)と conc(濃度) を書き入れる。
3. Fe に Ir を1%入れた場合の例を下記に示す。
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Fe          2        1        0.0     2     26    99.0
                                                            77      1.0
c------------------------------------------------------------

□ CPAでの計算
G~ = g0[1-t~g0]-1
GA=G~[1-(tA-t~)G~]-1
GB=G~[1-(tB-t~)G~]-1
cGA+(1-c)GB=G~
※ t~ は有効媒質(メディウム)と呼ばれる。

□ LMDでの計算
0. 強磁性でのSCF計算をする(通常のSCF計算でよい)。
1. cd util
2. case.fmg というファイルを作成して下記のように記述する。
  下記は4成分有る場合の例で、1と3のスピンを反対にしている。元の成分の数だけ数字を記述する。
  ../data/case 1 2 3 4
  ../data/case_lmd 1 -1 2 3 -1 4
  最初のものが通常のSCF計算で得たポテンシャルの情報で、二行目ではスピンを反対にしたい成分を記述した数値の直ぐ後に-1を記述する。この例では、1と3の原子のスピンを反対にする。同じファイル名に書き換えても良い。
3. ./fmg case.fmg
4. cd ..
5. cp in/case in/case_in_lmd
6. 入力ファイルとなる case_in_lmd で ポテンシャルファイル data/case の指定を先ほど作成したdata/case_lmd に下記のように 書き換える。
c----------------------case----------------------------------
      go   data/case
c------------------------------------------------------------
を下記のようにする。
c----------------------case----------------------------------
      go   data/case_lmd
c------------------------------------------------------------
7. specx < in/case_lmd > out/case_out_lmd
※ cd util とした後に、gfortran -o fmg fmg.f などと入力して fmg の実行ファイルを作成する。これは一度だけ行えばよい。
※ cat data/case.info を見ると良い。データは、格子定数、全エネルギー、全スピンモーメントの順で記述される。

■ DOS 計算
0. SCF計算をする。 
1. SCFで作成した入力ファイルをコピーする。
  例 cp in/case in/case_in_dos
2. 入力ファイル case_in_dos での calctyp を go から dos にする。そして、bzqlty を 8〜10 など、SCF計算のときよりも大きくする。
3. export OMP_NUM_THREADS=1
4. specx < in/case_in_dos> out/case_out_dos
5. case_out_dos に結果が出力される。spin upとdownを考慮した計算をした場合には最初にup、次にdownのデータが表示される。
6. case_out_dos ではcase_in_dosで記述した各成分順に下記のように出力される。※
  エネルギー、s軌道、p軌道、d軌道......
7. ./gpd out/case_out_dos と入力するとTotal DOSが表示される。
※ px, py, pz など細かな成分に分けるには、コードを書き換えれば可能になる。
  ドープした系の計算をすると、色々な散乱が多くなるので、DOSで鋭いピークが見られなっていても間違いではない。
※ upとdownのDOSが経験的な結果と逆の場合は、もう一度SCFを計算すると upとdownの入れ替えてくれる。プロットしたときに up→down、down→up と記述し直ても良い。
※ WIEN2kのTotal DOSに形状を合わせると、Akai-KKRではGGA91 を用いてmsex=1001, edelt=0.0036で比較的良く一致するかもしれない(gpd.f における msex も変えてmakeしないと./gpdは動作しないので注意)。

■ バンド分散
0. SCF計算をする。
1. cp in/case in/case_spc
2. case_spc で calctp をgoから spc に書き換える。
3. case_spc の最後にk点を書き加えてSAVEする。
  k点入力例は参考文献[1, 6] を参照すると良い。
4. specx < in/case_spc
5. cd data
6. case_spc_up.spc や case_spc_dn.spc など .spc が付いた名称で結果が記述される(入力ファイルに書いた data/の後の文字に _up.spc や _dn.spc と付く)
  データは下記のように記述される。A(E,k) はブロッホスペクトル関数。
  k点、エネルギー、A(E,K)
7. A(E,k) は系の散乱の程度により、幾らかの広がりを持っている。
  これらのデータをプロットソフトなどでプロットすればよい。
・gnuplotを用いた簡便なプロット方法
1. gnuplot
2. set pm3d map
3. splot "case_spc_up.spc"
・さらに詳細なgnuplotの使い方
1. gnuplot
2. set xrange[0.0:1.1]
3. set yrange[-0.3:0.3]
4. set ylabel "Energy / eV"
5. set title "case"
6. set pm3d map
7. splot "case_spc_up.spc"
モノクロにしたい場合:set palette gray
スケールを元に戻したい場合:set autoscale
コマンドを入力した後に再描画したい場合:replot
簡単な数値計算をしたい場合: 1: ($2*13.6) や ($2 * $3)など、$を付けて()内に式を書くと、計算した値を用いることができる。
終了する場合:q
gnuplot での参考文献:
http://www.cse.kyoto-su.ac.jp/~oomoto/lecture/program/gnuplot/gnuplot.html
http://www-tlab.math.ryukoku.ac.jp/wiki/pukiwiki.php?gnuplot 
http://www.cn.kagawa-nct.ac.jp/~kusama/standard/linux/li_gnuplot_tgif.pdf 
http://kkr.phys.sci.osaka-u.ac.jp/bbs/view_msg.cgi?id=66
※各成分での分散を描きたい場合は、bzqltyを0にして、bzmesh.f での50番の行に関連する項目を書き換える。そうすることで、vkp(r,1)で指定したk点でのPDOSが得られる。

■ K-edge XMCD 計算[6]
1. calctyp を go で、reltyp をsralsにして、input_go などの名称でSAVE
2. specx < in/input_go
3. calctyp の go を mcd に書き換えて、input_mcd などの名称で SAVE
4. source/cemesr.f での data ref/0.75d0/ を data ref/0.0d0/ に書き換えて SAVEし、 make と入力してコンパイルする。※1
5. specx < in/input_mcd > out/output_mcd
6. 得られたoutput_mcd の中で xmd of Component という行が出てくる。各Componentは入力したinput_mcd での順に出力される。
   各列は、エネルギー、吸収、吸収(up)、吸収(down) の順番で出力される。
7. 必要なComponentのデータをコピー&ペーストして、xmd_inputという名称でSAVE
8. xmdcnv を Akai Group から頂いている場合は xmdcnv と Terminal で入力して、最初に xmd_input 、次に xmd_output と入力する。※2
9. gnuplot と入力して、gnuplot を起動させる。
10. plot "xmd_output" と入力すると結果が出力される。
11. 終了させる場合はウインドウを閉じて、Terminal 上で quit と入力する。
12. cemesr.f の data ref/0.0d0/ を data ref/0.75d0/ に戻してSAVEし、make と入力してコンパイルする。
※ 1. ref は 内殻から非占有準位までの範囲を指定するewidth の 内殻からEFまでの比率を指定する。refがディフォルトの0.75だと、内殻からEFまでは ewidth * ref となる。ref を0 にすると EFから非占有準位までの範囲がewidthになる。
※ 2. gfortran -o xmdcnv xmdcnv.f と Terminal 上で入力すると xmdcnv という名称の実行ファイルが作成される。xmdcnv は得られた結果をローレンチアンでブロードニングする。このプログラムは自作しても良い。

■ 電荷密度分布
  概念上可能であるが「電荷密度分布」を出力させるのは green というプログラムの部分を書き換える必要があり難しい。マフィンティンポテンシャルを用いているので、フルポテンシャルではないという点には注意が必要となる。
  このHPの管理者の個人的な考えですが、Valenceの限られた範囲であれば、バンド図で描くk点の数を大きくして、valenceから非占有準位までのブロッホスペクトル関数をプロットしていけば描画ができるかもしれない。EF近傍や非占有準位を描けばMEMとの対比や電気伝導性やSTMを議論できるかもしれませんね。[7]

参考文献
[1] 『計算機マテリアルデザイン入門』、大阪大学出版会
[2] http://pmt.sakura.ne.jp/wiki/index.php?title=Main_Page での Doument for AkaiKKR
[3] http://kkr.phys.sci.osaka-u.ac.jp/jp/document.cgi 
[4] http://www.dyn.ap.eng.osaka-u.ac.jp/CMD/CMD4/CMD4_TEXT/CMD4_MACHIKANEYAMA2000_1.pdf 
[5] http://olymp.cup.uni-muenchen.de/index.php?option=com_remository&Itemid=20&func=startdown&id=51&lang=en 
[6] http://www.ghfecher.de/Run_AKAI-KKR.pdf
[7] http://kkr.phys.sci.osaka-u.ac.jp/bbs/thread.cgi?id=256
[8] http://gomisai.blog75.fc2.com/?tag=AkaiKKR#container 
------------------------------------------------------------------------------
■ 格納されているデータの説明 
A. ls と入力すると下の様なディレクトリ名が表示される(この他にも表示される場合がある)。
   data in makefile out source util
B. これらの意味は下記の通り。
   data: ポテンシャルファイルやバンド図を描くためのcase_up.spcやcase_dn.spc などのデータがある。
   in: 入力ファイルを入れるディレクトリ
   makefile: sourceにあるコードをコンパイルするためのファイル
   out: 出力ファイルを入れるディレクトリ
   source: ソースファイル。エラーが出る場合にはコードを書き換えて対処するのに必要。
   util: 反強磁性用のポテンシャルファイルを作るのに利用する。
------------------------------------------------------------------------------
■ 入力ファイルの例
・ Fe (BCC)
c----------------------Fe------------------------------------
     go   data/fe
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     bcc      5.27  ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.0         nrl        mjw       mag        2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      10        50     0.023
c------------------------------------------------------------
c    ntyp
      1
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Fe          1        1        0.0      2      26    100
c------------------------------------------------------------
c   natm
     1
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Fe
c------------------------------------------------------------

・ Fe(BCC、不純物 Irがある場合)
c----------------------Fe------------------------------------
     go   data/fe_ir
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     bcc      5.27  ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.0          sra      mjw        mag        2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4           50     0.023
c------------------------------------------------------------
c    ntyp
      1
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
      Fe           2        1       0.0       2      26    100
                                                              77      0
c------------------------------------------------------------
c   natm
     1
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0       Fe
c------------------------------------------------------------

・ Ni (FCC)
c--------------------Ni--------------------------------------
     go    data/Ni
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     fcc      6.6594  ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
    0.001     1.0       sra      gga91    mag      2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4       50    0.023
c------------------------------------------------------------
c    ntyp
      1
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
     Ni      1       1      0.0     2
                                          28     100
c------------------------------------------------------------
c   natm
     1
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Ni
c------------------------------------------------------------

・ Co (HCP)
atmicxでは1/3a などと記述しても良い。
c----------------------Co------------------------------------
     go   data/co
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     hcp      4.74 , 1.6215 ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.0        nrl         mjw        mag       2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4           50     0.023
c------------------------------------------------------------
c    ntyp
      1
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Co          1        1       0.0       2      27    100
c------------------------------------------------------------
c   natm
     2
c------------------------------------------------------------
c   atmicx                        atmtyp
     0a         0b         0c       Co
     1/3a       2/3b       1/2c     Co
c------------------------------------------------------------

・ CrAs
CrAsは岩塩構造でFCCと FCCの原点を 1/4 1/4 1/4 にした構造になる。しかし、1/2 1/2 1/2 や 3/4 3/4 3/4 を原点としたFCCの部分に大きな隙間がある。KKR法ではこのような隙間にVacancyを入れると精度が良くなると言われている(下記の例では Vc1とVc2が該当)。
c----------------------CrAs----------------------------------
     go   data/cras
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     fcc     10.68 ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.5          sra       mjw       mag        2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4          50     0.015
c------------------------------------------------------------
c    ntyp
      4
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Cr          1        1        0.0      2      24    100
       As          1        1        0.0      2      33    100
       Vc1        1        1        0.0      0       0    100
       Vc2        1        1        0.0      0       0    100
c------------------------------------------------------------
c   natm
     4
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Cr
     0.25       0.25       0.25     As
     0.5        0.5        0.5      Vc1
     0.75       0.75       0.75     Vc2
c------------------------------------------------------------

・ ZrNiSn
c----------------------XYZ--------------------------------
     go      data/zrnisn
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     fcc     11.5462 ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.5         sra        mjw       mag        2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4          50     0.015
c------------------------------------------------------------
c    ntyp
      4
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Sn         1         1       0.0       2     50    100
       Ni          1         1       0.0       2     28    100
       Zr          1         1       0.0       2      40    100
       Vac        1         1       0.0       2        0    100
c------------------------------------------------------------
c   natm
     4
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Ni
     0.25       0.25       0.25     Sn
     0.5        0.5        0.5      Vac
     0.75       0.75       0.75     Zr
c------------------------------------------------------------

・ ZrNi0953Ni0085Sn
c----------------------XYZ--------------------------------
     go      data/ZrNi0953Ni0085Sn
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     fcc     11.8890 ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
    0.001     1.5       sra      pbe      mag        2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4       200     0.015
c------------------------------------------------------------
c    ntyp
      4
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
     Sn       1      1      0.0     2     50    100
     Ni       2      1      0.0     2     28     95.3
                                           0      4.7
     Zr       1      1      0.0     2     40    100
     Vac      2      1      0.0     2      0     91.5
                                          28      8.5
c------------------------------------------------------------
c   natm
     4
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Ni
     0.25       0.25       0.25     Sn
     0.5        0.5        0.5      Vac
     0.75       0.75       0.75     Zr
c------------------------------------------------------------

・ Fe2VAl
c----------------------X2YZ-------------------------------
     go    data/fe2val
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     fcc     10.8962 ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.5         sra        mjw        mag       2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4          50     0.015
c------------------------------------------------------------
c    ntyp
      3
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Fe          1        1        0.0       2     26    100
       V            1        1        0.0       2     23    100
       Al           1        1        0.0       2     13    100
c------------------------------------------------------------
c   natm
     4
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Fe
     0.25       0.25       0.25     Al
     0.5        0.5        0.5      Fe
     0.75       0.75       0.75     V
c------------------------------------------------------------

・ Fe2V09Ti01Al
c----------------------X2YZ---------------------------------
     go    data/Fe2V09Ti01Al
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     fcc     10.8962 ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.5         sra        mjw      mag         2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4          500     0.015
c------------------------------------------------------------
c    ntyp
      3
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Fe          1        1       0.0       2     26    100
       V            2        1       0.0       2     23     90
                                                             22     10
       Al           1        1       0.0       2     13    100
c------------------------------------------------------------
c   natm
     4
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Fe
     0.25       0.25       0.25     Al
     0.5        0.5        0.5      Fe
     0.75       0.75       0.75     V
c------------------------------------------------------------

・ TiNiSi(No.62 Pnma)
c---------------------------------- 0 ----------------------------------
    go    data/NiSiTi-pbe
c------------------------------------------------------------
c   brvtyp     a       c/a     b/a   alpha   beta   gamma
     so     11.618,   1.141,  0.596,      ,      ,        ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
    0.001      1.2      sra       pbe     mag       2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4       500    0.020
c------------------------------------------------------------
c    ntyp
     3
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
    Ni      1       0      0.0     2    28      100
    Si      1       0      0.0     2    14      100
    Ti      1       0      0.0     2    22      100
c------------------------------------------------------------
c   natm
    12
c------------------------------------------------------------
c   atmicx                        atmtyp
    0.14200000a    0.25000000b    0.43910000c  Ni
    0.85800000a    0.75000000b    0.56090000c  Ni
    0.64200000a    0.25000000b    0.06090000c  Ni
    0.35800000a    0.75000000b    0.93910000c  Ni
    0.76510000a    0.25000000b    0.37710000c  Si
    0.23490000a    0.75000000b    0.62290000c  Si
    0.26510000a    0.25000000b    0.12290000c  Si
    0.73490000a    0.75000000b    0.87710000c  Si
    0.02120000a    0.25000000b    0.81970000c  Ti
    0.97880000a    0.75000000b    0.18030000c  Ti
    0.52120000a    0.25000000b    0.68030000c  Ti
    0.47880000a    0.75000000b    0.31970000c  Ti
c------------------------------------------------------------

・ SrTiO3
c----------------------STO-----------------------------------
     go      data/sto
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     sc      7.38 ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.5         sra        mjw       mag       2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      8          200     0.015
c------------------------------------------------------------
c    ntyp
      3
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
      Sr           1        1        0.0      2      38    100
      Ti            1        1        0.0      2      22    100
      O             1        1        0.0      2        8    100
c------------------------------------------------------------
c   natm
     5
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Sr
     0.5        0.5        0.5      Ti
     0.5        0.5        0.0      O
     0.5        0.0        0.5      O
     0.0        0.5        0.5      O
c------------------------------------------------------------

・Ni1%/Ta2%-STO
※ 「ねがてぃぶろぐ」さんのHPでご指摘いただいて気づきました。Ti 98になっていたので、Ti 97に修正しました。Akai-KKRでは占有率の合計が自動的に100%に調整さるようになっています。修正前の記述ではTiサイトで「Ti 98 : Ni 1 : Ta 2」となっていましたので、そのときのTiサイトの占有率は「Ti 970297... : Ni 0.99 : Ta 1.98」となります。
c----------------------STO-----------------------------------
     go      data/nitasto
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     sc      7.38 ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.5         sra        mjw       mag        2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4          200      0.015
c------------------------------------------------------------
c    ntyp
      3
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Sr          1        1       0.0       2      38    100
       Ti           3        1       0.0       2      22     97
                                                              73      2
                                                              28      1
        O          1        1       0.0       2         8    100
c------------------------------------------------------------
c   natm
     5
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Sr
     0.5        0.5        0.5      Ti
     0.5        0.5        0.0      O
     0.5        0.0        0.5      O
     0.0        0.5        0.5      O
c------------------------------------------------------------

・ LiCoO2
c----------------------LiCoO2--------------------------------
      go  data/licoo2
c------------------------------------------------------------
c   brvtyp     a         c/a   b/a   alpha   beta   gamma
     rhb    9.373   ,      ,    ,   32.97 ,     ,    ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.7         sra       mjwasa     mag       2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update     4           30        0.02
c------------------------------------------------------------
c    ntyp
      3
c------------------------------------------------------------
c   type     ncmp    rmt    field   mxl  anclr   conc
       Li            1        1       0.000   2     3     100
       Co           1        1       0.000   2    27    100
       O             1       1        0.000   2     8     100
c------------------------------------------------------------
c   natm
     4
c------------------------------------------------------------
c   atmicx                   atmtyp
     0           ,  0           ,  0       , Co
     0.5a        ,  0.5b        ,  0.5c    , Li
     0.26a       ,  0.26b       ,  0.26c   , O
     0.74a       ,  0.74b       ,  0.74c   , O
c------------------------------------------------------------

・ Li
c----------------------Li------------------------------------
     go   data/li
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     bcc      6.633  ,      ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
     0.001     1.2         sra      mjwasa    nmag      2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      5          50     0.023
c------------------------------------------------------------
c    ntyp
      1
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
       Li           1        1        0.0      2        3    100
c------------------------------------------------------------
c   natm
     1
c------------------------------------------------------------
c   atmicx                        atmtyp
     0          0          0        Li
c------------------------------------------------------------

・ Zr
c----------------------Co------------------------------------
     go   data/Zr
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     hcp      6.1076, 1.5925 ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
    0.001     1.0       sra      gga91      mag      2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4        50    0.023
c------------------------------------------------------------
c    ntyp
      1
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
    Zr       1       1      0.0     2
                                          40    100
c------------------------------------------------------------
c   natm
     2
c------------------------------------------------------------
c   atmicx                        atmtyp
     0a         0b         0c       Zr
     1/3a       2/3b       1/2c     Zr
c------------------------------------------------------------

・ Sn
c--------------------Sn--------------------------------------
     go    data/Sn
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     st      11.0204,  0.5456,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
    0.001     1.0       sra      gga91    mag      2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update       4       50    0.023
c------------------------------------------------------------
c    ntyp
      1
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
     Sn      1       1      0.0     2
                                          50     100
c------------------------------------------------------------
c   natm
     4
c------------------------------------------------------------
c   atmicx                        atmtyp
     0a        0b        0c         Sn
     1/2a      0b        1/2c       Sn
     1/2a      1/2b      1/2c       Sn
     0a        1/2b      1/2c       Sn
c------------------------------------------------------------
------------------------------------------------------------------------------
■ バンド図を描くための入力のテンプレート(spcでのk点入力)
下記はWIEN2kと同じk点の入力ファイルをAkaiKKR用に書き下したものである。空間群の情報は「http://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-table?from=kv」を参照すると良い。
下記に空間群の説明がある。
http://www.ghfecher.de/Run_AKAI-KKR.pdf
※ 出力データに基本ベクトルとなる ga, gb, gc が出力されている。http://www.ghfecher.de/Run_AKAI-KKR.pdf に書かれている点の値にga, gb, gc を掛けて、kx, ky, kz での値を算出する。
 FCCでX点(1/2, 0, 1/2)の場合には下記のようになる。
(1/2)*ga = (1/2)*(-1, 1, 1) = (-1/2, 1/2, 1/2)
0*gb = 0*(1, -1, 1) = (0, 0, 0)
(1/2)*gc = (1/2)*(1, 1, -1) =  (1/2, 1/2, -1/2)
上記を合計すると、 (-1/2, 1/2, 1/2) +  (0, 0, 0) +  (1/2, 1/2, -1/2) = (0, 1, 0) となる。(調査中)

□ Simple Cubic (R→LAMBDA→GAMMA→DELTA→X→Z→M→SIGMA→GAMMA)
0.500 0.500 0.500
0.475 0.475 0.475
0.450 0.450 0.450
0.425 0.425 0.425
0.400 0.400 0.400
0.375 0.375 0.375
0.350 0.350 0.350
0.325 0.325 0.325
0.300 0.300 0.300
0.275 0.275 0.275
0.250 0.250 0.250
0.225 0.225 0.225
0.200 0.200 0.200
0.175 0.175 0.175
0.150 0.150 0.150
0.125 0.125 0.125
0.100 0.100 0.100
0.075 0.075 0.075
0.050 0.050 0.050
0.025 0.025 0.025
0.000 0.000 0.000
0.050 0.000 0.000
0.100 0.000 0.000
0.150 0.000 0.000
0.200 0.000 0.000
0.250 0.000 0.000
0.300 0.000 0.000
0.350 0.000 0.000
0.400 0.000 0.000
0.450 0.000 0.000
0.500 0.000 0.000
0.500 0.050 0.000
0.500 0.100 0.000
0.500 0.150 0.000
0.500 0.200 0.000
0.500 0.250 0.000
0.500 0.300 0.000
0.500 0.350 0.000
0.500 0.400 0.000
0.500 0.450 0.000
0.500 0.500 0.000
0.450 0.450 0.000
0.400 0.400 0.000
0.350 0.350 0.000
0.300 0.300 0.000
0.250 0.250 0.000
0.200 0.200 0.000
0.150 0.150 0.000
0.100 0.100 0.000
0.050 0.050 0.000
0.000 0.000 0.000

□ FCC (GAMMA→DELTA→X) (http://www.ghfecher.de/Run_AKAI-KKR.pdfで下記のように記載されている)
#- k-vectors for spc -----------------------------------------#
#- fcc Gamma -- Delta -- X -----------------------------------#
0.00000000 0.00000000 0.00000000
0.00000000 0.02500000 0.02500000
0.00000000 0.05000000 0.05000000
0.00000000 0.07500000 0.07500000
0.00000000 0.10000000 0.10000000
0.00000000 0.12500000 0.12500000
0.00000000 0.15000000 0.15000000
0.00000000 0.17500000 0.17500000
0.00000000 0.20000000 0.20000000
0.00000000 0.22500000 0.22500000
0.00000000 0.25000000 0.25000000
0.00000000 0.27500000 0.27500000
0.00000000 0.30000000 0.30000000
0.00000000 0.32500000 0.32500000
0.00000000 0.35000000 0.35000000
0.00000000 0.37500000 0.37500000
0.00000000 0.40000000 0.40000000
0.00000000 0.42500000 0.42500000
0.00000000 0.45000000 0.45000000
0.00000000 0.47500000 0.47500000
0.00000000 0.50000000 0.50000000
#----------------------- end input data -----------------------#

New type input
◇ FCC ( check OK)
--------------------
c--- same line k point  as WIEN2k_FCC.klist_band, Conventional basis
110
c--- W-point
1.000 0.500 0.000
c--- L-point
0.500 0.500 0.500
c--- Gamma point
0.000 0.000 0.000
c--- X-point
1.000 0.000 0.000
c--- W-point
1.000 0.500 0.000
c--- K-point
0.750 0.750 0.000
--------------------

◇ BCC ( check OK)
--------------------
c--- same line k point as WIEN2k_BCC.klist_band, Conventional basis
120
c--- Gamma-point
0.000 0.000 0.000
c--- H-point
0.000 1.000 0.000
c--- N-point
0.500 0.500 0.000
c--- Gamma-point
0.000 0.000 0.000
c--- P-point
0.500 0.500 0.500
--------------------

◇ HCP (No check)
--------------------
c--- same line k point as WIEN2k_HCP.klist_band, Conventional basis
240
c--- Gamma-point
0.000 0.000 0.000
c--- M-point
0.500 0.000 0.000
c--- K-point
0.333 0.333 0.000
c--- Gamma-point
0.000 0.000 0.000
c--- A-point
0.000 0.000 0.500
--------------------

◇ SC (simple cubic and super cell)
--------------------
c--- same line k point as WIEN2k_SC.klist_band, Conventional basis
60
c--- R-point
0.500 0.500 0.500
c--- Gamma-point
0.000 0.000 0.000
c--- X-point
0.500 0.000 0.000
c--- M-point
0.500 0.500 0.000
c--- Gamma-point
0.500 0.000 0.000
--------------------
------------------------------------------------------------------------------
■ 物性の計算方法
・格子定数最適化
  構造最適化機能はないので、格子定数を徐々に変えて、全エネルギーが最も低くなる結果を採用させる。
1. 入力ファイルの最後に、格子定数を変えた入力ファイルをコピー&ペーストして貼り付ける。
2. 上記を格子定数を変えたいだけ行う。
3. specx < in/case_latt > out/case_latt で計算させる。
4. grep "total energy=" out/case_latt > out/case_latt_te とすると全エネルギーが入力ファイルで記述した順にcase_latt_teに記述される。
※ Murnaghan の状態方程式でフィットさせると、置換量と格子定数の変化が見やすくなる。
Igorでは、解析 > 回帰分析 > 新規回帰関数
回帰係数: B, BP, V0, E
独立変数: x
f(x) = (B*x^3)/(BP*(BP-1)) * (BP*(1-V0/x^3)+(V0/x^3)^BP - 1 ) + E
かなり正確に、予測される最適な体積V0とエネルギーEを入力することが必要になる。
Bは放物線での傾き(Bが大きいほど傾きが急)、BPはU型になる線を全体的にどの程度角度を変化(BPが小さいとU型は左回転)させるかに効いてくる。
ほぼ手作業で回帰係数を決めることになることも多々あるが、最適な格子定数から2%内であれば、私の場合は、綺麗に点の上に線がのった。

・ 平均場近似でのキュリー温度Tcの計算
  Tc = ( 2 /3 ) * ( Elmd - Eferro ) / kB
  ここで ElmdはLMD(50%をup、残りの50%をdownにした場合での計算)での全エネルギー/単位胞あたりの磁性に強く関係する原子数、Eferroはupに揃えた場合の全エネルギー/単位胞あたりの磁性に強く関係する原子数、kbはボルツマン定数で 6.3336*10^-6 (Ry/K) である。結果は実験結果よりも数百Kほど大きめになる。 1mRy ≒ 100 K と覚えておくと良い。
※ 全体でのspin momentが0にならないは、全ての原子に対してupとdownをそれぞれ50%ずつにする。

・ 起電力の計算
  -V = Total Energy (LiCoO2) - Total Energy (CoO2) - Total Energy (Li)
結果は 13.6 * Ry = eV となるが、1e あたりにすると、Vとなる。

・形成エネルギー(=生成エンタルピー)
 =Total Energy (化合物ABC) - Total Energy (原料A) - Total Energy (原料B) - Total Energy (原料C)
原子数は揃える。
 形成エネルギー > 0: 化学反応でいう吸熱反応で絶対零度ではこの反応は起きない。有限温度では多少なりとも入る。不純物原子の乱雑な配置によるエントロピー増加が効いた効果である。その濃度は、
 濃度 c = exp ( -1 * 形成エネルギー / (kB * T) )
で与えられる。[TF1] この形成エネルギーの評価は0.1 eVくらいの誤差で予測できる。[TF1]
[TF1] 赤井純久、白井光雲、「密度汎関数法の発展」、シュプリンガー・ジャパン
※ この形成エネルギーおよび濃度による不純物の固溶限の評価に対して、局所的な原子位置の緩和を導入しない場合には、その予測精度が落ちると思われる。構造緩和をしていない場合は、構造緩和していないときよりもエネルギー的に不利である。そのため固溶限が実験値よりも小さく出ることが考えられることを忘れてはいけない。もちろん、有限温度では、エントロピー(配置、そして、フォノン)の影響も忘れてはいけない。
※ ボルツマン定数 kB = 1.3807 * 10-23 J/K = 8.6171 * 10-5 eV/K

・自由エネルギー
 G = Total Energy + pV  - S * T = Total Energy + pV - (-kB * Σi (xi * ln (xi)) ) * T
ここで、xi は濃度で、Σi xi = 1の関係がある。Sは配置のエントロピーとフォノンでのエントロピーになるが、フォノンでのエントロピーは計算が大変なので、最右辺では配置のエントロピー[- (-kB * Σi (xi * ln (xi)) ) * T]のみにしている。
G: http://en.wikipedia.org/wiki/Gibbs_free_energy
S: http://en.wikipedia.org/wiki/Entropy
分子系(振動のエントロピーの計算方法):http://www2.itc.nagoya-u.ac.jp/pub/pdf/pdf/vol06_03/282_296kouza.pdf
※ 蒸気圧の低い母材(金属など)を溶かして試料を作る場合、大部分が気体になって仕事をしているというわけではないので、pVによる仕事を無視してpV=0とする。
http://www.escoltd.co.jp/page_technicalinf/target_blank/meltpoint_win.html

・ hyper fine (NMRでよく見る。Coreを横切る直線(z軸方向で考えると分かりやすい)でのup spin と down spin の差)
  出力ファイルに出力される。

※ te = Total energy [Ry / (atoms in the unit cell)]
 出力ファイルに[atoms in the unit cell]が記述されている。さらに、[atoms in the unit cell]は、入力ファイルで指定した原子位置の記述の部分に対応している(単位胞中の原子数は、原子位置の記述の行数に対応している)。
------------------------------------------------------------------------------
■ formation energy: calculation test
・ sdftyp: mjw, bzqlty: 4, magtyp: mag, reltyp: sra
Al(FCC, a=7.65 bohr):  -483.8260411 Ry/atom
V (BCC, a=5.72 bohr): -1894.7477552 Ry/atom
Fe(BCC, a=5.41 bohr): -2541.0206618 Ry/atom
Fe2VAl(FCC, a=10.8962): -7460.7728034 Ry/4atom
Formation energy = Fe2VAl - Fe*2 - V - Al
= -7460.7728034 - ( -7460.6151199 ) = -0.1576835 Ry/4atom
= -0.1576835 Ry/4atom * 13.6058 = -2.14541 eV/4atom

-----------------------------------
・ sdftyp: gga91, bzqlty: 4, magtyp: mag, reltyp: sra
Al(FCC):  -485.8619040 Ry/atom
V (BCC): -1899.1734716 Ry/atom
Fe(BCC): -2546.2273215 Ry/atom
Fe2VAl(FCC): -7477.6160146 Ry/4atom
Formation energy = Fe2VAl - Fe*2 - V - Al
= -7477.6160146 - ( -7477.4900186 ) = -0.125996 Ry
-0.125996 Ry * 13.6058 = -1.71428 eV

・ sdftyp: gga91, bzqlty: 4, magtyp: mag, reltyp: sra
Fe2V103Al097(FCC): -7520.0129793 Ry
-----
G = Total Energy + pV - S * T = Total Energy + pV - (-kB * Σi (xi * log (xi)) ) * T
kB = 8.6173324(78)×10^−5 eV / K
※ ここで、Σi xi = 1

S = -8.6173324(78) * 10^(-5) * ( [(0.03)*ln(0.03) + (0.97)*ln(0.97)] )
= -8.6173324(78)×10^(-5) *[ -0.134722 ] = 1.161128 * 10^(-5) eV / K
S(300K) = 0.0000116112 eV/K * 300 K = 0.0034834 eV ( 0.0002560 Ry )
S(600K) = 0.0000116112 eV/K * 600 K = 0.0069667 eV ( 0.0005120 Ry )
※ ここでは、エントロピーは4dサイトのAlに過剰なVが一部ランダムに置換するとして計算した。
FCC
Fe 4a 0 0 0
V 4c 1/4 1/4 1/4
Fe 4b 1/2 1/2 1/2
Al 4d 3/4 3/4 3/4
or
Fm-3m
Al 4a 0 0 0
V 4b 1/2 1/2 1/2
Fe 8c 1/4 1/4 1/4
-----
Formation energy (0 K) = Fe2V103Al097 - ( Fe*2 - V*1.03 - Al*0.97 )
= -7520.0129793 -( -5092.454643 + -1956.148675748 + -471.28604688 ) = -7520.0129793 + 7519.889365628
= -0.123613672 Ry/4atom = -1.68186 eV
-----
F = E - ST
G = E + pV - ST = -7520.0129793 + pV - (0.0000116112 / 13.6058)  * T [Ry]

Formation energy (300K) = E - ST - ( Fe*2 - V*1.03 - Al*0.97 )
= -7520.0129793 - 0.0002560 -( -5092.454643 + -1956.148675748 + -471.28604688 ) = -7520.013235 + 7519.889365628
= -0.1238694 Ry/4atom = -1.68534 eV

Formation energy (600K) = -1.68882 eV
-----
・ sdftyp: pbe, bzqlty: 4, magtyp: mag, reltyp: sra
Zr (hcp, a=6.187, c/a=1.5925): -14395.1440238 Ry/2atom = -7197.5720119 Ry/atom
Ni (fcc, a=6.7260): -3041.5743821 Ry/atom
Sn (st, a=11.6938, c/a=0.540144): -49418.3963520 Ry/4atom = -12354.599088 Ry/atom
ZrNiSn (fcc, a=11.7713): -22593.8780373 Ry/3atom
Formation energy = ZrNiSn - (Zr + Ni + Sn)
= -22593.8780373 - ( -7197.5720119 + -3041.5743821 + -12354.599088 ) =  -0.1325 Ry
-0.1325 Ry * 13.6058 =  -1.80 eV

ZrNi0953Ni0085Sn (fcc, a=11.8890): -22700.3046071 Ry
Formation energy = ZrNi0953Ni0085Sn - (Zr + Ni*(0.953+0.085) + Sn)
= -22709.4302936 - ( -7197.5720119 + -3041.5743821*(0.953+0.085) + -12354.599088 ) =  -0.1050 Ry
 -0.1050 Ry * 13.6058 = -1.43  eV

ZrNi2Sn (fcc, a=12.0067): -25635.4903589 Ry/4atom
Formation energy = ZrNi2Sn - (Zr + Ni*2 + Sn)
= -25635.4903589 - ( -7197.5720119 + -3041.5743821*2 + -12354.599088 ) = -0.1597  Ry
-0.1597 Ry * 13.6058 = -2.17 eV

Cu (fcc, a=6.9679): -3309.9567849 Ry/atom
ZrNiSnCu007 (Cu/Vac, a=11.8890) = -22825.5606435 Ry
  Formation energy = ZrNiSnCu007 - (Zr + Ni + Sn + Cu*0.07)
  = -22825.5606435 - ( -7197.5720119 + -3041.5743821 + -12354.599088 +  -3309.9567849*0.07 )
  = -0.1181865 Ry = -1.608023 eV
  S = -8.6173324 * 10^(-5) * ( (0.07)*ln(0.07) + (0.93)*ln(0.93))
  = 0.000021857 eV
  F = E - ST = -1.608023 - 0.000021857*300 (@300K) = -1.614580 eV (@300K)
※ Vacuume site をCuが占有したときのエントロピー
  Vacuume site
Cu: 0.07 0.07
Vacuume: 0.93 0.93
S = -8.6173324 * 10^(-5) * ( (0.07)*ln(0.07) + (0.93)*ln(0.93) = 0.000021857 eV

ZrNiSnCu007 (Cu/Ni, Ni/Vac, a=11.8890) = -22825.5596185 Ry
  Formation energy = ZrNiSnCu007 - (Zr + Ni + Sn + Cu*0.07)
  = -22825.5596185 - ( -7197.5720119 + -3041.5743821 + -12354.599088 +  -3309.9567849*0.07 )
  = -0.117162 Ry = -1.594077 eV
  S = -8.6173324 * 10^(-5) * ( (0.07*0.07)*ln(0.07*0.07) + (0.93*0.07)*ln(0.93*0.07) + (0.07*0.93)*ln(0.07*0.93) + (0.93*0.93)*ln(0.93*0.93) )
  = 0.000043714 eV
  F = E - ST = -1.594077 - 0.000043714*300 (@300K) = -1.60719 eV (@300K)
※ Cu/Ni で残りのNiがVacuumeサイトを占有するときのエントロピー
  Cu: 0.07 Ni: 0.93
Ni: 0.07 0.07 * 0.07 0.07 * 0.93
Vacuume: 0.93 0.93 * 0.07 0.93 * 0.93
S = -8.6173324 * 10^(-5) * ( (0.07*0.07)*ln(0.07*0.07) + (0.93*0.07)*ln(0.93*0.07) + (0.07*0.93)*ln(0.07*0.93) + (0.93*0.93)*ln(0.93*0.93) ) = 0.000043714 eV
エントロピーは乱雑(無秩序)の方が大きくなる。例えば、確率が0.5であると一番大きな値になる。

dE = -1.608023 - -1.594077 = -0.013946 eV = 162 K (1 eV = 11604 K)
dF(@300K) = -1.614580 - -1.60719 = 0.00739 eV
dF(@638K) = 0 eV
-----
・ sdftyp: pbe, bzqlty: 8, magtyp: mag, reltyp: sra
Zr (hcp, a=6.1687, c/a=1.5925): -14395.1432687 Ry/2atom = -7197.57163435 Ry/atom
Ni (fcc, a=6.7260): -3041.5743375 Ry/atom
Sn (st, a=11.8400, c/a=0.5333922): -49418.3970222 Ry/4atom = -12354.59925555 Ry/atom
ZrNiSn (fcc, a=11.7713): -22593.8778914 Ry/3atom
Formation energy = ZrNiSn - (Zr + Ni + Sn)
= -22593.8778914 - ( -7197.57163435 + -3041.5743375 + -12354.59925555 )
= -0.132664 Ry = -1.804998 eV

ZrNi0953Ni0085Sn (fcc, a=11.8890): -22709.4302586 Ry
Formation energy = ZrNi0953Ni0085Sn - (Zr + Ni*(0.953+0.085) + Sn)
= -22709.4302586 - ( -7197.57163435 + -3041.5743375*(0.953+0.085) + -12354.59925555 )
= -0.105206375 Ry = -1.431417  eV

ZrNi2Sn (fcc, a=12.0067): -25635.4903307 Ry/4atom
Formation energy = ZrNi2Sn - (Zr + Ni*2 + Sn)
= -25635.4903307 - ( -7197.57163435 + -3041.5743375*2 + -12354.59925555 )
= -0.1707658  Ry = -2.323405 eV

Cu (fcc, a=6.9679): -3309.9588276 Ry/atom
ZrNiSnCu007 (Cu/Vac, a=11.8890) = -22825.5606594 Ry
  Formation energy = ZrNiSnCu007 - (Zr + Ni + Sn + Cu*0.07)
  = -22825.5606594 - ( -7197.57163435 + -3041.5743375 + -12354.59925555 +  -3309.9588276*0.07 )
  = -0.118314068 Ry = -1.609757546 eV
  S = -8.6173324 * 10^(-5) * ( (0.07)*ln(0.07) + (0.93)*ln(0.93))
  = 0.000021857 eV
  F = E - ST = -1.609757546 - 0.000021857*300 (@300K) = -1.616314646 eV (@300 K)
ZrNiSnCu007 (Cu/Ni, Ni/Vac, a=11.8890) = -22825.5596477 Ry
  Formation energy = ZrNiSnCu007 - (Zr + Ni + Sn + Cu*0.07)
  = -22825.5596477 - ( -7197.57163435 + -3041.5743375 + -12354.59925555 +  -3309.9588276*0.07 )
  = -0.117302368 Ry = -1.595992559 eV
  S = -8.6173324 * 10^(-5) * ( (0.07*0.07)*ln(0.07*0.07) + (0.93*0.07)*ln(0.93*0.07) + (0.07*0.93)*ln(0.07*0.93) + (0.93*0.93)*ln(0.93*0.93) )
  = 0.000043714 eV
  F = E - ST = -1.595992559 - 0.000043714*300 (@300K) = -1.609106759 eV (@300 K)
dE = -1.609757546 - -1.595992559 = -0.013764987 eV = 159.73 K (1 eV = 11604 K)
dF(@300 K) = -1.616314646 - -1.609106759 = 0.007207887 eV
dF(@629.77 K) = 0 eV

平衡状態において、O KでZrNiSnCu007 (Cu/Vac)が100%
有限温度では
ZrNiSnCu007 (Cu/Ni, Ni/Vac)の濃度 c = exp ( -1 * 形成エネルギー / (kB * T) ) or
exp ( -1 * 自由エネルギー / (kB * T) )
= exp ( -1 * dF / (kB * T) )
= EXP(-1*((-1.595992559-0.000043714*T)-(-1.609757546-0.000021857*T))/(8.6173324*10^(-5)*T))
300 K: 0.756 = 0.756*100 [%] = 75.7%
600 K: 0.987 = 0.987*100 [%] = 98.7%
※ 形成エネルギーの評価は0.1 eVくらいの誤差で予測できると教科書に書かれているので参考にしかならないが、理論では”平衡状態”で600 Kにおいて98.7%がZrNiSnCu007 (Cu/Ni, Ni/Vac)である。
※ 他の材料だと、歪取焼鈍(1273K, 1h),規則化焼鈍(673K, 4h)を行なって結晶性の高い試料を得たりする(873 Kで粘性かなり高くなる?)。
※ 形成エネルギーを比較して安定な構造を決めてしまうのはまずい。この例のように有限温度では安定な構造が別のものになることがあるので注意して欲しい。
https://www.jstage.jst.go.jp/article/tetsutohagane1955/78/3/78_3_399/_pdf
http://www.geocities.jp/e_kamasai/shiryou/d-base/material.html#3
-----
TiNiM0.05Sn (M=Co, Cu), 0 K
compsition a Ry moment
Co/Ni & Ni/Vac 11.2189  -22733.1703740 0.00000
Co/Vac 11.2189  -22733.1675953 0.03664
Cu/Ni & Ni/Vac 11.2240 -22759.3209509 0.00000
Cu/Vac 11.2240 -22759.3218961 0.00000
T = (-22759.3209509 - -22759.3218961)*13.6058/(0.000043714 - 0.000021857) = 588.38 K
Stable
Co 4c site
Cu 4d site (=< 588.38 K), Cu 4c site (>588.38 K)

Co (4c site) [*1]
Cu (4d site) [*1]
[*1] R. A. Dowine et al., Chem. Mater., 27 (2015) 2449-2459.: http://pubs.acs.org/doi/abs/10.1021/cm5045682
-----
------------------------------------------------------------------------------
■ エラーの対処法
・配列の数が少ないと表示された場合の対処法
1. error の記述をメモっておく。
2. source/specx.f を開いて、error で記述された配列や引数に近い数値を大きくしてSAVE
3. make と入力すればコンパイルされる。
  よく起こるエラーの例としては、 ncmpmx で 成分の数が多くなって計算が出来ないという場合である。この他、mxlmxなどが該当する。

・spc でエネルギー分解能を高くしたい場合
1. kwrite ./source/specx.f または emacs ./source/specx.f
2. msex= の後の数値を大きくする。
3. SAVEして閉じる。
4. make
※ msex=1001 としても正常に動作していた。

・spc でk点の入力数を増やしたい場合 または nk3 に関するエラーが出た場合
1. kwrite ./source/specx.f または emacs ./source/specx.f
2. nk3x= の後の数値を大きくする。
3. SAVEして閉じる。
4. make
※ nk3x=350 としても正常に動作していた。
------------------------------------------------------------------------------
■ 入力ファイルの詳細
http://kkr.phys.sci.osaka-u.ac.jp/pdf/machikaneyama2000_memo.pdf を読んでおくと良い。
・ edelt: 0.001 のままで良い。0.0001なども可能。pdosを計算するときに用いられる。
・ sdftyp: 交換相関ポテンシャルの指定。ガスモデルのmjwで良い(vbhもvwnもLDAでガスモデル。GGAが良いならgga91とする)。
・ magtyp: 磁性ではmag, 非磁性ではnmag にする。
・ record:
  init: ポテンシャルファイルを一から新規に作る。
  2nd: 最新のポテンシャルデータ。
  1st: 2番目に新しいデータ。
  dataにポテンシャルファイルが無い場合は init になる。
・ outtyp: update だと更新する。quitは更新しない。
・ maxitr: 最大イタレーション回数。
・ pmix: Vin と Vout の混ぜ具合。大きくすると早く収束するが、発散する場合もある。
・ rmt: マフィンティン半径。0を指定したり、非現実的な値だと自動的に調整される。0だと原子半径比、1だと全元素で一定の値。マニュアル(http://kkr.phys.sci.osaka-u.ac.jp/pdf/akaikkr_j.pdf)には「与えられたマフィンティン半径では球どうしが重なってしまう場合、マフィンティン半径は与えられた半径比で球が接するように設定し直される。0を与えた時または省略した時は半径比として原子半径比がとられる」とある。
・ mxl: マフィンティン球内で計算するlの最大値。最低限の値で良い。入力に入れない最大値を超えるlはlocal orbitalで処理。
 散乱に取り入れる各運動量の最大値がmxlであり、他は平面波となる(無限個のカットオフ?)。
・ 最初に#やcと付いた行はコメントとして無視される。空白も無視される。
・ 改行は関係無い。
※ spinの量子化軸は常にz軸方向になっている。
■ 出力ファイルの詳細
  だいたい下記の順に出力される。
・ 入力ファイルの詳細が記述される。
・ 次にメッシュ(計算するエネルギー経路)が記述される。
・ ポテンシャルファイルを作成
・ rmt のチェック
・ 格子定数の修正
・ ポテンシャルファイルの修正
 幾つかの記述
・ 最初の方にあるitrは各成分の原子の計算の繰り返し回数
・ state で数値が記述されているのがCore
・ neu:  -だと電子が足らない、+だと電子が多い、との記述。
・ spin moment: 磁気モーメント。単位は[μB/atom]。
・ orbital moment: スピン軌道相互作用(srals)を入れた場合に右回りと左周りで値が異なるために0の値にならなくなる。
・ intc: マフィンティン球同士の間の隙間の電荷
・ ints: マフィンティン球同士の間の隙間の磁気モーメント
・ rms err: up、downの順に記述
・ ef: up、downの順の記述
・ dep または def: up dos (EF), down dos (EF)
・ valence charge: 計算の妥当性を確認できる。
・ meshr: 原子に対するメッシュ数
・ mse: 入力ファイルでのk点(bzqlty)を入れると自動的に出力される。
・ ng : 内挿のプログラム
※ 電位勾配や残留抵抗(0 Kでの抵抗値)の出力には対応していない?
[1] http://kkr.phys.sci.osaka-u.ac.jp/bbs/thread.cgi?id=56
[2] http://kkr.phys.sci.osaka-u.ac.jp/bbs/view_msg.cgi?id=168
[3] http://kkr.phys.sci.osaka-u.ac.jp/bbs/view_msg.cgi?id=37
[4] http://kkr.phys.sci.osaka-u.ac.jp/bbs/view_msg.cgi?id=258
[5] http://kkr.phys.sci.osaka-u.ac.jp/bbs/view_msg.cgi?id=68
[6] http://kkr.phys.sci.osaka-u.ac.jp/bbs/thread.cgi?id=38
------------------------------------------------------------------------------
■ 交換相関項
  source/potenv.f や source/potenv2.f などに交換相関ポテンシャルの記述がある。
 具体的なポテンシャルへは exczzz.f (zzzはxcポテンシャルの略称)
  コードを解読できていないので申し訳ないが、コードを書き換える能力のある方は、交換相関項を PBE, DFT+UやSIC、PBE0などにして色々とチャレンジして頂けるとありがたい。cpa2002v009c では下記のものが実装されている(多体問題は粒子密度がしっかりと求まっていれば良いことも覚えておくとよい)。
c     +---------------------------------------------------+
c        sdf parametrization by vbh     ... vbh
C                            by mjw     ... mjw
c                            by vwn     ... vwn
c     Langreth, Mehl non-local + mjw    ... lmmjw
c                   non-local + vbh     ... lmvbh  missing
c     Perdew, Yue    non-local + mjw    ... pymjw
c                   non-local + vbh     ... pyvbh  missing
c                   non-local + vwn     ... pyvwn
C     Perdew, Wang   GGA91              ... gga91
c     Engel, Vosko   GGA exchange-only  ... ev
c     +---------------------------------------------------+

sdftyp type of the parametrization of the exchange energy. isdf [XC1]
vbh v. Barth and Hedin [1], 1
mjw Moruzzi, Janak, and Williams [2], 2
vwn Vosko, Wilk, and Nussair [3, 4], 3
lmmjw non-local + vbh, Langreth and Mehl [5, 6], 4
lmvbh not implemented, 5
pymjw non-local + mjw, Perdew et al [7, 8, 9], 6
pyvbh not implemeted, 7
pyvwn non-local + vwn, Perdew et al [7, 8, 9], 8
gga91 GGA, Perdew et al [10], 9
ev GGA only, Engel and Vosko [11]. 10

[XC1] http://www.ghfecher.de/Run_AKAI-KKR.pdf
[10] J. P. Perdew, et al., Phys. Rev. B, 46 (1992) 6671.: http://journals.aps.org/prb/pdf/10.1103/PhysRevB.46.6671,  PW91 と呼ばれる(wikipediaより)。
------------------------------------------------------------------------------
■ LDA+U
ΣHF = i * G * v だが、Veff = Veff + U * n(i) - double-counting の方が容易。
n(i) は chrdnc.f にあるので、そこからデータを potenv.f に持ってくる。
ro = ρ= ni なので、Veff = VLDAiUi(1/2 - ni ) は素人でも出来そうだが、特定の原子にUを入れられるが、特定の軌道には入れられない。F = E - ST = E - 0*T = E = -1.80 eV
rom = ρm= nim とできれば、Veff = VLDA + ΣimUim(1/2 - nim) までは出来るかもしれない。しかし、|Φim><Φim| がない状態での計算となる。

http://arxiv.org/pdf/1309.3355v1.pdf
Veff = VLDA + ΣimUi(1/2 - nim)|Φim><Φim|
上記の論文だと、Uは特定の原子iの特定の軌道mに入る。そして、ポテンシャルを計算する時に特定の原子iの軌道mを計算するときのみにUが入るようになっている。
※ 素人が見ると、この形式では、vを軌道毎にも分けなければならない。Akai-KKRでのポテンシャルはv(メッシュ、原子{成分}、スピン)となっており、軌道で分けられてない。LDA+Uの導入は数行書き換えだけでは済まなそうである。ある原子の特定の軌道にUをいれたいなら、vはpotenv.f に持ってきてはダメ。
------------------------------------------------------------------------------
■ ハートリーフォック近似の交換ポテンシャル
ΣHF = i * G(r,r') * v(r-r')
v: 裸のクーロン相互作用
------------------------------------------------------------------------------
■ ハイブリッド
Veff = VLDAxc + a*(ΣHF - VLDAx)
Akai-KKRでは potenv.f -> excg91.f において、upとdownスピンに対する交換相関ポテンシャル vxcupとvxcdn が計算される。行番号100の前で交換ポテンシャル vx が計算され、その後に相関ポテンシャル vcupとvcdnが vxcupとvxcdn に追加されているように見える。そのため、VLDAxc = vxcup and/or vxcdn,  VLDAx = vx として Veff = v としてしまえば良いように思えるが確認はしていない。
------------------------------------------------------------------------------
■ SIC
SIC: 自身の軌道でのVxcを差し引く。
------------------------------------------------------------------------------
■ 電場勾配(getefg.f)
d2V/dr2 = - (1/π) * Im GLL'
ここで、LL' は 02, 11, 22。GLL'はgreen.f にある。GLL'は通常φ0φ0, φ1φ1, φ2φ2などと計算されている。
電荷 ρ(r) = ΣLρL(r) * YL(r) にて L =2 が二回微分に対応。

下記は調査中
< do 10 m=1, mmx
< mm=mmxl*(m-1)+m
< 10 g(m,ji,k)=g(m,ji,k) + t(m,ji,k)*(1d0-t(m,ji,k)*wk(mm,2))*pf(m,m,ji,k)
---
> do 10 m=1, mmx
> do 10 m'=1, mmx
> mm'=mmxl*(m'-1)+m'
> 10 g(m, m', ji,k) = g(m, m',ji,k) + t(m',ji,k)*(1d0-t(m',ji,k)*wk(mm',2))*pf(m',m,ji,k)
------------------------------------------------------------------------------
■ 引数・配列の意味(初心者向け)
□ cpa2002v009c/source/specx.f の最初の部分にある配列
・natmmx: maximum number of atoms in a unit cell. (原子の数の最大値)
・ncnpmx: maximum number of types of atoms in a unit cell. (原子の種類の数の最大値)
・msizmx: maximum size of the KKR matrix. (Green 関数行列の最大サイズ)
・mxlmx: l ≥mxlmx are truncated. (mxl-1 の角運動量までで展開する。mxl の最大値)
・nk1mx,nk3mx: nk1mx+nk3mx is maximum number of k-points in the Brillouin zone. (nkmx = nk1mx+nk3mx がk 点の最大数)
  SCF や DOSでは nk1mxが、spcではnk3mxが用いられる。nk3mxを大きくしたければ、nk1mxを小さくする。
・msex: mesh points on the energy contour for go=dos, spc. (エネルギーメッシュの点数の最大値。go=dos や spc ではこの点数を使う)
・ngmx: The chebyshev expansion is performed up to ngmx-th order. (チェビシェフ展開の次数)
・nrpmx: maximum number of lattice points used in the Ewalt’s sum. (実空間の格子点の最大数)
・ngpmx: maximum number of reciprocal lattice points used in the Ewalt’s sum. (逆空間の格子点の最大数)
・msr: radial mesh points. (動径座標のメッシュの点数)
・mse0: (data) mesh points on the energy contour for go=go. ((deta 文)goでのエネルギーメッシュの点数)
・tol: (data) tolerance of convergence ((deta 文)ポテンシャルがtol まで収束したら計算を終了する)
------------------------------------------------------------------------------
■ 引数・配列の意味
・ids: 初期で0が代入される。dspで2, dosで1, mcdやxmdで3, spcで4, fsmで5が代入される。
・meshr: msr=動径座標のメッシュの点数の値がdata文で宣言される
・xr: メッシュ。イメージとして、メッシュ番号に対して指数関数的な感じでrが設定される。しかし、指数関数ではない。
  xr(k)=c*x**a*(x+g)**b+d (rmesha.f)
・dr: メッシュの微分
  dr(k)=(xr(k)-d)*(a/x+b/(x+g))/(dble(meshr-1)-e) (rmesha.f)
 ※ xr と dr があることで d/dr = dr/dk * d/dr が計算できる。
・e: 複素数でe(msex,2)と宣言される。msexはエネルギーメッシュの点数の最大値
  ref で価電子帯の割合を指定する。通常 ref = 0.75 としている。
  tfp.f(Thomas-Fermi potential)ではeは上記と異なる。
・de: ewidth/dble(kmx-1)
  kmx=mse=msex=エネルギーメッシュの点数の最大値, er=ewidth, ed=edelt
・di: dcmplx(0d0,edelt)
・ebtm: dble(e(1))で価電子帯の底のエネルギー。eb=ebtm
・a: mixing parameter
・sm: ∫n(e)*チェビシェフ多項式(e) de
・d: density
・intc: interstitial charge
・ints: intersititial spin
・is: スピン up は1, スピン down は2
・ro: 4 * pi * charge density (poisnb.f)
  ro(動径座標のメッシュの点数,単位胞内の原子の成分)
・z: solution v has the form v = -z(x)/x (poisnb.f) 
・anclr: point charge at the origin neutralizing the system (poisnb.f)
・a, b, xr: radial mesh must have a form xr(s)=a*(exp(b*s)-1) (poisnb.f)
・msr: s = 1,2,3,...,msr (poisnb.f)
・ns1: array size for r0,z,xr (poisnb.f)
・ncmpx: 単位胞内の原子の成分
・cm: fのチェビシェフ展開係数
・tchc, tchs: Tchebycheff expansion (prntcs.f)
  Tchebycheff expansion coeficients tchc,tchs and the shift constants fcs (cpshift.f)
・t, t(l): t-matrix (or its inverse) (cpshift.f)
・pf: phase factor fanction (cpshift.f)
・str: To facilitate the calculation of the phase of the Green's function this program also returns function 'str'. (cpshift.f)
・gs(l): the single potential Green's function. (cpshift.f)
・str(l): the quantity needed to calculate phase of the KKR Green's function. (cpshift.f)
・tc: Tchebycheff polinomials for a complex argument e. (cgntcd.f)
・td: Tchebycheff polinomials's derivatives for a complex argument e. (cgntcd.f)
・f: sum (1/(-1/t + b(struct) +i sqrt(e))) , ImGに相当するもの。
・tof: fの合計
・rstr: ΦL
・conc: 成分比
・cj: 
・ck: 
・ng: 出力ファイルの最初の方に値が記載されている。
・g: g(m,ji,k)=g(m,ji,k)+t(m,ji,k)*(1d0-t(m,ji,k)*wk(mm,2))*pf(m,m,ji,k)
  the local green's function for each site and component.
・ro: 電荷密度ρ(L=0)
・wk: work area (gsdatp.f)
・wkc: g と同じ。グリーン関数。wkc(方位量子数lのデータ, 原子, エネルギーメッシュ )
 方位量子数lのデータは、real spherical harmonics で格納されている。
  -----
  l=1: m=-1: px
  l=1: m=0 : pz
  l=1: m=1: py
  -----
  l=2: m=-2 : dxy
  l=2: m=-1: dyz
  l=2: m=0 : dz^2
  l=2: m=1: dxz
  l=2: m=2 : dx^2-y^2
  -----
  eg: dz^2 + dx^2-y^2
  t2g: dxy + dyz + dxy
  -----
  References
  [1] http://zurich.disneyresearch.com/~wjarosz/publications/dissertation/appendixB.pdf
  [2] http://www.sjbrown.co.uk/2004/10/16/spherical-harmonic-basis/
・w: ワーキングエリア( XCなどが入れられる)hp/calculations/page1
・v: ポテンシャル
・u: XC
・expsl: exchangeのスプリティング(↑- ↓)
・vint: interstitial site でのポテンシャルV
・ew: 計算するvalence stateの中心のエネルギー
・ez:  ew を中心にして、± ez の範囲を計算する
・emx: KKRで用いられるエネルギーの最大の絶対値(単位は{2π/a}^2) (etaopt.f)
・vc: 単位胞の体積 (etaopt.f)
・rmx: Ewald sum計算で用いるための実数部の最大長さ (etaopt.f)
・gmx: Ewald sum計算で用いるための逆格子ベクトルの最大長さ (etaopt.f)
・eta: 収束パラメーター (etaopt.f) (malng.f)
・atmicp: 単位胞内での原子の位置 (malng.f)
・rpt: 隣接する格子点 (malng.f)
・nrpt: rptで与えられる格子点の数 (malng.f)
・gpt: 逆格子ベクトル (malng.f)
・rmt: 各原子のマフィンティン半径 (malng.f)
・amdlng: 静電ポテンシャルを構築するために用いられる計算された定数 (malng.f)
・v(i): i番目のマフィンティンポテンシャルに対する静電ポテンシャル (malng.f)
・q(j): j番の原子の全電荷 (malng.f)
------------------------------------------------------------------------------
■ 使用されている組み込み関数
※ cmplex*16 で宣言された配列では、実数部はdble()を、虚数部はdcmplx()を用いる。
・ dble: 倍精度実数への変換
・ dcmplx: 倍精度複素数への変換
・ dimag: 倍精度の虚数成分
・ index:
------------------------------------------------------------------------------
■ F90以降への書き換え
 恐らく、COMMON文をUSE文に書き換えれば良いと考えられる。細かなエネルギーメッシュが欲しい場合には検討されると良い。
------------------------------------------------------------------------------
■ サブルーチン
※ サブルーチンへの引数の受け渡しはアドレス渡しである。配列は左側にある添字から処理される(メモリが順番に格納される)。
・gpdos.f: construct DOS data for cpa98 output.
・corconf.f: 内殻の電子数の設定(内殻に正孔を設定しても、全体の電子数は保たれるようになっている。
  up(i) に電子数が入る。iは入力した原子の順番。
・drvmsh: dosやmcd, xmd, spc のとき、cemesrのサブルーチンを呼ぶ。それ以外のとき cemeshを呼ぶ。
・chrasa: 動径波動関数と結合したsmデータから電荷密度を構築する。rmt球の内側の電荷は系の正しい電子数が得られるように再規格化される。
・chrdnc: 動径波動関数と結合したsmデータから電荷密度を構築する。
 ※ sm = ∫n(e)*チェビシェフ多項式(e) de
・spmain.f: pdos や tdosを出力する。Total DOSはIntegrated DOSを微分することで得ている。
  Total DOS と Partial DOSは計算方法が異なるので、Total DOS = ΣPartial DOS とはならない。また、PDOSの強度がTotal DOSを超えることも起こる。
 reltypを読み取る場合はlsがあればls=1を、nrlがあれば isr=0を、sraがあればisr=1となる。
・tchmta.f: 第一種チェビシェフ多項式での計算。メッシュ数 ng (ew±ezの範囲でng個に取ったサンプリングポイント)に対して、角度(ang)を2π/ng * k として、kをngまで1ずつ増やして計算する。elvl(k)=ew+ez*cosθ が計算される。
 Tn+1(x) = 2*x*Tn(x) - Tn-1(x) = 2 * T1(x) * Tn(x) - Tn-1(x) (wikipediaより)から、tm(l,k)=2d0*tm(2,k)*tm(l-1,k)-tm(l-2,k) となる。ここで、Tn(x) = cos(nt), x = cos(t) の定義からT1(x) = x である。この後、tm(1,k)は(1/ng)が掛けられ、lが2以上のtm(l,k)では(2/ng)が掛けられる。
・cgtable.f: cg.f によって得られる数がclksに格納される。c >= 0.0000001のときには、ncubにはmj(下記のsubscr.f のmlと同じ。mjは(2 * l l + 1)^2まである)が格納される。そうでなければ clks と ncub ともに 0 が入る。
・subscr.f: subscr(方位量子数、磁気量子数、ml) で、mlの番号から方位量子数と磁気量子数を取り出す。
 式は次の通り、方位量子数 = √(ml - 1) + 0.00001, 磁気量子数 = ml - l * ( l + 1) - 1となる。ただし、方位量子数の計算では ( ml-1) < 0 のときは (ml - 1)を0 とする。
・cg.f: クレブシュ-ゴルダン係数, gount number
・cgc.f: クレブシュ-ゴルダン係数の計算。Clebsh-Gordan coefficient
 ※ cg.f と cgc.f はどちらかがspherical harmonics、残りの一つは real harmonics である。
・etaopt.f: madlng.f で用いる収束パラメーターを計算? 
・genrpt.f: 3つの基本並進ベクトル(primitive translation vectors)を与えて、格子ベクトルのセットを作る。
・madlng.f: マーデルングポテンシャルに関する計算
・hexmtr.f: real harmonics に対する六方対称の回転行列を作成する?
・getrtm.f:
・realh.f: r^l * real spherical hermonics(yに出力される)を計算する。入力fはベクトル v と lmax。
・vrotat.f: 入力をv1とし、オイラー角による座標変換に対応する回転操作を行った結果をv2とする。
・drvlud.f:
・transp.f: ランク jmx の行列aを転置する。※ 入力も出力も同じ配列になるので注意。
・cubmtr.f:  real harmonics に対する立方対称の回転行列を作成する?
・uutrns.f:
・brvsym.f: ブラベー格子に対する対称操作を与える
 スピン軌道相互作用を考慮にいれると対称操作が変わる
・getnum.f:
・stfact.f: 
・atmrot.f:
・bzmesh.f: BZ内でのメッシュを作成する。
・rndkpt.f:
・gnpset.f:
・gsdatp.f: SCF計算における初期ポテンシャルを作成する。Mattheiss' prescriptionが用いられる。LSDによってセルフコンシステントに原子ポテンシャルが作成される。
・potenv.f: XCポテンシャルの計算
・finitn.f: 原子番号(anclr)からポテンシャル上での有限な核サイズ補正を作る
・phasea: スピン軌道相互作用を考慮しない場合、こちらがspmain.f から呼ばれる。
・phaseb: スピン軌道相互作用を考慮する場合、こちらがspmain.f から呼ばれる。
 cは光速度(Ry単位)。
 配列の期化を行った後、ポテンシャル v を微分して wk に格納する。
・phase.f: ClやSlを求める。t-matricsを計算する? 
・errtrp.f: エラー表示用のサブルーチン。エラーの種類によって表示を分け、ルーチンとメッセージも準備してある。
・clrarr.f: clarr(配列、配列数)は、配列数nである倍精度実数(real*8)配列の中身を0.0 (= 0d0 )に初期化する。
・diffn.f: diffn(関数、微分した関数、配列 xr、刻み、配列数)は、配列の最初(4点公式)と中央(5点公式)と最後(-1 * 4点公式)の部分に分けて微分の計算を行う。
 ※ 配列 xr は微分の計算では使われていない。diffn に入れる引数は最後のものを除いて、最後の引数の配列数である。
・drivtv.f: 微分値を計算する。ルーチンの中では5点近似公式などが用いられている。
・mdbsl.f: コメント欄に詳細な説明がある。変形ベッセル関数とノイマン関数、それらの微分に関するサブルーチン?
・radial.f: Adams-Moulton法が使われている(常微分方程式の数値解法である予測子修正子法の一種。簡単には y' =f(x,y) で、yを求める方法と理解してみるとよい)。任意のメッシュで良くてスピンオービタルを入れた計算を比較的簡単に記述できる。
 ※ 類似の方法としては、Hamming法がある。d2L/dX2 = y' として dL/dX = y を得る。dL/dX = y' として L = y を得る。
・poisnb.f: 各引数の説明がある。
・reduce.f: sm行列から各エネルギー点に対する重みを生成する
・rrconv.f: r(up,down) --> r(charge,spin) conversion 
・sbtime.f: Count cpu time used in a specific subroutine.
・ifkey.f: ifkey(key, 文字列)は、組み込み関数 index を用いて、入力ファイルから読み取った文字の中にkeyとなる文字があるかを調べる。あればtrueを返す。
・leftfil.f: lerftil(文字列)は、組み込み関数lenで文字列の文字数を調べ、左側の空白を除いて左に詰める。詰めて残った部分は空白になる。
・inijip.f: jip マッピング(j,i)-> ji を与える。jは成分、iはtype。
------------------------------------------------------------------------------
■ ファイル番号
unit 5: 入力ファイル
------------------------------------------------------------------------------
■ 状態密度の計算
・ 部分状態密度を計算するには、Green 関数を複素平面上で実軸に沿って等間隔にメッシュを取る。虚数の値はedeltで一定。edeltはローレンチアンでの半値幅と概念的に考えることができる。PDOSでは同じローレンチアンでなましていることになる。
・ DOSが負になることがあるのは、iがあるためである。
・ 複素空間での上半平面で積分し、始点と終点が同じなら、積分経路上には極はないものとすると、積分経路によらず結果は同じとなる。
※ dosやpdos, spcはEFを0として出力されるようになっている。
□ drvmsh.f が cemesr.f を呼んで、複素空間上で、実軸のメッシュを切る。
------------------------------------------------------------------------------
■ 数値計算法
・4点公式
・5点公式
・Adams-Moulton法
 任意のメッシュで計算でき、スピンオービットを入れた計算の記述が比較的容易である。
・第一種チェビシェフ多項式
 理学上、工学上そのほかの分野で利用頻度の多い関数に直交多項式または特殊多項式とよばれるものがある。その中でもよく使われる多項式にはルジャンドル(Legendre)多項式、ラゲール(laguerre)多項式、エルミート(Hermite)の多項式、チェビシェフ(Tchebycheff)の多項式などがある。
 これらの関数はべき乗項の有限項の多項式であるので数式をFortranで書き出せばそのまま関数プログラムとなると考えるかもしれない。しかしこれらの多項式は各項ごとに符号が正負交互に交代する交代多項式であるために、項数が多くなり引数の絶対値が大きくなると、各項の相殺(キャンセレイション)が著しくなり、精度が出ず正しい数値を与えない。このためにべき乗の有限項多項式であるが専用の関数副プログラムが必要となる。(渡部力ら、「Fortran77による数値計算ソフトウェア」、丸善 より)
・クレブシュ-ゴルダン係数
・Mattheiss' prescription
------------------------------------------------------------------------------
■ アルゴリズム
------------------------------------------------------------------------------
■ 理論背景
[1] http://pmt.sakura.ne.jp/wiki/images/Cmdtxt1.pdf
[2] http://kkr.phys.sci.osaka-u.ac.jp/pdf/akaikkr_j.pdf
[3] http://kkr.phys.sci.osaka-u.ac.jp/pdf/machikaneyama2000_memo.pdf
[4] http://www.psi-k.org/Psik-training/Ebert.pdf
------------------------------------------------------------------------------
■ 書き換え例
※ 参考にして頂ければと思いますが、結果については自己責任でお願い致します。

■ PDOS (spmain.f)
※ is で up spinとdown spinを指定しますが、pdosの引数であるwkcにはそれが無いように見えます。しかし、下記の部分が存在しており、up で wkc を得て出力した後に、downでもwkcを得て出力するようにしています。
c     --- perform kkr --
      call clrarr(sm,2*ng*mmxl*ncmpx)
      do 70 is=1,nspin
  (省略)
  PDOSなどの出力
   70 continue
下記は方位量子数を分けて出力するための書き換え方の一例です。

before
---------
c     --- print partial and the total DOS if required.
      if(ids .eq. 1 .or. ids .eq. 2 .or. ids .eq. 3) then
      estep=dble(e(2,is))-dble(e(1,is))
      do 69 i=1,ncmpx
      write(*,'(//1x,a,i2,a,i2)')'DOS of component',i
      do 69 k=1,kk
      xmd(i,k,1,is)=-dimag(wkc(2,i,k))/pi
cc    &    (-dimag(wkc(4,i,k))/pi)+(-dimag(wkc(2,i,k))/pi)
      xmd(i,k,2,is)=-dimag(wkc(4,i,k))/pi
cc    &    (-dimag(wkc(4,i,k))/pi)-(-dimag(wkc(2,i,k))/pi)
c  69 write(*,'(1x,f7.4,3x,9f8.4)') dble(e(k,is))-ef(is)
c    &      ,( -dimag(wkc(l,i,k))/pi,l=1,mxl**2)
      do 160 l=1,mxlcmp(i)
c     do 160 l=1,2
      do 160 m=1,2*(l-1)
  160 wkc(l**2,i,k)=wkc(l**2,i,k)+wkc(l**2-m,i,k)
c     wkc(5,i,k)=wkc(5,i,k)+wkc(6,i,k)+wkc(8,i,k)
c     wkc(7,i,k)=wkc(7,i,k)+wkc(9,i,k)
   69 write(*,'(1x,f7.4,3x,4f10.4)') dble(e(k,is))-ef(is)
     &      ,(-dimag(wkc(l**2,i,k))/pi,l=1,mxlcmp(i))
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'total DOS',(dble(e(k,is))-estep/2d0-ef(is)
     &      ,dimag((detl(k,is)-detl(k-1,is))/(e(k,is)-e(k-1,is)))
     &        ,k=2,kk)
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'integrated DOS',(dble(e(k,is))-ef(is)
     &       ,dimag(detl(k,is)),k=1,kk)
---------

after
---------
c     --- print partial and the total DOS if required.
      if(ids .eq. 1 .or. ids .eq. 2 .or. ids .eq. 3) then
      estep=dble(e(2,is))-dble(e(1,is))
      write(*,'(///a)')
     & '(PDOS DATA: Ry, s, px, pz, py, dxy, dyz, dz^2, dxz, dx^2-y^2)'
      do 69 i=1,ncmpx
      write(*,'(//1x,a,i2,a,i2)')'DOS of component',i
      do 69 k=1,kk
      xmd(i,k,1,is)=-dimag(wkc(2,i,k))/pi
cc    &    (-dimag(wkc(4,i,k))/pi)+(-dimag(wkc(2,i,k))/pi)
      xmd(i,k,2,is)=-dimag(wkc(4,i,k))/pi
cc    &    (-dimag(wkc(4,i,k))/pi)-(-dimag(wkc(2,i,k))/pi)
   69 write(*,'(1x,f7.4,3x,9f8.4)') dble(e(k,is))-ef(is)
     &      ,( -dimag(wkc(l,i,k))/pi,l=1,mxl**2)
c      do 160 l=1,mxlcmp(i)
c     do 160 l=1,2
c      do 160 m=1,2*(l-1)
c  160 wkc(l**2,i,k)=wkc(l**2,i,k)+wkc(l**2-m,i,k)
c     wkc(5,i,k)=wkc(5,i,k)+wkc(6,i,k)+wkc(8,i,k)
c     wkc(7,i,k)=wkc(7,i,k)+wkc(9,i,k)
c   69 write(*,'(1x,f7.4,3x,4f10.4)') dble(e(k,is))-ef(is)
c     &      ,(-dimag(wkc(l**2,i,k))/pi,l=1,mxlcmp(i))
      if(is .eq. 1) then
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'total_up TDOS_up',(dble(e(k,is))-estep/2d0-ef(is)
     &      ,dimag((detl(k,is)-detl(k-1,is))/(e(k,is)-e(k-1,is)))
     &        ,k=2,kk)
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'integrated_up IDOS_up',(dble(e(k,is))-ef(is)
     &       ,dimag(detl(k,is)),k=1,kk)
      end if
      if(is .eq. 2) then
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'total_dn TDOS_dn',(dble(e(k,is))-estep/2d0-ef(is)
     &      ,dimag((detl(k,is)-detl(k-1,is))/(e(k,is)-e(k-1,is)))
     &        ,k=2,kk)
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'integrated_dn IDOS_dn',(dble(e(k,is))-ef(is)
     &       ,dimag(detl(k,is)),k=1,kk)
      end if
---------
------------------------------------------------------------------------------
■ SPC for PDOS at one k point (bzmesh.f)

before
---------
      do 50 i=1,3
   50 vkp(i,1)=zero
---------

after [e.g. k point (0.5 0.25 0.123)]
---------
c      do 50 i=1,3
c   50 vkp(i,1)=zero
      vkp(1,1)=0.5
      vkp(2,1)=0.25
      vkp(3,1)=0.125
---------
then, bzqlty = 0.
------------------------------------------------------------------------------
hybrid functionals (potenv.f) (under constructing)
------------------------------------------------------------------------------
QRコード
携帯用QRコード
アクセス数
ページビュー数
[無料でホームページを作成] [通報・削除依頼]