AkaiKKRのewidth その1

AkaiKKR(machikaneyama)の入力ファイルで指定するパラメータの中で、今回は ewidth と edelt について書きます。

akaikkr-handson-2-140815051743-phpapp02.png

Fig.1: 状態密度と複素数エネルギーメッシュの関係。第5回CCMSハンズオン(ソフトウェア講習会): AkaiKKRチュートリアル 2. AkaiKKRの実習より



go計算のewidthとedelt


AkaiKKR(machikaneyama)では、他の多くのDFTコードとは異なりKKR法をつかいます。したがって、計算に必要なパラメータが違います。例えばAkaiKKRには、平面波法で重要になるカットオフエネルギーのパラメータは存在しません。その一方でewidthやedeltといったパラメータを指定する必要が出てきます。今回はewidthとedeltについてまとめます。

Fig.1に示したewidthとedeltと状態密度の関係を表した図は、第5回CCMSハンズオン(ソフトウェア講習会): AkaiKKRチュートリアル 2. AkaiKKRの実習のP11に載っているものです。ewidth はフェルミエネルギーから負の方向にどのエネルギーまで計算するかを指定します。この値は、価電子バンドを完全に含むように選ばなければなりません。一方、ewidthの方は、今回は深入りしませんが、ほとんどの場合 0.001 Ry としておけば問題ないようです。

   complex energy mesh
1( -2.0000, 0.0000) 2( -1.9998, 0.0041) 3( -1.9989, 0.0095)
4( -1.9965, 0.0168) 5( -1.9913, 0.0263) 6( -1.9810, 0.0388)
7( -1.9620, 0.0546) 8( -1.9286, 0.0742) 9( -1.8730, 0.0976)
10( -1.7856, 0.1238) 11( -1.6574, 0.1507) 12( -1.4839, 0.1750)
13( -1.2692, 0.1926) 14( -1.0291, 0.1999) 15( -0.7873, 0.1954)
16( -0.5679, 0.1804) 17( -0.3876, 0.1581) 18( -0.2518, 0.1327)
19( -0.1569, 0.1076) 20( -0.0946, 0.0849) 21( -0.0555, 0.0657)
22( -0.0319, 0.0501) 23( -0.0181, 0.0378) 24( -0.0101, 0.0284)
25( -0.0056, 0.0212) 26( -0.0031, 0.0157) 27( -0.0017, 0.0116)
28( -0.0009, 0.0086) 29( -0.0005, 0.0064) 30( -0.0003, 0.0047)
31( -0.0001, 0.0035) 32( -0.0001, 0.0026) 33( -0.0000, 0.0019)
34( -0.0000, 0.0014) 35( -0.0000, 0.0010)


上に示したのがgo計算用に用いられた energy contour の例です。これは ewidth=2.0 Ry, edelt=0.001 Ry としたときの複素エネルギーの経路で、計算のアウトプットのかなりはじめの方(16行目ぐらいから)に出力されます。Fig.1のように(-ewidth, 0)から半楕円状に(0, edelt)までの経路をとっていることが分かります。go計算では、この経路上の35点を計算していますが、この35という分解能は source/specx.f の97行目ぐらいにある mse0 で指定されています。

contour.png
Fig.2: go計算の実際の energy contour (ewidth=2.0 Ry, edelt=0.001 Ry)


DOS計算のewidthとedelt


go計算のewidthは必ず価電子バンドを完全に含むように選ばなければならないと書きましたが、そのことを確認するためには、結局のところ状態密度(またはバンド構造)を計算してみなければなりません。この際の energy contour はgo計算の場合と少し異なります。
以下に示すのは、ewidth=2.67 Ry, edelt=0.001 Ry としたときのdos計算のときの energy contour の例です。go計算のときは半楕円状でしたが、dos計算の場合は、エネルギーの虚数成分は常に edelt の値で一定です。実数成分のエネルギー範囲も異なっており -0.75*ewidth から 0.25*ewidth の範囲となっています。(逆に言えばgo計算でのewidthの範囲とまったく同じ範囲の状態密度を描画しようとする場合、ewidth/0.75を新しいewidthとしなければならないという事です。本エントリのgo計算でのewidth=2.0なのでdosでは2.0/0.75=2.67としました。)この 0.75 という値は source/cemesr.f の ref の値を編集することにより変更することができます。エネルギーの分解能も201点となっており、これは source/specx.f の20行目ぐらいにある msex で指定されています。これらはdos計算でもspc計算でも共通です。

   complex energy mesh
1( -2.0025, 0.0010) 2( -1.9891, 0.0010) 3( -1.9758, 0.0010)
4( -1.9625, 0.0010) 5( -1.9491, 0.0010) 6( -1.9357, 0.0010)
7( -1.9224, 0.0010) 8( -1.9090, 0.0010) 9( -1.8957, 0.0010)
10( -1.8823, 0.0010) 11( -1.8690, 0.0010) 12( -1.8557, 0.0010)
13( -1.8423, 0.0010) 14( -1.8289, 0.0010) 15( -1.8156, 0.0010)
16( -1.8022, 0.0010) 17( -1.7889, 0.0010) 18( -1.7755, 0.0010)
19( -1.7622, 0.0010) 20( -1.7489, 0.0010) 21( -1.7355, 0.0010)
22( -1.7222, 0.0010) 23( -1.7088, 0.0010) 24( -1.6954, 0.0010)
25( -1.6821, 0.0010) 26( -1.6687, 0.0010) 27( -1.6554, 0.0010)
28( -1.6421, 0.0010) 29( -1.6287, 0.0010) 30( -1.6153, 0.0010)
31( -1.6020, 0.0010) 32( -1.5886, 0.0010) 33( -1.5753, 0.0010)
34( -1.5619, 0.0010) 35( -1.5486, 0.0010) 36( -1.5353, 0.0010)
37( -1.5219, 0.0010) 38( -1.5086, 0.0010) 39( -1.4952, 0.0010)
40( -1.4819, 0.0010) 41( -1.4685, 0.0010) 42( -1.4551, 0.0010)
43( -1.4418, 0.0010) 44( -1.4284, 0.0010) 45( -1.4151, 0.0010)
46( -1.4017, 0.0010) 47( -1.3884, 0.0010) 48( -1.3750, 0.0010)
49( -1.3617, 0.0010) 50( -1.3483, 0.0010) 51( -1.3350, 0.0010)
52( -1.3216, 0.0010) 53( -1.3083, 0.0010) 54( -1.2950, 0.0010)
55( -1.2816, 0.0010) 56( -1.2683, 0.0010) 57( -1.2549, 0.0010)
58( -1.2416, 0.0010) 59( -1.2282, 0.0010) 60( -1.2148, 0.0010)
61( -1.2015, 0.0010) 62( -1.1882, 0.0010) 63( -1.1748, 0.0010)
64( -1.1614, 0.0010) 65( -1.1481, 0.0010) 66( -1.1347, 0.0010)
67( -1.1214, 0.0010) 68( -1.1080, 0.0010) 69( -1.0947, 0.0010)
70( -1.0814, 0.0010) 71( -1.0680, 0.0010) 72( -1.0547, 0.0010)
73( -1.0413, 0.0010) 74( -1.0280, 0.0010) 75( -1.0146, 0.0010)
76( -1.0012, 0.0010) 77( -0.9879, 0.0010) 78( -0.9746, 0.0010)
79( -0.9612, 0.0010) 80( -0.9479, 0.0010) 81( -0.9345, 0.0010)
82( -0.9212, 0.0010) 83( -0.9078, 0.0010) 84( -0.8944, 0.0010)
85( -0.8811, 0.0010) 86( -0.8678, 0.0010) 87( -0.8544, 0.0010)
88( -0.8411, 0.0010) 89( -0.8277, 0.0010) 90( -0.8144, 0.0010)
91( -0.8010, 0.0010) 92( -0.7876, 0.0010) 93( -0.7743, 0.0010)
94( -0.7610, 0.0010) 95( -0.7476, 0.0010) 96( -0.7343, 0.0010)
97( -0.7209, 0.0010) 98( -0.7076, 0.0010) 99( -0.6942, 0.0010)
100( -0.6808, 0.0010) 101( -0.6675, 0.0010) 102( -0.6542, 0.0010)
103( -0.6408, 0.0010) 104( -0.6275, 0.0010) 105( -0.6141, 0.0010)
106( -0.6008, 0.0010) 107( -0.5874, 0.0010) 108( -0.5740, 0.0010)
109( -0.5607, 0.0010) 110( -0.5474, 0.0010) 111( -0.5340, 0.0010)
112( -0.5207, 0.0010) 113( -0.5073, 0.0010) 114( -0.4940, 0.0010)
115( -0.4806, 0.0010) 116( -0.4673, 0.0010) 117( -0.4539, 0.0010)
118( -0.4405, 0.0010) 119( -0.4272, 0.0010) 120( -0.4139, 0.0010)
121( -0.4005, 0.0010) 122( -0.3872, 0.0010) 123( -0.3738, 0.0010)
124( -0.3605, 0.0010) 125( -0.3471, 0.0010) 126( -0.3337, 0.0010)
127( -0.3204, 0.0010) 128( -0.3071, 0.0010) 129( -0.2937, 0.0010)
130( -0.2804, 0.0010) 131( -0.2670, 0.0010) 132( -0.2537, 0.0010)
133( -0.2403, 0.0010) 134( -0.2269, 0.0010) 135( -0.2136, 0.0010)
136( -0.2003, 0.0010) 137( -0.1869, 0.0010) 138( -0.1736, 0.0010)
139( -0.1602, 0.0010) 140( -0.1469, 0.0010) 141( -0.1335, 0.0010)
142( -0.1201, 0.0010) 143( -0.1068, 0.0010) 144( -0.0935, 0.0010)
145( -0.0801, 0.0010) 146( -0.0668, 0.0010) 147( -0.0534, 0.0010)
148( -0.0401, 0.0010) 149( -0.0267, 0.0010) 150( -0.0134, 0.0010)
151( 0.0000, 0.0010) 152( 0.0133, 0.0010) 153( 0.0267, 0.0010)
154( 0.0400, 0.0010) 155( 0.0534, 0.0010) 156( 0.0667, 0.0010)
157( 0.0801, 0.0010) 158( 0.0934, 0.0010) 159( 0.1068, 0.0010)
160( 0.1201, 0.0010) 161( 0.1335, 0.0010) 162( 0.1468, 0.0010)
163( 0.1602, 0.0010) 164( 0.1736, 0.0010) 165( 0.1869, 0.0010)
166( 0.2003, 0.0010) 167( 0.2136, 0.0010) 168( 0.2269, 0.0010)
169( 0.2403, 0.0010) 170( 0.2536, 0.0010) 171( 0.2670, 0.0010)
172( 0.2803, 0.0010) 173( 0.2937, 0.0010) 174( 0.3070, 0.0010)
175( 0.3204, 0.0010) 176( 0.3337, 0.0010) 177( 0.3471, 0.0010)
178( 0.3604, 0.0010) 179( 0.3738, 0.0010) 180( 0.3871, 0.0010)
181( 0.4005, 0.0010) 182( 0.4139, 0.0010) 183( 0.4272, 0.0010)
184( 0.4405, 0.0010) 185( 0.4539, 0.0010) 186( 0.4672, 0.0010)
187( 0.4806, 0.0010) 188( 0.4939, 0.0010) 189( 0.5073, 0.0010)
190( 0.5206, 0.0010) 191( 0.5340, 0.0010) 192( 0.5473, 0.0010)
193( 0.5607, 0.0010) 194( 0.5740, 0.0010) 195( 0.5874, 0.0010)
196( 0.6007, 0.0010) 197( 0.6141, 0.0010) 198( 0.6275, 0.0010)
199( 0.6408, 0.0010) 200( 0.6542, 0.0010) 201( 0.6675, 0.0010)


関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR ewidth edelt 状態密度 DOS 

AkaiKKRとecaljでCuGaTe2 その2

AkaiKKRとecaljでCuGaTe2 その1ではAkaiKKR(machikaneyama)を用いてCuGaTe2の状態密度を計算しました。今回は、さらにecaljを用いて、LDAとGW近似の計算を行いました。その結果、AkaiKKRがCuGaTe2のバンドギャップを過小評価するのは、LDA/GGAに共通する問題であろう事、また、ecaljを用いてGW近似を適用すると、バンドギャップが改善されることが分かりました。

AkaiKKRLDAGW
Bandgap~0.4 eV0.1586 eV1.0962 eV
table.1: AkaiKKR(GGA), ecalj(LDA), ecalj(GW近似)のそれぞれのバンドギャップ



CuGaTe2のバンドギャップ問題


AkaiKKRとecaljでCuGaTe2 その1では、AkaiKKR(machikaneyama)を用いてCuGaTe2の状態密度を計算しました。その結果、得られたバンドギャップの大きさは過小評価となりました。これがLDA/GGAに起因する問題ならば、ecaljのLDA計算でも過小評価となり、かつ、GW近似を適用することで改善が見られるはずです。そこで今回は、LDAとGW近似の両方でCuGaTe2の計算を行います。

ecaljの結晶構造ファイル


ecaljの実行手順(LDA計算)ecaljの実行手順(GW近似)で書いたとおりecaljのために最低限必要なのは結晶構造ファイル ctrls.cugate2 だけです。ecaljの結晶構造ファイルは、基本並進ベクトルで格子を指定し、直交座標系で原子位置を指定する必要があります。

AkaiKKRではブラベ格子の種類は「ブラベ格子のキーワード」を指定する方法と「基本並進ベクトル」を指定する方法の2種類があります(参考: AkaiKKRのブラベ格子AkaiKKRの基本並進ベクトル その1AkaiKKRの基本並進ベクトル その2)。また、原子位置も「直交座標系」で与える方法と「分率座標系(fractional coordinate)」で与える方法の2種類あります。

もしもAkaiKKRの入力ファイルが、ブラベ格子のキーワードを利用している、あるいは、原子位置が分率座標系で与えられている場合、ecaljの入力フォーマットに合うように直す必要があります。AkaiKKRの入力ファイルからecaljの入力ファイルを作る方法としては、いったんcifファイルなどを経由してVESTAを利用する方法が正攻法でしょう(参考: VESTAでAkaiKKRのための基本並進ベクトル)。

しかし、もしもすでにAkaiKKR用の入力ファイルがあるのなら、ワンステップだけとはいえ、余計な手間がかかっている感があります。幸いにして、AkaiKKRの出力には、基本並進ベクトルと直交座標系での原子位置が出力されるので、今回はこれらを利用してecaljのための入力ファイルを作ります。

具体的には、AkaiKKRとecaljでCuGaTe2 その1の出力結果から、以下の基本並進ベクトルと原子位置の部分を抜き出します。

   primitive translation vectors
a=( -0.50000 0.50000 0.99600)
b=( 0.50000 -0.50000 0.99600)
c=( 0.50000 0.50000 -0.99600)


   atoms in the unit cell
position= 0.23703000 0.25000000 0.24900000 type=Te
position= 0.76297000 0.75000000 0.24900000 type=Te
position= 0.75000000 0.23703000 0.74700000 type=Te
position= 0.25000000 0.76297000 0.74700000 type=Te
position= 0.50000000 0.50000000 0.00000000 type=Ga
position= 0.50000000 0.00000000 0.49800000 type=Ga
position= 0.00000000 0.00000000 0.00000000 type=Cu
position= 0.00000000 0.50000000 0.49800000 type=Cu
position= 0.75000000 0.25000000 0.24900000 type=Es1
position= 0.25000000 0.75000000 0.24900000 type=Es1
position= 0.75000000 0.75000000 0.74700000 type=Es1
position= 0.25000000 0.25000000 0.74700000 type=Es1
position= 0.00000000 0.00000000 0.49800000 type=Es2
position= 0.50000000 0.50000000 0.49800000 type=Es2
position= 0.00000000 0.50000000 0.00000000 type=Es2
position= 0.50000000 0.00000000 0.00000000 type=Es2


結局CuGaTe2の結晶構造ファイルは、以下のようになります。

STRUC   ALAT=11.5388
PLAT=-1/2 1/2 0.99600
1/2 -1/2 0.99600
1/2 1/2 -0.99600
SITE ATOM=Te POS=0.23703000 0.25000000 0.24900000
ATOM=Te POS=0.76297000 0.75000000 0.24900000
ATOM=Te POS=0.75000000 0.23703000 0.74700000
ATOM=Te POS=0.25000000 0.76297000 0.74700000
ATOM=Ga POS=0.50000000 0.50000000 0.00000000
ATOM=Ga POS=0.50000000 0.00000000 0.49800000
ATOM=Cu POS=0.00000000 0.00000000 0.00000000
ATOM=Cu POS=0.00000000 0.50000000 0.49800000


結果


LDAを用いた結果、ワンショットGW近似を用いた結果、および、前回計算したAkaiKKRのGGAの結果を示します。

CuGaTe2-GW.png
Fig.1: QSGWによるCuGaTe2の状態密度

CuGaTe2-LDA.png
Fig.2: LDAによるCuGaTe2の状態密度

CuGaTe2DOS.png
Fig.3: GGA(AkaiKKR)によるCuGaTe2の状態密度


まず、大まかに見れば状態密度の形状はどれもほとんど同じであることが分かります。ただし、バンドギャップができるはずの部分では、LDAやAkaiKKR(GGA)の結果は、GW近似の結果と比較してギャップエネルギーを過小評価していることが分かります。

ecaljのログファイルには、バンドギャップの値が出力されます。前回のAkaiKKRの値とあわせて比較するとTable.1のようになりました。

ただし、LDA/GGAはバンドギャップを過小評価するという以上のことを言おうと思うと、微妙な議論になるかもしれません。
結晶構造と格子定数を合わせてあるとはいえ、細かい設定は異なります。AkaiKKRではpbe(GGA)を、ecaljのLDAはvwnを使っている点。AkaiKKRはフルポテンシャルでない点、などです。実際、空孔を入れたり原子球近似(ASA)を使ったりすることでギャップの値が改善されるということは、バンドギャップの値がポテンシャルの形状に影響を受けるはずだということです(参考: AkaiKKRでダイヤモンド型構造半導体)。

LDA/GGAのバンドギャップ過小評価は、第一原理計算パッケージ全般についてまわる問題だけあって、私には結論が出せませんが、仮にバンドギャップが正しい値を示さなくても、バンド構造や全エネルギーに関して議論を行うことは、何らかの意味があるかもしれません。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR ecalj LDA GGA GW近似 QSGW 半導体 バンドギャップ 

AkaiKKRとecaljでCuGaTe2 その1

第一原理計算パッケージには、それぞれ特徴があり、計算したい物質によって適切に使い分ける必要に迫られることがあります。AkaiKKR(machikaneyama)は不規則系に適しており、ecaljは半導体のバンドギャップを求めるのに適しています。

例えば、不規則を含む半導体の計算をAkaiKKRで行いたいと考えたとき、不規則を含まない端成分の計算をecaljの結果と比較しておくことは有用です。今回はCuGaTe2を対象として、AkaiKKRで状態密度の計算をおこないました。

CuGaTe2DOS.png
Fig.1: CuGaTe2の状態密度



AkaiKKRとecaljの長所


AkaiKKR(machikaneyama)は、コヒーレントポテンシャル近似(CPA)を導入することによって、合金などの不規則性を扱うことが可能であるという特徴があります。
またecaljはGW近似を用いて、半導体のバンドギャップの見積もりを局所密度近似(LDA)から改善できる長所があります。

他にもさまざまな第一原理計算パッケージが、それぞれ特有の長所を持っています。このため、しばしば複数のコードでの計算結果を比較するということが起こります。

今回と次回では、AkaiKKRの掲示板に投稿された CuGaTe2 のバンドギャップをこれら二つのコードで計算し、バンドギャップと状態密度の比較を行います。今回はAkaiKKRでの計算です。

計算手法


入力ファイルはCannot reproduce the bandgap of CuGaTe2に投稿されているものとほとんど同じですが、少しだけ変更してあります。一つ目の変更点は、スピン軌道相互作用を(計算が重いので)はずした事。二つ目はewidthを小さくしたことです。

c--------------------CuGaTe2---------------------------------
go data/cugate2
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
bct 11.5388 1.992 1 90 90 90
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 0.7 sra pbe nmag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 4 500 0.015
c------------------------------------------------------------
c ntyp
5
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Cu 1 0 0.0 2 29 100
Ga 1 0 0.0 2 31 100
Te 1 0 0.0 2 52 100
Es1 1 0 0.0 0 0 100
Es2 1 0 0.0 0 0 100
c------------------------------------------------------------
c natm
16
c------------------------------------------------------------
c atmicx atmtyp
0.23703x 1/4y 1/8z Te
0.76297x 3/4y 1/8z Te
3/4x 0.23703y 3/8z Te
1/4x 0.76297y 3/8z Te
1/2x 1/2y 0.0z Ga
1/2x 0.0y 1/4z Ga
0.0x 0.0y 0.0z Cu
0.0x 1/2y 1/4z Cu
c
0.75x 1/4y 1/8z Es1
0.25x 3/4y 1/8z Es1
3/4x 0.75y 3/8z Es1
1/4x 0.25y 3/8z Es1
c
0.0x 0.0y 0.25z Es2
1/2x 1/2y 0.25z Es2
0.0x 1/2y 0.0z Es2
1/2x 0.0y 0.0z Es2
c------------------------------------------------------------


結果


Fig.1に状態密度を示します。
AkaiKKRでの状態密度やバンド構造(ブロッホスペクトル関数)のエネルギー分解能は source/specx.f の msex で指定することが可能で、デフォルトでは msex=201 となっています。したがって、状態密度を計算するために ewidth = 0.8 Ry とした場合の分解能は 4 mRy 程度になります。その結果、状態密度の図だけを見ると、バンドギャップが存在するか否かが微妙です。

AkaiKKRでバンドギャップの測り方では、バンドギャップを決める場合、状態密度から値を読むよりも、バンド構造から見るほうが良さそうであると書きました。CuGaTe2は、伝導帯の上端(CBM)と価電子帯の下端(VBM)が共にΓ点に存在する直接遷移型の半導体であるとの事なので、その付近のバンド構造をプロットしたのがFig.2です。

CuGaTe2band.png
Fig.2: Γ点周辺のCuGaTe2のバンド構造


GaAsの場合と異なり、CBMにフェルミ準位(というか計算上のエネルギー基準点)が張り付いてしまっていますが、電子の数を足し上げるときの数値計算上の誤差と思うので、いまは気にしないことにします。

ローレンツ関数へのフィッティングは、あまりきれいにいかなかったので、目視で読むと、バンドギャップの大きさはおよそ 30 mRy 程度でしょうか。換算すると 0.4 eV 程度となるので、Cannot reproduce the bandgap of CuGaTe2に書かれている通り 1 eV 程度存在するはずのバンドギャップから見ると過小評価です。

AkaiKKRに限らず密度汎関数理論(DFT)に局所密度近似(LDA)や一般化勾配近似(GGA)を組み合わせた第一原理計算パッケージは、バンドギャップを過小評価してしまう問題が広く知られています。
ecaljで利用できるGW近似は、この問題に対する回答のひとつです。AkaiKKRとecaljでCuGaTe2 その2では、ecaljを用いてCuGaTe2の状態密度とバンドギャップを計算します。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR CPA ecalj 半導体 バンドギャップ バンド構造 分散関係 GW近似 

AkaiKKRで凝集エネルギー(未完)

凝集エネルギーは(ゼロ点振動のエネルギーなどを無視すれば)、凝集状態の全エネルギーと孤立原子の全エネルギーの差で表すことができます。仮に、孤立原子のエネルギーを、結晶の格子定数を大きくしていった極限として見積もることが可能なら、通常の第一原理計算パッケージでも凝集エネルギーを計算できることになります。

今回はAkaiKKR(machikaneyama)を用いてダイヤモンド構造のシリコンの格子定数を大きくしていって、全エネルギーの変化を調べる事により、シリコンの凝集エネルギーを見積もることを試みました。しかしながら、格子定数を大きくするに従って収束させるのが難しくなり、凝集エネルギーとしてリーズナブルな値を得ることはできませんでした。

CohesiveSi.png

Fig.1: ダイヤモンド型構造シリコンの全エネルギーの格子定数依存性。



凝集エネルギー


全エネルギーって何だよ?では、第一原理計算パッケージの出力する全エネルギーと固体物理の熱力学で言うところの凝集エネルギーの関係について書きました。
話の結論としては、全エネルギーは、エネルギーの原点に物理的な意味がないので、凝集エネルギーとは直接比較することができないが、よく似た概念であるということでした。

そこで、凝集エネルギーの原点である孤立原子の全エネルギーを、第一原理計算パッケージで計算しておけば、凝集状態の全エネルギーとの差から、凝集エネルギーを計算できるのではないでしょうか?

CohesiveOsaka2002.png

固体中の原子を引き離していったら…(PDF)で計算されているシリコンの格子定数と全エネルギー


このような試みはosaka2002のテクニカルノートにも見つけることができます(参考: 固体中の原子を引き離していったら…(PDF))。そこで今回は、AkaiKKR(machikaneyama)を用いてダイヤモンド構造のシリコンの格子定数を大きくしていって、全エネルギーの変化を調べました。

計算方法


これまで通り、入力ファイルのテンプレートからシェルスクリプトを利用して連続的に入力ファイルを作成していきます。交換・相関ポテンシャルにはmjwを使いました。これはMoruzzi Janak and Williamsに原子のエネルギーと凝集エネルギーの値が載っている事を意識してのことでしたが、この中にシリコンの値はありませんでした。Osaka2002のテクニカルノートにはスピン分極に関する言及もあるので、スピン分極ありの計算(mag)としました。相対論の効果は考慮せず nrl を指定しています。

c------------------------------------------------------------
go data/CohesiveSi_ABOHR
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
fcc ABOHR , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 1.0 nrl mjw mag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 4 1000 0.01
c------------------------------------------------------------
c ntyp
2
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Si 1 1 0.0 2
14 100
Vc 1 1 0.0 2
0 100
c------------------------------------------------------------
c natm
4
c------------------------------------------------------------
c atmicx(in the unit of a) atmtyp
0 0 0 Si
1/4 1/4 1/4 Si
1/2 1/2 1/2 Vc
3/4 3/4 3/4 Vc
c------------------------------------------------------------


#!/bin/csh -f

## *** プロジェクト名 ***
set PROJECT="CohesiveSi"
## *** 格子定数 ***
set ABOHR_LIST=( 10.26 11.97 13.68 15.39 17.10 18.81 20.52 22.23 23.94 )

## *** 繰り返し計算 ****
foreach ABOHR ( ${ABOHR_LIST} )
## ポテンシャルが存在しなければコピー
if ( ! -e data/${PROJECT}_${ABOHR} ) then
if ( -e data/${PROJECT} ) then
cp data/${PROJECT} data/${PROJECT}_${ABOHR}
endif
endif
## 入力ファイルの作成
sed 's/'ABOHR'/'${ABOHR}'/g' template/${PROJECT}_template.in > in/${PROJECT}_${ABOHR}.in
## 第一原理計算
specx < in/${PROJECT}_${ABOHR}.in > out/${PROJECT}_${ABOHR}.out
## 次のポテンシャルをコピー
cp data/${PROJECT}_${ABOHR} data/${PROJECT}
end


格子定数を大きくするほどバンドの分散が小さくなるので、ewidthも格子定数を大きくするのに従って小さくするほうが収束しやすくなるのだと思いますが、すべての計算で同じ値としました。

結果


Fig.2-3に格子定数 a0 = 10.26 Bohr のときのバンド構造と a = 15.39 Bohr (a = 1.5 a0, a0/a = 0.67)のときのバンド構造を示します。格子定数が大きくなることによってバンドの分散が小さくなっていることが分かります。

CohesiveSi1026.png
Fig.2: ダイヤモンド構造シリコンのバンド構造(a0 = 10.26 Bohr)

CohesiveSi1539.png
Fig.3: ダイヤモンド構造シリコンのバンド構造(a = 15.39 Bohr)


Fig.1に a0 = 10.26 Bohr としたときの格子定数の変化に対する全エネルギーをプロットしたものを示します。格子定数を大きくするほど全エネルギーが大きくなりましたが a = 20.52 Bohr (a0/a = 0.5) よりも大きい辺りから収束させるのが難しくなっていきました。これはbzqlty=4としましたが、bzqlty=9としても同様です。

計算できなくなる直前の値が、ほぼ孤立原子と同じエネルギーになっていると考えると凝集エネルギーは Ecoh = 0.47 Ry (6.37 eV) 程度となります。実験値は Ecoh = 4.63 eV との事なので(スピン分極を考慮した計算であるにもかかわらず)過大評価をしています。

議論(反省)


実を言うと凝集エネルギーに関する議論は、AkaiKKRの掲示板に投稿があります(How can I calculate cohesive energy? (in Japanese))。その内容によると、孤立原子のエネルギーは計算できるようにする予定が、少なくとも当時は、あったようですが、現在のバージョンにも実装はされていないようです。

前述の通り、格子定数を大きくするとバンドの分散が小さくなるのでewidthを小さくするほうがよさそうな気がします。
その一方で、同じ基準でエネルギーを比較するためにはewidthを固定する必要があるようにも思います。

格子定数 a = 10.26 Bohr では ewidth=1.0 Ry 程度が必要です。通常のgo計算では ewidth の範囲を source/specx.f の mse0 で指定されたエネルギー分解能で計算するので、この値を大きくすることにより、大きな ewidth でも高分解能で起算できるかもしれません。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR 

デバイ模型

デバイ模型やデバイ温度といった言葉は、固体物理学の分野でたびたび目にします。
近似自体は、非常に単純なものなのですが、あらゆる物性に顔を出す反面、初歩的な教科書には比熱のことしか書いてなかったりするので、ときどき良く分からなくなったりします。
デバイ模型は、そもそも、フォノンの分散関係を近似的に表すための模型であり、そこから(格子比熱を含む)色々な物性を計算することができるようになります。
今回は、デバイ模型の仮定している近似についてまとめると共に、デバイ温度が現れる物性の計算式を列挙します。

001_201511080539137af.png

Fig.1: Quantum ESPRESSOで計算したダイヤモンドのフォノンの分散関係



フォノンの分散関係と物性値


結晶の格子振動(フォノン)にだけに依存する物性値は、フォノンの分散関係が分かればすべて計算できるということになっています。デバイ模型(Debye model)やアインシュタイン模型(Einstein model)は、単純に格子比熱を計算するためだけのモデルだと見做される事もありますが、本来はこの分散関係を近似的にあらわす模型です。したがって、デバイ模型を使えば格子比熱以外の物性も計算することができます。

フォノンの分散関係は、以下の式で表され、一般には複雑な形状をしています。

\begin{equation}
\omega = \omega(\vec{q})
\end{equation}

ここで$\omega$は角振動数で$\vec{q}$は波数ベクトルです。Fig.1に示したのは、Quantum ESPRESSOで計算したダイヤモンドのフォノンの分散関係です。音響フォノンと光学フォノンの両方を持っていること、異方性を持っていることなどが分かります。

デバイ模型の2つの仮定


そこで、デバイ模型では「分散関係をシンプルな形で表す」「異方性をなくす」という二つの方針で簡単な式に置き換えます。

仮定(1) 角振動数$\omega$は波数$q$に比例し、その比例定数を音速$v$とする
\begin{equation}
\omega = v q
\end{equation}

仮定(2) ブリュアンゾーン(BZ)を同じ体積の球で置き換える
\begin{equation}
q_D = \left( \frac{6 \pi^2 N}{\Omega} \right)^{\frac{1}{3}}
\end{equation}
ここで$q_D$はデバイ波数と呼ばれ、体積Ωとその中の原子の数Nのみで表すことができます。

この波数$q$は、$0 \leqq q \leqq q_D$の範囲に限られます。この$q_D$は、体積Ωと原子数Nだけで決まるので、デバイ模型を特徴付けるパラメータは、音速$v$のみということになります。しかしながら、音速$v$はこの用途にはあまり使われず、代わりにデバイ温度$\Theta_D$がよく使われます。音速$v$とデバイ温度$\Theta_D$は、いったんエネルギーの次元で考えた後、以下のように換算できます。
\begin{equation}
k_B \Theta_D = \hbar \omega_D = \hbar v q_D \\
\therefore v = \frac{k_B \Theta_D}{\hbar q_D}
\end{equation}

以降では、デバイ模型のパラメータであるデバイ温度 ΘD が出てくる物性値について列挙します。

熱力学


ヘルムホルツの自由エネルギー(ヘルムホルツエネルギー)は、以下のような式で、エネルギー E の項とエントロピー S の項から表されます。
\begin{equation}
F=E-TS
\end{equation}
エネルギーにせよエントロピーにせよ、色々な原因による項(磁性、配置のエントロピー、電子比熱によるエントロピーなど)がありますが、格子系のエネルギーとエントロピーをデバイ模型を用いて表します。(参考: AkaiKKRで金属の熱物性)

格子系のエネルギー
\begin{equation}
E_D = 3 N k_B T \left( \frac{T}{\Theta_D} \right)^3 \int_0^{\Theta_D / T} \frac{z^3}{\exp(z)-1} \mathrm{d}z + E_0
\end{equation}

ゼロ点エネルギー
\begin{equation}
E_0 = \frac{9}{8} N k_B \Theta_D
\end{equation}

格子系のエントロピー
\begin{equation}
S_D = 3 N k_B \left[ \frac{4}{3} \left( \frac{T}{\Theta_D} \right)^3 \int_0^{\Theta_D / T} \frac{z^3}{\exp(z) - 1} \mathrm{d}z \\ - \ln \left\{ 1 - \exp \left(- \frac{\Theta_D}{T} \right) \right\} \right]
\end{equation}

定積比熱は格子系のエネルギーを、温度で微分したものです。
\begin{equation}
C_V = \frac{\partial E_D}{\partial T} = 9 N k_B \left( \frac{T}{\Theta_D} \right)^3 \int_0^{\Theta_D / T} \frac{z^4 \exp(z)}{(\exp(z)-1)^2} \mathrm{d} z
\end{equation}

格子振動の大きさ


デバイワラー因子
\begin{equation}
M_D = \frac{6 h^2 T}{M k_B \Theta_D^2} \left( \frac{T}{\Theta_D} \int_0^{\Theta_D / T} \frac{z}{\exp(z)-1} \mathrm{d} z + \frac{\Theta_D}{4T} \right) \left( \frac{\sin \theta}{\lambda} \right)^2
\end{equation}

平均二乗変位
\begin{equation}
u_D^2 = \frac{3 \hbar T}{M k_B \Theta_D^2} \left( \frac{T}{\Theta_D} \int_0^{\Theta_D/T} \frac{z}{\exp(z)-1}\mathrm{d}z + \frac{\Theta_D}{4T} \right)
\end{equation}

原子の振幅が一定以上大きくなると、融解が起こるのではないかという考えから、振幅と融点の関係式が提案されました。当然ながら、かなりラフな推定ですが、少なくとも定性的には、デバイ温度の高い物質の融点が高いという関係が見られるようです。(参考: ザイマン 固体物性論の基礎)

リンデマンの融解公式
\begin{equation}
T_m = \frac{x_m^2}{9 \hbar} M k_B \Theta_D r_s^2
\end{equation}

電気抵抗率


金属の電気伝導は、格子系だけでは決まらないので、電子と格子の相互作用に起因するパラメータも必要となります。

電気抵抗率
\begin{equation}
\rho_D = A \left( \frac{T}{\Theta_D} \right)^5 \int_0^{\Theta_D / T} \frac{z^5}{(\exp(z)-1)(1-\exp(-z))} \mathrm{d} z
\end{equation}

BCS理論による超伝導転移温度
\begin{equation}
T_c = \Theta_D \exp \left( - \frac{1}{N(0)V} \right)
\end{equation}

つまりデバイ模型/デバイ温度とは何なのか?


結局のところ、デバイ温度とかデバイ模型とか言われるものは、いったい何なのでしょうか?
結晶の比熱を計算するための近似がデバイ模型であり、デバイ温度はその中のパラメータである、というのはひとつの答えです。
実際、ねがてぃぶろぐでもScilabで数値積分: 固体の比熱のエントリで比熱の計算をしています。

001_20130522024102.png
Fig.2: 銅の比熱の温度依存性(ΘD = 343.5 K)。青破線がアインシュタインモデル、赤実線がデバイモデルによる計算。


しかしながら、上記のように、比熱以外の色々な物理量の計算式に顔を出します。また、直接数式の中に出てこない場合でも、例えば「デバイ温度よりも充分高温のときだけ成り立つ」みたいな文脈で現れることもあります。あるいは「ダイヤモンドは固いのでデバイ温度が高い」のような記述が見つかるかもしれません。

温度がデバイ温度よりも高いかどうかを議論しているときは、格子振動を古典力学的に扱ってよいのか、量子力学的な効果を考慮しなければならないのかを問題にしている場合が多いです。例えばFig.2の銅の比熱の温度依存性では、低い温度では比熱が急激に変化していますが、高温ではほとんど変化しなくなります。高温での比熱の温度依存性がなくなることは、古典力学的から導かれ、デュロン・プティの法則として知られています。その一方で、低温の急激な温度依存性は、量子力学を考えなければ説明できません。
また、物質の硬さを議論しているときには、直感的に言えば、デバイ温度を原子同士をつなぐ「ばね定数」のようなものだと考えています。

そんなわけで、最初に書いたとおり「デバイ模型とは、結晶のフォノンの分散関係に対する最も簡単な近似のひとつ」で「デバイ温度は、デバイ模型の中で、その物質の個性を表す最も重要なパラメータ」が答えとなるわけですが、その近似が単純であるにもかかわらず、色々な応用が考えられるので、出くわすたびに、それまでもっていた理解とは、ほんの少し違う側面が見えることがあります。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: デバイ模型 フォノン Scilab quantumESPRESSO 

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

LTspiceAkaiKKRScilabmachikaneyamaKKRPSoCOPアンプPICCPA強磁性モンテカルロ解析常微分方程式トランジスタode状態密度インターフェースDOSPDS5022ecaljスイッチング回路定電流半導体シェルスクリプトレベルシフト乱数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

最新コメント
リンク

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