エレクトロニクス、材料科学や物性物理学では電気伝導度(電気抵抗の逆数)や誘電率、帯磁率といった応答係数を計算する方法がいくつかあり、よく知られた手法の一つが線形応答理論です。これを使うと、金属における物性値をバンド図から計算することができるようになります。材料の電気特性を調べるのによく使う基本的な方法なのですが、僕は理解するまでに苦労したのでここで紹介させていただきます。
記事が長くなってしまうので、ここではバンドとは何か分かる程度の量子力学の知識を前提とさせていただきます(おそらく物理学専攻で大学3年を終わった時点くらいの知識量です)。
この記事では実際にコードを書いて簡単な例を計算してみます。
理論パートは
を見てください。文字の表記などもこの記事と同じものを使います。
ここでは一番簡単な２次元正方格子模型を考えます。
タイトバインディング模型
正方格子は基本並進ベクトルが$(a_1,a_2) = (e_x,e_y)$で、$\Delta n_i = \pm1$のときに$t$が有限です。ほかのパラメーターは電子スピン$\sigma = \pm1$を考えます。
H(q) = \begin{pmatrix}
2t (\cos{q_1}+\cos{q_2}) & 0 \\
0 & 2t (\cos{q_1}+\cos{q_2})
\end{pmatrix} \\
J_i(q) = \begin{pmatrix}
-2t \sin{q_i} & 0 \\
0 & -2t \sin{q_i}
\end{pmatrix}
(i=x,y,z)
実装
$\delta$は微小量に置き換えて、最後に実部をとることとします。
簡単なモデルの場合でも、数値誤差が出やすく、$\delta$の取り方や数値計算的なテクニックによる工夫が必要なことが多いですので、気を付けてください。例えば以下のような原因があります。
- $\Phi$の分母の$E_{j_1} - E_{j_2}$の部分が$0$に近づいて発散する。
- フェルミ/ボーズ分布がオーバーフロー/アンダーフローする。
復習すると
\Phi_{j_1j_2}^{core}(\mathbf{q},i\omega_l;\mathbf{k}) = \frac{ f[E_{j_2}(\mathbf{k+q})] - f[E_{j_1}(\mathbf{k})] }{ E_{j_1}(\mathbf{k}) - E_{j_2} (\mathbf{k+q}) + i\omega_l } \\
\Phi_{i_1 \cdots i_4}^{loc}(\mathbf{q},i\omega_l;\mathbf{k}) = \sum_{j_1,j_2} U^*_{i_1j_1}(\mathbf{k})U_{i_3j_1}(\mathbf{k})U^*_{i_4j_2}(\mathbf{k+q})U_{i_2j_2}(\mathbf{k+q})\Phi_{j_1j_2}^{core}(\mathbf{q},i\omega_l;\mathbf{k}) \\
\Phi_{i_1 \cdots i_4}(\mathbf{q},i\omega_l) = \frac{1}{L^3}\sum_{i_1\cdots i_4}\sum_\mathbf{k}\sum_{j_1,j_2} \Phi_{i_1 \cdots i_4}^{loc}(\mathbf{q},i\omega_l;\mathbf{k}) \\
\chi_{BA}(\mathbf{q},i\omega_l) = \sum_{i_1\cdots i_4} B_{i_1i_2} A_{i_3i_4} \Phi_{i_1 \cdots i_4}(\mathbf{q},i\omega_l)
を$i\omega_l \rightarrow \omega + i\delta$として計算します。
ヘッダー
import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as plt
from joblib import Parallel, delayed
パラメーター
t = 1.
system_size = 32
dimension = 2
volume = system_size**dimension
temp_limit = 4.
temp_iteration = 128
temp_widtg = temp_limit/temp_iteration
temp_offset = 1e-4
帯磁率
ここでは帯磁率$\chi_{SS}$を計算してみましょう。磁化は$S=\sigma_z$(パウリ行列)で表されます。
まずフェルミ分布の関数をつくります。オーバーフロー、アンダーフローに気を付けます。
def calc_boltzmann_factor(energy, kbt):
inv_temp = 1/(kbt)
a = -inv_temp*(energy)
if a > 709: #prevent overflow
a = 709
elif a < -745: #prevent underflow
a = -745
return np.exp( a )
def fermi_distribution(kbt,energy,chemical_potential):
return 1/(calc_boltzmann_factor(chemical_potential-energy, kbt)+1)
次に$\Phi$の各成分を計算します。幸い、ここではハミルトニアンが最初から対角なので$\Phi$も対角成分しかありません。(1)$\Phi^{core}$、(2)$\Phi^{loc}$、(3)$\Phi$の順に計算していきます。
まず、計算する$k$と温度点のリストを作成します。
def extend_k_mesh(k_mesh_list, ki_mesh_list):
return [ k + ki for k in k_mesh_list for ki in ki_mesh_list ]
def make_k_mesh_list():
global system_size, dimension
k0 = 2*np.pi/system_size
k_mesh_list = [ [ n*k0 ] for n in range(system_size) ]
ki_mesh_list = copy.deepcopy( k_mesh_list )
for i in range( dimension-1 ):
k_mesh_list = extend_k_mesh(k_mesh_list, ki_mesh_list)
return k_mesh_list
def make_temperature_list():
global temp_limit, temp_iteration, temp_width, temp_offset
return [ n*temp_width + temp_offset for n in range(temp_iteration+1) ]
反強磁性に対応する$\mathbf{q}=(\pi,\pi)$で計算すると、
$\chi_{\uparrow\uparrow} = \chi_{\downarrow\downarrow}=\chi$、$\chi_{\uparrow\downarrow} = \chi_{\uparrow\downarrow}=0$なので、$\chi_{SS}=2\chi$を計算します。
#(1)
def culc_green_function_core( E_kq, E_k, omega, kbt, chemical_potential):
global volume
delta = 0.1/volume
numerator = fermi_distribution(kbt,E_kq,chemical_potential) - fermi_distribution(kbt,E_k,chemical_potential)
denominator = ( E_kq - E_k + omega - 1j*delta )/ ( ( E_kq - E_k + omega )**2 + delta**2 )
return -numerator*denominator
# eigenenergy
def hop(k):
global t
return -2*t*( np.cos(k[0]) + np.cos(k[1]) )
#(2)
def susceptibility_loc_k(kbt,omega,q,k):
E_k= hop(k)
E_kq = hop(k+q)
return culc_green_function_core( E_kq, E_k, omega, kbt, 0)
#(3)
def susceptibility( susceptibility_index, q, omega, kbt, chemical_potential, k ):
chi0 = (1/volume)*np.sum( np.array( [ test_susceptibility_loc_k(kbt,omega,q,k) for k in make_k_mesh_list() ] ) )
return 2*chi0
def calc_kbt_dependence_of_susceptibilities():
q = np.array( [np.pi, np.pi] )
omega = 0
return [ susceptibility(kbt,omega,q) for kbt in make_temperature_list() ]
では計算してみます。
temperature = make_temperature_list()
chi1 = calc_kbt_dependence_of_susceptibilities()
plt.style.context('classic')
plt.plot(temperature, chi1)
plt.xlim([0,temp_limit])
plt.xlabel("temperature")
plt.ylabel("susceptibility")
plt.ylim([-0.001,2.5])
plt.show()
結果は以下のようになります。
参考文献
関連