ねがてぃぶろぐ

cif2cell srtio3.cif -p vasp --vasp-cartesian --vasp-format=5
vasp2ctrl POSCAR
mv ctrls.POSCAR.vasp2ctrl ctrls.srtio3

ecaljの入力ファイル

ecaljには、結晶構造ファイル ctrls.* からコントロールファイル ctrl.* を半自動的に生成するスクリプト ctrlgenM1.py が付属しています。従って、コントロールファイル ctrl.* がecaljの実質的な入力ファイルでありながら、ユーザーは結晶構造ファイル ctrls.* さえ作ればよいので簡単です。(参考: ecaljの実行手順(LDA計算))

結晶構造ファイルの中身はシンプルなので、人間が手動で作成してもたいした手間ではありませんが、ケアレスミスの可能性もあるので、自動化できるに越したことはありません。そこで今回は、結晶構造を記述するファイル形式として有名なCIFから結晶構造ファイル ctrls.* を半自動的に生成する方法について書きます。
対象にした結晶は、正方晶(tetragonal)の SrTiO₃ です。この結晶構造ファイルは Crystallography Open Databaseからダウンロードすることが出来ます。このCIFをVESTAで描画したのがFig.1です。

cif2cellのインストール

SourceForgeからcif2cellをダウンロードし、ホームディレクトリに置きます。私がダウンロードした時点での最新版は cif2cell-1.2.10.tar.gz でした。これを展開します。

tar xzvf cif2cell-1.2.10.tar.gz

展開されたディレクトリへ移動し、インストールを行います。

cd cif2cell-1.2.10/
sudo python setup.py install

更に .bashrc に以下を追記してパスを設定しておきます。

export PATH=$PATH:$HOME/cif2cell-1.2.10

追記したら .bashrc を再読み込みさせます。

source .bashrc

cifから結晶構造ファイル ctrls.* の作成

まずcif形式のファイルを用意します。今回はCrystallography Open Databaseからダウンロードします。

wget http://www.crystallography.net/cod/cif/9/00/28/9002806.cif

次に cif2cell を用いてcif形式からvasp形式のファイルを作成します。

cif2cell 9002806.cif -p vasp --vasp-cartesian --vasp-format=5

このとき私の環境では、以下のような警告が出ますが、無視して進めます。

***Warning : Site occupancies not found, assuming all occupancies = 1.

すると POSCAR というファイルが出来ています。このファイルからecaljの結晶構造ファイルを作ります。

vasp2ctrl POSCAR

結晶構造ファイル ctrls.POSCAR.vasp2ctrl が出来ているはずです。 vasp2ctrl という名前ですが、作成されるのはコントロールファイル ctrl.* では無く、結晶構造ファイル ctrls.* です。

ecaljで強磁性鉄のスピン分極計算やecaljで反強磁性NiOなどの場合は ctrls.POSCAR.vasp2ctrl を基にして自生の設定のための編集が必要になりますが、今回は単純に名前だけ変更します。

mv ctrls.POSCAR.vasp2ctrl ctrls.srtio3

あとは通常通りにLDA計算を行うだけです。

ctrlgenM1.py srtio3
cp ctrlgenM1.ctrl.srtio3 ctrl.srtio3
lmfa srtio3
mpirun -np 2 lmf-MPIK srtio3

関連エントリ

• ecaljの実行手順(LDA計算)
• ecaljで強磁性鉄のスピン分極計算
• ecaljで反強磁性NiO

参考URL

• ecalj
• SourceForge
• Crystallography Open Database
• VESTA

フィードバック

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



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

結晶構造の最適化

多くの第一原理計算パッケージには、結晶構造の最適化機能があります。これは多くの場合、原子にかかる力を計算して、その方向に原子を動かすことによって実現されています。しかしながらAkaiKKR(machikaneyama)では原子にかかる力の計算が出来ません。AkaiKKRでも、全エネルギーが最小になるようにすることで結晶構造の最適化が出来るはずですが、手動ですべてのパラメータの探索を行うのは大変です。そこで最急降下法を使って全エネルギーが最小になる構造を探して見ます。
Scilabで最急降下法 その1とScilabで最急降下法 その2では、簡単な関数f(x,y)に対して、最小値を探すことで最急降下法のアルゴリズムを確認しました。今回は f(x,y)の代わりに第一原理計算を行います。今回はテスト計算として六方最密充填構造(hcp)のコバルトの計算をします。具体的には以下の3つのステップを行います。

1. マフィンティン半径の決定
2. エネルギーマップの作成
3. 最急降下法計算

マフィンティン半径の決定

AkaiKKRで全エネルギーを比較するときは、色々な計算条件が変化しないようにする必要があります。特にマフィンティン半径の設定には注意が必要です。AkaiKKRの入力ファイルでは、マフィンティン半径は格子定数 a で規格化されます。私の理解では「AkaiKKRでマフィンティン半径を固定する」というのは「原子の位置や計算セルのサイズを変更したときに、計算セルの体積に対する各マフィンティン球の体積の比が変化しないようにする」ということです。(これは他の計算手法、例えばAPW法などとは違うことに注意が必要です。)この条件を満たす範囲でマフィンティン半径を最大になるようにするのが最良であると理解しています。(AkaiKKRでBain機構 その1とその2も参照。)

そこでまず、全エネルギー最小を探索する格子定数の範囲を決めます。次に、その範囲でrmt=1としたときの実際のマフィンティン球の体積比が一番小さくなってしまう条件を探します。この条件でのマフィンティン半径を全ての計算に用いれば、マフィンティン球が重なることなく、かつ、マフィンティン半径を固定することが出来ます。
決められたマフィンティン半径は、格子体積の1/3乗で規格化しておくのが後々の事を考えると便利なはずです。具体的には下記のようなシェルスクリプトを作成し、計算しました。1.60 ≦ c/a ≦ 1.70 の範囲では r_MT/V^1/3 = 0.439518 ぐらいのようです。
このスクリプトを走らせる際には、第一原理計算を収束させる必要はないので bzqlty=0とmaxitr=0としておきます。

#!/bin/csh -f
#setenv GFORTRAN_UNBUFFERED_ALL y

## *** パラメーター範囲 ***
set COA_LIST=( 1.60 1.61 1.62 1.63 1.64 1.65 1.66 1.67 1.68 1.69 1.70 )
set OMEGA_LIST=( 160 158 156 154 152 150 148 146 144 142 140 )

# set COA_LIST=( 1.63 )
# set OMEGA_LIST=( 150 )

## *** マフィンティン半径の計算 ***
set OMEGA=${OMEGA_LIST[1]}
set COA0=`echo $COA_LIST[1]`
set A0=`echo "e((1/3)*l(2*${OMEGA}/(sqrt(3)*${COA0})))" | bc -l`
if (`echo "${COA0} > 2*sqrt(2)/sqrt(3)" | bc -l` == 1) then
  set RMTB0=`echo "${A0}/2" | bc -l`
else
  set RMTB0=`echo "${A0}*sqrt(1/3+(${COA0}^2)/4)/2" | bc -l`
endif
set COA1=`echo $COA_LIST[$#COA_LIST]`
set A1=`echo "e((1/3)*l(2*${OMEGA}/(sqrt(3)*${COA1})))" | bc -l`
if (`echo "${COA1} > 2*sqrt(2)/sqrt(3)" | bc -l` == 1) then
  set RMTB1=`echo "${A1}/2" | bc -l`
else
  set RMTB1=`echo "${A1}*sqrt(1/3+(${COA1}^2)/4)/2" | bc -l`
endif
if (`echo "${RMTB0} < ${RMTB1}" | bc -l` == 1) then
  set RMTB=`echo $RMTB0`
else
  set RMTB=`echo $RMTB1`
endif
set RMTOMEGA=`echo "${RMTB}/e((1/3)*l(${OMEGA}))" | bc -l | sed -e 's/^\./0./g'`
echo ${RMTOMEGA}

なおhcp構造の場合はAkaiKKRでコバルトのc/a その2の方法でマフィンティン半径を決めることが出来るのですが、今回は将来的により複雑な構造をやることも考えてこのようにしました。

エネルギーマップの作成

必要ないはずの手順ですが、今回は最急降下法のテストでもあるので、考える体積・軸比の範囲の全ての全エネルギーを計算しておきます。この結果がFig.1のカラーコンターです。

Cシェルスクリプトには関数を実装することが出来ませんが、下請けのシェルスクリプトに変数を渡すことによって、それっぽい事は出来ます。
今回は という下請けのスクリプトを作成しました。これは、体積・軸比・マフィンティン半径の3つを引数として受け取り、第一原理計算を行って全エネルギーを返り値にします。

#!/bin/csh -f
##setenv GFORTRAN_UNBUFFERED_ALL y

## *** プロジェクト名 ***
set PROJECT="hcpCo"

## *** 実行ファイル ***
set EXEC="~/kkr/cpa2002v009c/specx"

## *** 標準入力から値を読み取る ***
## 格子体積
set OMEGA=$1
## c/a
set COA=$2
## RMT
set RMTOMEGA=$3

## ファイル名
set INFILE="in/${PROJECT}_go_${OMEGA}_${COA}.in"
set OUTFILE="out/${PROJECT}_go_${OMEGA}_${COA}.out"
set POTFILE="data/${PROJECT}_${OMEGA}_${COA}"
set POTBACK="data/${PROJECT}_${OMEGA}"

## 格子定数 a の計算
set ABOHR=`echo "scale=7; e((1/3)*l(2*${OMEGA}/(sqrt(3)*${COA})))" | bc -l | sed -e 's/^\./0./g'`
## マフィンティン半径の計算
set RMTA=`echo "scale=7; ${RMTOMEGA}*e((1/3)*l(${OMEGA}))/${ABOHR}" | bc -l | sed -e 's/^\./0./g'`

## テンプレートから入力ファイルを作成
sed 's/'OMEGA'/'${OMEGA}'/g' template/${PROJECT}_go_Template.in | sed 's/'ABOHR'/'${ABOHR}'/g' | sed 's/'COA'/'${COA}'/g' | sed 's/'RMTA'/'${RMTA}'/g' > ${INFILE}

## ポテンシャルファイルのコピー
if ( ! -e ${POTFILE} ) then
  if ( -e ${POTBACK} ) then
    cp ${POTBACK} ${POTFILE}
  endif
endif

## 計算回数の初期化
set num=0
## 最大計算回数
set nummax=20
## 第一原理計算
${EXEC} < ${INFILE} > ${OUTFILE}
while ( ( ! { grep -q "err= -6." ${OUTFILE} } ) && ( $num < $nummax ) )
  ${EXEC} < ${INFILE} > ${OUTFILE}
  @ num++
end

## ポテンシャルのバックアップ
cp ${POTFILE} ${POTBACK}

set ENE=`grep "total energy" ${OUTFILE} | sed -e s/total//g  -e s/energy=//g`
echo ${ENE}

enemap.sh は下請けスクリプト dogo.sh を利用して 1.60 ≦ c/a ≦ 1.70, 140 ≦ V ≦ 160 の範囲で全エネルギーのマップを作ります。dogo.sh は最急降下法のシェルスクリプトでも利用します。

最急降下法

エネルギーのマップから最安定な格子定数が V = 150 Bohr³, c/a = 163 付近にあることが予想できます。とりあえず初期値を V = 156 Bohr³, c/a = 1.68 として計算してみます。その結果がFig.1上に白で示された経路です。最急降下法では、その関数の値の微分の方向に向かって次の入力パラメータを探します。今回の計算では、4回程度で格子定数の最適化が出来ている事が確認できます。

最急降下法のCシェルスクリプト steepest.sh は以下の通りです。

#!/bin/csh -f
#setenv GFORTRAN_UNBUFFERED_ALL y

## *** マフィンティン半径 ***
set RMTOMEGA=0.43951815598150528923

## *** 係数 ***
set KEISUU_OMEGA="10000.0"
set KEISUU_COA="0.5"

## *** 初期値 ***
set OMEGA=$1
set COA=$2

## *** 微分のステップ ***
set dOMEGA="1.0"
set dCOA="0.01"
set OMEGA_PLUS=`echo "${OMEGA}+${dOMEGA}" | bc -l`
set OMEGA_MINUS=`echo "${OMEGA}-${dOMEGA}" | bc -l`
set COA_PLUS=`echo "${COA}+${dCOA}" | bc -l`
set COA_MINUS=`echo "${COA}-${dCOA}" | bc -l`

## *** 第一原理計算 ***
set ENE=`./dogo.sh ${OMEGA} ${COA} ${RMTOMEGA}`
echo "Center energy:" ${ENE} "(Ry)"
set ENE_OMEGA_PLUS=`./dogo.sh ${OMEGA_PLUS} ${COA} ${RMTOMEGA}`
echo "Omega plus:   " ${ENE_OMEGA_PLUS} "(Ry)"
set ENE_OMEGA_MINUS=`./dogo.sh ${OMEGA_MINUS} ${COA} ${RMTOMEGA}`
echo "Omega minus:  " ${ENE_OMEGA_MINUS} "(Ry)"
set ENE_COA_PLUS=`./dogo.sh ${OMEGA} ${COA_PLUS} ${RMTOMEGA}`
echo "c/a plus:     " ${ENE_COA_PLUS} "(Ry)"
set ENE_COA_MINUS=`./dogo.sh ${OMEGA} ${COA_MINUS} ${RMTOMEGA}`
echo "c/a minus:    " ${ENE_COA_MINUS} "(Ry)"

## *** 数値微分(中心差分) ***
set dENEdOMEGA=`echo "(${ENE_OMEGA_PLUS}+(-1*${ENE_OMEGA_MINUS}))/(2*${dOMEGA})" | bc -l`
echo ${dENEdOMEGA}
set dENEdCOA=`echo "(${ENE_COA_PLUS}+(-1*${ENE_COA_MINUS}))/(2*${dCOA})" | bc -l`
echo ${dENEdCOA}

echo ${OMEGA} ${COA} ${ENE} >> analysis/steepest.txt

set OMEGA=`echo "scale=7; (${OMEGA}-${KEISUU_OMEGA}*${dENEdOMEGA})/1.0" | bc -l`
set COA=`echo "scale=7; (${COA}-${KEISUU_COA}*${dENEdCOA})/1.0" | bc -l`
echo "Next:"
echo "./steepest.sh" ${OMEGA} ${COA}

関連エントリ

• AkaiKKRの使い方
• Scilabで最急降下法 その1
• Scilabで最急降下法 その2
• AkaiKKRでBain機構 その1
• AkaiKKRでBain機構 その2
• AkaiKKRでコバルトのc/a その2

参考URL

• AkaiKKR(machikaneyama)

付録

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

• hcpCo_go_Template_in.txt
• dogo_sh.txt
• rmtomega_sh.txt
• enemap_sh.txt
• steepest_sh.txt

参考文献/使用機器

フィードバック

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



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

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

参考文献/使用機器

フィードバック

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



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