ねがてぃぶろぐ

入力ファイルのテンプレート

まずPWscfの入力作成補助で紹介したxlt2pw.pyを利用して、ZnOの入力ファイルの雛形を作ります。入力の元になるxtlファイルは、CIFからVESTAを通じて出力します。CIFはDFT計算用データベース MatNaviで紹介した無機材料データベース(AtomWork)などからダウンロードします。

xtl2pw.py ZnO.xtl vc-relax 1
cp ZnO.in ZnO.relax.in

擬ポテンシャルとカットオフエネルギー

PWscfの擬ポテンシャルを参考に入力ファイルを編集します。
擬ポテンシャルの種類には PAW ポテンシャルを選びました。

ATOMIC_SPECIES
  Zn    65.38   Zn.pbe-dnl-kjpaw_psl.1.0.0.UPF
   O   15.999   O.pbe-n-kjpaw_psl.0.1.UPF

カットオフエネルギーは、擬ポテンシャルファイルに書かれている推奨値を参考にします。波動関数のカットオフの推奨値は酸素の方が大きく、電荷密度の推奨値は亜鉛の方が大きいことに注意が必要です。それぞれ大きいほうの1.5倍ぐらいを選びました。

    ecutwfc = 75.0 ,
    ecutrho = 420.0 ,

圧力

圧力下の構造最適化を行うこともできます。その場合は press に圧力を設定します。単位は kbar です(1 GPa = 10 kbar)。

原子位置の制約

xtl2pw.py の最後の入力パラメータは、原子位置の制約の有無に関連しています。
0 を入力すると全ての原子について位置の最適化を行います。
1 を入力すると結晶学的に動かなさそうな原子を固定して計算します。

下記は xtl2pw.py の最後のパラメータに 1 をしていた場合の入力ファイルです。原子の座標の後ろに更に3つの数字が続いていますが、これが原子を動かすか否かのフラグになっていて 0 なら固定、1 なら動かせることを意味しています。今回の例では、酸素原子のz方向だけ原子を最適化するようになっていることが分かります。

ATOMIC_POSITIONS crystal
Zn     0.333333   0.666667   0.000000 0 0 0
Zn     0.666667   0.333333   0.500000 0 0 0
O      0.333333   0.666667   0.389010 0 0 1
O      0.666667   0.333333   0.889010 0 0 1

入力ファイル

結局、入力ファイルは以下のようにしました。

&control
    calculation='vc-relax' ,
    restart_mode='from_scratch' ,
    prefix='ZnO' ,
    outdir = './ZnO/' ,
    wfcdir = './ZnO/' ,
    pseudo_dir = './' ,
    disk_io='default' ,
    forc_conv_thr= 0.001 ,
    verbosity = 'default' ,
    nstep = 100 ,
 /
 &system
    ibrav= 4 ,
    celldm(1) = 6.14709011 ,
    celldm(3) = 1.60158627686 ,
    nat =  4 ,
    ntyp = 2 ,
    ecutwfc = 75.0 ,
    ecutrho = 420.0 ,
 /
 &electrons
    electron_maxstep = 100 ,
    mixing_beta = 0.7 ,
!   use smaller conv_thr for better results ,
    conv_thr = 1.0d-12 ,
 /
 &ions
    ion_dynamics='bfgs' ,
 /
 &CELL
    cell_dynamics = 'bfgs' ,
    press = 0.001,
    press_conv_thr = 0.05 ,
!   cell_dofree = 'xyz' ,
 /
ATOMIC_SPECIES
  Zn    65.38   Zn.pbe-dnl-kjpaw_psl.1.0.0.UPF
   O   15.999   O.pbe-n-kjpaw_psl.0.1.UPF
ATOMIC_POSITIONS crystal
Zn     0.333333   0.666667   0.000000 0 0 0
Zn     0.666667   0.333333   0.500000 0 0 0
O      0.333333   0.666667   0.389010 0 0 1
O      0.666667   0.333333   0.889010 0 0 1
K_POINTS automatic
 6 6 4 0 0 0

結果

圧力を変化させながら格子定数がどのように変化するかを調べた結果が Fig.1 です。40 GPaまでは正常に圧縮されていっていますが 50 GPa で異常が見られます。
pwout2xtl.pyをもちいて出力したxtlファイルをVESTAで描画してみると、酸素原子のz位置が大きく動いて、結晶全体がc軸方向につぶれたことが良く分かります。

関連エントリ

• PWscfの入力作成補助
• DFT計算用データベース MatNavi
• PWscfの擬ポテンシャル

参考URL

• xlt2pw.py
• 無機材料データベース(AtomWork)
• VESTA

参考文献/使用機器

フィードバック

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



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

具体的には二次元のガウス関数に-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

参考文献/使用機器

フィードバック

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



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