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

入力ファイルのQ&A

 ここでは、バンド計算の入力ファイル作成において、よく疑問となる事項についてまとめていきます。
------------------------------------------------------------------------------
■ 一般的な入力ファイルのQ&A
GAMESS-US, Gaussian -> Gabedit QC-GUI

Q. 周期律表はどこから手に入れればよいですか?
A. 文部科学省HPから周期律表が見られるようになっています。
http://www.mext.go.jp/a_menu/kagaku/week/shuki.htm 

Q. 波動関数のカットオフ(cut off)の値はいくつに設定すれば良いか?
A. ウルトラソフト擬ポテンシャルの場合は20〜25 Ry 以上に設定して下さい(ノルム保存型擬ポテンシャルの場合はこれより3〜4倍にしたりする)。詳細にこの値を決定するには、結晶の計算をして、格子定数や体積弾性率の実験値を再現するようにします[1]。この他、徐々に値を大きくして結果が収束したところで値を決定させることもよく用いられます(Phaseでは擬ポテンシャルファイルに予め値が入力されており、特に設定する必要がありません。一方、VASPでは推奨値が擬ポテンシャルファイル中に記述されています{詳しくはVASPのマニュアルを参照下さい})。WIEN2でのCutoff は 13.6 * (RKmax / Rmt ) ^2 [eV]で得られます(http://www.jim.or.jp/journal/j/pdf3/73/09/683.pdf とも比較してみると良いでしょう。ちなみに、Rmtは最小の値で計算します。計算してみると13.6*(9.0/2.0)^2=275.4 eVとなります)。
よく用いられる擬ポテンシャルの種類は下記の「 電荷密度分布のカットオフ」に関する記述を参照。
※ VASPでのPAWに対するカットオフについては下記の報告も参考になる。defaultから600 eVを採用すべきか?
※ 多くの先生方に伺うと 500 eV 〜 520 eV の方が多い。
※ WIEN2kとの比較がされていて十分な精度なので、私の場合、最近では推奨値を用いています(左欄の「各第一原理計算コードの特徴」を参照してください)。
[VPAW0] http://wolf.ifj.edu.pl/workshop/work2004/files/talks/kresse_vasp_phonon.pdf
[VPAW1] http://www.nuclear.jp/~kyodo/report/0022/F/18F-13.pdf
[VPAW2] https://www.jstage.jst.go.jp/article/msjmag/34/2/34_1003R004/_pdf
[VPAW3] http://www.toyotariken.jp/activities/65/S09_yuhara.pdf
[VPAW4] http://books.google.co.jp/books?id=yOznjeSms18C&pg=PA9&lpg=PA9&dq=PAW+%E3%82%AB%E3%83%83%E3%83%88%E3%82%AA%E3%83%95&source=bl&ots=fyrsls8w3s&sig=nWIw9-LW0_8asl0Pj6UlSrSMY98&hl=ja&sa=X&ei=b6CxUcv6H4uElAWmzYGwAw&ved=0CFkQ6AEwCTgK#v=onepage&q=PAW%20%E3%82%AB%E3%83%83%E3%83%88%E3%82%AA%E3%83%95&f=false
[VPAW5] http://nsed.jaea.go.jp/fme/group3/group3_i5.pdf
[VPAW6] http://kgur.kwansei.ac.jp/dspace/bitstream/10236/9722/1/M0349.pdf

下記のHPの資料(PWscf)にcut offの値を変えた場合の結果が示されています。
http://www.fisica.uniud.it/~giannozz/QE-Tutorial/tutorial_pwscf_ex.pdf 
2 Ry = 1 Hartree = 27.2113845 eV
1 bohr = 1 a.u. = 0.05291772108 nm
*カットオフの意味については、下記の中山先生のHPでのstudyを参照するとよい。
http://homepage3.nifty.com/mnakayama/index.htm

Q. 電荷密度分布のカットオフ(cut off)の値はいくつに設定すれば良いか?
A. 波動関数のカットオフの4倍以上[1]
・ノルム保存型擬ポテンシャル:波動関数のカットオフの4倍以上。この種類で用いられる擬ポテンシャルではTroullier-Martin(TM)(トウリエ-マーティン)が有名。
・ウルトラソフト擬ポテンシャル:波動関数のカットオフの4倍以上、典型的には8倍-12倍の値。この種類で用いられる擬ポテンシャルではVanderbilt(ヴァンデルビート)が有名。
・PAW:波動関数のカットオフの4倍以上
※ Phaseでは擬ポテンシャルファイルに予め値が入力されており、特に設定する必要がありません。

Q. qパラメーター(フォノンや吸収の計算で入力することになる)はどの値にすればよいの?
A. 私も良く分からなかったのですが、qは逆格子空間上でのメッシュの間隔に対応しており、その間隔を実空間に戻したときに(ディスプレイスメントが)計算される範囲(10Å以上あればよい)と考えれば良いようです。格子定数が大きければガンマ点(1点)のみでもよくなります。
※ フォノンの計算方法は幾つかありますが、基本的なものとしては、スーパーセルを作って、一部の原子を少し動かしたとき(ディスプレイスメント)の力(フォース)を利用したりします。そう考えるとqをどう考えれば良いかが分かりますなりますね!

Q. bandの数はどれくらいに設定すればよいの?
A. Phaseでは計算する系での価電子帯の全電子数の2割強程度増やした値(ノンスピンの場合はその半分)を入力してください[1]。ABINITでは"nband=占有+非占有のバンドの数"です。(PhaseやABINIT, PWscf, WIEN2kは自動的にデフォルト値が入力される。WIEN2k_11では、デフォルトで、エネルギー範囲内での単位胞あたりの電子数×2.0+5に設定されます。ノンスピンの場合はその半分の値)

Q. LDAとGGA、PBEはどんなときに使い分けるのか?
A. 一般的に PBE(が無ければGGA)を用いて下さい。しかし、希ガス固体とファンデルワールス(vdW)力が関係するグラファイトなどではLDAを採用する方が計算結果が実験に合います。最近の論文(2011年)でもグラフェン(グラファイトが一層になったもの)を扱った系ではLDAで計算されているものがあります。最近では、vdWが可能な第一原理計算コードが普及してきており、PBE(が無ければGGA)を用いてもvdWを入れた計算が可能なものがあります。余談ですが、CASTEPであればSiなどでのギャップの値がLDAよりも改善されるsx-LDAがあります。

Q. LDA+UやGGA+U、PBE+U はどんなときに使うの?
A. 3d 軌道などの局在する軌道に対して用います。そうすると、ギャップや電子構造が実験値と合う計算結果を得ることができます。SrTiO3でのTi 3d 軌道などを例として良く使われています。(effective U for WIEN2k)

Q. 3d 軌道などの局在軌道が無い系ではどうすればよいの?
A. 酸素の2pなどにUを入れている例も見掛けます。この他に、計算に時間が掛かりますが、PBE0, YS-PBE0, mBJ, sx-LDAなどを用いることも考慮してみては良いのではないでしょうか。


Q. スカラー相対論 (scalar relativestic calculation) とはなんですか?
A. 非相対論 (non relativestic calculation) に、質量・速度項 と ダーウィン項を加えた計算です。Full relativestic calculation は、スカラー相対論 に スピン・軌道相互作用項を含めたものです(WIEN2kでは Spin Orbit couplingを計算するために入れる項になります)。

Q. XPSなどでよく見る内殻の電子状態を示す2P1/2や2P3/2などは理論的にどのよう考えれば良いでしょうか?
A. スピン軌道相互作用の項を理解すると良いでしょう。スピン軌道相互作用の項は次式になります。
  ΔEnlj = (1/2) * Anl * (h/2π)^2 * { j * (j+1) - l * (l+1) - 3/4}
ここでAnlは複雑な式なので詳細を省略します(文献[SOC1-3]を参照して下さい)。j = l + s です。このエネルギー準位は2j+1重に縮退しています。内殻の電子の場合で、l = 1の場合の例を下記で見てみましょう(他のlの値でも同様です)。
  j = 3/2 の場合: P3/2、ΔEnlj =(1/2) * Anl * (h/2π)^2: 2j+1 = 4 (-3/2, -1/2, 1/2, 3/2 の4つ)
  j = 1/2 の場合: P1/2、ΔEnlj =-Anl * (h/2π)^2: 2j+1 = 2 (-1/2, 1/2 の2つ)
j = 1/2 の方が内殻に近いエネルギーになります。スピン軌道相互作用の式から、電子のエネルギーは P3/2 > P1/2 であると理解し、また、エネルギー準位が2j+1重に縮退していることからP3/2 とP1/2の強度比が、2j+1の式から、2*(3/2)+1 : 2*(1/2)+1= 4 : 2 = 2 : 1になることも覚えておくとよいでしょう。重要なこととして、upとdownがそれぞれの軌道に均等に入ります(lとsを負まで考慮に入れて j = l + s の表を作ってみてください{正確にはj に対するmの表、つまり、mj = ml + ms が正しいでしょう})。-3/2, -1/2, -1/2, 1/2, 1/2, 3/2 の6つの値が得られます。これらのどれがP3/2 とP1/2に属するかは ガシオロウィッツ、「量子力学U」、丸善株式会社など、色々な参考書を読んで勉学してみてください)。l=0や1, 2などでも同じ様に考えることができます。また、WIEN2kのcase.instを見ると内殻の電子はupとdownがそれぞれの軌道に均等に配置されていることが分かります。
[SOC1] 中嶋貞雄、「量子力学U」、岩波書店
[SOC2] http://202.11.2.126/iwaya/ryoushi/h19/12-2.pdf
[SOC3] http://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/96501/1/KJ00004707301.pdf

Q. PWscf での擬ポテンシャルファイルの読み方は?
A. 下記のHP[1,2]にて解説されています。求める計算の条件に合う擬ポテンシャルを選んで下さい。
[1] http://www.quantum-espresso.org/wiki/index.php/Naming_convention_for_the_pseudopotential
[2] http://www.quantum-espresso.org/pseudo-notes.php

Q. conventional unit cell と primitive unit cell はどう使い分ければいいの?
A. 下記のHP[1]にて conventional unit cell と primitive unit cell での入力する原子座標の違いが示されています。conventional unit cell は原子座標の入力が直観的に分かりやすいですが、一方、primitive unit cell は入力する原子の座標は少ないものの初学者には難しいと思います。 primitive unit cell では入力する原子座標の数が少ないので、計算時間は短くなると思われます。しかし、primitive unit cell に慣れるまでは、構造を確認し易い conventional unit cell 形式での入力ファイルの作成が良いと思います。ちなみに、cif 形式のデータを読み取って計算するWIEN2k では 私の場合 conventional unit cell での原子座標の入力になり易いです。この他にも、Osaka2k などでは conventional unit cell での原子座標の入力になっています。初めは、分かりやすい conventional unit cell で計算をしてみて、徐々に primitive unit cell に挑戦して、結果に違いが無いかを調べながら様々な計算を行っていくと良いと思います。(xcrysdenでconventional unit cell と primitive unit cell の表示を切り替えられます。そこから自動的に入力ファイルが作成できるかは未調査)
[1] http://mpdc.mae.cornell.edu/Courses/MAE715/HwsPDF/HW5/PWSCF-Presentation.pdf 

Q. k点数の設定はどうすればよいのか?
A. ユニットセル内の原子数が多くなるほど、k点は少なくて良い。論文の値などを参考にするもの良いでしょう。
SCF case
  Metal: 0.035 [1/Angstrom]
  Semi-conductor: 0.088 [1/Angstrom]
  Insulator: 0.088 [1/Angstrom]
DOS case
  Metal: 0.035*2 ~ 0.035*4 [1/Angstrom]
  Semi-conductor: 0.088*2 ~ 0.088*4 [1/Angstrom]
  Insulator: 0.088*2 ~ 0.088*4 [1/Angstrom]
a) Phaseの場合
自動的に処理してくれるので、特に気にする必要はない。
b) ABINITやPWscf の場合
格子定数が3.91オングストローム程度で4原子程度だと6x6x6の216点程度。
c) WIEN2kの場合
格子定数が13オングストローム程度で140原子程度だと、500点で十分なDOSを描くことができる。一般のAl金属などでは1000点。グラフェンなどでは5000点の例がある。8000点を超えたk点数を採用することはほとんどない(SCF計算後に[initialize calc.]のx kgenでk点を増やし、run SCFを実行させ、これらの計算結果が飽和するまで徐々にk点を増やしてみるのがよい。k点の数を多くしていき、3 cycleで収束していくようになれば、それ以上k点を増やさなくてもよい一つの目安になると思います)。

Q. スーパーセルでのk点はどれだけ取ればよいのか?
A. スーパーセルにする前の単位胞で採用したk点を、スーパーセルにするために拡大した単位胞の倍率分だけ割って下さい。例えば、2×2×1倍のスーパーセルにしたら、k点は1/4倍とすれば良いです。半導体や絶縁体は金属よりもk点を少なくしてよく、SrTiO3を3x3x3のスーパーセルとした系などではk点が8点でよい(スーパーセルの作り方は参考文献[9]を参照してください)。

Q. 表面の計算をするために真空層はどれだけ必要か?
A. 10オングストロームほど真空層(何もない空間)を作って下さい。論文でもだいたい10オングストロームが採用されています。
※ ただし、表面で分極が大きい場合は、20〜30オングストロームにすることもあります。

Q. 格子定数の決定はどうすればよいか?
A. 理論屋さんであれば、第一原理計算の範囲内で最適化により決定します。実験屋さんであれば、実験値を用いることもあるでしょう。実験値を用いて計算した結果が、インパクトファクターの高い論文に記載されていることも多々見かけます。

Q. 各種パラメータはどうやって決定すればよいか? 計算の手順を教えて下さい。
A. 論文があればそれをまず参考にして同じ結果が出力されるか調べる。論文が無ければ、格子定数に実験値を入力して固定し、k点→カットオフ→k点と繰り返しながら値を徐々に大きくしていって結果が収束するまで繰り返す。次に、格子定数の最適化を行う。対称性の関係などから(Γ点など)、k点を偶数か奇数で増やしていくかは注意深くみておく必要がある。これは形成エネルギーを検討する場合に条件を揃えるためにも有効である。


Q. 書籍にあるデータベースの構造を入力したのに計算できません。
A. NIMS物質・材料データベース(http://mits.nims.go.jp/)を利用してみて下さい。それ以外にも、左欄のデータベースを参照してみて下さい。

Q. WIEN2kでは -c が表示されたり、表示されなかったりします。これは何を意味しているのでしょうか?
A. 平面波に虚数成分が入って計算されているときに-cが表示されます。
単純な平面波の式を見てみましょう(結晶での計算では結晶の周期性を反映させるためにさらにブロッホ関数が掛けられる)。
exp(ikr) = A*cos(kr)+i*B*sin(kr)
例えば、[r = -r]で同じ原子がある場合、[exp(ikr) = exp(-ikr)]となるような結晶があるとします(これは時間反転対称性と呼ばれる)。
式をより細かく見てみると、
A*cos(kr)+i*B*sin(kr) = A*cos(-kr)+i*B*sin(-kr) = A*cos(kr)-i*B*sin(kr)
2*i*B*sin(kr) = 0 となります。
実数部分のAは任意の値で上記を満たします。しかし、虚数成分のBは0のときにどのようなkrでも上記を満たすことになります。
 このように時間反転対称性が成り立つ系の場合には、虚数成分が0になって計算の必要がなくなるために、WIEN2kなどでは実数のみでの計算がなされます。
 時間反転対称性が無い場合にはWIEN2kで計算中に -c と表示され虚数成分も計算されます。


Q. MPI では -np と OMP_NUM_THREADS の値をどのように設定すればよいの?
A. -np の後に入れる数字はCPUの型番を調べてコア数を記入して下さい。そして、OMP_NUM_THREADS では [(OS上で表示されているCPU数) / (-npで入力した数値)]の値を指定して下さい。PWscf では アウトプットファイルにMPIに関する記述が書き込まれますので確認してみて下さい。下記に core-i7 2700K  のCPU でPWscf を動作させる場合の例を示します。
1) export OMP_NUM_THREADS=2 を .bashrc に書き入れる。
2) /usr/local/openmpi/bin/mpirun -machinefile /home/*/hosts -np 4 /home/*/espresso-4.3.2/bin/pw.x atom.out
以上で、2 threads * 4 CPU = 8 processor cores の動作となります。
※ abinit などにおいてもほぼ似たようなコマンド形式になるので同様に設定すればよいです。
※ OpenMP だけの計算の場合は OMP_NUM_THREADS=OS上で表示されているCPU数 としても良いでしょう。
※ -machinefile のオプションで、アドレスを指定した hosts ファイルに MPI で計算させる PC の名称を記入します。その優先順位に添って計算がされていきます(同じPCの名称が何行続いてもOK。1台のPCで MPIの並列計算をする場合、筆者は core-i7 2700K が 4 core なので 4行同じ PC の名称を記入したりする)。
※ OMP_NUM_THREADS=1 とし、 -np 直後に入力する数値を実際のCore数-1とする方が計算が速いこともある。どの指定が計算が速いかは似た系で調査してみるしかない。

Q. SrTiO3 で Ti に対して 約1% だけ遷移金属をドープする場合に、PCのメモリ数はどれだけ準備しておけばよいの?
A. CASTEP では メモリが24GB実装されていれば、500原子を扱えると言われています。SrTiO3などでTiに対して約1% 置換する系を計算するには、スーパーセルで 4x4x4 の構造が必要になります。これは 320 原子の計算になりますから、24 GBのメモリを扱えるPCを準備する必要があります。(ちなみに、Ti に対して約4%置換する系ですと3x3x3 で 135 原子なので、8〜16 GB 必要になります。core-i7 2700K で メモリ16GBのPCを用いて2 threads * 4 CPU で 3x3x3 構造を PWscf で計算したところ、格子定数も含めた構造最適化をすると1日から8日程度掛かります。遷移金属のドープが多いほど収束までの日数は長くなりました。ここで、擬ポテンシャルはPAW以外のもので主にVanderbilt型を、交換相関項はPBEを用い、spinを考慮した計算としている)

Q. ギャップのある系でバンド分散を描くと、本来、直接遷移であるべきなのに、間接遷移の結果が表示されます。どのように対処すれば良いでしょうか?
A. 格子定数を僅かに変えると、直接遷移や間接遷移に変わると言われています。構造最適化した場合での結果と、XRDなどで得られた実験で計算した場合のどちらかが正しい場合は、その旨を論文に記して、正しい値が得られた方で議論していけば良いでしょう。この他、k点の数を増やしたりしてみます。もし、あなたが擬ポテンシャルを用いて計算している場合は、擬ポテンシャルを他のものに変えてみたりもしてください。また、FLAPW法(WIEN2kやexciting)で計算してみるのも良いでしょう。ちなみに、擬ポテンシャルの作成は理論計算で論文を書いているプロでも難しいので、擬ポテンシャルを作成する場合には、擬ポテンシャル作成のためのプログラムを書いている方やあたなたが計算している第一原理計算コードのプログラムを書いている方に相談してください。

Q. 変分原理から考えるとSCFサイクルの最初から最後へはエネルギーが小さくなるはずですが、そうならない結果もあります。どうしてなのでしょうか?
A. イフェクティブポテンシャル(Veff)がSCFの初期と最後で異なるためです。Veffが変化するとハミルトニアンも変化します。もし、SCF計算で初期から最後までハミルトニアンが固定されていた場合にはエネルギーが高→低になります。

Q. WIEN2kなどで計算するとギャップが開いた結果になりますが、実際に試料を作製した結果ではギャップがつぶれています。それは何故でしょうか?
A. 理想的な完全結晶を実験で作ることは困難です。そのため、欠陥や不純物、または非化学量論組成になることによる準位がフェルミ端近傍に形成されることがあります。KKR-CPAなどで検討してみてください。

------------------------------------------------------------------------------
■ PWscfでの入力ファイルのQ&A

Q. PAWとUSPPの両方を用いているのですが、SCF計算が走りません。
A. PWscfでのQ&Aでも擬ポテンシャルは混ぜて用いても良いことが書かれていますが、上手く動作しない経験を持った方は多いのではないでしょうか。参考までに筆者が成功した例を下記に示します。
[System] -> [Optional variables]
--- Noncolinear calculation ---
Perform noncolinear calculation (noncolin): Yes
Allow the use of a pseudopotential with spin-robit (lspinorb): Yes
上記の条件で計算可能になった擬ポテンシャルの例を下記に示します(注:VはUSPP)。
Fe.rel-pbe-kjpaw.UPF
V.pbe-sp-van-UPF
Al.rel-pbe-n-kjpaw_psl.0.2.2.UPF

Q. 続きから計算するためには?
A. PWgui → Control → Restart mode (restart_mode): from previous interruptied run を選択後、上の欄にあるRunを押せば続きから計算します。多くの場合、100 cycleで収束しなかったことによる計算の停止が多いと思います。SCFの計算回数の制限を多くするためには、Electrons → Optional Variables → Max. # iterations in SCF step (electron_maxstep): に大きな数値を入力しておくことも良いでしょう。 
------------------------------------------------------------------------------
■ ABINITでの入力ファイルのQ&A

Q. ABINITでの入力ファイルのvariablesのどれを選べばよいか分かりません。
A. 物質・材料研究機構の新井先生が書かれた「ABINITによる第一原理電子構造計算例」http://www.nims.go.jp/cmsc/staff/arai/memo/kyushu/example1.pdf及び「ABINITによる第一原理電子構造計算例:ABO3http://www.nims.go.jp/cmsc/staff/arai/memo/kyushu/example2.pdfを参考にすると、以下のようになります。

GWは http://www.powershow.com/view/f0f8f-NjNkY/Abinit_Workshop_flash_ppt_presentation を参照

------------------------------------------------------------------------------
case.out または case.log で、認識された空間群とVariableの設定状況の確認方法
------------------------------------------------------------------------------
入力された Variable の確認(下記の2通りの方法がある)
A) case. out で下記を確認
  1) Symmetries : space group 空間群 (# 空間群の番号)を調べる。
  2) Values of the parameters that define the memory need of the present rum
    直ぐ下にある This job should need less than の後の使用メモリも見ておくと良いだろう。
  3) -outvars: echo values of preprocessed input variables ------ 
  4) 以上で自動で設定された値も含めて、入力された variable の値が分かる。
   目的のものと違っていた場合は、case.in に書き入れる。
  5) 上記になかった variable は ctrl + F で検索して、どのように使われているか調べる。
B ) case.log で下記を確認
  1) symspgr : spgroup= で空間群の番号を調べる。
  2) Values of the parameters that define the memory need of the present rum
    直ぐ下にある This job should need less than の後の使用メモリも見ておくと良いだろう。
  3) -outvars: echo values of preprocessed input variables ------ 
  4) 以上で自動で設定された値も含めて、入力された variable の値が分かる。
    目的のものと違っていた場合は、case.in に書き入れる。
  5) 上記になかった variable は ctrl + F で検索して、どのように使われているか調べる。
  ※ mband = nband * nsppol
※ Fermi energy は case.out に記載されている。

case.filesには、以下を記述する。()内は記述しない。
./A/case.in とファイルAの中にある指定も可能。
--------------------------------------------------
case.in (パラメーターの入力ファイル名)
case.out (出力ファイル名)
casexi (その他の入力ファイルのヘッダー)
casexo (その他の出力ファイルのヘッダー)
temp (一時ファイル用)
06-C.LDA.fhi (擬ポテンシャルファイル名)
--------------------------------------------------

case.in
--------------------------------------------------
acell 7.38 7.38 7.38 #セルの長さの設定。 3*7.38 でもよい。(a.u.単位)
angdeg 90 90 90 #基本ベクトル間の角度(a2とa3の間の角) (a1とa3の間の角) (a1とa2の間の角)
ntypat 3 #原子種の数
znucl 38 22 8 #各原子種の原子番号
natom 5 #原子数
typat 1 2 3 3 3 #各原子の種類を指定 (例だと、Sr Ti O O O となる)
xred #各原子の相対座標 または xcart #各原子の座標
0.0 0.0 0.0
0.5 0.5 0.5
0.5 0.5 0.0
0.5 0.0 0.5
0.0 0.5 0.5 
ixc 11 #交換相関項の設定。GGAは11。擬ポテンシャルもGGA用が必要
nsppol 2 #スピン分極を考慮する場合は2
spinat #各原子の初期のスピン(ディフォルト値は0.0d0、h-bar/2単位)
0.0 0.0 0.0 #typat で指定した順なので原子番号38。x, y, z のスピンは全て0
0.0 0.0 0.0 #原子番号22で、x, y, z の各軸の方向で "upの数 - down の数" の値を入力。
0.0 0.0 0.0 #原子番号8で、x, y, z のスピンは全て0
0.0 0.0 0.0 #原子番号8で、x, y, z のスピンは全て0
0.0 0.0 0.0 #原子番号8で、x, y, z のスピンは全て0
nband #考慮する状態(バンド)数
68 68
ecut 30.0 # 平面波展開のカットオフエネルギー(Hartree単位)
occopt 7 # 7はGaussian smearing
tsmear 0.01 #フェルミ端のブロードニング(Hartree単位)
kptopt 1 # 対称性を考慮してirreducible Bz内でのk点を作成
ngkpt 6 6 6 #Monkhorst-Pack法によるk点(x, y, z方向のk点数)
nstep 40 # 最大イタレーション
toldfe 1.0d-6 # 収束判定条件(total energy) または toldff #収束判定条件(atomic forces)
--------------------------------------------------

case.in (複数の計算を行う場合 : multi dataset mode)
この方法は、原子位置を徐々に変えて計算したり、SCF計算の後にDOSやバンド図を計算させたりすることを設定することとができる。下記はカットオフを徐々に変えた例。
--------------------------------------------------
ndtset 3 # データセットの数
ecut1 30.0
ecut2 35.0
ecut3 40.0
または下記のようにする
ndtset 3 
ecut : 30.0
ecut+ 5.0
--------------------------------------------------

計算の実行は、abinit < case.files > case.log

Q. ABINITでDOSやband図、Gamma点のフォノン、M点での不安定モードに対応する原子変位を計算をする方法を教えて下さい。
A. SCF計算で出力されたcase0_DENファイルを複製して、casei_DENという名前に変えて下さい。

case.in(DOSを得るための計算)
--------------------------------------------------
iscf -2 # SCF計算をしない
prtdos 3 # Local Density of States inside a sphere centered on an atomを出力
nband 28 # 考慮する状態(バンド)数
tolwfr 1.0e-10 # wf許容誤差
--------------------------------------------------

case.in (バンド図のための計算)
--------------------------------------------------
iscf -2 # SCF計算をしない
nband 28 # 計算するバンドの数
enunit 1 # eV単位の固有エネルギーとなるように設定
tolwfr 1.0e-10 # 2|nk>  <= 1.0e-10
kptopt -5 # 負の場合はkptboundsとndivkが対応する
kptbounds #
0.0 0.0 0.0 # Gamma点
0.5 0.0 0.0 # X点
0.5 0.5 0.0 # M点
0.0 0.0 0.0 # Gamma点
0.5 0.5 0.5 # R点
0.5 0.5 0.0 # M点
ndivk 15 15 20 25 15 # k点の分割
-------------------------------------------------- 

case.in (Gamma点に対するフォノンの計算)
--------------------------------------------------
ndtset 2 # データセットの数
# キーワードの直後に1が書かれていると最初に計算される。
kptopt1 1 # 対称性を考慮してirreducible Bz内でのk点を作成
tolvrs1 1.0d-10 # SCF停止の基準
# キーワードの直後に2が書かれていると二番目に計算される。
rfphon2 1 # フォノンに関する応答関数
rfatpol2 1 5 # 1から5の原子が移動する
rfdir2 1 0 0 # x軸方向
nqpt2 1 # Q点の数
qpt2 0 0 0 # Gamma点
getwfk2 -1 # GET the wavefunctions from _WFK file
kptopt2 2# 時間反転対称
tolvrs2 1.0d-8 # SCF停止の基準
#以下は共通。キーワードの直後に数字が書かれない。
acell 7.267 7.267 7.267 # セルの長さの設定
ntypat 3 # 原子の種類
znucl 38 22 8 # 各原子種の原子番号
natom 5 # 原子数
typat 1 2 3 3 3 # Sr Ti O O O
xred # 各原子の相対座標
0.0 0.0 0.0
0.5 0.5 0.5
0.5 0.5 0.0
0.5 0.0 0.5
0.0 0.5 0.5 
ecut 40.0 # 平面波展開のカットオフエネルギー
ngkpt 6 6 6 # Monkhorst-Pack法によるk点(x, y, z方向のk点数)
--------------------------------------------------

case.files (M点での不安定モードに対応する原子位置変位の計算)
--------------------------------------------------
case.in
case.out
casexo_DS3_DDB
foo
--------------------------------------------------
case.in (M点での不安定モードに対応する原子位置変位の計算)
--------------------------------------------------
eivec 1 # flag to turn on the analysis of phonon eigenvectors
nph1l 1 # number of q-points for phonon interpolation
qph1l 1 1 0 2 # list of q-points for phonon interpolation
--------------------------------------------------

計算の実行は、anaddb < case.files > case.log

Q. ABINITで特定の原子を固定した場合の構造最適化の方法を教えて下さい。
A. 下記のようにnatfixとiatfixを用いて下さい。

case.in (特定原子を固定した構造最適化)
--------------------------------------------------
ionmov 3 # 修正Broydenアルゴリズムを使用(全エネルギー考慮)
ntime 10 # Broydenタイムステップの最大値
tolmxf 5.0d-4 # 構造最適化での収束判定条件(Hartree / Bohr)
nsym 1 # 空間群 No.1
acell 4.73565 4.73565 47.25
angdeg 90 90 120
natfix 2 # 固定する原子種の数
iatfix 6 7 # 固定する原子の番号
ntypat 2 # 原子種の数
znucl 6 27 # 各原子種の種類を指定
natom 7 # 原子数
typat 1 1 2 2 2 2 2 # C C Co Co Co Co Co
xred # 各原子の相対座標 または xcart # 各原子の座標
0.0 0.0 0.4
0.3333 0.6666 0.4
0.3333 0.6666 0.488
0.0 0.0 0.56938
0.3333 0.6666 0.65076
0.0 0.0 0.73214
0.3333 0.6666 0.81352
nband 36 # 考慮する状態(バンド)数
ixc 7 # 交換相関項の設定。Perdew-Wang 92 function
ecut 30.0 # 平面波展開のカットオフエネルギー(Hartree単位)
kptopt 1 # 1にすると手動でk点を設定することになる
ngkpt 6 6 2 # セル内を6 x 6 x 2で分割する
nstep 40 # 最大イタレーション
toldfe 1.0d-6 # 収束判定条件(total energy) または toldff # 収束判定条件(atomic forces)
--------------------------------------------------

Q. ABINITでチュートリアルの入力ファイルを書き直して実際の計算をする場合、どのくらいの精度に設定し直せばよいですか?
A. XtalEditのsample/ABINITの中には、入力ファイルの例が書かれています。この程度の精度になるように、チュートリアルの値を書き換えると良いと思います。(原子種や原子座標の設定は除く)
Ga2O3 : ionmov 2, ntime 200, acell 10.0531007914 10.0531007914 10.0531007914, tsmear 0.04, ixc 7, kptop 1, ngkpt 4 4 4, ecut 40.0, nstep 100, toldff 5.0d-5, tolmxf 5.0d-4, diemac 12.0
La2O3 : ionmov 2, ntime 200, optcel 0, tsmear 0.04, ixc 7, kptopt 1, ngkpt 4 4 2, ecut 40.0, nstep 100, toldff 5.0d-5, tolmx 5.0d-4, diemac 12.0
LaAlO3(Lattice optimize) : ionmov 1, dtion 100.0, vis 50.0, acell 11.322570056 16.039447611 11.401013107, ecut 50.0, kptopt 0 #(use nkpt.), nkpt 1, nstep 100, toldfe 1.0d-6, diemac 12.0
LaGaO3 (1000 k) : ionmov 2, ntime 200, acell 10.3651470474 14.6831709313 10.4369566349, tsmear 0.04, ixc 7, kptopt 1, ngkpt 3 3 3, ecut 45.0, nstep 100, toldff 5.0d-5, tomlxf 5.0d-4, diemac 12.0
LaGaO3(super cell) : ionmov 2, ntime 200, optcel 0, acell 14.7093964442 14.6831710000, 14.7093964442, tsmear 0.04, ixc 7, kptopt 1, ngkpt 2 2 2, ecut 45.0, nstep 100, toldff 5.0d-5, tolmx 5.0d-4
SrO : ionmov 2, ntime 200, optcel 2, acell 6.8949883953 6.8949883953 6.8949883953, ixc 7, tsmear 0.04, kptopt 1, ngkpt 4 4 4, ecut 40, nstep 100, toldff 5.0d-5, tolmx 5.0d-4, diemac 12.0

 

Q. TDDFTの入力ファイルを教えてください。

A. チュートリアルには下記のようにあります。(見たところ、全てのnbandとboxcenter, nbdbufを消して、diemacとcommon dataの部分を計算したい系の情報に書き換えればよいようだ{調査中})

 case.in
  --------------------------------------------------
  ndtset 2 (SCFとTDDFTの二つのプロセスがあるという意味)
  # data set 1 SCF
  iscf1 5 # SCFの計算の設定
  tolwfr1 1.0d-15 #2|nk>  <= 1.0e-15
  nband1 5 # 考慮する状態(バンド)数
  prtden1 1 # TDDFTの計算で用いるための電荷密度分布のデータを出力する
  getwfk1 0 # no use of previously computed output wavefunction file appended with _DSx_WFK is done.
  #data set 2 TDDFT
  iscf1 -1 # TDDFTの計算の設定
  tolwfr2 1.0d-9 #2|nk>  <= 1.0e-9
  nband2 12 # 考慮する状態(バンド)数、非占有準位の数も入れる
  getden2 1 # data set "1" からの電荷密度分布のデータを取得する
  getwfk2 1 #  data set "1" からの output wavefunction file を得る
  #common data
  acell 6 2*5 Angstrom # 6 x 5 x 5 の大きさの単位胞を設定
  boxcenter 3*0.0d0 # 中心位置を設定(TDDFTでの計算のみに用いられる)。ディフォルトは 0.5 0.5 0.5
  diemac 1.0d0 # For silicon, use 12.0 . A similar value is likely to work well for other semiconductors.  For wider gap insulators, use 2.0 ... 4.0. ディフォルトは 106 (metallic damping).
  diemix 0.5d0 # model dielectric mixing factor.  0.0 から1.0の間の値。ディフォルトは 1.0
  ecut 25 # 平面波展開のカットオフエネルギー(Hartree単位)
  ixc 7 # 交換相関項の設定。Perdew-Wang 92 function
  kptopt 0 # k点設定のオプション(この場合、1点計算になる)
  natom 2 # 単位胞内に2つの原子を用いる
  nbdbuf 0 # bufferとしてのバンドの数
  nstep 25 # 最大イタレーション
  ntypat 1 # 1つの原子種のみ用いる
  znucl 7 # 窒素原子Nを1番に設定
  typat 1 1 # N と N の座標を単位胞内で設定する
  xcart -0.54885 0 0 0.54885 0 0 Angstrom # 座標の設定
  --------------------------------------------------

 

Q. GW近似での入力ファイルを教えてください。

A. チュートリアルには下記のようにあります。擬ポテンシャルはLDA用のものが必要です。最後にpspncやfhiと付いたファイルや、Troullier-Martines pspやHamann-type, LDA CA Perdewangとファイルに書き込まれたものを用いて下さい。

 case.files 
  --------------------------------------------------
  case.in
  case.out
  casei
  caseo
  caset../../Psps_for_tests/14si.pspnc
-------------------------------------------------- 

 

 case.in
--------------------------------------------------
nstset 3 #
# Definition of parameters for the calculation of the KSS file
nbandkss1 -1 # KSS fileにおけるバンド数に関する設定
nband1 9 # 占有と非占有のバンド数
# Calculation of the screening
optdriver2 3 # Screening calculation
getkss2 -1 # 前のデータセットからのKSS fileを取得
nband2 17 # screening 計算で用いるバンド数
ecutwfn2 2.1 # 波動関数のためのカットオフエネルギーを設定, ecutwfn is smaller than ecut, It would be more convenient to keep the default ecut value.
ecuteps2 3.6 # 誘電関数のためのカットオフエネルギーを設定, A value of ecuteps between 5 and 10 Hartree often leads to converged results
ppmfrq2 16.7 eV # Note that, if the plasmon-pole approximation is good, then, the choice of ppmfrq should have no influence on the final result.  One should check whether this is the case.  In general, the plasmon frequencies of bulk solids are of the order of 0.5 Hartree.
awtr2 0 # 時間反転対称性無しにAdler-Wiser expressionを用いる
#Calculation of the Self-Energy matrix elements
optdriver3 4 # Self-Energy calculation
getkss3 -2 # KSS file をデータセット1から得る
getscr3 -1 # SCR fileを前のデータセットから得る
nband3 30 # Self-Energy 計算で用いるバンド数
ecutwfn3 5.0 # It would be more convenient to keep the default ecut value.
ecutsigx3 6.0 # It would be more convenient to keep the default ecut value.
nkptgw3 1 # GW補正のためのk点の数
kptgw3 # k点の指定
-0.125 0.0 0.0
bdgw3 4 5 # 4から5のバンドに対してGW補正を計算
#Data common to the three different datasets
acell 3*10.217 # 10.217 x 10.217 x 10.217のユニットセル
rprim # primitive vectors
0.0 0.5 0.05
0.5 0.0 0.5
0.5 0.5 0.0
ntypat 1 #
znucl 14 #
#Definition of the atoms
natom 2 #
typat 1 1
xred 0.0 0.0 0.0 #
0.25 0.25 0.25
#Definition of the k-point grid
ngkpt 4 4 4 #
nshiftk 4 #
shiftk 0.5 0.5 0.5 #
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5
istwfk *1 # This is mandatory in all the GW steps.
#Use only symmorphic operations
symmorphi 0 #
# Definition of the planewave basis set
ecut 8.0
# Definition of the SCF procedure
nstep 10 #
toldfe 1.0d-6
diemac 12.0 #
iscf 5
--------------------------------------------------

上記は下記のようにして、ディフォルトの値を用いるようにしても、0.001 eVでGapの値が求められる。

# ppmfrq2 16.7 eV # Note that, if the plasmon-pole approximation is good, then, the choice of ppmfrq should have no influence on the final result.  One should check whether this is the case.  In general, the plasmon frequencies of bulk solids are of the order of 0.5 Hartree.

#wtr2 0 # 時間反転対称性無しにAdler-Wiser expressionを用いる

 

case.files 
--------------------------------------------------
  case.in
  case.out
  casei
  caseo
  caset
  ../../Psps_for_tests/13al.981214.fhi
-------------------------------------------------- 

 

case.in
--------------------------------------------------
# Create the KSS file
kssform 1 # ディフォルトは1
nbandkss 30 # the KSS file でのバンドの数。GW計算で用いられる。
nband 6 # バンド数の設定
prtwf 0 # 波動関数のアウトプットは無し
# definition of occupation numbers
occopt 3 # フェルミディラック関数を用いる
tsmear 0.05 # フェルミ端のブロードニング(Hartree単位)
   # Definition of the unit cell
acell 3*7.652 # 7.652 x 7.652 x 7.652のユニットセル
rprim # FCC primitive vectors
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
# Definition of the atom types
ntypat 1 # 1つの原子種のみ用いる
znucl 13 # 原子番号13、つまりAlを"1"に設定
# Definition of the atoms
natom 1 # 単位胞内で設定する原子の数
typat 1 # "1"( つまりAl )を単位胞内で用いる
xred 0.0 0.0 0.0 # "1"( つまりAl )の座標を設定
# Definition of the planewave basis set
ecut 8 # 平面波展開のカットオフエネルギー(Hartree単位)
# Definition of the k-point grid
ngkpt 4 4 4 # セル内を4 x 4 x 4で分割する
#64 k point
nshiftk 1 # k点のグリッドをシフトさせる数
shiftk 0.0 0.0 0.0 # k点のシフト。
istwfk 8*1 #
#Definition of the SCF procedure
nstep 50 # 最大イタレーション
toldfe 1.0d-8 # 収束判定条件(total energy)
prtvol 5 # 結果を出力するボリュームを制御する
enunit 1 # eV単位の固有エネルギーとなるように設定
--------------------------------------------------

 

Q. DFT+Uでの入力ファイルを教えてください。

A. チュートリアルには下記のようにあります。
擬ポテンシャルは paw が必要です。(ldaの記述があるものも必要か調査中)

case.files
--------------------------------------------------
case.in
case.out
casei
caseo
casex
../../Psps_for_tests/Ni.atompaw.lda.paw
../../Psps_for_tests/08o-lda.atompaw 
--------------------------------------------------

case.in
--------------------------------------------------
#Spin
nsppol 1
nspden 2
nspinor 1
spinat 0 0 1 # 各原子の初期のスピン
0 0 -1
0 0 0
0 0 0
#Parameters
nstep 50 # 最大イタレーション
ecut 15 # 平面波展開のカットオフエネルギー
pawecutdg 30
iscf 17
toldfe 1.0d-5 # 収束判定条件(total energy)
nband 40 #考慮する状態(バンド)数
occopt 1
#Structural parameters
natom 4 # 原子数
ntypat 2 # 原子種の数
typat 1 1 2 2 # Ni Ni O O
znucl 28 8 # 各原子種の原子番号
xred 0 0 0 # 各原子の相対座標
0.0 0.0 0.5
0.5 0.5 0.25
0.5 0.5 0.75
acell 3*7.92 # セルの長さの設定
rprim 0.0 1/2 1/2
1/2 0.0 1/2
1.0 1.0 0.0
# Kpoint Grid
ngkpt 2 2 2 # セル内を2 x 2 x 2で分割する
chksymbreak 0 # The k point grid is not symmetric, but the calculations being for the ground-state, this is not a problem.
# LDA+U
usepawu 1 # FLL法だと1, AMF法だと2
lpawu 2 -1 # znucl での原子の順で入力。-1: U無し、2:d軌道にU, 3: f軌道にUを考慮
upawu 8.0 0.0 eV
jpawu 0.8 0.0 eV # upawu の10%の値を入力
--------------------------------------------------

Q. 単位胞内の電子数を変えるには?
A. chagreを持ちいて下さい。
nelect = zion - charge になります。

下記も参照してださい。
http://www.abinit.org/doc/helpfiles/for-v7.10/input_variables/vargs.html#charge
「To treat a system missing one electron per unit cell, set charge=+1.」
=「単位胞あたりで1電子欠けた系を取り扱うため、charge=+1を設定する」
------------------------------------------------------------------------------


------------------------------------------------------------------------------
■ WIEN2kでの入力ファイルのQ&A
http://www.wien2k.at/reg_user/faq/ も見てみると良い。

Q. 遠隔操作をする方法を教えて下さい。
A. Google Chron や Firefox、sufari の URLに WIEN2kがインストールされている計算機のIPアドレス:7890 と入力して下さい。例えば、XXX.Y.ZZ.AAA:7890 などです。これで上手く動作しない場合は WIEN2kがインストールされている計算機のTerminal でw2web を入力して得られる最後の4つの数字を 7890と取り替えて下さい。各種のデータの取得は、[Files>>]を開き、[show all files]で、得たいファイルの左側にあるフロッピーディスクのアイコンをクリックすると可能になります。余談ですが、RunSCFのところで、e-mailの行にチェックを入れて、e-mailアドレスを書き込むと計算が終了した時に連絡してくれるようになります。(どうしても動作しない場合は、WIEN2kがインストールされている計算機側からgoogleの検索画面などが見れるか調べて下さい)

Q. WIEN2kでcore electrons leakと出力されます。どうすれば良いでしょうか?
A. Separation energyの値を変更します。ここでは価電子帯(Valence)と準内殻(semi-core)を分けるエネルギー位置を入力することになります。余談ですが、準内殻(semi-core)と内殻(core)を分けるのは、先に設定した値から- 3Ryした値になります。

Q. WIEN2kでのRmtはどのように設定すれば良いでしょうか?
A. 隣接原子と重ならない限り、出来るだけ大きくとる。しかし、小さな値にしても精度的に論文で非難されるようなことはありませんが、計算時間が増大します。WIEN2k HPのPAPERSから値を参考にするのも手でしょう。ちなみに、Rmtの単位はa.u.です。(Ca : 2.0, Cu : 1.9, Ti : 1.9, Ge : 1.9, Sn : 2.2, V : 1.9, O : 1.6程度の値だったりします) 
  上記までで問題がある場合は、入力ファイルにおいて、結晶群、格子定数と角度、原子種、座標のみ入力し、あとはWIEN2k側で入力させるようにすると解決する場合があります。cifファイルを使って入力ファイルを作るのも手です。
※ Total DOSは変わりませんが、Partial DOSの大きさは、RMTで変わります。RMTを大きくすると、大きくした原子に対するPartial DOSが大きくなります。また、積算した電子数も大きくなります。

Q. WIEN2kでクーロンや交換相関エネルギーを得るにはどうすればよいの?
A. WIEN2k GUI (web画面)で左側にある[Files>>]で[SCF files]を選択します。case.scfのファイルを見ると、下記のデータが書かれています。(case=ファイル名)
 ENE : The total Energy.
 SUM : The band-structure energy including the lowlying semicore state.
 DEN : The potential density integrals are contains E-coul, Exc and kinetic E-correction.
となり、内殻のエネルギー core energy = ENE - SUM - DEN となります。
WIEN2k_10では、case.in0でKXCを用いると、
 DEN : coulomb, kinetic and exchange-correlation energy
を分離して得ることができます。

Q. バンド分散用のk点やCIFファイルの利用、InterstitialのDOS表示、フェルミ面の可視化はどうすればいいの?
A. 下記のアフィニティサイエンス WIEN2k FAQをご参照ください。
http://www.affinity-science.com/wien2k/faq.html 

Q. バンド分散を描くときに任意のk点でプロットする方法を教えて下さい。
A. [Tasks >>] → [Bandstructure] → Generate k-mesh using Xcrysden を押して、得たいk点を次々と押していきます。Total number of k-points along the path: の右の空欄に選択したk点の初めから終わりまでを分割する数(全体のk点数)を入力して(k点数は入力した値に近い点数で最適なものに自動で書き換えられます)、ファイル名をxcrysden.klistとしてください(Generate k-mesh using Xcrysden の右に書かれているものをコピー&ペーストするのが便利です)。create wien.klist_band の左を from xcrysden にして、create wien.klist_band を押します。再度 Bandstructureの画面に入り、create wien.klist_band の左を from xcrysden にしてからは、いつも通りにcreate wien.klist_band のある行の次から、上から順にボタンを(濃い灰色の部分は除いて)押していきます。(つまり、[Tasks >>] → [Bandstructure] → Generate k-mesh using Xcrysdenを押したあとは、http://www.affinity-science.com/wien2k/faq.html の3.から始めればよい)

Q. LDA+U またはGGA+Uの計算をしたいのですが、どのようにすればよいか教えてください。
A. 下記の文献[L1-L3]を見て下さい。3dを含んだ化合物に対して、3dに対してUを含んだ計算してみると良いでしょう。例えば文献[L4-L7]にあるような系で計算してみて、Uを導入することでギャップが大きくなるかを調べて、LDA+Uの計算に慣れてみるのも良いでしょう。
spを入れて非磁性の計算にし、SOとLDA+Uの計算をする場合のコマンド: runsp_c_lapw -so -orb
  runsp_c calculates only spin-up vectors and copies them to spin-dn. Thus spin up and dn are exactly the same and any spin-moment is gone. (P.Blaha)
  run_lapw -orb was NEVER possible (as far as I remember). You must use a spin-polarized setup and use runsp_c_lapw -orb (P.Blaha)
  In order to speed up the runsp_c calculations, in the new version lapw2 -dn is not executed, but the necessary files are copied from the corresponding "up" files.  Since lapwdm is usually so fast, this was not done for this step. However, in parallel mode it would not work this way.  Therefore one has to make the same "trick" as with lapw2: comment the spin-dn execution and add two "cp" lines:
total_exec lapwdm -up $para
# total_exec lapwdm -dn $para
cp $file.scfdmup $file.scfdmdn
cp $file.dmatup $file.dmatdn
(P.Blaha)
  Please note, we do not recommend to use -orbdu .  Did you run initso_lapw (with spin-polarization !!) and accepted eventually the newly generated struct file ???  SO + spin-plo. may lower symmetry. (The "ipass" message indicates to me that you still have very high symmetry) (P. Blaha)
[L1] http://www.apph.tohoku.ac.jp/kiso/files/0803_magne.pdf 
[L2] http://aperonya.blog115.fc2.com/blog-category-3.html
[L3] http://www.nims.go.jp/cmsc/staff/arai/wien/ldau.html
[L4] http://www.mizuho-ir.co.jp/solution/research/semiconductor/nano/titania2/index.html
[L5] http://www.fplo.de/workshop/ws2008/koepernik2_fplo2008.pdf
[L6] http://www.wien2k.at/reg_user/textbooks/WIEN2k_lecture-notes_2011/Exercises_11.pdf での Exercise 6: を試してみる。
[L7] http://www.wien2k.at/reg_user/textbooks/Constraint_U.pdf
[L8] http://www.wien2k.at/reg_user/textbooks/WIEN2k_lecture-notes_2013/ldau.pdf 
[L9] http://cometscome.blog76.fc2.com/blog-entry-1009.html
※ Memo(他者から収束したファイルを貰った場合)
1. Expert のモードで select spin-polarizedにチェックを入れて、下にある「CHECK BATCH VALUES」を押す
2. initso_lapw で x symmetso -c から始める

Q. LDA+U での計算が上手く収束しません。どうすればよいでしょうか?
A. LDA+Uは計算が不安定になりやすくなります。LDAでSCF計算を行った後に、Uを導入した計算をすると良いでしょう。それでも駄目なときは case.in1のLAPWのエネルギー値を変えます。
※ 幾つか計算した後に、k点数を変えたい場合は、[Utiles] -> [initso_lapw]を押して、[kgen] や [kgen -so]でk点数を変えて下さい。
※ Uを入れた計算はHF近似の計算も入ることになりますから、複数の解が得られることもあります。経験的計算にもなりますので、物理的な考察も大切になります。浜田典昭「固体物理 39」(2004) 743. は読んでおくと良いでしょう。
1. SCF calculation.
2. Copy from case.inorb and case.indm or case.indmc in SRC_templates to them in calculation directory.
3. change from case to calculation directory name.
4. rewrite them.
5. check a box for orbital pot(LDA+U), It-number 1 (1 shot calculation case) and  spinorbit (+SO case), then start SCF cycle.


Q. WIEN2kで電場や磁場を掛けた計算方法を教えて下さい。
A. zig-zag potentialの概念図は[10]を参考にするとよいでしょう。電場の場合は、case.in0の4行目を書き換えます。磁場の場合はcase.inorbを書き換えます。(編集中)
J. Stahn et al., Phys. Rev. B 63 (2001) 165205.
 : http://prb.aps.org/abstract/PRB/v63/i16/e165205 (電場)

Q. Spin orbital coupling での計算をしたいのですが、case.inso での NX のある行はどのようにすればよいのですか?
A. 通常はNX がある最初の行のNX を 0 にて、次の行を消してください。もし、semi-core の範囲(ディフォルトでは-6 Ry から -3 Ry)に p1/2の軌道がある場合は、最初の行のNXへそれに対応する原子の番号(StructGen や DOSと同じ)の数を入力して、NXの数だけ NX1の行をコピー&ペーストして、NX1に原子の番号(StructGen やDOS)を入力します。

Q. WIEN2kを用いてラシュバ効果の計算をしたいのですが、どのように設定すればよいでしょうか?
A. spinを入れた計算でもvalence bandではscalar relativesticの計算となるので、スラブモデル(真空層を入れたモデル)で計算し、upとdown spinでバンド構造が異なるかを調べてみてください。spin orbit couplingを入れて計算をしたい場合は、[run SCF]→Optionsでのspinorbitにチェックを入れます。(semi-core{ディフォルトでは-6 Ry から -3Ry}の範囲にp1/2があり、相対論用にp1/2を入れた結果を議論したい場合{結果がより実験値に近づく場合もあるかもしれない}は、case.insoでNXのところに計算に入れる数を入力します。次に、NX1に原子に付けられた番号を入れます(この行の部分を、NXで指定した数だけ増やします)[11-13])

Q. WIEN2kを用いてEELSの計算をしたいのですが、特定の方位での遷移確率を調べるにはどうすればよいでしょうか?
A. TELNESをExpert modeにします。InnesGenTM for TLNESがでますので、Orientation sensitive : のところにチェックを入れて、計算させたい方位を指定すればよいです。case.incで計算させたい原子の内殻をOCCUPを2から1にして内殻空孔を導入します。そして、case.in2のNEの値を一つ増やします。case.in1で1.5を2.5や3にします(非占有準位の計算範囲を増やす)。上記は[Remove 1 core electron, add 1 electron to conduction band]での計算となる。さらに、TELNESで x lapw1 -c の計算の後に、case.in2のNEの値を元に戻して、続けて計算する方法もある[14]。この他に、内殻空孔を入れない方法(sudden approximaiton)や、原子番号がZ+1である原子に入れ替えて計算する方法(Z+1 approximation)などもあります[15]。
Expert modeのさせ方:WIEN2kのGUIで左の欄にある[Configuration]を押してELNES expert mode にチェックを入れ、Save these values を押します。WIEN2kのソースファイルを解凍したフォルダで、SRC_lapw2のフォルダにあるmodules.F をkwriteやemacs、vi などで開いて LXDOS = 1 を LXDOS = 3 にして save します。./siteconfig_lapw で r を押して、s を選択し、lapw2 と入力します。

Q. スーパーセルで元素置換を行った系の計算をしたいのですが上手くいきません。
A.1.  文献[9]を参考にスーパーセルを作成して下さい。その後、置換する元素が、置換する前の元素よりも大きなZ番号を持つ場合は、格子定数を大きめに設定して計算して下さい(例えば、5%程度)。計算が収束したら、構造最適化を行って、最適な格子定数を決定します。

A.2.  下記に手順を示します。
  1) [ single prog. ] → [ Other programs ] → [ supercell ] にチェックを入れ[ Ececute !]を押す。
  2) Number of cells in x direction : に x 方向に拡張したい数を入力
    y, z も同様に行い、 一番下にある [Ececute !] を押す。
    ( [ Enter your target lattice type: ] も選択した方が良いかもしれない)
  3) case_super.strcut という名称のファイルができる。このファイル以外を削除し、case_super.strcut を case.struct という名称に変える。
  4) [StructGenTM]にて、置換したい部分の原子名とZを変える。
  5) [StructGenTM] → [set automatically RMT and continue editing] → [Resude RMTs by %] → [ do it ]
  6)  [ Save structure ] → [ save file and clean up ]
  7) 通常と同じように進めていき([ initialize calc. ]にて上から順にボタンを押していく)、[ Use new struct-file ? ] では No, [ Use strcut-file generated by sgroup ? ] で Yes を選択し、x symmetry に進む。
  8) 後は通常通り進めて行けばよい。
A.3. スーパーセルの構造を計算して、VESTAでP1にて構造を入力して、確認し、cif データとして書き出して、WIEN2kに読み取らせることもできる。この場合、3) の手順から進めればよい。


Q. Use struct-file generated by sgroup? (Usually No, unless WARNING appeared above)ではどのようにすればよいでしょうか?
A. ここでYesを選択すると、WIEN2k側でより対称性を良くした構造の入力ファイルを作成してくれます。しかし、適切な構造の入力ファイルとならないこともあるので、通常はNoを選択して下さい。(スーパーセルを用いる場合などでは、Yesを選択して、XCrySDenなどで構造を確認して問題が無ければ、計算を走らせてみることもお勧めします。上手く入力ファイルが作成されれば、計算に掛かる時間をかなり短くできます。この場合、スーパーセルで置換した系を計算するのに格子定数を実験値をより大きめに設定しなくてもエラー無く収束した例が報告されています)

Q. 構造最適化が上手く動きません。
A. [StructGenTM] → [set automatically RMT and continue editing] → [Resude RMTs by %]の部分を1から5程度の値を入れて計算させて下さい。optimizeで For option 1-4 の下の欄に体積変化率を入力します。edit optimize.job で必要とするコマンドの#を消します。磁性の有無や反転対称性の有無などです。最適化が終われば、E vs. volume から、lattice parameter がbohrとAng(オングストローム)で表示されます。

Q. 原子位置の最適化を行いたいのですがどうすればよいですか?
A. RMTを小さくし、k点を1にして、[mini.positions] を押して下さい。[Prepare commandline] のボタンを押します。[start it ...]のボタンを押します。
  「これでも時間が掛かりそうなときは、実験などから真の値に出来るだけ近い構造を初期値として入力し、SCF計算します。DOSとバンド図などを得ておきます。そして、新しくディレクトリを作り直し、RmtKmaxを5〜6の値にしたり(Cutoff =13.6 * (RmtKmax / Rmt ) ^2 [eV]となるので必要なカットオフのRmtKmaxにする )、x lstart での core と valence band の区切りを -6.0 Ry から -3.0 などにします。このとき x lstart に戻されますが、無理やり次のボタンに移ってください。前に得られた結果を吟味しながら、case.in1 での El やエネルギー刻みを設定します。case.in2 の最後にあるGmaxも変えます(GMAXは電荷ρとポテンシャルvを表現するのにどの程度まで精度が欲しいかに関連する。http://arxiv.org/pdf/1103.4466.pdfに詳細な式がある)。そして、k点を1にして、[mini.positions]で計算していきます。構造最適化の後にk点を大きくして、構造最適化前と後を比べます。物理的に意味のある結果になっているかしっかりと他の実験結果と吟味します」(←この記述は確認中)

Q. SCF計算を途中で止めるにはどうすればよいですか?
A. [Utils.>>]→[stop SCF]か[stop mini]で止めることが可能です。再度計算させるには、[Execution >>]→[run SCF]→[Remove the files case.broyd[1][2]]に対応する部分を選択して、再度[run SCF]を押して再び計算を始めることができます。計算の経過は[Utils.>>]→[show dayfile]で知ることができます。

Q. DOSやバンド図を描けないのですが、どこか間違っているのでしょうか?
A. 濃い灰色の部分にあるボタンを押さずにDOSやbandstructureなどの計算をしてみて下さい。ただし、bandstructure で各原子の軌道を計算したい場合は濃い灰色の部分の上から二番目のボタン[x lapw2 -band -qtl -p]を押して下さい。そして、edit case.insp でjatomとjtypeのある行の左から一番目と二番目をDOSのときと同じように表示させたい原子と軌道の番号を書き入れます。

Q. DOS を計算したときに'NaN'と表示されます。どうすればよいでしょうか?
A. 下記のHPに解決方法が記載されています。
http://zeus.theochem.tuwien.ac.at/pipermail/wien/2011-September/015419.html

Q. DOSの表示範囲を変えたいのですが、どのようにすればよいでしょうか?
A. case.in1_st または case.in1 でのEMAXの部分を大きくして run_SCF をして下さい。続いて、DOSの edit.int でもEMAXの部分があります。

Q. DOSの単位は何ですか?
A. case.structに記載されている原子の数だけ計算されます。代表的な原子座標を記述すると等価な原子が自動的に記載されていることが分かります。その数の合計と各原子でのDOSが一致します。そのため、case.strcutに記載されている等価な原子の数でDOSを割ることで、単位は[Density of states / eV. atom] になります。同様に、Total DOS でもcase.structに記載されている原子数で割り[Density of states / eV. atom]にします。[Density of states / eV. unitcell]ではないので注意してください。
 Total DOS = Σ Partial DOS + interstitial DOS (case.struct の Inewuivalent Atoms と edit case.int のAtomが対応)
 Total DOS = Total DOS (up) + Total DOS (down)
※ 単位胞内での原子数が分かっている場合は、Total DOSが下記のようにして出力されます。
1) BCC: 全原子*1/2個 (case.structで記述されている原子の数と同じになることが分かる)
 得られた Total DOS を 2 倍すると [Density of states / eV. unitcell] になる。
2) FCC: 全原子*1/4個 (case.structで記述されている原子の数と同じになることが分かる)
 得られた Total DOS を 4 倍すると [Density of states / eV. unitcell] になる。
3) HCP: 全原子(case.structで記述されている原子の数と同じになることが分かる)
 得られた Total DOS を 1 倍すると [Density of states / eV. unitcell] になる。
4) SC: 全原子(case.structで記述されている原子の数と同じになることが分かる)
 得られた Total DOS を 1 倍すると [Density of states / eV. unitcell] になる。
※ case.dos1(Ryでの出力。integrated DOSも出力される)
 Ry → eVへの変換: (Energy - EF)*13.6058, DOS / 13.6058 とする。
 EF値は、case.dos1の二行目の一番最初に記載されている。
※ case.dos1eV(eVでの出力。EFが0の位置に自動的に調整されている)

Q. バンド図を描く場合にエラーがでます。
A. Insert correct EFのところに、xxxxとある部分があります。そこに、右上に書かれている数値を入力してください。データはspaghetti_eneとして出力されます。OpenOfficeなどのspreadsheetでスペースなどで区切って、最後の2列(最初の行が空白の場合はそれを削除する)をigorなどのプロットソフトにコピー&ペーストして下さい。全く別の話ですが、DOSを計算した場合の結果はdos1eVとして出力されます(dos2eVやdos3eVが存在する場合はそこにもデータが存在しています)。ちなみに、PDOSの単位はDOS / eV.atomです。

Q. バンド図を描いた後に、DOSやSCF計算をするにはどうすればよいですか?
A. [Files >>]→[def files]→[lapw1.def]で、case.klist_band の部分を case.klist に書き換えます。

Q. 不要になったセッションを削除したいのですが。
A. [Session Mgmt. >>]→[delete session]で削除してください。また、データが格納されているファイルも残っているので、それも削除してください。余談ですが、[change session]で他のセッションに移動できます。

Q. 単位胞内の原子数を知りたいのですが。
A. case.struct の MULTI の総和を求めてください。その総和に係数を掛けます。case.structの二行目の最初に PやF, B, CXY, CYZ, CXZ, R, H などの記述があると思います。係数はP=1, F=4, B=2, CXY=1(調査中), CYZ=1(調査中), CXZ=1(調査中), R=3(調査中), H=1(調査中)となります。文献[17]に他のコードでの記載があります。

Q. OpenMP での並列計算ができないのですが……。
A. Linux 上でcd $HOME と入力して、 [view] → [Show Hidden Files] を選択すると .bashrc というファイルが見つかります。export OMP_NUM_THREADS=1 で数値を並列計算させたいスレッド数に変えてください。下記のk点並列やMPI並列の設定を各ディレクトリで設定するのが面倒な場合はOSが認識している最大のCPU数を設定しておけば良いでしょう。(慣れてきたらk点並列などを試してみましょう)

Q. 並列計算ができないのですが……。
A. 並列計算には、k点並列とMPI並列の二つがあります。k点並列の方が計算が早く済むことが多いですから、k点並列を用いることをお勧めします。
[Utiles>>] → [ edit .machines ] の最後を書き換えます。自動的にWIEN2k側である程度入力ファイルが作られます。それを参考に、k点並列とMPI並列で、下記のように書き換えます(k点並列の場合は特に大幅な変更はないでしょう)。
※ 簡単な系で、下記にあるk点並列の方法で、どれが計算速度が速いかを調べてみてください。
◇ OpenMP並列(4 core の場合)
.bashrc で  export OMP_NUM_THREADS=4 にする。
以上の設定をして、run SCF で計算します。
◇ k点並列計算の場合(4 core の場合)(最初の数値がk点を割り振る比率)
  1:localhost
  1:localhost
  1:localhost
  granularity:1
  extrafine:1
※ .bashrc で  export OMP_NUM_THREADS=1 か  export OMP_NUM_THREADS=2 にする。(未確認)
以上の設定をして、run SCF でのOptions でparallel1にチェックを入れます。
◇ k点並列計算の場合(4 core の場合)(最初の数値がk点を割り振る比率)(未確認)
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  granularity:1
  extrafine:1
※ .bashrc で  export OMP_NUM_THREADS=1 か  export OMP_NUM_THREADS=2 にする。
以上の設定をして、run SCF でのOptions でparallel1にチェックを入れます。
◇ k点並列計算の場合(4 core 8 スレッド の場合)(最初の数値がk点を割り振る比率)(未確認)
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  granularity:1
  extrafine:1
※ .bashrc で  export OMP_NUM_THREADS=1 とする。
以上の設定をして、run SCF でのOptions でparallel1にチェックを入れます。
◇ k点並列計算の場合(6 core 12 スレッド の場合)(最初の数値がk点を割り振る比率)(未確認)
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  1:localhost
  granularity:1
  extrafine:1
※ .bashrc で  export OMP_NUM_THREADS=1 とする。
以上の設定をして、run SCF でのOptions でparallel1にチェックを入れます。
◇ MPI並列計算の場合(4 core の場合)
  1: localhost:4
  lapw0: localhost:4
  granularity:1
  extrafine:1
以上の設定をして、run SCF でのOptions でparallel1にチェックを入れます。
文献:  http://www.hpc.co.jp/assets/files/Wien2kManual.pdf
※ BTOなどで購入したPCを複数用いて並列化させる場合は、NFSの設定が必要になります。
※ スレッド:参考文献:日本分光会、『分光測定入門シリーズ3 分光装置Q&A』、講談社サイエンティフィック

Q. 他のPCとでの並列計算ができないのですが……。
A. NFSの設定が必要です。そして、intel fortran でのコンパイルオプションで -i-static を忘れないように入れておいて下さい。

Q. XMCDの入力ファイルにおいて、WIEN2kのマニュアルと論文[16]での記述が異なっているのですが、どちらが正しいのですか?
A. Blaha教授に聞いたところ、論文に書かれている入力ファイルの方が正しいとのことです。コンパイルに不安がある場合は、論文と同じ結果が得られるかどうか試してみるとよいでしょう。

Q. WIEN2kのワークショップに参加できなかったのですが、独学する方法はありますか?
A. 下記のHPにおけるExercises_11.pdfを参照してみて下さい。
WIEN2k lecture note: http://www.wien2k.at/reg_user/textbooks/WIEN2k_lecture-notes_2011/
また、左欄にあるWIEN2kの中頃から後半を参照してみてください。

Q. 初期構造をLow spin配置にして計算する方法を教えてください。
A. Initial calc. での instgen_lapw にて、selected blowの欄にチェックを入れます。各原子に記入されているuをnに変えます。こうすると、case.inst にて、up と downが同じ数だけ配置されるようになります。もし、これで計算してもLow spinの状態にならないときには、その物質がhigh spin で安定であると考えられます。

Q. case.in1にあるEl や de の入力について説明してある Userguide にある Fig.7.1 はどのように解釈すれば良いのでしょうか?
A. 図中のul(r,El)はradial function(動径分布関数)です。例として、各軌道が混成し、軌道混成を形成すると、図7の右のDOSのような結果が得られます。エネルギーの高い反結合軌道エネルギー(Etop)をはじめに探して、少しずつ下に向かい、結合軌道のエネルギー(Ebottom)を見つけます。EtopとEbottomが見つかれば、その平均をとって、その混成軌道に前よりもより適した動径分布関数のエネルギー El を見つけます。
  WIEN2kではRmt を使いましたね。平面波だけで計算すると多くの平面波が必要になり、さらに時間も掛かるので、Rmt の距離よりも内側では内殻用の軌道(例えば、sやp軌道などの動径関数を用います)を用いることで(Rmtの外側で使われる)平面波を減らして計算時間を短縮します。Fig. 7の左でRmtまでが記述されているのは、動径分布関数がそこで終わるからです。その先は平面波で滑らかに繋がるようにしています。
 EtopとEbottomを探すエネルギー刻みが 0.000 の場合には、ElがEF-0.2 または EF+0.2 した値に設定されます。
 ゴーストバンドを避けたい場合には、LAPWではなくAPW(局在化している3dや4f電子だけをAPW)にしたりすることも考慮してみてください(ただし、時間が掛かるようになります)。
※ LAWPを使うときは RKmax=6〜10, APWを使うときはRKmax=5〜9でよい。[x1]
※ loとLOともにローカルオービタル。
※ APW+loでは interstital region で APW -> PW に、LAPW+LOでは interstital region でLOが消える。APW+loはLAPW+loと併用して使われたりする。
※ APW+lo: ゴーストバンド出ない(と言われる)。エネルギー的に非線形の基底関数になるので永年方程式を解くのが大変。
※ エネルギー的に線形化するようにしたのがLAPW+LO(よく使われる)。もし、ゴーストバンドが出現していたら、バンド分散を描くと直線になっており、さらに、ゴーストバンドに電子が詰められるので、本来よりもEFの位置が下がる。それでゴーストバンドが発生しているかを判断してみるとよい。正しいエネルギー位置を設定すると問題が解消される。
マフィンティンの内側では、これら双方とも局所的な球面関数で基本的には同じ手法で展開される。各原子での角運動量においては、動径部分が一つかいくつかの動径部分で展開される。

Q. case.inc とcase.in1 の違いや詳細を教えて下さい。
A. case.inc は 内殻の軌道に関する情報を記述します。case.in1 はRmt内でcoreよりも上のエネルギー、つまり、semi-core や valenceに関係する軌道を記述します。case.in1の最初の行ではEFの値がRy単位で記述されています。例えば、Tiの3sが-3.00と記述されていたら、-3.0 - EF = -3.0 - 0.5 = -3.5 Ry = -3.5*13.6 eV = -47.6 eV に Ti 3s の状態がDOSとして現われると初期では考えることになります。SCF計算のサイクル中に、上の質問に書いてある通り、EtopとEbottomを探して、その平均を Elにします。case.in1のEFは各サイクルで得られたEFの値が入れられます。計算値や実験値を参考にしてEFとElを近い値にしておくと収束が速くなる場合もあります。

Q. APW+lo の設定方法を教えてください。
A. case.in1 にあるCONTの右側の数値が1であればAPW(0だとLAPW)です。(L)APWでのloは第2項の値が負になっていることが多いです。ディフォルトの場合、APW+loになっていると思います。case.in1を確認してみて下さい。

Q. LOPW('LOPW' - Plane waves exhausted)のエラーの対処方法を教えて下さい
A. RKmaxを増やしてください。それでもだめな場合は、下記のように.in1でdバンドを除いてみます。

----------------------------------
WFFIL  EF= 0.50000       (WFFIL, WFPRI, ENFIL, SUPWF)
  7.00       10    4 (R-MT*K-MAX; MAX L IN WF, V-NMT
  0.30    4  0      (GLOBAL E-PARAMETER WITH n OTHER CHOICES, global APW/LAPW)
2   -1.83      0.002 CONT 0
2    0.30      0.000 CONT 0
0    0.30      0.000 CONT 1
1    0.30      0.000 CONT 1
----------------------------------
を下記のようにする。
----------------------------------
WFFIL  EF= 0.50000       (WFFIL, WFPRI, ENFIL, SUPWF)
7.00       10    4 (R-MT*K-MAX; MAX L IN WF, V-NMT
0.30    3  0      (GLOBAL E-PARAMETER WITH n OTHER CHOICES, global APW/LAPW)
2   -1.83      0.002 CONT 0
0    0.30      0.000 CONT 1
1    0.30      0.000 CONT 1
----------------------------------
http://www.mail-archive.com/wien%40zeus.theochem.tuwien.ac.at/msg03910.html

Q. ゴーストバンド( 'l2main' - QTL-B.GT.15., Ghostbands, check scf files )のエラーの対処方法を教えてください.
A. http://www.wien2k.at/reg_user/faq/qtlb.html に対処方法が記載されています。
※ 要点としては、case.in1 または case.in1c で Eパラメーターの値(ディフォルト 0.30 [Ry])をEFの値から0.2ほど小さい値に変えて再度計算させます。それでも駄目であれば、Eパラメーターを手作業で少しずつ変えていきます(それはエラーの生じたところから行っていけばよいと思います{私の経験})。LOを入れた部分でエラーが出ている場合は、どちらかの軌道の EパラメーターをEFに近づけるか、LOを消去してしまいます。
1. case.scf2 ファイルを見る。
  下記に似た表示があることが分かる。
:FER  : F E R M I - ENERGY(TETRAH.M.)=   1.01765
EFは 1.01765 Ry となったことが分かる。そして、
     QTL-B VALUE .EQ.  120.91788   in Band of energy    1.12346   ATOM=    1   L=  2
ここで、energy 1.12345 Ryで、ATOMが1番、軌道l=2で問題が生じていることを示している。
2. case.scf1 ファイルを見る。
  :E2_0001: E( 2)=    0.3000   E(BOTTOM)=   -0.066   E(TOP)= -200.000
などの記述がある部分ところを見る。これは、この軌道のE(BOTTOM)の軌道は見つかっているが、E(TOP)が見つかっていないことを示している(-200.000で計算をストップした)。
3.  case.in1 を修正する。
  case.scf2 で記述されいたATOMの設定を見ると下記に似た形になっている。1
0.30    6  0      (GLOBAL E-PARAMETER WITH n OTHER CHOICES, global APW/LAPW)
0    0.30      0.000 CONT 1
(省略)
2    0.80      0.000 CONT 1
  GLOBAL E-PARAMETER の行以降は次の順番で記述されている。
軌道l、RmtでのAPW/LAPWの軌道lのエネルギー(=El)、Rmtでの最適なエネルギーを探すためのエネルギー刻み (=de)、Rmtでの最適なエネルギーを探すのに必要なEtop(軌道lが関連するバンドの上端)かEbottom(バンドの下端)が見つからなくても続けるか?、1=APW か 0=LAPWか?
  対応するATOMと軌道lに対する Rmtでのエネルギーを EF に近い値にする。つまり、下記のようにする。
 E( 2) を正確な値に設定する(FEから0.2を引いた値にするか、試行錯誤する)。
 他の全ての軌道をディフォルトの 0.30 から、FEの値 から0.2 を引いた値の0.8などにしてEFの値に近づける。(※ 問題が生じた軌道の行は、QTL-Bの行で記述されている値 に近い値にする)
・筆者の経験
  ATOM=4 L=1 でエラーが出ており、case.in1 では下記のように書き換えた。
1    0.30      0.000 CONT 1 → 1    0.50      0.001 CONT 1
  最初のcycleで発生し、初期値が EF=0.5 なので、ATOM 4 の l=1軌道のエネルギーをEF近傍の値となる(というよりもほぼEFと同じ値) 0.5 に設定した。最適なAPW/LAPWのRMTでのエネルギー(El)を探すためエネルギーステップを 0.001 と小さい値にした。これらのパラメーターを書き換えた後は、RKmax 5.0から7.0までの小さい値にし、k点を1や8にして、cycleが走るかを調べ、上手く動作したら本格的にRKmaxやk点を大きくしていけば良い。上手くいかなければ、 El と de を少しずつ変えて試行錯誤する。
※ 上記の他、私の経験では、ゴーストバンドが表示された場合、 例えば、”QTL-B VALUE .EQ.  120.91788   in Band of energy    -2.12346   ATOM=    9   L=  2”(上記のゴーストバンドのエラーの対処方法を参照)と表示された場合に、case.in1または case.in1c でそれに対応する原子の軌道のEパラメーターに Band of energy に近い値を指定すると上手く計算が収束してくれたことがある。例えば、case.in1c で "2    0.30      0.005 CONT 1" ->  "2    -2.00      0.005 CONT 1" と変更したりする。それに加えて、単位胞の中で原子番号が似たものは全て同じ様に値を変更してみるとよい。
※ WIEN2k_13, ZrCoBiではCo l=0, 1でゴーストバンドのエラーが出た。対処方法は、case.in1cでEF=0.6にしてゴーストバンドに対処し、全エネルギーの収束が良くなってから発散している様子が見えていたのでmixingを0.05にしてcycle毎の変化を小さくした。case.in1cでCo l=0, 1のエネルギーを0.3から変えてもいた。しかし、どのような条件で上手くいったかは忘れてしまった。mixingとcase.in1cのパラメーターを変えた直後はinitialize calc.は2回以上行っている。

Q. lapw.error Cholesky INFO が生じるときとは?
A. case.struct(二つの原子が同じ位置を指定している)、case.in1(energy parameter の設定)が間違っているか、EKMAXが合理的ではないほどに大きい、などの理由が考えられます。

Q. 計算負荷を小さくしたい場合にはどこを変えればよいでしょうか?
A. 下記の点を考察してみて下さい。non-spinで計算して良い場合にはspin polarization を指定しない。プロだとここまで変えても正しい値が得られる。ただし、初心者が行うと収束が上手くいかず、計算時間が長くなるだけでなく、正しい結果が得られない。プロは十分に理解して書き換え、計算負荷を小さくして正しい値を得ることが出来ていると理解しておけば良いだろう。
※ 虚数を平面波展開での係数に含む場合には、入力ファイルの最後に c が付く。
※ 経験では RKmaxが小さいと、バンドが広がって見える。
※ 下記の条件で計算を行って、RKmaxを増加させたり、k点数を多くしたりしていくと良い。
■ 初心者
a) x lstart (case.in0)
  x lstart で LSDA を選択。GGA-PBEでも良い。
b) case.in1
  RKmax を 6.5 や 6.0 にする。どうしてもメモリが足りないなどの場合は、5.00±0.04に近い値にする。
 価電子帯(-15 eV 程度か?)局在した軌道があるかや Rmt の大きさにもよるので一概には言えないが、
 第四周期元素ならば経験上 5.0 くらいが良いか?
 重元素(原子番号50を超える)ならば経験上 6.0 くらいが良いか?
c) case.in2
  計算する系が金属の場合(フェルミ端の直ぐ上に状態がある場合)は、TETRAを TEMPS または TEMP にして、 その右の数値を  0.002 にする。 TEMPSやTEMPは計算する系が金属の場合、一般的に収束が速くなり易い。0.002 はディフォルト値。
d) case.klist
  x kgen で1桁か2桁の数字(例えば、8 など)を入れて、view klist で見たときにENDも含めて数行だけ書き込まれているよう x kgen での値を調整する。SCF収束後にDOSを見て、満足いかなければ x kgenで値を増やして再びrun SCF を行う。
※ 原子での電子密度分布を初期に用いるので、少ないk点数で結晶中での電子密度分布に近づけてから、k点を多くする方が結果的には速く収束する。私の周りの人だと、8 kや32 k, 64 kなどの小さいk点での計算の後に、256 k, 512 k, 1000 kにするなどしている。
※ 下記の中級者と上級者向けの変更は、通常のSCFで RKmax 5、1 k-point(または数点) として計算を行った後に行う。 
■ 中級者 (semi-core と valenceの範囲を変更)
a) x lstart (case.in0)
  x lstart で LSDA を選択。GGA-PBEでも良い。
b) x lstart
 -6.0 Ry から -3.0 Ry や -2.0 Ry(case.in1で 0のAPW+loを指定時) にする。
 ※ x lstart で入力した Ry に -3.0 Ry した値が case.in1 での emin に対応する。x lstart で -3.0 Ry などと入力して、 x lstart に 戻されても 0.97 から 0.9999 内に収まっていれば無視して次に進む(これで精度が良いのは、4fや5dがcoreに、5pや5dがsemi-coreにあると考えられる場合)。
c) case.in1( semi-core と valence での local orbital と LAPW, APW の指定など)
  RKmax を 6.0 や 5.0 にする。どうしてもメモリが足りないなどの場合は、5.00±0.04に近い値にする。
 価電子帯(-15 eV 程度か?)局在した軌道があるかや Rmt の大きさにもよるので一概には言えないが、
 第四周期元素ならば経験上 5.0 くらいが良いか?
 重元素(原子番号50を超える)ならば経験上 6.0 くらいが良いか?
 ※ semi-core と valence での local orbital と LAPW, APW の指定は x lstart で自動的に書き換えてくれる(ver. 12 で確認)。
 ※ 私の経験では、ゴーストバンドが表示された場合、 例えば、”QTL-B VALUE .EQ.  120.91788   in Band of energy    -2.12346   ATOM=    9   L=  2”(上記のゴーストバンドのエラーの対処方法を参照)と表示された場合に、case.in1または case.in1c でそれに対応する原子の軌道のEパラメーターに Band of energy に近い値を指定すると上手く計算が収束してくれたことがある。例えば、case.in1c で "2    0.30      0.005 CONT 1" ->  "2    -2.00      0.005 CONT 1" と変更したりする。それに加えて、単位胞の中で原子番号が似たものは全て同じ様に値を変更してみるとよい。
 ※ 一番右の 0はLAPWを使うことを指定している。APWでは1を記入。LAPWは、APWの計算式においてエネルギーの関数で対数微分する部分を線形に近似した方法のこと。APWの方が精度がよい。局在化している3d電子と4f電子だけをAPWに指定すれば良いことが多い。例えば、Znではd軌道まで、Mn, Alp軌道までloを入れて計算し、APW+loとしてみる。APW+lo法は、LAPW法と同じ結果に収束し、尚且つRkmak(平面波の最大波数の目安)の値が1小さくて済み、基底関数の数が最大50%少なくなり、計算時間は1桁短くなる可能性がある。
d) case.in2
  計算する系が金属の場合(フェルミ端の直ぐ上に状態がある場合)は、TETRAを TEMPS または TEMP にして、 その右の数値を  0.002 にする。 TEMPSやTEMPは計算する系が金属の場合、一般的に収束が速くなり易い。0.002 はディフォルト値。
e) case.klist
  x kgen で1桁か2桁の数字(例えば、8 など)を入れて、view klist で見たときにENDも含めて数行だけ書き込まれているよう x kgen での値を調整する。SCF収束後にDOSを見て、満足いかなければ x kgenで値を増やして再びrun SCF を行う。
f) case.inst(恐らく編集は必要ない)
 周期律表に記載されている[希ガスの原子]と他の軌道の記述を真似ればよい。
g) [run_SCF]にて、エネルギーの収束判定値を「0.0001 (Ry) * case.struct に記載されている原子の数(Posと書かれている行の数)」にする。1 meV / atom であれば十分である。※ 0.0001 Ry * 13.6058 = 0.0013 eV = 1.3 meVである。
■ 上級者 (semi-core にある軌道{ s や pなど}を core{case.inc} に指定する)(中級者のレベルでこれを行うと 1 cycle の時間は短いが、収束までの時間はあまり短くならない。私としては上記の中級者レベルまでの計算をオススメしたい)
下記は一例である。 (.lcore という名称の空のファイルを作成しておくと良いだろう)
a) x lstart (case.in0)
  x lstart で LSDA を選択。GGA-PBEでも良い。
b) x lstart
 -6.0 Ry から -3.0 Ry や -2.0 Ry(case.in1で 0のAPW+loを指定時) にする。
 ※ x lstart で入力した Ry に -3.0 Ry した値が emin に対応する。x lstart で -3.0 Ry などと入力して、 x lstart に 戻されても 0.97 から 0.9999 内に収まっていれば無視して次に進む(これで精度が良いのは、4fや5dがcoreに、5pや5dがsemi-coreにある考えられる場合)。
c) case.in0
  LSDA となるように TOT の次の数字を 5 にする。
  min IFFT-parameters: LSDA なので 0 0 0 とする。
  enhancement factor: 1.00 とする。
  ※ x lstart で LSDA を選べば同然ながら、TOT の次の数字は 5 になっている。
d) case.in1 (恐らく、x lstart でのsemi-coreとvalenceの範囲の変更には自動的に対応して書き換えてくれるだろう)
 case.in1から軌道を除いた場合は、case.inc に軌道を記述する。多くの場合、価電子帯にあるs軌道をcase.in1から除いた場合には、特にcase.incに記述する必要は無い(計算する平面波で十分な精度が得られることがある。ただし、その場合、RKmaxは6.0以上にしておくのが良いだろう)。-1.0 Ry 近傍にあるような軌道をcase.in1から除いた場合には、case.incに軌道を記述する。
  RKmax: 5.0 から 6.0 の値 ← ケースにより様々なので 中級者から上級者へとなってきた方は 6.0 以上としておくのが無難。
  emin: -6.0 
  emax: 2.0 
  ※ この他に、 coreに属するものは除いて、上記のcase.in1のQ&Aを参考に記述する。
  ※ x lstart で入力した Ry に -3.0 Ry した値が emin に対応する。x lstart で -3.0 Ry などと入力して、 x lstart に 戻されても 0.97 から 0.9999内に収まっていれば無視して次に進む(これで精度が良いのは、4fや5dがcoreに、5pや5dがsemi-coreと考えられる場合である。経験上は semi-core(x lstart で{入力したRy} から {入力した Ry -3.0 Ry})の範囲にdがある場合などでも精度よく計算可能である)。case.in1とcase.inc が連動して最適な値に変わっているかは確認が必要。
  ※ これまでの経験などで EF が分かる場合は、それに近い値を入力する。そして、Elの値も入力したEFを基準に、(XPSなどのデータブックを見るなどして)存在しそうなエネルギー(Ry)の値を書き込む。
 ※ SOを入れた計算をする場合には、EFから遠い方のDOSの中心部分を El とする(調査中)。理由:SOを入れてSCF計算した後、SO無しでDOSを描くと高束縛エネルギー側のDOSだけが主に描かれていたため。恐らく、4f などでは、SOを入れると低束縛エネルギー側にスプリットした準位が生じるためであると考えている。
 ※ 局在していないs 軌道(l=0で、EF近傍にあるようなもの)がLAPWやAPWとしてcase.in1に記述してある場合、それを除いたときは、RKmax を 6.0 以上にする(調査中)。
 ※ 私の経験では、ゴーストバンドが表示された場合、 例えば、”QTL-B VALUE .EQ.  120.91788   in Band of energy    -2.12346   ATOM=    9   L=  2”(上記のゴーストバンドのエラーの対処方法を参照)と表示された場合に、case.in1または case.in1c でそれに対応する原子の軌道のEパラメーターに Band of energy に近い値を指定すると上手く計算が収束してくれたことがある。例えば、case.in1c で "2    0.30      0.005 CONT 1" ->  "2    -2.00      0.005 CONT 1" と変更したりする。それに加えて、単位胞の中で原子番号が似たものは全て同じ様に値を変更してみるとよい。
 ※ 一番右の 0はLAPWを使うことを指定している。APWでは1を記入。LAPWは、APWの計算式においてエネルギーの関数で対数微分する部分を線形に近似した方法のこと。APWの方が精度がよい。局在化している3d電子と4f電子だけをAPWに指定すれば良いことが多い。例えば、Znではd軌道まで、Mn, Alp軌道までloを入れて計算し、APW+loとしている模様。APW+lo法は、LAPW法と同じ結果に収束し、尚且つRkmak(平面波の最大波数の目安)の値が1小さくて済み、基底関数の数が最大50%少なくなり、計算時間は1桁短くなる可能性がある。
e) case.inc (x lstart でsemi-coreとvalenceをシフトさせた分、coreに追記する軌道が増える。計算していてゴーストバンドのエラーが出ている場合、ここでの記述を忘れていないか確認するようにしておくとよい。加えて、case.in1で軌道をディフォルトから省いた場合も記述が必要になる)
  一般的な周期律表に、計算で扱う原子の[希ガスの原子]と他の軌道が書かれていると思うが、[希ガスの原子]のところまで記述を書き込む(追加する)。
  NUMBER OF ORBITALS: 追加する内殻の軌道の数を入力
  SHIFT: 0.00
  N, KAPPA, OCCUP を入力する。
  OCCUP は KAPPAの絶対値の2倍
  ※ case.inst は UPとDNを分けているので、KAPPAの絶対値
 ※ error case following,
      'CORE' - NSTOP= 362 positive eigenvalue for  5P  Atom:  20 Yb               
      'CORE' - Try to apply a potential shift in case.inc   
  change: [17 0.00  0  NUMBER OF ORBITALS (EXCLUDING SPIN), SHIFT, IPRINT] -> [17 1.00  0  NUMBER OF ORBITALS (EXCLUDING SPIN), SHIFT, IPRINT]
 ※ shift: shift of potential for “positive” eigenvalues (e.g. for 4f states as core states in lanthanides)
f) case.in2
 Emin: x lstart で入力したい値を記述する。中級者から上級者に移ってきたばかりの人はディフォルトのままでよい。
 ne: x lstart で入力した値よりも低束縛エネルギー側(=valence)の電子数(=「case.structまたは StrcutGenTMに記載されている原子の数」 * 「x lstartで入力した値よりも低束縛エネルギー」)となっている。
  計算する系が金属の場合(フェルミ端の直ぐ上に状態がある場合)は、TETRAを TEMPS または TEMP にして、 その右の数値を  0.002 にする。 TEMPSやTEMPは計算する系が金属の場合、一般的に収束が速くなり易い。0.002 はディフォルト値。
  Gmax: 12.00 ←12.00近辺を入力して Gmax > Gminとなるようにする。Userguideに推奨されている最小値は14.00。
 Gmaxが12.00 よりも小さな値でも大丈夫か調査中。
g) case.inm
  mixing FACTOR for BROYD/PRATT scheme: 0.10 ←ミキシングファクター(Userguide 7.10 MIXERに記載がある。新しく計算された電荷密度分布の割合を記述している。ディフォルトは0.2。ちなみにDVXa法SCATであれば0.3である。最初は大きな値にして数サイクル計算後に小さい値にしたり、初期値の電子密度分布を適切な値に設定していれば小さい値で始めたりする)
  PW and CLM-scaling factor: 0.40 1.00 ←電荷が振動する場合は1.0より小さくする(0.4とか0.5とか)。CLM-scaling factor はUserguideに1.00とするとある。
h) k点は行数が、並列化で用いるcore数-1 (spin有りの場合は、(並列core数-1)/2とするとよい?) になるようにする。
 (「並列計算するcoreの中で、1つだけは統括のような役目をして、計算にはあまり寄与されない?」だったかな?)
i) case.inst (恐らく、編集の必要は無い)
 周期律表に記載されている[希ガスの原子]と他の軌道の記述を真似ればよい。
※ 入力ファイルは、まず _st と付いたものが作られ、Prepare input file で _st の無いファイルにコピーされる。実際に計算に使われるのは _st の無いものである。虚数を平面波展開での係数に含む場合には、入力ファイルの最後に _st を取り去った後に、c が付く。
j) [run_SCF]にて、エネルギーの収束判定値を「0.0001 (Ry) * case.struct に記載されている原子の数(Posと書かれている行の数)」にする。1 meV / atom であれば十分である。※ 0.0001 Ry * 13.6058 = 0.0013 eV = 1.3 meV である。

Q. k点数はどうやって大きくすれば良いでしょうか?
A. SCF や DOS 計算の後に、k点を大きくしたい場合は下記のようにします。
通常: [initialize calc.] -> [x kgen] -> k点を入力 して、Execute ! を押す -> [run SCF]
SOの場合: [initso_lapw] -> [x kgen -so] -> k点を入力、Execute ! を押す -> [run SCF] 

Q. 本来、あり得ないと思われるところに非常にシャープなスパイク構造が見えるのですが・・・・・・。
A. k点数を多くしてみてください。k点数が少ないと、スパイク的な構造が見える場合があります。上記の[k点数はどうやって大きくすれば良いのでしょうか?] を参照してk点数を大きくしてみてください。

Q. ".lcore"という名称の空(0MB)のファイルはどういう意味があるのですか?
A.  ".lcore"はマフィンティン球(muffin-tin spheres)からリークした(染み出た)内殻電子の電荷密度に対するテール(裾)の補正を組み込ませます。例えば、彼方が、準内殻ではなく、内殻の状態としてGa-3d 電子を選択すれば、このオプションは必須になります。
英文:".lcore" includes the tail correction for the charge density of core electrons leaking out of the muffin-tin spheres. For example,  you select the Ga-3d electrons as core states (not semicore), so this option is essential. (Dr. Eeuwe S. Zijlstra)
"touch .lcore"コマンド(Linuxでのターミナル上で touch .lcore と入力すると .lcore という名称の空のファイルが作られる)は内殻の状態とした電子に対するテール(裾)の補正するプログラムを起動させます。これを用いないと、多くても 0.01の電子がディフォルトとしてマフィンティン球(muffin-tin spheres)からリークするでしょう。この補正を用いると、0.1よりも多くの電子が結果に影響することなしにリークできます。それ故、"touch .lcore"オプションを用いると、内殻に多くの電子準位を含ませることができ、大規模計算を若干経済的(計算コストの少ない割安)な実行形式にさせることができます。
英文:The "touch .lcore" command invokes tail corrections for the core states. Without it, at most 0.01 electrons may by default leak out of the muffin-tin spheres. With the corrections, as much as 0.1 electrons may leak out without impacting the results. Therefore, with the "touch .lcore" option, we can include more electronic levels in the core and make a big calculation a little bit cheaper to perform. (Dr. Eeuwe S. Zijlstra)
(書き換える部分は、case.inc に内殻とする電子を書き入れ、case.in1から内殻とした電子のlの記述を省きます)

Q. どのプログラムでどのファイルが出力されるか分かりますか?
A. Userguide に記載されています。
 program needs generates という欄があるテーブルがUserguide にあります。
 例えば、KGEN ( x kgen )を実行すると case.outputkgen(IBZ内で全てのk点に対する情報)や case.klist, case.kgenが出力されることが分かります。necessaryで分かるように、_stが付いたファイルを入力としているプログラムはSYMMETRYを除いてありません。直接ディレクトリにあるファイルを編集する場合には_stの付いていないものにすることが必要であることも理解しておくと良いでしょう。

Q. Spin orbit couplid の入力ファイルを initso_lapw で作成すると対称性が低下します。
A. case.structの一番下に対称操作の行列が記載されています。このデータは、SRC_symmetso/symho.f では transに格納されます。例えば、磁化の方向をz軸にした場合には、trans(1,1)*trans(2,2) - trans(1,2)*trans(2,1) を計算して、それが0.01の範囲で0に近ければ、対称操作から除き、1であればAと分類してその対称操作を残し、-1であればBと分類して対称操作(z軸に時間反転のものとなる)を残します。残された対称操作から空間群を探し出します。 
 SOの場合は、case.struct での空間群の記載があてになりません(ver.12)。case.struct の最下部に書かれている対称操作を見て、International tables A から対応する空間群を探すべきです。例えば、I23 (No.197)の場合には、International tables Aで [ 197 (I23) ] を押して、Maximal non-isomorphic subgroup にある空間群を調べます。Positions の x, y, z などと書かれているところと、case.structの最下部が対応しているかを調べます。そうすると、23 (I222) が対応していることが分かります。つまり、SOの計算のために対称性を低下させて、197 (I23) -> 23 (I222) になっていることが分かります。

Q. RKmaxを大きくしすぎると結果が正しくないような気がするのです。
A. ブラハ教授に聞きましたが「線形性が壊れる」からだそうです。 
他のコードで似たような事例としては下記があります。
「fcc, hcp, bccなどの結晶構造をもつ高密度のバルク系において、多くの基底関数を割り当てた場合には、基底関数が過完備になる可能性があります。そのような場合には、数値不安定性のために誤った固有値が生じる可能性があります。この過完備性を避けるために、最適化された基底関数を使用することが推奨されます」(OpenMX HP)

その他のQ&Aなどについては下記を参照してみてください。
WIEN2k-FAQ: How to set a good RKmax value?
 
Frequently asked Questions
 
WIEN2k software package

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

参考文献
[1] ナノ・物質・材料・マルチスケール機能シミュレーション PHASE ver.9.00 ユーザーマニュアル
[2] 山内淳 著、第一原理計算を理解するために : http://ci.nii.ac.jp/els/110004800210.pdf?id=ART0007532250&type=pdf&lang=jp&host=cinii&order_no=&ppv_type=0&lang_sw=&no=1302015226&cp=
[3] 赤木和人 著、企画の趣旨 : http://www.jstage.jst.go.jp/article/jsssj/28/3/128/_pdf/-char/ja/ (以下は表面科学での特集としてシリーズとなっている)http://www.jstage.jst.go.jp/browse/jsssj/28/3/_contents/-char/ja/ 
[4] 小林一昭 著、 バンド計算と実験との距離 : http://www.jstage.jst.go.jp/article/jsssj/28/3/129/_pdf/-char/ja/
[5] 山内淳 著、第一原理計算の諸条件 : http://www.jstage.jst.go.jp/article/jsssj/28/3/135/_pdf/-char/ja/
[6] 吉本芳英 著、第一原理計算による構造最適化とその周辺 : http://www.jstage.jst.go.jp/article/jsssj/28/3/144/_pdf/-char/ja/
[7] 中井浩巳 著、表面-分子相互作用系の量子化学計算に関する最近の動向 : http://www.jstage.jst.go.jp/article/jsssj/28/3/150/_pdf/-char/ja/
[8] 赤木和人 著、第一原理計算パッケージを使ってみる : http://www.jstage.jst.go.jp/article/jsssj/28/3/160/_pdf/-char/ja/
[9] Tsuchiura et al., A "Beginners guide" to WIEN2k :
http://www.apph.tohoku.ac.jp/sakuma-lab/docs/wien2k.pdf
[10] M. Gmitra et al., Phys. Rev. B 80 (2009) 235431.
  http://prb.aps.org/pdf/PRB/v80/i23/e235431
[11] http://www.wien2k.at/events/ws2008/talks/Laskowski-SO-NCM.pdf
[12] http://www.mri.psu.edu/conferences/WIEN/ProgramFiles/printversie_relativity.pdf
[13] http://www.wien2k.at/events/ws2006/magnetism-so.pdf
[14] http://interface.t.u-tokyo.ac.jp/home/teru/wien2k.html
[15] G. Bertoni et al., Micron 37 (2006) 486.
  http://www.emat.ua.ac.be/pdf/1387.pdf
[16] L. Pardini et al., arXiv:1104.1558v1, (2011).
  http://arxiv.org/PS_cache/arxiv/pdf/1104/1104.1558v1.pdf 
[17] 笠井秀明ら著、計算機マテリアルデザイン入門、大阪大学出版会
[xx] http://www.mri.psu.edu/conferences/WIEN/ProgramFiles/WIEN2k-getting_started-07.pdf
[X1] 和光システム研究所、固体の中の電子、WSL
------------------------------------------------------------------------------

QRコード
携帯用QRコード
アクセス数
ページビュー数