ねがてぃぶろぐ

点欠陥の形成エネルギー

結晶中に点欠陥を作るために必要な欠陥の形成エネルギーは、スーパーセル法を用いた第一原理計算から計算することができます。
具体的には、まず、スーパーセルで完全結晶の全エネルギーを計算します(E₀)。次にスーパーセルから1つの原子を取り除いた系の計算を行い、全エネルギー(E₁)を求めます。スーパーセルの原子数が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個のアルミニウム原子を置いた完全結晶の全エネルギー(E₀)とスーパーセルからひとつの原子を点欠陥に置き換えたスーパーセルの全エネルギー(E₁)の計算を行いました。得られた全エネルギーは、以下のようになりました。

E₀ = -15482.440171667 (Ry)
E₁ = -14998.335584355 (Ry)

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

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

関連エントリ

• AkaiKKRでスーパーセル その1
• AkaiKKRとUbuntu 12.04 のスワップ領域

参考URL

• AkaiKKR (machikaneyama)
• Wang et al. (2004) (PDF)

付録

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

• SuperAl_in.txt
• SuperAlVc_in.txt

参考文献/使用機器

フィードバック

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



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

不純物の計算

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

他方、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でニッケル・鉄・コバルトなどの純鉄の状態密度と一致することが確認できます。



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



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

関連エントリ

• AkaiKKRの角運動量(方位量子数)のメモ
• AkaiKKRでニッケル・鉄・コバルト

参考URL

• AkaiKKR
• how to run a system of over 30 atoms?
• スーパーセル(気象)

付録

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

• FeSuper_in.txt
• FeLiSuper_in.txt
• FeLi_in.txt

参考文献/使用機器

フィードバック

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



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