xyzzyでAkaiKKR-mode

xyzzyAkaiKKRの入力ファイル用のメジャーモードを作りました。と、いってもキーワードとコメントに色が付くだけですが。


001_20130115163802.png

Fig.1: AkaiKKR-modeのスクリーンショット



インストール


  1. AkaiKKR.txtAkaiKKR(拡張子なし)にリネームしてetc/へ置く(参考:xyzzyのキーワード設定方法)
  2. AkaiKKR-mode_l.txtAkaiKKR-mode.lにリネームしてsite-lisp/へ置く
  3. 下記を.xyzzyファイルかsiteinit.lへ書く
  4. xyzzy再起動&再ダンプ

;;; AkaiKKR-mode
(require "AkaiKKR-mode")
;;(pushnew '("\\.in$" . AkaiKKR-mode) *auto-mode-alist* :test 'equal)


3行目のコメントアウトをはずすと、拡張子が.inのファイルを開いたときに自動的にAkaiKKR-modeにするように出来ます。

キーワードファイルについて


キーワードファイルに入力したキーワードのリストはInstalling and Runnig AKAIKKRを参考にしました。このPDFでは、本家であるKKR-Green関数法によるバンド計算Machikaneyama2000の使用に関するメモ、あるいは計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)で紹介されていないようなものも含まれているようです。

私は、これらのオプションがどの程度機能するのかは確認していません。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR CPA xyzzy 

CygwinでAkaiKKR(Machikaneyama)

このエントリの内容は、やや古くなったため新しいエントリを書き直しました。→ AkaiKKRのインストール

AkaiKKR(Machikaneyama)は、大阪大学の赤井久純先生が公開されている第一原理計算パッケージで、KKR-CPA法を採用することで不規則合金のバンド計算が可能です。


NiDOS.png

Fig.1: Niの状態密度


本エントリでは、Cygwinのg77コンパイラを利用してWindows上でAkaiKKRを使う方法について書きます。

Cygwin上のAkaiKKRの利用に関しては、AkaiKKRのBBSの0038
How can I calculate cohesive energy? (in Japanese)
にmakefileの編集だけでコンパイルできるとの報告がありますが、最新版(009c)ではエラーが出てしまいます。そこでmakefileの編集に加えてspecx.fの編集が必要になります。(AkaiKKRのヴァージョンを下げると解決するという報告もあるようですが、おそらくspecx.fの編集の方が正攻法のはずです。)

Cygwin上でのAkaiKKR利用に必要な手順は下記の通りです。
  1. Cygwinのインストール
  2. AkaiKKRのダウンロードと展開
  3. makefileの編集
  4. specx.fの編集とコンパイル
  5. テスト計算(Niの状態密度)



1. Cgywinのインストール


AkaiKKRをコンパイルするためには、FORTRAN77コンパイラが必要です。Windows上で利用できるFORTRAN77コンパイラはいろいろな種類があると思うのですが、今回は無料で利用できるCygwin上のg77コンパイラを利用します。
Cygwinのインストールの方法自体は、ウエブ上にたくさん解説記事があると思うので書きませんが、必ずg77をインストールしてください。デフォルトではインストールされません。


002_20130114153227.png

Fig.2: Cygwinのインストールでかならずg77がインストールされるようにする。全パッケージのインストールを選ぶとg77もインストールされる。


私は面倒くさいので全てインストール(All Install)としました(Fig.2)。これを選ぶとインストールにかなりの時間がかかります。
夜寝る前にインストール作業を開始して、朝起きたらインストールが終わっているか確認するぐらいがいいと思います。

Cygwinのインストールが完了したら、Cygwinのターミナル上で
g77 --version

とタイプしてみます。コンパイラのバージョン情報が表示されればインストール完了です。

2. AkaiKKRのダウンロードと展開


AkaiKKRのダウンロードには『計算コード利用許諾契約』への同意と利用登録が必要です。
利用登録に必要な情報は
  • 名前
  • メールアドレス
  • 所属
の3つだけです。

利用登録が完了すると入力したメールアドレスにパスワードが送信されてきます。
このパスワードを利用して、このエントリを書いている段階での最新版であるcpa2002v009c.tar.gzをダウンロードします。

AkaiKKRのページでは、プログラム本体のほかにマニュアルも配布しています。KKR-Green関数法によるバンド計算Machikaneyama2000の使用に関するメモは印刷して手元においておくとよいと思います。

また市販の書籍として計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)密度汎関数法の発展 -マテリアルデザインへの応用にAkaiKKRに関する言及があります。特に前者は、第一原理計算を"良くわからないけど、とにかくやってみる"という立場の人(私とか)にはよい本だと思います。



KKR-Green関数法によるバンド計算の『プログラムの入手、コンパイルと実行』にならった方法でインストールを進めます。

まず、Cygwinのホームディレクトリにkkrというディレクトリを作成します。Cygwinのホームディレクトリは、普通にインストールを行っていれば『C:\cygwin\home\ユーザー名』となっているはずです。
次に、kkrディレクトリに、ダウンロードしたcpa2002v009c.tar.gzを解凍します。(解凍はマニュアルにあるとおりコマンドを打ち込んでも、解凍レンジなどWindowsの解凍ツールでも、お好きなほうで。)
すると~/kkr/cpa2002v009cというディレクトリができます。

3. makefileの編集


さて、~/kkr/cpa2002v009cのなかをみるとmakefileというファイルがあります。
このファイルをメモ帳などのテキストエディタで開きます。すると5行目に
fort = ifort

と書いてあるのが分かります。これはIntel Fortran Compilerを使ってコンパイルする設定です。今回はg77を使うのでこの行を
fort = g77

と、書き換えます。makefileの編集はこれだけでコンパイルできました。

4. specx.fの編集とコンパイル


KKR-Green関数法によるバンド計算、あるいは計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)のチュートリアルの通りに進めるのならば、ここでソースコードのコンパイルを行います。

コンパイルの方法は簡単でCygwinのターミナル上で~/kkr/cpa2002v009cへ移動し
make

とタイプするだけです。

しかしながら、このままだと下記のようなエラーメッセージが表示されてコンパイルできません。

g77 -O3 -o source/specx.o -c source/specx.f
In file included from source/specx.f:0:
source/specx.f:42: error: size of variable '_BLNK__' is too large
makefile:206: recipe for target `source/specx.o' failed
make: *** [source/specx.o] Error 1


このエラーメッセージの意味するところは、計算に使うメモリの量が多すぎるということです。そこでspecx.fの編集を行います。
~/kkr/cpa2002v009c/source/にあるspecx.fをテキストエディタで開きます。
すると19-20行目に以下の記述があります。

     & (natmmx=22, ncmpmx=12, msizmx=198, mxlmx=3, nk1x=550, nk3x=21,
& msex=201, ngmx=15, nrpmx=650, ngpmx=650, npmx=350, msr=400)


これらのパラメータは『配列の大きさを指定するパラメータ』です。詳しくはKKR-Green関数法によるバンド計算のP29の表5を見てください。

この中でnatmmxncmpmxの2つのパラメータに注目します。これらはそれぞれ計算するブラベー格子の中の『原子の個数』と『原子の種類』の最大値をあらわします。
例えば体心立方構造(bcc)や面心立方構造(fcc)の『原子の個数』は1個、六方最密充填構造(hcp)の場合は2個です。また『原子の種類』の方は、純金属で1種類、二元合金で2種類、三元合金で3種類です。

もともとのspecx.fをみるとブラベー格子に12種類22個の原子を含む結晶の計算ができるような条件となっていますが、本ブログの扱おうとしている対象からは明らかに過剰です。そこで19-20行目を下記の通り書き直します。

     & (natmmx=2, ncmpmx=3, msizmx=198, mxlmx=3, nk1x=550, nk3x=21,
& msex=201, ngmx=15, nrpmx=650, ngpmx=650, npmx=350, msr=400)

これはbcc, fcc, hcp構造の三元合金まで計算できる設定です。

前述の通りCygwinのターミナル上で~/kkr/cpa2002v009cへ移動し
make

とタイプすると、実行ファイルであるspecx.exeが生成されます。

今回は、単位構造の中に3種類2個までの原子を含む設定でコンパイルしましたが、これは逆に言うと、計算したい物質の構造や組成によってspecx.fを編集しなおしてその都度コンパイルを行わなければならないことを意味しています。たとえば、入力ファイルのサンプルとして付いているgamn6mn0as0asは、今回の設定では計算できません。natmmxncmpmxの値をもっと大きくとる必要があるでしょう。

5. テスト計算(Niの状態密度)


それでは、早速計算を行ってみましょう。
最早ここまで来ればKKR-Green関数法によるバンド計算の通りサンプルファイルの計算を行ってもよいのですが、サンプルファイルは自己無撞着計算だけしかやってくれません。そこで標準のニッケル入力ファイルから、状態密度(DOS)計算まで一気にやってしまう入力ファイルを用意しました。

c ***********************************************
c AkaiKKR(Machikaneyama) input file
c for fcc Nickel
c ***********************************************

c *** Self-consistent calculation ***
go data/ni
fcc 6.67 , , , , , ,
0.001 1.2 nrl mjw mag init
update 4 100 0.02
c update 13 100 0.02
1
Ni 1 0 0 2 28 100
1
0 0 0 Ni

c *** Density of States (DOS) calculation ***
dos data/ni
fcc 6.67 , , , , , ,
0.001 1.2 nrl mjw mag 2nd
update 8 100 0.02
c update 13 100 0.02
1
Ni 1 0 0 2 28 100
1
0 0 0 Ni


上記のNi.txtをダウンロードして~/kkr/cpa2002v009c/inへ置きます。Cygwinのターミナル上で~/kkr/cpa2002v009cへ移動して
specx.exe <in/Ni.txt >out/Niout.txt

とタイプすると計算を開始します。計算が完了するまでにはしばらく時間がかかります。

計算が完了すると~/kkr/cpa2002v009c/outNiout.txtが出来るはずです。

計算結果(Niout.txt)にはいろいろな情報が出力されています。
出力された情報の詳細は計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)KKR-Green関数法によるバンド計算を参照していただくとして、今回は計の全状態密度(total DOS)をプロットしました(Fig.1)。磁性を考慮した(mag)計算をしているので、状態密度は上向きと下向きの二つのスピンのものが出力されています。Fig.1の赤のラインと緑のラインがそれぞれに対応します。

その他


たくさんの入力ファイルを順番に計算させるときにはバッチファイルを作ると楽です。Pathを通しておくとよいでしょう。

計算に使うPCの性能は、そこそこのもので充分です。私が使っているのは、いわゆる鼻毛鯖にメモリを8GBまで増設したものです。
自分が寝ている間に計算させるなら、うるさくない所におくことも重要です。というか、PCの性能自体よりも、メインのマシンと別個の物を用意できることのほうが大事だと思います。

OSなしの廉価PCを使う場合ubuntu+ifortでも計算できました(と、いうか私もどちらかというとこっちを使ってる)。この場合もspecx.fの編集は必須です。

テキストエディタはある程度以上高機能なものを使うほうがよいです。ある程度以上の高機能とは、ファイル内の検索が強力なもの(grep, 正規表現)で、大きなサイズのファイルが問題なく開ける必要があります。AkaiKKRの計算結果を出力したテキストファイルは、それなりの容量になり、自分にとって必要な部分を検索機能で探すことになるからです。

私はxyzzyを使っていますが、秀丸サクラエディタとかでもいいと思います。いずれにせよ標準のメモ帳では力不足だと思います。

データの解析には、差し当たりMicrosoft OfficeLibreOfficeのようなOfficeスイートの表計算ソフトがあるとよいでしょう。なおNiの状態密度(Fig.1)はgnuplotで描画しました。



お詫びと注意事項


このエントリは、少し以前に書き溜めておいたものをブログの予約投稿機能をつかって週一更新となるように公開しています。

【2013年7月29日追記】下記の件、AkaiKKRのウエブページが復活していました!

自由電子近似を使った(つまり第一原理計算パッケージを使わなくても計算できる)話を書いてから公開しようと思ってグズグズしていたところ本家であるAkaiKKR(Machikaneyama)のページが公開を停止されたようです。日本バンド(計算)屋さんマップによると(5/17 以降、2013、アクセス不能を確認)との事です。

そのため本エントリにたくさんのリンク切れが出来てしまい申し訳ありません。

AkaiKKR(Machikaneyama)の入手に関してはCMDワークショップへ参加することが最良の手段だと思われますが、計算機ナノマテリアルデザイン手法の開発にて以前のバージョン(cpa2002v007)が公開されているようです。私はcpa2002v007のコンパイルも実行も試していませんので、このバージョンでエントリの内容を再現できるかも確認していません。

なお、この「お詫びと注意事項」を書いている時点とエントリが公開状態になるまでにもタイムラグがあるので、上記の内容が公開時点では的外れになっているかもしれませんが、もしそうなっていたらごめんなさい。

また私は、AkaiKKRの開発者である赤井先生、及び、大阪大学の固体電子論グループとは何の関係も無い一個人であり、自分の固体物理の勉強のためにAkaiKKRを利用させていただこうとしている身です。更に言うなら、バンド計算屋さんではなくただの電子工作趣味人です。従って、このブログの内容の正しさに関しては一切保障できませんので、真似をされる方は自己責任でお願いします。おそらくブログ内の記述には、間違っている点もたくさんあると思います。間違いを見つけられた方は、このブログのコメント欄にてお知らせください。

関連エントリ(予定)




参考URL




付録


このエントリで使用したAkaiKKR(Machikaneyama)の入力ファイルと出力結果を収めたファイルを添付します。


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR CPA cygwin 

Scilabで金属の電子比熱

Scilabで金属の化学ポテンシャルでは非線形方程式ソルバfsolveを使って金属の化学ポテンシャルを温度の関数として求めました。これを利用すると電子系の内部エネルギーが温度ごとに計算できます。

今回は、この電子系のエネルギーを温度で数値微分して電子比熱を計算しました。

001_20130616223024.png


電子比熱


電子比熱は、電子系のエネルギーue(T)を温度Tで微分することによって得られます。

c_e(T) = \frac{\mathrm{d}u_e(T)}{\mathrm{d}T}

電子系の内部エネルギーは

u_e(T) = \int_{-\infty}^{\infty}\epsilon f(\epsilon,T) D(\epsilon) \mathrm{d}\epsilon

f(\epsilon,T) = \frac{1}{\exp{\left(\frac{\epsilon - \mu(T)}{k_B - T}\right)}+1}

として表されるため前回数値データとして得られた化学ポテンシャルμ(T)を代入することによって、温度の関数としてue(T)の数値データを得ることが出来ます。

今回は、微分を数値的に行うことによってue(T)からce(T)を求めます。

数値微分


数値的な微分とは、実際には数値差分です。
微分の定義は

f^{'}(x) = \lim_{\Delta x \to 0} \frac{f(x+\Delta x) - f(x)}{\Delta x}

なのでΔxが充分小さければ

f^{'}(x) \simeq \frac{f(x+\Delta x)-f(x)}{\Delta x}

と引き算(差分)と割り算で計算できます。

プログラミング


電子比熱の微分を差分に置き換えると

c_e(T) = \frac{u_e (T+\Delta T)-u_e(T)}{\Delta T}

差分はベクトルueの中の隣り合うようその引き算からもとまりますが、これを一気にやってくれるScilab関数がdiffです。
今回のソースコードでは、化学ポテンシャルμ(T)を計算するループの中でついでに内部エネルギーue(T)を計算した後
// 数値計算による電子比熱
dUdT = diff(Uenum) / tstep;
Cenum = [0,dUdT];

として差分を求めています。
差分を求める際に、ベクトルのようその数がひとつ減ってしまうので、先頭にT=0のときの値ce=0を補っています。

作成したソースコードはElectronicSpecificHeat_sce.txtです。

結果と考察


電子比熱の温度依存性のグラフをFig.1に示します。
これまでの計算ではリュードベリ原子単位系を使ってきましたが、実験データと直接比較できるように縦軸がモル比熱となるようにしました。


001_20130616223024.png
Fig.1: 自由電子近似を用いたアルミニウムの電子比熱の温度依存性。青実線はゾンマーフェルト展開、赤破線は数値計算による結果。


前回書いたとおりゾンマーフェルト展開から

c_e(T) = \frac{\pi^2}{3} k_B^2 D(\epsilon_F)T

またScilabで差し込みグラフ:金属の比熱で書いたとおり金属の低温での電子比熱は、電子比熱係数γを用いて

c_e(T) = \gamma T

と書ける事が実験から知られています。

これらを比較すると

\gamma = \frac{\pi^2}{3} k_B^2 D(\epsilon_F)

となります。

実験から求められた電子比熱係数はγ=1.35(mJ/mol/K^2)です。これに対して「自由電子近似+ゾンマーフェルト展開」で求めたものはγ=0.91(mJ/mol/K-2)となり、自由電子近似はアルミニウムに対しておおよそ妥当といったところだと思います。(参考:Lin et al. 2008 Phys. Rev. B)

数値微分に関してなのですが、今回は温度の刻み幅を一定にしてしまったため、本当は低温側で刻み幅が温度そのものと比較して小さくならない様になってしまっています。このことは微分の結果をおかしくしてしまう筈なのですが、そもそもほとんど直線になるような式であったため、その効果が見えていないようです。もっと複雑な式を微分するときには注意しなければなりません。

数値計算とゾンマーフェルト展開の比較は、前回と同様です。
自由電子近似を用いたアルミニウムの電子状態では、融点以上までゾンマーフェルト展開は妥当であるといえそうです。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 比熱 電子比熱 数値微分 

Scilabで金属の化学ポテンシャル

Scilabの非線形方程式ソルバfsolveを用いて、積分を含む方程式である金属の化学ポテンシャルを計算しました。
対象とする系は、自由電子近似が良く成り立つはずのアルミニウムとし、数値計算による化学ポテンシャルの計算結果とゾンマーフェルト展開を用いた近似値の比較を行いました。

その結果、この系に関しては、数値計算とゾンマーフェルト展開を用いた計算結果がアルミニウムの融点以上までほとんど同じ結果を示すことがわかりました。

001_20130609181459.pngFig.1: アルミニウムの化学ポテンシャルの温度依存性、ゾンマーフェルト展開を用いた近似(青実線)とfsolveを利用した数値計算(赤点線)。インセットは自由電子近似による電子の状態密度(青)、低温時(10K,緑)と高温時(赤)の電子の分布は絶対零度の状態密度に各温度のフェルミ分布関数を掛けたもの、エネルギーの原点はフェルミ準位。



金属電子論と数値計算


Scilabで差し込みグラフ:金属の比熱では実験的に求められたアルミニウムの電子比熱係数γ=1.35(mJ/mol/K^2)という値をもとに電子比熱の温度依存性をグラフにプロットしました。

これに対して、次は理論的にアルミニウムの電子比熱を予想してみます。
差し当たり今回は、電子比熱を計算する前段階として、アルミニウムの化学ポテンシャルを自由電子近似と数値計算から求めます。

金属の電子構造を計算するに当たり、次の2つの近似がよく使われます。
  • 自由電子近似(free electron approximation)
  • ゾンマーフェルト展開(Sommerfeld expansion)

自由電子近似はナトリウムやアルミニウムといった比較的単純な金属に対しては良い近似になることが知られています。また、ゾンマーフェルト展開は温度が充分低い条件で成り立ちます。

ゾンマーフェルト展開は、複雑な複雑な積分の計算を、比較的計算しやすい解析的な式に変形するための級数展開なのですが、Scilabなどで数値計算を行うことでこの近似に頼らずに電子比熱が計算できます。
今回は「自由電子近似+ゾンマーフェルト展開」と「自由電子近似+数値計算」の両方を行い比較します。

金属の化学ポテンシャル


金属の電子比熱は、以下の式で定義されます。

c_e = \frac{\mathrm{d} u_e}{\mathrm{d} T}

ここで、電子系の内部エネルギーは

u_e = \int^{\infty}_{-\infty}\epsilon f(\epsilon) D(\epsilon)\rm{d}\epsilon

さらにD(ε)は電子の状態密度(density of states)でf(ε)はフェルミ・ディラック分布関数です。

f(\epsilon) = \frac{1}{\exp{\left(\frac{\epsilon-\mu}{k_B - T}\right)}+1}

まずD(ε)ですが、これは金属の個性を現す重要なパラメータなのですが、自由電子近似と金属の体積、価電子数から求めることが出来ます(後述)。
一方で、フェルミ・ディラック分布関数内の化学ポテンシャルμは、温度の関数となり少々厄介です。本エントリの主題は、この化学ポテンシャルをScilabの非線形方程式ソルバーを用いて数値的に計算し、ゾンマーフェルト展開によって近似的に求められた化学ポテンシャルと比較を行うことです。

自由電子近似


詳細には立ち入りませんが、自由電子近似による金属の状態密度とフェルミエネルギー(絶対零度における化学ポテンシャル)は、リュードベリ原子単位系では以下の様に表されます。

D(\epsilon) = \frac{V}{2 \pi^2} \epsilon^{\frac{1}{2}}

\epsilon_F = \left( \frac{3 \pi^2 n}{V} \right)^{\frac{2}{3}}

ここでVは単位物質量あたりの体積です。具体的には格子体積やモル体積を用います。今回は差し当たり単位原子辺りの体積(atomic volume)を用いて計算します。nは電子の数で、今回は単位原子辺りの電子の数なので価電子数と一致します。(参考:金持 徹著固体電子論)

自由電子近似は、ナトリウムやアルミニウムでは良く成り立ちますが、遷移金属などには適用できません。

ゾンマーフェルト展開


これも詳細には立ち入らず結論だけ書きます。
ゾンマーフェルト展開の結果、内部エネルギーue、化学ポテンシャルμ、電子比熱ceは

u_e = u_{e0} + \frac{\pi^2}{6}k_B^2 D(\epsilon_F)T^2

\mu = \epsilon_F - \frac{\pi^2}{6}k_B^2 \frac{D^{'}(\epsilon_F)}{D(\epsilon_F)}T^2

c_e = \frac{\pi^2}{3}k_B^2 D(\epsilon_F)T

となります。(参考:無機化学特論I ゾンマーフェルト展開について)

ゾンマーフェルト展開は充分低温のみで適用できます。(参考:Sommerfeld 展開の導出- 武智公平)

自由電子近似+ゾンマーフェルト展開


ゾンマーフェルト展開に自由電子近似の結果を代入すると、金属の化学ポテンシャルや電子比熱の値が具体的に計算できます。

\mu = \epsilon_F \frac{\pi^2 k_B}{12 \epsilon_F}T^2

c_e = \frac{\pi^2 k_B^2 n}{2 \epsilon_F}T

押さえておかなくてはならない重要な点は、「自由電子近似」と「ゾンマーフェルト展開」は相互に全く関係の無い別の近似であり、適用できる条件もそれぞれ異なるということです。

数値計算


さて、本題です。
ここまでで「自由電子近似+ゾンマーフェルト展開」によって化学ポテンシャル(と電子比熱)を計算する準備が整いました。
ここからは、ゾンマーフェルト展開を使わずに化学ポテンシャルを計算することを考えます。

価電子数は温度に依存しないので

n = \int_{-\infty}^{\infty} D(\epsilon)f(\epsilon)\mathrm{d}\epsilon

をμに付いて解く事によって、化学ポテンシャルμをそれぞれの温度で計算することが出来ます。
この式は、解析的には解けないため、Scilabの非線形ソルバfsolveを利用します。

作成したソースコードはChemicalPotential_sce.txtです。

計算結果(自由電子近似+数値計算)


Fig.1に化学ポテンシャルの計算結果を示します。自由電子近似を用いたアルミニウムの場合、ゾンマーフェルト展開を用いた近似と数値積分を利用して求めた結果はほとんど変わらない結果を示します。(高温で若干ずれているのは、本当にゾンマーフェルト展開が妥当でなくなっていっているのか、数値計算の精度の問題かはちょっと自信がないです。計算するエネルギー範囲が狭いせいではなさそう・・・。)

アルミニウムの融点は660.32 ℃(933.47 K)なので、固体の間はゾンマーフェルト展開が破綻することは無いであろう事がわかります。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 非線形方程式ソルバ 電子比熱 数値積分 

Scilabで差し込みグラフ:金属の比熱

Scilabで数値積分:固体の比熱では、積分を含む方程式としてあらわされるデバイの比熱式を計算しました。

今回は、アルミニウムの比熱に関して、デバイの比熱式に加えて金属の電子比熱の計算も行い、ほとんどの温度領域においてはデバイの比熱式で表される格子振動の寄与が支配的であり、しかしながら、数ケルビン程度の低温では電子比熱の影響が大きくなることを確認しました。


001_20130606070707.png
Fig.1: アルミニウムの格子比熱・電子比熱の温度依存性のシミュレーション。極めて低温では電子比熱(赤破線)が支配的になるものの、それ以外では格子比熱(青実線)に比べ電子比熱は無視できるほど小さい。



金属の比熱


Scilabで数値積分:固体の比熱では、積分を含む方程式としてあらわされるデバイの比熱式を計算しました。これは、固体の原子の熱振動に起因する熱容量で、格子比熱と呼びます。固体金属の比熱もまた、デバイの比熱式でほとんど問題なく計算できます。

しかしながら、金属の場合は、わずかながら伝導電子に起因する電子比熱も存在します。これは常温では、格子比熱と比較して無視できるほど小さいのですが、数ケルビン程度の極めて低い温度では格子比熱が急速に小さくなるため、電子比熱の寄与が相対的に大きくなります。

今回は、このアルミニウムに関して格子比熱と電子比熱の大きさをScilabで計算し、差し込みグラフ(インセットグラフ)を描画してみます。

プログラミング


格子比熱の計算には、Scilabで数値積分:固体の比熱で計算を行ったデバイモデルを用います。

C_l(T) = 9 R \left( \frac{T}{\Theta_D} \right)^3 \int^{\Theta_D / T}_{0}\frac{x^4 e^x}{(e^x - 1)^2}{\rm d}x

アルミニウムのデバイ温度は428Kとします。(参考:デバイ模型:wikipedia)

電子比熱は、低温では単純に温度に比例することが知られていて、電子比熱係数γを用いて以下のようにあらわします。

C_e(T) = \gamma T

アルミニウムの電子比熱係数は1.35mJ/mol/K^2です。(参考:第4講- 金属の基本物性の電子論-(PDF):志賀@高槻)

// アルミニウムの電子比熱係数(J/mol/K^2)
egamma = 1.35e-3;
// アルミニウムのデバイ温度
dt = 428;
// 気体定数 (J/K/mol)
r = 8.314

// 格子比熱 (Debye model)
function Cl = Cl(T)
Cl = 9 * r * ((T ./ dt) .^ 3) .* integrate('(x .^ 4) .* exp(x) ./ ((exp(x) - 1) .^ 2)','x',0,dt ./ T);
endfunction
// 電子比熱
function Ce = Ce(T)
Ce = egamma .* T;
endfunction

// 高温までのプロット
// 温度ベクトル
T = [1:1:500];
// 格子比熱のプロット
// 絶対零度の計算は出来ないので後から補う
plot([0,T],[0,Cl(T)],'-b');
// 電子比熱のプロット
plot([0,T],[0,Ce(T)],'--r');
legend(['Lattice specific heat';'Electronic specific heat'],2);
xlabel("Temperature (K)");
ylabel("Specific heat (J/K/mol)");

// 低温部分のプロット
xsetech([0.4,0.32,0.5,0.5]);
// 温度ベクトル
T = [1:0.1:10];
// 格子比熱のプロット
// 絶対零度の計算は出来ないので後から補う
plot([0,T],[0,Cl(T)],'-b');
// 電子比熱のプロット
plot([0,T],[0,Ce(T)],'--r');
xlabel("Temperature (K)");
ylabel("Specific heat (J/K/mol)");


格子比熱や電子比熱は、高温までの計算と低温のみの計算の2回の計算を行うため、あらかじめfunctionを用いて関数化してあります。(参考:Scilab入門―電気電子工学で学ぶ数値計算ツール)

差込グラフ(インセットグラフ)は、差し込むほうのグラフを後から小さいサイズで上書きすることで作成することが出来ます。
差し込むグラフの位置とサイズはxsetechで指定することが出来ます。(参考:コマンドxsetech([x座標始点,y座標始点,幅,高さ])を繰り返し使うことで一つのウィンドウに複数のグラフを記述することができます.及びxsetech - プロットのためのグラフィックスウィンドウの サブウィンドウを設定)

また、複数のパネルに分割する場合はsubplotを利用するほうが簡単かもしれません。

// 高温までのプロット
subplot(2,1,1);
// 温度ベクトル
T = [1:1:500];
// 格子比熱のプロット
// 絶対零度の計算は出来ないので後から補う
plot([0,T],[0,Cl(T)],'-b');
// 電子比熱のプロット
plot([0,T],[0,Ce(T)],'--r');
legend(['Lattice specific heat';'Electronic specific heat'],2);
xlabel("Temperature (K)");
ylabel("Specific heat (J/K/mol)");

// 低温部分のプロット
subplot(2,1,2);
// 温度ベクトル
T = [1:0.1:10];
// 格子比熱のプロット
// 絶対零度の計算は出来ないので後から補う
plot([0,T],[0,Cl(T)],'-b');
// 電子比熱のプロット
plot([0,T],[0,Ce(T)],'--r');
xlabel("Temperature (K)");
ylabel("Specific heat (J/K/mol)");


プロット部分を上記のものに差し替えた場合、グラフは以下のようになります。


002_20130606074346.png

Fig.2: subplotで分割した場合


ただし、下のパネルの横軸ラベルが切れてしまうので、きれいな図を描きたい場合は、やはりxsetechのほうがオススメです。

またScilabで関数フィッティング:金属の電気抵抗のように複数のグラフウインドウをひとつのプログラムから立ち上げたい場合は、subplotの代わりにscfでグラフ番号を指定します。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 差し込みグラフ インセットグラフ 比熱 電子比熱 格子比熱 デバイモデル 

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

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

最新コメント
リンク

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