スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

DFT計算用データベース MatNavi

AkaiKKR(machikaneyama)ecalj, Quantum ESPRESSO(PWscf)といった第一原理計算パッケージを用いるには、色々なデータベースがあると便利です。

物質・材料研究機構(NIMS)が、この目的に適したデータベースを公開しています。
利用にはユーザー登録が必要ですが無料です。ユーザー登録自体はとても簡単です。今回は、雰囲気を紹介するためにいくつかのデータを見てみます。


状態図


実際に計算する物質が組成や温度によってどのような相を撮るのかを確認しておくことは重要です。計算状態図データベース(CPDDB)には、熱力学パラメータから計算された状態図が収録されています。例えば、以下に示すのは Fe-Cr の二元系状態図です。

crfe_lee.png
Fig.1: 計算状態図データベース(CPDDB)で見ることが可能な状態図の一例、 Fe-Cr 二元系状態図。


結晶構造


実際の第一原理計算には、格子定数や原子位置などの詳細な結晶構造データが必要です。

無機材料データベース(AtomWork)では、これらのパラメータやX線回折パターンを表示することが出来ます。

2018y06m10d_005628909.jpg
Fig.2: 無機材料データベース(AtomWork)で見ることの出来る結晶構造の一例、MgAl2O4の結晶構造。ただし、実際にはCIF形式のファイルをダウンロードしてVESTAなどで描画するほうが便利なことが多い。


更に結晶構造を記述するCIF形式のファイルをダウンロードすることも可能です。色々な第一原理計算パッケージがCIFから入力ファイルを作成する機能を持っています。たとえば本ブログでもCIFからecalj入力の作成, CIFからecalj入力の作成 その2です。

またCIFから直接入力ファイルを作る機能が無いDFTパッケージであってもVESTAなどを経由して比較的に簡単に入力ファイルを作成することが出来る場合も多いです。例えばVESTAでAkaiKKRのための基本並進ベクトルPWscfの入力作成補助などです。

バンド構造


第一原理計算を実行する際、既知のバンド構造の計算結果と比較して、答え合わせを行っておくことはとても重要です。電子構造計算データベース (CompES-X)を利用してバンド構造や状態密度などを確認することが出来ます。

2018y06m09d_234430344.png
Fig.3: 電子構造計算データベース (CompES-X)で確認できるバンド構造の一例、PbTeのバンド構造。


関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


 ↓ この記事が面白かった方は「拍手」をお願いします。
スポンサーサイト

tag: AkaiKKR machikaneyama ecalj QuantumESPRESSO 状態図 分散関係 状態密度  

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)でスピネル型の結晶構造を持つ磁鉄鉱(マグネタイト)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 強磁性 ハーフメタル 状態密度 分散関係 

ecaljのgetsyml.py

ecaljには、バンド分散を描画するときの対称性のよいパスを自動的に生成してくれる getsyml.py というpythonスクリプトが ~/ecalj/GetSyml/ にあります。そのディレクトリの README に、以下のようなインストール方法が書いてあるのですが、どういう訳か私の環境では上手く行きませんでした。どうやら私のubuntuのpython環境に何かかの問題があるようでした。

===========================
Requirement and Install:

1.seekpath
>git clone https://github.com/giovannipizzi/seekpath/
>python setup.py install

2.matplotlib for 3D plot
> python -m pip install --update pip #pip update
> pip install matplotlib

3.spglib for crystal structure symmetry
>git clone https://github.com/atztogo/spglib.git
>python setup.py install --user
--user install it locally.


そこでAnaconda で Python 環境をインストールするを参考にして、わたしのubuntuにPython 2.7をインストールしたところ getsyml.py が使えるようになりました。


Anaconda の Python 2.7 のセットアップ


まずAnaconda の Python 2.7 をインストールします(3.xではありません)。Anacondaのダウンロードページから、インストールスクリプトをダウンロードし、実行します。

cd ~
wget https://repo.continuum.io/archive/Anaconda2-5.0.1-Linux-x86_64.sh
bash Anaconda2-5.0.1-Linux-x86_64.sh


すると端末上に、対話型のインストーラーが表示されるので、言われるがままに進めます。最後にAnacondaのpythonをPATHに加えるか聞かれるので yes と答えます。
この段階だと、単純に .bashrc に追記しただけなので source コマンドで .bashrc を再読み込みさせた後 python のバージョンを確認します。

source ~/.bashrc
python --version


以下のように Ananaconda でインストールされたものが表示されていれば成功です。

Python 2.7.14 :: Anaconda, Inc.


seekpath のセットアップ


Python 2.7 のセットアップが完了したら、次に seekpath のセットアップをします。
以下のコマンドを順番に端末に入力します。

cd ~
git clone https://github.com/giovannipizzi/seekpath/
cd seekpath/
python setup.py install


matplotlib のセットアップ


私の環境では特に何もしなくても大丈夫でした。Anacondaではデフォルトでmatplotlibが入ってる?

spglib のセットアップ


以下のコマンドを順番に端末に入力します。

cd ~
git clone https://github.com/atztogo/spglib.git
cd spglib/python/
python setup.py install --user


getsyml.pyの場所をパスに追加


私は ~/ecalj/GetSyml/ をパスに追加しました。
~/.bashrc に以下を追記します。

export PATH="$HOME/ecalj/GetSyml:$PATH"


テスト計算


CIFからecalj入力の作成CIFからecalj入力の作成 その2のセットアップが完了しているという前提で、シリコンのCIFファイルからバンド計算まで一気にやってみます。適当なディレクトリ、例えば ~/ecalj/project/Si-GetSyml/ で以下の順に実行します。

cp ~/cif2cell-1.2.10/cifs/Si.cif si.cif
cif2ctrl.sh si
getsyml.py si
lmfa si
mpirun -np 2 lmf-MPIK si
job_band si -np 2


getsyml.py si を実行すると以下のようなグラフィカルなウインドウが立ち上がります。

Screenshot from 2017-11-15 003A553A41

Fig.1: getsyml.py で得られるブリルアンゾーンの図


最後に job_band si -np 2を実行するとバンド分散の図が得られます。

Screenshot from 2017-11-15 003A583A51

Fig.2: シリコンのバンド構造


ecaljでシリコンのバンド構造(LDA計算)で得られたものと同じバンド分散結果が得られていることが分かります。

関連エントリ




参考URL




フィードバック



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

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


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


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

tag: ecalj 分散関係 

AkaiKKRのspc計算の入出力フォーマット その2

AkaiKKRのspc計算の入出力フォーマットでは、AkaiKKR(machikaneyama)の26 Augsut 2015のバージョンを利用して、それよりも以前の 22 May 2015 のバージョンと同じspcの入出力フォーマットに戻す方法を書きました。入力フォーマットの変更は簡単だったのですが、出力フォーマットの指定方法の変更は、少し面倒くさい感じでした。最新版の 22 July 2016 バージョンのAkaiKKRでは、この出力フォーマットの指定方法も簡単になりました。

具体的には source/spmain.firdfmt の値を変更することによって入力のフォーマットを、また、 iwrtfmt の値を変更することによって出力のフォーマットを変更できます。22 May 2015 のバージョンと同じspcの入出力フォーマットに戻すには、両方とも 1 を指定します。デフォルトの 3 を指定しておくほうがバンド構造のテスト計算には便利かもしれません(参考:AkaiKKRでバンド構造(分散関係))。しかしながら、ねがてぃぶろぐでは、今後とも基本的に22 May 2015 のバージョンと同じ入出力フォーマットで計算をすることにします。

関連エントリ




参考URL




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama 分散関係 

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

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

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。