Kelvinの公式でSeebeck係数

AkaiKKR(machikaneyama)で計算された状態密度から Kelvin の公式を利用して、遷移金属のゼーベック係数を計算しました。
\begin{equation}
S = - \frac{1}{|e|}\frac{\mathrm{d}\mu}{\mathrm{d}T}
\end{equation}
ここで S はゼーベック係数、 e は電気素量、 μ は化学ポテンシャル、T は絶対温度です。結果はそこそこ良く実験値を再現しました。

Kelvin.png
Fig.1: Kelvinの公式で計算されたゼーベック係数(実線)と実験から得られた文献値(丸シンボル)。パラジウム(赤)、プラチナ(青)、タングステン(緑)、モリブデン(黒)。



ゼーベック係数


ゼーベック係数(熱電能)は、以下の式で表されます。

\begin{equation}
S = \frac{1}{eT}\frac{K_1}{K_0}
\end{equation}
ここでKn
\begin{equation}
K_n = \int_{-\infty}^{\infty}\sigma(\epsilon)(\epsilon - \mu)^{n} \left( - \frac{\mathrm{d}f}{\mathrm{d}\epsilon}\right) \mathrm{d}\epsilon
\end{equation}

この式の中の σ(ε) はエネルギーに依存する電気伝導度とでも呼ぶべきもので、これを具体的に計算するには、通常のバンド計算から得られる電子の群速度に加えて、電子がどのように散乱されるかを表す散乱時間も必要になります。これは不可能ではありませんが、結構大変です。

これに対して熱電材料の物質科学―熱力学・物性物理学・ナノ科学 (物質・材料テキストシリーズ)では、とても簡単な近似式として Kelvin の公式を示しています。
\begin{equation}
S = - \frac{1}{e}\frac{\mathrm{d}\mu}{\mathrm{d}T}
\end{equation}
Scilabで金属の化学ポテンシャルに書いたとおり、金属の状態密度さえ分かっていれば化学ポテンシャルは計算できるので、バンド計算的には、とても簡単な手法です。(というか、これで精度よくゼーベック係数が計算できるのなら、群速度とか散乱時間とか一体なんだったのという感じ)

計算手順


まずAkaiKKRで状態密度を計算しました。2500 K 程度なら化学ポテンシャルの大きさもたいしたこと無いはずなので ewidth を小さくして計算すべきですが、雑な計算ということで価電子を全て含むエネルギーとしました。

c------------------------------------------------------------
go data/Pt
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
fcc 7.41 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 0.9 sra mjw nmag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 4 200 0.023
c------------------------------------------------------------
c ntyp
1
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Pt 1 1 0.0 2
78 100
c------------------------------------------------------------
c natm
1
c------------------------------------------------------------
c atmicx atmtyp
0 0 0 Pt
c------------------------------------------------------------

c------------------------------------------------------------
dos data/Pt
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
fcc 7.41 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 1.2 sra mjw nmag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 20 200 0.023
c------------------------------------------------------------
c ntyp
1
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Pt 1 1 0.0 2
78 100
c------------------------------------------------------------
c natm
1
c------------------------------------------------------------
c atmicx atmtyp
0 0 0 Pt
c------------------------------------------------------------


Pt-DOS.png
Fig.2: プラチナの状態密度


次に計算された状態密度から、電子の数密度neを計算します。
\begin{equation}
n_e = \int_{-\infty}^{\infty}D(\epsilon)f(\epsilon, T)\mathrm{d}\epsilon
\end{equation}
ここでフェルミ分布関数は以下のようになります。
\begin{equation}
f(\epsilon, T) = \frac{1}{\exp \left(\frac{\epsilon - \mu(T)}{k_B T} \right) + 1}
\end{equation}

電子の数密度neは温度に関わらず一定なので、絶対零度 T = 0 (K) のときの電子数密度 ne0 をまず計算します。フェルミ分布関数はこのとき ε < εF で f(ε, 0) = 1, ε > εF で f(ε, 0) = 0 なので
\begin{equation}
n_{e0} = \int_{-\infty}^{0} D(\epsilon) \mathrm{d}\epsilon
\end{equation}
です。

あとは以下の条件を満たすように非線型方程式ソルバで、化学ポテンシャルμを求めます。
\begin{equation}
\int_{-\infty}^{\infty}D(\epsilon)f(\epsilon, T)\mathrm{d}\epsilon - n_{e0} = 0
\end{equation}

化学ポテンシャルが計算できたら、これを数値微分します。
\begin{equation}
\frac{\mathrm{d}\mu}{\mathrm{d}T} \sim \frac{\mu(T+\Delta T) - \mu(T - \Delta T)}{2\Delta T}
\end{equation}

clear;

// *** 物理定数 ***
// アボガドロ数 (/mol)
na = 6.0221413E23;
// 1 (Ry) = 2.179872E-18 (J)
eRy = 2.179872E-18; //(J)
// リュードベリ原子単位系でのボルツマン定数
// Boltzmann constant kB = 1.3806488E-23 (J/K)
kB = 1.3806488E-23 / eRy; // (Ry/K)
// 電気素量
chage = 1.60217662E-19;

// *** 状態密度の読み出し ***
X = fscanfMat("Kelvin-Pt-calc.txt");
Edat = X(:,1);
Ddat = 2 * X(:,2);

// *** 計算用 ***
// エネルギー
E = linspace(min(Edat), max(Edat), 10000);
// 温度
tstart = 10; tstep = 10; tend = 2500;
T = [tstart:tstep:tend];
// 状態密度
D = interp1(Edat, Ddat, E, "linear");

// *** フェルミ分布関数 ***
function fermi = fermi(mu, energy, temp)
fermi = 1 ./ (exp((energy - mu) ./ (kB * temp)) + 1)
endfunction

// *** フェルミ分布関数 ***
n = intsplin(E, (D .* fermi(0, E, 0)));
function y = f1(x, temp)
y = intsplin(E, D .* fermi(x, E, temp)) - n
endfunction

// *** 化学ポテンシャル ***
Snum = ones(T);
for i = 1:length(T) do
temp = 1.01 * T(i);
mu1 = fsolve(0, f1);
temp = 0.99 * T(i);
mu2 = fsolve(0, f1);
// 数値計算による電子比熱
Snum(i) = - eRy * (mu1 - mu2) / (0.02 * T(i)) / chage;
end

// 数値計算による電子比熱
plot(T, 1E6 * Snum, "-b");
Y = fscanfMat("Kelvin-Pt-lit.txt");
plot(Y(:,1), Y(:,2), ".b");

// *** グラフの装飾 ***
xlabel("Temperature (K)");
ylabel("Seebeck coefficient (uV/K)");
xgrid(color("gray"));


Sommerfeld展開


Scilabで金属の化学ポテンシャルでは Sommerfeld 展開から得られた化学ポテンシャルが以下のように書かれるとしています。
\begin{equation}
\mu = \epsilon_F - \frac{\pi^2}{6}k_B^2 \frac{D'(\epsilon_F)}{D(\epsilon_F)}T^2
\end{equation}
したがって、その温度微分は以下のようになります。
\begin{equation}
\frac{\mathrm{d}\mu}{\mathrm{d}T} = - \frac{\pi^2}{3}k_B^2 \frac{D'(\epsilon_F)}{D(\epsilon_F)}T
\end{equation}
状態密度が鋭く変化している(D'(εF)が大きい)ほど大きなゼーベック係数を持つことが分かります。

関連エントリ




参考URL




付録


このエントリで使用したファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: AkaiKKR KKR ゼーベック係数 状態密度 分散関係 

AkaiKKRで窒化ガリウム

AkaiKKR(machikaneyama)ウルツ型結晶構造の窒化ガリウムGaNの計算を行いました。

hexGaN.png
Fig.1: ウルツ型構造の窒化ガリウムの状態密度


ウルツ型構造の結晶構造はFig.2のような六方晶です。このCIFはCrystallography Open Databaseからダウンロードすることが出来ます。

2018y04m17d_140952920.png
Fig.2: ウルツ型窒化ガリウムの結晶構造


AkaiKKRの入力ファイルは、以下のようになりました。原子位置(atomicx)はVESTAでAkaiKKRのための基本並進ベクトルの方法で出力したものを利用しました。

c----------------------GaN-----------------------------------
go data/GaN
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
hex 6.01 , 1.6245 , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 1.8 sra mjw nmag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 4 200 0.03
c------------------------------------------------------------
c ntyp
2
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Ga 1 1 0.0 2
31 100
N 1 1 0.0 2
7 100
c------------------------------------------------------------
c natm
4
c------------------------------------------------------------
c atmicx atmtyp
0.333333343a 0.666666687b 0.000000000c Ga
0.666666627a 0.333333313b 0.500000000c Ga
0.333333343a 0.666666687b 0.384999990c N
0.666666627a 0.333333313b 0.884999990c N
c------------------------------------------------------------


関連エントリ




参考URL




付録


このエントリで使用したファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: AkaiKKR machikaneyama KKR 状態密度 

AkaiKKRで磁鉄鉱

AkaiKKR(machikaneyama)でスピネル型の結晶構造を持つ磁鉄鉱(マグネタイト)Fe3O4の計算を行いました。この計算はYanase and Siratori (1984, J. Phys. Soc. Jpn)で報告されている通り、基底状態ではハーフメタリックな強磁性体になります。(ハーフメタルに関してはAkaiKKRでハーフメタルを参照。)

fccMagnetitie-dos.png
Fig.1: 磁鉄鉱の状態密度。基底状態ではハーフメタルになる。


入力ファイル


磁鉄鉱の結晶構造は以下のように図示できます。

2017y09m18d_154337081.png
Fig.2: 磁鉄鉱の結晶構造


一見すると複雑な結晶構造ですが、fcc格子を持つと考えると計算セルの中の原子の数は14個まで減ります。鉄のAサイトの位置は ±(1/8, 1/8, 1/8) で、Bサイトの位置は (0, 0, 1/2), (1/4, 0, 3/4), (0, 1/4, 3/4), (1/4, 3/4, 0) です。酸素のサイトは ±(-u, u±1/4, -u), ±(-u, -u, u±1/4) で磁鉄鉱の場合 u=0.0048 です。結局、入力ファイルは以下のようになります。

c------------------------------------------------------------
go data/fccMagnetite
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
fcc 15.87 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 0.9 sra mjw mag 2nd
c 0.001 1.2 sra mjw mag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 4 200 0.035
c------------------------------------------------------------
c ntyp
3
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Fe1 1 1 0.0 2
26 100
Fe2 1 1 0.0 2
26 100
O 1 1 0.0 2
8 100
c------------------------------------------------------------
c natm
14
c------------------------------------------------------------
c atmicx atmtyp
1/8 1/8 1/8 Fe1
-1/8 -1/8 -1/8 Fe1
0 0 1/2 Fe2
1/4 0 3/4 Fe2
0 1/4 3/4 Fe2
1/4 3/4 0 Fe2
0.2548 0.2548 0.2548 O
-0.2548 -0.2548 -0.2548 O
0.2548 -0.0048 -0.0048 O
-0.2548 0.0048 0.0048 O
-0.0048 0.2548 -0.0048 O
0.0048 -0.2548 0.0048 O
-0.0048 -0.0048 0.2548 O
0.0048 0.0048 -0.2548 O
c------------------------------------------------------------


バンド構造


以下にブロッホスペクトル関数(バンド構造)を示します。バンド構造から見てもハーフメタルになっていることが分かります。

Magnetite-up.png
Magnetite-dn.png
Fig.3-4: バンド構造


関連エントリ




参考URL




付録


このエントリで使用したファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: AkaiKKR machikaneyama KKR 強磁性 ハーフメタル 状態密度 分散関係 

Ubuntuでシンボリックリンクの作成

AkaiKKR(machikaneyama)などのソフトウエアパッケージをインストールする際に、以前のバージョンも残しておきたいと思うことが良くあります。一方で、最新版には簡単にアクセスできるようにもしておきたいです。そこで最新版のディレクトリにシンボリックリンクを作成すると便利です。Windowsで言うところのショートカットです。現在の最新版である January 17, 2018 を ~/kkr/20180117/cpa2002v009c/ にインストールしてあるとして、以下のようにコマンドを打ち込みます。

ln -s ~/kkr/20180117/cpa2002v009c/ ~/kkr/cpa2002v009c


このようにする事によって、あたかも ~/kkr/cpa2002v009c/ に最新版があるかのように操作することが出来ます。新しいバージョンが公開されたら、シンボリックリンクだけを作り直せばOKです。

参考URL




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Ubuntu AkaiKKR machikaneyama 

AkaiKKRのMT球の充填率 その2

AkaiKKRのMT球の充填率 その1では、面心立方構造と体心立方構造の金属のマフィンティン半径を変えながら全エネルギーを計算しました。その結果、マフィンティン球をできるだけ大きくするとき、全エネルギーが低くなることが確認できました。
これに対して、プリミティブセルに2種類以上の原子を含む場合、それぞれの原子のMT半径の比をどのように取るのがよいのかという問題が発生します。今回は、MgOのマグネシウムと酸素のマフィンティン球半径の比を、マフィンティン球が接する条件を保ちながら変化させたときに全エネルギーがどのように変わるかを確認しました。


計算手法


AkaiKKRのMT球の充填率 その1と同様に入力ファイルのテンプレートとシェルスクリプトを用いて、酸素とマグネシウムの比を変化させながら全エネルギーを計算します。

c------------------------------------------------------------
go data/MgO_RMT
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
fcc 7.9582 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 2.0 sra pbe nmag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 4 200 0.035
c------------------------------------------------------------
c ntyp
2
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Mg 1 2 0.0 2
12 100
O 1 RMT 0.0 2
8 100
c------------------------------------------------------------
c natm
2
c------------------------------------------------------------
c atmicx atmtyp
0 0 0 Mg
1/2 1/2 1/2 O
c------------------------------------------------------------


#!/bin/csh -f

## *** 実行ファイル ***
#setenv GFORTRAN_UNBUFFERED_ALL=y
#set EXEC="~/kkr/20170222/cpa2002v009c/specx"
set EXEC="specx"

## *** プロジェクト名 ***
set PREFIX="MgO"

## *** 酸素のMT半径設定値 ***
set RATIO_LIST=( 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 )

## *** ポテンシャルの再利用フラグ ***
## (0: 利用しない, 1: 利用する)
set POTCOPY=0

## *** 計算結果の出力先 ***
set ANALYSIS="./analysis/${PREFIX}.txt"
if ( -e ${ANALYSIS} ) then
cat ${ANALYSIS} >> ${ANALYSIS}.back
endif
echo "MT_ratio Filling(%) Total_energy(Ry)" > ${ANALYSIS}

## *** 繰り返し計算 ***
foreach RATIO ( ${RATIO_LIST} )
set RMT=`echo "${RATIO}*2" | bc -l | sed -e 's/^\./0./g'`
## *** ファイル名 ***
set TEMPLATE="./template/${PREFIX}_Template.in"
set KKRINP="./in/${PREFIX}_${RMT}.in"
set KKROUT="./out/${PREFIX}_${RMT}.out"
set POTFILE="./data/${PREFIX}_${RMT}"
set POTBACK="./data/${PREFIX}"

## 前回のポテンシャルが存在すれば利用
if (( -e ${POTBACK} ) && ( ${POTCOPY} == 1 )) then
cp ${POTBACK} ${POTFILE}
endif

## *** 入力ファイルの作成 ***
sed 's/'RMT'/'${RMT}'/g' ${TEMPLATE} > ${KKRINP}

## *** 第一原理計算の実行 ***
## 最大計算回数
set nummax=20
## 計算回数の初期化
set num=0
## 最初の第一原理計算
${EXEC} < ${KKRINP} > ${KKROUT}
## *** 繰り返し計算 ***
while ( ( ! { grep -q "err= -6." ${KKROUT} } ) && ( $num < $nummax ) )
${EXEC} < ${KKRINP} > ${KKROUT}
@ num++
end

## 前回のポテンシャルが存在すればバックアップ
if ( ${POTCOPY} == 1 ) then
cp ${POTFILE} ${POTBACK}
endif

## *** 結果の出力 ***
set ENE=`grep "total energy" ${KKROUT} | sed -e s/total//g -e s/energy=//g`
set FIL=`grep "volume filling=" ${KKROUT} | sed -e s/volume//g -e s/filling=//g -e s/%//g`
echo ${RATIO} ${FIL} ${ENE} >> ${ANALYSIS}
end


計算結果


計算結果を以下に示します。

MgO.png

Fig.1


酸素とマグネシウムのマフィンティン半径の比rO/rMgが1.2~1.3ぐらいで全エネルギーが最小になっています。
ちなみにマフィンティン球による充填率は、rO/rMg=1のとき最も低くなります。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: AkaiKKR machikaneyama KKR マフィンティン半径 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性PICCPAOPアンプecalj常微分方程式モンテカルロ解析状態密度odeトランジスタインターフェースDOS定電流スイッチング回路PDS5022半導体シェルスクリプト分散関係レベルシフト乱数HP6632Aトランジスタ技術ブレッドボード可変抵抗温度解析I2CR6452A反強磁性バンドギャップ確率論数値積分セミナー偏微分方程式絶縁バンド構造熱設計非線形方程式ソルバシュミットトリガISO-I2CLEDマフィンティン半径GW近似三端子レギュレータLM358A/DコンバータカオスフォトカプラUSBPC817C直流動作点解析サーボ74HC4053アナログスイッチTL431発振回路カレントミラー数値微分単振り子量子力学開発環境補間2ちゃんねるチョッパアンプbzqltyFFT電子負荷アセンブラBSchLDA標準ロジックパラメトリック解析ブラべ格子基本並進ベクトルイジング模型VESTAVCAMaximaSMPewidthGGA仮想結晶近似FET位相図キュリー温度QSGWTLP621ランダムウォーク不規則合金gfortranコバルト相対論失敗談抵抗状態方程式スレーターポーリング曲線ラプラス方程式スピン軌道相互作用スイッチト・キャパシタ六方最密充填構造熱伝導繰り返しcygwinTLP552条件分岐TLP521NE555LM555マントル詰め回路MCUテスタFXA-7020ZR三角波過渡解析ガイガー管自動計測Writer509UPSQNAPダイヤモンドデータロガー格子比熱熱力学平均場近似OpenMPブラウン運動スーパーセルUbuntuフェルミ面差し込みグラフubuntuハーフメタルfsolve最適化第一原理計算固有値問題シュレディンガー方程式最小値awk起電力井戸型ポテンシャルCIFxcrysden最大値結晶磁気異方性PGATeX非線型方程式ソルバ2SC1815等高線OPA2277面心立方構造初期値FSM正規分布interp1ウィグナーザイツ胞フィルタfccL10構造合金BaOウルツ鉱構造CapSense岩塩構造ルチル構造ZnO二相共存磁気モーメント不純物問題電荷密度重積分SICスワップ領域リジッドバンド模型multiplotジバニャン方程式gnuplotc/a全エネルギー半金属デバイ模型edeltquantumESPRESSOノコギリ波フォノン固定スピンモーメントspecx.f等価回路モデル円周率パラメータ・モデルヒストグラム不規則局所モーメントTS-112TS-110直流解析PCExcelシンボルGimp日本語最小二乗法フラクタルマンデルブロ集合縮退クーロン散乱三次元ゼーベック係数キーボード入出力関数フィッティング文字列疎行列Realforceトラックボール線種EAGLE連立一次方程式MBECrank-Nicolson法AACircuit負帰還安定性ナイキスト線図マテリアルデザインP-10化学反応ifort境界条件陰解法熱拡散方程式MAS830LCK1026グラフの分割軸ラベル凡例片対数グラフトランスHiLAPW両対数グラフLMC662PIC16F785ヒストグラム確率論

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ