マテリアル計算科学入門 — 目次 第V部 分子から固体へ / 第13章

第13章van der Waals力とMadelungエネルギー

第11章では永年方程式から共有結合とイオン結合を導き，第12章では対称性（点群と指標表）を使って，どの原子軌道どうしが結合をつくれるかを定性的に判定できるようになった．ここまでは主役がつねに「電子の波動関数の重なり」であった．しかし現実の物質を凝集させている力は，それだけではない．無極性の希ガス原子ですら，十分に冷やせば液体になり固体になる．電子雲がほとんど重ならない距離でも働く弱い引力——van der Waals力（van der Waals force）——が存在するからである．一方，NaCl のようなイオン結晶では，逆に電子の授受が完全に起きてしまい，あとに残るのは点電荷の集合体である．この場合，結晶を安定化させているのは無限に続くクーロン和，すなわちMadelungエネルギー（Madelung energy）である．本章では，この「もっとも弱い結合」と「もっとも古典的で強い結合」の両極端を，いずれも数式を一行も飛ばさずに導く．前者では量子力学の零点振動が主役になり，後者では純粋に幾何学的な定数（Madelung定数）が主役になる．最後に，実際の研究で Madelung エネルギーがどう使われているかを，担当教員自身の研究例で紹介する．

この章で学ぶこと
  • 中性・無極性の原子どうしが引き合う理由——「電気的つり合いが瞬間的に破れる」とはどういうことか
  • 2個の連成振動子モデルからの van der Waals 力の完全な導出：ハミルトニアン $\Ham=\Ham_0+\Ham_1$ の設定，$R\gg\abs{x_1},\abs{x_2}$ での展開，双極子–双極子相互作用 $-2\ee^2x_1x_2/(4\pi\eps R^3)$ の出現
  • 基準座標（45度回転）$x_s=(x_1+x_2)/\sqrt{2}$, $x_a=(x_1-x_2)/\sqrt{2}$ による対角化と，2つの独立な調和振動子への分離
  • 零点エネルギーの差が $\Delta U=-A/R^6$ を生むこと．量子力学がなければ van der Waals 力は存在しないという帰結，および Lennard-Jones ポテンシャル
  • イオン結晶の静電エネルギーと Madelung 定数 $A_\mathrm{M}$ の定義，1次元イオン結晶での $A_\mathrm{M}=2\ln 2$ の手計算
  • 3次元の Madelung 定数（NaCl型 1.7476 など）と，Born–Mayer 型斥力を加えた凝集エネルギー・体積弾性率の見積もり
  • アルカリハライドの化学的傾向：$R_0$ が小さいほど硬く，バラバラにしにくい
  • pymatgen（Ewald法）による Madelung エネルギーの実計算と，研究における応用例（構造相転移・表面再構成）
前提：調和振動子の量子力学（$E_n=(n+\tfrac12)\hbar\omega$，学部2年「量子力学」），Born–Oppenheimer近似（第8章），共有結合とイオン結合の永年方程式（第11章），Taylor展開・Maclaurin展開と等比級数（付録A「数学の道具箱」）．線形代数の直交変換（回転行列）も使うが，必要な事項はその都度示す．

13.1 分子から固体へ — 無極性分子さえ凝集させる力

13.1.1 4種類の化学結合

物質を凝集させる相互作用は，伝統的に次の4つに分類される．第11章・第12章で扱ったのは，このうち上の2つであった．

表13.1 化学結合の4分類と結合エネルギーの目安（1原子対あたり）．数値は代表的な値であり，物質によって大きく変わる．
種類本質結合エネルギーの目安代表例
共有結合原子軌道の重なりによる電子の共有（第11・12章）200〜1000 kJ/molダイヤモンド，Si，H$_2$
イオン結合電子の完全な授受と点電荷間のクーロン力600〜1000 kJ/molNaCl，MgO
金属結合非局在化した電子の海（第14章で扱う）100〜800 kJ/molNa，Fe，Cu
van der Waals結合電気的つり合いの瞬間的な破れ0.3〜10 kJ/molAr結晶，分子性結晶，黒鉛の層間

van der Waals結合のエネルギーは，他の3つに比べて2〜3桁も小さい．それにもかかわらず本章でわざわざ丁寧に扱うのは，この力が例外なくすべての原子・分子の間に働くからである．共有結合もイオン結合もつくれない希ガス原子ですら，Ar は 84 K で融け，そのわずかな引力で結晶をつくる．生体高分子の折りたたみ，黒鉛の層間のはがれやすさ，ヤモリが壁に張りつく仕組みまで，van der Waals力が支配している現象は多い．

なぜ？：中性で無極性なのに，なぜ引き合うのか

高校の化学では「無極性分子さえ凝集させるのが van der Waals 力である」と習う．しかし，共有結合とイオン結合を数式で理解したあとにこの説明を読み返すと，かえって分からなくなる．Ar 原子は電気的に中性で，電子雲は球対称だから永久双極子モーメントも持たない．静電気学の常識では，中性で無極性の2物体の間にクーロン相互作用は働かないはずである．それなのに，なぜ引力が生じるのか．

結論を先に述べておく．中性原子の電気的つり合いは，「時間平均で」成り立っているにすぎない．電子は原子核のまわりを量子力学的にゆらいでおり，ある瞬間を切り取れば電子雲の重心は核の位置からずれている．つまり瞬間双極子（instantaneous dipole）が絶えず生まれては消えている．そしてこの瞬間双極子は，隣の原子に双極子を誘起する．誘起された双極子は，もとの双極子と必ず引き合う向きになる．したがって，時間平均をとっても引力だけが残る．本章の13.2〜13.5節は，この定性的な話を，量子力学の厳密解が使えるモデルで最後まで計算しきる作業である．

13.1.2 van der Waals力の2つの成分

「電気的つり合いが破れる」原因は，大きく2つに分けられる．

  1. 電子の量子論的なふるまいによる自発的な分極．外から何もしなくても，電子の位置は不確定性原理のためにゆらぐ．このゆらぎが隣の原子と相関することで生じる引力をLondon分散力（London dispersion force）と呼ぶ．無極性分子どうしに働くのはこの成分だけである．
  2. 外部電荷や永久双極子による分極．近くにイオンや極性分子があると，その電場によって双極子が誘起される．これを誘起双極子（induced dipole）による相互作用と呼ぶ．

本章で導出するのは 1 の London 分散力である．これがもっとも普遍的で，かつもっとも「量子力学らしい」からである．

(a) 時間平均：無極性 電子雲の重心＝核の位置 双極子モーメント = 0 (b) ある瞬間：ゆらぐ δ− δ+ 瞬間双極子 電子雲の重心が核からずれる (c) 隣の原子が誘起分極 δ+ δ− 引力（頭と尾が向き合う）
図13.1 London分散力の描像．(a) 時間平均すれば中性原子は無極性である．(b) しかし電子の量子的なゆらぎのため，ある瞬間には電子雲の重心が核からずれ，瞬間双極子が生じる．(c) その電場が隣の原子を分極させる．誘起された双極子は必ず「頭と尾が向き合う」配置になるので，平均をとっても引力だけが残る．

記法：本書で使う静電気の定数

本書はSI単位系で通す（仕様どおり原子単位系は使わない）．電気素量を $\ee$，真空の誘電率を $\eps$ と書き，2つの点電荷 $q_1,q_2$ が距離 $r$ にあるときの静電ポテンシャルエネルギーを

$$ U=\frac{1}{4\pi\eps}\frac{q_1q_2}{r} $$

と書く．同符号なら $U\gt 0$（斥力），異符号なら $U\lt 0$（引力）である．以後，$\dfrac{\ee^2}{4\pi\eps}$ という組み合わせが繰り返し現れる．13.8節でこの値を実用的な単位で与える（$14.400\ \mathrm{eV}\,\text{Å}$）ので，覚えておくと計算が速い．

13.2 連成振動子モデルを立てる

13.2.1 モデルの設定

「電子雲がゆらぐ」を数式にするために，原子をできるだけ単純化する．もっとも簡単な模型は，核にバネでつながれた1個の電子である．これを2組用意し，距離 $R$ だけ離して直線上に並べる．

定義：van der Waals力の連成振動子モデル

直線（$x$ 軸）上に4個の粒子を置く．粒子1と粒子3は陽子（電荷 $+\ee$），粒子2と粒子4は電子（電荷 $-\ee$，質量 $m$）である．粒子1と2，粒子3と4がそれぞれバネ定数 $C$，自然長 $L$ のバネで結ばれている．Born–Oppenheimer近似（第8章）により，陽子は重いので静止していると考え，その間隔を $R$ に固定する．

つり合いの位置からの電子の変位を，粒子2について $x_2$，粒子4について $x_1$ とする．粒子の座標は

$$ \text{陽子1}:0,\qquad \text{電子2}:L+x_2,\qquad \text{陽子3}:R,\qquad \text{電子4}:R+L+x_1 $$

である．すなわち2つの「原子」は同じ向きを向いており，双極子モーメントはともに $-x$ 方向を向く．

x バネ定数 C バネ定数 C + − + − R（陽子間距離：固定） L + x₂ L + x₁ 陽子1 電子2 陽子3 電子4 2–3間：R − (L + x₂) 2–4間：R + L + x₁ − (L + x₂) 1–4間：R + L + x₁
図13.2 van der Waals力を導く連成振動子モデル．陽子（赤）はBorn–Oppenheimer近似により固定され，電子（青）だけがバネにつながれて変位 $x_1$, $x_2$ で振動する．4通りの粒子対（1–3, 2–4, 1–4, 2–3）にクーロン相互作用が働く．

なぜ？：どうして原子を「バネ」で表してよいのか

電子は本当はバネにつながれてなどいない．しかし，電子が核のまわりのつり合い位置から少しだけずれたときの復元力は，どんな束縛ポテンシャル $V(x)$ でも，つり合い点 $x=0$ のまわりで Taylor 展開すれば

$$ V(x)=V(0)+\underbrace{V'(0)}_{=0}x+\frac{1}{2}V''(0)x^2+\cdots\simeq \text{const.}+\frac{1}{2}Cx^2 $$

となり，必ず調和振動子になる（$C\equiv V''(0)$）．つり合い点では $V'(0)=0$ だから1次の項が消えるのがポイントである．すなわち「小さなゆらぎは必ず調和振動子」であり，これは物理学のあらゆる場面で使われる最強の近似である．しかも調和振動子は量子力学で厳密に解ける数少ない系のひとつだから，あとの計算がすべて手で追える．バネを持ち出すのは手抜きではなく，正しい第一近似なのである．

13.2.2 ハミルトニアンを書き下す

系の全ハミルトニアンを，バネの部分 $\Ham_0$ と，2つの「原子」の間のクーロン相互作用 $\Ham_1$ に分けて書く．陽子は固定されているので運動エネルギーを持たない．電子2の運動量を $p_2$，電子4の運動量を $p_1$ と書くと，

$$ \begin{equation} \Ham_0=\frac{p_1^2}{2m}+\frac{1}{2}Cx_1^2+\frac{p_2^2}{2m}+\frac{1}{2}Cx_2^2 . \label{eq:13-h0} \end{equation} $$

これは独立な2つの調和振動子である．角振動数はどちらも

$$ \begin{equation} \omega_0=\sqrt{\frac{C}{m}} \label{eq:13-omega0} \end{equation} $$

であり，基底状態のエネルギーは $\tfrac12\hbar\omega_0$ ずつ，合わせて $\hbar\omega_0$ である．

次に，2つの原子の間に働くクーロン相互作用を書く．図13.2から，粒子の組み合わせは 1–3（陽子どうし），2–4（電子どうし），1–4，2–3 の4通りである（1–2 と 3–4 はバネの中に繰り込まれているので数えない）．それぞれの距離は図13.2に示したとおりだから，

$$ \begin{equation} \Ham_1=\frac{1}{4\pi\eps}\left[ \underbrace{\frac{\ee^2}{R}}_{1\text{–}3} +\underbrace{\frac{\ee^2}{R+L+x_1-(L+x_2)}}_{2\text{–}4} -\underbrace{\frac{\ee^2}{R+L+x_1}}_{1\text{–}4} -\underbrace{\frac{\ee^2}{R-(L+x_2)}}_{2\text{–}3} \right]. \label{eq:13-h1} \end{equation} $$

符号に注意してほしい．1–3 は陽子どうし（$+\ee$ と $+\ee$），2–4 は電子どうし（$-\ee$ と $-\ee$）なので積は $+\ee^2$ で斥力，1–4 と 2–3 は陽子と電子なので積は $-\ee^2$ で引力である．全ハミルトニアンは

$$ \Ham=\Ham_0+\Ham_1 $$

である．この $\Ham_1$ が $x_1$ と $x_2$ を結びつける（連成させる）ことが，以下のすべての鍵になる．

数学ノート：1次元調和振動子の厳密解（復習）

質量 $m$，バネ定数 $k$ の1次元調和振動子

$$ \hat{H}=\frac{\hat{p}^2}{2m}+\frac{1}{2}k\hat{x}^2 $$

のエネルギー固有値は，角振動数 $\omega=\sqrt{k/m}$ を用いて

$$ E_n=\left(n+\frac{1}{2}\right)\hbar\omega,\qquad n=0,1,2,\dots $$

で与えられる．$n=0$ でもエネルギーは $0$ にならず $\tfrac12\hbar\omega$ が残る．これを零点エネルギー（zero-point energy）と呼ぶ．古典力学では「静止していれば運動エネルギーも位置エネルギーも $0$」だが，量子力学では不確定性関係 $\Delta x\,\Delta p\ge\hbar/2$ のため，$x$ と $p$ を同時に $0$ にできない．本章の主役はまさにこの $\tfrac12\hbar\omega$ である．

13.3 $R\gg\abs{x_1},\abs{x_2}$ での展開 — 双極子–双極子相互作用の出現

13.3.1 まず $L$ を落とす

式 \eqref{eq:13-h1} は分母に $x_1,x_2$ が入っているので，このままでは扱えない．原子どうしが十分離れている状況，すなわち $R$ が結合長 $L$ よりずっと大きい場合を考える．$R\gg L$ ならば

$$ R+L+x_1\simeq R+x_1,\qquad R-(L+x_2)\simeq R-x_2 $$

としてよい．2–4 間の距離は $L$ がもともと相殺して $R+x_1-x_2$ であるから，結局

$$ \begin{equation} \Ham_1\simeq\frac{\ee^2}{4\pi\eps}\left[\frac{1}{R}+\frac{1}{R+x_1-x_2}-\frac{1}{R+x_1}-\frac{1}{R-x_2}\right] \label{eq:13-h1-nol} \end{equation} $$

となる．以下ではこの形から出発する．

注意：$L$ を落とすことの物理的意味

「$R\gg L$ だから $L$ を落とす」という操作は，数学的には $L/R$ の1次の項を捨てることに相当する．捨てられた項の中には，自然長 $L$ そのものがつくる永久双極子どうしの相互作用 $-2\ee^2L^2/(4\pi\eps R^3)$ が含まれている．しかし実際の中性原子では，電子雲の重心は時間平均すると核の位置にあり，永久双極子は存在しない（図13.1(a)）．$L$ は「バネの自然長」というモデル上の便宜であって，永久双極子を意味するものではない．物理的に意味があるのはつり合いからのゆらぎ $x_1,x_2$ だけであり，それゆえ $L$ を落とした式 \eqref{eq:13-h1-nol} が，我々が知りたい London 分散力を正しく記述する．永久双極子を持つ極性分子どうしの相互作用（Keesom力）を扱いたいときは，$L$ の項を残せばよい．

13.3.2 二項展開を実行する

$R$ に比べて変位が小さい，すなわち $R\gg\abs{x_1},\abs{x_2}$ とし，$\ee^2/(4\pi\eps R)$ をくくり出す．

$$ \Ham_1=\frac{\ee^2}{4\pi\eps R}\left[1+\frac{1}{1+\dfrac{x_1-x_2}{R}}-\frac{1}{1+\dfrac{x_1}{R}}-\frac{1}{1-\dfrac{x_2}{R}}\right] $$

数学ノート：$(1+u)^{-1}$ の展開

$\abs{u}\lt 1$ のとき，等比級数の和の公式から

$$ \frac{1}{1+u}=1-u+u^2-u^3+\cdots $$

が成り立つ．実際，右辺を $S$ とおくと $S-(-u)S=1$ すなわち $(1+u)S=1$ である．より一般には二項展開

$$ (1+u)^{\alpha}=1+\alpha u+\frac{\alpha(\alpha-1)}{2!}u^2+\cdots $$

で $\alpha=-1$ とおけば $1-u+u^2-\cdots$ が得られる．本節では $u$ の2次まで残す．2次まで残すのが本質であることは，すぐ後で分かる．

$u_1\equiv x_1/R$, $u_2\equiv x_2/R$ とおき，それぞれを2次まで展開する．

$$ \begin{align} \frac{1}{1+(u_1-u_2)}&=1-(u_1-u_2)+(u_1-u_2)^2+\cdots \notag\\ \frac{1}{1+u_1}&=1-u_1+u_1^2+\cdots \notag\\ \frac{1}{1-u_2}&=1+u_2+u_2^2+\cdots \notag \end{align} $$

これらを代入して，角括弧の中身を次数ごとに整理する．

$$ \begin{align} \left[\ \cdots\ \right]&=1+\bigl\{1-(u_1-u_2)+(u_1-u_2)^2\bigr\}-\bigl\{1-u_1+u_1^2\bigr\}-\bigl\{1+u_2+u_2^2\bigr\}\notag\\[4pt] &=\underbrace{(1+1-1-1)}_{0\text{次}}+\underbrace{\bigl\{-(u_1-u_2)+u_1-u_2\bigr\}}_{1\text{次}}+\underbrace{\bigl\{(u_1-u_2)^2-u_1^2-u_2^2\bigr\}}_{2\text{次}} .\notag \end{align} $$

各次数を個別に見よう．まず0次は

$$ 1+1-1-1=0 $$

で消える．これは当然で，$x_1=x_2=0$ のとき2つの中性原子の間には（距離 $R$ の点電荷対として見れば）正味の相互作用がないことを意味する．次に1次は

$$ -(u_1-u_2)+u_1-u_2=-u_1+u_2+u_1-u_2=0 $$

となり，これも完全に相殺する．最後に2次は，$(u_1-u_2)^2=u_1^2-2u_1u_2+u_2^2$ を使って

$$ (u_1^2-2u_1u_2+u_2^2)-u_1^2-u_2^2=-2u_1u_2 $$

となり，$u_1u_2$ の交差項だけが生き残る．したがって

$$ \Ham_1=\frac{\ee^2}{4\pi\eps R}\cdot\left(-\frac{2x_1x_2}{R^2}\right) $$

すなわち

$$ \begin{equation} \Ham_1=-\frac{1}{4\pi\eps}\frac{2\ee^2x_1x_2}{R^3}. \label{eq:13-dipole} \end{equation} $$

13.3.3 これは双極子–双極子相互作用である

式 \eqref{eq:13-dipole} の意味をはっきりさせよう．原子1（陽子1＋電子2）の瞬間双極子モーメントは，つり合いからのずれの分だけ考えれば $p_2=-\ee x_2$，原子2（陽子3＋電子4）のそれは $p_1=-\ee x_1$ である．すると \eqref{eq:13-dipole} は

$$ \Ham_1=-\frac{1}{4\pi\eps}\frac{2p_1p_2}{R^3} $$

と書ける．これは電磁気学で学ぶ，同一直線上に並んだ2つの双極子（縦配置）の相互作用エネルギーそのものである．距離の3乗に反比例することと，係数が $2$ であることが特徴である．

物理的意味：$1/R^3$ と係数 $2$ はどこから来るか

双極子モーメント $p_1$ が距離 $R$ の点につくる電場は，双極子軸方向（縦配置）で $E=\dfrac{1}{4\pi\eps}\dfrac{2p_1}{R^3}$ である．この電場の中に双極子 $p_2$ が同じ向きに置かれると，相互作用エネルギーは $U=-p_2E=-\dfrac{1}{4\pi\eps}\dfrac{2p_1p_2}{R^3}$ となる．係数 $2$ は「双極子の軸上では電場が横方向の2倍強い」ことに由来する．$1/R^3$ は，双極子が正負の電荷の対であるために，遠方で $1/R$ のクーロン項が相殺して1つ次数が落ち，さらにもう一方も双極子なのでもう1つ落ちる，と理解すればよい．1次の項がすべて相殺したのは，両方の原子が中性だからである．片方がイオンなら1次の項が残り，$-1/R^4$ 型のイオン–誘起双極子相互作用になる．

例題13.1　$\Ham_1$ の大きさを見積もる

$R=4$ Å，$x_1=x_2=0.1$ Å のとき，式 \eqref{eq:13-dipole} の値を eV 単位で求めよ．ただし $\dfrac{\ee^2}{4\pi\eps}=14.400\ \mathrm{eV}\,\text{Å}$ を使ってよい．

［解］式 \eqref{eq:13-dipole} に代入すると

$$ \Ham_1=-\frac{2\times 14.400\ \mathrm{eV\,\text{Å}}\times(0.1\ \text{Å})\times(0.1\ \text{Å})}{(4\ \text{Å})^3} =-\frac{2\times 14.400\times 0.01}{64}\ \mathrm{eV}=-4.5\times10^{-3}\ \mathrm{eV}. $$

約 $-4.5$ meV，モルあたりに直すと $-4.5\times10^{-3}\times96.485=-0.43$ kJ/mol である．表13.1の van der Waals 結合のエネルギー（0.3〜10 kJ/mol）と同じ桁になっており，モデルが妥当な大きさを与えていることが分かる．ただしこれは「たまたま $x_1$ と $x_2$ が同符号で $0.1$ Å ずれていた瞬間」の値であり，これを本当のエネルギー変化にするには，次節以降の量子力学的な扱いが必要である．

注意：ここで話を終えてはいけない

式 \eqref{eq:13-dipole} は $x_1x_2$ に比例するから，$x_1$ と $x_2$ の符号が同じなら引力，逆なら斥力である．もし2つの原子のゆらぎが無相関なら，$\braket{x_1x_2}=\braket{x_1}\braket{x_2}=0$ となって，平均のエネルギー変化はゼロになってしまう．「瞬間双極子ができるから引力」という高校的説明では，実はここで話が止まってしまうのである．本当の答えは，$\Ham_1$ が2つの振動子を連成させ，基底状態そのものを変えてしまうところにある．これを見るために，次節で座標変換を行う．

13.4 基準座標への変換（45度回転）

13.4.1 なぜ座標変換するのか

全ハミルトニアン

$$ \Ham=\frac{p_1^2}{2m}+\frac{1}{2}Cx_1^2+\frac{p_2^2}{2m}+\frac{1}{2}Cx_2^2-\frac{2\ee^2}{4\pi\eps R^3}x_1x_2 $$

は，$x_1x_2$ という交差項があるために「$x_1$ だけの部分」と「$x_2$ だけの部分」に分けられない．すなわち変数分離ができず，そのままでは Schrödinger 方程式が解けない．

ところが，この形は連成振り子や2原子分子の振動でおなじみの構造をしている．こういうときは，基準座標（normal coordinate，基準振動座標）に移れば必ず対角化できる．今の場合，2つの振動子は完全に対等（$C$ も $m$ も同じ）だから，答えは「和」と「差」である．

定義：基準座標（対称・反対称座標）

$$ \begin{equation} x_s=\frac{1}{\sqrt2}(x_1+x_2),\qquad x_a=\frac{1}{\sqrt2}(x_1-x_2) \label{eq:13-normal-coord} \end{equation} $$

運動量についても同じ変換を行う．

$$ \begin{equation} p_s=\frac{1}{\sqrt2}(p_1+p_2),\qquad p_a=\frac{1}{\sqrt2}(p_1-p_2) \label{eq:13-normal-mom} \end{equation} $$

逆変換は

$$ x_1=\frac{1}{\sqrt2}(x_s+x_a),\qquad x_2=\frac{1}{\sqrt2}(x_s-x_a),\qquad p_1=\frac{1}{\sqrt2}(p_s+p_a),\qquad p_2=\frac{1}{\sqrt2}(p_s-p_a) $$

である．添字 $s$ は symmetric（対称），$a$ は antisymmetric（反対称）を表す．

数学ノート：これは $(x_1,x_2)$ 平面の45度回転である

式 \eqref{eq:13-normal-coord} を行列で書けば

$$ \begin{pmatrix}x_s\\x_a\end{pmatrix} =\frac{1}{\sqrt2}\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix} =\begin{pmatrix}\cos\theta&\sin\theta\\ \sin\theta&-\cos\theta\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix},\qquad \theta=\frac{\pi}{4} $$

である．この行列は直交行列（転置＝逆行列）であり，$(x_1,x_2)$ 平面を45度傾けた新しい直交軸をとることに対応する．直交変換なので長さの2乗が保たれ，

$$ x_1^2+x_2^2=x_s^2+x_a^2,\qquad p_1^2+p_2^2=p_s^2+p_a^2 $$

が成り立つ．また，正準交換関係も保たれる．実際

$$ [\hat{x}_s,\hat{p}_s]=\frac{1}{2}\bigl([\hat{x}_1,\hat{p}_1]+[\hat{x}_1,\hat{p}_2]+[\hat{x}_2,\hat{p}_1]+[\hat{x}_2,\hat{p}_2]\bigr)=\frac{1}{2}(i\hbar+0+0+i\hbar)=i\hbar, $$ $$ [\hat{x}_s,\hat{p}_a]=\frac{1}{2}\bigl([\hat{x}_1,\hat{p}_1]-[\hat{x}_1,\hat{p}_2]+[\hat{x}_2,\hat{p}_1]-[\hat{x}_2,\hat{p}_2]\bigr)=\frac{1}{2}(i\hbar-i\hbar)=0 $$

となるので，$(x_s,p_s)$ と $(x_a,p_a)$ はそれぞれ独立な正準変数の組であり，それぞれを普通の調和振動子として量子化してよい．ここが議論の要である．

13.4.2 変換を実行する

まず $\Ham_0$ を書き換える．逆変換を代入して丁寧に展開すると

$$ \begin{align} p_1^2+p_2^2&=\frac{1}{2}(p_s+p_a)^2+\frac{1}{2}(p_s-p_a)^2\notag\\ &=\frac{1}{2}(p_s^2+2p_sp_a+p_a^2)+\frac{1}{2}(p_s^2-2p_sp_a+p_a^2)\notag\\ &=p_s^2+p_a^2 , \label{eq:13-p-inv} \end{align} $$

同様に

$$ \begin{align} x_1^2+x_2^2&=\frac{1}{2}(x_s+x_a)^2+\frac{1}{2}(x_s-x_a)^2=x_s^2+x_a^2 . \label{eq:13-x-inv} \end{align} $$

交差項が相殺するのが気持ちよい．したがって

$$ \begin{equation} \Ham_0=\frac{p_s^2}{2m}+\frac{1}{2}Cx_s^2+\frac{p_a^2}{2m}+\frac{1}{2}Cx_a^2 . \label{eq:13-h0-normal} \end{equation} $$

次に $\Ham_1$ である．積 $x_1x_2$ は

$$ x_1x_2=\frac{1}{\sqrt2}(x_s+x_a)\cdot\frac{1}{\sqrt2}(x_s-x_a)=\frac{1}{2}(x_s^2-x_a^2) $$

だから，式 \eqref{eq:13-dipole} より

$$ \begin{equation} \Ham_1=-\frac{2\ee^2}{4\pi\eps R^3}\cdot\frac{1}{2}(x_s^2-x_a^2) =-\frac{\ee^2}{4\pi\eps R^3}\left(x_s^2-x_a^2\right). \label{eq:13-h1-normal} \end{equation} $$

驚くべきことに，交差項だった $\Ham_1$ が，基準座標では2乗の項だけになった．つまりバネ定数の補正として吸収できる．両者を足すと

$$ \begin{equation} \Ham=\frac{p_s^2}{2m}+\frac{1}{2}\left(C-\frac{2\ee^2}{4\pi\eps R^3}\right)x_s^2 +\frac{p_a^2}{2m}+\frac{1}{2}\left(C+\frac{2\ee^2}{4\pi\eps R^3}\right)x_a^2 . \label{eq:13-h-normal} \end{equation} $$

ここで，$\tfrac12Cx_s^2-\dfrac{\ee^2}{4\pi\eps R^3}x_s^2=\tfrac12\left(C-\dfrac{2\ee^2}{4\pi\eps R^3}\right)x_s^2$ とまとめたことに注意してほしい（係数 $\tfrac12$ をくくり出すので $2$ が現れる）．

式 \eqref{eq:13-h-normal} は完全に独立な2つの調和振動子の和である．座標変換（45度回転）ひとつでハミルトニアンがここまで綺麗になり，調和振動子の厳密解がそのまま使えるようになった．

13.4.3 2つの基準モードの物理的意味

それぞれのモードが何を表すかを見ておこう．

対称モード xs（x₁ = x₂） + − + − 双極子が同じ向き（頭と尾が向き合う） 静電的に有利 → バネが柔らかくなる C − 2e₀²/(4πε₀R³)　→　ωs ＜ ω₀ 反対称モード xa（x₁ = −x₂） + − + − 双極子が逆向き（頭と頭が向き合う） 静電的に不利 → バネが硬くなる C + 2e₀²/(4πε₀R³)　→　ωa ＞ ω₀
図13.3 2つの基準モード．赤い矢印は各原子の瞬間双極子モーメントの向き（正電荷から見た変位の向き）を表す．対称モードでは双極子が一直線に揃うので静電的に安定化し，実効バネ定数が減る．反対称モードではその逆になる．

なぜ？：結局これは第11章の永年方程式と同じことをしている

第11章では，2つの同じ原子軌道 $\phi_A,\phi_B$ が共鳴積分 $\beta$ で結合すると，$\tfrac{1}{\sqrt2}(\phi_A+\phi_B)$（結合性）と $\tfrac{1}{\sqrt2}(\phi_A-\phi_B)$（反結合性）に分裂し，エネルギーが $\alpha\pm\beta$ になることを学んだ．本節でやっているのはまったく同じ構造の問題である．「同じものが2つあって，それらが相互作用でつながっている」とき，答えは必ず和と差であり，対称性が高い方が安定化し，低い方が不安定化する．第12章の対称性の言葉で言えば，$x_s$ と $x_a$ はそれぞれ，2原子系の反転操作に対する対称・反対称既約表現に属する基底である．違いは，第11章では分裂したエネルギー準位を電子が片方だけ占有できたのに対し，ここでは2つの振動子モードが両方とも基底状態にあることである．だからこそ，次節で見るように「相殺しきらずに残るわずかな差」が主役になる．

13.5 零点エネルギーの差が引力を生む

13.5.1 2つのモードの角振動数

式 \eqref{eq:13-h-normal} は，バネ定数がそれぞれ $C\mp\dfrac{2\ee^2}{4\pi\eps R^3}$ である2つの独立な調和振動子だから，角振動数は $\omega=\sqrt{k/m}$ より

$$ \begin{equation} \omega_s=\sqrt{\frac{1}{m}\left(C-\frac{2\ee^2}{4\pi\eps R^3}\right)},\qquad \omega_a=\sqrt{\frac{1}{m}\left(C+\frac{2\ee^2}{4\pi\eps R^3}\right)} . \label{eq:13-omega-sa} \end{equation} $$

ここで，無次元の小さな量

$$ \begin{equation} \epsilon\equiv\frac{2\ee^2}{4\pi\eps CR^3}\qquad(\abs{\epsilon}\ll1) \label{eq:13-eps-def} \end{equation} $$

を導入すると，$\omega_0=\sqrt{C/m}$（式 \eqref{eq:13-omega0}）を使って

$$ \omega_s=\omega_0\sqrt{1-\epsilon},\qquad \omega_a=\omega_0\sqrt{1+\epsilon} $$

と簡潔に書ける．$R$ が大きいほど $\epsilon$ は小さく，$R\to\infty$ で $\omega_s=\omega_a=\omega_0$ に戻る（相互作用なし）．

数学ノート：$(1+u)^{1/2}$ の展開

二項展開 $(1+u)^\alpha=1+\alpha u+\dfrac{\alpha(\alpha-1)}{2!}u^2+\cdots$ で $\alpha=\tfrac12$ とおくと

$$ \sqrt{1+u}=1+\frac{1}{2}u+\frac{\frac12\left(\frac12-1\right)}{2}u^2+\cdots =1+\frac{1}{2}u-\frac{1}{8}u^2+\cdots $$

を得る．$u\to-u$ とすれば

$$ \sqrt{1-u}=1-\frac{1}{2}u-\frac{1}{8}u^2+\cdots $$

である．2次の項の符号が両方とも負であることが，この節の結論を決める．

この公式を $\epsilon$ の2次まで適用すると

$$ \begin{align} \omega_s&=\omega_0\left(1-\frac{1}{2}\epsilon-\frac{1}{8}\epsilon^2+\cdots\right), \label{eq:13-omega-s-exp}\\ \omega_a&=\omega_0\left(1+\frac{1}{2}\epsilon-\frac{1}{8}\epsilon^2+\cdots\right). \label{eq:13-omega-a-exp} \end{align} $$

1次の項は符号が逆であるが，2次の項はどちらも同じ符号（負）である．ここが決定的である．

13.5.2 零点エネルギーの変化

2つの基準モードはそれぞれ独立な調和振動子だから，系の基底状態エネルギー（零点エネルギー）は

$$ \begin{equation} E_0=\frac{1}{2}\hbar\omega_s+\frac{1}{2}\hbar\omega_a \label{eq:13-e0-def} \end{equation} $$

である．式 \eqref{eq:13-omega-s-exp}, \eqref{eq:13-omega-a-exp} を代入して和をとると，$\epsilon$ の1次の項が相殺して

$$ \begin{align} \omega_s+\omega_a&=\omega_0\left[\left(1-\frac{\epsilon}{2}-\frac{\epsilon^2}{8}\right)+\left(1+\frac{\epsilon}{2}-\frac{\epsilon^2}{8}\right)\right]\notag\\ &=\omega_0\left[2-\frac{\epsilon^2}{4}\right] \label{eq:13-omega-sum} \end{align} $$

となる．したがって

$$ E_0=\frac{\hbar}{2}\left(\omega_s+\omega_a\right)=\frac{\hbar\omega_0}{2}\left(2-\frac{\epsilon^2}{4}\right)=\hbar\omega_0-\frac{\hbar\omega_0}{8}\epsilon^2 . $$

一方，クーロン相互作用がまったく存在しない（$R\to\infty$）ときの零点エネルギーは

$$ E_0^{(0)}=\frac{1}{2}\hbar\omega_0+\frac{1}{2}\hbar\omega_0=\hbar\omega_0 $$

である．両者の差が，クーロン相互作用によって生じるエネルギー変化 $\Delta U$ である．

$$ \begin{equation} \Delta U=E_0-E_0^{(0)}=-\frac{\hbar\omega_0}{8}\epsilon^2 =-\frac{\hbar\omega_0}{8}\left(\frac{2\ee^2}{4\pi\eps CR^3}\right)^2 =-\frac{A}{R^6}, \label{eq:13-deltaU} \end{equation} $$

ここで

$$ \begin{equation} A=\frac{\hbar\omega_0}{8}\left(\frac{2\ee^2}{4\pi\eps C}\right)^2=\frac{\hbar\omega_0}{2}\left(\frac{\ee^2}{4\pi\eps C}\right)^2\ \ (\gt 0) \label{eq:13-A-def} \end{equation} $$

である．2つの中性原子の間には，距離の6乗に反比例する引力ポテンシャルが働くことが示された．これが London 分散力である．

零点エネルギーの変化の図．相互作用がないときは2つの振動子の準位が同じ高さ（½ħω₀）にあり，相互作用があると対称モードの準位 ω_s は下がり反対称モードの準位 ω_a は上がるが，下がる幅のほうが大きく，2つの和 E₀ は ħω₀ より小さくなる．変化は ε = 0.5 として誇張して描いてある．
図13.4 零点エネルギーの変化．対称モードの角振動数は下がり，反対称モードは上がる．1次の効果ではこの2つがちょうど打ち消し合うが，平方根の凸性（2次の項が両方とも負）のために，下がる幅がわずかに大きい．その差 $-\hbar\omega_0\epsilon^2/8$ が van der Waals 引力である．

物理的意味：量子力学がなければ van der Waals 力はない

式 \eqref{eq:13-deltaU} には $\hbar$ が露わに含まれている．$\hbar\to0$（古典極限）とすれば $\Delta U\to0$ である．これは偶然ではない．古典力学では，絶対零度において2つの振動子は $x_1=x_2=0$，$p_1=p_2=0$ で完全に静止できる．そのとき $\Ham_1=-2\ee^2x_1x_2/(4\pi\eps R^3)=0$ となり，相互作用エネルギーはぴったりゼロで，引力は生じない．

量子力学では不確定性関係により電子を静止させることができず，必ず零点振動が残る．この「止まれないこと」が2つの原子の間で相関し，その相関の分だけエネルギーが下がる．すなわちvan der Waals 力（London分散力）は，純粋に量子力学的な起源をもつ力である．高校で習った「瞬間的に双極子ができる」の「瞬間的なゆらぎ」の正体は，まさにこの零点振動なのである．

注意：1次の項が消えることこそが本質

式 \eqref{eq:13-omega-sum} で $\epsilon$ の1次の項が相殺したことを軽く見てはいけない．もし1次の項が残っていれば，$\Delta U\propto\epsilon\propto1/R^3$ となり，van der Waals 力は $1/R^3$ の距離依存性を持つはずである．実際には1次が消えるので $\epsilon^2\propto1/R^6$ となる．2次の摂動でしか効かないという点は，第4章の摂動論で学んだ「2次のエネルギー補正は基底状態を必ず下げる」という定理と同じ構造である．実際，$\Ham_1$ を摂動として第4章の2次摂動論を適用しても，同じ $-A/R^6$ が得られる．

13.5.3 分極率との関係と数値的な検証

式 \eqref{eq:13-A-def} の $A$ は，モデルのパラメータ $C$ と $\omega_0$ で書かれていて分かりにくい．しかしこれは，測定可能な量である分極率（polarizability）で書き直せる．

このモデルの原子に一様電場 $E$ をかけると，電子には力 $-\ee E$ が働き，バネの復元力とつり合う位置は $Cx=-\ee E$ すなわち $x=-\ee E/C$ である．生じる双極子モーメントは $p=-\ee x=\ee^2E/C$ だから，分極率は $\alpha=\ee^2/C$ となる．体積の次元をもつ分極率体積 $\alpha'\equiv\alpha/(4\pi\eps)=\ee^2/(4\pi\eps C)$ を使えば，式 \eqref{eq:13-A-def} はきわめて簡潔になる．

$$ \begin{equation} A=\frac{\hbar\omega_0}{2}\,\alpha'^2,\qquad \Delta U=-\frac{\hbar\omega_0}{2}\frac{\alpha'^2}{R^6}. \label{eq:13-A-alpha} \end{equation} $$

すなわち「分極しやすい原子ほど強く引き合う」．ヨウ素や Xe のように電子雲がふわふわした大きな原子ほど van der Waals 力が強く，沸点が高いことは，この式で説明できる．

数学ノート：3次元にすると係数は $3/4$ になる（Londonの公式）

本節のモデルは電子が $x$ 軸方向にしか動けない1次元モデルだった．実際の原子では電子は3方向に動ける．$y$ 方向・$z$ 方向の変位（横配置の双極子）については，双極子–双極子相互作用の係数が $-2$ ではなく $+1$ になる（軸上と横方向で電場の強さが2倍違うため）．したがって横モードでは $\epsilon_\perp=\epsilon/2$ となり，その寄与は

$$ -\frac{\hbar\omega_0}{8}\left(\frac{\epsilon}{2}\right)^2\times 2\ (\text{2方向})=-\frac{\hbar\omega_0}{16}\epsilon^2 $$

である．縦モードの寄与 $-\dfrac{\hbar\omega_0}{8}\epsilon^2$ と足すと

$$ \Delta U=-\left(\frac{1}{8}+\frac{1}{16}\right)\hbar\omega_0\epsilon^2=-\frac{3}{16}\hbar\omega_0\epsilon^2 =-\frac{3}{4}\frac{\hbar\omega_0\alpha'^2}{R^6} $$

となり，Londonが1930年に導いた有名な公式（文献[6]）に一致する．本書の1次元モデルはその係数を $1/2$ と見積もる（縦モードのみ）．桁は正しく，係数だけが $2/3$ 倍という良い近似である．

例題13.2　Ar 二量体の結合エネルギーを見積もる

Ar 原子の分極率体積は $\alpha'=1.64\ \text{Å}^3$，特性励起エネルギー（第1イオン化エネルギーで代用）は $\hbar\omega_0\simeq15.8$ eV である．Ar$_2$ の平衡核間距離 $R=3.76$ Å における分散力のエネルギーを，Londonの公式（係数 $3/4$）で見積もり，実測のポテンシャル井戸の深さ $12.3$ meV と比べよ．

［解］$R^6=(3.76)^6$ を計算する．$3.76^2=14.14$，$3.76^3=53.16$ なので $R^6=(53.16)^2=2.826\times10^{3}\ \text{Å}^6$．したがって

$$ \Delta U=-\frac{3}{4}\times\frac{15.8\ \mathrm{eV}\times(1.64\ \text{Å}^3)^2}{2.826\times10^3\ \text{Å}^6} =-\frac{0.75\times15.8\times2.690}{2.826\times10^3}\ \mathrm{eV}=-1.13\times10^{-2}\ \mathrm{eV}. $$

すなわち約 $-11$ meV である．実測値 $-12.3$ meV と 10 % 程度の一致であり，たった1個のバネで原子を表す粗いモデルとしては驚くほど良い．モルあたりに直すと $-11\times10^{-3}\times96.485=-1.1$ kJ/mol となり，表13.1の値とも整合する．なお本書の1次元モデル（係数 $1/2$）では $-7.5$ meV となり，こちらもオーダーは合っている．

13.5.4 斥力と Lennard-Jones ポテンシャル

式 \eqref{eq:13-deltaU} は引力しか与えない．しかし原子どうしを無限に近づけられるはずがない．近づきすぎると，電子雲どうしが重なり，Pauli の排他律（第7章）によって激しい斥力が生じる．この斥力を加えないと，原子は $R\to0$ に落ち込んでしまう．

斥力の距離依存性は，本来は波動関数の重なりが指数関数的に減衰することから $\exp(-R/\rho)$ 型になる（13.8節で再登場する）．しかし分子動力学計算などでは，計算の便宜上

$$ \begin{equation} U(R)=-\frac{A}{R^6}+\frac{B}{R^{12}} =4\varepsilon\left[\left(\frac{\sigma}{R}\right)^{12}-\left(\frac{\sigma}{R}\right)^{6}\right] \label{eq:13-lj} \end{equation} $$

という形がよく使われる．これをLennard-Jonesポテンシャル（Lennard-Jones potential）と呼ぶ．$\varepsilon$ は井戸の深さ，$\sigma$ は $U=0$ となる距離である．

例題13.3　Lennard-Jonesポテンシャルの最小点

式 \eqref{eq:13-lj} の右辺が最小となる $R$ と，そのときの $U$ を求めよ．

［解］$U(R)=4\varepsilon\left(\sigma^{12}R^{-12}-\sigma^6R^{-6}\right)$ を $R$ で微分すると

$$ \frac{\dd U}{\dd R}=4\varepsilon\left(-12\sigma^{12}R^{-13}+6\sigma^6R^{-7}\right) =\frac{24\varepsilon}{R^7}\left(\sigma^6-\frac{2\sigma^{12}}{R^6}\right). $$

これが $0$ になるのは $R^6=2\sigma^6$，すなわち $R_{\min}=2^{1/6}\sigma\simeq1.122\,\sigma$ のときである．このとき $(\sigma/R)^6=1/2$ だから

$$ U(R_{\min})=4\varepsilon\left[\left(\frac12\right)^2-\frac12\right]=4\varepsilon\left(\frac14-\frac12\right)=-\varepsilon . $$

したがって井戸の深さはちょうど $\varepsilon$ である．Ar では $\varepsilon=12.3$ meV，$\sigma=3.35$ Å（したがって $R_{\min}=3.76$ Å）が標準的な値である．

Lennard-Jonesポテンシャルのグラフ．赤の破線の引力項 −A/R⁶ と青の破線の斥力項 +B/R¹² の和が緑の実線で，R = 2^(1/6)σ で最小値 −ε をとる．
図13.5 Lennard-Jonesポテンシャル．赤の破線が本章で導出した引力項 $-A/R^6$，青の破線が経験的な斥力項 $+B/R^{12}$，緑の実線がその和である．$R=2^{1/6}\sigma$ で深さ $\varepsilon$ の井戸をつくる．曲線は式 \eqref{eq:13-lj} を実際に計算して描いた．

注意：$R^{-12}$ に物理的根拠はない

引力項の $R^{-6}$ は本章で導出したとおり明確な物理的根拠を持つ．しかし斥力項の $R^{-12}$ には理論的根拠がない．単に「$R^{-6}$ の2乗だから計算機で速く計算できる」（$R^{-6}$ を一度計算すればその2乗で済む）という便宜的な理由で選ばれたものである．物理的には $\exp(-R/\rho)$ 型（Born–Mayer型）の方が正しく，精密な分子間ポテンシャルではそちらが使われる．イオン結晶を扱う13.8節では Born–Mayer 型を用いる．

13.6 イオン結晶の静電エネルギーとMadelung定数

13.6.1 もう一方の極端：完全な電子移動

ここまでは，電子雲がほとんど重ならない極端（van der Waals）を扱った．本節からは正反対の極端に移る．第11章で見たように，2つの原子の電気陰性度の差が大きいと，永年方程式の解は片方の原子軌道にほぼ完全に偏る．極限では電子が一方から他方へ完全に移り，系は陽イオンと陰イオンの集まりになる．NaCl はその代表例で，Na$^+$ と Cl$^-$ がそれぞれ閉殻構造をとって岩塩型構造に並ぶ．

こうなると話は一気に古典的になる．イオンを点電荷とみなし，Coulomb の法則で全静電エネルギーを足し上げればよい——ように見える．ところが，この「足し上げ」がまったく簡単ではない．ここに Madelung の問題がある．

NaCl型（岩塩型）構造 紫の丸：原点にとった陽イオン 陽イオン（Na⁺） 陰イオン（Cl⁻） 原点の陽イオンから見た近接殻 殻 距離 個数 N 符号 1R₀6− (異符号) 2√2 R₀12+ (同符号) 3√3 R₀8− (異符号) 42 R₀6+ (同符号) 5√5 R₀24− (異符号) 6√6 R₀24+ (同符号) 7√8 R₀12+ (同符号) 83 R₀30− (異符号) ※ √7 R₀ の殻は存在しない（7 は3個の平方数の和で表せない）．
図13.6 NaCl型構造と，原点にとった陽イオンから見た近接殻．距離が $\sqrt{n}R_0$ の殻に含まれるイオン数 $N$ は，整数の組 $(n_1,n_2,n_3)$ で $n_1^2+n_2^2+n_3^2=n$ となるものの個数に等しい．$n_1+n_2+n_3$ が奇数なら異符号，偶数なら同符号のイオンである．

13.6.2 静電エネルギーの和を書き下す

結晶中のある1個の陽イオンを原点にとる（原点はどこでもよいが，結晶の中央あたりの陽イオンを選ぶと考えやすい）．このイオンが他のすべてのイオンから受ける静電ポテンシャルエネルギーを足し上げる．同じ距離にあるイオンをひとまとめにして殻（shell）と呼び，$j$ 番目の殻には $N_j$ 個のイオンが距離 $r_j$ にあるとする．図13.6の例では $N_1=6$（最近接の Cl$^-$），$N_2=12$（第2近接の Na$^+$），$N_3=8$（第3近接の Cl$^-$）である．

陽イオンの価数を $z_+$，陰イオンの価数を $z_-$ とすると（いずれも符号を含む数．NaCl なら $z_+=+1$, $z_-=-1$），第1近接からの寄与は

$$ N_1\frac{z_+z_-\ee^2}{4\pi\eps}\frac{1}{r_1}, $$

第2近接からの寄与は（相手が陽イオンなので $z_+z_+$）

$$ N_2\frac{z_+z_+\ee^2}{4\pi\eps}\frac{1}{r_2}, $$

となり，これが遠方までどこまでも続く．$\abs{z_+}=\abs{z_-}$ の場合（アルカリハライドや MgO など，多くのイオン結晶がそうである）には $z_+z_+=-z_+z_-$ なので，すべての殻を $z_+z_-$ でくくり出すことができる．第 $j$ 殻が異符号のイオンからなるとき $\sigma_j=+1$，同符号のとき $\sigma_j=-1$ とおけば，

$$ \begin{equation} U_\mathrm{M}(R)=\sum_{j=1}^{\infty}\sigma_jN_j\,\frac{z_+z_-\ee^2}{4\pi\eps}\frac{1}{r_j} \label{eq:13-um-sum} \end{equation} $$

と書ける．これがMadelungポテンシャルエネルギー（Madelung potential energy）である．$U_\mathrm{M}$ は「イオン1個あたり」のエネルギーであることに注意してほしい．

13.6.3 Madelung定数の定義

ここで，殻の半径 $r_j$ を最近接距離 $R_0$ で測り直す．$r_j$ は $R_0$ の定数倍だから

$$ \begin{equation} r_j=c_jR_0\qquad(c_j\ \text{は無次元の幾何因子},\ c_1=1) \label{eq:13-cj} \end{equation} $$

と書ける．図13.6の場合，$c_1=1,\ c_2=\sqrt2,\ c_3=\sqrt3,\ c_4=2,\ c_5=\sqrt5,\dots$ である．これを式 \eqref{eq:13-um-sum} に代入すると，$R_0$ が和の外にくくり出せる．

$$ U_\mathrm{M}(R_0)=\sum_{j=1}^{\infty}\sigma_jN_j\frac{z_+z_-\ee^2}{4\pi\eps}\frac{1}{c_jR_0} =\left[\sum_{j=1}^{\infty}\frac{\sigma_jN_j}{c_j}\right]\frac{z_+z_-\ee^2}{4\pi\eps}\frac{1}{R_0} $$

角括弧の中身は，もはや物質の種類にも大きさにもよらず，結晶構造の幾何だけで決まる純粋な数である．これを Madelung 定数と呼ぶ．

定義：Madelung定数と Madelung エネルギー

$$ \begin{equation} A_\mathrm{M}\equiv\sum_{j=1}^{\infty}\frac{\sigma_jN_j}{c_j}, \qquad \sigma_j=\begin{cases}+1&\text{第 }j\text{ 殻が原点と異符号のイオンのとき}\\[2pt]-1&\text{第 }j\text{ 殻が原点と同符号のイオンのとき}\end{cases} \label{eq:13-am-def} \end{equation} $$

この $A_\mathrm{M}$ をMadelung定数（Madelung constant）と呼ぶ．$\abs{z_+}$, $\abs{z_-}$ をイオン価数の絶対値と書けば $z_+z_-=-\abs{z_+}\abs{z_-}$ なので，Madelung エネルギーは

$$ \begin{equation} U_\mathrm{M}(R_0)=-A_\mathrm{M}\frac{\abs{z_+}\abs{z_-}\ee^2}{4\pi\eps}\frac{1}{R_0} \label{eq:13-um-final} \end{equation} $$

と表される．$A_\mathrm{M}\gt 0$ なので $U_\mathrm{M}\lt 0$，すなわちイオン結晶は静電的に安定である．$A_\mathrm{M}$ は結晶構造だけで決まり，NaCl 型なら（Na でも K でも Cl でも I でも）つねに $1.7476$ である．

注意：Madelung定数の符号の置き方は本によって違う

講義スライドや多くの教科書では，$N_j$ に符号を含めてしまい（異符号の殻を正，同符号の殻を負とする），$A_\mathrm{M}=-\sum_jN_j/c_j$ と書いたり，逆に $A_\mathrm{M}=\sum_jN_j/c_j$ と書いたりする．マイナス符号を $A_\mathrm{M}$ の定義に含めるか，$U_\mathrm{M}$ の前に出すかの違いにすぎないが，両方にマイナスを付けてしまうと符号が反転して「イオン結晶は不安定」というおかしな結論になる．最終的に $U_\mathrm{M}\lt 0$（安定）になっているかを必ず確認すること．本書では式 \eqref{eq:13-am-def} のように $A_\mathrm{M}\gt 0$ と約束し，マイナスは式 \eqref{eq:13-um-final} の前に置く．これは Kittel [1] をはじめ多くの標準的教科書の流儀である．

また，$A_\mathrm{M}$ の値はどの長さを基準にとるかにも依存する．本書では最近接イオン間距離 $R_0$ を基準にとる．格子定数 $a$ を基準にした値を載せている本もあるので（NaCl型なら $a=2R_0$ なので値が2倍になる），表を引用するときは必ず定義を確認してほしい．

13.6.4 素朴に足すと収束しない

ではさっそく NaCl 型の $A_\mathrm{M}$ を計算してみよう．図13.6の表から，殻の順に足していくと

$$ A_\mathrm{M}=\frac{6}{1}-\frac{12}{\sqrt2}+\frac{8}{\sqrt3}-\frac{6}{2}+\frac{24}{\sqrt5}-\frac{24}{\sqrt6}-\frac{12}{\sqrt8}+\frac{30}{3}-\cdots $$

部分和を順に計算すると次のようになる．

表13.2 NaCl型構造のMadelung定数を，近接殻の順に素朴に足し上げたときの部分和．正解は $1.7476$ である．
殻まで12345678
部分和6.000−2.4852.134−0.8679.8670.069−4.1745.826

まったく収束する気配がない．これは計算ミスではなく，この和が絶対収束しないことの表れである．殻の個数 $N_j$ は距離の2乗に比例して増える（球面上のイオン数だから）のに対し，各項は距離の1乗に反比例して減るだけなので，$\sum\abs{\sigma_jN_j/c_j}$ は発散してしまう．符号が交代するおかげでかろうじて有限の値に落ち着くだけである（条件収束，conditional convergence）．

物理的意味：なぜ収束が悪いのか，どう直すのか

球殻ごとに足すやり方の問題は，殻を切ったところで電荷の中性が破れることにある．たとえば第1殻まで（6個の陰イオン）で打ち切ると，原点の $+1$ に対して $-6$ の電荷が余分に残り，遠方から見ると大きな正味電荷が残った状態になる．この人工的な電荷が $1/R$ という長距離のポテンシャルをつくるので，和が暴れるのである．

解決策は2つある．ひとつは Evjen 法で，立方体の表面にあるイオンの電荷を面上なら $1/2$，稜上なら $1/4$，頂点なら $1/8$ と分割して，各段階で電気的に中性なブロックを積み上げる方法である．もうひとつが Ewald 法（文献[7]）で，各点電荷を「鋭い点電荷 − なめらかなGauss電荷」と「なめらかなGauss電荷」に分解し，前者を実空間で，後者を逆格子空間（Fourier空間）で足す．実空間側も逆格子空間側も指数関数的に収束するので，実用的な精度が短時間で得られる．13.10節で使う pymatgen の EwaldSummation は，まさにこの Ewald 法を実装したものである．

13.7 1次元イオン結晶のMadelung定数

13.7.1 手で計算できる唯一の例

3次元では手計算ができないが，1次元なら最後まで手で解ける．陽イオンと陰イオンが間隔 $R_0$ で交互に無限に並んだ1次元結晶を考えよう．

原点を陽イオンにとる．原点から見ると，

ある．したがって式 \eqref{eq:13-um-sum} は

$$ \begin{align} U_\mathrm{M}(R_0)&=2\cdot\frac{(+1)(-1)\ee^2}{4\pi\eps}\frac{1}{1\cdot R_0} +2\cdot\frac{(+1)(+1)\ee^2}{4\pi\eps}\frac{1}{2R_0} +2\cdot\frac{(+1)(-1)\ee^2}{4\pi\eps}\frac{1}{3R_0}+\cdots\notag\\[4pt] &=2\left(1-\frac12+\frac13-\frac14+\cdots\right)\times\frac{(+1)(-1)\ee^2}{4\pi\eps}\frac{1}{R_0} \label{eq:13-1d-sum} \end{align} $$

となる．2行目では，$(+1)(-1)$ を共通因子としてくくり出したので，陽イオンの殻（第2, 4, …）の項の符号が反転して $+\tfrac12$, $+\tfrac14$ が $-$ に変わったことに注意してほしい．式 \eqref{eq:13-am-def} と見比べれば，丸括弧の中身の2倍がまさに Madelung 定数である．

$$ A_\mathrm{M}=2\left(1-\frac12+\frac13-\frac14+\cdots\right) $$

数学ノート：$\ln(1+x)$ のMaclaurin展開

等比級数の和

$$ \frac{1}{1+t}=1-t+t^2-t^3+\cdots\qquad(\abs{t}\lt1) $$

の両辺を $t=0$ から $t=x$ まで積分する．左辺は $\displaystyle\int_0^x\frac{\dd t}{1+t}=\bigl[\ln(1+t)\bigr]_0^x=\ln(1+x)$，右辺は項別に積分して

$$ \int_0^x(1-t+t^2-\cdots)\dd t=x-\frac{x^2}{2}+\frac{x^3}{3}-\frac{x^4}{4}+\cdots $$

となる．したがって

$$ \begin{equation} \ln(1+x)=x-\frac{x^2}{2}+\frac{x^3}{3}-\frac{x^4}{4}+\cdots \label{eq:13-ln-series} \end{equation} $$

である．この級数の収束半径は $1$ であるが，境界の $x=1$ では交代級数となり（Leibniz の判定法により）収束する．すなわち $x=1$ を代入してよく，

$$ \ln 2=1-\frac12+\frac13-\frac14+\cdots $$

が成り立つ．この級数を交代調和級数（alternating harmonic series）と呼ぶ．

式 \eqref{eq:13-ln-series} に $x=1$ を代入すれば，ただちに

$$ \begin{equation} A_\mathrm{M}^{(1\mathrm{D})}=2\left(1-\frac12+\frac13-\frac14+\cdots\right)=2\ln 2\simeq1.386 \label{eq:13-am-1d} \end{equation} $$

を得る．1次元イオン結晶の Madelung 定数が，自然対数という解析的な形で書けてしまうのは実に美しい．

13.7.2 収束は非常に遅い

ただし，この級数の収束はきわめて遅い．第 $n$ 項までの部分和を実際に計算すると次のようになる．

表13.3 交代調和級数による1次元Madelung定数の部分和．極限値は $2\ln2=1.386294$．
$n$12410501001000$\infty$
部分和2.00001.00001.16671.29131.36651.37631.38531.3863

1000項も足してようやく小数第2位が確定する程度である．誤差は $1/(2n)$ 程度でしか減らない．3次元ではこれよりさらに悪く（表13.2），Ewald 法のような特別な工夫が不可欠になる理由がここにある．

(a) 1次元イオン結晶 + − + − + − + 原点 R₀ 2R₀ 左右対称なので，どの殻もイオンは2個ずつ AM = 2(1 − 1/2 + 1/3 − 1/4 + …) = 2 ln 2 ≒ 1.386 交代調和級数の部分和 S_n の項数 n への依存．n = 1 から 30 まで，緑の折れ線が極限値 2 ln 2 = 1.3863（赤の破線）をまたいで振動しながら，ゆっくりと近づく．
図13.7 (a) 1次元イオン結晶と，その Madelung 定数 $2\ln2$．(b) 交代調和級数の部分和の振る舞い．項を1つ足すごとに極限値をまたいで振動し，収束が非常に遅い．曲線は実際に部分和を計算して描いた．

注意：条件収束級数は「足す順番」を変えてはいけない

交代調和級数は条件収束するが絶対収束しない（$\sum1/n$ は発散する）．Riemann の再配列定理によれば，条件収束する級数は，項の順序を入れ替えることで任意の値に収束させられる．たとえば $1+\tfrac13-\tfrac12+\tfrac15+\tfrac17-\tfrac14+\cdots$ と並べ替えると $\tfrac32\ln2$ に収束してしまう．したがって Madelung 和では「近い殻から順に，電気的に中性なまとまりで足す」という手順が本質的に重要である．順番を変えると答えが変わってしまうのだから，これは数学的な形式論ではなく，物理的に正しい答えを得るための必須条件である．

例題13.4　原点を陰イオンにとっても同じ答えになるか

1次元イオン結晶で，原点を陰イオンにとって $A_\mathrm{M}$ を計算し，式 \eqref{eq:13-am-1d} と一致することを確かめよ．

［解］原点を陰イオン（価数 $-1$）にとると，距離 $R_0$ には陽イオンが2個，距離 $2R_0$ には陰イオンが2個，… と並ぶ．第1殻の寄与は

$$ 2\cdot\frac{(-1)(+1)\ee^2}{4\pi\eps}\frac{1}{R_0}, $$

第2殻の寄与は $2\cdot\dfrac{(-1)(-1)\ee^2}{4\pi\eps}\dfrac{1}{2R_0}$ である．$(+1)(-1)=(-1)(+1)$ でくくり出すと，括弧の中身は前と同じ $2(1-\tfrac12+\tfrac13-\cdots)$ になる．よって $A_\mathrm{M}=2\ln2$ で一致する．

これは当然の結果である．$A_\mathrm{M}$ は式 \eqref{eq:13-am-def} のとおり「異符号か同符号か」だけで決まり，どちらを原点にとるかによらない．ただし，$\abs{z_+}\neq\abs{z_-}$ の物質（CaF$_2$ など）では陽イオン基準と陰イオン基準で $U_\mathrm{M}$ の値が異なるので，「$A_\mathrm{M}$ を何単位あたりで定義したか」に注意が必要である．

13.8 3次元のMadelung定数と凝集エネルギー

13.8.1 代表的な構造のMadelung定数

Ewald 法などで正しく計算した結果を表にまとめる．いずれも最近接イオン間距離 $R_0$ を基準にとった値である．

表13.4 代表的なイオン結晶構造のMadelung定数（最近接イオン間距離 $R_0$ を基準）．Kittel『固体物理学入門』[1] などによる．
構造型代表物質陽イオンの配位数$A_\mathrm{M}$
NaCl型（岩塩型）NaCl，KCl，MgO61.747565
CsCl型CsCl，CsBr81.762675
閃亜鉛鉱型（zinc blende）ZnS，GaAs41.638055
ウルツ鉱型（wurtzite）ZnO，GaN41.6413
蛍石型（fluorite）CaF$_2$，ZrO$_2$82.5194

注目すべきことが2つある．ひとつは，配位数が 4 から 8 まで大きく変わるのに，$A_\mathrm{M}$ は $1.64$ から $1.76$ までの狭い範囲にしか収まっていないことである（蛍石型は $\abs{z_+}\abs{z_-}=2$ で定義されているので別扱い）．配位数が増えれば近くのイオンは増えるが，その分だけ同符号のイオンも近づくので，正味の利得はほとんど変わらないのである．

もうひとつは，CsCl 型（1.7627）が NaCl 型（1.7476）よりわずかに大きいことである．それでも NaCl が岩塩型をとるのは，静電エネルギーだけで構造が決まるわけではないからである．イオン半径比によって決まる斥力（次項）や，分極の効果が効いてくる．

数学ノート：暗算に便利な定数 $\ee^2/(4\pi\eps)=14.400\ \mathrm{eV}\,\text{Å}$

SI単位系のままだと数値計算が面倒なので，次の値を覚えておくとよい．

$$ \frac{\ee^2}{4\pi\eps}=14.400\ \mathrm{eV}\cdot\text{Å} $$

すなわち「$1$ Å 離れた $+1$ と $-1$ の点電荷の静電エネルギーは $-14.4$ eV」である．実際，$\ee=1.602\times10^{-19}$ C，$1/(4\pi\eps)=8.988\times10^{9}\ \mathrm{N\,m^2/C^2}$ から

$$ \frac{\ee^2}{4\pi\eps}=(1.602\times10^{-19})^2\times8.988\times10^{9}=2.307\times10^{-28}\ \mathrm{J\,m} $$ $$ =\frac{2.307\times10^{-28}}{1.602\times10^{-19}}\ \mathrm{eV\,m}=1.440\times10^{-9}\ \mathrm{eV\,m}=14.40\ \mathrm{eV}\,\text{Å} $$

となる（$1\ \text{Å}=10^{-10}$ m）．ついでに，エネルギー換算 $1\ \mathrm{eV}=96.485\ \mathrm{kJ/mol}$，圧力換算 $1\ \mathrm{eV}/\text{Å}^3=160.22\ \mathrm{GPa}$ も覚えておくと本節の計算がすべて手でできる．

例題13.5　NaCl のMadelungエネルギー

NaCl の最近接 Na–Cl 距離は $R_0=2.82$ Å である．イオン対（Na$^+$ 1個と Cl$^-$ 1個）あたりの Madelung エネルギーを eV 単位で求めよ．

［解］式 \eqref{eq:13-um-final} に $A_\mathrm{M}=1.7476$，$\abs{z_+}=\abs{z_-}=1$ を代入する．

$$ U_\mathrm{M}=-1.7476\times\frac{14.400\ \mathrm{eV}\,\text{Å}}{2.82\ \text{Å}} =-1.7476\times5.106\ \mathrm{eV}=-8.92\ \mathrm{eV}. $$

すなわちイオン対あたり $-8.92$ eV である．モルあたりに直すと $-8.92\times96.485=-861$ kJ/mol となる．

［注意］式 \eqref{eq:13-um-final} は「イオン1個が感じるエネルギー」だが，$N$ 個のイオン対からなる結晶の全静電エネルギーを求めるときは，$2N$ 個のイオンそれぞれについて $U_\mathrm{M}$ を足すと各ペアを2回数えてしまうので $\tfrac12$ を掛ける必要がある．結果として全エネルギーは $\tfrac12\times2N\times U_\mathrm{M}=NU_\mathrm{M}$，つまりイオン対あたりちょうど $U_\mathrm{M}$ になる．式 \eqref{eq:13-um-final} は「イオン対あたりの値」と読んでよい，というのがこの計算の結論である．

13.8.2 斥力を加える：Born–Mayer型ポテンシャル

$U_\mathrm{M}\propto-1/R$ だけでは，$R\to0$ でエネルギーが $-\infty$ に落ち込んでしまい，結晶は潰れてしまう．実際にはイオンの閉殻電子雲どうしが重なると，Pauli の排他律による強い斥力が働く．この斥力は電子雲の重なりに比例し，原子軌道が指数関数的に減衰する（第5章）ので，距離とともに指数関数的に減る．そこで

$$ \lambda\exp\!\left(-\frac{R}{\rho}\right) $$

という形（Born–Mayer型斥力，Born–Mayer repulsion）を仮定する．$\lambda$ は強さ，$\rho$ は斥力の到達距離を表すパラメータで，実験値（体積弾性率など）から決める．斥力は近接イオンとの間にしか働かないので，最近接の $z$ 個（NaCl 型なら $z=6$）だけを数えればよい．イオン対あたりの全ポテンシャルエネルギーは

$$ \begin{equation} U(R)=-A_\mathrm{M}\frac{\abs{z_+}\abs{z_-}\ee^2}{4\pi\eps}\frac{1}{R}+z\lambda\exp\!\left(-\frac{R}{\rho}\right) \label{eq:13-born-mayer} \end{equation} $$

となる．以下，簡単のため $\abs{z_+}=\abs{z_-}=1$ とし，$K\equiv A_\mathrm{M}\dfrac{\ee^2}{4\pi\eps}$ と略記する．

13.8.3 平衡距離と凝集エネルギー

平衡距離 $R_0$ は $U(R)$ が最小になる点，すなわち $\dd U/\dd R=0$ を満たす点である．式 \eqref{eq:13-born-mayer} を微分すると

$$ \frac{\dd U}{\dd R}=\frac{K}{R^2}-\frac{z\lambda}{\rho}\exp\!\left(-\frac{R}{\rho}\right) $$

だから，$R=R_0$ で

$$ \begin{equation} z\lambda\exp\!\left(-\frac{R_0}{\rho}\right)=\frac{\rho K}{R_0^2} \label{eq:13-equil} \end{equation} $$

が成り立つ．この関係を使うと，平衡点でのエネルギーから未知パラメータ $\lambda$ を消去できる．実際，式 \eqref{eq:13-born-mayer} に $R=R_0$ を代入し，\eqref{eq:13-equil} を使うと

$$ \begin{align} U(R_0)&=-\frac{K}{R_0}+\frac{\rho K}{R_0^2} =-\frac{K}{R_0}\left(1-\frac{\rho}{R_0}\right)\notag \end{align} $$

すなわち

$$ \begin{equation} U(R_0)=-A_\mathrm{M}\frac{\ee^2}{4\pi\eps R_0}\left(1-\frac{\rho}{R_0}\right) \label{eq:13-ucoh} \end{equation} $$

を得る．実に美しい結果である．凝集エネルギーは Madelung エネルギーに $(1-\rho/R_0)$ を掛けただけで得られる．$\rho$ は典型的には $0.3$ Å 程度，$R_0$ は $2$〜$4$ Å 程度だから，$\rho/R_0\simeq0.1$，つまり斥力は Madelung エネルギーの安定化を1割ほど食いつぶすにすぎない．イオン結晶の凝集エネルギーの9割は Madelung エネルギーであると言ってよい．

例題13.6　NaCl の凝集エネルギー

NaCl について $R_0=2.820$ Å，$\rho=0.321$ Å である（Kittel [1] の表による）．式 \eqref{eq:13-ucoh} からイオン対あたりの凝集エネルギーを求め，実験値 $182.6$ kcal/mol と比較せよ．

［解］例題13.5 より $-K/R_0=-8.92$ eV である．補正因子は

$$ 1-\frac{\rho}{R_0}=1-\frac{0.321}{2.820}=1-0.1138=0.8862 . $$

したがって

$$ U(R_0)=-8.92\times0.8862=-7.91\ \mathrm{eV}\ (\text{イオン対あたり}). $$

モルあたりに直すと $-7.91\times96.485=-763$ kJ/mol である．実験値は $182.6\ \mathrm{kcal/mol}\times4.184=764$ kJ/mol（の安定化）だから，誤差 0.2 % 未満という驚異的な一致である．「点電荷を並べて足し算するだけ」というモデルが，これほど定量的に成功するのは，イオン結合が本当にほぼ純粋な静電相互作用であることを意味している．

13.8.4 体積弾性率まで出る

式 \eqref{eq:13-born-mayer} からは，エネルギーだけでなく硬さも計算できる．体積弾性率（bulk modulus）$B$ は，体積 $V$ に対するエネルギーの2階微分

$$ B=V\frac{\dd^2U_{\text{全}}}{\dd V^2}\bigg|_{V=V_0} $$

で定義される．NaCl 型構造では，立方格子定数が $a=2R$ で，その中に4組のイオン対が入る．したがってイオン対あたりの体積は $V/N=(2R)^3/4=2R^3$，$N$ 対では $V=2NR^3$ である．ここから

$$ \frac{\dd V}{\dd R}=6NR^2\ \Longrightarrow\ \frac{\dd}{\dd V}=\frac{1}{6NR^2}\frac{\dd}{\dd R} $$

となる．$U_{\text{全}}=NU(R)$ だから，まず1階微分は

$$ \frac{\dd U_{\text{全}}}{\dd V}=\frac{1}{6NR^2}\cdot N\frac{\dd U}{\dd R}=\frac{U'(R)}{6R^2} $$

となり，もう一度微分して

$$ \frac{\dd^2U_{\text{全}}}{\dd V^2}=\frac{1}{6NR^2}\frac{\dd}{\dd R}\left[\frac{U'(R)}{6R^2}\right] =\frac{1}{6NR^2}\left[\frac{U''(R)}{6R^2}-\frac{2U'(R)}{6R^3}\right] $$

を得る．平衡点では $U'(R_0)=0$ なので第2項が落ち，$V_0=2NR_0^3$ を掛けると

$$ \begin{equation} B=2NR_0^3\cdot\frac{U''(R_0)}{36NR_0^4}=\frac{U''(R_0)}{18R_0}. \label{eq:13-bulk} \end{equation} $$

あとは $U''$ を計算すればよい．式 \eqref{eq:13-born-mayer} をもう一度微分すると

$$ U''(R)=-\frac{2K}{R^3}+\frac{z\lambda}{\rho^2}\exp\!\left(-\frac{R}{\rho}\right) $$

で，平衡条件 \eqref{eq:13-equil} を代入して $\lambda$ を消去すると

$$ U''(R_0)=-\frac{2K}{R_0^3}+\frac{1}{\rho^2}\cdot\frac{\rho K}{R_0^2} =\frac{K}{R_0^2}\left(\frac{1}{\rho}-\frac{2}{R_0}\right) =\frac{K}{R_0^3}\left(\frac{R_0}{\rho}-2\right) $$

となる．よって

$$ \begin{equation} B=\frac{A_\mathrm{M}\ee^2}{4\pi\eps}\cdot\frac{1}{18R_0^4}\left(\frac{R_0}{\rho}-2\right). \label{eq:13-bulk-final} \end{equation} $$

例題13.7　NaCl の体積弾性率

式 \eqref{eq:13-bulk-final} に $A_\mathrm{M}=1.7476$，$R_0=2.820$ Å，$\rho=0.321$ Å を代入して $B$ を GPa 単位で求め，実験値 $24$ GPa と比較せよ．

［解］まず $K=A_\mathrm{M}\ee^2/(4\pi\eps)=1.7476\times14.400=25.165\ \mathrm{eV}\,\text{Å}$．次に

$$ \frac{R_0}{\rho}-2=\frac{2.820}{0.321}-2=8.785-2=6.785, $$ $$ R_0^4=(2.820)^4=(7.952)^2=63.24\ \text{Å}^4 . $$

したがって

$$ B=\frac{25.165\times6.785}{18\times63.24}\ \frac{\mathrm{eV}}{\text{Å}^3} =\frac{170.7}{1138.3}=0.1500\ \mathrm{eV}/\text{Å}^3 . $$

$1\ \mathrm{eV}/\text{Å}^3=160.22$ GPa を使うと

$$ B=0.1500\times160.22=24.0\ \mathrm{GPa}. $$

実験値 $2.40\times10^{10}\ \mathrm{N/m^2}=24.0$ GPa と一致する．（もっとも，Kittel の表の $\rho$ はもともと $B$ の実験値から決めたパラメータなので，これは「うまく合った」というより「$\rho$ の決め方が自己無撞着である」ことの確認である．しかし，こうして1つのモデルからエネルギーと弾性率が同時に出てくること自体が重要である．）

注意：体積弾性率の単位の混乱

古い教科書では体積弾性率を CGS 単位系の $\mathrm{dyn/cm^2}$ で表している．換算は $1\ \mathrm{dyn/cm^2}=0.1\ \mathrm{N/m^2}=0.1$ Pa なので，

$$ 1\times10^{11}\ \mathrm{dyn/cm^2}=1\times10^{10}\ \mathrm{N/m^2}=10\ \mathrm{GPa} $$

である．Kittel の表に「NaCl: 2.40（$\times10^{11}\ \mathrm{dyn/cm^2}$）」とあるのは $24.0$ GPa の意味であって，$2.4$ GPa ではない．単位の桁を1つ間違えると結論が変わるので，必ず換算を確認すること．ちなみに $24$ GPa は，鉄（$170$ GPa）やダイヤモンド（$443$ GPa）よりはずっと柔らかく，アルミニウム（$76$ GPa）よりも柔らかい．岩塩が爪で傷つけられるほど柔らかいことと整合する．

13.9 アルカリハライドの化学的傾向

13.9.1 データを並べる

イオン結晶の代表格であるアルカリハライド $AX$（$A=$ Li, Na, K, Rb; $X=$ F, Cl, Br, I）は，いずれも NaCl 型構造をとる（したがって $A_\mathrm{M}$ はすべて $1.7476$ で共通である）．違うのは最近接距離 $R_0$ だけだから，式 \eqref{eq:13-ucoh} と \eqref{eq:13-bulk-final} が正しければ，すべての性質が $R_0$ の関数として整理できるはずである．

表13.5 NaCl型構造をもつアルカリハライドの諸性質．$R_0$：最近接距離，$\rho$：斥力の到達距離，$B$：体積弾性率，$U_{\exp}$：実験による凝集エネルギー（自由イオンを基準）．すべて Kittel『固体物理学入門』[1] の表（原データは M. P. Tosi [8]）による．$U_{\text{計算}}$ は本書の式 \eqref{eq:13-ucoh} で計算した値である．
$R_0$ (Å)$\rho$ (Å)$B$ (GPa)$U_{\exp}$ (kJ/mol)$U_{\text{計算}}$ (kJ/mol)
LiF2.0140.29167.110141031
LiCl2.5700.33029.8832824
LiBr2.7510.34023.8794773
LiI3.0000.36617.1744711
NaF2.3170.29046.5897917
NaCl2.8200.32124.0764763
NaBr2.9890.32819.9726723
NaI3.2370.34515.1683670
KF2.6740.29830.5794807
KCl3.1470.32617.4694692
KBr3.2980.33614.8663661
KI3.5330.34811.7627619
RbF2.8150.30126.2759770
RbCl3.2910.32315.6667665
RbBr3.4450.33813.0639636
RbI3.6710.34810.6606599

計算値と実験値は，どの物質でも数 % 以内で一致している．これは驚くべきことである．$A$ が Li から Rb へ，$X$ が F から I へ変わっても，必要な情報は $R_0$ と $\rho$ の2つだけなのだ．

13.9.2 散布図が語ること

アルカリハライド16種について，最近接距離 R₀ に対する体積弾性率 B（左，緑の点）と凝集エネルギーの絶対値 |U|（右，赤の点）の散布図．右の破線は |U| ∝ 1/R₀ の曲線（ρ/R₀ = 0.11 で固定）で，データ点はこの曲線の近くに並ぶ．
図13.8 表13.5のデータの散布図．(a) 最近接距離 $R_0$ に対する体積弾性率 $B$．(b) $R_0$ に対する凝集エネルギーの絶対値．破線は式 \eqref{eq:13-ucoh} が予言する $\abs{U}\propto1/R_0$ の曲線（$\rho/R_0=0.11$ で固定）で，データ点がよく乗っていることが分かる．プロットはすべて表13.5の実データによる．

図13.8から次のことが読み取れる．

物理的意味：バラバラにするのが大変な固体は，硬い

凝集エネルギーは「結晶を自由なイオンにバラバラにするのに必要なエネルギー」，体積弾性率は「押し縮めにくさ」である．一見別の量だが，どちらも同じポテンシャル曲線 $U(R)$ から出てくる．凝集エネルギーは井戸の深さ，体積弾性率は井戸底の曲率である．深い井戸は幅も狭い（$1/R$ 型のポテンシャルはそういう性質を持つ）ので，深さと曲率は自動的に相関する．

「バラバラにするのが大変な固体は硬い」というのは，日常的な直感とも一致する．ダイヤモンドは凝集エネルギーが大きく（$7.4$ eV/原子），体積弾性率も最大級（$443$ GPa）である．逆に固体アルゴンは凝集エネルギーが $0.08$ eV/原子しかなく，体積弾性率も $1$ GPa 程度しかない．本章の前半で扱った van der Waals 結晶が柔らかいのは，まさにこのためである．

例題13.8　MgO はなぜ硬いのか

MgO は NaCl 型構造をとり，$R_0=2.11$ Å，$\abs{z_+}=\abs{z_-}=2$ である．Madelung エネルギーを求め，NaCl と比べよ．

［解］式 \eqref{eq:13-um-final} より

$$ U_\mathrm{M}=-1.7476\times\frac{2\times2\times14.400}{2.11}\ \mathrm{eV}=-1.7476\times\frac{57.60}{2.11}=-47.7\ \mathrm{eV}. $$

NaCl の $-8.92$ eV の実に 5.3 倍である．内訳は，価数の積による因子 $4$ と，距離が短いことによる因子 $2.82/2.11=1.34$ の積（$4\times1.34=5.35$）である．MgO の体積弾性率は $160$ GPa と NaCl の7倍近くあり，融点も $2852\ ^\circ$C（NaCl は $801\ ^\circ$C）と非常に高い．価数を1つ上げる効果は絶大であることが分かる．耐火物やセラミックスに2価・3価のイオン結晶が使われる理由がここにある．

13.10 pymatgenでMadelungエネルギーを計算する

13.10.1 Materials Genome Initiative と pymatgen

2011年，米国は Materials Genome Initiative（MGI，材料ゲノム計画）という国家プロジェクトを立ち上げた．「新材料の開発期間と費用を半分にする」ことを目標に，実験と計算とデータ科学を統合しようという構想である．ここで開発された Python モジュールが pymatgen（Python Materials Genomics，文献[9]）であり，結晶構造の生成・変換・解析から，第一原理計算の入出力処理，相図の作成まで，材料計算に必要な道具がひととおり揃っている．中心メンバーは G. Ceder，K. A. Persson，S. P. Ong，A. Jain らである．

pymatgen と対になるのが Materials Project（文献[10]）で，十数万物質の第一原理計算結果（構造，生成エネルギー，バンド構造，状態密度など）が無料で公開されている．マテリアルズ・インフォマティクスを始めるなら，まず pymatgen で Materials Project のデータをいじってみるのが定石である．生成AIのない時代にこれだけのソフトウェア基盤を人力で作り上げたのは，率直に言って驚異的である．

本節では，Materials Project から NaCl の構造情報（POSCAR 形式）をダウンロードし，pymatgen の EwaldSummation クラスで Madelung エネルギーと Madelung 定数を計算する．第10章で PySCF を使ったときと同様，Google Colab で実行してもローカルの Ubuntu で実行してもよい．

!pip3 install pymatgen

コード13.1 pymatgen のインストール（Colab の場合）．エラー ModuleNotFoundError: No module named 'pymatgen' が出たらこれを実行する．

13.10.2 計算スクリプト

# calc_madelung_constant.py
import argparse
from pymatgen.core import Structure
from pymatgen.analysis.ewald import EwaldSummation

KE = 14.399645  # e0^2 / (4 pi eps0) [eV * Angstrom]

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--poscar", required=True, help="POSCAR ファイル名")
    args = parser.parse_args()

    # 1) 結晶構造を読み込む
    st = Structure.from_file(args.poscar)
    print(st)
    print("POSCAR:", args.poscar)

    # 2) 形式電荷(酸化数)を推定して各サイトに与える
    #    Ewald 和は点電荷の集まりに対する計算なので，電荷が必要
    st.add_oxidation_state_by_guess()
    print("Oxidation states: guessed")

    # 3) Ewald 法で静電エネルギーを求める
    ewald = EwaldSummation(st)
    e_cell = ewald.total_energy          # セルあたりの値 [eV]

    reduced, z = st.composition.get_reduced_composition_and_factor()
    print("Formula:", st.composition.formula)
    print("Reduced formula:", reduced.formula)
    print("Formula units per cell (Z):", int(z))
    print("Madelung energy (cell):", e_cell, "eV")
    print("Madelung energy (per formula unit):", e_cell / z, "eV")

    # 4) 基準にする cation-anion ペアから Madelung 定数を逆算する
    cations = [i for i, s in enumerate(st) if s.specie.oxi_state > 0]
    anions  = [i for i, s in enumerate(st) if s.specie.oxi_state < 0]
    i0 = cations[0]
    j0 = min(anions, key=lambda j: st.get_distance(i0, j))
    d_min = st.get_distance(i0, j0)
    zz = abs(st[i0].specie.oxi_state * st[j0].specie.oxi_state)
    e_pair = -KE * zz / d_min            # 1 ペア分の静電エネルギー [eV]
    a_madelung = (e_cell / z) / e_pair   # 比をとれば Madelung 定数

    print()
    print("Reference cation-anion pairs:")
    print(f"  {st[i0].specie} - {st[j0].specie}: d_min = {d_min:.6f} A,"
          f" E_pair = {e_pair:.9f} eV, |z+ z-| = {zz:.0f},"
          f" M_eff = {a_madelung:.9f}")
    print()
    print(f"Madelung constant: {a_madelung:.9f}"
          f" (reference pair: {st[i0].specie} - {st[j0].specie})")

if __name__ == "__main__":
    main()

コード13.2 POSCAR を読み込んで Madelung エネルギーと Madelung 定数を計算するスクリプト．

!python3 calc_madelung_constant.py --poscar POSCAR_NaCl_Fm-3m

コード13.3 実行方法．

13.10.3 出力の読み方

Full Formula (Na4 Cl4)
Reduced Formula: NaCl
abc   :   5.691694   5.691694   5.691694
angles:  90.000000  90.000000  90.000000
Sites (8)
  #  SP       a    b    c
---  ----  ----  ---  ---
  0  Na+   0     0    0
  1  Na+   0.5   0.5  0
  2  Na+   0.5   0    0.5
  3  Na+   0     0.5  0.5
  4  Cl-   0     0    0.5
  5  Cl-   0.5   0    0
  6  Cl-   0     0.5  0
  7  Cl-   0.5   0.5  0.5
POSCAR: POSCAR_NaCl_Fm-3m
Oxidation states: guessed
Formula: Na4 Cl4
Reduced formula: NaCl
Formula units per cell (Z): 4
Madelung energy (cell): -35.369871389878 eV
Madelung energy (per formula unit): -8.842467847469 eV

Reference cation-anion pairs:
  Na+ - Cl-: d_min = 2.845847 A, E_pair = -5.059880404206 eV, |z+ z-| = 1, M_eff = 1.747564594633

Madelung constant: 1.747564594633 (reference pair: Na+ - Cl-)

コード13.4 NaCl に対する出力例．

出力を上から確認していこう．

注意：「セルあたり」と「組成式あたり」を取り違えない

計算結果を論文や実験値と比べるとき，もっとも多い間違いが単位の取り違えである．「Madelung エネルギーは $-35.4$ eV だ」と言ってしまうと，それは NaCl 4組ぶんの値であって，教科書の $-8.9$ eV と4倍ずれる．第一原理計算のソフトウェアはほとんどが「セルあたり」の全エネルギーを出力するので，必ず Z（セル中の組成式の数）で割ってから比較する習慣をつけてほしい．逆に，実験の凝集エネルギーは通常「モルあたり（＝組成式あたり）」で与えられるので，$\times96.485$ で kJ/mol に直せば直接比較できる．

例題13.9　出力から手計算で検算する

コード13.4の出力を使って，Madelung 定数が $1.7476$ になることを手で確かめよ．

［解］最近接距離は $d_{\min}=2.845847$ Å である．1組の Na$^+$–Cl$^-$ ペアの静電エネルギーは

$$ E_{\text{pair}}=-\frac{\ee^2}{4\pi\eps d_{\min}}=-\frac{14.3996\ \mathrm{eV}\,\text{Å}}{2.845847\ \text{Å}}=-5.0599\ \mathrm{eV}, $$

これは出力の E_pair = -5.059880 eV と一致する．一方，組成式あたりの Madelung エネルギーは $-8.842468$ eV だから，その比は

$$ A_\mathrm{M}=\frac{-8.842468}{-5.059880}=1.747565 . $$

表13.4の NaCl 型の値と完全に一致する．Madelung 定数とは「最近接1ペアぶんの静電エネルギーの何倍か」を表す数であると理解すればよい．NaCl では最近接の Cl$^-$ が6個あるのだから素朴には $6$ 倍になりそうなものだが，第2近接の Na$^+$ からの斥力などが効いて $1.75$ 倍まで減ってしまう．それでも符号は負のままで，結晶は安定である．

13.11 研究例：Madelungエネルギーが駆動する現象

Madelung エネルギーは，教科書の中だけの量ではない．第一原理計算で得られた全エネルギーの変化を「なぜそうなったのか」と解釈するとき，全エネルギーを静電項（Madelung 項）とそれ以外に分解して比べる，という手法が実際の研究で使われている．本節では，担当教員が関わった3つの研究例を紹介する．

13.11.1 逆ペロブスカイトの構造相転移

ペロブスカイト $ABX_3$ は，$B$ を中心とする $BX_6$ 八面体が頂点を共有してつながり，その隙間に $A$ が入った構造である．通常の酸化物ペロブスカイト（SrTiO$_3$ など）では $A$, $B$ が陽イオン，$X$ が陰イオンだが，逆ペロブスカイト（antiperovskite，inverse perovskite）$A_3BX$ ではこれが逆転し，陽イオンと陰イオンの役割が入れ替わっている（例：Ba$_3$PN，Ca$_3$AsN）．

ペロブスカイトでは，イオン半径の比が理想からずれると $BX_6$ 八面体が傾いて（八面体回転ひずみ），立方晶 $Pm\bar{3}m$ から直方晶 $Pbnm$ へ構造相転移する．文献[11]では，$A_3BX$ 型逆ペロブスカイト8種類について第一原理計算を行い，相転移エネルギー $E_{Pbnm}-E_{Pm\bar{3}m}$ と，その Madelung 項だけを取り出した $E^{\text{Madelung}}_{Pbnm}-E^{\text{Madelung}}_{Pm\bar{3}m}$ を比べた．結果，この2つの量がきれいに相関することが分かった．

八面体回転ひずみによる構造相転移 立方晶 Pm-3m → 直方晶 Pbnm 八面体が交互に傾き，イオン間距離が最適化される → Madelungエネルギーが下がる 相転移エネルギーとの相関（模式） Madelungエネルギー差 (meV/f.u.) ΔE (meV/f.u.) −1600 −800 0 −600 −400 −200 0 Ba₃PN Sr₃PN Ca₃PN Ba₃AsN A₃PN A₃AsN
図13.9 逆ペロブスカイト $A_3BX$ の立方晶 $Pm\bar{3}m$ から直方晶 $Pbnm$ への構造相転移（左）と，相転移による全エネルギーの低下量と Madelung エネルギーの低下量の相関（右，文献[11]の結果をもとにした模式図）．八面体回転が大きい物質ほど両者がともに大きく下がる．

この結果が面白いのは，酸化物ペロブスカイトではこの相関が成り立たないことである．SrTiO$_3$ や CaTiO$_3$ のバルクの八面体回転は，Madelung エネルギーの減少ではなく，共有結合性の変化（Ti–O の混成）によって駆動されると考えられている．逆ペロブスカイトは，より「イオン性が高い固体」であるがゆえに，静電エネルギーが素直に主役を務めるのだろう，というのが文献[11]の解釈である．

13.11.2 ペロブスカイト酸化物の表面再構成

ところが，同じ酸化物ペロブスカイトでも表面では話が変わる．文献[12]では，$A^{2+}B^{4+}$O$_3$，$A^{+}B^{5+}$O$_3$，$A^{3+}B^{3+}$O$_3$ の3系列について，(001) 表面がとりうる再構成（cation exchange，checker，stripe，zigzag の4種類）を系統的に第一原理計算で調べた．その結果，どの系列でも表面エネルギーと「Madelung 表面エネルギー」が正の相関を示すことが分かった．すなわち，静電エネルギーを下げる再構成ほど，実際に安定な表面構造になっている．

物理的意味：バルクでは脇役，表面では主役

同じ物質でも，バルクの構造相転移では Madelung エネルギーが主因ではないのに，表面再構成では主因になる．この対比は示唆的である．バルクでは，どの原子も対称的な環境にあり，Madelung ポテンシャルは全体としてよく釣り合っている．そこから少し構造をひずませても静電エネルギーの変化は小さく，共有結合性のような「細かい」効果が勝つ．

一方，表面では対称性が根本から破れ，イオンが片側にしか隣人を持たない．この静電的なアンバランスは非常に大きく，それを解消する方向に原子が並び替わるのが表面再構成である．つまり，対称性が破れているところほど静電エネルギーがものを言う．表面，界面，粒界，点欠陥のまわり——材料科学の面白い現象が起きる場所は，たいていここである．

13.11.3 表面ランプリング：BaO の例

もうひとつの例が表面ランプリング（surface rumpling）である．イオン結晶の (001) 表面では，最表面の陽イオンと陰イオンが同じ平面に留まらず，垂直方向に少しずれる．一般には，表面の陽イオンは配位数が減ってしまうため，失った配位を稼ごうとしてバルク側にめり込む（陰イオンより内側に沈む）と説明される．

ところが文献[13]の第一原理計算によれば，BaO ではこの常識が破れ，Ba が真空側に飛び出す（逆向きのランプリング）ことが予測された．BaO の Ba$^{2+}$ はイオン半径が大きく配位数も高いはずなのに，である．この逆転もまた，Madelung エネルギーが減少する方向にランプリングが起きる，として説明された．

(a) 一般的なイオン結晶（例：BeO, MgO） 真空 理想表面の位置 + − + − 陽イオンがバルク側にめり込む （失った配位数を稼ぐため） (b) BaO：逆向きのランプリング 真空 + − + − Ba が真空側に飛び出す （Madelungエネルギーが下がるため）
図13.10 表面ランプリングの模式図．(a) 一般には表面の陽イオンがバルク側にめり込む．(b) BaO では逆に Ba が真空側に変位すると第一原理計算で予測された（文献[13]）．この逆転は Madelung エネルギーの減少で説明される．

この3つの例に共通するのは，「第一原理計算で得られた全エネルギーの差」という，それだけでは解釈しにくい数値を，Madelung エネルギーという古典的で分かりやすい量に投影して理解するという戦略である．本章で学んだ点電荷の足し算は，最先端の計算材料科学でも現役の道具なのである．

13.12 まとめと演習

13.12.1 この章のまとめ

  1. van der Waals力の起源．中性・無極性の原子どうしが引き合うのは，電子の量子的なゆらぎによって電気的つり合いが瞬間的に破れ，生じた瞬間双極子が隣の原子に双極子を誘起するからである．この成分を London 分散力と呼ぶ．
  2. 連成振動子モデル．核にバネでつながれた電子を2組並べ，Born–Oppenheimer 近似で核を固定する．$\Ham_0$ は独立な2つの調和振動子（式 \eqref{eq:13-h0}），$\Ham_1$ は4つの粒子対のクーロン相互作用（式 \eqref{eq:13-h1}）である．
  3. 展開．$R\gg\abs{x_1},\abs{x_2}$ で $(1+u)^{-1}=1-u+u^2-\cdots$ を使って展開すると，0次と1次の項がすべて相殺し，2次の交差項だけが残って $\Ham_1=-\dfrac{2\ee^2x_1x_2}{4\pi\eps R^3}$（式 \eqref{eq:13-dipole}）となる．これは縦配置の双極子–双極子相互作用である．
  4. 基準座標．45度回転 $x_s=(x_1+x_2)/\sqrt2$, $x_a=(x_1-x_2)/\sqrt2$ により，ハミルトニアンは実効バネ定数 $C\mp\dfrac{2\ee^2}{4\pi\eps R^3}$ をもつ2つの独立な調和振動子に分離する（式 \eqref{eq:13-h-normal}）．対称モードは柔らかく，反対称モードは硬くなる．
  5. 零点エネルギー．$\sqrt{1\pm\epsilon}$ を2次まで展開して和をとると，1次の項が相殺し，2次の項だけが残る．その結果 $\Delta U=-\dfrac{\hbar\omega_0}{8}\epsilon^2=-\dfrac{A}{R^6}$（式 \eqref{eq:13-deltaU}）．$\hbar$ が消えないので，量子力学の零点振動がなければ van der Waals 力は存在しない．斥力を加えると Lennard-Jones ポテンシャル（式 \eqref{eq:13-lj}）になる．
  6. Madelungエネルギー．イオン結晶の静電エネルギーは $U_\mathrm{M}=-A_\mathrm{M}\dfrac{\abs{z_+}\abs{z_-}\ee^2}{4\pi\eps R_0}$（式 \eqref{eq:13-um-final}）と書ける．$A_\mathrm{M}=\sum_j\sigma_jN_j/c_j$ は結晶構造だけで決まる純粋に幾何学的な数である．
  7. 1次元の厳密解．$A_\mathrm{M}^{(1\mathrm{D})}=2\ln2\simeq1.386$（式 \eqref{eq:13-am-1d}）．交代調和級数は条件収束であり，収束が非常に遅く順番も変えられない．3次元では Ewald 法が必須である．
  8. 凝集エネルギーと硬さ．Born–Mayer 斥力を加えると，$U(R_0)=-A_\mathrm{M}\dfrac{\ee^2}{4\pi\eps R_0}\left(1-\dfrac{\rho}{R_0}\right)$（式 \eqref{eq:13-ucoh}），$B=\dfrac{A_\mathrm{M}\ee^2}{4\pi\eps}\dfrac{1}{18R_0^4}\left(\dfrac{R_0}{\rho}-2\right)$（式 \eqref{eq:13-bulk-final}）．NaCl では凝集エネルギー $763$ kJ/mol（実験 $764$），$B=24.0$ GPa（実験 $24.0$）という驚異的な一致が得られる．
  9. 化学的傾向．$R_0$ が小さいほど Madelung エネルギーの絶対値が大きく，凝集エネルギーも体積弾性率も大きい．「バラバラにするのが大変な固体は硬い」．
  10. 実際の道具として．pymatgen の EwaldSummation で任意の結晶の Madelung エネルギーが計算できる．研究では，第一原理計算の全エネルギー差を Madelung 項に投影して現象を解釈する手法が使われている（逆ペロブスカイトの構造相転移，ペロブスカイト表面の再構成とランプリング）．

次章（第14章）では，イオン結晶を「点電荷の集まり」ではなく「電子の波」の側から見直す．本章で計算した Madelung ポテンシャルが，実は各イオンの軌道エネルギーを上下にずらす働きをしており，それがイオン結晶のバンドギャップを決めていることが分かる．本章の $A_\mathrm{M}$ は，第14章でそのまま再登場する．

13.12.2 演習問題

演習13.1　展開の1次項が消えることの確認

式 \eqref{eq:13-h1-nol} を $u_1=x_1/R$, $u_2=x_2/R$ の3次まで展開し，3次の項を求めよ．さらに，$\Ham_1$ の展開で0次と1次の項が消える理由を，「両方の原子が中性である」という事実から説明せよ．

ヒント：$(1+u)^{-1}=1-u+u^2-u^3+\cdots$ を使う．3次の項は $-(u_1-u_2)^3+u_1^3-u_2^3$ を計算すればよい．$(u_1-u_2)^3=u_1^3-3u_1^2u_2+3u_1u_2^2-u_2^3$ を展開して整理すると，$u_1^2u_2$ と $u_1u_2^2$ の項だけが残るはずである．中性であることの意味については，片方をイオン（正味電荷 $+\ee$，すなわち電子を取り去った状態）に置き換えたときに0次・1次の項がどうなるかを考えてみるとよい．

演習13.2　零点エネルギーの符号

(1) 式 \eqref{eq:13-omega-sum} で $\epsilon$ の1次の項が相殺することを，$\sqrt{1+u}$ の展開を書き下して確かめよ．
(2) 一般に，$f(u)=\sqrt{1+u}$ について $f(\epsilon)+f(-\epsilon)\lt 2f(0)$ が成り立つことを，$f$ が上に凸（$f''\lt0$）であることから説明せよ．この事実が「van der Waals力は必ず引力である」ことを保証していることを述べよ．
(3) $\Delta U$ が $\hbar$ に比例することから，古典極限で van der Waals 力が消えることを説明せよ．

ヒント：(2) は Jensen の不等式そのものである．上に凸な関数では，2点での値の平均は中点での値より小さい．$\tfrac12[f(\epsilon)+f(-\epsilon)]\lt f(0)$ を図示してみるとよい．(3) は「絶対零度で古典的な振動子はどこにいるか，そのとき $\Ham_1$ はいくらか」を考える．

演習13.3　2次元正方イオン格子のMadelung定数

陽イオンと陰イオンが間隔 $R_0$ で交互に並ぶ2次元正方格子を考える．原点の陽イオンから見て，第1近接（4個の陰イオン，距離 $R_0$），第2近接（4個の陽イオン，距離 $\sqrt2R_0$），第3近接（4個の陰イオン，距離 $2R_0$），第4近接（8個の陽イオン，距離 $\sqrt5R_0$）までの寄与を足し，部分和を順に求めよ．収束しているといえるか．

ヒント：$A_\mathrm{M}$ の部分和は $\dfrac41-\dfrac{4}{\sqrt2}+\dfrac42-\dfrac{8}{\sqrt5}+\cdots$ である．数値を入れて順に足すこと．$(n_1,n_2)$ で表される格子点までの距離は $\sqrt{n_1^2+n_2^2}\,R_0$，符号は $n_1+n_2$ の偶奇で決まる．第4近接が8個になるのは $(2,1),(1,2),(-1,2),\dots$ の8通りがあるからである．正しい極限値は $A_\mathrm{M}^{(2\mathrm{D})}=1.6155$ であり，表13.2と同様に部分和は大きく振動する．

演習13.4　KBr の凝集エネルギーと体積弾性率

表13.5から KBr の $R_0=3.298$ Å，$\rho=0.336$ Å を読み取り，(1) Madelung エネルギー，(2) 式 \eqref{eq:13-ucoh} による凝集エネルギー（eV/イオン対 と kJ/mol の両方），(3) 式 \eqref{eq:13-bulk-final} による体積弾性率（GPa）を計算し，表の実験値と比較せよ．

ヒント：$\ee^2/(4\pi\eps)=14.400\ \mathrm{eV}\,\text{Å}$，$A_\mathrm{M}=1.7476$，$1\ \mathrm{eV}=96.485$ kJ/mol，$1\ \mathrm{eV}/\text{Å}^3=160.22$ GPa を使う．例題13.5〜13.7と同じ手順を踏めばよい．$R_0^4=(3.298)^4$ の計算では，まず $3.298^2=10.877$ を求めてからその2乗をとると楽である．

演習13.5　なぜ NaCl は CsCl 型をとらないのか

表13.4によれば CsCl 型の Madelung 定数（1.7627）は NaCl 型（1.7476）より約 0.9 % 大きい．それにもかかわらず NaCl が NaCl 型構造をとる理由を，次の観点から論じよ．
(1) CsCl 型では最近接の配位数が 6 から 8 に増える．式 \eqref{eq:13-born-mayer} の斥力項はどう変わり，平衡距離 $R_0$ はどうなるか．
(2) 陽イオンと陰イオンの半径比 $r_+/r_-$ が小さいとき，8配位は幾何学的に可能か．

ヒント：(1) 斥力項は $z\lambda\exp(-R/\rho)$ で，$z$ は配位数である．配位数が増えれば斥力も増えるので平衡距離 $R_0$ が伸びる．式 \eqref{eq:13-ucoh} は $1/R_0$ に比例するから，$A_\mathrm{M}$ の 0.9 % の利得は $R_0$ が 0.9 % 伸びるだけで帳消しになってしまう．(2) 8配位（立方体の頂点に陰イオン，中心に陽イオン）が安定に存在できる半径比の下限は，立方体の体対角線を考えると $r_+/r_-=\sqrt3-1=0.732$ である．Na$^+$/Cl$^-$ の半径比は $1.02/1.81=0.56$ であり，この条件を満たさない．

演習13.6　pymatgen で MgO を計算する

Materials Project から MgO（NaCl 型，$a\simeq4.26$ Å）の POSCAR をダウンロードし，コード13.2 を使って Madelung エネルギーと Madelung 定数を計算せよ．得られた「組成式あたりの Madelung エネルギー」を，式 \eqref{eq:13-um-final} による手計算の値と比べよ．さらに，この値が例題13.8で求めた値とどう違うか，その理由を述べよ．

ヒント：MgO は $\abs{z_+}=\abs{z_-}=2$ であることに注意する．$R_0=a/2$．Madelung 定数そのものは構造が同じなので NaCl と同じ 1.7476 になるはずである．例題13.8では実験値 $R_0=2.11$ Å を使ったが，Materials Project の計算値は実験値と数 % 違うことがある．その差がエネルギーにどう伝わるかを考えるとよい．もし add_oxidation_state_by_guess() がうまく推定しない場合は，st.add_oxidation_state_by_element({"Mg": 2, "O": -2}) のように明示的に与えればよい．

13.12.3 参考文献

  1. C. Kittel『固体物理学入門（上）』第8版，宇野良清ほか訳，丸善出版（2005）第3章「固体の結合エネルギー」．本章の表13.4・表13.5のデータ，および13.2〜13.5節の導出の流儀はこの本による．
  2. 原田義也『量子化学（上巻）』裳華房（2007）．調和振動子の厳密解と，分子間力の量子化学的な扱い．
  3. P. A. Cox『固体の電子構造と化学』魚崎浩平ほか訳，技報堂出版（1989）（原著：The Electronic Structure and Chemistry of Solids, Oxford University Press, 1987）．イオン結晶のMadelungポテンシャルと電子状態の関係（第14章の準備）．
  4. シュライバー・アトキンス『無機化学（上）』第6版，田中勝久ほか訳，東京化学同人（2016）．イオン半径比則と格子エンタルピー（Born–Haberサイクル）．
  5. P. Atkins, J. de Paula, Atkins' Physical Chemistry, 11th ed., Oxford University Press (2018), 分子間相互作用の章．London分散力・Keesom力・Debye力の分類．
  6. F. London, "Zur Theorie und Systematik der Molekularkräfte", Zeitschrift für Physik 63, 245 (1930). 分散力の原論文．式 \eqref{eq:13-A-alpha} の係数 $3/4$ はここに由来する．
  7. P. P. Ewald, "Die Berechnung optischer und elektrostatischer Gitterpotentiale", Annalen der Physik 369, 253 (1921). Ewald法の原論文．
  8. M. P. Tosi, "Cohesion of Ionic Solids in the Born Model", Solid State Physics 16, 1 (1964). 表13.5の原データ．
  9. S. P. Ong et al., "Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis", Computational Materials Science 68, 314 (2013).
  10. A. Jain et al., "Commentary: The Materials Project: A materials genome approach to accelerating materials innovation", APL Materials 1, 011002 (2013).
  11. Y. Mochizuki et al., Physical Review Materials 4, 044601 (2020). 逆ペロブスカイトの構造相転移とMadelungエネルギー（13.11.1節）．
  12. Y. Mochizuki et al., Chemistry of Materials 35, 2047 (2023). ペロブスカイト酸化物の表面再構成とMadelung表面エネルギー（13.11.2節）．
  13. Y. Mochizuki et al., Physical Review Materials 2, 124603 (2018). BaOの逆向き表面ランプリング（13.11.3節）．