ねがてぃぶろぐ

AkaiKKRの部分状態密度

AkaiKKR(machikaneyama)ではgoモードでセルフコンシステント計算を行った後、出力されたポテンシャルファイルからdosモードでspecxを実行することによって状態密度を計算することが出来ます。
状態密度は全状態密度とコンポーネントごとに分かれた部分状態密度の両方が出力されます。デフォルトでは部分状態密度は s, p, d, f, ...軌道に分かれて出力されます。第一原理計算入門 AkaiKKRのページではt_2gやe_gに分解した部分状態密度を計算する方法が書かれています。
今回はこれにしたがって source/spmain.f を編集して鉄の部分状態密度を計算しました。

部分状態密度計算用の実行ファイル

第一原理計算入門 AkaiKKRのページに書かれている通りですが、手順を書きます。具体的には source/spmain.f を編集して再び make するだけです。

まず source/spmain.f と実行ファイルの specx のバックアップを取ります。

cp source/spmain.f source/spmain.f.back
mv specx specx.back

次に source/spmain.f の該当部分を第一原理計算入門 AkaiKKRのページに書かれている通りに編集します。(以下のコマンドでは emacs で編集していますが vi でも gedit でもお好きなエディタでどうぞ)

emacs -nw source/spmain.f

まず下記に該当する部分を探します。

c     --- print partial and the total DOS if required.
      if(ids .eq. 1 .or. ids .eq. 2 .or. ids .eq. 3) then
      estep=dble(e(2,is))-dble(e(1,is))
      do 69 i=1,ncmpx
      write(*,'(//1x,a,i2,a,i2)')'DOS of component',i
      do 69 k=1,kk
      xmd(i,k,1,is)=-dimag(wkc(2,i,k))/pi
cc    &    (-dimag(wkc(4,i,k))/pi)+(-dimag(wkc(2,i,k))/pi)
      xmd(i,k,2,is)=-dimag(wkc(4,i,k))/pi
cc    &    (-dimag(wkc(4,i,k))/pi)-(-dimag(wkc(2,i,k))/pi)
c  69 write(*,'(1x,f7.4,3x,9f8.4)') dble(e(k,is))-ef(is)
c    &      ,( -dimag(wkc(l,i,k))/pi,l=1,mxl**2)
      do 160 l=1,mxlcmp(i)
c     do 160 l=1,2
      do 160 m=1,2*(l-1)
  160 wkc(l**2,i,k)=wkc(l**2,i,k)+wkc(l**2-m,i,k)
c     wkc(5,i,k)=wkc(5,i,k)+wkc(6,i,k)+wkc(8,i,k)
c     wkc(7,i,k)=wkc(7,i,k)+wkc(9,i,k)
   69 write(*,'(1x,f7.4,3x,4f10.4)') dble(e(k,is))-ef(is)
     &      ,(-dimag(wkc(l**2,i,k))/pi,l=1,mxlcmp(i))
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'total DOS',(dble(e(k,is))-estep/2d0-ef(is)
     &      ,dimag((detl(k,is)-detl(k-1,is))/(e(k,is)-e(k-1,is)))
     &        ,k=2,kk)
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'integrated DOS',(dble(e(k,is))-ef(is)
     &       ,dimag(detl(k,is)),k=1,kk)

上記の部分をごっそりと削除して、下記の記述に置き換えます。

c     --- print partial and the total DOS if required.
      if(ids .eq. 1 .or. ids .eq. 2 .or. ids .eq. 3) then
      estep=dble(e(2,is))-dble(e(1,is))
      write(*,'(///a)')
     & '(PDOS DATA: Ry, s, px, pz, py, dxy, dyz, dz^2, dxz, dx^2-y^2)'
      do 69 i=1,ncmpx
      write(*,'(//1x,a,i2,a,i2)')'DOS of component',i
      do 69 k=1,kk
      xmd(i,k,1,is)=-dimag(wkc(2,i,k))/pi
cc    &    (-dimag(wkc(4,i,k))/pi)+(-dimag(wkc(2,i,k))/pi)
      xmd(i,k,2,is)=-dimag(wkc(4,i,k))/pi
cc    &    (-dimag(wkc(4,i,k))/pi)-(-dimag(wkc(2,i,k))/pi)
   69 write(*,'(1x,f7.4,3x,9f8.4)') dble(e(k,is))-ef(is)
     &      ,( -dimag(wkc(l,i,k))/pi,l=1,mxl**2)
c      do 160 l=1,mxlcmp(i)
c     do 160 l=1,2
c      do 160 m=1,2*(l-1)
c  160 wkc(l**2,i,k)=wkc(l**2,i,k)+wkc(l**2-m,i,k)
c     wkc(5,i,k)=wkc(5,i,k)+wkc(6,i,k)+wkc(8,i,k)
c     wkc(7,i,k)=wkc(7,i,k)+wkc(9,i,k)
c   69 write(*,'(1x,f7.4,3x,4f10.4)') dble(e(k,is))-ef(is)
c     &      ,(-dimag(wkc(l**2,i,k))/pi,l=1,mxlcmp(i))
      if(is .eq. 1) then
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'total_up TDOS_up',(dble(e(k,is))-estep/2d0-ef(is)
     &      ,dimag((detl(k,is)-detl(k-1,is))/(e(k,is)-e(k-1,is)))
     &        ,k=2,kk)
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'integrated_up IDOS_up',(dble(e(k,is))-ef(is)
     &       ,dimag(detl(k,is)),k=1,kk)
      end if
      if(is .eq. 2) then
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'total_dn TDOS_dn',(dble(e(k,is))-estep/2d0-ef(is)
     &      ,dimag((detl(k,is)-detl(k-1,is))/(e(k,is)-e(k-1,is)))
     &        ,k=2,kk)
      write(*,'(//1x,a/(1x,f12.7,f13.5))')
     &      'integrated_dn IDOS_dn',(dble(e(k,is))-ef(is)
     &       ,dimag(detl(k,is)),k=1,kk)
      end if

編集したら make を実行します。

make

出来たファイルはpdos用に別名保存するようにしました。ついでにバックアップしておいた通常の specx を復帰させておきます。

mv specx specx.pdos
mv specx.back specx

また source/spmain.f も元に戻しておいたほうがいいでしょう。代わりに編集したバージョンのバックアップを取っておきます。

mv source/spmain.f source/spmain.f.pdos
mv source/spmain.f.back source/spmain.f

鉄の計算

入力ファイルは通常の鉄の計算用のものと基本的には変わりませんが、今回の specx.pdos はd電子までしか計算できないので l=2 とします。また部分状態密度は全状態密度よりもギザギザになりやすいので bzqlty はかなり高めにしました。

c----------------------Fe------------------------------------
    dos   data/fe
c------------------------------------------------------------
c   brvtyp     a        c/a   b/a   alpha   beta   gamma
    bcc      5.27  ,     ,      ,      ,       ,      ,
c------------------------------------------------------------
c   edelt    ewidth    reltyp   sdftyp   magtyp   record
    0.001     1.0       nrl      mjw      mag      2nd
c------------------------------------------------------------
c   outtyp    bzqlty   maxitr   pmix
    update     18        100   0.035
c------------------------------------------------------------
c    ntyp
     1
c------------------------------------------------------------
c   type    ncmp    rmt    field   mxl  anclr   conc
    Fe       1       1      0.0     2
                                          26    100
c------------------------------------------------------------
c   natm
     1
c------------------------------------------------------------
c   atmicx(in the unit of a)     atmtyp
     0         0          0        Fe
c------------------------------------------------------------

結果として得られる部分状態密度は s, p, p_y, p_z,
d_xy, d_yz, d_z2, d_xz, d_x2-z2の順に出力されます。
今回は e_g=d_z2+d_x2-y2 と t_2g = d_xy+d_yz+d_xz についてプロットしました。

関連エントリ

• AkaiKKRの使い方

参考URL

• AkaiKKR(machikaneyama)
• 第一原理計算入門 AkaiKKR

付録

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

• FePDOS_plt.txt
• FePDOS_up.txt
• FePDOS_dn.txt

参考文献/使用機器

フィードバック

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



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

具体的には二次元のガウス関数に-1を掛けた関数を f(x,y) とします。
\begin{equation}
f(x,y)=\frac{1}{2 \pi \sigma_1 \sigma_2 \sqrt{1-\rho^2}}\exp\left\{-\frac{1}{2(1-\rho^2)}\\
\times\left(\frac{(x-\mu_1)^2}{\sigma_1^2}-2\rho\frac{(x-\mu_1)(y-\mu_2)}{\sigma_1\sigma_2} +\frac{(y-\mu_2)^2}{\sigma_2^2}\right) \right\}
\end{equation}
ここで
\begin{equation}
\rho=\frac{\sigma_{12}}{\sigma_1\sigma_2}
\end{equation}

最急降下法の値の更新は二次元の場合は以下のように行います。
\begin{equation}
\begin{pmatrix}
x^{(k+1)}\\
y^{(k+1)}
\end{pmatrix}
=
\begin{pmatrix}
x^{(k)}\\
y^{(k)}
\end{pmatrix}
-\alpha
\begin{pmatrix}
\frac{\partial f(x^{(k)}, y^{(k)})}{\partial x}\\
\frac{\partial f(x^{(k)}, y^{(k)})}{\partial y}
\end{pmatrix}
\end{equation}
停止条件は ∂f/∂x < ε かつ ∂f/∂y < ε とすればよいと思います。

clear;

// *** 二次元ガウス分布 ***
mu1 = 3; mu2 = 1;
mu = [mu1; mu2];
sigma1 = sqrt(10); sigma2 = sqrt(10); sigma12 = 5;
SIGMA = [sigma1^2, sigma12; sigma12, sigma2^2];
rho = sigma12 / (sigma1 * sigma2);
function z = func(x, y)
  z = -1 ./ (2 * %pi * sigma1 * sigma2 * sqrt(1 - rho)) ..
  .* exp(-1 / (2 .* (1 - rho^2)) ..
  .* ((x - mu1) .^ 2 ./ (sigma1 ^ 2) - 2 .* rho .* (x - mu1) .* (y - mu2) ./ (sigma1 * sigma2) + ((y - mu2) .^ 2) ./ (sigma2 ^ 2)))
endfunction


// *** 数値微分 ***
h = 1E-3;
function dzx = dfx(x, y)
  dzx = (func((x+h), y) - func((x-h), y)) ./ (2 * h)
endfunction

function dzy = dfy(x, y)
  dzy = (func(x, (y+h)) - func(x, (y-h))) ./ (2 * h)
endfunction

// *** グラフのプロット ***
x = linspace(-10,10);
y = linspace(-10,10);
[X,Y] = ndgrid(x,y);
Z = func(X,Y);
// 色の設定
//xset("colormap",jetcolormap(64))
// xset("colormap",graycolormap(64))
// Sgrayplot(x,y,Z)
xset("fpf"," ");		// 等高線に値を表示しない
contour2d(x,y,Z,10);
xmin = min(X); xmax = max(X); ymin = min(Y); ymax = max(Y);
zoom_rect([xmin,ymin,xmax,ymax]);




// *** 最小値の計算 ***
// 停止条件
err = 1E-5;
a = 200;
// 初期値
x = -4; y = 4;
z = func(x, y);
dx = dfx(x, y);
dy = dfy(x, y);
plot(x, y, "xk");
// 最小値の計算
while ((abs(dx) > err) | (abs(dy) > err))
  x = x - a * dx;
  y = y - a * dy;
  dx = dfx(x, y);
  dy = dfy(x, y);
  plot(x, y, ".r");
end
// 計算結果
x, y
plot(x, y, "xk");

関連エントリ

• Scilabで最急降下法 その1

参考URL

• 半塲 滋 非線形計画法
• Wikipedia 最急降下法

付録

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

• steepest2d_sce.txt

参考文献/使用機器

フィードバック

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



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

最小値を求める関数

さて、実際に最小値を求める関数 f(x) ですが、今回は単純にガウス関数に -1 を掛けたものにします。
\begin{equation}
f(x) = - \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{(x-\mu)^2}{2\sigma^2}\right)
\end{equation}
当然ながら、f(x)が最小になるのは x = μ のときです。

最急降下法

高校の数学で習ったとおり f(x) が最大値や最小値(や極値)をとるときその微分は f'(x) = df(x)/dx = 0 となります。最急降下法は、関数の微分を計算しその傾きが大きいほうへ f'(x)=0 となる x を探すアルゴリズムです。具体的には以下の手続きを繰り返します。

1. x の初期値 x⁽⁰⁾ を決める
2. f'(x) < ε なら終了
3. x^(k+1) = x^(k) - αf'(x^(k))
4. 2.に戻る

実際には α や ε を上手に決めておく必要があります。αは勾配の方向にどの程度進むかを決めるパラメータ(下記Scilabスクリプトではa)で、εは計算の終了条件を決めるパラメータ(下記Scilabスクリプトではerr)です。

Scilabスクリプト

clear;

// *** 一次元ガウス分布 ***
function y = func(x)
  mu = 3;
  sigma = 1;
  y = -1 / sqrt(2*%pi*sigma^2) * exp(-1 * ((x - mu) .^ 2) ./ (2*sigma^2))
//  y = cos(x)
endfunction

// *** 数値微分 ***
function y = dfunc(x)
  h = 1E-4;
  y = (func(x+h) - func(x-h)) ./ (2*h)
endfunction

// *** グラフのプロット ***
X = linspace(0,6);
Y = func(X);
dY = dfunc(X);
subplot(2,1,1);
plot(X, Y);
subplot(2,1,2);
plot(X, dY);

// *** 最小値の計算 ***
// 停止条件
err = 1E-3;
a = 0.5;
// 初期値
x = 1;
y = func(x);
dx = dfunc(x);
subplot(2,1,1);
plot(x, y, ".r");
subplot(2,1,2);
plot(x, dx, ".r");

// *** 最小値の計算 ***
while abs(dx) > err
  x = x - a * dx;
  dx = dfunc(x);
  y = func(x);
  subplot(2,1,1);
  plot(x, y, ".r");
  subplot(2,1,2);
  plot(x, dx, ".r");
end

// *** 計算結果 ***
x
subplot(2,1,1);
plot(x, y, "xk");

参考URL

• 半塲 滋 非線形計画法
• Wikipedia 最急降下法

付録

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

• steepest1d_sce.txt

参考文献/使用機器

フィードバック

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



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

積分範囲が長方形領域の場合

積分範囲が長方形の領域の場合、すなわち以下のような式で表すことができる場合は、簡単に数値積分できます。

\begin{equation}
\int_{x_0}^{x_1}\int_{y_0}^{y_1}f(x,y)\mathrm{d}x\mathrm{d}y
\end{equation}

Scilabの int2d では長方形領域を2つの三角形のパッチワークとして与えます。
積分範囲を int2d に渡すために行列XとYを用意します。それぞれ2つの三角形の頂点のx座標とy座標を与えます。

\begin{equation}
X =
\begin{pmatrix}
x_{0} & x_{0} \\
x_{1} & x_{1} \\
x_{1} & x_{0} \\
\end{pmatrix},
Y =
\begin{pmatrix}
y_{0} & y_{0} \\
y_{0} & y_{1} \\
y_{1} & y_{1} \\
\end{pmatrix}
\end{equation}



実際に以下の積分を計算して見ます。

\begin{equation}
\int_{0}^{\pi}\int_{0}^{R}x^4 \sin(y)\mathrm{d}x\mathrm{d}y
\end{equation}

clear;

r = 1;

// *** 積分する関数の定義 ***
function z = func(x,y)
  z = (x .^ 4) * sin(y)
endfunction
// 積分範囲
x0 = 0; x1 = r;
y0 = 0; y1 = %pi;

// *** 二重積分 ***
X = [x0, x0;
     x1, x1;
     x1, x0];
Y = [y0, y0;
     y0, y1;
     y1, y1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = 2*(r^5)/5

数値化解と解析解が同じ値になることが確認できます。

積分範囲が三角形の組み合わせで表せる場合

積分範囲が長方形の場合は2つの三角形の組み合わせで表現されますが、より複雑な形状の場合も任意の個数の三角形の組み合わせで表現できるはずです。今回は逆に簡単になってしまいますが、1個の三角形で表現できる例を計算します。

\begin{equation}
\int \int_D xy^2 \mathrm{d}x\mathrm{y}
\end{equation}



積分領域が三角形ひとつ分なので、与える行列は3行1列になります。

\begin{equation}
X =
\begin{pmatrix}
x_{0} \\
x_{1} \\
x_{1} \\
\end{pmatrix},
Y =
\begin{pmatrix}
y_{1} \\
y_{0} \\
y_{1} \\
\end{pmatrix}
\end{equation}

この計算を行うScilabスクリプトは以下のようになります。

clear;

// *** 積分する関数の定義 ***
function z = func(x,y)
  z = x .* (y .^ 2)
endfunction

// *** 二重積分 ***
X = [0;
     1;
     1];
Y = [1;
     0;
     1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = 3/20

このスクリプトも数値解と解析解が同じに値になることが分かります。

同様にしてN個の三角形の組み合わせで表現される積分範囲の場合3行N列の行列で指定することができます。

更に複雑な積分領域の場合

どんなに複雑な積分領域の形状であっても三角形のパッチワークで表現できるはずですが、現実的には大変です。そこでOctaveの精義―フリーの高機能数値計算ツールを使いこなすで紹介されている方法を試してみましたが、現状うまく行っていません。上手く行っていませんがとりあえず方法だけは紹介します。
具体的にはScilabの論理演算で条件分岐の考え方を使って積分領域外では値がゼロになるように被積分関数の定義を行います。

\begin{equation}
\int\int_D -\frac{1}{(2x + y + 1)^2}\mathrm{d}x\mathrm{d}y
\end{equation}

clear;

// *** 積分する関数の定義 ***
function z = func(x,y)
  region = y >= x .^ 2
  z = - 1 ./ ((2 * x + y + 1) .^ 2) .* region
endfunction

// *** 二重積分 ***
X = [0;
     1;
     1];
Y = [0;
     0;
     1];
// 数値解
I = int2d(X, Y, func)
// 解析解
A = (1/3) * log(4) - 1/2

原理的にはこのスクリプトでよいはずですが、実際には正しく計算してくれません。Scilab 6.0ではエラーで停止します。Scilab 5.5.2ではそれっぽい値を返しますが、解析解の値とはかなりずれた値となっており、不正確です。

関連エントリ

• 数値積分タグ
• Scilabの論理演算で条件分岐

参考URL

• 高校数学の美しい物語

付録

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

• int2dsample1_sce.txt
• int2dsample2_sce.txt
• int2dsample3_sce.txt

参考文献/使用機器

フィードバック

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



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