PWscfでZnOの構造緩和

PWscf(Quantum ESPRESSO)を用いてウルツ鉱構造のZnOの高圧力下での構造最適化を行いました。

003_201807261443018d6.png

Fig.1: ウルツ鉱構造ZnOの格子定数の圧力変化



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


まず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軸方向につぶれたことが良く分かります。

001_20180726144258855.jpg
002_20180726144259aba.jpg

Fig.2-3: ZnOの結晶構造の変化


関連エントリ




参考URL




参考文献/使用機器




フィードバック



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

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


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


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

tag: PWscf QuantumESPRESSO 最適化 ZnO 

Scilabで最急降下法 その2

Scilabで最急降下法 その1では1変数の関数 f(x) に対して最急降下法を用いて最小値を求めました。今回は2変数の関数 f(x,y) について最小値を求めます。

001_20170509143114675.png
Fig.1: 最急降下法による最小値探査



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


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 最適化 最小値 最大値 

Scilabで最急降下法 その1

Scilabで何らかの関数 f(x) の最小値(あるいは最大値)を計算することを考えます。関数の値を計算するのが簡単な場合は x の定義域全体で f(x) を計算した後 minmax を使うという方法もあります。しかしながら f(x) の計算にそれなりの時間がかかる場合や f(x, y) といったように引数がたくさんある場合は効率的ではないと思います。

そこで今回は最急降下法のアルゴリズムを利用して f(x) の最小値を求めるということをやってみます。

001_20170507020049254.png

Fig.1: 最急降下法での最小値探索。上が関数f(x)の値、下が微分値f'(x)



最小値を求める関数


さて、実際に最小値を求める関数 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(0) を決める
  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




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 最適化 最小値 最大値 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAPICOPアンプecalj状態密度常微分方程式モンテカルロ解析odeトランジスタインターフェースDOSPDS5022スイッチング回路定電流半導体分散関係シェルスクリプト乱数レベルシフトHP6632A可変抵抗温度解析トランジスタ技術ブレッドボードR6452AI2C確率論セミナー数値積分反強磁性バンドギャップ熱設計非線形方程式ソルバ絶縁バンド構造偏微分方程式三端子レギュレータフォトカプラカオスマフィンティン半径ISO-I2CGW近似LM358A/DコンバータシュミットトリガLEDUSB数値微分サーボアナログスイッチ補間発振回路カレントミラー直流動作点解析TL43174HC4053PC817C単振り子FFTVESTA開発環境bzqlty電子負荷量子力学基本並進ベクトルパラメトリック解析標準ロジックチョッパアンプBSchLDAアセンブラブラべ格子2ちゃんねるイジング模型PWscf状態方程式仮想結晶近似キュリー温度Quantum_ESPRESSO熱伝導VCAスイッチト・キャパシタewidth最適化QSGWTLP621GGASMPMaxima失敗談位相図六方最密充填構造繰り返しスピン軌道相互作用相対論ランダムウォークFETgfortranコバルトスレーターポーリング曲線ラプラス方程式抵抗cygwin不規則合金格子比熱熱力学マントル条件分岐MCU井戸型ポテンシャルダイヤモンドQNAPUPS固有値問題シュレディンガー方程式自動計測ガイガー管詰め回路OpenMPTLP521ハーフメタルLM555ubuntufsolveブラウン運動平均場近似NE555ZnOTLP552QuantumESPRESSOxcrysdenCIF最小値最大値awkフェルミ面テスタ第一原理計算Ubuntu差し込みグラフFXA-7020ZR三角波過渡解析Writer509データロガースーパーセル起電力CK1026AACircuitMAS830LフィルタMBEP-10PGAトランスナイキスト線図ノコギリ波負帰還安定性EAGLEOPA2277PIC16F785CapSenseLMC6622SC1815入出力固定スピンモーメントFSMTeX結晶磁気異方性全エネルギーc/a合金multiplotgnuplot非線型方程式ソルバL10構造正規分布等高線ジバニャン方程式初期値interp1fcc面心立方構造ウィグナーザイツ胞半金属デバイ模型磁気モーメント電荷密度重積分SIC不純物問題ゼーベック係数cif2cellPWgui擬ポテンシャル二相共存ウルツ鉱構造edeltquantumESPRESSOフォノンリジッドバンド模型スワップ領域BaO岩塩構造ルチル構造ヒストグラム確率論マテリアルデザインフラクタルマンデルブロ集合キーボードRealforceクーロン散乱三次元疎行列縮退化学反応関数フィッティング最小二乗法Excel直流解析PCTS-110TS-112日本語パラメータ・モデル等価回路モデル文字列状態図陰解法熱拡散方程式HiLAPW両対数グラフCrank-Nicolson法連立一次方程式specx.fifort境界条件片対数グラフグラフの分割円周率ヒストグラム不規則局所モーメントGimpシンボル軸ラベル凡例線種トラックボール

最新コメント
リンク

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