熱膨張率の理論計算¶
これまでの説明では、原子間ポテンシャルの温度変化が無視されていたが、実際には温度によって変化する。
より精密に熱膨張を計算するには、この温度変化を考慮する必要がある。
この温度変化した原子間ポテンシャルの極小値から、平衡原子間距離\(r(T)\)と ヘルムホルツ自由エネルギー :math:`F(T)` が求められる。
原子間ポテンシャルの温度変化は格子振動の効果によって現れる。格子振動の効果を計算によって求めるには、
原子間に働く力
ダイナミカルマトリックスを求める。ダイナミカルマトリックスを対角化し、固有値 \(\omega(k)\) （フォノン分散）を得る。フォノン分散を波数空間で積分し、フォノンDOSを得る。フォノンDOSとボース分布関数から、有限温度における格子振動状態を見積もる。
このようにさまざまな段階を踏んで計算を行わなければならず、まともに計算をするのはかなり大変である。
ここではMoruzziらが考案したデバイ－グリュナイゼン近似により、熱膨張について考えてみる。デバイ－グリュナイゼン近似は、原子間ポテンシャルから、ヘルムホルツ自由エネルギーを計算する方法である。
課題¶
下記の例に倣い、熱膨張を計算するpythonコードを完成させて提出してください。
提出ファイルはJupyter notebookを推奨するが、それ以外でも可。ただし実行できる形で提出すること。
デバイ－グリュナイゼン近似¶
熱膨張を電子論的に計算する手法として，デバイ－グリュナイゼン近似[1],[2]による方法を説明する．
この方法としてはまず，電子状態計算により基底状態における化合物の体積と全エネルギーの関係を計算する．このときエネルギーが極小値をとる体積 \(V_0\) を含む範囲で体積を変化させる．
化合物の体積と全エネルギーのフィッティングを以下のモース型ポテンシャルにより行う．
\(r\) は原子間距離，\(r_0\) は平衡原子間距離，\(A\), \(D\), \(\lambda\)はフィッティングパラメータである．
これらのパラメータを用いると、以下の式により体積弾性率\(B(r_0)\)が求められる．
つぎにこの式を平衡体積におけるデバイ温度の式に導入する．
\(h\)はプランク定数，\(\omega_{D0}\) はデバイ振動数，\(k_B\) はボルツマン定数，\(\hslash\) はディラック定数，\(m\)は化合物A\(_x\)B\(_y\) の有効原子質量であり全質量の対数平均で定義される．
ここで\(k(v)\) はポアソン比\(\nu\)を用いて以下の式で表される．
大抵の固体はポアソン比が0.2~0.3であるため，ここでは\(\nu\) =0.2として取り扱うこととする．
デバイ温度体積の関係はグリュナイゼン定数\(\gamma\)用いて以下の式のように定義した．
また，グリュナイゼン定数は
で表される．グリュナイゼン定数の定義としてはSlater 近似[3]、Dugdale-MacDonald 近似 [4] などがあるが，ここでは広い温度域での再現性が良好な後者を使用している．
次に各体積でのヘルムホルツ自由エネルギーを算出する．
デバイ関数 \(f_D\) は以下の様に定義される。
各温度についてエネルギーの底 \(r_0\) を取り，線膨張係数\(\alpha\)を以下の式を用いて算出する．
[1] X. -G. Lu, M. Selleby and B. Sudman, Acta Materialia, 55(2007), 1215
[2] X. -G. Lu, Ma. Selleby and B. Sudman, Acta Materialia, 53(2005), 2259
[3] J.C.Slater, Introduction to chemical physics. New York (NY): McGraw-Hill; 1939
[4] J.S.Dugdale, D. K. C. MacDonald, Phys. Rev., 89(1953), 832
エネルギーの体積依存性¶
物質の全エネルギーの体積依存性について見てみよう。２つの原子を近づけていくと、ある距離までエネルギーが減少する。これは２つの原子に電子が共有されることが原因である。一方で近づきすぎると、電子同士の反発する力が強くなるためエネルギーは急上昇する。
このような全エネルギーの体積依存性は、第一原理計算と呼ばれる計算手法によって求めることが可能である。ここでは、Materials Projectに収録されている計算データを利用する。FCC構造のCuについて、全エネルギーと体積との関係を図示してみる。
[1]:
from pylab import * #this includes numpy as np!
from scipy.optimize import leastsq
plt.rcParams["figure.figsize"] = (5,4)
必要なライブラリをインポートする。
[2]:
#Cu(mp-30)
v = np.array([8.6664, 8.9586, 9.5623, 9.8741, 10.5179, 9.2571, 10.8500, 11.5351, 12.9905, 12.2484, 13.7620, 12.6158, 14.5635, 11.1890, 14.9757, 15.8231, 10.1926, 15.3955, 11.8882, 14.1590, 13.3726])
r = np.cbrt(v)
e = np.array([-3.3243, -3.4941, -3.7562, -3.8536, -3.9921, -3.6372, -4.0371, -4.0873, -4.0660, -4.0937, -4.0115, -4.0836, -3.9364, -4.0681, -3.8927, -3.7961, -3.9314, -3.8457, -4.0956, -3.9763, -4.0416])
plot(r,e,'ro')
xlabel('r ($\AA$)')
ylabel('Energy (eV)')
[2]:
Text(0, 0.5, 'Energy (eV)')
エネルギーが最安定となる\(r_0 \sim 2.3\) の周りで原子間距離を変化させたエネルギーの変化を示している
Morse ポテンシャルでフィッティング¶
この様な固体のミクロの振舞いを記述するモデルを考える。
固体の原子間に働く力を数学的に記述する2体間ポテンシャル (pair potential) にはさまざまなものがあるが，簡単なものとしてモース (Morse) 型がある。
Morseポテンシャルは，2原子分子の原子間距離\(r\)に依存するポテンシャルである.
ここで\(r\)は原子間距離。\(r_0\)は平衡原子間距離。\(A\)、\(D\)、\(\lambda\)はフィッティングパラメータである。
Cuのエネルギーと体積の関係に、モース型ポテンシャルをフィッティングしてみよう。
*****here*****
に正しい式を記入して実行してください。
[4]:
#difinition of the functions
def Morse(parameters,r):
A = parameters[0]
D = parameters[1]
r0 = parameters[2]
la = parameters[3]
#*****here*****
return E
def objective(pars,y,x):
#we will minimize this function
err = y - Morse(pars,x)
return err
[7]:
#now here are our initial guesses.
A = -4.1
D = 1
r0 = 2.3
la = 0.1
x0 = [A, D, r0, la]
rfit = np.linspace(min(r),max(r),100)
morpars, ier = leastsq(objective, x0, args=(e,r)) #this is from scipy
plot(r,e,'ro')
plot(rfit, Morse(morpars,rfit), label='Morse fit')
xlabel('r ($\AA$)')
ylabel('Energy (eV)')
legend(loc='best')
#add some text to the figure in figure coordinates
ax = gca()
text(0.4,0.5,'A = %1.2f eV' % morpars[0],transform = ax.transAxes)
text(0.4,0.4,'D = %1.2f eV' % morpars[1],transform = ax.transAxes)
text(0.4,0.3,'r0 = %1.2f $\AA$' % morpars[2],transform = ax.transAxes)
text(0.4,0.2,'lambda = %1.2f $\AA^{-1}$' % morpars[3],transform = ax.transAxes)
savefig('a-eos.png')
show()
/var/folders/qx/bp1cy1rx2g3b3fqmxf9knzph0000gn/T/ipykernel_14877/2984365570.py:10: RuntimeWarning: overflow encountered in exp
E = A + D*np.exp(-2*la*(r-r0)) -2*D*np.exp(-la*(r-r0))
/var/folders/qx/bp1cy1rx2g3b3fqmxf9knzph0000gn/T/ipykernel_14877/2984365570.py:10: RuntimeWarning: invalid value encountered in subtract
E = A + D*np.exp(-2*la*(r-r0)) -2*D*np.exp(-la*(r-r0))
二次関数の時よりもフィッティング結果が良好であることがわかる。ここで、
\(r_0\)の平衡原子間距離は2.29 Å
\(A\)=-1.31 eV、\(D\)=2.79 eV、\(\lambda\)=1.80 Å\(^{-1}\)
である。
体積弾性率B¶
これらのパラメータを使うことで、熱膨張を求める際に必要な 体積弾性率B(bulk modulus) を計算することができる。体積弾性率とは、物体に圧力を加えたときの変形のしにくさを示す指標である。体積弾性率\(B(r_0)\)は以下の式で表される。
*****here*****
に正しい式を記入して実行してください。
なお、得られた体積弾性率の単位は eV/Å\(^3\) なのでkbarへ単位変換をしてください（eVA3_to_kbar）。
[8]:
D = morpars[1]
r0 = morpars[2]
la = morpars[3]
#transformation
eVA3_to_kbar=160.218*10
x0=np.exp(-la*r0)
#*****here*****
print("B(r_0)=",B0,"kbar")
B(r_0)= 1749.9896494141262 kbar
デバイ温度¶
次にデバイ温度 \(\Theta_D\)を求める。デバイ温度とは格子振動や原子間力と深くかかわり合いを持つ温度のことを指し、物質の硬さの指標となる。一般的に物質のデバイ温度は、硬い物質ほど高く、逆に柔らかい物質ほど低い値を取る。
平衡体積におけるデバイ温度は以下の式で表される。
\(r_0\)は平衡原子間距離、Bは体積弾性率、Mは原子質量である。
デバイ温度体積の関係はグリュナイゼン定数\(\gamma\)用いて以下の式のように定義した．
またグリュナイゼン定数は
で表される。グリュナイゼン定数はSlater近似やDugdale-MacDonald近似などがあるが、ここでは広い温度域での再現性が良好なDugdale-MacDonald近似を利用している。
*****here*****
に平衡体積におけるデバイ温度の式を記入して実行してください。
[9]:
mass=63.546 #Cu
#*****here*****
Gru=la*r0/2
print("ThetaD0=",TD0,"K")
print("Gru=",Gru)
def DebyeTemp(r):
#return 67.48*sqrt(r*B/mass)
return TD0*(r0/r)**(3*Gru)
ThetaD0= 330.55437406219244 K
Gru= 2.0548457352477607
ヘルムホルツ自由エネルギーの体積依存性¶
これらの結果を利用して、次式で表されるヘルムホルツ（Helmholtz）の自由エネルギーの体積依存性を0～3000 Kまでの範囲で算出する。
デバイ関数\(f_D\)は以下の様に定義される。
[10]:
from scipy import integrate
def integrand(T):
if T==0:
return 0
return (TD/T)**3/(np.exp(TD/T)-1)
def Debye(T):
if T==0:
return 0
y, err = integrate.quad(integrand, 0, TD/T)
return (3.)*T * (T/TD)**3 * y
def FreeEnergy(r,T):
#kB=1.380662e-23 #J/K
kB=1/1.16048e4 #eV/K
F = Morse(morpars,r) - kB*T*(Debye(T)-3*np.log(1-np.exp(-TD/T))) + 9/8*kB*TD
#print(Debye(T))
return F
[12]:
#make a vector to evaluate fits on with a lot of points so it looks smooth
T_list = np.linspace(0,3000,7)
Nr=rfit.shape[0]
NT=T_list.shape[0]
def BulkModulus(r):
x=np.exp(-la*r)
return -D*r*la**3/(6*np.pi*np.log(x))*eVA3_to_kbar
with open('out.csv', 'a') as f:
for T in T_list:
F=[]
for r in rfit:
#print(T)
#F = Debye(T)
#plt.scatter(T,F)
B=BulkModulus(r)
TD=DebyeTemp(r)
#print(TD)
Tfit=np.full(Nr,T)
F.append(FreeEnergy(r,T))
#print(F)
#csvファイルとして保存
data = np.c_[Tfit,rfit,F]
np.savetxt(f,data,delimiter='\t')
#plt.figure(figsize=(16,10))
plt.scatter(rfit,F, label=T)
#plt.plot(x,F_model,label=PlotT); # 折れ線グラフのプロット
#plt.title(PlotT) # x軸ラベルの設定
plt.xlabel('r ($\AA$)') # x軸ラベルの設定
plt.ylabel('Free energy (eV/atom)'); # y軸ラベルの設定
plt.legend()
plt.show()
/var/folders/qx/bp1cy1rx2g3b3fqmxf9knzph0000gn/T/ipykernel_14877/2017970456.py:17: RuntimeWarning: divide by zero encountered in scalar divide
F = Morse(morpars,r) - kB*T*(Debye(T)-3*np.log(1-np.exp(-TD/T))) + 9/8*kB*TD
/var/folders/qx/bp1cy1rx2g3b3fqmxf9knzph0000gn/T/ipykernel_14877/2017970456.py:6: RuntimeWarning: overflow encountered in exp
return (TD/T)**3/(np.exp(TD/T)-1)
熱膨張¶
ここまでで求めたヘルムホルツ自由エネルギーから、各温度における自由エネルギーの最小値をとる原子間距離\(r_0\)を求めることができる。これをつなげていくと、\(r_0\)を\(T\)の関数で表すことができる。つまり、各温度における平衡原子間距離\(r_0\)を求めることができれば、熱膨張係数が計算できることになる。
各温度におけるヘルムホルツ自由エネルギーの最小値をとるには、pandasで最大値・最小値の行名・列名を取得するidxmax, idxminを使います。
https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.idxmin.html
[13]:
import pandas as pd
df = pd.read_csv('out.csv', delimiter='\\t', dtype='float', names=('T', 'r', 'F'))
print("Temperature list")
Tlist=df['T'].unique()
print(Tlist)
r0_T=[]
for i in Tlist:
PlotT = i
#print(\"T=\",PlotT,\"K\")
df2 = df[df['T'] == PlotT ]
x=df2['T']
r=df2['r']
F=df2['F']
#print(df2['r'][df2['F'].idxmin()])
r0_T.append(df2['r'][df2['F'].idxmin()])
plt.scatter(Tlist,r0_T)
#plt.title(PlotT) # x軸ラベルの設定
plt.xlabel('T') # x軸ラベルの設定
plt.ylabel('$r_0 (\AA)$'); # y軸ラベルの設定
plt.legend()
plt.show()
/var/folders/qx/bp1cy1rx2g3b3fqmxf9knzph0000gn/T/ipykernel_14877/2585520366.py:3: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support regex separators (separators > 1 char and different from '\s+' are interpreted as regex); you can avoid this warning by specifying engine='python'.
df = pd.read_csv('out.csv', delimiter='\\t', dtype='float', names=('T', 'r', 'F'))
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
Temperature list
[ 0. 500. 1000. 1500. 2000. 2500. 3000.]
線熱膨張係数¶
線熱膨張系数\(\alpha\)は以下の式で定義される。
*****here*****
に正しい式を記入して実行してください。
[15]:
#*****here*****
print("α=" ,alpha, "(1/K)")
α= 1.2060301913129738e-05 (1/K)
（参考）極小値の求め方¶
[16]:
import pandas as pd
from pylab import * #this includes numpy as np!
from scipy.signal import argrelmin, argrelmax
df = pd.read_csv('out.csv', delimiter='\\t', dtype='float', names=('T', 'r', 'F'))
#print("Temperature list")
Tlist=df['T'].unique()
#print(Tlist)
PlotT = 500
#print(\"T=\",PlotT,\"K\")
df2 = df[df['T'] == PlotT ]
x=df2['T']
r=df2['r']
F=df2['F']
print(df2['r'][df2['F'].idxmin()])
plt.scatter(r,F)
#plt.title(PlotT) # x軸ラベルの設定
plt.xlabel('r') # x軸ラベルの設定
plt.ylabel('Free energy (J/mol)'); # y軸ラベルの設定
plt.legend()
plt.show()
/var/folders/qx/bp1cy1rx2g3b3fqmxf9knzph0000gn/T/ipykernel_14877/1048669438.py:5: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support regex separators (separators > 1 char and different from '\s+' are interpreted as regex); you can avoid this warning by specifying engine='python'.
df = pd.read_csv('out.csv', delimiter='\\t', dtype='float', names=('T', 'r', 'F'))
No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
2.307654929948991
[ ]: