体積変化による自由エネルギーの様子

phonopy,MedeA,Quasi-harmonic近似法を用いてAlが300Kのときの最安定な格子定数 を算出した．図5.2は縦軸に自由エネルギー[KJ/mol],横軸に格子定数を体積変化させた 倍率を示している．Calculation-dataはMedeAとphonopyの計算値であり，実線は計算 値に対してfittingを行ったものである．MedeA,Quasi-harmonic近似法は極小値を持つが

phonopyは極小値がなく値は温度が上がるにつれて自由エネルギーの値は低下している．

このことからphonopyで算出した自由エネルギーは最安定をとることが出来ないという ことが明らかとなった．最安定がとれなければ熱膨張係数を算出することが出来ず，これ 以上の計算は不可能である．従ってphonopyの算出する結果は正しく再現出来ているの は絶対値のみであり，熱膨張の影響は考慮できていないということが分かった．

図 5.2: 各手法における極小値の有無.

第 6 _{章 総括}

本研究ではphonopyとMedeAを用いてAl,NaClの有限温度における自由エネルギーを算 出し，お互いの計算結果を比較することでphonopyの計算精度を検証した．研究開始当初は phonopyとMedeAの結果が一致しなかったため,Moruzzi らが提案したQuasi-harmonic 近似法を用いてAlの自由エネルギーを算出することで，どちらの計算が正しいかを見 積もった. その結果,MedeAとQuasi-harmonic近似法の算出した値は一致し，phonopy の計算結果を再度検討することになった．計算を見直し，再検討した結果，phonopyと MedeA,Quasi-harmonic近似法の自由エネルギーの絶対値はほぼ一致した．phonopyと

MedeAの結果が違った原因としてphonopyは自由エネルギー値を１ユニットセル毎に算

出しているため，1原子あたりに算出しているMedeAと比較する際はphonopyの値を1 原子あたりのエネルギーに変更しなければいけないことが分かった．しかし，同一温度で 各格子定数ごとの自由エネルギーの様子を比較すると，MedeA,Quasi-harmonic近似法で は最小値を得られたが，phonopyでは値は軸比が上がる毎に自由エネルギーが低下し，極 小値は得られなかった．この結果から熱膨張係数などの熱膨張の影響を考慮した物性値の 計算に移行することは困難であり，phonopyはMedeA同等の計算結果を再現できないと いう事が明らかとなった．

謝辞

本研究を遂行していくにあたり，終始多大なる有益なご指導，また私の就職後の事につ いても親身になって考えてくださいました西谷滋人教授に深い感謝の意を表します．研究 を進めるにあたり，多くの援助，親身に相談にのって頂きました清原資之先輩をはじめ，

さまざまな知識の供給などの御協力を頂いた西谷研究室の皆様，そして私が無事に卒業で きるまで，支えてくださった両親に感謝の意を表します．最後になりましたが，この場を 借りて心から深く御礼申し上げます．

付録

phonopy の導入と使用方法

インストール方法

以下の手順で環境を構築し，phonopyのインストールを行う．VirtualBox,Vagrant,UbuntuOS のインストールはphonopyをインストールする際に使用するapt-getツールを導入するた めに行うので，既存でAPTがある場合は省略してもかまわない．

VirtualBox

VirtualBoxは無料で提供されている仮想化アプリケーションである．仮想化とはホス

トOS上にソフトウェアで仮想PCを作成して、その中で別のゲストOSを実行するア プリケーションである．本研究ではゲストOSとしてUbuntuを導入する際に使用する．

https://www.virtualbox.org/wiki/Downloads [1]にアクセスし，VirtualBox version for

OS X hostsを選択し，ソフトウェアをダウンロード，GUIを用いて指示される手順の通

りにインストールする．詳細はVirtualBoxをダウンロードする際にUserManualもダウ ンロードされるのでそちらを参照してほしい

Vagrant

Vagrantは仮想環境の構築から設定までを自動的に行うことができる．つまり，本研究

ではVirtualBoxとUbuntuのつなぎ目の役割を果たしている．

http://www.vagrantup.com/downloads.html [2]にアクセス，MAC OS Xを選択し，

Vagrantをダウンロードする．インストールはmacのターミナル上で以下の通りに行う．

[yanase-no-MacBook-Pro:~/.vagrant.d] yanase%

vagrant plugin install vagrant-vbguest

Vagrantの起動,入り方

virtualBoxを起動させ，Vagrantを立ち上げるコマンドを以下に記す．これにより，phonopy コマンドを用いることが出来る．

[yanase-no-MacBook?-Pro:~/Vagrant] yanase% vagrant up

Vagrantにsshでログインする．

yanase% vagrant ssh

UbuntuOS

UbuntuOSはコミュニティーにより開発されているLinux系列のオペレーティングシス

テムであり，無償で提供されている．無償にも関わらずセキュリティーアップデートが頻 繁に行われ，常に最新のパージョンを容易に手に入れることが出来る．

[yanase-no-MacBook-Pro:~/.vagrant.d] yanase%

vagrant box add ubuntu http://opscode-vm-bento.s3.amazonaws.com/vagrant/

virtualbox/opscode_ubuntu-13.10_chef-provisionerless.box

[yanase-no-MacBook-Pro:~/.vagrant.d] yanase% vagrant init ubuntu

phonopy

http://sourceforge.net/projects/phonopy [4] からphonopyをダウンロード，以下の

ようにvagrant上でインストールを行う．インストールはPythonのパスを通すことが必

須となっており，Pythonのライブラリはapt-getを使用する．つまりVagrant上での作業

1. Pythonの環境を整えるためにpythonのライブラリをインストール 2. apt-getが最新のバージョンか確認

3. ホームディレクトリにphonopyファイルを移動 4. phonopyのパッケージを解凍

5. pythonのPATHを通す．このパスはそれぞれ異なるため，自分のpythonのファイ ルを探してパスを通す

6. phonopyのファイルsetup.pyがあるディレクトリ上でインストールを行う．

1 sudo apt-get install python-dev python-numpy

 python-matplotlib python-tk python-lxml python-yaml 2 sudo apt-get update --fix-missing

3 mv /vagrant/phonopy-1.8.4.2-rc3.tar.gz . 4 tar xvfz phonopy-1.8.4.2-rc3.tar.gz 5 export PYTHONPATH=/usr/lib/python2.7/

6 sudo python setup.py install --home=.

上記のパスでpythonのパスが上手く通らないときは以下のようにPATHを変更する必 要がある．

＃PATHの確認 echo $PATH 

＃パスの編集

PATH=$PATH:/home/vagrant/phonopy-1.8.4.2-rc3/lib/python export PATH

setup.pyが上手くいかない場合は，インストールの実行を２回に分けて行う．

sudo python setup.py build sudo python setup.py install

phonopy Tutorial

1. Pre-process : モデルの構造方法

2. 力計算 : VASPで得られたデータから力定数を計算する方法

3. Post-process : Phonon-DOS,FreeEnergy,Phonon-dispersionなどの計算方法

Pre-process

前段階として，スーパーセルの変位は結晶の対称性を考慮し，単位結晶から作成され る．例として，２×２×２構造のスーパーセルを得るためにphonopyを起動する

phonopy -d --dim="2 2 2"

結果として，以下のようなファイルが生成される

% ls

disp.yaml POSCAR POSCAR-001 POSCAR-002 POSCAR-003 SPOSCAR

SPOSCARは完全なスーパーセル構造であり，disp.yamlには変位に関する情報が含まれ

ている．また，POSCAR-(001,002,003)はatom の変位のスーパーセルである.POSCAR-numberはdisp.yamlに書かれているそれぞれ異なったatomの変位と一致する．

力計算

Force constants（力定数）は構造のファイルであるPOSCAR-(001,002,003)を用いて計算 する．VASPの場合，有限変位法の計算は，単にVASP計算のPOSCARとして POSCAR-(001,002,003)を使用することができる 以下がVASPの計算で用いたINCARの例である

PREC = Accurate IBRION = -1

ENCUT = 500 EDIFF = 1.0e-08

ISMEAR = 0; SIGMA = 0.01 IALGO = 38

LREAL = .FALSE.

ADDGRID = .TRUE.

LWAVE = .FALSE.

LCHARG = .FALSE.

構造を緩和（リラックス）させないように注意すること．その後，VASPインターフェー スを使ってFORCE-SETSファイルを作成する． 微小変位したファイルそれぞれをvaspに かけること．disp-00...とはPOSCAR-numberと対応している．つまりそれぞれのvasprun ファイルを一緒に計算にかける

% phonopy -f disp-001/vasprun.xml disp-002/

vasprun.xml disp-003/vasprun.xml

Post-process mesh.conf

phonon-DOS,自由エネルギーの計算に必要なmesh.confファイルをemacsを用いて作

成する

ATOM_NAME = Al DIM = 2 2 3 MP = 8 8 8

Post-process phonon-DOS計算

POSCAR,FORCE-SET,mesh.confファイルを用意し，以下のコマンドを実行する．計 算のグラフは出力できないので値をmac上のgnuplotで表示する．

phonopy -p mesh.conf

Post-process FreeEnergy計算

POSCAR,FORCE-SET,mesh.confファイルを用意し，以下のコマンドを実行すると自 由エネルギーが算出される．F[KJ/mol]が自由エネルギー，S[J/K/mol]がエントロピー

である

vagrant:~/phonopy-1.8.4.2-rc3/Al/Al-test$ phonopy -t -p mesh.conf __ | |__ ___ _ __ ___ _ __ _ _

| ’_ \| ’_ \ / _ \| ’_ \ / _ \ | ’_ \| | | |

| |_) | | | | (_) | | | | (_) || |_) | |_| |

| .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, |

|_| |_| |___/

1.8.4.2-rc3 Mesh sampling mode

Settings:

Sampling mesh: [8 8 8]

Supercell: [2 2 2]

Spacegroup: Fm-3m (225)

Calculating force constants...

Number of irreducible q-points: 20 Calculating thermal properties...

# T [K] F [kJ/mol] S [J/K/mol] C_v [J/K/mol] E [kJ/mol]

0.000 14.0264543 0.0000000 0.0000000 14.0264543 10.000 14.0262949 0.0701190 0.2109276 14.0269961

...

Quasi-harmonic プログラム

結合エネルギー箇所

restart;with(plots):with(plottools):with(stats):with(LinearAlgebra):

with(linalg):with(ListTools):with(combinat,permute):with(Statistics):

with(StringTools):with(plots):with(LinearAlgebra):with(stats):

# E-V Curveの作成

# 結合エネルギーはユニットセルの結合エネルギーを2×2×2倍した値 p1:=[[0.95,-116.341475],[0.98,-119.249719],[0.99,-119.636394], [1.00,-119.751101],[1.01,-119.658924],[1.02,-119.361763], [1.03,-118.880567],[1.04,-118.234993],[1.05,-117.449962], [1.06,-116.543482],[1.07,-115.526365],[1.08,-114.416418]];

# eVをRyに単位変換 (8はSuperCellをUnitCellに換算) for i from 1 to nops(p1) do

p1[i][2]:=p1[i][2]*(2/27.2)/8;

end do;

# 各軸比×Al格子定数=原子が体積膨張した時の格子定数 for i from 1 to nops(p1) do

p1[i][1]:=p1[i][1]*8.0986795425999993;

end do;

# 格子定数の単位をÅからa.u.に変換 for i from 1 to nops(p1) do p1[i][1]:=p1[i][1]/0.529177;

end do;

# data11=[[各格子定数(a.u.)],[結合エネルギー(Ry)]]

# Alはfcc構造なので原子間距離は格子定数のsqrt(2)/2倍である data11:=convert(transpose(convert(p1,array)),listlist);

data1:=[(data11[1]*evalf((sqrt(2)/2))/2),data11[2]];

q1:=convert(transpose(convert(data1,array)),listlist);

# 上記の結合エネルギーの点に関数をfitting，プロットする

fit1:=fit[leastsquare[[x,y], y=a+b*x+c*x^2+d*x^3+e*x^4+f*x^5]](data1);

f2:=unapply(rhs(fit1),x);

pp1:=pointplot(q1);

pp2:=plot(f2(r),r=5..6);

display(pp2,pp1,labels=["AtomicDistance[a.u.]","BindingEnergy[Ry]"], labeldirections=[HORIZONTAL,VERTICAL]);

# Morse 型ポテンシャルの結合エネルギー

f1:=(a,b,c,d,r)->a+b*exp(-d*r)+c*exp(-2*d*r);

# 平衡原子間距離とその結合エネルギー x0:=fsolve(diff(f2(x),x),x=5..6);

y0:=f2(x0);

# fittingした関数を微分したものに平衡原子間距離を代入している.y2は２階微分

y1:=subs(x=x0,diff(f2(x),x));

y2:=subs(x=x0,diff(f2(x),x,x));

y3:=subs(x=x0,diff(f2(x),x,x,x));

# 結合エネルギーの関数f1にa,b,c,d,x0の値を代入した．まだa,b,c,dは不明 y0=f1(a,b,c,d,x0);

# 上の式をx0まわりで微分する．y2はその２階微分 subs(x=x0,diff(f1(a,b,c,d,x),x))=0;

y2=subs(x=x0,diff(f1(a,b,c,d,x),x,x));

y3=subs(x=x0,diff(f1(a,b,c,d,x),x,x,x));

# 微分から得た式を用いて方程式を解く

eqs:={y0=f1(a,b,c,d,x0),subs(x=x0,diff(f1(a,b,c,d,x),x))=0,y2=subs(x=x0, diff(f1(a,b,c,d,x),x,x)),y3=subs(x=x0,diff(f1(a,b,c,d,x),x,x,x))};

sol1:=solve(eqs,{a,b,c,d});

# 結合エネルギーの関数f1に求めたa,b,c,dを代入．新たにf3を作る f3:=unapply(subs(sol1,f1(a,b,c,d,x)),x);

pp3:=plot(f3(r),r=5..6,color=blue);

display(pp1,pp2,pp3,labels=["AtomicDistance[a.u.]","[Ry]"], labeldirections=[HORIZONTAL,VERTICAL]);

変数早見

P1: [[倍率a/a0],[基底状態エネルギー(eV)]]

data11: [[格子定数(a.u.)],[結合エネルギー(Ry)]]

data1: [[原子間距離(a.u.)],[結合エネルギー(Ry)]]

q1: data1をlist型にコンバート    fit1: fitting関数

f2: fit1の右辺を関数定義       pp1,pp2: プロット変数

f1: Morse型ポテンシャル関数    x0,y0: 平衡原子間距離,そのときのエネルギー y0,y1,y2,y3: f1の微分，0,1,2,3はその微分回数

eqs: y0,y1,y2,y3の微分を方程式としてまとめる sol1: eqs方程式を解く,a,b,c,dが算出される

f3: f1にa,b,c,dを代入して関数化   pp3: プロット関数

体積弾性率箇所

# B1×B2=体積弾性率の式

B1:=unapply((-lambda^3*exp(-lambda*r))/(12*Pi*ln(exp(-lambda*r))),r);

B2:=unapply((b+4*c*exp(-lambda*r))-(2/ln(exp(-lambda*r)))*

(b+2*c*exp(-lambda*r)),r);

# 単位変換

B3(r):=B1(r)*B2(r)*((27.2/2)/(0.529177^3))*160.218*10;

# 関数化,プロット B:=unapply(B3(r),r);

evalf(B(x0));r0:=x0;

plot(B(r),r=6..8.5,color=black,labels=["AtomicDistance[a.u.]",

"BulkModulus\UTF{008E}\UTF{0087}[Kbar]"], labeldirections=[HORIZONTAL,VERTICAL]);

変数早見

B1: 式2.18の前部 B2: 式2.18の後部

B3: B1,B2の単位変換 160.218: ev/ÅをGPa

※その他の数字に関しては赤本の冒頭にある単位に関するノートを参照 B: B1,B2の乗算

r0: 平衡原子間距離=x0

ドキュメント内 phonopy phonon-dispersion phonon-dos MedeA phonopy phonopy, phonopy MedeA phonopy MedeA Moruzzi Quasi-harmonic MedeA,Quasi-harmonic phonopy phonopy. (ページ 35-52)