Phonon¶

1. Introduction¶

1.1 フォノン (phonon)とは¶

フォノンは、一言で表すと周期境界条件をもつ結晶における振動モードです。

詳細の説明に入る前に、Phononを目で見て感覚を掴むためこちらのサイトを見てましょう。

ecae7acec1e04a9c9a751f480373726e — フォノンのアニメーション表示。http://henriquemiranda.github.io/phononwebsite/phonon.html より。

上記Websiteへ行く -> 左上のmaterials という表から Si を選択 -> 右図の“フォノン分散図”から適当に点をクリックしてみましょう。すると、真ん中のシリコン結晶におけるシリコン原子それぞれが周期的な動きをしていることが確認できます。そして、全体として片側から反対側に波が伝わるような動きとなっていることが見て取れます。

このような結晶における振動を格子振動と呼びます。格子振動は振動数が非常に大きくなるため、通常は量子力学的な効果を考慮して取り扱われます。格子振動に量子効果を含めたものにはフォノンという名前がついています。ただし、以下の説明では複雑さを避けるため、格子振動とフォノンは特別に用語を区別せずに解説しています。より詳しくは末尾のAppendixや固体物理学の教科書を参考にしてください。

波の伝播する方向は波数ベクトル $\vec q$ で表されます。フォノンの解析を行う際に重要な問題は、ある波数ベクトル $\vec q$ に対し、どのような振動モードが許されているのか、その振動モードにおける原子の振動数がどうなるかを求めることです。

フォノンは準粒子であり、エネルギー $E$と運動量 $\mathbf{p}$をもちます。固体中のフォノンは熱物性や音波特性に深い関わりがあります。例えば、フォノンの励起は内部エネルギー $U$の増加に寄与したり、フォノンの輸送は音や熱の伝導を引き起こします。このように、フォノン解析は固体の性質を探る上で重要なシミュレーションツールとなります。

1.2 フォノン振動数の計算方法¶

前章で紹介した振動解析と同様、フォノンも調和振動子近似を用いて計算されます。ただし違いとして、各原子はUnit cellをまたいで向こうにいる原子とも相互作用をします。周期境界条件を考慮して、結晶構造から微小変位をさせた系のポテンシャルエネルギーは以下のように表されます。

\[V(\mathbf{r}) \approx V(\mathbf{r_0}) +　\frac{1}{2} \sum_{ij} \frac{\partial^2 V(\mathbf{r_0})}{\partial r_i^a \partial r_j^b} \Delta \mathbf{r}_i^a \Delta \mathbf{r}_j^b\]

ここで、$r_i^a$ $a$番目のスーパーセルに属する、$i$番目の原子の座標です。結晶における力定数マトリクス (前章で紹介) は$3NM \times 3NM$ 行列となります。 $N$ はユニットセルの原子数、 $M$ はユニットセルの数です。このため、力定数マトリクスはとても大きな行列となります。

特定の振動モードにおける、各原子の平衡点からの変位 $u$は次の平面波のような形を持つと仮定します。

\[u(\vec R, t) = A e^{ \vec{q} \vec{R} -i\omega_{m} t}\]

詳細は省きますが、この仮定のもとで実空間における離散フーリエ変換を適用することで、　大きな力定数マトリクスを$3N \times 3N$のdynamical matrixに変換することができます。 Dynamical matrixは波数ベクトル$\vec q$に依存しており、その固有値と固有ベクトルがフォノンの振動数と振動モードに相当します。詳細な導出については、末尾のappendixを参照ください。

次章では、フォノンの計算をMatlantis上で提供しているライブラリであるmatlantis-featuresを用いて行いながら、フォノンの重要なコンセプトを説明します。

2. Siでのフォノン計算¶

2.1 構造準備¶

ここでは、フォノン計算の例としてシリコン結晶を用います。まずは構造緩和を行います。フォノン分散を計算する時は、その結晶構造の基本単位胞 (primitive cell)を用いることが多いです。セル内を最小原子で構成できるためです。

[1]:

import numpy as np
np.set_printoptions(precision=3)

from ase.build import bulk
from ase.optimize import BFGS
from ase.constraints import ExpCellFilter

import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode

estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v3.0.0")
calculator = ASECalculator(estimator)
atoms = bulk("Si")
atoms.calc = calculator

exp = ExpCellFilter(atoms)
BFGS(exp).run(fmax=0.01)

      Step     Time          Energy         fmax
BFGS:    0 07:07:45       -9.108037        0.4967
BFGS:    1 07:07:45       -9.113260        0.0002

[1]:

True

2.2 力定数マトリクスの計算¶

振動解析同様、フォノンの振動数計算には力定数マトリクスが必要です。この力定数マトリクスは、ForceConstantFeatureで計算されます。

Siの基本単位胞は2原子が含まれています。この基本単位胞をたくさん繰り返したスーパーセル構造を作成することが精度の良い結果を得るために必要です。力定数マトリクスの各成分 $\frac{\partial^2 V(\mathbf{r_0})}{\partial r_i^a \partial r_j^b}$ は、$a$番目のユニットセルの$i$番目の原子と $b$番目のユニットセルの$j$番目の原子との距離で減衰していきますが、スーパーセルが十分大きい場合にのみ力定数マトリクスで無視できないような要素は含まれているとみなせます。

本チュートリアルでは、 Si結晶を 10x10x10 に拡張して計算を行います。

matlantis-featuresのForceConstantFeatureを用いて計算します。

[2]:

from matlantis_features.features.phonon import ForceConstantFeature

fc = ForceConstantFeature(
    supercell = (10,10,10),
    delta = 0.1,
)
force_constant = fc(atoms)
print(force_constant.force_constant.shape)

(1000, 6, 6)

計算結果として得られる力定数マトリクス $C$ は 1000x6x6 のshapeを持つtensorとなります。

[3]:

print(force_constant.force_constant[1])

[[-0.196 -0.195 -0.115 -0.044 -0.026  0.031]
 [-0.195 -0.196 -0.115 -0.026 -0.044  0.031]
 [ 0.115  0.115  0.43   0.031  0.031 -0.136]
 [-3.125 -2.088  2.088 -0.196 -0.195  0.115]
 [-2.088 -3.125  2.088 -0.195 -0.196  0.115]
 [ 2.088  2.088 -3.125 -0.115 -0.115  0.43 ]]

force_constant.force_constant[1] は 6x6 のtensorです。これは、1つめと2つめのユニットセルにおける力定数 $C^{1,2}$、すなわち $:nbsphinx-math:frac{partial^2 V(mathbf{r_0})}{partial r_i^1 partial r_j^2} $です。前章の振動解析同様、力定数 $C^{1,2}$は有限差分法で近似計算されています。

\[C^{1,2}[i,j]=\frac{\partial^2 V(\mathbf{r_0})}{\partial r_i^1 \partial r_j^2} \approx \frac{F(\mathbf{r_0} + \Delta r_i^1)_j^2 - F(\mathbf{r_0})_j^2}{|\Delta{r_i^1}|}\]

今回は、10x10x10 のスーパーセルを用いているため、ユニットセルの組み合わせである1000x1000 の数だけ 6x6 matrix が得られるべきです。しかし、結晶の並進対称性を考慮し、$C^{a,b} = C^{0, b-a}$ となることを利用することで、力定数マトリクスを 1000x6x6 tensorで表すことができます。

2.3 Dynamical matrixの計算¶

分子の振動数は力定数マトリクス$C$の対角化で得られましたが、結晶の振動数は実空間の力定数マトリクスを離散フーリエ変換したdynamical matrix $D_{\vec q}$ の固有値として得られます。 Dynamical matrixへの変換の際に、波数ベクトル $\vec q$ というパラメータが追加で導入されます。

ここでは、matlantis-featuresのPhononFreqencyを用いて、シリコン結晶の波数ベクトル $\vec q = [0.1, 0.0, 0.0]$ に対応するdynamical matrix を計算してみます。

[4]:

from matlantis_features.features.phonon.utils import PhononFrequency

phonon_frequency = PhononFrequency(
    force_constant.force_constant, force_constant.supercell, force_constant.unit_cell_atoms
)
dynamic_matrix = phonon_frequency._dynamic_matrix_q(np.array([0.1, 0.0, 0.0])).real
print("Dynamical matrix of q=[0.1, 0.0, 0.0]:")
print(dynamic_matrix)

Dynamical matrix of q=[0.1, 0.0, 0.0]:
[[ 0.456 -0.001 -0.001 -0.433 -0.013 -0.013]
 [-0.001  0.456  0.001 -0.013 -0.433  0.013]
 [-0.001  0.001  0.456 -0.013  0.013 -0.433]
 [-0.433 -0.013 -0.013  0.456 -0.001 -0.001]
 [-0.013 -0.433  0.013 -0.001  0.456  0.001]
 [-0.013  0.013 -0.433 -0.001  0.001  0.456]]

Dynamical matrixは 6x6 の対称行列となっています。これは、力定数マトリクス $C$ から以下の式の計算により得られています。

\[\tilde{D}_{\vec{q}}[k\alpha, k'\beta] = \frac{1}{\sqrt{M_k M_{k'}}} \tilde{C}(\vec{q})[k\alpha, k'\beta]\]

\[\tilde{C}_{\vec{q}}[k\alpha, k'\beta] = \sum_{b} {\left( \frac{\partial ^{2} E}{\partial{u_{k}^{a\alpha}} \partial{u_{k'}^{b\beta}}} \right) e^{i\vec{q}(\vec{R_b} - \vec{R_a})}} = \sum_{b} {C^{0,b}[k\alpha, k'\beta] e^{i\vec{q}(\vec{R_b} - \vec{R_a})}}\]

その後、この行列を対角化することで固有値・固有ベクトルを得ます。分子の振動解析と似たように、固有値が振動数に対応し、固有ベクトルはそのフォノンモードにおける原子の変位と関連付けられます。式の導出など、詳細はappendixを参照ください。

それでは、シリコン結晶の $\vec q = [0.1, 0.0, 0.0]$ におけるフォノンの振動数を計算してみましょう。

[5]:

frequency, eigenvector = phonon_frequency._frequency_mode_q(np.array([0.1, 0.0, 0.0]), return_mode=True)
for i in range(6):
    print(f"Band {i}")
    print(f"    Frequency: {frequency[i]:.3f} meV")
    print(f"    Eigenvector:")
    print(eigenvector[i].real)

Band 0
    Frequency: 5.910 meV
    Eigenvector:
[[-0.075  0.028 -0.103]
 [-0.074  0.029 -0.103]]
Band 1
    Frequency: 5.912 meV
    Eigenvector:
[[ 0.052  0.087 -0.035]
 [ 0.048  0.083 -0.035]]
Band 2
    Frequency: 11.188 meV
    Eigenvector:
[[ 0.073 -0.073 -0.073]
 [ 0.077 -0.077 -0.077]]
Band 3
    Frequency: 60.955 meV
    Eigenvector:
[[-0.073  0.073  0.073]
 [ 0.077 -0.077 -0.077]]
Band 4
    Frequency: 61.389 meV
    Eigenvector:
[[-0.09  -0.014 -0.077]
 [ 0.088  0.011  0.077]]
Band 5
    Frequency: 61.390 meV
    Eigenvector:
[[-0.013 -0.09   0.077]
 [ 0.01   0.088 -0.077]]

波数ベクトル :math:`vec q`

フォノンの計算では、波数ベクトル $\vec q$ が導入されました。この波数ベクトルの物理的な意味は、その振動モードの伝播する方向と関連します。また、理論的に波数ベクトル $\vec q$ は第一ブリルアンゾーン (first Brillouin zone)の中に制限できることが示せます (詳細はappendixを参照)。

ここでは、get_brillouin_zone_3d関数を使い、シリコン結晶における第一ブリルアンゾーンを可視化してみます。 Wikipedia上の以下の画像と比較してみてください

347940aa14b146b8b19b0c2b54a01bd9 — 対称性ラベルを付した面心立方格子構造(FCC構造)の第一ブリルアンゾーン。 https://en.wikipedia.org/wiki/Brillouin_zone より。

[6]:

from matlantis_features.utils.visual_utils.brillouin_zone import get_brillouin_zone_3d

fig = get_brillouin_zone_3d(atoms.cell, show_reciprocal_basis=True, show_special_k_points=True)
fig.update_layout(width=600, height=500)