この記事では,この記事や,この記事でこれまで解説してきた密度汎関数理論(Density Functional Theory, DFT)に,これまでの議論では欠けていた(特殊)相対性理論による補正を取り込む方法について述べたいと思います。なお,この記事はGPT-5 ThinkingなどのLLMを積極的に援用しています。LLMが苦手な方はご注意ください。
序論 なぜ(特殊)相対性理論による補正が必要なのか
重元素では相対論的効果が無視できません。例えば,水素様ランタン($La^{56+}$)は,ボーア模型の見積もり$ v/c \approx Z\alpha $により,光速の約42%で運動します($Zα \simeq 57/137 \approx 0.416 $)。この速度域では,相対論的量子力学(Dirac方程式)に基づく補正が不可欠です。実際,金(Au)が金色を示すこと,水銀(Hg)が常温常圧で液体であること,白金(Pt)の高い触媒活性などは,6s軌道の相対論的収縮や5d準位の相対論的変化に起因する効果としてよく説明されます。さらに,ランタノイド収縮は主として4f電子の遮蔽の弱さ(有効核電荷の増大)に起因しますが,相対論的効果も副次的に寄与します。加えて,軽元素であっても微細構造(質量速度補正(運動エネルギー補正)項,Darwin項,スピン-軌道相互作用(Spin Orbit Coupling; SOC)など)を議論する際には,相対論的補正が欠かせません。したがって,第一原理計算・量子化学計算でも,対象や要求精度に応じて,Kohn–Sham方程式に相対論的効果を組み込む必要があります。以下では,まず球対称ポテンシャルにおけるDirac方程式を概説します。
球対称ポテンシャルにおけるDirac方程式
相対論的効果は内殻の深いところで生じるので,Dirac方程式は球対称ポテンシャルを仮定して解くことができます。内殻領域では球対称近似がよく成り立ち,その結果は擬ポテンシャル(後述)生成などに広く活用されます。Dirac方程式の復習から始めましょう。Diracは1928年に,Schrödinger方程式を相対論的共変形式に一般化した式
i \hbar \frac{\partial}{\partial t} \Psi=\left(c \boldsymbol{\alpha} \cdot \mathbf{\hat{p}}+\beta m c^2+V\right) \Psi=H \Psi \quad(1)
を提案しました。ここで$\Psi$は4個の成分(スピノル)をもつ1粒子波動関数で,スピン1/2を持つ粒子を表し,$V$は外部ポテンシャルです。$V = -e^2/(4\pi\epsilon_0 r)$である水素原子の場合には,Dirac方程式は解析的に解け,束縛状態のエネルギー固有値を厳密に求めることができます(例えば,adharaさんの記事,相対論的水素原子に潜むN = 2超対称性などを参照)。$\mathbf{\hat{p}}=-i\hbar \boldsymbol{\nabla}$は通常の運動量演算子であり,$(4 \times 4)$の行列 $\alpha_i$と$\beta$はPauli行列を使って,
\alpha_i=\left(\begin{array}{ll}
0 & \sigma_i \\
\sigma_i & 0
\end{array}\right), \quad \beta=\left(\begin{array}{rr}
\mathbf{1} & 0 \\
0 & \mathbf{- 1}
\end{array}\right) \quad(2)
と書けます。ここで$\sigma_i$は$(2 \times 2)$のPauliスピン行列で,
\sigma_1=\left(\begin{array}{ll}
0 & 1 \\
1 & 0
\end{array}\right), \quad \sigma_2=\left(\begin{array}{rr}
0 & -i \\
i & 0
\end{array}\right), \quad \sigma_3=\left(\begin{array}{rr}
1 & 0 \\
0 & -1
\end{array}\right) \quad(3)
であり,$\beta$の中の単位要素は$(2 \times 2)$の単位行列です。この方程式の解は次式
\Psi\left(x^\mu\right)=\mathrm{e}^{-\mathrm{i} \varepsilon t / \hbar}\binom{\phi(\mathbf{r})}{\chi(\mathbf{r})} \quad(4)
のように書くと便利です。ここで$\phi(\mathbf{r})$と$\chi(\mathbf{r})$は時間に依存しない2成分のスピノルで,空間とスピン1/2の自由度を表します。このようにすれば,Diracの方程式は$\phi$と$\chi$についての連立方程式
\begin{aligned}
& c(\sigma \cdot \mathbf{\hat{p}}) \chi=\left(\varepsilon-V-m c^2\right) \phi \quad(5a) \\
& c(\sigma \cdot \mathbf{\hat{p}}) \phi=\left(\varepsilon-V+m c^2\right) \chi \quad(5b)
\end{aligned}
となります。電子(正のエネルギーを持つ解)については$\phi$は大きく,$\chi$は相対的に小さく(概ね$v/c \sim |\mathbf{p}|/(mc)$ のオーダー)なります。$\phi$をLarge成分と呼び電子に対する波動関数,$\chi$をSmall成分と呼び陽電子運動との相互作用が電子運動に対する影響を指すと説明されることがありますが,実際には両成分ともその混合であり,混合の度合いは重元素になるほど増します。外部ポテンシャル$V(r)$が球対称の場合には,パリティと量子数$j m$によって表される全角運動量の保存則を使うことができます。各主量子数$n$に対する波動関数は動径関数と角度-スピン関数(非相対論での球面調和関数に対応)を使って書くことができ,
\psi_{n j m}^l=\binom{g_{n j}(r) \varphi_{j m}^l}{-i f_{n j}(r) \frac{\sigma \cdot \mathbf{r}}{r} \varphi_{j m}^l} \quad(6)
と表せます。この式は同じ$j m$を持ちますが,2個のとり得る値$j=l \pm \frac{1}{2}$に対する反対のパリティを持つ2つの関数を表しています。2個の成分を持つ関数 $\varphi_{j m}^l$ の詳しい形は下記のようになります。
\begin{aligned}
& j=l+\frac{1}{2}: \\
& \varphi_{j m}^l=\sqrt{\frac{l+\frac{1}{2}+m}{2 l+1}} Y_l^{m-\frac{1}{2}} \chi_{\frac{1}{2}+\frac{1}{2}}+\sqrt{\frac{l+\frac{1}{2}-m}{2 l+1}} Y_l^{m+\frac{1}{2}} \chi_{\frac{1}{2}-\frac{1}{2}} \quad(7a) \\
& j=l-\frac{1}{2}: \\
& \varphi_{j m}^l= \sqrt{\frac{l+\frac{1}{2}-m}{2 l+1}} Y_l^{m-\frac{1}{2}} \chi_{\frac{1}{2}+\frac{1}{2}}-\sqrt{\frac{l+\frac{1}{2}+m}{2 l+1}} Y_l^{m+\frac{1}{2}} \chi_{\frac{1}{2}-\frac{1}{2}} \quad(7b)
\end{aligned}
ここで$Y_l^{m}$は球面調和関数です。ここで,動径関数に対する式は,エネルギーを
\varepsilon^{\prime}=\varepsilon-m c^2 \quad(8)
とし,動径方向に変化する質量を
M(r)=m+\frac{\varepsilon^{\prime}-V(r)}{2 c^2} \quad(9)
とし,量子数$\kappa$を
\kappa= \mp\left(j+\frac{1}{2}\right) \quad\left\{\begin{array}{l}
j=l+\frac{1}{2} \text { のときは }\Rightarrow \kappa=-(l+1), \\
j=l-\frac{1}{2} \text { のときは }\Rightarrow \kappa=+l
\end{array}\right\} \quad(10)
とおけば簡単になります。ここで $\kappa(\kappa+1)=l(l+1)$ が両方の場合について成り立つことに注意しましょう。このようにすれば連立方程式は動径方程式の形で書くことができ,
\begin{aligned}
& -\frac{\hbar^2}{2 M} \frac{1}{r^2} \frac{\mathrm{~d}}{\mathrm{~d} r}\left(r^2 \frac{\mathrm{~d} g_{n \kappa}}{\mathrm{~d} r}\right)+\left[V+\frac{\hbar^2}{2 M} \frac{l(l+1)}{r^2}\right] g_{n \kappa} \\
& -\frac{\hbar^2}{4 M^2 c^2} \frac{\mathrm{~d} V}{\mathrm{~d} r} \frac{\mathrm{~d} g_{n \kappa}}{\mathrm{~d} r}-\frac{\hbar^2}{4 M^2 c^2} \frac{\mathrm{~d} V}{\mathrm{~d} r} \frac{(1+\kappa)}{r} g_{n \kappa}=\varepsilon^{\prime} g_{n \kappa}
\end{aligned} \quad(11)
また,
\frac{\mathrm{d} f_{n \kappa}}{\mathrm{~d} r}=\frac{1}{\hbar c}\left(V-\varepsilon^{\prime}\right) g_{n \kappa}+\frac{(\kappa-1)}{r} f_{n \kappa} \quad(12)
となります。これらの方程式は球対称ポテンシャルに対する一般形であり,これまでのところどのような近似も使っていません。すると(11)式は質量$M$が動径の関数であり,左辺に2つの付加項があるということ以外は普通のSchrödinger方程式と同じです。ところでこの付加項はそれぞれDarwin項とSOCの項です。前者は原子核の有効ポテンシャルを変え,これは電子と原子核の静電相互作用が,電子のジグザグ運動(Zitterbewegung)や高速量子振動によって乱されていると解釈することができます。後者は次の関係式
\frac{2}{\hbar^2} \mathbf{L} \cdot \mathbf{S} \varphi_{\kappa m}=-(1+\kappa) \varphi_{\kappa m} \quad(13)
とスピン演算子$\mathbf{S}(=\hbar \boldsymbol{\sigma} / 2)$を陽に使って書くことができます(ここで,$\mathbf{L}$は軌道角運動量演算子です)。(13)式中の$\varphi_{\kappa m}$は(10)式に示したように,$\kappa$によって(7a)式および(7b)式の適切な$\varphi_{j m}^l$に対応します。
もしSOC項が小さいと仮定すると,$g$と$f$の動径方程式ではそれを省略でき,摂動として扱うことができます。そうすれば(11)式と(12)式は主量子数$n$と軌道角運動量$l$のみに依存することになり,近似関数$\tilde{g}_{n l}$と$\tilde{f}_{n l}$を使って,
-\frac{\hbar^2}{2 M} \frac{1}{r^2} \frac{\mathrm{~d}}{\mathrm{~d} r}\left(r^2 \frac{\mathrm{~d} \tilde{g}_{n l}}{\mathrm{~d} r}\right)+\left[V+\frac{\hbar^2}{2 M} \frac{l(l+1)}{r^2}\right] \tilde{g}_{n l}-\frac{\hbar^2}{4 M^2 c^2} \frac{\mathrm{~d} V}{\mathrm{~d} r} \frac{\mathrm{~d} \tilde{g}_{n l}}{\mathrm{~d} r}=\varepsilon^{\prime} \tilde{g}_{n l} \quad(14)
と
\tilde{f}_{n l}=\frac{\hbar}{2 M c} \frac{\mathrm{~d} \tilde{g}_{n l}}{\mathrm{~d} r} \quad(15)
のように書くことができ,これに規格化条件
\int\left(\tilde{g}_{n l}^2+\tilde{f}_{n l}^2\right) r^2 \mathrm{~d} r=1 \quad(16)
が付きます。(14)式はスカラーの相対論的動径方程式であり,通常の非相対論方程式の場合と同じ方法で解くことができます。その後で他の方程式は動径方向のグリッド上で簡単に扱うことができます。最後にSOC項はMacDonaldら[1]に従って扱うことができます。(13)式と共に波動関数の大きい成分をつなぐスピン-軌道Hamiltonianは
\hat{H}_{\mathrm{SO}}=\frac{1}{2 M^2 c^2} \frac{1}{r} \frac{d V}{d r} \mathbf{L} \cdot \mathbf{S} \quad(17)
で与えられ,この項は小さい摂動として扱うことができる場合も多いです。この項は$\frac{1}{r} \frac{\mathrm{~d} V}{\mathrm{~d} r}$が大きくなる原子核に近い深い内殻で生じるので,今考えているSOC項の導出は,原子だけでなく固体や分子にそのまま適用できます。ここで,相対論的補正を明示的な演算子として解釈するにはFoldy–Wouthuysen (FW)変換をとるのが通例なので,こちらも紹介しておきます。球対称ポテンシャル$V(\mathbf{r})$(外部磁場なし)を仮定すると,(5a)式および(5b)式のFW変換による$1 / c^2$展開は,
\left[\frac{\hat{\mathbf{p}}^2}{2 m}+V-\frac{\hat{\mathbf{p}}^4}{8 m^3 c^2}+\frac{1}{2 m^2 c^2} \frac{1}{r} \frac{d V}{d r} \mathbf{L} \cdot \mathbf{S}+\frac{\hbar^2}{8 m^2 c^2} \nabla^2 V\right] \phi=\epsilon_k \phi \quad(18)
と書き下せます。(18)式の左辺カギ括弧の中は一電子Breit–Pauli Hamiltonianと呼ばれ,Large成分の2成分スピノルしか必要としません。左辺カギ括弧内の第3項は速度が質量に与える効果である質量速度補正項,第4項はすでに紹介したSOC項((17)式と微妙に一致しないのは,$M(r)$に対して非相対論的な展開を行っていないためです),第5項もすでに紹介したDarwin項です。ここで,(11)式で明示的に質量速度補正項が現れなかったのは,$M(r)$に対して非相対論的な展開を行っていないためであり,$M(r)$に対してこのような展開を行えば,質量速度補正項がきちんと現れます。
相対論的密度汎関数理論
ここまでは,球対称ポテンシャルを仮定した原子系における相対論的補正について述べてきました。しかし,多電子系を第一原理計算・量子化学計算で扱うには,DFTを相対論的に拡張する必要があります。本節ではその概念とその際に現れるDirac-Kohn-Sham方程式について述べます。
相対論的電磁気学にならい,外部電磁場を4成分ポテンシャル
A^\mu(\mathbf{r})=(\phi(\mathbf{r}), \mathbf{A}(\mathbf{r})) \quad(19)
と書きます。ここで$\phi$は非相対論的な外部静電ポテンシャル,$\mathbf{A}(\mathbf{r})$はベクトルポテンシャルです。電子の4成分電流密度は
j^\mu(\mathbf{r})=(cn(\mathbf{r}), \mathbf{j}(\mathbf{r})) \quad(20)
と書きます。ここで,$n(\mathbf{r})$は電子密度,$\mathbf{j}(\mathbf{r})$は電流密度です。ここで,相対論的なHohenberg–Kohn(HK)の定理であるRajagopal-Callawayの定理は,非縮退基底状態に対して,次の2つが成り立つとします。
- 第1定理:ゲージの自由度を除いて,4成分外部ポテンシャル(核電子静電ポテンシャルやベクトルポテンシャルなど1電子ポテンシャル)$A^\mu$は4成分電流密度$j^\mu$で決定される。
- 第2定理:表現可能な4成分電流密度$j^\mu$について,常にエネルギーの変分原理が成り立ち,エネルギー汎関数は下に有界で,基底状態で極小値をとる。
第1定理は,非相対論的なHK定理と同様に,相対論的なHamiltonian演算子が4成分電流密度$j^\mu$で一意的に表わせることを意味します。この定理は,同じ4成分電流密度$j^\mu$に対応する異なる2つの4成分外部ポテンシャル$A^\mu$が存在していると仮定し,その仮定が4成分スピノル波動関数に対する変分原理と矛盾することを背理法で示すことで証明されています。第1定理について気をつけたいのは,外部ポテンシャルは4成分ポテンシャル$A^\mu=(\phi, \mathbf{A})$で記述しますが,ここで磁場は $\mathbf{B}=\nabla \times \mathbf{A}$ で決まるため,
\mathbf{A} \rightarrow \mathbf{A}^{\prime}=\mathbf{A}+\nabla \Lambda \quad(21)
のようなゲージ変換をしても$\mathbf{B}$は変わらないという点です(定常状態問題では外部静電ポテンシャル$\phi$の加法定数は物理量に影響しません)。このため,見かけ上ちがう4成分ポテンシャル$A^\mu$と$A^{\prime \mu}$ が,同じ磁場$\mathbf{B}$をつくることがあります。すると,4成分ポテンシャルの違い(ゲージの違い)が基底状態の4成分スピノル波動関数$\Psi$にまで同じ影響を与えるとは限らず, $A^\mu$と$\Psi$の一対一対応は保証できません。しかし,Rajagopal-Callawayの定理が必要としている一意性は,$A^\mu$と$\Psi$の対応ではなく,4成分電流密度$j^\mu$と4成分外部ポテンシャル$A^\mu$の対応(ゲージの自由度を除く)です。定理の変分原理も,基本変数としての4成分電流密度$j^\mu= (cn, \mathbf{j})$を用いて書かれており,$A^\mu \leftrightarrow \Psi$の厳密な一対一対応は一般には要請されません。ただし,外部磁場$\mathbf{B} = \mathbf{0}$の定常状態問題では,非縮退かつ適切なゲージ固定のもとで,実質的に$A^\mu \leftrightarrow \Psi$の一意性が回復します。また第1定理について重要なのは,HK定理の場合と同様に,任意の4成分電流密度$j^\mu$が,ある4成分外部ポテンシャル$A^\mu$のもとで基底状態から必ず実現できるとは限らないという$V$-表示可能性問題が存在していることです。この問題については,4成分外部ポテンシャルのゲージ変換の自由度を除けば,4成分外部ポテンシャルから導かれた4成分波動関数と4成分電流密度には一対一対応関係が存在することが証明されています。ゲージ変換については,静電ポテンシャルの場合はHelmholtzの定理($|\mathbf{r}| \rightarrow \infty$で$V(\mathbf{r}) \rightarrow 0$)により固定されて変換の自由度がなくなることが知られています。以上の前提(ゲージ固定・非縮退など)のもとで,第1定理の主張は成り立ちます。
第2定理は,4成分電流密度で表されたHamiltonian演算子により求められたエネルギーは最低解をもつとする変分原理です。この定理については実はまだ厳格な証明が存在していません。というのは,Rajagopal-Callawayの定理に必要な量子電磁力学(Quantum electrodynamics; QED)が,繰り込み操作なしには変分原理が成り立たないからです。しかし,この繰り込み操作が正しいと仮定するのであれば,HK定理と同じように背理法によって,変分原理が縮退していない基底状態について成立することが証明できます。またこれに関連して,4成分電流密度で表わされたエネルギーに導関数が存在することが変分原理の前提条件ですが,これは保証されていません。
ここで,この定理を数学的に定式化するには,Legendre変換で普遍汎関数を導くのが便利です。まず
W\left[A^\mu\right]=\inf _{\Psi}\langle\Psi| \hat{H}\left[A^\mu\right]|\Psi\rangle, \quad j^\mu(\mathbf{r})=\frac{\delta W\left[A^\mu\right]}{\delta A_\mu(\mathbf{r})} \quad(22)
とおき,Legendre-Fenchel変換で
\mathcal{F}\left[j^\mu\right]=\sup _{A^\mu}\left\{\int d^3 r j^\mu(\mathbf{r}) A_\mu(\mathbf{r})-W\left[A^\mu\right]\right\} \quad(23)
を定義します。すると,外部ポテンシャル $A_{\mathrm{ext}}^\mu$ のもとでの基底エネルギーは
E\left[A_{\mathrm{ext}}^\mu\right]=\inf _{j^\mu \in \mathcal{J}_N}\left\{\mathcal{F}\left[j^\mu\right]-\int d^3 r j^\mu(\mathbf{r}) A_\mu^{\mathrm{ext}}(\mathbf{r})\right\} \quad(24)
となります($\mathcal{J}_N$は$N$表現可能な4成分電流密度の集合)。これがRajagopal-Callawayの定理に対応する変分原理の形です。前述のとおり,ゲージの固定と$V$-表示可能性に注意が必要です。後者の回避には,非相対論の場合と同様にLevyの制限付き探索の相対論拡張を用いて, $\mathcal{J}_N$上で普遍汎関数$\mathcal{F}$を定義・最小化するのが有効です(詳細は本稿では割愛します)。
Dirac-Kohn-Sham方程式
前節で述べた,Rajagopal-Callaway定理の第2定理の変分原理にもとづき,(1)式のDirac方程式を利用した相対論的DFTを考えると,次のDirac-Kohn-Sham方程式が得られます。
\left[c \boldsymbol{\alpha} \cdot \hat{\boldsymbol{p}}+\boldsymbol{\beta} m c^2+\alpha^\mu V^\mu\right] \psi_k(\boldsymbol{r})=\epsilon_k \psi_k(\boldsymbol{r}) \quad(25)
ここで,
V^\mu \equiv\left(-e \phi+v_{\mathrm{H}}+v_{\mathrm{xc}}, c\left(e \mathbf{A}+\mathbf{A}_{\mathrm{xc}}\right)\right) \quad(26)
ととれば(ここで$v_{\mathrm{H}}$はHartreeポテンシャル,$v_{\mathrm{xc}}$は交換相関ポテンシャルです。詳しくは後述します),
\alpha^\mu V_\mu=\left(-e \phi+v_{\mathrm{H}}+v_{\mathrm{xc}}\right)+c \boldsymbol{\alpha} \cdot (e\mathbf{A}+ \mathbf{A}_{\mathrm{xc}}) \quad(27)
となり(ここで$\alpha^\mu V_\mu$は,$V^0+\boldsymbol{\alpha} \cdot \mathbf{V}$の略記でありMinkowski計量での4-ベクトル縮約ではありません),最小結合
\hat{\mathbf{p}} \rightarrow \hat{\mathbf{p}}+e \mathbf{A}+\mathbf{A}_{\mathrm{xc}} \quad(28)
の形と一致します(静電問題では空間成分は$0$)。なお,本記事では$\mathbf{A}_{\mathrm{xc}}$は$e \mathbf{A}$と同じ次元に規格化しています。ただし,$\alpha^\mu$は4成分のDiracの$\alpha$行列
\alpha^\mu=(1, \boldsymbol{\alpha}) \quad(29)
です。4成分ポテンシャル$V^\mu$には4成分外部ポテンシャル$V_0^\mu$以外に多電子間相互作用ポテンシャルであるHartree4成分ポテンシャル$J^\mu$(定常状態系では時間成分のみ)と,4成分交換相関ポテンシャル$V_{\mathrm{xc}}^\mu$も含まれ,すべて4成分電流密度$j^\mu$で表わされます。すなわち,
V^\mu=V_0^\mu+J^\mu+V_{\mathrm{xc}}^\mu \quad(30)
です。ここで,
V_0^\mu(\mathbf{r})=(-e \phi(\mathbf{r}), ce \mathbf{A}(\mathbf{r})) \quad(31)
であり,
\begin{aligned}
& J^\mu(\mathbf{r})=\left(v_{\mathrm{H}}(\mathbf{r}), \mathbf{0}\right) \quad(32) \\
& v_{\mathrm{H}}(\mathbf{r})=\frac{\delta E_{\mathrm{H}}[n]}{\delta n(\mathbf{r})} \quad(33) \\
& E_{\mathrm{H}}[n]=\frac{e^2}{8 \pi \varepsilon_0} \iint \frac{n(\mathbf{r}) n\left(\mathbf{r}^{\prime}\right)}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|} d^3 r d^3 r^{\prime} \quad(34)
\end{aligned}
です。ここで,静的なCoulomb相互作用のみを仮定する通常のDFTでは,磁気自己場を解かないためベクトル成分は取りません。また,
\begin{aligned}
& V_{\mathrm{xc}}^\mu(\mathbf{r})=\left(v_{\mathrm{xc}}(\mathbf{r}), c\mathbf{A}_{\mathrm{xc}}(\mathbf{r})\right) \quad(35) \\
& v_{\mathrm{xc}}(\mathbf{r})=\frac{\delta E_{\mathrm{xc}}[n, \mathbf{j}]}{\delta n(\mathbf{r})} \quad(36) \\
& \mathbf{A}_{\mathrm{xc}}(\mathbf{r})=\frac{\delta E_{\mathrm{xc}}[n, \mathbf{j}]}{\delta \mathbf{j}(\mathbf{r})} \quad(37)
\end{aligned}
です。ここで,$v_{\mathrm{xc}}$には,非相対論的なKohn-Sham方程式のもの(LDAやGGAなど)がそのまま使われ,また通常計算では$\mathbf{A}_{\mathrm{xc}}=\mathbf{0}$とする電流非依存汎関数を採るのが一般的です。ここで電子密度は
n(\mathbf{r})=\sum_{k}^{occ} f_k \psi_k^{\dagger}(\mathbf{r}) \psi_k(\mathbf{r}) \quad(38)
で定義されます($f_k$は占有数です)。重要なのは,Dirac-Kohn-Sham方程式の解である(25)式の $\psi_k$が4成分スピノルであることです。ここでは,この$\psi_k$を軌道スピノルと呼び,$\epsilon_k$も軌道スピノルエネルギーと呼びます。Dirac-Kohn-Sham波動関数は通常,Kohn-Sham方程式と同様に,軌道スピノルで表わされたSlater行列式で表わされます。
\Psi\left(\boldsymbol{r}_1, \boldsymbol{r}_2, \cdots, \boldsymbol{r}_N\right)=\frac{1}{\sqrt{N!}} \operatorname{det}\left|\psi_1\left(\boldsymbol{r}_1\right) \psi_2\left(\boldsymbol{r}_2\right) \cdots \psi_N\left(\boldsymbol{r}_N\right)\right| \quad(39)
ここで$N$は電子数,$\boldsymbol{r}_i$は電子$i$の座標ベクトルです。
第一原理計算における相対論的補正
ここまでは,球対称ポテンシャルを仮定した原子系における相対論的補正と相対論的DFTについて述べてきました。この節では,周期系(固体)の第一原理計算,とくに平面波基底を前提に,相対論的補正の取り込み方を簡潔に説明します。実務では,相対論的効果を原子問題(擬ポテンシャル)側に埋め込み,計算時にはスカラー相対論(Scalar Relativistic; SR)を既定とし,必要に応じてノンコリニア(noncollinear)DFT によるSOC項を付加するのが主流です。ここで,擬ポテンシャル(PseudoPotential; PP)について詳しく説明しようとすると,それだけで記事が一つ出来上がってしまうので,ごく簡単に説明すると,内殻の効果を有効ポテンシャルに吸収し,価電子だけを明示的に扱う手法です。この前提のもとでSRを説明します。SRでは,PP生成の段階で原子のスカラー相対論方程式((14)式に相当するKoelling–Harmon型)を解き,質量速度補正項とDarwin項の効果をPPに取り込みます。これにより,固体で解くKohn–Sham方程式はSchrödinger型のまま,平面波基底で扱えます。SOCが必要な場合は,このSRを土台にノンコリニアDFTでSOCを付加します(ノンコリニアDFTに関しては,気が向いたら記事を執筆する予定です)。実際,定常状態・磁場なし($\mathbf{B} = \mathbf{0}$)かつ電流依存交換相関汎関数(これも詳しく解説しようとすると電流密度DFTの記事が別途必要になります)を用いない場合は$\mathbf{A}=\mathbf{A}_{\mathrm{xc}} = \mathbf{0}$とでき,Dirac-Kohn-Sham方程式の有効ポテンシャルは(26)式から分かるとおり,時間成分$-e \phi+v_{\mathrm{H}}+v_{\mathrm{xc}}$のみになります。したがって外部静電場があってもSRの枠組みで扱えます。この方針はVASPなどの平面波系コードで一般的であり,私が主に用いている数値原子軌道(Numerical Atomic Orbital; NAO)ベースのOpenMXにおいても基本的に同じ方針です。第一原理計算で相対論効果を一電子(Dirac–Coulomb)レベルで厳密に取り込むには,前節で述べたDirac–Kohn–Sham方程式をそのまま4成分(four-component; 4c)形式で解く必要があります。しかし,周期系での4c計算は,4成分交換相関ポテンシャル$V_{\mathrm{xc}}^\mu$の計算コストの高さに加え,負エネルギー状態の分離(no-pair近似,後述)や二電子相対論項の扱いなど実装上の困難が大きく,ほとんど用いられていません。
量子化学計算における相対論的補正
量子化学計算では,周期系の第一原理計算と比べて,より多様な相対論的手法が用いられます。4c Dirac–Kohn–Shamを直接解く方法に加えて,2成分(two-component; 2c)の有効HamiltonianであるZORA(Zero-Order Regular Approximation),DKH(Douglas–Kroll–Hess)法,X2C(exact two-component)法などが広く利用されます(周期系でも使われることはありますが,平面波系ではSR+SOCが主流で,2c計算の普及は相対的に限定的です)。また,二電子相互作用について,非相対論の瞬時Coulomb相互作用をQEDの$1/c^2$近似で補正するBreit相互作用を取り込む場合もあります。Breitの磁気項に相当するのがGaunt項であり,完全なBreitはこのGaunt項に加えて遅延(retardation)補正を含みます。これら(二電子相対論補正)は計算コストや実装上の理由から,分子系でも適用は比較的まれです。以下では,ZORA,DKH法,X2C法について述べたいと思います。
ZORA(Zero-Order Regular Approximation)
膨大な計算時間のかかる4c Dirac-Kohn-Shamに代わる相対論的補正として利用されるのが,2c法です。Dirac-Kohn-Sham計算には,計算時間の割にさほど重要ではないLarge成分とSmall成分とのpair項が存在します。この項を取り除く近似が前述したno-pair近似であり,このno-pair近似が2c法の基礎をなします。まず,改めてDirac-Kohn-Sham方程式の4成分スピノル波動関数を書くと,
\Psi_k=\left(\begin{array}{l}
\Psi_\alpha^{\mathrm{L}} \\
\Psi_\beta^{\mathrm{L}} \\
\Psi_\alpha^{\mathrm{S}} \\
\Psi_\beta^{\mathrm{S}}
\end{array}\right) \quad(40)
ですが,これを(25)式に代入し,(5a)式と(5b)式を使って整理すると,
\begin{aligned}
&c(\boldsymbol{\sigma} \cdot \hat{\boldsymbol{p}}) \psi_k^{\mathrm{S}}+V \psi_k^{\mathrm{L}}=\epsilon_k \psi_k^{\mathrm{L}} \quad(41a) \\
&c(\boldsymbol{\sigma} \cdot \hat{\boldsymbol{p}}) \psi_k^{\mathrm{L}}+\left(-2 m c^2+V\right) \psi_k^{\mathrm{S}}=\epsilon_k \psi_k^{\mathrm{S}} \quad(41b)
\end{aligned}
と書けます。ここで,$\psi_k^{\mathrm{L}}$,$\psi_k^{\mathrm{S}}$は2成分の軌道スピノルです(それぞれ(5)式の$\phi$と$\chi$に相当)。こうすると,Small成分軌道スピノルはLarge成分軌道スピノルで表わせ,
\psi_k^{\mathrm{S}}=\left(\epsilon_k+2 m c^2-V\right)^{-1} c(\boldsymbol{\sigma} \cdot \hat{\boldsymbol{p}}) \psi_k^{\mathrm{L}}=\hat{K} \frac{\boldsymbol{\sigma} \cdot \hat{\boldsymbol{p}}}{2 m c} \psi_k^{\mathrm{L}} \quad(42)
と書けます。$\hat{K}$ は,
\hat{K}=\left(1+\frac{\epsilon_k-V}{2 m c^2}\right)^{-1} \quad(43)
です。電子運動が光速より十分に遅い場合,この $\hat{K}$ は,
\hat{K}=\left(1+\frac{\epsilon_k-V}{2 m c^2}\right)^{-1} \approx 1-\frac{\epsilon_k-V}{2 m c^2}+\cdots \quad(44)
と展開できます。(44)式の$\hat{K}$の展開において問題になるのは,ポテンシャルが発散($V \rightarrow \infty$)する原子核近傍です。この問題を回避するため,
\begin{aligned}
\left(\epsilon_k+2 m c^2-V\right)^{-1} & =\left(2 m c^2-V\right)^{-1}\left(1+\frac{\epsilon_k}{2 m c^2-V}\right)^{-1} \\
& =\left(2 m c^2-V\right)^{-1} \hat{K}^{\prime}
\end{aligned} \quad(45)
と定義される$\hat{K}^{\prime}$で置き換える方法が使われます。この置き換えにより,$\epsilon_k /\left(2 m c^2-V\right) \ll 1$なので,展開がいつも良い近似になります。さらには$\hat{K}^{\prime}=1$とおくこともかなり良い近似です。この近似をZORAと呼びます。ZORAに基づいた相対論的Kohn-Sham方程式は,
[(\sigma \cdot \hat{\boldsymbol{p}}) \frac{c^2}{2 m c^2-V}(\sigma \cdot \hat{\boldsymbol{p}})+V ]\psi_i=\epsilon_i^{\mathrm{ZORA}} \psi_i \quad(46)
となります。ここで,この(46)式の左辺かぎ括弧の中の第1項を$1/c^2$で展開すると,
(\boldsymbol{\sigma} \cdot \hat{\boldsymbol{p}}) \frac{c^2}{2 m c^2-V}(\boldsymbol{\sigma} \cdot \hat{\boldsymbol{p}}) \simeq \frac{\hat{\mathbf{p}}^2}{2 m}+\frac{1}{8 m^2 c^2}\left\{V, \hat{\mathbf{p}}^2\right\}+\frac{\hbar^2}{8 m^2 c^2} \nabla^2 V+\frac{1}{4 m^2 c^2} \boldsymbol{\sigma} \cdot(\nabla V \times \hat{\mathbf{p}})+\mathcal{O}\left(c^{-4}\right) \quad(47)
となり,(47)式を(18)式と比較すると,Darwin項とSOC項は再現できますが,質量速度補正項は与えられません(代わりに$\frac{1}{8 m^2 c^2}\lbrace V, \hat{p}^2\rbrace$という項を得ます)。ここで,ZORA-Kohn-Sham法で非常によく利用されるのは,スカラーZORAを使う2c法です。この方法では,Pauliスピン行列を単位行列に置き換えて$\sigma \cdot \hat{\boldsymbol{p}}$を$\hat{\boldsymbol{p}}$とすることにより,Kohn-Sham計算で律速にならないほど劇的に計算時間を減らすことができます。 この置き換えは,(47)式の右辺第4項である,SOC項を無視することと同等です。それにもかかわらず,エネルギーに与える相対論的効果の大部分を占める質量と速度による効果を取り込むことができるので,きわめて効率的な近似であるといえます。しかし,高スピン系やきわめて重い原子を含む系の電子物性はSOCを無視して定量的に評価することはできないので,スカラーZORAは問題となります。
DKH(Douglas–Kroll–Hess)法
前節で2c法の一つであるZORAを説明しましたが,本節ではDKH法について説明します。DKH法はFW変換の拡張として捉えることができ,FW変換での展開パラメ-タである光速$c$のかわりに外部ポテンシャル$V$を展開パラメ-タとします。外部ポテンシャルを展開パラメ-タとして用いることでFW変換の持つ強い特異性を避けることができますし,変分的に安定です。まず,相対論的1電子Dirac Hamiltonian(電流や外部磁場はとりあえず無視します)は,Large/Small成分で$2 \times 2$ブロックに分けると
\hat{H}_D=\left(\begin{array}{ll}
\hat{H}_{L L} & \hat{H}_{L S} \\
\hat{H}_{S L} & \hat{H}_{S S}
\end{array}\right)=\beta m c^2+\hat{E}+\hat{O} \quad(48)
と書けます。ここで偶演算子(even)$\hat{E}$は,$[\beta, \hat{E}]=0$(ブロック対角),奇演算子(odd)$\hat{O}$は${\beta, \hat{O}}=0$(ブロック非対角)で,自由粒子なら$\hat{O}=c \boldsymbol{\alpha} \cdot \hat{\boldsymbol{p}} ,\hat{E}=V(\mathbf{r})$です。ここでの目的は,ユニタリ変換で奇演算子を消してブロック対角化し,Largeだけの有効2成分Hamiltonianを得ることです。すなわち,
U_{\mathrm{FW}}^{\dagger} \hat{H}_D U_{\mathrm{FW}}=\left(\begin{array}{cc}
\hat{H}^L & 0 \\
0 & \hat{H}^S
\end{array}\right) \quad(49)
とすることです。ここで,ブロック対角化後の$\hat{H}^L$が電子用の2c有効Hamiltonianです。前述したFW変換は,Coulomb場のような強い外場では高次($1/c^4$以降)に核近傍で特異な($\delta$型等の)項や高階微分が多数現れ,打切り近似の数値安定性・収束性の面で実用上扱いにくくなります。この問題を段階的ユニタリ変換の積で回避し,奇成分を摂動次数ごとに消していくよう定式化したのがDouglas-Kroll(のちHessにより改良)によるDKH法です。変換が有限回ならDKH2,DKH3,…と呼称します。DKHの逐次ブロック対角化(2c化)の具体的手順を以下で説明します。Dirac Hamiltonianを,正負エネルギー(Large/Small)成分の混合を表す「偶(even)成分」と「奇(odd)成分」に分解するところから始めます。記号$\beta$はDirac行列, ${\cdot, \cdot}$は反交換子です。 $\beta$と可換なものを偶,反可換なものを奇と呼ぶ約束にすると,
\hat{H}_D=\underbrace{\beta m c^2+V(\mathbf{r})}_{\hat{E}}+\underbrace{c \boldsymbol{\alpha} \cdot \hat{\boldsymbol{p}}}_{\hat{O}}, \quad[\hat{E}, \beta]=0, \quad\{\hat{O}, \beta\}=0 \quad(50)
ここで$\hat{E}$はLarge/Smallブロック対角(偶),$\hat{O}$はそれらを結ぶブロック非対角(奇)です。目的は,ユニタリ変換$\hat{U}=\exp (\hat{S})$を段階的に構成して$\hat{H}^{\prime}=\hat{U} \hat{H}_D \hat{U}^{\dagger}$をブロック対角化(奇成分の逐次消去)することです。生成子 $\hat{S}$は反エルミ-トかつ奇($\hat{S}^{\dagger}=-\hat{S},{\hat{S}, \beta}=0$)にとるのが肝要で,この選び方により変換がブロック構造を保ったまま奇成分を弱めていきます。変換はBaker-Campbell-Hausdorff(BCH)展開で,
\hat{H}^{\prime}=\hat{H}_D+\left[\hat{S}, \hat{H}_D\right]+\frac{1}{2}\left[\hat{S},\left[\hat{S}, \hat{H}_D\right]\right]+\cdots \quad(51)
と書けます。DKHは摂動パラメータ$1 / c$で秩序付けを行い,各段で$\hat{S}$を選び直してその段の精度に現れる奇成分を消すようにします。
ここで,まず最小の近似である,一次デカップリング(DKH1)について説明します。最初の$\hat{S}_1$ は,自由粒子の大きなギャップ$2 m c^2$を使って最も強い混合$\hat{O}$を打ち消すよう
\hat{S}_1=\frac{\beta \hat{O}}{2 m c^2} \quad(52)
と選びます($\beta \hat{O}$ は反エルミ-トなので$\hat{S}_1$は反エルミート)。この選択で一次の奇成分が最小化され,変換後のHamiltonianは,
\hat{H}^{(1)}=\hat{E}+\beta \frac{\hat{O}^2}{2 m c^2}+\underbrace{\frac{1}{2 m c^2} \beta[\hat{O}, \hat{E}]}_{\text {奇 }}+\cdots \quad(53)
となります。ここで $\hat{O}^2=c^2 \hat{\boldsymbol{p}}^2$ を用いると,Largeブロック($\beta \rightarrow+1$)には,
m c^2+\frac{\hat{\boldsymbol{p}}^2}{2 m}+V(\mathbf{r}) \quad(54)
が現れ,DKH1(Large ブロック)は素直に非相対論Hamiltonian(静止質量を差し引けば $\hat{\boldsymbol{p}}^2 / 2 m+V$ )に一致します(実用上はこれは無意味であり,次段のDKH2以降が使われます)。一方,まだ $[\hat{O}, \hat{E}]$ 由来の残りの奇成分が残っています。これが次段での消去対象です。
第2段(DKH2)では, $\hat{H}^{(1)}$ に残った奇成分 $\hat{O}^{(1)}=\frac{1}{2 m c^2} \beta[\hat{O}, \hat{E}]$ を消すように,
\hat{S}_2=\frac{\beta \hat{O}^{(1)}}{2 m c^2}=\frac{1}{4 m^2 c^4}[\hat{O}, \hat{E}] \quad(55)
を選び, $\hat{H}^{(2)}=\exp \left(\hat{S}_2\right) \hat{H}^{(1)} \exp \left(-\hat{S}_2\right)$ を作ります。BCHを二重交換子まで評価すると,Largeブロックの偶成分に
\hat{H}_L^{(\mathrm{DKH} 2)}=m c^2+V+\frac{\hat{\boldsymbol{p}}^2}{2 m}-\frac{\hat{\boldsymbol{p}}^4}{8 m^3 c^2}-\frac{1}{8 m^2 c^2}[\hat{O},[\hat{O}, V]] \quad(56)
が現れます。最後の二重交換子を$\hat{O}=c \boldsymbol{\alpha} \cdot \hat{\boldsymbol{p}}$で具体的に計算すると
-\frac{1}{8 m^2 c^2}[\hat{O},[\hat{O}, V]]=\frac{\hbar^2}{8 m^2 c^2} \nabla^2 V+\frac{1}{4 m^2 c^2} \boldsymbol{\sigma} \cdot(\nabla V \times \hat{\boldsymbol{p}})
\quad(57)
となります。すなわちDarwin項とSOC項が得られ,前述したFW変換での$1 / c^2$補正の(18)式と一致します(球対称ポテンシャルでは(57)式の右辺第2項は$\frac{1}{2 m^2 c^2} \frac{1}{r} \frac{d V}{d r} \mathbf{L} \cdot \mathbf{S}$となります)。これがDKH2と呼ばれる有効2成分Hamiltonianです。第3段(DKH3)以降も同じ考え方を繰り返します。一般に, $(k-1)$ 段まで終えたHamiltonianを
\hat{H}^{(k-1)}=\hat{E}^{(k-1)}+\hat{O}^{(k-1)}
\quad(58)
(ただし$\hat{E}^{(k-1)}$偶,$\hat{O}^{(k-1)}$奇)と書けば,その段の生成子は
\hat{S}_k=\frac{\beta \hat{O}^{(k-1)}}{2 m c^2}
\quad(59)
と選べばよく,$\hat{H}^{(k)}=\exp \left(\hat{S}_k\right) \hat{H}^{(k-1)} \exp \left(-\hat{S}_k\right)$によりその段の秩序に現れる奇成分が消え,代わりにより高次(より小さい)偶成分補正が積み上がることになります。こうして$k$を上げるほど,Largeブロックの2成分有効Hamiltonianに $1 / c^4, 1 / c^6, \ldots$の高次運動エネルギー補正や高次スピン結合(例えばThomas項の高次修正,$\left\lbrace \hat{\boldsymbol{p}}^2, \nabla^2 V \right\rbrace$と$\hat{\boldsymbol{p}}$のより高階の交換子に由来する項)が系統的に付加されていきます。有限次数で打ち切れば,その次数までブロック対角化(奇成分の消去)が厳密に満たされ,残りの奇成分はより高次の小さな量に追いやられています。
通常の化学の問題においては,DKH2法は満足のいく結果を与えます。MolzbergerとSchwarz[2]による数値的な解析によると,DKH2法は$Z^6 / c^4$のオーダーまでエネルギーをリカバ-し,高次の相対論効果の大部分を考慮することができます。しかし,非常に重い元素や核近接性質にはDKH3以上が有利です。
X2C(exact two-component)法
X2C法は,Dirac-Kohn-Sham方程式を有限基底上で一度の手続きで厳密にブロック対角化して,電子(正エネルギー)部分の2成分有効Hamiltonianを得る方法です。DKHが$1 / c$の摂動級数として奇成分を順次消していくのに対し,X2Cは有限基底での正エネルギー固有ベクトルから直接デカップリング行列を構成し,事実上DKHの無限次(DKH $\infty$)に相当する2成分化を一発で行います。しかし,X2Cはあくまで行列表現(採用基底)内で厳密であることに注意してください。核近傍での発散や高階微分による数値不安定を避けられ,変分的にも安定です。実用上はスピンフリーX2C(SFX2C)でスカラー相対論を構築し,必要なら一電子/二電子のSOCを(厳密または平均場近似で)付加する構成が広く使われます。
一電子Hamiltonianを,Large/Small成分で$2 \times 2$ブロックに書くと(外部磁場なし)
\hat{H}_D=\left(\begin{array}{ll}
\hat{H}_{L L} & \hat{H}_{L S} \\
\hat{H}_{S L} & \hat{H}_{S S}
\end{array}\right) \quad(60)
となります。ここで $\hat{H}_{L L}=V(\mathbf{r})+m c^2, \hat{H}_{S S}=V(\mathbf{r})-m c^2$,$\hat{H}_{L S}=\hat{H}_{S L}=c \boldsymbol{\alpha} \cdot \hat{\boldsymbol{p}}$の像(後述の基底化で厳密には行列)です。目的は,ユニタリ変換$\hat{U}$により,
\hat{U}^{\dagger} \hat{H}_D \hat{U}=\left(\begin{array}{cc}
\hat{H}^{(+)} & 0 \\
0 & \hat{H}^{(-)}
\end{array}\right) \quad(61)
とブロック対角化し,正エネルギーブロック$\hat{H}^{(+)}$を2成分有効Hamiltonianとして用いることです。X2Cは,Small成分をLarge成分で表すデカップリング演算子$\hat{X}$(行列)を
\psi^{\mathrm{S}}=\hat{X} \psi^{\mathrm{L}} \quad(62)
と仮定して構成します。これを満たす$\hat{X}$を正エネルギー固有ベクトルから決定し,非直交化(正規化)行列$\hat{R}$を
\hat{R}=\left(\mathbf{1}+\hat{X}^{\dagger} \hat{X}\right)^{-1 / 2} \quad(63)
と定めると,電子用2成分有効Hamiltonianは厳密に
\hat{H}^{\mathrm{X} 2 \mathrm{C}}=\hat{R}^{\dagger}\left(\hat{H}_{L L}+\hat{H}_{L S} \hat{X}+\hat{X}^{\dagger} \hat{H}_{S L}+\hat{X}^{\dagger} \hat{H}_{S S} \hat{X}\right) \hat{R} \quad(64)
で与えられます。通常は静止質量を差し引いた
\hat{H}_{\mathrm{el}}^{\mathrm{X} 2 \mathrm{C}}=\hat{H}^{\mathrm{X} 2 \mathrm{C}}-m c^2 \quad(65)
を2成分電子Hamiltonianとして用います。(64)式は基底に依らない一般形で,有限基底に投影した行列として実装されます。なおSFX2Cとは,X2C法のスカラー相対論版であり,4成分Dirac(またはDirac-Kohn-Sham)方程式を有限基底上でブロック対角化して電子(正エネルギー)だけの2成分有効Hamiltonianを作るとき,スピン依存($\boldsymbol{\sigma}$ を含む)項,すなわちSOC項を最初は入れないでデカップリングする流儀を指します。SFX2Cは,質量速度補正項とDarwin項を含むスカラー相対論Hamiltonianが得られます。そこに,一電子SOC項,すなわち(57)式の右辺第2項と全く同じ項
\hat{H}_{\mathrm{SOC}}^{(1)}=\frac{1}{4 m^2 c^2} \boldsymbol{\sigma} \cdot(\nabla V \times \hat{\boldsymbol{p}}) \quad(66)
に加え,二電子SOC(Breit相互作用のGaunt相互作用由来)を平均場近似(AMFI/MMFなど)で補う項を加えて2c(ノンコリニア)計算を行うのが通例です。相対論的DFTでは二電子SOCを完全に2cで扱うのは重く,原子平均場(AMFI)や分子平均場(MMF)近似が広く使われます。
ZORA,DKH法,X2C法の比較
以下では,ZORA,DKH法,X2C法を,主に計算精度と計算コストの面から簡単に比較したいと思います。
1.物理内容と近似の性格
- ZORA
Small成分の消去で現れる$\left(\epsilon+2 m c^2-V\right)^{-1}$を$\left(2 m c^2-V\right)^{-1}$に置き換えて正則化し,さらに必要なら $1 / c^2$ でテイラ-展開して2c化します。
これにより核近傍の特異性を回避しつつ,Darwin項とSOC項は正しく再現しますが,FW変換/DKHに出る質量速度補正項$-\hat{\mathbf{p}}^4 /\left(8 m^3 c^2\right)$は$\frac{1}{4 m c^2}\left\lbrace V, \hat{\mathbf{p}}^2\right\rbrace$のような形に置き替わります(近似の仕方の違いによるもの)。スカラーZORAではSOC項を落とした2cスカラー化が得られます。
長所:定式化が簡潔,行列構築が軽い,核近傍の数値安定性は良好。
注意:厳密なブロック対角化ではないため,演算子の表現変換をZORA形で整合的に定義する必要があり,磁場や電流を扱う場合はゲージ依存に配慮が要ります。ここで,演算子の表現変換とは,演算子の表現変換(picture change)とは,4c Dirac表現から2c(有効電子空間)表現へ基底を変えたときに,物理量の演算子も同じユニタリ変換で写像し直すことを指します。 - DKH法
生成子$\hat{S}$を用いたユニタリ変換の逐次ブロック対角化で奇ブロックを消していく方法。DKH2はFW変換の$1 / c^2$展開と等価で,質量速度補正項・Darwin項・SOC項を正しい形で与えます。DKH3以降は$1 / c^4, 1 / c^6, \ldots$の高次補正を系統的に取り込め,核近傍の発散的な挙動を避けるよう正則化された演算子が得られます。
長所:次数を上げれば精度を単調に改善できます(有限次数で厳密にその次数までブロック対角)。一電子演算子の表現変換も一貫しています。
注意:次数を上げるほど交換子列の評価や積分変換の手間が増えます。実用上はDKH2とDKH3がコストパフォ-マンスに優れ,「定番」とされます。 - X2C法
有限基底上で4c—電子行列を対角化し,正エネルギー固有ベクトルからデカップリング行列$X$と正規化$R$を構成することで,一段のユニタリ変換で厳密(=その基底内で無近似)に2c化します。結果は「無限次DKH を一気にやった」のと等価で,スピンフリーのSFX2CやSOCを含むSO-X2Cに分けて使えます。SFX2C-1e(スカラー)で土台を作り,SOCは1e-SOCと,必要に応じて2e-SOCもAMFI/MMFで付加します。
長所:高精度かつ安定,一電子演算子の表現変換も自然に付随。通常はDKH3以上に匹敵する精度をDKH2並みの手間で得られます。
注意:二電子SOCは平均場(AMFIなど)で補うのが一般的で,基底はRKB等の相対論的にバランスが取れた基底が望ましいとされます。
2.数値安定性と実装上の挙動
- 核近傍の特異性:ZORA,DKH法,X2C法はいずれも4c直接解より扱いやすいとされます。FW変換の高次打切りで問題化する$\delta$型・高階微分の扱いを,DKH法は正則化,X2C法は基底内厳密化で回避します。
- ゲージ問題:外部磁場や電流密度を扱う場合,ZORAはゲージ依存性に敏感。DKH法・X2C法はRKBやGIAO併用で良好です。
- 基底の相性:DKH法・X2C法用に再収縮・最適化された基底(DK-,x2c-系など)があるとベストです。ZORAもバランスの取れた基底が無難とされます。
3.計算コスト(相対比較)
- ZORA:最軽量で,1回の有効運動エネルギー演算子の組み立ておよび通常の2c SCFのみで可能です。平面波やNAOでも負担が小さいとされます。
- DKH2:ZORAよりやや重いですが軽量級です。交換子の評価と一電子行列の変換が主コストです。
- DKH3/4:次数とともに上昇します(交換子鎖が増える)。それでも4c直解よりは大幅に軽いです。
- X2C:初回に4cの固有値問題を解きますが,その後の計算を2cに集中できるので実用上は軽いです。計算コストはDKH2と同程度からやや上で,DKH3より軽いことが多いです。
-参考順位(概念図):
$$
\mathrm{ZORA}<\mathrm{DKH} 2 \approx \mathrm{X} 2 \mathrm{C}<\mathrm{DKH} 3 / 4 \ll 4 \mathrm{c} \text { (Dirac 直解) }
$$
Breit/Gaunt相互作用の取り扱い
非相対論の電子-電子相互作用はCoulomb$1 / r_{12}$のみですが,これはLorentz共変ではないので,力が瞬時に働くと仮定しており,厳密には量子電磁力学(QED)で記述される遅延(retardation)や磁気的相互作用の効果を落としています。量子化学計算では,QEDの厳密式を使わずに微細構造定数($1 / c$)展開で導かれるBreit相互作用を,Coulombに対する相対論補正として加えます。核-電子のQED放射補正(自己エネルギーや真空分極)は高次で小さく無視できる一方,電子-電子相互作用の相対論補正(Breit相互作用)は$1 / c^2$次で主寄与となる,というのが実用上の出発点です。
二電子Hamiltonianの有効相互作用を,CoulombにBreit相互作用を足した形で書くと,
\hat{H}_{12}^{(\mathrm{DCB})}=\frac{e^2}{4 \pi \varepsilon_0 r_{12}}+\hat{B}_{12}, \quad \hat{B}_{12}=-\frac{e^2}{8 \pi \varepsilon_0 r_{12}}\left[\boldsymbol{\alpha}_1 \cdot \boldsymbol{\alpha}_2+\frac{\left(\boldsymbol{\alpha}_1 \cdot \mathbf{r}_{12}\right)\left(\boldsymbol{\alpha}_2 \cdot \mathbf{r}_{12}\right)}{r_{12}^2}\right] \quad(67)
となります(Coulombゲージ)。ここで$\boldsymbol{\alpha}$はDiracの$\alpha$行列,$\mathbf{r}_{12}=\mathbf{r}_1-\mathbf{r}_2$です。かぎ括弧の第1項が磁気的なGaunt相互作用項$\left(-\frac{1}{2 r_{12}} \boldsymbol{\alpha}_1 \cdot \boldsymbol{\alpha}_2\right)$ ,第2項が遅延(retardation)項に相当します(全体が$1 / c^2$次)。それぞれの項の物理的意味は,Gaunt相互作用項が電流(スピン・運動)同士の相互作用で,スピン-軌道・スピン-スピン等の二電子スピン結合の一部と同じ起源であり,EPR物性・超微細相互作用・重元素のスピン依存物性に寄与します。また,遅延項は光子交換に伴う有限伝播速度の効果を表す補正を表し,一般にGaunt項よりさらに小さくなります。いずれの項も$1 / c^2$ 次(非相対論の上に乗る相対論補正)で,Coulombにくらべ総じて小さく,通常の結合エネルギーや構造に対する影響はかなり限定的,という経験則があります。
多電子系の一電子+二電子をまとめると,よく使われる形は,
\hat{H}_{\mathrm{DCB}}=\sum_i\left[c \boldsymbol{\alpha}_i \cdot \hat{\mathbf{p}}_i+\beta_i m c^2+V_{\mathrm{ne}}\left(\mathbf{r}_i\right)\right]+\sum_{i<j}\left[\frac{e^2}{4 \pi \varepsilon_0 r_{i j}}-\frac{e^2}{8 \pi \varepsilon_0 r_{i j}}\left(\boldsymbol{\alpha}_i \cdot \boldsymbol{\alpha}_j+\frac{\left(\boldsymbol{\alpha}_i \cdot \mathbf{r}_{i j}\right)\left(\boldsymbol{\alpha}_j \cdot \mathbf{r}_{i j}\right)}{r_{i j}^2}\right)\right] \quad(68)
であり,(68)式をDirac–Coulomb–Breit Hamiltonianと呼びます。(68)式から,スピン依存(SO)とスピンフリー(SF)成分の厳密分離を行う定式化(Dyallの分離)も確立しており,後述のX2C/AMFなどの実装で活用されます。実用計算では,Breit相互作用のうちGaunt相互作用だけを取り込むDirac–Coulomb–Gaunt(DCG)近似がよく使われます(Breit相互作用を完全に取り込んだ計算はDCBと呼ばれます) 。常田は,Breit相互作用の標準形とGaunt近似の位置づけをまとめ,化学では二電子ポテンシャルの変形(=Breit相互作用起因)の寄与は概して小さいことを指摘しています。
Breit/Gauntを取り込むと,スピンに依存する二電子スピン-軌道(2e–SO)効果が平均場を通じて一電子 SOCを「遮蔽」する(核のみから生じるSOCを弱める)方向に働きます。これを経験的に模倣するSNSO(Screened-Nuclear Spin–Orbit)のような簡便法もありますが,厳密性・拡張性の観点からは二電子演算子そのものを二成分枠組みへ正しく写像して扱うのが理想です。歴史的にはAMFI(Atomic Mean-Field Integral)に代表される平均場近似が開発され,スカラー相対論ハミルトニアン(たとえばDKH2)に1電子SOCと2電子SOCの平均場項を組み合わせる手順が実装されてきました。AMFIは4cに比べて大幅に軽く,スピン-軌道分裂や励起エネルギー,分光定数などで良好な性能を示すことが報告されています。
2c理論,とくにX2Cにおける二電子相対論の扱いでは,「演算子の表現変換」をどう近似するかが要点になります。4成分空間で定義されたHamiltonianやCoulomb積分は,2成分の有効電子空間へユニタリに写し替えられねばなりません。前述したように,X2Cはこの写像(デカップリング)を有限基底上で厳密に構成できるが,二電子部分の完全な表現変換(2e-pc)は分子計算では計算負荷が大きいとされます。
そこで,「原子ごとに4成分の平均場を解き,その情報だけを分子計算に持ち込む」というAMF (Atomic Mean-Field)の発想が生まれ,2e-pcを原子ブロックの差分Fock行列として再現するamfX2Cが提案されました。amfX2Cは4cで得た原子Fock行列を正確に2cへブロック対角化し,分子では未変換の2e積分との差分を原子ブロックに埋め戻すことで2e-pcを近似的に再現します。評価は原子ごとで済むため,系サイズに対して線形(あるいは原子タイプの重複が多いと亜線形)にスケールします。
Breit/Gaunt相互作用を含むより汎用的な平均場としては,X2CAMFが整備されています。X2CAMFは4c DFTを原子で実行し,X2Cデカップリングで2cへ変換した分子平均場(MMF)を構築する厳密版(X2CMMF)を出発点に,分子ではコストの重い2e相対論項を原子由来のモデルで置き換えます。X2CAMFはCoulomb相互作用のスピン依存部分(2e-SO)を厳密に分離し,それをAMF積分として取り入れるとともに,Gaunt相互作用(場合によっては完全なBreit相互作用)もAMFに含められるため,2e-SOとGaunt相互作用を同一枠組みで効率的に扱えます。化学的に重要な性質に対しては,重原子を含む分子でもGaunt相互作用の寄与は小さく,適切な基底(拡散関数を過度に含めないVTZ級など)を選べば,X2CAMF (DCG) は結合長や調和振動数へのGaunt相互作用の補正を化学精度で再現できると報告されています。一方,スカラー2e-pcを1中心近似で足し込むamfX2Cは,拡散関数の扱いによって数値安定性の問題が生じ得ることが指摘されていて(拡散関数ブロックからのAMF積分の除外で安定化する実務策もあります),2e-SOのように局在性の高い項に比べてスカラー2e-pcの1中心近似は慎重な検証が必要です。
まとめると,Breit相互作用は二電子相互作用の相対論補正であり,Gaunt相互作用がその磁気(スピン依存)主要部を担います。化学量への定量的影響は一般に小さいとはいえ,超重元素や高精度の分光定数,超微細結合定数,磁気物性などでは無視できない局面があります。4cでDCBを直接解くのが最も直接的ですが,実用上は2cのX2Cを土台に,二電子演算子の表現変換をAMFI/MMFあるいはモデルポテンシャルで効率よく近似し,必要に応じてGaunt相互作用(場合により完全なBreit相互作用)を加える設計が,精度と計算量の両面でバランスが取れています。常田が強調するように,化学では二電子ポテンシャルの相対論的変形そのものは小さいですが,その取り込み方(Gaunt相互作用の近似にとどめるのか,どの平均場スキームで2e-pcを補うのか,拡散関数の扱いをどう安定化するか)が,最終的な信頼性を左右します。
参考文献
[1] A. H. MacDonald, W. E. Pickett, and D. Koelling, "A linearised relativistic augmented-plane-wave method utilising approximate pure spin basis functions," J. Phys. C: Solid State Phys. 13:2675-2683, 1980.
[2] K. Molzberger, W. H. E. Schwarz, Theor. Chim. Acta, 94, 213 (1996).
[3] R.M.マーチン 『物質の電子状態 上』シュプリンガー・ジャパン株式会社(2010)
[4] T. Tsuneda, J. Comput. Chem. Jpn. 13, 71–82 (2014).
[5] T. Nakajima, J. Comput. Chem. Jpn. 13, 50–70 (2014).
[6] A. Wodyński, and M. Kaupp, J. Chem. Theory Comput. 16, 314–325 (2020).
[7] X. Wang, C. Zhang, J. Liu, and L. Cheng, arXiv:2504.19479 (2025).
[8] S. Knecht, M. Repisky, H. J. A. Jensen, and T. Saue, J. Chem. Phys. 157, 114106 (2022).
[9] D. D. Koelling and B. N. Harmon, J. Phys. C 10, 3107 (1977).
[10] E. van Lenthe, E. J. Baerends, and J. G. Snijders, J. Chem. Phys. 99, 4597 (1993).
[11] M. Douglas and N. M. Kroll, Ann. Phys. (N.Y.) 82, 89 (1974).
[12] B. A. Hess, Phys. Rev. A 33, 3742 (1986).
[13] G. Breit, Phys. Rev. 34, 553 (1929).
[14] J. A. Gaunt, Proc. R. Soc. Lond. A 122, 513 (1929).
[15] J. Sucher, Phys. Rev. A 22, 348 (1980).
[16] L. L. Foldy and S. A. Wouthuysen, Phys. Rev. 78, 29 (1950).