Pmt method

From Ecal

Jump to: navigation, search

Contents

バンド計算の方法

ポイントになるのは、

  • 波動関数
  • 電子密度
  • 一体ポテンシャル(local)

がどうやって表現されているか?です。


PMT波動関数

PMT

PMTはaugmenationの方法であり、APW(PWをaugmentしたもの)とMTO(smooth Hankel関数のブロッホ和)をaugmentしたものを用いて対角化する方法である。 全領域をMuffin-tin(MT)内とinterstitial region(MT外)に分ける。以下、ノーテーションに注意する。

まず対角化のための基底関数を作るためにenvelope関数のセットを用意する。これを\{{F}_{0i}({\mathbf r})\}と書く。 iの数が対角化の次元(ハミルトニアンの次元に対応する。ここでは波数ベクトルは略する)。 上述のようにenvelope関数として、PWとsmooth Hankelがある。 MTOに関しては、「vMTO」と「cMTO(セミコア軌道)」がある。cMTOでは、

  1. セミコアのエネルギーレベルでradialシュレディンガー方程式を解き(この解をphi^cと呼ぶ)、MT端での波動関数のradial微分、二回微分と値を求める。
  2. それらの値にマッチするようにsmooth Hankelのパラメーターを選ぶ。
  3. EHは自動的に求められる。(要確認)

というようにして作っている。cMTOの中心原子では、phi^cでaugmentすることになる。

smooth HankelとはHankel関数を修正したものである。単純にはヘルムホルツ型方程式で原点にデルタ関数的なソースをおいた解(exp(-kappa r)/r)
のかわりに,Gaussian的にひろがったソースをおいた解を考えるものである。

index iを持つ基底関数を F_{0i}({\mathbf r}) が全域でnonzero、 F_{1i}({\mathbf r}) F2i(r)がMT内でnonzeroに定義された関数。 F_{2i}({\mathbf r}) がMT内でnonzeroでsmoothな成分で、F_{0i}({r})とMT内では同じである。

{F}_i({\mathbf r}) = F_{0i}({\mathbf r}) \oplus \{ F_{1i,a}({\mathbf r}) \} \ominus \{ F_{2i,a}({\mathbf r}) \}
F_{1i,a}(r)= \sum_{akL} C_{akl}^i \tilde P_{akL}(r)
F_{2i,a}(r)= \sum_{akL} C_{akl}^i P_{akL}(r)

と構成する。 ここでaは原子index,kは後出するLaguerre多項式(PakL(r))のindex、Lは角運動量を意味する。 (紛らわしいがここでのkは波数ではない。) \tilde P_{akL}(r)はあるenergyのall electronの波動関数φ\dot \phiを用いてPakL(r)rMTでスムーズにつながるようにした関数になる。

F0,F1,F2はeq(PAW)のそれぞれの3つの項と対応している。 F2i(r)F0i(r)をMT内へprojectionをかけた結果であり、 PAWの説明と同じくMT内で完全系からprojectionを作っていればF0 = F2である。 PMTではF0i(r)のinterstitialの基底、index i、はplane waveと多くの場合は原子を中心としたsmoothed Hankel関数eq(smoothed_Hankel)を用いる。 ((F)LAPWはplane waveだけを用い,(fp)LMTOは(smoothed) Hankelだけを用いる。) 実はbulkにおいて原子の性質が残っているtailの部分をplane wave展開するのはかなり非効率的であり、 PMTでは両者を用いることによりplane waveの数を大幅に減らすことが可能となる。

  • {F}_i({\mathbf r})の表式を額面どうり受けとるのは注意する必要がある。計算において、

3つの成分を実際に加算することは全くない。むしろF_iが3つの成分をもつと考えるべき。

smoothed HankelのF

F0i(r)を張る smoothed Hankelはatomの波動関数のtailをfitして作成する。smoothed Hankelのheadはtailから自動的に決まる。 実際はheadは後でAE波動関数のF1iと差し替えるだけであるからある程度smoothな扱いやすい関数なら何でもいいのである。

訂正: F0i(r)はatomの波動関数のtailからつくるのではなくEH、EH2という2つのenergyをもつsmoothed Hankelで生成する。 PlanewaveとEH,EH2の係数がsecular equationを解いて求まり、F2の自分の原子の部分(head)とまわりの原子からの寄与(tail)をMT内でなめらかに繋ぐという条件からF2の係数C_kLが(解析的に)求まる。その係数C_klを用いてF1が作られ、F0+F1-F2としてall electronの波動関数が求められる。

smoothed Hankelは原子を中心に作成すると F2i(r)は中心原子のheadと他原子からのtailからなる。(この部分はfpLMTOと同じである。) どちらもMT内をF1iでaugumentし全系のAE基底であるFiを作成する。

通常はatomのvalence orbital {nL}全てのnLのsmoothed Hankelを使用して基底を張るが、overcompletenessの問題が生ずるsmoothed HankelはF0の基底から除く必要がある。 (このovercompletenessはPMT固有の問題ではなく同じ形式を持つFLAPW,PAWにも発生する。)

他サイト展開

F2の形であるが、MT内で角運動量展開する。

planewave

 eq(pw-sbessel):

\exp(i\vec r\cdot \vec k)=
\exp(irk\cos\theta)=\sum_l (i)^l(2l+1)j_l(kr)P_l(\cos\theta)
 = \sum_{lm} (i)^l 4\pi (2l+1)j_l(kr) Y_{lm}^*(\hat k) Y_{lm}(\hat r)

\hat x=\vec x/|x|、 と球ベッセル関数jlで展開できることはよく知られているが、ここではLaguerre polynomial Lで展開する。J.Mathの文献。 eq(exp-Laguerre)に従い、


\exp(iqr)= \sum_{akL} C^{q}_{akL} P_{akL}(r) 

と展開しF2の境界条件を求めるのに用いる。(紛らわしいが)summationのLは角運動量L=(l,m)で、kはLaguerre polynomialの足である。

ソース

bstrux.F, pauggp, kmax,napw, b=C_kL^i


smoothed Hankelのtail関数

他のMTからの寄与をLaguerre polynomial L (eq(Laguerre-poly))を用い


P_{kl}(r)=\frac{ (-2)^k k! (2l+1)!! }{(2k+2l+1)!!} L_k^{l+1/2}(\frac{r^2}{r_s^2}) 
\frac{r^l}{r_s^l} Y_L(\hat r)

F_{2i}(r)= \sum_{akL}C^{i}_{akL}{P}_{akL}(r)

とLaguerre polynomialを用いて展開できる。C^{i}_{akL}も解析的に計算できる。(see eq(Hankel_expansion))

F1

planewaveもsmoothed Hankelも F_{1i}(r)= \sum_{aL}\tilde{F}_{aL}(r)の方は logderiもしくはenuを指定してMT内で解いたφnLとそれに線形独立な\dot\phi_{nL}と の線型結合 c_{1nL}\phi_{nL}+c_{2nL}\dot\phi_{nL}F0(もしくは同等であるF2)と 角運動量展開の範囲内で一回微分まで滑らかに接続する。ここでのc1,c2に不定性はない。 enuは{nL}のpartial DOSの重心と合わせる。

(F0がplanewaveのみでphi,phidotを毎回計算しないでphi,phidotをatomの値に固定するとPAWと等価なformalismになる。)

H,S

HamiltonianHij、overlap matrixSij はinterstitialのenvelope関数の次元を持つ。 PMTはPAWと同様にtailの足持つという意味でtailの部分を決めるとheadが決まる構造になっている。 ((fp)LMTOはHamiltonianが原子の足を持っているという意味でMT内のheadを決めるとからtailが決まる。 単なる見方の問題でもちろん両者は等価である。)

overcompleteにならないよう注意がいるが interstitial基底の張りかたがかなり効率的になっており同じ精度を出すとして bulkでは1/10程度と行列の次元が大幅に減るためflapwより高速に対角化できる。 overcompletenessを議論する際に原子の足とplanewaveの足を持っていたが、原子の足の方はsmoothed Hankelの足である。


実際の計算では

  • smoothed Hankelのtail, planewaveの展開のsummationは

\sum_{k=0}^{\rm KMXA}\sum_{l=0}^{\rm LMXA}と上限を定めている。

  • planewaveの数はgmax or napwで与えられる。
  • schreodinder equationは各l(各運動量)ごとに解く。

論文では各lでDOSの重心位置で解くある。

Comment:

DOSの重心で解くということなっているがMarkさんのoriginal コードにはバグがありDOSの中心にはなっていなかった。現在は修正されており、正しくDOSの中心にする該当OPTIONはOPTIONS_PFLOAT=1がある。しかし、test suitsには過去の結果との数値チェックのため=0になっているものもある。 なお、DOSの中心に最適化されていなくてもφ\dot\phiの線形結合で展開されているためか現実的に問題はほぼないことは確認している。

PMT charge

重なり積分の類は


F_i(r)F_j(r)= F_{0i}^*(r)F_{0j}(r) + F_{1i}^*(r)F_{1j}(r) - F_{2i}^*(r)F_{2j}(r) 

と書く。ただしF1,F2はMT内のみで定義されており MT内が完全系でprojectされていればF0i = F2iであるから MT内でF_{0i}^*(r)F_{0j}(r)=F_{2i}^*(r)F_{2j}(r)が成り立つ。

F1,F2はMT内のみで定義されているので 重なり積分は


O_{ij}=
\langle F_i|F_j\rangle =
\int\! d^3 r F_i(r)^*F_j(r)= \int\! d^3 r \; F_{0ia}(r)^*F_{0ja}(r) + 
\sum_a \int_{r<R_{MT,a}}\! d^3 r \; \left( F_{1ia}(r)^*F_{1ja}(r) - F_{2ia}(r)^*F_{2ja}(r) \right)

となる。ここでFxia(r)は原子aでの値である。

Hamiltonian matrixも問題なく定義できる。 secular equationは

Hαp = εpOijαp

であり、ここでαpは固有値εpの固有ベクトルである。 対応する波動関数はψp = αpFである。

total charge densityは重なり積分の書き方にならい、 
n(r)=\sum^{\rm occup.} \psi_p^*(r) \psi_p(r) = n_0(r) + n_1(r) - n_2(r)
と書く。

Force

全エネルギーの表式は、\sum_{i,j} rho_ij F^*_i(r)F_j(r)の汎関数であるが、原子位置が固定されているなら rho_ijによる変分はゼロであるので,原子をうごかしたときのF^*_i(r)F_j(r)の変化分による項 (Pulay項)と静電相互作用に入っている原子位置からくる項のみを考えればよい。

このとき、PakLは、固定されているし、\tilde P_{akL}の変化分を 考えないならC_{akL}^iの変化分と、F^_i(r)の第0成分の変化分を考えればよい。 これがPulay補正である。

core charge

(sectionが一段階上かもしれない。) frozen core近似=coreのchargeは変化しない、通常はatomのもの、を用いる。主に2つの扱いがある。

  • LFOCA=0だとMT内に閉じ込めるようにcoreの波動関数を求めて、そこからcore chargeを出す。
n(r) = nMT,T
  • LFOCA=1だとatomのものを用いる。必ずしもMT内に局在するとは限らない。全域とMT内に分ける。
n(r) = nSH(r) + nTrue,MT(r) − nSH,MT

LFOCA=1が正確な取扱いである。 nTrue(r)を波動関数から求めた真のcore chargeだとして、nTrue,MT(r) はそのMT内の成分である。 ここでnSH(r) はsmoothed Hankelでcore charge tailをfitした関数であり、 Kinoが見たところn_{SH}(R_{MT})\ne n_{True,MT}(R_{MT})で、tailの値がほぼ一致するのは現在の決めかたではMTの外側である。ここから誤差が生じている。

小谷reply:n_{SH}(R_{MT})\ne n_{True,MT}(R_{MT})であるために mtopw9j.pdfの多重極変換してから静電相互作用を計算する方法だとすこしは誤差を生じます。ですがそれほど深刻ではないとおもう。そのあたりの評価もある程度書きました。

interstitial regionを求めるための外部から見て等価な電荷への変換

全電荷は nT(r) = ncore(r) + nvalence(r) として与えられる。

interstitial regionのCoulomb potentialとG空間のchargeを求めるためには等価な電荷分布と置き換え、 nTのフーリェ変換の収束を速くする。(実際はフーリェ変換しているのはinterstitial regionはn0の部分である。) これはcore chargeはpeakyであるし、valenceにおいてもノードを持った波動関数もcore程でないがまたpeakyだからである。以下、足したり引いたりして変換を行うがn0-n2まで全部足せばn=\tilde nである。

(planewaveだけを用いると波動関数のcutoffがGなら電荷のcutoffはO(G^2)であるがこの場合のFFTのcutoffは自明ではない。smoothed Hankel tailがあるのでG^2より大きく必要、というよりはsmoothed Hankelの展開がうまく行えるという方の条件で決まっているはず。このあたりはfpLMTOでの経験を使える。)

ここでのn \rightarrow \tilde nの変換は、「smoothにするため」ではなくて「外側からみて等価な電荷」にするための変化です。静電ポテンシャルを計算するとき、「格子上の点であらわされた電荷」はFFTで解きますが,「Gaussianであらわされた電荷」は格子に落とさずに扱います。

「外側からみて等価な電荷」にするのは手段で目的ではありません。 その目的=意味の一つは\tilde n_0をfourier変換してinterstitialのCoulomb potentialを求めていますが、そのG cutoffを小さくすることです。--Kino 12:38, 15 January 2010 (JST) \langle \tilde n_0(site1) | v | \tilde n_0(site2) \rangleという項を求めるときに...

core

coreの場合は電荷は球対称である。MT内の電荷と同じtotal電荷を持つgaussianの和で置き換える演算子\mathcal Gを定義する。 Gaussの定理によりn(r)でも{\mathcal G}(n(r))でもMTの外では両者は同じポテンシャルを与える。

MT内の電荷と同じtotal電荷を持っていれば各項の外部からの見えかたは同じわけだが、 どういうわけか以下の様に分割している。


\tilde n(r)=\left( n_{SH}(r) + {\mathcal G}(-Z\delta(r) + n_{True,MT}(r)- n_{SH,MT} \right)


+ \left(-Z\delta(r) + n_{True,MT}(r))  \right)


+ \left( n_{SH,MT} + {\mathcal G}(-Z\delta(r) + n_{True,MT}(r)- n_{SH,MT} \right)

ここで\tilde n\mathcal Gにより変換されたn0,n1,n2の電荷分布を持つnを示す。 Zはこの原子の電荷で、各項がn0,n1,n2相当である。 十分遠方からみるとtotal charge-core charge=valence charge分の電荷が見える。次のvalence chargeでその分が足される。 LFCOA=0ではn_SH=0として上の式を用いている。 上の議論は電荷を多重極子モーメント(角運動量展開)と置き換えても成り立つ。

valence

valence部分もcoreと同様に各項の多重極子が変化しないように\mathcal Gのa siteへの変換op\mathcal G_aを用い変換したnを定義する。ここでの\mathcal Gは多重極子モーメント(角運動量毎の電荷量)への変換である。 全体としてn=\tilde nだが各項n0,n1,n2が違う。


\tilde n(r)= \left( n_0(r) + \sum_a \mathcal G_a( n_1(r-R_a) - n_2(r-R_a) ) \right)


+ n1(rRa)
a


 - \sum_a \left( n_2(r-R_a) + \mathcal G_a( n_1(r-R_a) - n_2(r-R_a) \right)

これらがvalenceのn0,n1,n2となる。

全電荷

\mathcal G_aにより変換された全電荷は 
\tilde n^T = \tilde n^{valence} + \tilde n^{core}
となる。 
\tilde n^Tもn0,n1,n2に分けられている。 \tilde n_0はもともとの分布よりsmoothになっており、gaussianとsmoothed Hankelなので解析的にフーリェ変換できる。n1,n2の方はMT内で動系方向、L毎に数値的に方程式を解くことによりpotentialを求める。(境界条件はL毎のpotentialの値。)Poisson eq.は線形なので、重ね合わせて最終的なpotentialが得られる。

上では3つに分けたが、実際は

  • smooth density\tilde n_0からpotentialのsmoothed partを解き、
  • \tilde n_1+ \tilde n_2からrapid partを解く。rapid partはMT上のL-channel毎のpotential vval(ilm)を境界条件として各サイトで動径方向へ方程式を解いている。

と2つに分けて解いている。 FFTを行いG空間で値がでるものは\tilde n_0(G)、とinterstitial regionのV0(G)である。




ソース解説

lm7K解説の最後にもある程度ヒントはある。ただ古くなっている。 まずpi,sigma,tau積分の定義は lm7K/MarkOriginalDoc/nfp-doc.ps.gz (M. Methfessel, M. van Schilfgaarde, and R. A. Casali)の式(21),(23),(28),(29)にある。

For charge and kinetic energies,

  • (21) \sigma_{k,k',l}=\int_S \{ \tilde P_{kL} \tilde P_{k'L} - P_{kL}P_{k'L} \}
dr
  • (23) \tau_{k,k',l}=\int_S \{ \tilde P_{kL} [-\Delta]\tilde P_{k'L} - P_{kL}[-\Delta]P_{k'L} \}
dr

For potential energy integral,

  • (24) n(r) = n0(r) + {ρ1(r) − ρ2(r)}
  • (25) \tilde n_0(r) = n_0 + \sum_M q_M G_M(r), where G_M is a localized gaussian of moment unity with angular momemtum M.
  • (27)  Q_{kk'LL'M}=\int_S \{ \tilde P_{kL} \tilde P_{k'L} - P_{kL}  P_{k'L} \} r^{m}Y_M(\hat r) dr
  • (28) \pi^{\rm mesh}_{kk'LL'}=\sum_M Q_{kk'LL'M} \int \tilde V_0 G_M dr
  • (29) \pi^{\rm local}_{kk'LL'} = \int_S \{
\tilde P_{kL} V_1 \tilde P_{k'L'} -P_{kL} \tilde V_2 P_{k'L'} \} dr - \sum_M Q_{kk'LL'M} \int_S \tilde V_2 G_M dr

注意:コードの中核はbndfp.Fにある。このルーチンで、与えられた電子密度から(最初は原子ファイルを読み込んで重ね合わせる)、 新たな電子密度を作る(mixingもする。mix.*ファイルにmixing情報は書き出す)。 bndfp.Fはかなり冗長になっている。lmfgwの新しいモードを試験的に組み込むのが途上になっているのもある(これはPMTでQSGWをするのを目指している)。 subs/switch.Fでnewsigmaswitch()はFになっているが、いずれTにする予定。sigmamodeは自己エネルギーを読み込まないときはF。

Note: smoothed Hankel functionの展開

  • まず、head partとtail partの定義をする。ある原子を中心としたMT(muffin-tin)領域Sが定義される。この領域S内にあるsmoothed Hankelの成分はこの原子を中心としたもの(head)と隣の原子を中心としたもの(tail)と2通りある。

(21)式は理論上はplanewaveもsmoothed Hankel(headとtail) もLaguerre polynomialに展開して計算する式になっているが、実際はhead partのsmoothed HankelをLaguerre polynomialに展開すると(形式は綺麗だが)収束がのろくて効率が悪いことが分かった。そこで tail partの他サイト展開の部分はLaguerre polynonialに一度直すが、head partのsmoothed HankelはLaguerre polynomialに展開せずに計算を行っている。


具体的にはMTを中心とするsmoothed Hankel自分自身と他サイト展開からきているLaguerre polynomialからF2を書くことになる。つまり、あるMT内の領域Sで定義するF2は


F_2=\sum_m H_{mL}+ \sum_{k=0}^{\rm KMXA}C_{kL} P_{kL}

と書く。 ここでHmLはEH,EH2で指定したvMTOのsmoothed Hankel、もしくはcMTOのsmoothed Hankel関数のhead部分のみ。 係数C_kLは周りのsmoothed Hankel tailをLaguerre polynomialで展開した成分のみ、となる。F2は一意に与えられる。(ここでindex kは波数ではなくLaguerre functionのindexであることを思い出して欲しい。)


この表式を用いたF2F2の計算は、

  • HmLHm'L, HH=sig(3)
  • (HmLPi'L)Ci'L, H(m)P(i)=sig(2)
  • CiL(PiLPi'L)Ci'L, P(i)P(i')=sig(1)

とCの入り方の違いから内部では3通り考えている。 なお、radial方向のこれらのすべての積分はPkL(r)を数値arrayとして与えてgaugm()内で数値積分を行っている。 smoothed Hankel関数 H(r)を与えているのはfradhd(), Laguerre polynomial PkL(r)を与えているのがfradpk()である。 この定義拡張に応じて(20)式の実際の計算も少し異なるものになる。

(
F_2=\sum_{k=-2}^{\rm KMXA}C_{kL} P_{kL}
, where P_kL (k<0)はsmoothed Hankel, k=-1はEHのsmoothed Hankel, k=-2はEH2のsmoothed Hankel、C_kL=1 (k<0)と思ってもいい。こう書くと見掛け上は他の式を変えなくてもいい。)

ポテンシャル生成:mkpot(bndfp:L818)

bndfp.Fはこのmkpotからスタートする。これは与えられた電子密度(osmrhoが第0成分,orhoatが第1成分,第2成分である。ただしorhoatは第3成分としてcore成分rhoc=n^c_a(r)を含む)から、 ポテンシャルosmpotを生成するルーチンである。注意すべきは、sigma,tau,pi積分は、このルーチン内で球対称ポテンシャルをつくり計算されてしまうのでradial関数(波動関数を書くのに必要)は、このルーチンの外に出ない(bndfpのレベルではやりとりされない)。またosmpotに対応するポテンシャルは陽には出力されない(ただし次項参照)。

  • 対角化のあと、orhoatを再計算する必要がある;このときにもradial関数は必要である。これは、bndfp->mkroutで行われている。そこでは、rv_p_ov0 => ssite(ib)%rv_p_ov0によりradial関数生成のための球対称ポテンシャルを引用している(よくない引用方法。あたえるのは、mkpot-locpot:L377以下のdpscopだと思われる。こういうトランシーバー的な方法はいかん)。
  • sigma,tau,pi積分は、mkpot-locpot-augmat-gaugmで得られる。nfp-doc.ps.gzは正確でない:pkl関数(augmentationで差し引く部分を展開するためのもの)にsmHankelそのものが入っているのが略されている(これはsmHankelをGaussian polynomialだけでは展開しきれないため)。それで(PP-P~P~)タイプの積分が3種類になる。すなわち、(HH-H~H~)、(PP-P~P~)、(HP-H~P~)などとなる。gaugmが3回呼ばれているのはそのため(augmat.F:L390のC ... Hsm*Hsm、L401のPkl*Pkl,L412のHsm*Pklである)。まずノートを確立し、indexを拡張すれば、もっと単純化できるはず。

原則としてはこのsigma積分のところからradial関数をとりだすことができる

ハミルトニアンのindex

ndimh(ハミルトニアンの次元=基底の数)=nlmto+napwである。 nlmtoはMTO基底の数。napwはAPW基底の数。napwはqごとに変わる。

MTO部分

bndfp.Fにおいてgen_hamindex(bndfp.F:L1246)を通るようにしたのでMTO基底のindexが表示される。 (2011Apr26; https://github.com/tkotani/ecalj/commit/b49c04911fb50cbd9acfec02bb5e57c44dc95512) たとえば(ちょっと物理的には違和感のある例だが)、lmf実行時、

 gen_hamindex: --- Hamiltonian index ---
 ib l  k  offl(iorb)+1  offl(iorb)+2*l+1   trim(spec)
 1  0  1      1         1            Fe    # k=1はEH,iorb=1
 1  1  1      2         4            Fe      # iorb=2
 1  2  1      5         9            Fe      # iorb=3
 1  3  1     10        16            Fe    # iorb=4
 1  0  2     17        17            Fe      # k=2はEH2, iorb=5
 1  1  2     18        20            Fe      # iorb=6
 1  2  2     21        25            Fe      # iorb=7
 1  3  2     26        32            Fe      # iorb=8
 1  1  3     33        35            Fe      # k=3はcMTO, iorb=9
 2  0  1     36        36            Z2    # k=1はEH, iorb=10
 2  1  1     37        39            Z2      # iorb=11
 2  0  2     40        40            Z2      # k=2はEH2, iorb=12
 2  1  2     41        43            Z2      # iorb=13
 2  0  3     44        44            Z2    # k=3はcMTO, iorb=14

(上で#の後はここで追記したコメント) などと表示される。 これはFe,Z2の二種の原子がそれぞれセル内で1siteある場合でnlmto=44の例である。

hamiltonianのindexの並びと、ib,l,kの対応を表示します。たとえば、上の表ではhamiltonianのindexの37から39まではibas=2(spec Z2)のl=1,EHチャンネル(k=1がEH,k=2がEH2、k=3がcMTOです。)(kl=と書いてあったがk=に修正。)を意味します。iorbが軌道チャンネルindexで、上の場合、14チャンネルあると数えます。この場合、hamiltonianのindex 45以上は、APWに対応。 l=2,l=3などの場合の例えばoffl=2から4,offl=5から9の実球面調和関数の並びは、GWmanualなどにあります。

この出力を出す部分は2つある。

  1. fp/bndfp.F中のcall gen_hamindex()直後。(コメントアウトされていないか確認。宣言部分の初期化でonlyonce=.true.とされていることも確認。)
  2. subs/suham.F,gen_hamindex()中(らしい。)

非常に面倒なインデックスシステムとなっている。

またポインターを用いて通常では考えられないようなアクセスをしている部分もあるようでので変更の際は注意がいる。

APW部分

APW部分の基底に関しては、igv2xにGが格納されている。これはgvlst2で作られている。その他のG(電子密度展開のとき)も同じようなルーチンで作られている。 ただ注意すべきは、一ヶ所で作られて引用されるのでなく、on the flyで作られている場合が多いことである。例えば原子のSCFを毎iteration計算している。

napwには、APW基底の数が入っている。

Hamiltonian,overlap行列の生成

Hamiltonian行列とoverlap行列は、bndfp.Fのhamblで作られる。

 if (newsigmasw()) then
     call hambl(...)
 else 
     call hambls(...)
 endif

newsigmasw()=F(subs/switch.F)としてあり、hambls()の方を呼ぶがhambls()内部でhambl()を呼んでいる。 (newsigmasw()=Tの方はテスト段階。)

hamm(ndimh,ndimh,nspec),ovlm(ndimh,ndimh,nspec)となる。 nspecはLS couplingを入れる(HAM_SO=1)か入れないか(HAM_SO=0)で変わる。 HAM_SO=0なら(NSPIN=2であってもnspec=1、 HAM_SO=1ならスピンについての非対角項を加えnspec=4となる;このとき対角化の次元は当然2倍でスピンは混ざることになる。ソースの変更時にはこのあたりの場合分けも考慮しなければいけない。

Hamiltonian,overlap行列のonsiteパート(1,2成分)

これは、pi,sigma,tau積分に、Laguerre多項式の係数であったC^i_akLを乗じて足し合わせることで得られる。bndfp->hambl->augmblがこれを行う。

bndfp()

Cl   lcplxp:0 if ppi is real; 1 if ppi is complex
     lcplxp = 0
      lekkl = int(sctrl%pfloat)
      if (lso .ne. 0) lcplxp = 1
      if (isum(nbas,lldau,1) .ne. 0) lcplxp = 1

  call dfaugm(...,sv_p_osigx,...)  # allocate variables
  call mkpot( ...,sv_p_osigx,...) # call locpot call augmat
  call hambl()

mkpot()

  cal locpot( ,sv_p_osig, )

locpot()

  call augmat( ,sv_p_osig, )

augmat()

C ... Hsm*Hsm
   call gaugm(,  sv_p_osig(3)%v, )
C ... Pkl*Pkl
   call gaugm(,  sv_p_osig(1)%v, )
C ... Hsm*Pkl
   call gaugm(,  sv_p_osig(2)%v, )


augmbl( ,sv_p_osig, )

     call bstrux( ,b, )     #  b=CkL matrix
C   --- Add 1-center and 2-center terms ---
        if (lcplxp .eq. 0) then
            call augq12 (  , sv_p_osig ( 3 , ia ) %v , sv_p_oppi( 3 , ia )%v 
     .      , sv_p_osig ( 2 , ia ) %v , , b , ) # sighh=osig(3), sighp=osig(2), osig(1): Pkl-Pkl overlap matrix

        else
            call augq2z (  , sv_p_osig ( 3 , ia ) %v , sv_p_oppi( 3 , ia )%v 
     .      , sv_p_osig ( 2 , ia ) %v , s , b , )
        endif
C   --- Add B+ sig B to S and B+ ppi B to H ---
        call augqs3 (  , sv_p_osig( 1 , ia )%v    , b ,  ) 
 

sugqs3()

              g(k1,ilm) = g(k1,ilm) + sig(k1,k2,l,isp)*b(k2,ilm,i2)

          s(i1,i2) = s(i1,i2) + sum( dconjg(b(:,:,i1))*g(:,:) )

augq12

     .          s(i1,i2) = s(i1,i2) + sighh(ik1,ik2,l1,isp)

                cadd = sighp(ik1,k,l1,isp)*b(k,ilm1,j)
                s(i1,j) = s(i1,j) + cadd


この式(20)の係数C^(i)_akLはaugmbl()->bstrux()により返された変数b()に入っている。

  • (20) \int \tilde F_i^* \tilde F_j = \int F_i^* F_j dr + \sum_{kk'L} C_{kL}^{(i)*} \sigma_{kk'l}  C_{k'L}^{(j)}

ただし,実際にはF2はH(r)とPkL(r)からなっているため実際に評価している式は(20)式そのものではなく(20')式になる。

  • (20') \int \tilde F_i^* \tilde F_j = \int F_i^* F_j dr + \sum_{kk'L} C_{kL}^{(i)*} \sigma_{kk'l}^{(1)}  C_{k'L}^{(j)}
+ (\sum_{kk'L} \sigma_{kk'l}^{(i,2)}  C_{k'L}^{(j)} + h.c.)
+ \sigma_{kk'l}^{(i,j,3)}\delta_{ij}

ちゃんと書いていないが、sigma2はhead(H)-tail(PkL), sigma3はhead(H)-head(H)の寄与であるためこうなっている。

(ソースコードにおいてsigma2,sigma3の前は+でいいのだろうか。要チェック)

対角化ルーチン

bndfp.F

  call zhev_tk(..., t_zv, ...)
  call addrbl(..., t_zv, ...)
  • hamm,ovlmがzhev_tkへ送られ対角化される。対角化の結果、evl(固有値)、t_zv(固有ベクトル)が得られる。この固有ベクトル×基底関数の和で波動関数は作られる。

t_zvは以下のaddrblへ渡されevecとしてアクセスされる。

電子密度のvalence第ゼロ成分の生成

対角化で得られた固有ベクトルt_zvはaddrbl内ではevec変数でアクセスされaddrbl内で出力電子密度が計算される。データはaddrbl -> rsibl -> rsibl1,rsiblpと流れていく。これらでMTOとAPWのq+G平面波展開での値が計算されている。psi(ng,nspec,nblock)の並びはngがGの数(これは実空間メッシュのうちの半分位の数である。)、nspecは紛らわしいがスピン、nblockはバンドインデックスをブロック化したものとなっている。rsiblpに置いてはAPWがq+Gそのものなので簡単である。このルーチンは、「napw個のGの並びと、実空間メッシュのGの対応テーブル」を用いている。

変数nlmto,napw,ndimhがそれぞれsmoothed Hankel基底(MTO)の数、平面波基底の数、全基底の数である。

「rsibl1,rsiblp」の直後のrsibl2で実空間で波動関数の積の和をとり電子密度を作っている。この部分から波動関数の第0成分を取り出すことができる。

bndfp内のaddrblの出力変数srout_zvに、第0成分の電子密度が入る(実空間メッシュ)。

--TakaoKotani 01:38, 27 April 2011 (JST)


Fermi levelを切るバンドまでしか固有ベクトルが入っていない。計算しているのはmkewgt()。

addrbl()

  if (lwtkb .eq. 0) then
    call mkewgt(...,nevec,...)  # output nevec is number of bands to the EF
  else
    ...
  endif

      if (lrout .gt. 0) then
         call rsibl ( ssite , sspec , slat , lfrce , nbas , isp , q ,
     .   iq , ndimh , nspc , napw , igapw , iv_p_oidxsh , numq , nevec
     .   , evec , ewgt , k1 , k2 , k3 , smpot , smrho , f )

fp/lmfp.F

      ndham=sham%ndham

      call bndfp ( nbas , nsp , nlibu , lmaxu , lldau_iv , ssite 
     ., sspec , slat , sctrl , sham , spot , sbz , sstrn , ndham , 
     . leks , lrout , lfrce , lpnu , dmxp , iter , maxit  , !evl_rv, 
     . frc_rv , dmatu, vorbdmat, llmfgw )

ndhamはどこから?

suham()

      call iinit(ldham,16)
      if (mod(pwmode,10) .eq. 2) then
        call info0(20,0,0,' suham: LMTO basis will be excluded')
      else
           call makidx ( nl , nkaph , 1 , nbasp , 0 , sspec , iv_p_oips 
     .     , iv_p_offh , iv_p_oidxsh , ldham ) 

     sham%ndham=ldham(1)

意味不明...。

--bandの時

--bandでこれを回避できるだろうか。

      if (cmdopt('--band',6,0,strn)) then
        plbnd = 1
        plbopt = strn(7:)
        lrout = 0
        lfrce = 0
        nkp = 0
        numq = 1
        allocate(ifbls_iv(ndham*nspc*2))
        ifbls_iv(:)=0

      allocate( evlall(ndham,nsp,nkp))  # nkp=ではまずいかも

wkp.{name}からfermi level (ef0)を読む。

      #nkp=0で呼ぶ。
           call suqlst ( plbopt , iopq , ndhamx , ef0 , i , iwdummy , nfbn
     .     , ifbls_iv , nkp , qp , onesp ) 
        call mpibc1(nkp,1,2,mlog,'bndfp','nkp')
... 他にも色々ある。
      allocate(qplist(3,nkp))
      do  iq = 1, nkp
        if (plbnd .ne. 0) then
          i = nsp
          if (onesp .ne. 0 .or. nspc .eq. 2) i = 1
           call suqlst ( plbopt , 0 , ndhamx , ef0 , i , iwdummy , nfbn 
     .     , ifbls_iv , nkp , qp , onesp ) 
        else
           call dpscop ( rv_p_oqp , qp , 3 , 3 * iq - 2 , 1 , 1d0 ) 
        endif
        qplist(:,iq)=qp
      enddo

      print *,'-------- qplist --------' 
      do iq=1,nkp
        write(6,"(i5,3f8.3)")iq,qplist(:,iq)  
      enddo
  # nkp, qplist はここから取れる。
 
            nmx = min(nevmx,ndimhx)
            nmx = max(nmx,nevec) !
bndfp 0: nevmx,ndimhx,nmx,nevec=    0   22    0    0
...
bndfp 1: nevl,nmx,nev=   22   -1    0  
# nmx==0なのでeigen vectorは求めない。
# call zhev_tk()
# 求めたeigen vectorの数(nev)は当然0
bndfp 2: nevl,nmx,nev=   22   -1    0

fat band

C       Need all eigenvalues if 'fat bands' plotting mode
            if (nfbn(1) .gt. 0) then
              nmx = ndimhx
              efmax = 99999
            endif

colorを指定すると'fat band'になる。

lmf --band~col=1:8~fn=syml fe 
bndfp 0: nevmx,ndimhx,nmx,nevec=    0   23    0    0
bndfp 1: nevl,nmx,nev=   23   23   23
bndfp 2: nevl,nmx,nev=   23   23   23

このモードでeigen vectorがすべて求まる。 しかし、plbnd==0だとaddrbl()に行かないので別にコードが必要。

rslbl()の構造

hsibl1()にてFT of smoothed Hankel functionが計算される。 H_L(q)=\frac{-4 \pi}{\epsilon - q^2} \mathcal Y_L(-iq) e^{\gamma(\epsilon - q^2)}
。 変数he,hrが返される。 he=1/(e-q^2) and hr=exp(-gamma q^2), g=(1/4)*R_{MT}^2. \mathcal Y_L(r)=r^l Y_L(\hat r). exp(gamma e)と\mathcal Y_L(-iq)はrsibl5で計算され掛けられる。

rsibl()

    if (nlmto .gt. 0) then
         call hsibl1 () # calc. factors for G space smoothed Hankel
    endif
    call rsibl1() # calc wavefunction in the G space, add smoothed Hankel in the G space 
    call rsiblp() # add PW part in the G space
    call rsibl2() # FFT to get realspace wavefunction, and calc. <psi|psi> to get the R space charge

rsibl1()

    call dpzero(psi, 2*ng*nspc*nvec)
     do io=1,norb
       call rsibl5() # calc. factors of smoothed Hankel in the G space
     enddo
     call rsibl6() # add PW wavefunction

rsibl5()

      fac = -4d0*pi*dexp(e*rsm*rsm*0.25d0)/vol   # -4 pi exp(e gamma)
      cc = (0d0,1d0)*fac                         # phase factor of {\mathcal Y}_L(-iq)
      do  l = 0, lmax
        cc = cc*(0d0,-1d0)
        do...
          cfac(ilm) = cc
        enddo
      enddo
      do  ispc = 1, nspc
        do  ilm = nlm1, nlm2
          ii = ilm-nlm1+ioff+1
          do  iv = 1, nvec
            evp(ilm,ispc,iv) = evec(ii,ispc,iv)*cfac(ilm)
          enddo
        enddo
        do  ilm = nlm1, nlm2
          do  iv = 1, nvec
            do  i = 1, ncut
              psi(i,ispc,iv) = psi(i,ispc,iv)+ylw(i,ilm)*evp(ilm,ispc,iv)
            enddo
          enddo
        enddo
      enddo

rsibl6()

      do  iv = 1, nvec
        do  ispc = 1, nspc
          do  i = 1, ng
            psi(i,ispc,iv) = psi(i,ispc,iv)
     .      + psi0(i,ispc,iv)*dcmplx(cosgp(i),singp(i))  # phase factor exp(ikR_spc)
          enddo
        enddo
      enddo

rsiblp()

      do  ispc = 1, nspc
        do  i = 1, nvec
          do  igv = 1, napw
            psi(ivp(igv),ispc,i) = psi(ivp(igv),ispc,i)
     .      + evec(nlmto+igv,ispc,i)/sqv
          enddo
        enddo
      enddo

rsibl2()

      do  ispc = 1, nspc
        do  i = 1, nev
          call gvputf(ng,1,kv,k1,k2,k3,psi(1,ispc,i),f)
          call fftz3(f,n1,n2,n3,k1,k2,k3,1,0,1)

                  smrho(i1,i2,i3,iq,ispc) = smrho(i1,i2,i3,iq,ispc) +
     .            wgt1*dconjg(f(i1,i2,i3))*f(i1,i2,i3)

        enddo
      enddo

GGAでNaNがでる部分

smooth rhoeps =         NaN   rhomu =         NaN  avg vxc =         NaN 

をだしているのはsmvxc2()


smvxcm: negative smrho_w number,min(smrho_w)=       31422 -2.23817234277943872E-006

を出しているのはsmvxc()



 avg es pot at rmt=-1.572817  avg sphere pot= 0.000099  vconst= 1.572817

 smooth rhoves      6.601484   charge     1.989995
jwe0245i-e In DLOG(dx) or LOG(dx) or DLOG10(dx) or LOG10(dx) or DLOG2(dx) or LOG2(dx), dx.le.0.0 (dx=-nan                   ).
 error occurs at m_xcpbe.xcpbe_                  line 992 loc 00000000006835b6 offset 00000000000089e6 
 m_xcpbe.xcpbe_                      at loc 000000000067abd0 called from loc 00000000005057b5 in vxcnlm_      line 160 
 vxcnlm_      at loc 00000000005032a0 called from loc 00000000005015a7 in smvxc2_      line 506 
 smvxc2_      at loc 00000000004ffc60 called from loc 00000000004a6fa8 in mkpot_       line 406 
 mkpot_       at loc 00000000004a6080 called from loc 000000000042ac23 in bndfp_       line 840 
 bndfp_       at loc 00000000004269a0 called from loc 000000000040a55e in lmfp_        line 856 
 lmfp_        at loc 0000000000405c00 called from loc 0000000000404ccb in MAIN__       line 439 
 MAIN__       at loc 0000000000402260 called from o.s.  
taken to (standard) corrective action, execution continuing.
jwe0245i-e In DLOG(dx) or LOG(dx) or DLOG10(dx) or LOG10(dx) or DLOG2(dx) or LOG2(dx), dx.le.0.0 (dx=nan                    ).
 error occurs at m_xcpbe.xcpbe_                  line 1067 loc 00000000006838f9 offset 0000000000008d29 
 m_xcpbe.xcpbe_                      at loc 000000000067abd0 called from loc 00000000005057b5 in vxcnlm_      line 160 

iterationの構造

fp/lmfp.F:

   integer:: iter  # iteration
   integer:: maxit # max. iteration

      iter = 1
    5 continue    #------------- loop start

     call bndfp(,iter,)

     if (cmdopt('--wden',6,0,strn)) then
       call ioden() # ioden is called everytime!
     endif

     call iors()

     call nwit(,iter,...,lsc)  # check SCF is converged or not, return=lsc

        iter = iter+1
C     Continue iterations toward self-consistency
        if (lsc .gt. 2) then
          irs(1) = 3
          goto 5          #----------- loop end
        endif

nwit

Co   lsc   :0 self-consistency achieved (diffe<=etol, qdiff<=qtol)
Co         :1 if not self-consistent, but encountered max. no. iter.
Co         :2 Harris energy from overlap of free atoms (iter=1 and lhf=t)
Co         :3 otherwise

smrho関係メモ

(構造体やソースが古いものを元にしている。)

smrho=> pot~osmrhoである。

  • bndfp-mkpot: \tilde nからpotentialをつくる。smrhoはvalence chargeのsmoothed partのはずで、mkpotにinputとして渡ってきている。smrhoはどこで値をいれているのか・・・。
  • bndfp-mkpot-rhomom: 多重極を計算
    • rhomom-corprm: coreのG,Hを計算,Zを取る。
    • rhomom-pvrhom: valenceのqmom=multiple moment=各LのMT内の電荷量を計算。L=0はcoreの寄与が入る。
  • bndfp-mkpot-smves: \tilde n_0からFFTを用いてCoulom potentialを求める。MT上のpotentialの値vval(ilm)はplanewaveをspherical Besselで展開する式eq(pw-sbessel)を用い求める。
  • sm*が\tilde n_0に関する量。
  • locpot-locpt2: 各サイトのpotentailを境界条件はvval(ilm)とcharge zとしてPoission eq.を解く。


  • fp/smvxcm

w(ocgh1)=from smcorm (foca hankel heads)
smrho = smooth valence density on uniform mesh
smrhoの規格化はcharge=sum_i1,i2,i3 smrho(i1,i2,i3)*vol/(n1*n2*n3)

w(osmrho)=w(ocgh1) -> FFT
L113足しているL115でcopyしている。 n_0^core +n_0 = w(osmrho)= w(osmrho)+smrho (係数は1であっている)
ここでn_0^core=n_sH^core_0

L150: smvxc2(w(osmrho))

  • fp/smvxc-smvxc2

L310: rho=2*smrho L326: evxcv(rho,...)


  • fp/locpt2

rho1: local valence density = sum_ilm rho1(ilm) * r**2 Y_L(ilm), on radial mesh
rho2:local smoothed valence density, defined as rho1, local atomic valence density is rho1-rho2
rhoc:core density times 4*pi*r*r

y0=q/sqrt(4pi)
L637:
n_1+n^c_1 = rhol1 = rho1 + y0*rhoc

L811: dpadd(rho2(1,1,isp),rhochs,1,nr,y0/nsp)
rho2= rho2 + y0/nsp*rhochs
としてから L813: call vxcnsp()としてpotentialへ。
その後
L816: call dpadd(rho2, rhochs, ..., -y0/nsp)
この後はrho2=rho2-y0/nsp*rhochs
同じ変数で足したり引いたりして使いまわしているので注意。


  • fp/vxcnsp(rl), rl:full charge density * r**2

L106: vxcns3(rl), rl=rl/ri**2
L113: vxcns2(rl)

  • fp/vxcnsp-vxcns2(rl)

L291: dgemm('N','T',nr,np,nlm...), rps(spin)=rl(spin)*yl
L292: daxpy, rp=rp+rps(spin) (rp=spin sum)
L354: evxcv(rp=rho spin sum,rps=rho per spin), smooth partと同じroutineを呼んでいるのでsmrho+factor*(rho1-rho2)のfactorは1
(rl(nr,nlm,nsp), yl(np,nlm), rps(nr,np,nsp))


fp/ioden.F

 47 Ci         :   core=#   specifies how local rho is to be included
 48 Ci         :            #=0 include core densities - nuclear charge
 49 Ci         :            #=1 include core densities
 50 Ci         :            #=2 (default) exclude core densities
 51 Ci         :            #=-1 no local densities to be included
 52 Ci         :            #=-2 true local density, no smoothed part
 53 Ci         :            #=-3 istl-local sm densities, no true part
     # mode=i
         call rhgcmp ( i , 1 , nbas , ssite , sspec , slat , sv_p_orhat
     .   , kmax , ng , cn_zv )

rhgcmp()

Ckino corprm() returns parameters for smooth core+nucleus representation
Ckino    the pseudocore density is
Ckino      cofg*g0(rg;r)*Y0 + cofh*h0(rfoca;r)*Y0        (1)

        call corprm(sspec,is,qcorg,qcorh,qsc,cofg,cofh,ceh,lfoc,rfoc,z)
Ckino Makes exp(-i p* (q+G)) for a list of reciprocal lattice vectors
Ckino qlat(3) * iv(ng,3) = g(ng,3)
        call suphas ( q0 , tau , ng , iv_iv , n1 , n2 , n3 , qlat 
     .  , cs_rv , sn_rv )
Ckino G_kL expansion of valence sphere densities
        call rhogkl ( ib , ib , nsp , modgkl , ssite , sspec , sv_p_orhoat
     .  , kmax , qkl_rv )
Ckino Convert G_kL expansion of function centered at a site to PW's
Ckino add Gaussian(G) contribution to cg
Ckino and, if lcor= mod ( mode  / 10 , 10 ) .ge.2, add H(G) contribution
        call rhgcm2 ( vol , rg , rfoc , ceh , cofh , kmax , mod ( mode 
     .  / 10 , 10 ) .ge.2 , qkl_rv , nlm , ng , g2_rv , yl_rv 
     .  , cs_rv , sn_rv , cg )
Ckino add Z * delta(r) in the G space
        if ( mode .ge. 200 ) call rhgcm3 ( - z , vol , ng , cs_rv 
     .  , sn_rv , cg )

rhogkl()

  #--- parent subroutine calls corprm. Calling the same subroutine again...
  call corprm(sspec,is,qcorg,qcorh,qsc,cofg,cofh,ceh,lfoc,rfoc,z)
  call pvrgkl ( mode , kmax , nlml , nr , nsp , rofi_rv , rwgt_rv 
     .   , sv_p_orhoat( 1 , ib )%v , sv_p_orhoat( 2 , ib )%v , sv_p_orhoat( 3 , ib )%v 
     .   , h_rv , cofh , rg , ceh , rfoc , z , qkl ( 0 , j1 ) ) 

rhgcm2()

      pi = 4d0*datan(1d0)
      y0 = 1d0/dsqrt(4d0*pi)
      gam = 0.25d0*rg*rg
      gamf = 0.25d0*rfoc*rfoc
      cvol = 1d0/vol
      cfoc = -4d0*pi*y0/vol
      do  i = 1, ng
        phase = dcmplx(cs(i),sn(i))
        aa = dexp(-gam*g2(i))*cvol
        cc = aa*phase*(0d0,1d0)
        ilm = 0
        do  l = 0, lmxl
          cc = cc*(0d0,-1d0)
          do m = -l,l
            ilm = ilm+1
            fac = 1
            sqkl = 0
            do  k = 0, kmax
              sqkl = sqkl + qkl(k,ilm)*fac
              fac = -g2(i)*fac
            enddo
         #------ This part calculates Gaussian(G)*exp(i pos G)
            cg(i) = cg(i) + sqkl*cc*yl(i,ilm)
          enddo
        enddo

        if (lcor) then
           #----- this calculates H(G)*exp(i pos G)
          aa = cfoc*dexp(gamf*(ceh-g2(i)))/(ceh-g2(i))
          cg(i) = cg(i) + cofh*aa*phase
        endif

      enddo

corprm

Cr   qcorg and qcorh are the charges in the gaussian and hankel parts.
Cr   The hankel part is used when the core is allowed to spill out of
Cr   the augmentation sphere.
Cr
Cr   cofg and cofh are the coefficients in front of the standard
Cr   gaussian and smoothed hankel functions for l=0.
Cr   That is: the pseudocore density is
Cr      cofg*g0(rg;r)*Y0 + cofh*h0(rfoca;r)*Y0        (1)
Cr   ceh and rfoc are the energy and sm.-radius for the hankel part.
Cr   cofg is set so that qc = integral of eq. 1 above.
Cr
Cr   For lfoc=0 there is no Hankel part; qc carried entirely by Gausian
Cr   For lfoc>0 there is no Hankel part; Gaussian carries difference
Cr              between qc and charge in Hankel part.

 どういうアルゴリズムで決めているのだ?

pvrgkl

Ci   mode  : a compound of digits specifying what is to be included
Ci         : in the expansion coefficients
Ci         : 1s   digit = 1 include local density rho1-rho2
Ci         :              2 include local density rho1
Ci         :              3 include local density rho2
Ci         : 10s  digit = 1 include core density rhoc
Ci                        2 include -1 * core density from sm-hankel
Ci                        3 combination 1+2
Ci         : 100s digit = 1 add -1 * nuclear density Y0 Z delta(r)

#valence part?
                  qkl(k,ilm) = qkl(k,ilm) + rwgt(i) * pkl(i,k,l) *
     .            (f1*rho1(i,ilm,isp) + f2*rho2(i,ilm,isp))  # if mod(mode,10)==1, f1=1,f2=-1
#core part
C         Case 1 or 3: add core density
            if (mod(mod(mode/10,10),2) .ne. 0) then
              do  i = 1, nr
                qkl(k,1) = qkl(k,1) + y0*rwgt(i)*rhoc(i,isp)*pkl(i,k,0)
              enddo
            endif
C         Case 2 or 3: subtract core density from sm. Hankel
            if (mod(mode/10,10) .ge. 2) then
              do  i = 1, nr
                call hansmr(rofi(i),ceh,1/rfoc,xi,1)
                smrch = cofh*xi(0)*rofi(i)**2
                qkl(k,1) = qkl(k,1) - rwgt(i)*smrch*pkl(i,k,0)
              enddo

smshft

  call pvsms1()

pvsms1()

 if (mod(job,10) .eq. 1) then

   do 14 i=1,nsp
    #  H(G) を求めてcg(ig,i)に足している。
   14 continue
   do 16 i=1,nsp
    # Gaussian(G)を求めてg(ig,i)に足している。
   16 continue

 elseif (job .eq. 12) then

           call rhgcmp ( 131 , ib , ib , ssite , sspec , slat , sv_p_orhoat 
     .     , kmax , ng , cwk ) 
      # ? 同じことを2度やっているように見える。

fp/ioden,古いコードを元にしているメモ

G空間でMT chargeを足しているのでvalence,coreのpeakyなchargeは表すことができない。そのためcore optionを変えても結果がcoreで変化がないように見えることがある。 また、下で書くがcore<0の場合は期待した動作をしない。

 47 Ci         :   core=#   specifies how local rho is to be included
 48 Ci         :            #=0 include core densities - nuclear charge
 49 Ci         :            #=1 include core densities
 50 Ci         :            #=2 (default) exclude core densities
 51 Ci         :            #=-1 no local densities to be included
 52 Ci         :            #=-2 true local density, no smoothed part
 53 Ci         :            #=-3 istl-local sm densities, no true part

L147-L149:
input: smrho(r)
local variable: w(opsrho)(r)= smrho(r)のコピーをFFTしたもの(gvgetfしていない)。
local variable: w(ocn)(G)

L160: core=-2: w(ocn)(G)=0
L161: core=-3: w(ocn)(G)=-w(ocn)(G)

L162: core<=-2 or core>=0: rhgcmpのcore modeの動作。 w(ocn)(G)に加える。G空間で足している!

  • core=0:add -1 * nuclear density Z delta(r) (Z is smoothed into the G_kL), add core density rhoc, add -1 * core density from sm-hankel in the local density, restoring it by adding the sm-hankel to the FT mesh, add local density rho1-rho2
  • core=1: add core density rhoc, add -1 * core density from sm-hankel in the local density, restoring it by adding the sm-hankel to the FT mesh, add local density rho1-rho2
  • core=2: add local density rho1-rho2
  • core=-2: add local density rho1
  • core=-3: add local density rho2

L163: core=-3: 再びw(ocn)(G)=-w(ocn)(G)

L169: core>=0: gvputfでw(ocn)(G)をw(opsrho)に(FFTの前処理)
L188: fftz3(), w(opsrho)(G)をFFT, w(opsrho)(r)に

L190: rlse(ocn)

L195: ioden2( w(opsrho) )

最後に使っているのはw(opsrho)。 core<0ではw(ocn)にデータを入れているがFFTの前にw(opsrho)へコピーしていないためデータを捨てているので期待する動作をしていない。(バグだろう。) core<-3にするとsmrhoそれ自体を出力するはず。

GWの波動関数

GWではoverlapと 
\langle i|v(q)|j\rangle
を用いる。 LDA基底,GW基底間の変換matrixをつくり LDA基底での値をGW基底へと変換してvの行列要素を求めている。

GWの側では波動関数をinterstitial領域はinterstitial planewave (IPW)で、MT内は\phi,\dot \phiとで表す。

IPW

IPWは IPW= PW - P(PW) と定義される。PはMT内での球Bessel関数とYLを用いたplanewave(PW)のprojectionとする。 IPWはMT内では0となる。(MT内のprojectionで見かけ出てくるangular momentum展開はmatrix elementを計算するとL一つだけ残るのでantular momentum展開の最大値は意識しなくてよい。)

PMTの側ではinterstitial regionはsmoothed HankelとPWを用いて表されていたため、 smoothed Hankelの方はPWに直す式 (<smoothed Hankel | PW >) があるのでtail,headともそれを用いてIPWに変換する。

これらを用いて 
\langle IPW|IPW \rangle, 
\langle IPW|{\rm smoothed Hankel}\rangle のmatrix elelementを計算し変換matrixを決めておく。 PMTではMTが重なっていても原理的に問題なかったが、IPWを使うGWの側ではMTは重なってはいけない。

ソース

gwd/pwmat.F, pwmatが変換を行っている。

smoothed Hankel functionの公式

Ref.[1] に記述がありますが、J. Math. Phys.はとりにくいのでsummaryを書いておきます。

spherical harmimonics

定義式が多いのですが・・・。 Y_L(\hat r)をreal valued spherical harmonicsとして


{\mathcal Y}_L(r) = r^l Y_L(\hat r)


{\mathcal Y}_L(r)
はrealのx,y,zの関数になる。


Y_K(\hat r)Y_L(\hat r) = \sum_M C_{KLM} Y_M(\hat r)

ここで、 real Gaunt or Clebsh-Gordan係数CKLMであり


 C_{KLM} = \int_0^{2\pi}\! d\theta sin(\theta)\int_0^{\pi}\! d\phi
 Y_K(\hat r) Y_L(\hat r) Y_{M}(\hat r)

と計算される。

{\mathcal Y}_L(r)の方は


{\mathcal Y}_K( r) {\mathcal Y}_L(r) = \sum_M C_{KLM} (x^2+y^2+z^2)^{k+l-m} {\mathcal Y}_M(r)

となる。

generalized Gaussian

ε = − κ2としてgeneralized Gaussianを


G_L(r)={\mathcal Y}_L(-\nabla) g(r)

と定義する。


g(r)=\left(\frac{a^2}{\pi}\right)^{3/2} \exp(\gamma \epsilon) \exp( -a^2 r^2)

であり、


G_L(r)=\left(\frac{a^2}{\pi}\right)^{3/2} \exp(\gamma \epsilon) 
(2 a^2)^l
\exp( -a^2 r^2) r^l Y_L(\hat r)

となる。

exp(γε)は Ref[1]で便宜のため導入している。

unsmoothed and smoothed Hankel

unsmoothed hankel {\check H}_L


(\Delta + \epsilon){\check H}_L(\epsilon; r) = -4\pi D_L(r)


 D_L(r)={\mathcal Y}_L(-\nabla) \delta(r)

で定義され smoothed Hankel HL(r)

eq(smoothed_Hankel):
(Δ + ε)HL(a,ε;r) =  − 4πGL(a,ε;r)

で定義される。なおHLは誤差関数を用いて解析的に表現可能である。

{\check H}_L は第一種Hankel関数を用いて

{\check H}_L(r)
= (-1)^l \kappa^{l+1} h_l^{(1)}(i\kappa r) Y_L(\hat r)

と書ける。 ここでε = − κ2 である。

overlap integral


 \int \! d^3 r
 H_{L_1}(\gamma_1,\epsilon_1;r-R_1)^*  
 G_{L_2}(\gamma_2,\epsilon_2;r-R_2) =
 \exp(\gamma_2 (\epsilon_2-  \epsilon_1 )) (-1)^{l_1}
\sum_M C_{L_1 L_2 M} H_{(l_1+l_2-m)/2,M}(\gamma_1+\gamma_2,\epsilon_1;R_1-R_2)

expansion of a smoothed Hankel function around a site

HL'を別の中心で展開する。

定義が多いが、


G_{kL}(\gamma;r) = \phi_{kl}(r) \left( \frac{a^2}{\pi} \right)^{3/2}
  \exp(-a^2 r^2) r^l Y_L(\hat r)


 \phi_{kl}(r) = (-1)^k ((2a)^2)^{k+l} 2^k k! L_k^{l+1/2} (a^2 r^2)

Laguerre polynomialは

 eq(Laguerre-poly)

 L_n^a(x)= \frac{e^x x^{-a}}{n!} \frac{d^n}{dx^n} ( e^{-x} x^{n+a} ) 
 =\sum_r^{n} (-1)^n\left( \begin{array}{c} n+a \\ n-r \end{array}\right) \frac{x^r}{r!}

と定義される。 PkL


 P_{kL}(\gamma;r)=p_{kl}(\gamma;r) Y_L(\hat r)


 p_{kl}(\gamma;r) = \phi_{kl}(r) \frac{ a^l (2l+1)!!}{ (2a^2)^{k+l} (2k+2l+1)!!} r^l

と書ける。 これを用いて

HL'(γ',ε';rR') = akL(γ,RR)PkL(γ;rR)
kL

を展開する。

展開係数は

eq(Hankel_expansion):

 a_{kL}(\gamma,R-R') =
\frac{4\pi}{ (2a^2)^k a^l 2^k k! (2l+1)!!} \int\! d^3 r
  G_{kL}(\gamma;r-R) H_{L'}(\gamma',\epsilon';r-R') 


 = \frac{4\pi}{ (2a^2)^k a^l 2^k k! (2l+1)!!}
 \exp(-\gamma \epsilon') (-1)^l 
\sum_M C_{LL'M} H_{(l+l'-m)/2+k,M}(\gamma+\gamma',\epsilon';R-R')

となる。

G_kLとP_kL

PkLGkLとは互いにdual basisである。

 \int \!d^3 r G_{kL}(\gamma;r) P_{k'L'}(\gamma;r) =
 \frac{1}{4\pi} (2a^2)^2 a^l 2^k k! (2l+1)!! \delta_{kk'} \delta_{LL'}

eq(exp-Laguerre):

exp(iqr) = ckLPkL(r)
kL

という展開も解析的に求められる。 展開係数C_kLは


\int \!d^3 r \exp(iqr) G_{kL}(r) = c_{kL} \int\! d^3 r P_{kL}(r)G_{kL}(r)

なので左辺はG_kLのフーリェ変換\hat G_{pL}(q)であり、解析的に


\hat G_{pL}(q)=(-q^2)^p{\mathcal Y}_L(-iq) \exp(\gamma(\epsilon - q^2))

と求められる。

[1] E. Bott et al. J. Math. Phys. 39, 3393 (1998).

Personal tools