BoltzTraP

  BoltzTraPはボルツマン方程式に従って輸送係数を求めるプログラムである。そこでは緩和時間を入力ファイルに記述しない場合、エネルギーに依存せず一定であるとした緩和時間近似が用いられている。温度によるフェルミ準位の変化は電子状態密度分布に対してフェルミディラック分布関数や他の式を用いて予めユーザー側が求めておくことが望まれる。
※ 化学ポテンシャルやゼーベック係数が計算可能。緩和時間が実験的に分かれば、電気伝導度の予測も可能になる。
------------------------------------------------------------------------------
■ Compiling
1. tar xvf BoltzTraP.tar.bz2
2. cd boltztrap-1.2.3/src
3. emacs Makefile
   FC=ifort
   FOPT = -FR -mp1 -w -prec_div -pc80 -pad -align -DINTEL_VML -traceback -axAVX,SSE4.2,SSE4.1,SSSE3,SSE3
   LDFLAGS = $(FOPT) -L$(MKLROOT)/lib/$(MKL_TARGET_ARCH) -pthread
   LIBS = -lmkl_lapack95_lp64 -lmkl_blas95_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread
4. make

■ SCF計算でのテクニック(WIEN2k)
1. 通常のk点でSCF計算を行う
2. [initialize_calc.]で "x kgen"のボタンを押して、k点数を2〜3万点に増やす
3. Terminalを開く
4. cd 計算しているディレクトリ
5. run_lapw -i 1
spinを入れた計算をしている場合は、 runsp_lapw -i 1 とする。

■ 計算方法
□  WIEN2k
SCF計算では格子定数が6Å程度でk点が20000〜30000点となるようにする。また、k点数は格子定数の長さに”逆比例”させる。
1) 計算させたい系の名前のファイルを作る。
2) 新しく作ったそのファイルをカレントディレクトリにする。
3) WIEN2kからcase.energy, case.struct コピー&ペーストする。
  case.energy は case.energyup, case.energydn, case.energyso などの場合もある。
  この他にcase.intrans が必要(README.pdfに詳細がある)。
  case.intransはCoSb3のものをコピーするのもよい。その場合は下記のようにする。
  A) EFは case.dos1 または case.output2(femiで検索すると良い)で書かれている値を入力。
  B) 電子数は case.in2 での NE (Number of electron)の値を入力。
  output2やoutput2up, output2dnにも同じ値で Number of electronが書かれている。どれでも良い。
  C) NOCALC を CALC にする。
4) $HOME/ boltztrap-1.2.3/src/x_trans BoltzTraP と入力。
  スピンの場合は、x_trans BoltzTraP -up や x_trans BoltzTraP -dn と入力)
※ CoSb3, Bi2Te3, Alについてはtestsに入力ファイルの例が記してある。ReferenceにあるCoSb3.trace や Bi2Te3.trace が良い。データを Efが同じところで比較してみると良い。大体似た結果が出ていればよい。
※ Fe2VAl を計算してみると、例えばWIEN2kでは、単位胞中に価電子帯の電子はXX個の設定になる。EFは9 eVとなる。しかし、実際の試料では、アーク炉での試料作成中にAlが0.3%程蒸発して無くなる可能性もある。そこで、Alが3価だからX電子程度が無くなったと考えてEFをシフトさせてみると、ゼーベック係数の傾向は実測値と比較的よく合うようになる(更に検証中)。
※ Fe2VAlの場合、交換相関項をハイブリッド汎関数のB1-CWにして、 Eg=0 eV にし、300Kでの抵抗の値をBoltzTraPの値で割って緩和時間を求め、改めて計算すると、実験の傾向と合うようになるかもしれない[S1]。mBJ[S2, S4]でもパラメーターが存在する。
※ WIEN2kでは交換相関項でmBJでの計算が可能[S2]、ただし、B1-CWではEg=0のバンド分散図が示されておらず、mBJ[S2]ではゼーベック係数が緩和時間を考慮して計算されていない。
※ 汎関数の傾向の違いなどを見るには文献[S3]が役に立つ。
[S1] I. Bilc and P. Ghosez, Phys. Rev. B 83 (2011) 205204.;
[S2] F. Tran and P. Blaha, Phys. Rev. Lett. 102 (2009) 226401.
[S3] D. Do et al., Phys. Rev. B 84 (2011) 125104.
[S4] M. Meinert, Phys. Rev. B 87 (2013) 045103.

□ Abinit
Abinit では BoltzTraP用のファイルを出力する prtbltztrap がある。SCF計算を1st dataset とすると、続いて下記のように2nd datasetで多くのk点での計算を行わせると良い。
prtbltztrp 1 # boltztrap output
iscf2 -2
getden2 -1
kptopt2 1 # Opiton for the automatic generation of k points, taking into acount the symmetry.
nshiftk2 4
shiftk2 0.5 0.5 0.5 # These shifts will be the same for all grids
0.5 0.0 0.0
0.0 0.5 0.0
0.0 0.0 0.5
ngkpt2 100 100 100
tolwfr2 1.0d-12
これによって、case_BLZTRAP_GEOM と case.BLZTRAP_EIGENが生成される。
References
[AB1] [CORRECTED 6.6] BoltzTrap output - prtbltztrp variable; http://forum.abinit.org/viewtopic.php?f=8&t=906

□ PWscf (PWscf の出力ファイルを pw2bgw.x を用いてBerkeleyGW に変換して計算する。BerkleyGWはUser accountで登録が必要になる)
botztrap-1.2.3のutilで、試験用のファイル Co4Sb12.nscf.out が容易されている。
1. ファイルに PWscf で出力された case.pw.out ファイルを準備する。
2. qe2boltz.py case.pw.out pw 1.0e7 0
Usage: C:\Users\\boltztrap-1.2.3\util\qe2boltz.py prefix format efermi nbnd_exclude [fn_pw [fn_energy]]

※ qe2boltz.py を wordpad などで開いて解読すると下記のようになっている。下記のargv[5]とargv[6]はargv[2]でpwを指定した場合、特に入力しなくても良い。
qe2boltz.py argv[1] argv[2] argv[3] argv[4] argv[5] argv[6]
argv[1]: prefixで case.pw.out を入力する(ファイルは case.nscf.out として、prifix に case を指定すればよいだろう)
argv[2]: pw, bands or inteqp のいずれかを選択するために指定。通常 pw、BerkeleyGWでは inteqp を指定。
argv[3]: efermi (eV) に対応する。入力した値をプログラム中でrydberg で割って処理している。1.0e6を超える場合は、prefix で指定したファイルから読み込まれる。「the Fermi energy is」または「highest occupied, lowest unoccupied level (ev): 」の行に書かれている値を用いる。後者の場合は (highest +lowest)/2/Ry としてプログラム中で処理される。
argv[4]: nband_excludeに対応する。0を指定すればよい。もし、計算時間を短くしたい場合は、バンドの番号を入力すると、低いエネルギーから入力したバンドの番号までが除かれる。
argv[5]: fname_pw または fn_pw で ディフォルトは prefix.nscf.out になっている。
argv[6]: fname_energy またはfn_energy
出力ファイルは、argv[1].intarns などとして出力される。

※pythonプログラムの変更方法
下記の部分の値を変えれば出力がそれに対応するようになる。
    deltae = 0.0005
    ecut = 0.4
    lpfac = 5
    efcut = 0.15
    tmax = 800.0
    deltat = 50.0
    ecut2 = -1.0
    dosmethod = 'TETRA'

※ τモデルへの書き換え方法(wordpadで行うと文字化けするので注意)(調査中)
f = open(fname_intrans, 'w')の前に下記を追加するとboltztrap-1.2.3のフォーマットに近づけられる?
    f_intrans += '32.0 2.0 ' + str(efermi) + ' 300  # tau-model, tauref(lifetime (fs)), tauexp(scattering parameter "r": 0 -> acoustic phonons, 2 -> ionized impurities), taurefen(Ry), taureftemp(K)\n'
    f_intrans += '0                         #number of fixed dopings\n'
    f_intrans += '1E20 -1E20          #fixed doping levels in cm3\n'
[PB1] http://www.quantum-espresso.org/wp-content/uploads/2013/06/tutorial_gwl.pdf (GWL)

□ VASP
util での vasp2boltz.txt を読んで vasp2boltz.py を動作させる。ASEの導入も必要。
■ boltztrap
-----
・python
1. yum install python (for CentOS)
1. sudo apt-get install python (for Ubuntu)
-----
・spglib-1.6.0 (for CentOS)
1.  download: http://spglib.sourceforge.net/
2. tar zxvf spglib*
3. cd spg*
4. ./confiugre
5. make
6. su
7. input password
8. make install
-----
・spglib-1.6.0 (for Ubuntu)
1.  download: http://spglib.sourceforge.net/
2. tar zxvf spglib*
3. cd spg*
4. ./configure
5. make
6. sudo make install
-----
・pyspglib-1.8.3.1 (for Ubuntu)
1.  download: https://pypi.python.org/pypi/pyspglib
2. tar zxvf pyspglib*
3. cd spg*
4. sudo python setup.py install
-----
・ASE (for Ubuntu)
1. sudo apt-get install python-ase
-----
・case.py
1. make new case.py file
------------------
from ase import io
from ase.lattice.spacegroup import Spacegroup
import vasp2boltz

ao = io.read('POSCAR')
# If you do not have spglib installed, you need to add the spacegroup no.
# e.g. sg = Spacegroup(216) for the 'F -4 3 m' space group
# If you use the conventional cell rather than the primitive unit cell in your
# POSCAR file for e.g. I or F spacegroups, you need to add this information.

# Example 1:
#sg = Spacegroup(1)
#ao.info = {'spacegroup': sg}

# Example 2. Here the unit cell in the POSCAR file is cubic (conventional):
#sg = Spacegroup(225)
#ao.info = {'spacegroup': sg, 'unit_cell': 'conventional'}

# Read number of electrons from OUTCAR.
foundnelect = False
try:
    for line in open('OUTCAR', 'r'):
        if 'NELECT' in line:
            nelect = float(line.split()[2])
        foundnelect = True
except:
    pass
if not foundnelect:
    print 'Number of electrons not found. Please set the number manually in hte.intrans.'
    nelect = 1.0

# The remaining part takes care of writing the files required by BoltzTraP.
bs = vasp2boltz.get_vasp_bandstructure()
vasp2boltz.write_bandstructure_boltztrap(bs)
vasp2boltz.write_structure_boltztrap(ao)
vasp2boltz.write_intrans_boltztrap(n_electrons = nelect)
------------------
2. put case.py and vasp2botz.py into files
3. python case.py
***************
None
get_kspace_operations(): atoms object has no space group information
-----
・boltztrap
1. mv energies.boltztrap case.energy
2. mv hte.intrans case.intrans
3. mv hte.struct case.struct
4. x_trans BoltzTraP
------
・plot
1. grep '  300.0000' '~/boltztrap-1.2.5/tests/case/case.trace' > caseout.txt
2. gnuplot
3. DOS vs. eV
    1) plot '~/boltztrap-1.2.5/tests/case/caseout.txt' u ($1*13.602):4 w l
    2) set zeroaxis
    3) replot
4. seebeck coefficient (microV/K vs. eV)
    1) plot '~/boltztrap-1.2.5/tests/case/caseout.txt' u 1:($5*1000000) w l
    2) set zeroaxis
    3) replot
5. conductivity/lifetime vs. eV
    1) plot '~/boltztrap-1.2.5/tests/case/caseout.txt' u ($1*13.602):6 w l
    2) set zeroaxis
    3) replot
I get same results for ZrCoSb comparing with http://www.phonon.t.u-tokyo.ac.jp/resource/yamamoto_b2012.pdf
動画:https://www.youtube.com/watch?v=OUauG1LM0QA
http://cms.mpi.univie.ac.at/vasp-forum/forum_viewtopic.php?4.9660 
http://www.phonon.t.u-tokyo.ac.jp/resource/yamamoto_b2012.pdf
https://sites.google.com/site/georgeyumnam/using-with-vasp
---------information----------
None is ktrafo. This case operate kp=enk['kpoint'].

Add "print bandstructure['E_Fermi']" at line 100 in vasp2boltz.py
vasp2boltz.py get fermi energy from OUTCAR,
vasp2boltz.py shift fermi energy and electron energy to zero.
I found that E_Fermi_zero is True at Line 141.
Line 142: e_up=e_up-bandstructure['E_Fermi']

pyspglib case must show 'get_kspace_operations(): atoms object has no space group information'
because vasp2boltz.py check atoms_info (spacegroup) at first.

energy unit is Ry.
    # output of bandstructure in a format which can be used by boltztrap code
    # input: bandstrcuture as dictionary like in get_vasp_bandstructure
    # default is unit conversion to Ry; if yscale is given this is taken as conversion factor
---------

■ spin を考慮した計算
conductivity or σ' = σ/ τで平均しなければならない。文献[4]には下記のようにある。
S = ( L11(↑) * S(↑) + L11(↓) * S(↓) ) / ( L11(↑) + L11(↓) )
where, S(↑) and S(↓) are calculated using the spin up and down band structure, respectively. In other words, the total
thermopower can be viewed as the average of spin up and down thermopowers weighted by the corresponding conductivity or σ' = σ/τ.
we note that L11 is, in fact, the conductivity .
L11 = L11(↑) + L11(↓)
where L↑ and L↓ are computed using the spin up and down band structures, respectively. Direct coupling of the thermopower to magnetic field via an energy dependent normal magnetoresistance is not included.

■ スピントロニクスに関する研究について
1)  他のオプションは、x_transを開けばわかる。スピントロニクスに関する研究を行うならば、-up や -dn のオプションを利用して輸送現象やゼーベック係数を検討してみるのも良いだろう。
2) 余裕があれば、他のcodeに対しても同様のプログラムを作成してみるとよい(例えば、TDDFTが利用可能なexcitingや開発中のPWscfなど)。
3) 電流 Je = J↑ + J↓ (J↑ :上向きスピンの電子の流れ。J↓:下向きのスピンのそれ)と、スピン流 Js = J↑ - J↓ を様々な物質で比較検討すればよい。
4) 「スピンホール効果:電流を流すと垂直方向にスピン流れが生じる」「逆スピンホール効果:スピン流を流すと垂直方向に電流が生じる」「スピンゼーベック効果:温度差によりスピン圧が生じる。スピン依存電気化学ポテンシャルの差『μ↑-μ↓』を計算すれば良い)。K. Uchida et al., Nature 455 (2008) 778.の論文を読んで検討してみるとよい」(磁性絶縁体でも検討してみられたい。加えて、日本物理学会誌,Vol 65, 2010. も参考にされるとよいだろう)

■ 生成するキャリア数について
生成するホールと電子のキャリア数は同数であるため下記のような関係になる。
hole carrier concentration =  electron carrier concentration

■ Pisarenko type plot
Seebeck coefficient vs electron count としたグラフも役に立つかもしれない。

■ 電気抵抗
ρ(T) /ρ(300 K)= a + b * T + exp(ΔE/T) [ER1]
ρ(T) = ρ0 + ρ1*ln(T) + ρ2 * (T/θD)n * Σ∫θD/T0  xn * e-kx dx [ER2]
  The first term ρ0 is the residual resistivity independent of temperature, the second term represents temperature dependent spin scattering, and the last term arises from the electron–phonon interaction. θD is the Debye temperature.
σ(T) = σ + σ1 * T1/2 + σ2 * exp(-EG/2kBT) [ER3]
ρ(T) = ρimp + a * T2 + b * T5 + c * T [ER4]
  Rimpは不純物によって決まる最低の電気抵抗で温度に依存しない。電子と電子の散乱がT2, フォノンによる電子の散乱がT5に比例する。
ρ(T) = ρ(0) + A * (T/θD)n * ∫θD/T0  xn / [(ex - 1)(1- e-x)] dx [ER5]
  ρ(0) is the residual resistivity due to defect scattering, A is a constant that depends on the velocity of electrons at the Fermi surface, the Debye radius and the number density of electrons in the metal. n is an integer that depends upon the nature of interaction:
n=5 implies that the resistance is due to scattering of electrons by phonons (as it is for simple metals)
n=3 implies that the resistance is due to s-d electron scattering (as is the case for transition metals)
n=2 implies that the resistance is due to electron–electron interaction.
[ER1] A. Slebarski et al., J. Phys.: condens. Matter 18 (2006) 10319.
[ER2] M. Vasundhara and V. Srinivas, Phys. Rev. B 77 (2008) 224415.
[ER3]
[ER4] 電気抵抗, wikipedia (Japanese)
[ER5] Electrical resistivity and conductivity, wikipedia (English); http://en.wikipedia.org/wiki/Electrical_resistivity_and_conductivity
[ER6] M. Vasundhara et al., J. Phys.: Condens. Matter 17 (2005) 6025.
[ER7] Heat transfer physics, wikpedia(English); http://en.wikipedia.org/wiki/Heat_transfer_physics
------------------------------------------------------------------------------

■ Theory
・τ-model[T4,T5]
the energy–power-dependent relaxation time: τ(E) = τ0*[(hk)2/2m*)]^(r-1/2) [T5]
k is the wave vector and m* the effective mass. r is the scattering parameter.
we estimated the Seebeck coefficients with various scattering parameters r as shown in Table I and found that r=2 was most suitable to explain the value of S at room temperature.
※ [T5]の論文では、scattering parameters r を変えたときのゼーベック係数の計算値 Srと実験値を比較し、r=2であると計算値が実験値に近い結果を与えることを報告している。最初の300 Kで比較し、後に300 K以下での温度変化でも上手く説明できることを示している。

the relaxation time: τ(E) = τ0(T)*[E/(kB*T)]^(r-1/2) [T4]
The values of r are not negative, and do not exceed the value 4 for known scattering processes.Since the exponent r takes on the values r=3/2 (the 2D case) and r=2 (the 3D case) for acoustic-phonon scattering, the asymptotic high-temperature dependence of S resulting from the Ziman variational formalism (25) corresponds to the linear dependence following from the relaxation-time approximation. The same concerns the 3D case considered in Refs. 2 and 3. [T4]
※ [T4]の論文では、弾性散乱のみを考慮した緩和時間(Eq. 31)と the parabolic (electron m>0) bandを考慮することで、上記の式を得ることができるとしている。

・ BoltzTraP
τ(E) = [τ0/T]*[E/(kB*T)]^(r-1/2) → τ(E) = [τ0*(T0/T)]*[(E-E0)/(kB*T)]^(r-1/2)
→τ(E) = [τ0*(T0/T)]*[(E-E0)/(kB*T)]*exp(r-1/2) ? → τ(E) = [τ0*(T0/T)]*exp{(r-0.5d0)*log([(E-E0)/(kB*T)]}
→τ(E) = τ0*exp{(r-0.5d0)*log([(E-E0)/(kB*T)]*[T0/T])}
※ 上記のような感じで式が作られていると思われるが、私には色々と分からない。

[T1] 竹内恒博, 日本金属学会誌第69巻第5 号(2005) 403-412 解説論文.; http://www.jim.or.jp/journal/j/pdf3/69/05/403.pdf 
[T2] Georg K. H. Madsen, David J. Singh, arXiv:cond-mat/0602203v1, 2008; http://arxiv.org/pdf/cond-mat/0602203v1.pdf 
[T3] Uichiro Mizutani, "Introduction to the Electron Theory of Metals", (Cambridge University Press, 2001)
[T4] K. Durczewski and M. Ausloos, Phys. Rev. B 61 (2000) 5303.
[T5] T. Okuda, et al., Phys. Rev. B 63 (2001) 113104.
[T6] 寺崎一郎, Netsu Sokutei 31(4)164-171; http://www.netsu.org/j+/Jour_J/pdf/31/31-4-164.pdf
[T7] http://kats.issp.u-tokyo.ac.jp/kats/semiconII/note4.pdf
[T8] http://www.phy.saitama-u.ac.jp/~saso/lectures/Lecture02.pdf

■ boltztrap
-------------------------------------------------------------------
□ 入出力ファイルの関係
計算時に入力するコマンド x_trans は *.def ファイルを作成する。
*.def ファイルはファイル読み込み番号が1になっている。
*.def はwriteやread文で読み込む番号との対応関係を記載している。

5: case.intrans
  入力ファイル。
・ setgap = 0 でない場合に shiftgap (gapchange) の値でギャップが形成される。(温度によるギャップの値の変化を考慮したプログラムになっているかは筆者自身が分かっていない)
・ HISTOの次の行には"tauref, tauexp, taurefen, tempexp, taureftemp"が読まれる。未入力の場合は tauref=1.0, tauexp=0.0, taurefen=0.0, taureftemp=-1.0となる。
  tauref: Reference lifetime (femtoseconds)
  tauexp: Scattering parameter "r" for exponential dependency of lifetime on energy
  Please, see PRB 63 113104 Okuda or PRB 61 5303 Durczewski for example
  taurefen: Reference energy point for energy dependent lifetime (Ry)
  taureftemp: Temperature for reference lifetime (K)
  taureffact = 1.d0
  もし、 taureftemp <= 1.0*10^-12ならば、taureffact = 1.d0 として計算される。
・ 上記のパラメーターはfermiintegrals.F90 で下記のように計算される。
  if (taureftemp > 1.d-12) taureffact = taureftemp/temp
  if (abs(tauexp) > 1.d-10) then
      ! the argument is evaluated in atomic units Ha/Ha
      !  temp is in Kelvin so temp*BOLTZMANN in Ha
      !  tauref is reference value for lifeteime (femtoseconds)
      !  taurefen is reference energy (rydberg), e.g. Fermi level, at which tauref is calculated
      !  taureftemp is reference temperature (Kelvin) at which tauref is given
      !  taureffact = taureftemp/temp is unitless (see just above)
      !  tauexp is the scattering parameter (unitless)
      !    0 -> acoustic phonons
      !    2 -> ionized impurities
      !  ene is in Hartree
      lifetime = (tauref/SECOND/1.e-15) * &
        &  exp((tauexp-0.5d0) * &
        &      log((abs(ene-taurefen*RYDBERG)/(temp*BOLTZMANN) &
        &           * taureffact)) )
  end if
  このようにして、温度とエネルギーに依存したlifetimeが計算される(何故logの部分をexpから出して前に置いて、論文の式の様にしていないのかは不明。数値計算での誤差を防ぐためであろうか?)。cond (the conductivity tensors: Eqs. 12 in [1])やnu (Eqs. 13 in [1]), kappa  (the electronic part of the thermal conductivity: Eqs. 14 in [1]), sigxyz  (Eqs. 15 in [1]) に liftimeの影響が考慮される。当然ながらゼーベック係数にも影響を与える。
・ 更に次の行では 0 よりも大きい数値が記述されている場合に、その次の行でドープ 量が読み取られる。単位は下記となる。
  doping: Doping levels to be output for, in carriers / cm^3
-------------------------------------------------------------------
case.intrans (for Al)  for BoltzTraP-1.2.2
WIEN                      # Format of DOS                                         
0 0 0 0.0                 # iskip (not presently used) idebug setgap shiftgap                         
0.49695 0.0005 0.4   9.   # Fermilevel (Ry), energygrid, energy span around Fermilevel, number of electrons             
CALC                      # CALC (calculate expansion coeff), NOCALC read from file                     
3                         # lpfac, number of latt-points per k-point                                    
BOLTZ                     # run mode (only BOLTZ is supported)                                 
.15                       # (efcut) energy range of chemical potential                             
300. 10.                  # Tmax, temperature grid                                     
-1                        # energyrange of bands given individual DOS output sig_xxx and dos_xxx (xxx is band number)
HISTO                  #scheme to obtain DOS. HISTO/TETRA: histogram/thetrahedron[4] sampling
32.0 2.0 0.49695 300      #τ-model. tauref(Reference lifetime (femtoseconds)), tauexp(scattering parameter "r": 0 -> acoustic phonons, 2 -> ionized impurities), taurefen(Ry), taureftemp(K)
0                         #number of fixed dopings
1E20 -1E20          #fixed doping levels in cm3
-------------------------------------------------------------------

6: case.outputtrans
  出力ファイル。case.intransからの計算条件なども表示される。

20: case.struct
  WIEN2kからの出力ファイル

10: case.energy + so or updn
  WIEN2kからの出力ファイル。

48: case.engre

49: case.transdos
  dosの結果が格納される。iが0からnpointsまでデータが出力される(elの低いエネルギーから高いエネルギーまでに対応)。
  一行目が el でエネルギー(Ry)。二行目が dos1(i) で電子状態密度分布。三行目がdosint(i) で電子数の積算。

50: case.sigxx
  e1,dos_sigxy(1:3,1:3,i) を出力する。

51: case.sigxxx
  e1,dos_sigxyz(1:3,1:3,1:3,i) を出力する。

21: case..trace
  # Ef[Ry] T [K] N DOS(Ef) S s/t R_H kappa0 c chi の結果を出力する。

22: case.condtens

24: case.halltens

※ 11 はdefの定義には無いが、存在すると taumodel が .TRUE. になる。サブルーチン generic_tau (genetic_lifetime.F90)とfite4がBoltzTraP.F90で呼ばれる。
file 11 の構成
1行目: タイトル
2行目: 読み取るk点数('Number of lifetime kpts in IBZ: ',nkpt_lt)
3行目: 1回目は S, T, Z, neup (動的メモリ配置の大きさを決めるために、ファイルにある最大値が調べられてnumeに格納される)
  2回目では k(1,k), k(2,k), kz(3,k), neup で読み込まれる。
4行目: neup だけ"life time"(寿命)の値を読み込む(1回目 e1で、2回目にlifetime(ii, k)で読み込まれる。ここでは、iiがneupまで読み込まれる)。
※ 1回目の読み込みが終了すると、サブルーチン init_ltmod が呼ばれ、動的メモリが宣言されて、ZERO(0.0d0)に初期化する。2回目に正式にlifetimeのデータが読み込まれる。
※ File number 11 with k and e dependent lifetimes found であるため、k点とエネルギーでの lifetimeを読み込んだプログラムであると考えられる。
-------------------------------------------------------------------
□ BoltzTraP.F90
メインルーチン。プログラムの全体の流れを理解することにも役立つ。以下が大まかな流れ。

USE reallocate: reallocate.F90 において module reallocate と宣言されている。
  reallocate.F90 で定義されたデータを用いる。
USE defs: gmlib2.F90 において MODULE defs が宣言されている。
基本的な定数について引数の定義がなされている。
USE input: m_input.F90 において MODULE input と宣言されている。
  intrans に記述されているデータの宣言に対応する。
  #of electrons は nval に対応している。
USE bandstructure: m_bandstructure.F90 で MODULE bandstructure が宣言されている。
USE ltmod: generic_lifetime.F90 で MODULE ltmod が宣言されている。
USE lattice_points: gen_lattpoints.F90 で MODULE lattice_points が宣言されている。

CALL gtfnam: extract the command-line argument
  入力ファイル名のチェックをしているように見える。

BoltzTraP.def ファイルを開いて、中に書かれている番号とファイル名などを対応付ける。
上から順番にファイルを開いていき、入力ファイルが全て揃っているかをチェック。
エラーがあれば表示される。
入出力でのファイル番号はディフォルトで6がディスプレイであるが、
BoltzTraP.defから読み込むファイル(case.outputtrans)にwrite(6,*)での記述が書き込まれることになる。

入力ファイル名で energyso がある場合は、spinorbit=ONEに設定する。(通常、spinorbit=TWO)

CALL read_input: m_input.F90 で subroutine read_input が宣言されている。

読み込むファイル形式を選択して、データを取得している。

CALL add_inv: 
! Add center off symmetry for non-centrosymmetric lattices

case.intrans で Run type: に CALCと記入した場合、modus1 は CALC となっており、egap や nwave が計算される。

CALL bandana: case.intrans で Run type: に CALCと記入した場合に呼ばれる。
バンドギャップを広げる指定をした場合は、INT(nval/spinorbit)+1以上の番号のバンド(case.energyに書かれている)のエネルギー(Ry) を shiftgap (プログラム中はgapchange) だけ加算する。
icut1には計算する範囲で最も低いエネルギーのバンドの番号が格納される。
icut2には計算する範囲で最も高いエネルギーのバンドの番号が格納される。

CALL fite4: case.intrans で Run type: に CALCと記入した場合に呼ばれる。
Interpolation scheme (PRB 38 p.2721)
The calculation of the expansion parameters, cR,i, are carried out in the subroutine FITE4.

boltztrap-1.2.2ではtaumodel = .FALSE.なので CALL generic_tau などは計算されない。
※ (File number 11 with k and e dependent lifetimes found)とした計算も可能なように改良してみるのも良いだろう。
もし、taumodelが.TRUE.であるならば、fite4で lifetimeを考慮した計算がなされて結果がtaugreに格納される。

run modeはBOLTZのみだけサポートされているので、modus2はBOLTZとなる。
ecut2(Range around Ef where bands are given individual output (Ry) ) > ZERO(0.0d0)の場合に CALL bandana 呼ばれる。

CALL dos: 電子状態密度分布(DOS)が計算されて、No.49-51に出力される。

CALL fermiintegrals:  
case.trace にある # Ef[Ry] T [K] N DOS(Ef) S s/t R_H kappa0 c chi という文字を出力する。
各温度に対するその数値はサブルーチンであるfermiint_fix_ef_Tで出力される。

'TRANSPORT END BoltzTrap calculation' と記述。

END PROGRAM BoltzTrap
-------------------------------------------------------------------
■ 引数の意味
□ ecut: energy span around Fermilevel
  計算に考慮に入れるバンド(のエネルギー)範囲。論文[1]での式12-15におけるエネルギーの積分範囲。
  icut1とicut2が決定される。icut1とicut2はcase.outputtransで Bands range: ',icut1 ,' - ',icut2 として表示される。',icut1 ,' - ',icut2 に対応するエネルギーはその左のEnergy rangeに書かれている。
  ※ バンド1本が1電子に対応する。あるエネルギー範囲にあるバンドの数が、そのエネルギー範囲にある電子数になる。
□ efcut: (efcut) energy range of chemical potential
  化学ポテンシャル μ を変えて計算する範囲
  range of μ in which the integrations should be performed
□ nval: number of  electrons, nval = NE in case.in2_st1
□ bandana.F90: egap=emin(INT(nval/spinorbit)+1)-emax(INT(nval/spinorbit))
 バンドギャップを広げる場合(setgapが0以外)にnvalが用いられる。
□ fermiintegrals.F90: sumelec=nval-(icut1-1)*spinorbit-(sumelec*deltaef*volume)
□ deltaef: energy step (step size), npoints = (ebmax-ebmin)/deltaef
  整数となる ceiling 関数 を用いてnpoints = ceiling((ebmax-ebmin)/deltaef となるかチェックが入る。
□ sumelec: 正孔のキャリア数=電子数(nvalの値) - (icut1までの電子数-1) *倍率(energysoを用いた場合1、それ以外2) - (フェルミディラック分布を掛けた電子密度分布をicut1からicut2まで合計*エネルギー刻み*体積)
  ※ 電子のキャリア数 = ∫D(ε) * FD(ε) dε から計算できる[6]。ここで、Dは電子状態密度分布、FDはフェルミディラック関数。εはエネルギー。
□ phon_band.F90: nval=nband/TWO
※ spinorbit は energyso を用いた計算の場合(ONE=1.0d0)、それ以外の場合はTWO=2.0d0 となっている。
------------------------------------------------------------------------------
□ 開発者からのコメント(著者の2人のうち、ボスの方からのみ返信が頂けた)
・ The chemical potential does shift with temperature but is not directly printed. The program instead prints things as a function of T and energy. However from the variation of electron number with energy you can get the chemical potential i.e. as the energy that gives the desired electron count.
Usually, I use a plotting script to get what I need. For example using the electron count rather than the energy to construct the x axis for a Pisarenko type plot. The program includes source, so you can change the formats to allow the desired number of digits.
以上の記述から、温度でフェルミ準位が大きく移動するような系の場合は、化学ポテンシャルと温度の組を選び出さなくてはならない。さらに、文献[5]では温度によるバンドギャップの変化を考慮することで、緩和時間近似を用いていても、Be2Te3において実験値と比較的というよりもかなり良い一致を示すことが分かる。
------------------------------------------------------------------------------
BoltzTraP を改良すべき点を挙げる
1. フォノンドラッグの項を入れる[C1]
2. 磁気での補正を入れる[C2-C4]
3. (スピン)エントロピーを考慮に入れる[C5] (NaxCoO2は既存のBoltzTraPでかなり一致するので必要ないかもしれない)
4. 温度による格子定数の変化を考慮に入れる
  格子定数はDebye lattice vibration model で計算できる。aL=KL * γL * CL / (3*V)
ここで CLは比熱、KLは等温圧縮率、γLはGruneissen定数。KLとγLは温度に依存しないと仮定。KL * γL / (3*V) = aexp / CL = A を常温で計算する。常温以下の温度ではaL=A*CL(T)として計算する。CL(T)はリートベルトから得られたデバイ温度を用いてDebyeモデルで計算する。 デバイ温度は「Using Debye-Waller factors Bi ,i5Fe, Ti, or Sn, obtained from the fit, the Debye temperatures θDi at T=300 K were calculated from the relation Bi(T) = (6h2T/kB * mi * θDi2 ) * [φ(θDi/T)+θDi/T], φ(θDi/T) =φ(x)=(1/x)∫0 x (ydy/ey-1), where mi is the atomic mass. The value of the Debye temperature」から得られる。 [C6]
[C1] http://iroha.scitech.lib.keio.ac.jp:8080/sigma/bitstream/handle/10721/667/document.pdf?sequence=1
[C2] G. Sundaram et al., Phys. Rev. B 59 (1999) 14915.; http://prb.aps.org/pdf/PRB/v59/i23/p14915_1
[C3] Di Xiao et al., Phys. Rev. Lett. 97 (2006) 026603.; http://prl.aps.org/pdf/PRL/v97/i2/e026603
[C4] Y. Hasegawa et al., Jpn. J. Appl. Phys. 43 (2004) pp. 35-42.; http://jjap.jsap.jp/link?JJAP/43/35/
[C5] 寺崎一郎、Netsu Sokutei 31 (4) 164-171.; http://www.netsu.org/j+/Jour_J/pdf/31/31-4-164.pdf
[C6] A. Slebarshki et al., Phys. Rev. B 62 (2000) 3296.; http://prb.aps.org/abstract/PRB/v62/i5/p3296_1
------------------------------------------------------------------------------
■ case.energy
  At the top are the energy-parameters E_l for each atom for the APW-expansion of u_l(r,E_l) (first line) and for the LOs (second line). Then comes the first k-point (0,0,0). It had a basis size of 219 and 44 eigenvalues (in Ry) follow. Then the next k-point, .... (P. Blaha)
200.60479200.60479200.67530  0.60479  0.60479  0.60479  0.60479  0.60479 0.60479  0.60479  0.60479  0.60479  0.60479  0.00000
   0.60479  0.60479  0.67530999.00000997.00000 -3.44000997.00000999.00000999.00000999.00000999.00000999.00000
  0.000000000000E+00 0.000000000000E+00 0.000000000000E+00         1   219  44  1.0
            1  -3.42926916127842
            2  -3.42926916127841
           ・・・・・・
          44   1.93420513481993
  next k-point kx ky kz (basis size) (# of eigenvalues) (weight for this k point)
        (# of eigenvalue) (eigenvalue)
------------------------------------------------------------------------------
理論値よりも実験値の方がゼーベック係数が大きい場合
1. 低温においては、フォノンの影響によるフォノンドラッグによってゼーベック係数が大きい結果が得られる。
2. エントロピー(磁性)の影響によりゼーベック係数が大きい結果が得られる。
3. 温度上昇に伴って格子定数が増大する場合。
[ES1] 竹内恒博, 日本金属学会誌第69巻第5 号(2005) 403-412 解説論文.; http://www.jim.or.jp/journal/j/pdf3/69/05/403.pdf 
[ES2] Y. Wang et al., Nature 423 (2003) 425-428.; http://www.nature.com/nature/journal/v423/n6938/abs/nature01639.html
[ES3] I. Terasaki et al., Phys. Rev. B 65 (2002) 195106.; http://prb.aps.org/abstract/PRB/v65/i19/e195106
[ES4] 寺崎一郎, Netsu Sokutei 31(4)164-171; http://www.netsu.org/j+/Jour_J/pdf/31/31-4-164.pdf
------------------------------------------------------------------------------
□ References
[1] Georg K. H. Madsen, David J. Singh, arXiv:cond-mat/0602203v1, 2008; http://arxiv.org/pdf/cond-mat/0602203v1.pdf
[2] Lijun Zhang and D. J. Singh, Phys. Rev. B 80 (2009) 075117.
[3] D. J. Singh and I. I. Mazin, Phys.  Rev. B 56 (1997) R1650.
[4] H. J. Xiang and D. J. Singh et al., Phys. Rev. B 76 (2007) 195111.
[5] http://books.google.co.jp/books?id=AoRHN_Ah5WkC&pg=PA125&lpg=PA125&dq=boltztrap&source=bl&ots=uxqtKrEcmt&sig=NkK-mmWCuQORnm9MrKuvKvRk2hM&hl=ja&sa=X&ei=xsHvUYLDLsytkgWDu4CwCQ&ved=0CFQQ6AEwBTgK#v=onepage&q=boltztrap&f=false
[6] D. F. Zou1, S. H. Xie1, Y. Y. Liu1, J. G. Lin1, and J. Y. Li2 , J. Appl. Phys. 113, 193705 (2013); http://dx.doi.org/10.1063/1.4804939 (7 pages); "Electronic structure and thermoelectric properties of half-Heusler Zr05Hf05NiSn by first-principles calculations "
[7] http://www.research.kobe-u.ac.jp/eng-nanoelectronics/Japanese/tsuchiya/suributuri/suributuri.pdf
------------------------------------------------------------------------------
■ 開発中(温度によるフェルミ準位の移動を考慮してゼーベック係数のデータを抜き出すIgorマクロ)
金属や擬ギャップ系のみ(半導体や絶縁体はFD分布の数値が広いエネルギー範囲で小さいので計算方法を変える予定)
□ 使い方
1. case.output2またはdos1, dos1upファイルの中に書かれているFermi準位をコピーして、マクロのEfにペーストしてEnterを押す。
2. BoltzTraP から case.trace を取り出して、Igor に入れる。
3. change name を押すと名称を wave0 などからマクロで処理する名称に書き換わる。
4. Temperature [K] を指定すれば、S-μ relation ボタンで図が表示される。series は全ての温度に対して描画される。
□ Igor マクロ

#pragma rtGlobals=1  // Use modern global access method.
// ver.1.00

Macro seebeck_panel()
NewPanel/W=(0, 0, 180, 190)
SetDrawEnv fillfgc= (48896,65280,48896)
DrawRect 4,2,176,188

Button button0 proc=ButtonProc_1
Button button0 title="change name",proc=ButtonProc_1
SetDrawEnv fsize= 14
Button button0 size={150,20}
Button button0 pos={15,3}

SetVariable setvar0 proc=SetVarProc,limits={-inf,inf,0},fSize=14
SetVariable setvar0 value= _NUM:0
SetVariable setvar0 size={70,20}
SetVariable setvar0 pos={55,30}

DrawText 25,45,"Ef ="
SetDrawEnv fsize= 14

DrawText 130,45,"Ry"
SetDrawEnv fsize= 14

PopupMenu popup0 proc=PopMenuProc
PopupMenu popup0 value="50;100;150;200;250;300;350;400;450;500;550;600;650;700;750;800;850;900;950;1000"
PopupMenu popup0 fSize=14
PopupMenu popup0 pos={120,55}

DrawText 15,70,"Temperature [K]"
SetDrawEnv fsize= 14

Button button1 title="S - μ relation",proc=ButtonProc
SetDrawEnv fsize= 14
Button button1 size={100,30}
Button button1 pos={60,80}

Button button2 title="series",proc=ButtonProc_2
SetDrawEnv fsize= 14
Button button2 size={150,20}
Button button2 pos={15,115}

SetDrawEnv fsize= 14
DrawText 25,155,"μ ="

SetVariable setvar1 proc=SetVarProc_1,limits={-inf,inf,0}
SetVariable setvar1 value= _NUM:0
SetVariable setvar1 size={50,15}
SetVariable setvar1 pos={60,140}

Button button3 title="plot",proc=ButtonProc_3
Button button3 size={50,20}
Button button3 pos={115,138}

Button button4 proc=ButtonProc_4
Button button4 title="series plot"
Button button4 size={70,20}
Button button4 pos={55,165}
Button button4 fSize=14

End

// convert Ry to eV and change name
Function ButtonProc_1(ba) : ButtonControl
 STRUCT WMButtonAction &ba

 switch( ba.eventCode )
  case 2: // mouse up
   // click code here
   WAVE wave0, wave1, wave2, wave3, wave4, wave5, wave6, wave7, wave8, wave9
   // rename
   Rename wave0,Ef_Ry; Rename wave1,T_K; Rename wave2,N; Rename wave3,DOS_Ef
   Rename wave4,Seebeck; Rename wave5,S_t; Rename wave6,R_H; Rename wave7,kappa0
   Rename wave8,c; Rename wave9,chi
   break
  case -1: // control being killed
   break
 endswitch

 return 0
End

// read non_dope_EF
Function SetVarProc(sva) : SetVariableControl
 STRUCT WMSetVariableAction &sva

 switch( sva.eventCode )
  case 1: // mouse up
  case 2: // Enter key
  case 3: // Live update
   Variable dval = sva.dval
   String sval = sva.sval
   Variable/G Ef_non_dope
   Variable Num_of_Wave
   WAVE Ef_Ry
   Num_of_Wave = DimSize(Ef_Ry,0)
   Make/O/N=(Num_of_Wave) Ef_eV
   Ef_non_dope = dval
   Ef_eV = (Ef_Ry - Ef_non_dope)*13.602
   break
  case -1: // control being killed
   break
 endswitch

 return 0
End

// set extacting temperature
Function PopMenuProc(pa) : PopupMenuControl
 STRUCT WMPopupAction &pa

 switch( pa.eventCode )
  case 2: // mouse up
   Variable popNum = pa.popNum
   String popStr = pa.popStr
   Variable/G temp
   temp = str2num (popStr)
   break
  case -1: // control being killed
   break
 endswitch

 return 0
End


// run bottun
Function ButtonProc(ba) : ButtonControl
 STRUCT WMButtonAction &ba

 switch( ba.eventCode )
  case 2: // mouse up
   // click code here
   NVAR temp
   print  "extract the", temp,"K data from the case.trace"
   sub_extract(temp)
   ModifyGraph width=0,height=0
   break
  case -1: // control being killed
   break
 endswitch

 return 0
End


// main routine
Function sub_extract(temp_K)
  Variable temp_K
  WAVE Ef_eV, T_K, N, DOS_Ef, Seebeck, S_t, R_H, kappa0, c, chi
 
  Variable Num_of_Wave, i, j
  Num_of_Wave = DimSize(T_K,0)
  Make/N=(Num_of_Wave) Shifted_Ef_tempo, Seebeck_at_Shifted_Ef_tempo
 
  i = 0; j = 0
  do
    if ( T_K[i] == temp_K )
      Shifted_Ef_tempo[j] = Ef_eV[i]
      Seebeck_at_Shifted_Ef_tempo[j] = Seebeck[i]
      j = j + 1
    endif
    i = i +1
  while ( i < DimSize(T_K,0) )
 
  String name_of_wave
  name_of_wave =  "Seebeck_at_" + num2str(temp_K) + "_K"
  Make/O/N=(j) Shifted_Ef, $name_of_wave
  Wave Seebeck_at_Shifted_Ef = $name_of_wave

  Shifted_Ef = Shifted_Ef_tempo
  Seebeck_at_Shifted_Ef = Seebeck_at_Shifted_Ef_tempo * 10^6
 
  Display  Seebeck_at_Shifted_Ef vs Shifted_Ef
 
  //setting graph
  ModifyGraph tick=2,mirror=1,fSize=14,standoff=0,font="Arial"
  Label left "\\Z14\rSeebeck coefficients / \\F'Symbol'm\\F'Arial'V\\S-1"
  ModifyGraph tick(bottom)=2,mirror(bottom)=1,fSize(bottom)=14,standoff(bottom)=0
  ModifyGraph font(bottom)="Arial"
  Label bottom "\\Z14\\F'Symbol'm\\F'Arial' / eV"
  SetAxis bottom -2,2
  ModifyGraph margin(left)=57
  ModifyGraph margin(bottom)=43
  ModifyGraph width=340.157,height=226.772
  Legend/C/N=text0/F=0/B=1/A=MC/X=30.00/Y=40.00
 
  KillWaves Shifted_Ef_tempo, Seebeck_at_Shifted_Ef_tempo
 
End

// series
Function ButtonProc_2(ba) : ButtonControl
 STRUCT WMButtonAction &ba

 switch( ba.eventCode )
  case 2: // mouse up
   // click code here
     Variable temp_K
     temp_K = 100
     sub_extract(temp_K)
     for ( temp_K=200; temp_K<=800; temp_K+=100)
       sub2_extract(temp_K)
     endfor
     ModifyGraph rgb(Seebeck_at_200_K)=(65280,21760,0)
     ModifyGraph rgb(Seebeck_at_400_K)=(40960,65280,16384)
     ModifyGraph rgb(Seebeck_at_500_K)=(0,52224,0)
     ModifyGraph rgb(Seebeck_at_600_K)=(0,65280,65280)
     ModifyGraph rgb(Seebeck_at_700_K)=(0,0,65280)
     ModifyGraph rgb(Seebeck_at_800_K)=(29440,0,58880)
     Legend/C/N=text0/F=0/B=1/A=MC/X=30.00/Y=25.00
     ModifyGraph width=0,height=0
   break
  case -1: // control being killed
   break
 endswitch

 return 0
End

//second main
Function sub2_extract(temp_K)
  Variable temp_K
  WAVE Ef_eV, T_K, N, DOS_Ef, Seebeck, S_t, R_H, kappa0, c, chi
 
  Variable Num_of_Wave, i, j
  Num_of_Wave = DimSize(T_K,0)
  Make/N=(Num_of_Wave) Shifted_Ef_tempo, Seebeck_at_Shifted_Ef_tempo
 
  i = 0; j = 0
  do
    if ( T_K[i] == temp_K )
      Shifted_Ef_tempo[j] = Ef_eV[i]
      Seebeck_at_Shifted_Ef_tempo[j] = Seebeck[i]
      j = j + 1
    endif
    i = i +1
  while ( i < DimSize(T_K,0) )
 
  String name_of_wave
  name_of_wave =  "Seebeck_at_" + num2str(temp_K) + "_K"
  Make/O/N=(j) Shifted_Ef, $name_of_wave
  Wave Seebeck_at_Shifted_Ef = $name_of_wave

  Shifted_Ef = Shifted_Ef_tempo
  Seebeck_at_Shifted_Ef = Seebeck_at_Shifted_Ef_tempo * 10^6
 
  AppendToGraph  Seebeck_at_Shifted_Ef vs Shifted_Ef
 
  KillWaves Shifted_Ef_tempo, Seebeck_at_Shifted_Ef_tempo
End

//
Function SetVarProc_1(sva) : SetVariableControl
 STRUCT WMSetVariableAction &sva

 switch( sva.eventCode )
  case 1: // mouse up
  case 2: // Enter key
  case 3: // Live update
   Variable dval = sva.dval
   String sval = sva.sval
   Variable/G mu =  dval
   break
  case -1: // control being killed
   break
 endswitch

 return 0
End

// S - K plot
Function ButtonProc_3(ba) : ButtonControl
 STRUCT WMButtonAction &ba

 switch( ba.eventCode )
  case 2: // mouse up
   // click code here
   NVAR mu
   sub3_extract(mu)
   ModifyGraph width=0,height=0
   break
  case -1: // control being killed
   break
 endswitch

 return 0
End

//third main
Function sub3_extract(mu)
  Variable mu
  WAVE Ef_eV, T_K, N, DOS_Ef, Seebeck, S_t, R_H, kappa0, c, chi
 
  Variable Num_of_Wave, i, j, k, num_of_temp
  Variable mu_low, mu_high
  Num_of_Wave = DimSize(T_K,0)
  Make/N=(Num_of_Wave) K_tempo, Seebeck_for_K_tempo
 
  i = 0; j = 0; k = 0
  do
    if ( Ef_eV[i] <= mu && mu <= Ef_eV[i+1] )
      mu_low = Ef_eV[i]
      mu_high = Ef_eV[i+1]
      do
        k = k + 1
      while ( T_K[i-k] < T_K[i-k+1] )
      num_of_temp = k
      for(k=0;k<num_of_temp;k+=1)
        K_tempo[k] = T_K[i+k+1]
        Seebeck_for_K_tempo[k] = (Seebeck[i+k+1] - Seebeck[i-num_of_temp+k+1]) * (mu - mu_low) / ( mu_high - mu_low ) +  Seebeck[i-num_of_temp+k+1]
        //print  T_K[i+k+1], T_K[i-num_of_temp+k+1]
        //print  Seebeck[i+k+1], Seebeck[i-num_of_temp+k+1]
      endfor
      break
    endif
    i = i +1
  while ( i < (DimSize(T_K,0) -1) )
 
  String name_of_wave
  name_of_wave = "Seebeck_at_" + num2str(mu) + "_eV"
  Make/O/N=(num_of_temp) temp_K, $name_of_wave
  Wave Seebeck_for_K = $name_of_wave

  temp_K = K_tempo
  Seebeck_for_K = Seebeck_for_K_tempo * 10^6
 
  Display  Seebeck_for_K vs temp_K
 
  //setting graph
  ModifyGraph mode=4,marker=19,msize=4
  ModifyGraph tick=2,mirror=1,fSize=14,standoff=0,font="Arial"
  Label left "\\Z14\rSeebeck coefficients / \\F'Symbol'm\\F'Arial'V\\S-1"
  ModifyGraph tick(bottom)=2,mirror(bottom)=1,fSize(bottom)=14,standoff(bottom)=0
  ModifyGraph font(bottom)="Arial"
  Label bottom "\\Z14\\F'Arial'Temperature / K"
  SetAxis bottom 50,800
  ModifyGraph margin(left)=57
  ModifyGraph margin(bottom)=43
  ModifyGraph width=340.157,height=226.772
  Legend/C/N=text0/F=0/B=1/A=MC/X=30.00/Y=40.00
 
  KillWaves K_tempo, Seebeck_for_K_tempo
 
End

Function ButtonProc_4(ba) : ButtonControl
 STRUCT WMButtonAction &ba

 switch( ba.eventCode )
  case 2: // mouse up
   // click code here
     NVAR mu
     Variable mu_shift
     Display
     for ( mu_shift=-0.15; mu_shift<=0.151; mu_shift+=mu)
       sub4_extract(mu_shift)
     endfor
     Legend/C/N=text0/F=0/B=1/A=MC/X=30.00/Y=25.00
     ModifyGraph width=0,height=0
   break
  case -1: // control being killed
   break
 endswitch

 return 0
End

//fourth main
Function sub4_extract(mu)
  Variable mu
  WAVE Ef_eV, T_K, N, DOS_Ef, Seebeck, S_t, R_H, kappa0, c, chi
 
  Variable Num_of_Wave, i, j, k, num_of_temp
  Variable mu_low, mu_high
  Num_of_Wave = DimSize(T_K,0)
  Make/N=(Num_of_Wave) K_tempo, Seebeck_for_K_tempo
 
  i = 0; j = 0; k = 0
  do
    if ( Ef_eV[i] <= mu && mu <= Ef_eV[i+1] )
      mu_low = Ef_eV[i]
      mu_high = Ef_eV[i+1]
      do
        k = k + 1
      while ( T_K[i-k] < T_K[i-k+1] )
      num_of_temp = k
      for(k=0;k<num_of_temp;k+=1)
        K_tempo[k] = T_K[i+k+1]
        Seebeck_for_K_tempo[k] = (Seebeck[i+k+1] - Seebeck[i-num_of_temp+k+1]) * (mu - mu_low) / ( mu_high - mu_low ) +  Seebeck[i-num_of_temp+k+1]
        //print  T_K[i+k+1], T_K[i-num_of_temp+k+1]
        //print  Seebeck[i+k+1], Seebeck[i-num_of_temp+k+1]
      endfor
      break
    endif
    i = i +1
  while ( i < (DimSize(T_K,0) -1) )
 
  String name_of_wave
  name_of_wave = "Seebeck_at_" + num2str(mu) + "_eV"
  Make/O/N=(num_of_temp) temp_K, $name_of_wave
  Wave Seebeck_for_K = $name_of_wave

  temp_K = K_tempo
  Seebeck_for_K = Seebeck_for_K_tempo * 10^6
 
  AppendToGraph  Seebeck_for_K vs temp_K
 
  //setting graph
  ModifyGraph mode=4,marker=19,msize=4
  ModifyGraph tick=2,mirror=1,fSize=14,standoff=0,font="Arial"
  Label left "\\Z14\rSeebeck coefficients / \\F'Symbol'm\\F'Arial'V\\S-1"
  ModifyGraph tick(bottom)=2,mirror(bottom)=1,fSize(bottom)=14,standoff(bottom)=0
  ModifyGraph font(bottom)="Arial"
  Label bottom "\\Z14\\F'Arial'Temperature / K"
  SetAxis bottom 0,2000
  ModifyGraph margin(left)=57
  ModifyGraph margin(bottom)=43
  ModifyGraph width=340.157,height=226.772
  Legend/C/N=text0/F=0/B=1/A=MC/X=30.00/Y=40.00
 
  KillWaves K_tempo, Seebeck_for_K_tempo
 
End


 ------------------------------------------------------------------------------
QRコード
携帯用QRコード
アクセス数
ページビュー数
[無料でホームページを作成] [通報・削除依頼]