ねがてぃぶろぐ

同じ元素に対しても擬ポテンシャルファイルには色々な種類があるようです。

交換相関汎関数

AkaiKKRやecaljの場合と同様に、交換相関汎関数にどのタイプを使うかを選択できます。計算結果は、どの汎関数を利用するかによって変わってきます(参考: AkaiKKRでLDAとGGA その1, その2)。
どれを使うかのが一番よいのかは、ケースバイケースなのでしょう。とりあえず PBE から試すのが最近の流行なのでしょうか。

擬ポテンシャルファイルの種類

擬ポテンシャルファイルには、以下の3種類の分類があるようです。

• ノルム保存型
• ウルトラソフト
• PAW(Projector Augmented-Wave)

それぞれに特徴があるので、目的によってどのタイプを使うのかを決めます。ただし擬ポテンシャルの種類によって選択できる元素の種類が限定されます。擬ポテンシャルファイルを自作することも可能らしいですが、そうでなければ計算したい化学組成についてどのタイプの擬ポテンシャルが使えるのかはよく確認しておく必要がありそうです。元素ごとに異なるタイプの擬ポテンシャルファイルを混在させることも可能ですが、出来れば避けたほうが良いだろうと思います。

PAWポテンシャルは、おそらく3つの中で最も新しい方法で、計算の高速さと正確さの両方を同時にあげることが出来るらしいです。ただし新しい手法なので(?)選べる元素の種類は少なそう(?)です。

ノルム保存型は3つの中ではおそらく最も古く、そのため選べる元素の種類も多そうです。

ウルトラソフト擬ポテンシャルは、ノルム保存型よりもカットオフエネルギーを小さく出来る、したがって計算速度を上げることが出来るメリットがあるようです。

カットオフエネルギー

カットオフエネルギーは、計算精度に関するパラメーターです。大きくするほど計算結果が正確になりますが、計算時間が掛かるようになります。したがって、必要充分な値を設定する必要があります。この値は擬ポテンシャルファイルの先頭部分に書いてあることが多いようです。以下は Siの擬ポテンシャルファイル の冒頭部分です。

Generated using "atomic" code by A. Dal Corso  v.5.0.2 svn rev. 9415
Author: ADC
Generation date: 11Sep2012
Pseudopotential type: PAW
Element: Si
Functional:  SLA  PW   PBX  PBC

Suggested minimum cutoff for wavefunctions:  38. Ry
Suggested minimum cutoff for charge density: 151. Ry
The Pseudo was generated with a Scalar-Relativistic Calculation
L component and cutoff radius for Local Potential:  2   2.0000

波動関数のカットオフエネルギーの最小値は 38 Ry で、電荷密度のカットオフエネルギーの最小値は 151 Ry であると書かれています。利用する全ての擬ポテンシャルファイルの中で、一番大きな値の 1.5 倍程度が目安な様です。

関連エントリ

• PWscfのインストール
• AkaiKKRでLDAとGGA その1, その2

参考URL

• PSLIBRARY

参考文献/使用機器

フィードバック

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



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

wget http://qe-forge.org/gf/download/frsrelease/247/1132/qe-6.2.1.tar.gz
tar xzvf qe-6.2.1.tar.gz
cd qe-6.2.1/
./configure
make all

Quantum ESPRESSO(PWscf)のウエブページを見ると、最新版はGitLabあるいはqe-forgeからソースコードをダウンロードするように書かれています。今回は wget を使って qe-6.2.tar.gz をダウンロードします。

wget http://qe-forge.org/gf/download/frsrelease/247/1132/qe-6.2.1.tar.gz

次にダウンロードされたファイルを展開します。

tar xzvf qe-6.2.1.tar.gz

ファイルが正しく展開された後は、アーカイブは削除してしまって構いません。どこかに保存しておいたほうがいいかもしれませんが。

rm qe-6.2.1.tar.gz

新しく作られたフォルダへ移動し、コンフィグを実行します。

cd qe-6.2.1/
./configure

configure: success と表示されれば成功です。Intel Fortran Compiler または gfortran がインストールされていれば、基本的には上手く行くはずですが･･･

いよいよ make します。単純に make とだけ端末でタイプすると make のオプションが表示されます。2. quantum-ESPRESSOのインストールでは make pw pp ph で充分と書かれているのですが、私はmake allとしてしまいました。

make all

途中でいくつか warning がでますが bin/ ディレクトリを見てみるとちゃんとバイナリファイルが出来ているようです。

関連エントリ

• AkaiKKRのインストール
• ecaljのインストール(Ubuntu + gfortran)
• Ubuntu 12.04 64bitでHiLAPW

参考URL

• Quantum ESPRESSO(PWscf)
• quantum ESPRESSO tutorial

参考文献/使用機器

フィードバック

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



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

ゼーベック係数

ゼーベック係数(熱電能)は、以下の式で表されます。

\begin{equation}
S = \frac{1}{eT}\frac{K_1}{K_0}
\end{equation}
ここでK_nは
\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------------------------------------------------------------

次に計算された状態密度から、電子の数密度n_eを計算します。
\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}

電子の数密度n_eは温度に関わらず一定なので、絶対零度 T = 0 (K) のときの電子数密度 n_e0 をまず計算します。フェルミ分布関数はこのとき ε < ε_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)が大きい)ほど大きなゼーベック係数を持つことが分かります。

関連エントリ

• AkaiKKRの使い方
• Scilabで金属の化学ポテンシャル

参考URL

• AkaiKKR(machikaneyama)

付録

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

• Kelvin_sce.txt
• Kelvin-Pt-calc.txt
• Kelvin-Pt-lit.txt

参考文献/使用機器

フィードバック

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



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

ウルツ型構造の結晶構造はFig.2のような六方晶です。このCIFはCrystallography Open Databaseからダウンロードすることが出来ます。



AkaiKKRの入力ファイルは、以下のようになりました。原子位置(atomicx)はVESTAでAkaiKKRのための基本並進ベクトルの方法で出力したものを利用しました。

c----------------------GaN-----------------------------------
     go   data/GaN
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
     hex      6.01 , 1.6245 ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
    0.001     1.8       sra      mjw      nmag     2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update      4       200    0.03
c------------------------------------------------------------
c    ntyp
      2
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
    Ga       1       1      0.0     2
                                          31    100
    N        1       1      0.0     2
                                          7     100
c------------------------------------------------------------
c   natm
     4
c------------------------------------------------------------
c   atmicx                        atmtyp
    0.333333343a   0.666666687b   0.000000000c  Ga
    0.666666627a   0.333333313b   0.500000000c  Ga
    0.333333343a   0.666666687b   0.384999990c  N
    0.666666627a   0.333333313b   0.884999990c  N
c------------------------------------------------------------

関連エントリ

• VESTAでAkaiKKRのための基本並進ベクトル

参考URL

• AkaiKKR(machikaneyama)
• Crystallography Open Database

付録

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

• GaN_in.txt

参考文献/使用機器

フィードバック

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



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

入力ファイル

磁鉄鉱の結晶構造は以下のように図示できます。



一見すると複雑な結晶構造ですが、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------------------------------------------------------------

バンド構造

以下にブロッホスペクトル関数(バンド構造)を示します。バンド構造から見てもハーフメタルになっていることが分かります。

関連エントリ

• AkaiKKRでハーフメタル

参考URL

• AkaiKKR(machikaneyama)
• Yanase and Siratori (1984, J. Phys. Soc. Jpn)

付録

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

• fccMagnetite_in.txt
  

参考文献/使用機器

フィードバック

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



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