AkaiKKRで点欠陥の形成エネルギー(未完)

Wang et al. (2004) (PDF)を参考に、AkaiKKR (machikaneyama)を用いて面心立方構造(fcc)のアルミニウムの点欠陥の生成エネルギーをスーパーセル法を用いて計算しました。
結果は、Wang et al. (2004) (PDF)の計算結果に対して、一桁程度の過大評価となってしまいました。


点欠陥の形成エネルギー


結晶中に点欠陥を作るために必要な欠陥の形成エネルギーは、スーパーセル法を用いた第一原理計算から計算することができます。
具体的には、まず、スーパーセルで完全結晶の全エネルギーを計算します(E0)。次にスーパーセルから1つの原子を取り除いた系の計算を行い、全エネルギー(E1)を求めます。スーパーセルの原子数がN個のとき、欠陥形成エネルギーは、以下の式から求めることができます(参考: Wang et al. (2004) (PDF))。
\begin{equation}
\Delta E_{f} = E_{1} - \frac{N-1}{N} E_{0}
\end{equation}

同様の計算は、スーパーセルではなくコヒーレントポテンシャル近似(CPA)でもできる可能性はあるのでしょうか?このような質問がAkaiKKR (machikaneyama)の掲示板に投稿されて・・・いましたが、現在は存在していないようです(#6678)。質問者の計算は、うまく行っておらず、形成エネルギーが一桁程度過大評価されているようです。

そこで今回は、AkaiKKR (machikaneyama)でもスーパーセル法を用いて点欠陥の形成エネルギーを計算してみます。

計算手法


Wang et al. (2004) (PDF)では、面心立方構造(fcc)のアルミニウムとニッケル、体心立方構造(bcc)のモリブデンとタンタルの点欠陥の形成エネルギーの計算が行われています。今回は、この中で面心立方構造のアルミニウムに関して計算を行います。スーパーセルのサイズは 2*2*2=32 原子としました。
Wang et al. (2004) (PDF)によると、交換相関汎関数はGGAよりもLDAのほうが良いであるとか、構造緩和はしないほうがむしろ良いだとか、色々議論があるようです。とりあえず今回は、LDA(mjw)で構造緩和もしない(というか大変なのでやりたくない)という方針で行きます。

c--------------------Al--------------------------------------
go data/SuperAlVc
c------------------------------------------------------------
c brvtyp a c/a b/a alpha beta gamma
sc 15.3 , , , , , ,
c------------------------------------------------------------
c edelt ewidth reltyp sdftyp magtyp record
0.001 1.0 sra mjw nmag 2nd
c------------------------------------------------------------
c outtyp bzqlty maxitr pmix
update 4 200 0.035
c------------------------------------------------------------
c ntyp
32
c------------------------------------------------------------
c type ncmp rmt field mxl anclr conc
Al1 1 1 0.0 2 0 100
Al2 1 1 0.0 2 13 100
Al3 1 1 0.0 2 13 100
Al4 1 1 0.0 2 13 100
Al5 1 1 0.0 2 13 100
Al6 1 1 0.0 2 13 100
Al7 1 1 0.0 2 13 100
Al8 1 1 0.0 2 13 100
Al9 1 1 0.0 2 13 100
Al10 1 1 0.0 2 13 100
Al11 1 1 0.0 2 13 100
Al12 1 1 0.0 2 13 100
Al13 1 1 0.0 2 13 100
Al14 1 1 0.0 2 13 100
Al15 1 1 0.0 2 13 100
Al16 1 1 0.0 2 13 100
Al17 1 1 0.0 2 13 100
Al18 1 1 0.0 2 13 100
Al19 1 1 0.0 2 13 100
Al20 1 1 0.0 2 13 100
Al21 1 1 0.0 2 13 100
Al22 1 1 0.0 2 13 100
Al23 1 1 0.0 2 13 100
Al24 1 1 0.0 2 13 100
Al25 1 1 0.0 2 13 100
Al26 1 1 0.0 2 13 100
Al27 1 1 0.0 2 13 100
Al28 1 1 0.0 2 13 100
Al29 1 1 0.0 2 13 100
Al30 1 1 0.0 2 13 100
Al31 1 1 0.0 2 13 100
Al32 1 1 0.0 2 13 100
c------------------------------------------------------------
c natm
32
c------------------------------------------------------------
c atmicx atmtyp
0 0 0 Al1
1/4 1/4 0 Al2
1/4 0 1/4 Al3
0 1/4 1/4 Al4

1/2 0 0 Al5
3/4 1/4 0 Al6
3/4 0 1/4 Al7
1/2 1/4 1/4 Al8

0 1/2 0 Al9
1/4 3/4 0 Al10
1/4 1/2 1/4 Al11
0 3/4 1/4 Al12

0 0 1/2 Al13
1/4 1/4 1/2 Al14
1/4 0 3/4 Al15
0 1/4 3/4 Al16

1/2 1/2 0 Al17
3/4 3/4 0 Al18
3/4 1/2 1/4 Al19
1/2 3/4 1/4 Al20

1/2 0 1/2 Al21
3/4 1/4 1/2 Al22
3/4 0 3/4 Al23
1/2 1/4 3/4 Al24

0 1/2 1/2 Al25
1/4 3/4 1/2 Al26
1/4 1/2 3/4 Al27
0 3/4 3/4 Al28

1/2 1/2 1/2 Al29
3/4 3/4 1/2 Al30
3/4 1/2 3/4 Al31
1/2 3/4 3/4 Al32
c------------------------------------------------------------


AkaiKKRでスーパーセル その1で書いたとおり、AkaiKKRでスーパーセルの計算を行うためには、それに適したパラメータをspecx.fに設定して再コンパイルする必要があります。計算するコンピュータのメモリが少ない場合、スワップ領域を使う必要があるかもしれません。今回はAkaiKKRとUbuntu 12.04 のスワップ領域で指定した下記のパラーメータをspecx.fに設定しました。

     & (natmmx=32, ncmpmx=32, msizmx=288, mxlmx=3, nk1x=500, nk3x=701,


結果と議論


面心立方構造のアルミニウムの点欠陥の形成エネルギーを計算するために、計算セルの中に32個のアルミニウム原子を置いた完全結晶の全エネルギー(E0)とスーパーセルからひとつの原子を点欠陥に置き換えたスーパーセルの全エネルギー(E1)の計算を行いました。得られた全エネルギーは、以下のようになりました。

E0 = -15482.440171667 (Ry)
E1 = -14998.335584355 (Ry)

よって点欠陥の生成エネルギーは
ΔEf = 0.2783319 (Ry) = 3.786906 (eV)
となりました。

この値はWang et al. (2004) (PDF)で報告されている 0.568 (eV) @ N=31, 0.511 (eV) @ N=108 とくらべて一桁程度の過大評価となってしまいました。したがって、AkaiKKRのCPA計算で点欠陥の生成エネルギーを過大評価してしまうのは、必ずしもCPAの問題ではないかもしれません。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR スーパーセル CPA 

AkaiKKRでスーパーセル その1

AkaiKKRではコヒーレントポテンシャル近似(CPA)を用いて不純物を含んだ系の第一原理計算を行うことが出来ます。しかしながら、これは(他のコードでしばしば行われる)スーパーセルの計算が出来ないことを意味しません。
今回は計算機マテリアルデザイン入門 に載っている例であるFe87.5Li12.5合金の計算をスーパーセル法で行いました。

Chaparral_Supercell_520.jpg

Fig.1: スーパーセル(気象)。なお今回のエントリとは関係ない。



不純物の計算


多くの第一原理計算パッケージでは、不純物や格子欠陥の計算を行う際に、単位格子を複数個並べてその中の少数の原子を他の組成や格子欠陥に置き換えると言うことを行います。このように単位格子をたくさん並べたものをスーパーセルと、スーパーセルを使う計算手法をスーパーセル法と呼びます。

他方、AkaiKKRではコヒーレントポテンシャル近似(CPA)を用いることで、単位格子を用いながら不純物や格子欠陥の計算ができます。
CPAでは不規則性の効果によりバンド構造がにじむ(電子が有限の寿命を持つ)効果がはっきりと確認できます。
これに対して、スーパーセル法では、不純物周りの構造緩和を取り入れることが出来ます。

AkaiKKRはCPAが使えるという点が特徴的であるのですが、当然ながら、スーパーセルの計算ができないと言うわけではありません。(実際に構造緩和までやろうとすると大変だとは思いますが。)

今回は計算機マテリアルデザイン入門 に載っている体心立方構造(bcc)の鉄を2×2×2のスーパーセルとして、その中の1原子をリチウムに置き換えた計算を行います。

入力ファイル


以下に8個の鉄原子のうち1個をリチウムに置換したスーパーセルの入力ファイルを示します。
この他にスーパーセル法のテストのために純鉄の計算もスーパーセルで行います。

 go    data/fe_li_super
bcc 10.85, , , , , ,
0.001 1.2 nrl mjw mag 2nd
update 10 200 0.024
3
Fe1 1 0 0 2 26 100
Fe2 1 0 0 2 26 100
Li 1 0 0 1 3 100
8
0.00000 0.00000 0.00000 Li
0.50000 0.00000 0.00000 Fe1
0.00000 0.50000 0.00000 Fe1
0.00000 0.00000 0.50000 Fe1
0.25000 0.25000 0.25000 Fe2
0.25000 0.25000 0.75000 Fe2
0.25000 0.75000 0.25000 Fe2
0.75000 0.25000 0.25000 Fe2


比較のためにCPAを用いた不規則合金の計算も行います。

 go    data/fe_li
bcc 5.43, , , , , ,
0.001 1.2 nrl mjw mag 2nd
update 10 200 0.024
1
FeLi 2 0 0 2 26 87.5
3 12.5
1
0.00000 0.00000 0.00000 FeLi


大きなセルでのspecx.fの設定


今回のスーパーセルは2×2×2なので、原子の数は高々8個なのですが、計算する原子の数が増えるとspecx.fを編集して再makeする事が必要になる可能性があります。

パラメータの設定にはhow to run a system of over 30 atoms?AkaiKKRの角運動量(方位量子数)のメモが参考になるはずです。
以下に、赤井先生の回答とそれを私が意訳したものを示します。

It is still possible to run for a system with over 30 atoms.
Take care about the following:

1) Decrease the number of k-points in order that the memory addressing does not overflow (in the case of 32bit addressing system, such as Pentium 4 with 32bit addressing). This can be done by decreasing "nk1x". The test run should be done bu "bzqlty=0" or "bzqlty=1" for which nk1x=1~20 would be enough.

2) Increase the number aoms that can be handled:
For that increase "natmmx", "ncmpmx", and "msizmx".

For example if you have 30 atoms and for each atom you would like to calculate up to d states (l=2), all atoms have different types, each site has a single component (i.e. not a random alloy) the following parameters are suitable:
natmmx=30, ncmpmx=30, msizmx=270, mxlmx=3

30原子以上の系に関しても計算を行うことが出来ます。以下の点に注意してください。

1) メモリのオーバーフローを避けるためk点の数を減らします。そのためには"nk1x"を小さくします。"bzqlty=0"や"bzqlty=1"のテスト計算ではnk1x=1~20で充分です。

2) 原子数を増やします。これには"natmmx", "ncmpmx" 及び "msizmx"を大きくします。

例えば30原子をd状態までの計算で、全ての原子が異なる組成で、各サイト(原子位置)には1種類の組成しか持たない(つまり不規則合金でない)場合は、下記のパラメータとなります。
natmmx=30, ncmpmx=30, msizmx=270, mxlmx=3


計算結果


まずテスト計算としての純鉄の状態密度をスーパーセル法で計算したものを示します。
AkaiKKRでニッケル・鉄・コバルトなどの純鉄の状態密度と一致することが確認できます。

FeSuper.png

Fig.2: スーパーセル法により計算した鉄の状態密度


次にスーパーセル法を用いたFeLi合金とコヒーレントポテンシャル近似(CPA)を用いて計算したFeLi合金の状態密度を示します。

FeLiCPA.png

Fig.3: スーパーセル法(赤)とコヒーレントポテンシャル近似(緑)によるFe87.5Li12.5合金の状態密度


スーパーセル法による結果とCPAによる結果は、状態密度のおおよその形状に関しては非常に良く似ています。しかしながらCPAでは不規則性の効果で状態密度の鋭さが減じているのに対して、スーパーセル法ではとがったままです。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR CPA スーパーセル 

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

最新コメント
リンク

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