Wannier function

From Ecal

Jump to: navigation, search

このページの内容はKinoのメモの段階です。 programはMiyakeさん & Kotaniさんにより作られています。

Contents

Definitions

ψk: Bloch function. Blochの定理によりψnk(r + T) = exp(ikTnk(r), where Tはtranslational vector.

uk: periodic part of ψk. unk(r) = exp( − ikrnk(r)

FBZ: first Brillouin zone. 総数はnqbz, vectorはqbz(1:3,iqbz)。

bb: =b。総数はnbb, そのvectorをbbv(1:3,ibb)とも書く。

tight binding modelでのdispersionを書くk点: <QPNT>(もしくは、最新版ではSYMLを優先)で指定する。q(1:3,1:nq)

rcut: wannier関数を作った後のtight binidng modelで rnn'(nn'間の距離) > rcutの場合Hnn' = 0とする。

exp(ibr)の評価

  • in MTs: 波動関数が原子を中心に展開しているため、各原子を中心にexp(ibr) = sum of Bessel functionsとして各MT内毎に計算。
  • interstitial region: < k | exp(ibr) | k + b > を直接計算

動かしかた

GWの途中

...
### Valence part of the self-energy
echo 0|$nfpgw/hbasfp0 >lbas

hbasfp0(mode0)の後で

### maxloc
echo 1|$nfpgw/hmaxloc   >lmaxloc1
$nfpgw/hpsig            >lpsig
echo 2|$nfpgw/huumat    >luumat2
echo 2|$nfpgw/hmaxloc   >lmaxloc2

として実行する。

GWinput 090910 version

https://github.com/nim-hrkn/maxloc090910

<MLWF>
3          # gaussian, nwf
1 1 1  # nphi(1:nwf)
9 14  2.0 1.0 phi,phidot and lamda(angs) of gaussian 1, #iphi(j,i),iphidot(j,i),
r0g(j,i),wphi(j,i)
10 15 2.0 1.0 phi,phidot and lamda(angs) of gaussian 2
12 17 2.0 1.0 phi,phidot and lamda(angs) of gaussian 3
</MLWF>

parameter in GWinput

これではわかりにくいですが、ないよりましなので書きます。

  • wan_out_ewin: boolean, energy outer window ?
  • wan_in_ewin: boolean, energy inner window?
  • wan_in_bwin: boolean, band idnex upper inner window?

if wan_out_ewin=true

  • wan_out_emin: real, energy lower outer window
  • wan_out_emax: real, energy upper outer window

else

  • wan_out_bmin: integer, energy lower band index window
  • wan_out_bmax: integer, energy upper band index window

endif

if wan_in_ewin=true

  • wan_in_emin: real, energy lower inner window
  • wan_in_emax: real, energy upper inner window

else if wan_in_bwin=true

  • wan_in_bmin: real, energy lower band index window
  • wan_in_bmax: real, energy upper band index window

endif

select Hilbert space [2]

  • wan_maxit_1st: real,default=100,
  • wan_conv_1st: real, default=1e-5
  • wan_mix_1st: real, default=0.1, α

find max. localized wannier function [1]

  • wan_maxit_2nd: real,default=100,
  • wan_conv_2nd: real, default=1e-5
  • wan_mix_2nd: real, default=0.1, α

ixc=2

tight binding Hamiltonian dispersion, output= bnds.maxloc.(up|dn)

ixc=2

tight binding Hamiltonian dispersion, output=bnds.tb.(up|dn)

  • wan_tb_cut: real, defaul=1.01d0, H_{nn'}(R)=0 if R>thisvalue in alat unit

ixc=3

tight binding Hamiltonian dispersion, output=bnds.cut.(up|dn)

  • wan_tbcut_rcut: real, default=wan_tb_cut, H_{nn'}(R)=0 if R>thisvalue in alat unit
  • wan_tbcut_heps: real, default=0, H_{nn'}(R)=0 if |H|< thisvalue

misc

tight binding Hamiltonian dispersion, output=bnds.tb.(up|dn)


  • wan_nb_below: real, default=0
  • wan_nb_above: real, default=0


tight binding Hamiltonianを書くときのk-path

  1. SYML を読む (lmfで用いたのと同じformat)
  2. SYMLがなければQPNTを読む

old GWinput

<MLOC>
3          # gaussian, nwf
1 1 1  # nphi(1:nwf)
9 14  2.0 1.0 phi,phidot and lamda(angs) of gaussian 1, #iphi(j,i),iphidot(j,i),
r0g(j,i),wphi(j,i)
10 15 2.0 1.0 phi,phidot and lamda(angs) of gaussian 2
12 17 2.0 1.0 phi,phidot and lamda(angs) of gaussian 3
0 16 18      ieo_swt,itin_i,itin_f
-10.0 5.0 eomin,eomax
0 0 0         iei_swt, itin_i,itin_f
-0.0 0.0   eimin,eimax
1000        nsc1
1d-5       conv1
0.7        alpha1
1000        nsc2
1d-12      conv2
1.1        alpha2
2.5        rcut (in units of alat)
#GaussianHead on
TailTruncate on
</MLOC>

nwf: wannier orbital の数 = バンドの数

nphi: 各wanneir orbitalの初期値のatomic orbitalの数 \Psi = \sum_i wphi \; \psi_i, \psi_i=c_1 \; phi + c_2 \; phidotψi(r = MT)においてtail のgaussian Nexp( − (r / r0g)2)(Nは規格化定数) にあわせるという接続条件からc1,c2が自動的に決まる。 (phi=φ,phidot=\dot \phi、上ではprogramの表示に合わせている。)

itin_i,itin_f: outer band のindex. eomin,eomaxとenergyとしても指定可能。ieo_swtでどちらか指定。ieo_swt=1(energy window), =0(band index window)

itin_i,itin_f: inner bandのindex. eimin,eimaxとenergyとしても指定可能。iei_swtでどちらか指定。iei_swt=1 (innder energy window) =2 (inner band index window), otherwise (off)

nsc1, conv1, alpha1 : hilbert spaceを求めるloop 1の最大回数と収束条件とα (see [2])

ncv2, conv2, alpha2 : wannier robitalを求めるloop 2の最大回数と収束条件とα (see [1])

rcut: tight binding Hamiltonian を書くときのR cutoff

GaussianHead on|off : initial orbitalとしてMT内をつかうかどうか。使わなければtailだけになる。defaultでは「使う」。

TailTruncate on|off: initial orbitalとしてMT外を使うかどうか. on=使わない。

他にもある。

コードの中身

GWinputの<QPNT>のk点は[2]のeq.(26)のk'点に相当。つまりWannier関数で'fit'したdispersionを描くときのk点。rotateした全Hnn'を用いると全く同じになるはずだがHnn'のrangeをtruncateすると全く同じにはならない。(最新版ではDFTでdispersionを書かせる時に用いたSYMLファイルがあるときはSYMLがこのk'点として読まれる。)

wannier functionを作るときのk点はこのk'点とは違う。wannier functionを作るときはGWinputのn1n2n3が用いられ、FBZのk点が自動的に読み込まれる。

FBZの点

BZDATAに書いてあり、読み込むのはread_BZDATA(gwsrc/rwbzdata.F)

FBZの総数はnqbzでdata consistencyをcheckするためにいくつかのファイルに書かれている。 例えば、hbe.dの5つめ、BBVECの2行目。BZDATAの1行目。 (data構造が繁雑なので簡単にするよう提案中。)--Kino 11:26, 13 October 2009 (JST)

prepare bb vector, hmaxloc.F, mode=1

  1. getbb(maxloc0.F)でbbの計算。
  2. writebb(maxloc0.F)でfile=BBVECにbbの情報を書き出す。

< ψk | g > , hpsig.F

g=initial wannier center.

k点はhpsig.Fにおいて

     do 1070 iqbz = 1,nqbz
           q1(:) = qbz(:,iqbz)

として回っている


gの初期値

iphi, iphidot, g0r, wphiを与える。 MT内で \psi_i(r)= c_1 \phi_{iphi}(r)+ c_2 \dot\phi_{iphidot}(r) としてindexを与える。これで原子位置と量子数を決める。(原子とindexの関係はどこで見ればいい?)

call getc2c2(...)

c1,c2の決め方はr0g(atomic unit)を用いて gaussian(r) = Nexp( − (r / r0g)2) (Nは規格化定数) というgaussianのradial directionのひろがりからMT上rMTψi(rMT) = gaussian(rMT)とその一回微分が連続という条件で 計算される。ここでのgaussianの中心はiphi,iphidotで指定された原子上である。

g(r) = wphiiψi(r)
i

call qggmat(...)

MT外は動径方向はgaussian(r)そのもの。(角度依存性はiphiのものが用いられる。)

input parameter

GaussianHead boolean (default=false) true=MT内を使わない。
TailTruncate boolean (default=false) true=MT外を使わない。

軌道の向きは 1 s, 2 p_y, 3 d_z, 4 p_x, 5 d_xy, 6 d_yz, 7 d_z^2, 8 d_xz, 9 d_x^2-y^2

さらにindexへの対応は

c l    0     1                 2
c n    1  2  1        2        1               2
c m    0  0 -1  0 -1 -1  0  1 -2 -1  0  1  2  -2  -1   0   1   2
c ind  1  2  3  4  5  6  7  8  9  10 11 12 13 14  15  16  17  18

phiは9=d_xy, 10=d_yz, 11=d_z^2, 12=d_xz, 13=d_x^2-y^2ということ。 phidotは14=d_xy, 15=d_yz,...。

< uk | uk + b > , huumat.F, mode=2

wannier connection  M_{mn}^{k,b}=<u_{m,k}|u_{n,k+b}> の評価

readbb(hpsig.F)でBBVECからbbを読む。

     call kbbindx(qbz,ginv,bb,
    d             nqbz,nbb,
    o             ikbidx,ku,kbu)

でk+bのindexを作る。 k,k+bのloopは

     do 1070 iqbz = 1,nqbz
       ...
       do 1080 ibb = 1,nbbloop ! nbbloop=nbb
           iqb = ikbidx(ibb,iqbz)
           q1(:) = qbz(:,iqbz)
           q2(:) = q1(:) + bbv(:,ibb)

として回る。k+bのindexがiqbで与えられる。uはperiodic partなのでFBZ内の点を選べばよい。

construct wannier function, hmaxloc.F, mode=2

GWinput<MLOC>のパラメタを読み込む。

wannier spread 
\Omega=\sum_n \left[\langle r^2 \rangle_n -  {r_n}^2 \right]
=\Omega=\Omega_I + \tilde \Omega
。 entangled bandのつなぎかたのambiguityを解消するため ΩIをHilbert spaceを(各kでどの{n}を)選択するために最小化する。[2] そのHilbert space内で\tilde \Omegaを最小化し (ΩIはHilbert space内のrotationでは不変。) wannier functionを局在化させる。[1]

step 1, energy window内でHilbert spaceを選択する

entangled bandの処理を行う。 Ref.[2]による。

energy windowとindexの変換。

        call ewindow(is,ieo_swt,iei_swt,itout_i,itout_f,itin_i,itin_f,...

input parameter (see Fig. 6 in [2])

  • outer energy (band index) window: eomin, eomax (itout_i,itout_f)
  • inner energy (band index) window: eimin, eimax (itin_i,itin_f)
entangled bandの処理をしない場合はouter windowのみ指定し、inner windowは用いない。
ieo_swt=0 (band index), =1(energy window)。
iei_swt=0 (inner window off), =1 (energy window), =2 (band index)。

output

  • outer energy window index: iko_ix, iko_fx, iko_i(1:nqbz), iko_f(1:nqbz)
  • inner energy window index: iki_i(1:nqbz), iki_f(1:nqbz)


loop

       do isc = 1,nsc1
  1. (1-2) <u_mk | P_k+b | u_nk>, where P_k=\sum_{n}|u_{nk}\rangle\langle u_{kn}|,call getupu(), mixing parameter: alpha1を使う
  2. (1-3) Zmn(k) -> phi,eval, call getzmn, eq.(21) in [2]
  3. (1-4) w_I(k), call get_omgik(), eq.(18) in [2]
  4. (1-5) w_I(k) -> Omaga_I eq.(11) in [2]
  5. (1-6) check self-consistency

step 2, localize Wannier fn.

Ref.[1]による。

       do isc = 1, nsc2
  1. (2-0) construct initlal u~ from u
  2. (2-1) initial: uumat -> Mmn
  3. (2-2) initial U
  4. <r_n> (eq.31 in [1])
  5. (2-3) A[R] matrix
  6. (2-4) S[T] matrix
  7. DW(k) (eq.57 in [1]), α = alpha2を用いる。
  8. update Mmn (eq.61 in [1])
  9. (2-6) Omeg_I, Omega_D and Omega_OD (eq.34,35,36 in [1])
  10. check self-consistency

How to update

initial guess gn

projection onto ψn space, |\phi_n^{(0)} \rangle = \sum_m |\psi_m\rangle\langle \psi_m | g_n \rangle =
\sum_m A_{mn} | \psi_m \rangle but, \phi_n^{(0)} is not orthogonal.

redefine first guess.  | \psi^{(0)} \rangle = S^{-1/2} | \phi^{(0)} \rangle = A S^{-1/2} | \psi_n \rangle , where  S=\langle \phi^{(0)} | \phi^{(0)} \rangle . ((22)-(23) in [2].)

|\psi_{n}^{(i)}\rangle = \sum_m | \psi_{m}^{(0)}\rangle c_{mnk}, where cmnk is used in the program.


eq.(57)-eq.(61) in [1].

M_{mn}^{(0)(k,b)}=\langle u_{nk}^{(0)}|u_{m,k+b}^{(0)}\rangle

|u_{n,k}\rangle = \sum_{m}  | u_{mk}^{(0)} \rangle U_{mn}^{k}

\Delta W^{k} = \frac{\alpha}{w} \sum_b w_b ( A[R] - S[T] )

U^{k} \rightarrow U^{k}\exp(\Delta W^{k} )

Mk,b = UkM(0)k,bUk + b

step 3, tight binding model(Hnn')を作る

input parameter, rcutでcutoffをする。dispersionを書く。do i=1,iq; q(1:3,iq)

結果の変換

nwf: band 数。

get_hrotr_wsでhrotk(nwf,nwf,nqbz)をhrotr(nwf,nwf,nrws)に変換。 eq.(25) in [2], H_{mn}(R)=\sum_k \exp(-i k R) H_{mn}(k)/N_{kp}=\langle w_{m0} | H | w_{nR} \rangle

\langle w_{m0} | H | w_{nR} \rangleの距離Rの情報はwigner_seitz()で計算して*rwsに入っている。

             rws(1:3,nrws) = rr(1:3,n) ! vector
             drws(nrws) = dd(n)  ! length
             irws(nrws) = ndegen  ! degeneracy

hrotr(nwf,nwf,nrws): rotated Hamiltonian (full range, H_nn')

get_hrotkp_wsでk空間へ変換->hrotkp

  • ifbnd= bnds.maxloc.up (or .dn), hrotkpを対角化, eval(:)に固有値が入る。
  • iftb=bnds.tb.up (or .dn), get_hrotkp_tb_wsでrcut以上のH_nn'を0にする。hrotrからhrotkpをつくる。対角化して固有値をevalに入れる。
  • ifsh=bnds.sh.up (or .dn). smalHamlitonian: hrotkpsをhrotkp(nsh1:nsh2,nsh1:nsh2)として切り取ったもの。対角化してevalsに固有値が入る。


binary dataは

  • ifmlw(1:2)= MLWU or (...D), q点とその基底のrotation matrix
  • ifmlwe(1:2)= LMWEU or (...D), q点とそのeigenvector, eigenvalue

に入る。

結果

収束性

ΩOD はすぐ収束する。ほんとかいな。

距離依存性

SVO3 (a=3.84Ang), Fermi levelを切っている3つのバンド(Vのd_xy, d_yz, d_zx相当)のH(1,1,R)のR依存性。 n_k=2(2x2x2)からn_k=10(10x10x10)まで変化させた。線は見やすくするために入れているもので意味はありません。

H(1,1,R=0)はorbitral energy。 nk=4は必要なのだろう。

別のwannier functionの作り方も可能なはずで今後の研究課題。

dispersion

赤: 全て使う(band計算結果), 緑:wannier function, H(i,j,R)=0 if (R>2.5a)

赤: 全て使う(band計算結果), 緑:wannier function, H(i,j,R)=0 if (R>3a), 青:wannier function H(i,j,R)=0 if (R>3a or abs(H)<1e-3).

赤: 全て使う(band計算結果), 緑:wannier function, H(i,j,R)=0 if (R>3a), 青:wannier function H(i,j,R)=0 if (R>3a or abs(H)<5e-4).

E=-0.08eV付近でO-X lineが斜めになっているところは5E-4以下の小さいtransfer integralが効いている。

R=3までは必要らしい。

wannier spread

それぞれのばんど(全て同じ)wannier spreadΩのk点数依存性。 k点はn1**3。

hrotrは収束しているようだが、 k点数に関してはwannier spreadが収束しているように見えない。 やや気持ち悪いがtable1, table2 [1]もkmeshの収束が遅いことを書いてある。

fplmtoでは4x4x4で1.844 Ang^2 (Miyakeさんより)なのでPMTの方がやや大きい。

TODO list

  • ψk(r) = ψk(r) * を使いwannier orbitalを必ずrealに ( (65) in [1] の説明。).

初期値の与え方の改良

rotdlmm が 座標回転によるYlmの変換行列を与える。

Y'(lm')rot = Y(lm)dlmm(m,m',l)
m

オイラー角による回転 Rz( − θ)Ry(φ)により (0.0.1)が (rcos(θ)cos(φ),rcos(θ)sin(φ),rsin(θ))になる。

input format

id=1 lm= dxy rgauss=2.2 centersite= 2  axis3= (cart 1 1 1) axis1=(cart 1 0 0)
id=2 lm= dyz rgauss=2.2 centersite= 2  axis3= (cart 1 1 1)  axis1=(cart 1 0 0)
id=3 lm= dz2 rgauss=2.2 centersite= 2  axis3= (cart 1 1 1)  axis1=(cart 1 0 0)
id=4 lm= dxz rgauss=2.2 centersite= 2  axis3= (cart 1 1 1)  axis1=(cart 1 0 0)
id=5 lm= dx2y2 rgauss=2.2 centersite= 2  axis3= (cart 1 1 1)  axis1=(cart 1 0 0)
  • id=: int
  • lm=: string
  • rgauss= real, Nexp( − (r / rgauss)2)
  • centersite= int
  • axis3=,axis1= ( (cart|abc) real real real )
  • coef=: real, 同じidがあると ({\rm coef})_1\psi_1 +  ({\rm coef})_2\psi_2+ \cdotsと軌道が作られる。

axis?n=normalized axis?, crot = axis3n, arot = axis1n − (axis1n,axis3n)axis3n, b_{\rm rot}=c_{\rm rot} \times a_{rot}として (a,b,c)で作られる軸に向けてlm=で指定される軌道を傾けて初期値をつくる。

default axis3=(cart 0 0 1), axis1=(cart 1 0 0)

@MNLA_CPHIを読んでcentersite,lmで指定される軌道を選択する。\phi_i,\dot \phi_iのindexが書いてある。

  • outside MT: Nexp( − (r / rgauss)2)
  • inside MT: c_1 \phi_i+ c_2 \dot\phi_i

としてMT半径で外と中をスムーズにつなぐようにc1,c2が内部で決められる。 規格化は全空間でNexp( − (r / rgauss)2)として取られている。

将来は rgauss=の代わりにn= intを指定して

  • outside MT: Nexp( − (r / rgauss)2)
  • inside MT: φnlm

というoptionも作りたい。

bond centerにpzをつくる場合

(cart 1 1 1)がbondの方向,site1 からsite2の方向だとする。

id=1 lm= s rgauss=.5 centersite= 1  axis3= (cart 1 1 1) axis1=(cart 1 0 0) coef=.7
id=1 lm= s rgauss=.5 centersite= 2  axis3= (cart 1 1 1)  axis1=(cart 1 0 0) coef=-.7

現状

今のところ別プログラムで提供。initial orbitalを書かせることも可能。

program sylm

It rotates orbitals.

in main3e.F

      axisc=(/1.0d0,1.0d0,0.5d0/)
      axisa=(/1.0d0,0.0d0,.0d0/)
      call axisac2rot(axisa,axisc,rot)

File:axisrotate.png

  • input: axics=\vec{c} and axisa=\vec{a'}, which has nonzero projection to \vec{a}. (\vec{a'} is a rough guess to \vec{a}.) The program calculates \vec{b} automatically.
  • output:rotation matrix(3x3)
      call rotdlmm(rot, 1, nl, dlmm) 

where nl=(Max. of L) +1. If you wantt to calculate upto d orbitals, nl=3. Upto p, nl=2.

  • input: rot and nl
  • output: dlmm

Rotated orbitals are calculated as

Ylm = Ylm'dlmm(m',m,l)
m'

The order of Ylm is 's'(l=0), 'py','pz','px'(l=1), and 'dxy','dyz','dz2','dxz','dx2y2'(l=2).


Examples

      call writemesh2(nl,lmax,lmx,dlmm,'c.cube')

which will generate cube files, from 11c.cube to 19c.cube.


  • original frame

Roughly it is expressed as File:frame_axis.png.

's','py','pz','px','dxy','dyz','dz2','dxz','dx2y2'

File:original_frame_orbitrals.jpg


  • rotated frame
      axisc=(/1.0d0,1.0d0,0.5d0/)
      axisa=(/1.0d0,0.0d0,.0d0/)

axisc is almost perpendicular to the screen.

File:rotated_frame_orbitals.jpg

References

[1] Marzari, Vanderbilt, PRB 56, 12847 (1997).

[2] Souza, Marzari, Vanderbilt, PRB 65, 035109 (2001).

Personal tools