ForceConstantFeature

はじめに#

ForceConstantFeature
は結晶の力定数を求めるのに用います。
力定数の値から、フォノンの分散や状態密度、いくつかの熱的な特性などを求めることができます。

計算手法#

ForceConstantFeature では VibrationFeature と同じく、結晶の力定数を調和近似により計算します。
結晶の周期性を取り入れるため、ユニットセルを3方向に繰り返したしたスーパーセルを作ります。(たとえば (4 times 4 times 4) 倍など)
このスーパーセルのポテンシャルエネルギーは次のように近似されます。

[V(mathbf{r}) approx V(mathbf{r_0}) + frac{1}{2} sum_{ijabalphabeta} frac{partial^2 V(mathbf{r_0})}{partial r_i^{a alpha} partial r_j^{b beta}} Delta mathbf{r}_i^{aalpha} Delta mathbf{r}_j^{bbeta}]

ここで、 (alpha) , (beta) はx, y, zを、 (i) , (j) は原子のインデックスを、
そして (a) , (b) はユニットセルのインデックスを表します。
(r_0) は原子の平衡位置です。
(r_i^{aalpha}) は (i) 番目の原子の (a) 番目のユニットセルでの (alpha) 方向の座標を表します。

上式のポテンシャルエネルギーの2階微分 (frac{partial^2 V(mathbf{r_0})}{partial r_i^{a alpha} partial r_j^{b beta}}) は力定数と呼ばれ、これを (C) で表わします。
matlantis-features では、力の微分を有限差分法から計算します。

[C=frac{partial^2 V(mathbf{r_0})}{partial r_i^{a alpha} partial r_j^{b beta}} = -frac {F_i^{aalpha}(Delta r_j^{bbeta}) – F_i^{aalpha}(-Delta r_j^{bbeta})} {2Delta r_j^{bbeta}}]

ここで、 (F_i^{aalpha}(Delta r_j^{bbeta})) は、
セル (b) 内の原子 (j) がその平衡位置から (beta) 方向に距離 (Delta r_j^{bbeta}) だけ動いたときの、
セル (a) 内の原子 (i) に作用する力の (alpha) 方向成分を示しています。

結晶の力定数は (N) をユニットセルあたりの原子数、
(M) をユニットセルの総数とすると、
(M times M times 3N times 3N) のテンソルになりますが、
結晶の並進対称性から (M times 3N times 3N) だけ考慮すればよいことがわかります。
たとえば、ユニットセルを (4 times 5 times 7) 倍した場合、
得られる力定数のサイズは (140 times 3N times 3N) となります。
最初のスライス、 (C[0]) はセル (0,0,0) の自分自身との相互作用を表します。
第2スライス、 (C[1]) はセル (0,0,0) とセル (0,0,1) との相互作用を表します。

力定数と実験から得られる物性値を直接比較することはできませんが、Post-calculation Featuresを用いて力定数をさらに解析することで、フォノンのバンドや状態密度および一部の熱力学特性を得ることができます。

Tips#

ForceConstantFeature の入力として使用できるのは安定構造のみです。
原子位置をあらかじめ最適化してからご使用ください。
入力構造にはprimitive cellの使用が推奨されます。
conventional cellでも動作しますが、結果から得られるフォノンのバンド構造が煩雑になってしまいます。

primitive_matrix 引数を指定することで手動または自動的に入力構造をprimitive cellに変換することができます。
```
# Manually set the basis vector of the primitive cell
force_constant = fc(
  atoms_Si,
  primitive_matrix=np.array(
      [
          [0.00, 0.25, 0.25],
          [0.25, 0.00, 0.25],
          [0.50, 0.50, 0.00],
      ]
  ),
  prec=1e-3,
)

# automatically detect the primitive cell
force_constant = fc(atoms_Si, primitive_matrix="auto", prec=1e-3)
```
この例の前半では、primitive cellの基本格子ベクトルを

[begin{split}a^{prime}=0 times a + 0.25 times b + 0.5 times c \
b^{prime}=0.25 times a + 0 times b + 0.5 times c \
c^{prime}=0.25 times a + 0.25 times b + 0 times cend{split}]

として指定しています。
ここで、 (a), (b), (c) は入力構造の基本格子ベクトルです。注意点として、手動で設定する場合には primitive_matrix を正確に入力するようにしてください。
もともとのcellにおいて折りたたまれる対象となる原子は、折りたたみ後の座標がprimitive cell内の原子のどれかと重なっている(距離が prec 以下である)必要があります。そうでない場合、例外 ValueError が生じます。

一方、 primitive_matrix="auto" を指定した場合、 spglib により対称性から自動的にprimitive cellを判定します。
パラメータ supercell を指定することで、ユニットセルを繰り返してスーパーセルを作成することができます。
スーパーセルは、力定数における無視できないセル間相互作用をすべて考慮できるように、十分に大きくする必要があります。
とるべきスーパーセルの大きさは材料に強く依存するため決定するのは困難ですが、
大きさを変化させながら収束性の確認を行なうことが望ましいです。
パラメータ delta により有限差分における変位の大きさを指定することができます。
PFPの出力に含まれる微細なノイズの影響を抑えるため、比較的大きな値（0.1程度）に設定することを推奨します。