銅の EXAFS デバイワラー因子と熱振動効果の第一原理計算

EXAFS（広域X線吸収微細構造）は、吸収原子のまわりの局所構造（配位数・結合距離・乱れ）を元素選択的に測定できる強力な手法で、触媒・電池材料・ナノ粒子・合金など結晶性の低い系の構造解析に広く使われています。EXAFS 振動は原子の熱振動によって $e^{-2\sigma^2 k^2}$ の形で減衰し、この減衰因子 $\sigma^2$（EXAFS デバイワラー因子）は同一ボンド両端の相対変位のボンド方向成分（MSRD: mean-square relative displacement）で決まります。隣接原子は音響フォノンでほぼ同位相に動くため、$\sigma^2$ は X線回折の B 因子（独立変位）から単純には求められず、格子力学の情報そのものを含んでいます。実測スペクトルの解析では $\sigma^2$ は配位数と強く相関するフィットパラメータであり、これを第一原理から与えられれば解析の定量性が大きく向上します。本事例では、第一原理計算ソフトウェア Advance/PHASE を力・エネルギーの評価エンジンとし、力定数フィッティングに ALAMODE、X線吸収スペクトル計算に FDMNES を連携させて、fcc 銅の Cu K 吸収端 EXAFS を「熱振動スナップショットの配置平均」として算出しました。さらに、$\sigma^2$ を、フォノンモード和・配置統計・スペクトル解析という 3 つの独立な方法で求めて相互検証を行い、Cu 箔の実験値と比較します。

Keywords: 第一原理計算 (DFT), EXAFS, XAFS, デバイワラー因子, MSRD, フォノン, 準調和近似 (QHA), 零点振動, FDMNES, ALAMODE, 銅

1. 計算方法：フォノンから EXAFS スペクトルまで

1.1 EXAFS デバイワラー因子 = ボンド射影の相対変位

第 1 配位圏の EXAFS 減衰を支配する $\sigma^2$ は、吸収原子 0 と散乱原子 $j$ の変位 $\mathbf{u}$ を用いて

$$\sigma^2 = \big\langle [(\mathbf{u}_j-\mathbf{u}_0)\cdot\hat{R}]^2 \big\rangle = 2\langle u_\parallel^2\rangle - 2\langle u_{0\parallel} u_{j\parallel}\rangle$$

と定義されます [1, 2]。調和力定数からスーパーセルの固有モード $(\omega_s, \mathbf{e}_s)$ を求めれば、$\sigma^2$ は量子統計（零点振動込み）のモード和

$$\sigma^2(T) = \sum_s \langle A_s^2\rangle\, \big|(\mathbf{e}_{s,j}-\mathbf{e}_{s,0})\cdot\hat{R}\big|^2 /m, \qquad \langle A_s^2\rangle = \frac{\hbar}{2\omega_s}\coth\!\frac{\hbar\omega_s}{2k_BT}$$

で厳密に評価できます（古典極限では $\langle A_s^2\rangle = k_BT/\omega_s^2$）。

1.2 準調和近似と熱振動スナップショット

格子定数の温度依存 $a(T)$ は、ALAMODE [3] を用いて3つの体積における自由エネルギー $F(V,T)=E_{\rm BM}(V)+F_{\rm ph}(V,T)$ の最小化（準調和近似, QHA）で決めます。その $a(300\,\mathrm{K})$ の格子上で、固有モードの正規座標に量子統計の振幅を与えた乱数変位

$$\mathbf{u} = \sum_s \mathbf{e}_s\, \xi_s \sqrt{\langle A_s^2\rangle}/\sqrt{m}, \qquad \xi_s \sim N(0,1)$$

で熱振動スナップショットを生成します。零点振動を含む量子分布からのサンプリングであり、古典 MD では失われる低温の量子効果を正しく取り込めます。

1.3 FDMNES による X 線吸収スペクトルと χ(k) 解析

各スナップショットの原子構造に対し、FDMNES [4, 5] の Green 関数（多重散乱）モードで Cu K 端の吸収断面積 $\mu(E)$ を光電子エネルギー約 800 eV（$k \lesssim 14.5\ \mathrm{Å}^{-1}$）まで計算します。FDMNES の Xan_atom 機能で同時出力される原子吸収 $\mu_{\rm at}(E)$ を用いれば、バックグラウンド関数のフィットなしに $\chi = (\mu-\mu_{\rm at})/\mu_{\rm at}$ が直接得られます。配置平均したスペクトルにコアホール寿命と光電子減衰のローレンツ幅 $\Gamma(E)$（Cu K殻 $\Gamma_{\rm hole}=1.55$ eV [6]）を畳み込み、$k^2\chi(k)$ のフーリエ変換・第 1 殻フィルタを経て $\sigma^2$ を抽出します。抽出には、無振動の静的参照スペクトルに $e^{-2\sigma^2k^2}$ を掛けたモデルへ実測と同一の窓関数処理を適用して合わせ込む「窓一貫フィット」を用います（従来の ln 比法にある R 窓切り出しの系統誤差が相殺されます）。全体のワークフローを図1に示します。

図1.　Cu K端 EXAFS デバイワラー因子の算出ワークフロー: Advance/PHASE × ALAMODE × FDMNES。①状態方程式: 4 原子 fcc セルの E–V 曲線 9 点から $a_0, B_0$ を決定。②フォノン: 3 体積の 32 原子スーパーセルで alm suggest の変位構造に Advance/PHASE の力を与え、alm optimize で調和力定数を取得。③QHA・サンプリング: $a(T)$・モード和 $\sigma^2(T)$ を計算し、量子統計の熱振動スナップショットを生成。④EXAFS: FDMNES（Green モード）で配置ごとの $\mu(E)$ を計算・平均し、$\chi(k)$ から $\sigma^2$ を抽出。

2. 計算モデルと計算条件

状態方程式には fcc 慣用セル（4 原子）、フォノンとスナップショットには 2×2×2 スーパーセル（32 原子）を用いました。交換相関汎関数は GGA-PBE、擬ポテンシャルは PAW です。主な計算条件を表1に示します。FDMNES 計算は静的参照 2 件 + スナップショット 32 件の計 34 件です。

表1.　計算条件の概要

項目	設定
対象	Cu（fcc, 空間群 \(Fm\bar{3}m\)）, K 吸収端 (8979 eV)
状態方程式用セル	fcc 慣用セル 4 原子
フォノン用スーパーセル	2×2×2 慣用セル（32 原子）× 3 体積（0.995 / 1.000 / 1.005）
交換相関汎関数 / 擬ポテンシャル	GGA-PBE / PAW
カットオフ（波動関数 / 電荷密度）	25 / 230 Rydberg
k 点メッシュ（EOS / スーパーセル）	8×8×8 / 2×2×2
変位量（有限変位法・調和）	0.02 Å（3 体積で共通）
熱振動スナップショット	16 配置 × 吸収サイト 2（量子統計, 300 K, \(a(T)\) 格子）+ 静的参照
FDMNES	Green（多重散乱）モード, クラスタ半径 6.0→5.0 Å（エネルギー依存）, −20〜800 eV, Xan_atom
スペクトル解析	\(\Gamma_{\rm hole}=1.55\) eV + arctan \(\Gamma(E)\), k 窓 3–12 Å⁻¹（k² 重み）, R 窓 1.6–3.0 Å, フィット 4–9 Å⁻¹

3. 計算結果と考察

3.1 状態方程式と準調和近似による熱膨張

E–V 曲線 9 点の Birch–Murnaghan フィットから $a_0=3.660$ Å、$B_0=136.2$ GPa、$B_0'=4.86$ を得ました（図2）。体積弾性率は実験値 137 GPa [7] と 1% 未満の差で一致します。格子定数は 295 K の実験値 3.6147 Å [7] より +1.25% 大きく、これは GGA-PBE の既知の傾向です。この差は後述のとおりフォノンをやや軟化させ、$\sigma^2$ を 1 割程度押し上げる方向に働きます。

図2.　Cu のエネルギー–体積曲線と Birch–Murnaghan（3 次）フィット。赤点が Advance/PHASE（GGA-PBE）による 9 体積の全エネルギー。平衡点で $a_0=3.6599$ Å、$B_0=136.2$ GPa、$B_0'=4.86$。緑破線は実験格子定数（295 K）。

3 体積の振動自由エネルギーから QHA で求めた $a(T)$ と線膨張係数 $\alpha(T)$ を図3に示します。零点振動により $a(0\,\mathrm{K})$ は静的平衡値より +0.19% 膨らみ（3.667 Å）、300 K で 3.680 Å に達します。$\alpha(300\,\mathrm{K}) = 16.1\times10^{-6}\ \mathrm{K}^{-1}$ は実験値 $16.6\times10^{-6}\ \mathrm{K}^{-1}$ [8] とよく一致し、フォノンの体積依存（グリューナイゼン物理）が正しく捉えられていることを示します。以降のスナップショット計算はこの $a(300\,\mathrm{K})$ の格子上で行います。

図3.　準調和近似による Cu の (a) 格子定数 $a(T)$ と (b) 線膨張係数 $\alpha(T)$。水平破線は零点振動を含まない静的 $a_0$。星印は300 K の実験値。

3.2 デバイワラー因子 σ²(T)：零点振動と変位相関

スーパーセル固有モード和による第 1 殻 MSRD $\sigma^2(T)$ を図4に示します。この図には、本事例における物理的な知見が凝縮されています。第一に、量子統計のモード和（青実線）は $T\to0$ で零点振動の極限値 $\sigma^2 = 3.4\times10^{-3}\ \mathrm{Å}^2$ に収束し、実験の外挿値 $\approx3.0\times10^{-3}\ \mathrm{Å}^2$ [9] を再現します。古典統計（青破線）は $T\to0$ で 0 になり、デバイ温度（$\theta_D\approx343$ K）以下では量子効果が無視できないことがわかります。第二に、非相関極限 $2\langle u_x^2\rangle$（紫点線）は $\sigma^2$ を約 4 割過大評価し、隣接原子の変位相関を正しく扱うことが EXAFS デバイワラー因子の定量に不可欠であることを示します。第三に、EXAFS 解析で慣用される相関デバイモデル（$\theta_{cD}=315$ K [1, 9], 緑一点鎖線）は実測値とともに本計算のモード和ともよく重なります。そして紫のひし形 ―― FDMNES 配置平均スペクトルの解析から抽出した $\sigma^2$ ―― がモード和の曲線上（配置統計と 1% 差）に載ることが、フォノン計算からスペクトル解析までのワークフロー全体の内部整合性を実証しています。

図4.　Cu 第 1 配位圏の EXAFS デバイワラー因子 $\sigma^2(T)$。フォノンモード和（量子/古典）、相関デバイモデル（$\theta_{cD}=315$ K）、非相関極限 $2\langle u_x^2\rangle$、FDMNES 配置平均スペクトルからの抽出値（ひし形）、実験値（星: Cu 箔 295 K, 四角: 零点極限 [9]）。低温の平坦部が零点振動。

3.3 配置平均 EXAFS スペクトル

静的格子と熱振動配置平均の $k^2\chi(k)$ を図5に、そのフーリエ変換を図6に示します。静的スペクトル（破線）は減衰せず高 k 領域まで振動が成長するのに対し、わずか 32 配置の平均（実線）で EXAFS 特有の $e^{-2\sigma^2k^2}$ 減衰が現れます。R 空間では第 1 殻ピークの振幅が約 1/3.4 に低下します。第 2 殻以遠のピークも、変位相関の弱さ（＝大きな $\sigma^2$）のため大きく減衰します。実測の室温 Cu 箔スペクトルでは、第 2〜第 4 殻（とくに直線配置の多重散乱によるフォーカシング効果で増強される第 4 殻）のピークが分離して観測されます。本計算でもこれらのピークは残存しており、静的格子に比べて高次殻の寄与が熱振動により強く抑圧される様子が再現されています。

図5.　FDMNES による Cu K端 EXAFS $k^2\chi(k)$。破線: 無振動の静的格子（$a(300\,\mathrm{K})$）、実線: 熱振動スナップショット 32 スペクトルの配置平均。配置平均だけで実測同様の熱減衰が現れる。

図6.　$k^2\chi(k)$ のフーリエ変換振幅（位相補正なし）。点線は fcc の配位殻の結晶学的距離 $R_1$–$R_4$。熱振動により第 1 殻振幅は約 1/3.4 に低下し、遠距離殻も大きく減衰する。黄色帯は $\sigma^2$ 抽出に用いた第 1 殻フィルタ窓。

スペクトルから窓一貫フィットで抽出した $\sigma^2$ を、他の 3 手法および実験と表2で比較します。モード和・配置統計・スペクトル解析の 3 者は 1–5% の差で内部整合し、実験値との差 +7〜12% は GGA の格子軟化（$a$ +1.3%）と整合する系統差です。なお、非調和性の効果を取り込む場合には、本ワークフローを第一原理 MD スナップショット（同一の枠組みで実行可能）に置き換えることで評価が可能です。

表2.　Cu 第 1 殻デバイワラー因子 σ²（300 K）の比較

方法	\(\sigma^2\) [10⁻³ Ų]	備考
フォノンモード和（量子, 零点込み）	9.95	32 原子スーパーセル, 300 K
フォノンモード和（古典）	9.53	300 K での量子補正は +4%
スナップショット配置統計	9.58	16 配置 × 12 ボンド × 2 サイトの実空間分散
FDMNES スペクトル解析（窓一貫フィット）	9.47	配置平均 \(\chi(k)\), フィット範囲 k = 4–9 Å⁻¹
相関デバイモデル（\(\theta_{cD}=315\) K, 295 K）	8.80	EXAFS 解析の慣用モデル [1]
実験（Cu 箔, 295 K）	8.85 ± 0.5	標準試料の代表値 [9]
実験（零点振動極限, T→0）	2.98 ± 0.4	本計算（量子モード和）: 3.37

4. まとめ

Advance/PHASE（力・エネルギー）、ALAMODE（力定数）、FDMNES（X線吸収スペクトル）を連携させ、fcc 銅の Cu K端 EXAFS を熱振動スナップショットの配置平均として第一原理計算しました。準調和近似による熱膨張 $\alpha(300\,\mathrm{K})$ は実験と 3% の差で一致し、EXAFS デバイワラー因子 $\sigma^2$ はフォノンモード和・配置統計・スペクトル解析の 3 つの独立な経路が 1–5% の差で内部整合したうえで、実験値を 1 割弱の系統差（交換相関汎関数の格子定数精度に対応）で再現しました。零点振動・隣接原子の変位相関・遠距離殻の選択的な減衰といった EXAFS 熱効果の物理を、フォノン計算とスペクトル計算の橋渡しによって定量的に扱えることを示しています。実測 EXAFS 解析への $\sigma^2$ の第一原理拘束（配位数決定の精度向上）、ナノ粒子・合金・表面吸着種への展開、また同一の枠組みで第一原理 MD スナップショットを用いる非調和効果（第 3 キュムラント・高温融解前駆現象）の評価など、局所構造解析と熱物性計算をつなぐ応用が期待できます。

本解析の詳細や、研究への適用可能性に関するご相談はこちら

お問い合わせ

参考文献

1. E. Sevillano, H. Meuth, and J. J. Rehr, "Extended x-ray absorption fine structure Debye-Waller factors. I. Monatomic crystals", Phys. Rev. B 20, 4908 (1979).
2. G. Beni and P. M. Platzman, "Temperature and polarization dependence of extended x-ray absorption fine-structure spectra", Phys. Rev. B 14, 1514 (1976).
3. T. Tadano, Y. Gohda, and S. Tsuneyuki, "Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations", J. Phys.: Condens. Matter 26, 225402 (2014).
4. O. Bunău and Y. Joly, "Self-consistent aspects of x-ray absorption calculations", J. Phys.: Condens. Matter 21, 345501 (2009).
5. J. D. Bourke, C. T. Chantler, and Y. Joly, "FDMX: extended X-ray absorption fine structure calculations using the finite difference method", J. Synchrotron Rad. 23, 551 (2016).
6. M. O. Krause and J. H. Oliver, "Natural widths of atomic K and L levels, Kα X-ray lines and several KLL Auger lines", J. Phys. Chem. Ref. Data 8, 329 (1979).
7. H. M. Ledbetter and E. R. Naimon, "Elastic Properties of Metals and Alloys. II. Copper", J. Phys. Chem. Ref. Data 3, 897 (1974).
8. G. K. White, "Thermal expansion of reference materials: copper, silica and silicon", J. Phys. D: Appl. Phys. 6, 2070 (1973).
9. S. a Beccara, G. Dalba, P. Fornasini, R. Grisenti, F. Pederiva, A. Sanson, D. Diop, and F. Rocca, "Local thermal expansion in copper: Extended x-ray-absorption fine-structure measurements and path-integral Monte Carlo calculations", Phys. Rev. B 68, 140301(R) (2003).

関連ページ

• 第一原理計算ソフトウェア Advance/PHASE
• 解析分野：ナノ・バイオ
• 産業分野：材料・化学