量子力学の数値計算の難しさ
量子力学で解析的に (つまり紙と鉛筆で) その性質を調べられる物理系といえば、水素原子とか調和振動子、あとは簡単なスピンモデルくらいで、他のものに関しては全て数値計算によって調べる必要がある。とは言っても、一般にこの数値計算は、考える系に含まれる粒子数が多くなるととても難しくなってしまう。
この難しさは主に、
- 波動関数を明示的に保持するために必要なメモリ量が指数的に増えること。
- 物理量を求めるための積分や離散和に指数的な時間がかかってしまうこと。
の 2 点から生まれるものだ。
例を挙げてみよう。
例えば「\(N\) 個のスピン 1/2 粒子の一般の状態を表現する」というタスクを考えよう。一般のスピン1/2粒子の状態はパウリ行列 \(\sigma_z\) の固有値 \(+1, -1\) に対応する固有状態をそれぞれ \(\ket{0}, \ket{1}\) として
\begin{align}
\ket{\psi} = \sum_{b_i \in \{0,1\}} c_{b_1\cdots b_N} \ket{b_1} \otimes \ket{b_2} \otimes \cdots \otimes \ket{b_N}
\end{align}
と書ける。係数 \(c_{b_1\cdots b_N}\) の添字 \(b_1\cdots b_N\) には \(2^N\) 通りの組み合わせがあるので、このスピン系の状態を一般に表すには、\(2^N\) 個の係数 \(c_{b_1\cdots b_N}\) を保持しなくてはならなくなる。これは \(30\) スピンくらいで、通常のコンピュータのメモリ量を超えてしまう計算だ。もちろん、この状態に対して適当な物理量 \(O\) の期待値 \(\bra{\psi}O\ket{\psi}\) を評価しようとすれば、さらに計算量は増えるだろう。
また、他にも「\(N\) 電子系の適当な状態のエネルギー期待値を求める」という数値計算を考えてよう。 \(\b{r}_i\) とスピン \(\sigma_i\) を一緒くたにしたものを \(\b{x}_i\) と書けば、電子の波動関数は \(\psi(\b{x}_1, \b{x}_2, \cdots, \b{x}_N)\) という形の関数となる。通常電子のハミルトニアンは
\begin{align}
H = \sum_i \left(-\frac{\hbar^2}{2m}\nabla_i^2 + U(\b{x}_i)\right) + \sum_{ij} U_{\text{int}}(\b{x}_i, \b{x}_j)
\end{align}
の形である。このハミルトニアンを用いて、波動関数 \(\psi(\b{x}_1, \b{x}_2, \cdots, \b{x}_N)\) のエネルギー期待値 \(E_{\psi}\) は
\begin{align}
E_{\psi} = \int d\b{x}_1 \cdots d\b{x}_N \psi^*(\b{x}_1, \b{x}_2, \cdots, \b{x}_N) H \psi(\b{x}_1, \b{x}_2, \cdots, \b{x}_N)
\end{align}
となる。ここでは簡単のため、ある与えられた座標 \((\b{x}_1, \b{x}_2, \cdots, \b{x}_N)\) に対して、\(\psi(\b{x}_1, \b{x}_2, \cdots, \b{x}_N)\) やそれにハミルトニアンを作用させた関数 \(H\psi(\b{x}_1, \b{x}_2, \cdots, \b{x}_N)\) が効率的に計算できるような波動関数 \(\psi(\b{x}_1, \b{x}_2, \cdots, \b{x}_N)\) が与えられているとしよう。その場合であっても、上記の積分を数値的に評価するために、例えば各軸を \(M\) 等分して計算しようとすれば、グリッドの数は \(M^N\) となり、これも指数回の計算が必要になってしまう。
このような計算量の指数関数的爆発を、乱数の力を使って改善しようと試みるのが
量子モンテカルロ法 (quantum Monte-Carlo, QMC)
である。
量子モンテカルロ法
さて、上にも述べた通り、
量子モンテカルロ法とは、量子力学に従う物理系の性質を、乱数の使用による確率的な計算によって求める手法の総称である。「総称」と書いたのは、量子モンテカルロには多種多様なバリエーションがあって、それぞれに〇〇量子モンテカルロ法といった名前がついているからだ。例えば経路積分量子モンテカルロ法、変分量子モンテカルロ法などがある。単に量子モンテカルロ法といった場合、上の説明で全てを語ることができてしまっていて面白く無いので、以下では各量子モンテカルロ法について、その概要を説明していこうと思う。
というのも、僕がひょんなことから知った量子モンテカルロというキーワードを調べたときには、スピンモデルの研究とか、量子化学の研究とか、色々なものが混じって出てきてしまって、いったいぜんたい量子モンテカルロとは何なんだ、という気持ちになったから。初学者に全然優しくなかった。そこで記事を書いて、(自分のためにも) まとめて整理しておこうと思ったのである。
最近では量子アニーリング研究の発展と相まって、もともと量子系の計算手法として使われていた手法を転用して、古典的な組み合わせ最適化を行うという方向性もある。このような場合も、そのまま名前だけ受け継がれて量子モンテカルロ法と呼ばれることがあるので注意したい。
変分量子モンテカルロ
まずは変分量子モンテカルロ、もしくは変分モンテカルロ (Variational Monte-Carlo,
VMC) と呼ばれる手法について説明しよう。VMC は、量子状態 \(\ket{\psi}\) に対してある特定の形 \(\ket{\psi} = \ket{\psi(\b{\lambda})}\) (\(\b{\lambda}\) はパラメータ) を仮定して、変分法によってハミルトニアン \(H\) の基底状態を求める手法である。
具体的には、VMC は、パラメータ付き量子状態 \(\ket{\psi{\b{\lambda}}}\) のエネルギー期待値
\begin{align}
E(\b{\lambda}) := \bra{\psi(\b{\lambda})} H \ket{\psi(\b{\lambda})}
\end{align}
を最小化するパラメータ \(\b{\lambda}^*\) をみつけることによって、近似的な基底状態 \(\ket{\psi(\b{\lambda}^*)}\) を構成する。
これだけなら例えば
Hartree-Fock 法 などのモンテカルロ的でない方法と同じなのだが、VMC では \(\bra{\psi(\b{\lambda})} H \ket{\psi(\b{\lambda})}\) の値をモンテカルロ積分で評価するところが異なる。モンテカルロ的では無い変分法では、\(\bra{\psi(\b{\lambda})} H \ket{\psi(\b{\lambda})}\) を直接効率的に計算できる必要があるので、そのような条件が満たされる \(\ket{\psi(\b{\lambda})}\) しか扱えなかった。例えば Hartree-Fock 法では、\(\ket{\psi(\b{\lambda})}\) に \(\braket{\b{x}_1\cdots \b{x}_N}{\psi({\b{\lambda}})}\) が Slater 行列式で表せる、という制限をつけて、その中で探索を行う。
一方で VMC では、モンテカルロ積分によってエネルギー期待値 \(\bra{\psi(\b{\lambda})} H \ket{\psi(\b{\lambda})}\) が評価できる量子状態 \(\ket{\psi(\b{\lambda})}\) ならばどんな \(\ket{\psi(\b{\lambda})}\) を使っても良い。これによって表現できる量子状態の幅が格段に向上し、モンテカルロ的でない方法では探索できない \(\ket{\psi(\b{\lambda})}\) を構成できるようになる。
普通、モンテカルロ積分で \(\bra{\psi(\b{\lambda})} H \ket{\psi(\b{\lambda})}\) を評価するには次のようにする。まず、適当な基底 \(\ket{i}\) でこの式を展開し
\begin{align}
\bra{\psi(\b{\lambda})} H \ket{\psi(\b{\lambda})} &= \bra{\psi(\b{\lambda})}\left(\sum_i\ket{i}\bra{i}\right) H \ket{\psi(\b{\lambda})} \\
&= \sum_{i} \braket{\psi(\b{\lambda})}{i} \bra{i} H \ket{\psi(\b{\lambda})} \\
\end{align}
とする。さらに
\begin{align}
&= \sum_{i} \braket{\psi(\b{\lambda})}{i} \bra{i}H\ket{\psi(\b{\lambda})} \frac{\braket{i}{\psi(\b{\lambda})}}{\braket{i}{\psi(\b{\lambda})}} \\
&= \sum_{i} |\braket{i}{\psi(\b{\lambda})}|^2\frac{\bra{i}H\ket{\psi(\b{\lambda})}}{\braket{i}{\psi(\b{\lambda})}}
\end{align}
と展開する。この和は、\(p_i = |\braket{i}{\psi(\b{\lambda})}|^2\) という確率分布から \(S\) 個のインデックス \(i_{s}\) (\(s=1,\cdots,S\)) をサンプルして、そのサンプルについて \(\bra{i_s}H\ket{\psi(\b{\lambda})}/\braket{i_s}{\psi(\b{\lambda})}\) の平均を取れば計算できる。つまり \(S\) が十分大きいとき、
\begin{align}
\frac{1}{S}\sum_s \frac{\bra{i_s}H\ket{\psi(\b{\lambda})}}{\braket{i_s}{\psi(\b{\lambda})}} &\approx \sum_{i} |\braket{i}{\psi(\b{\lambda})}|^2\frac{\bra{i}H\ket{\psi(\b{\lambda})}}{\braket{i}{\psi(\b{\lambda})}}\\
& = \bra{\psi(\b{\lambda})} H \ket{\psi(\b{\lambda})}
\end{align}
である。
上の具体的計算法を見ればわかるように、普通の変分量子モンテカルロでは、
- 確率分布 \(p_i\) からのサンプリングが効率的にできる。
- \(\bra{i_s}H\ket{\psi(\b{\lambda})}/\braket{i_s}{\psi(\b{\lambda})}\) が効率的に計算できる。
の2点が達成されるような \(\ket{\psi(\b{\lambda})}\) を使う必要がある。(もちろん、上の方法以外の方法によって \(\bra{\psi(\b{\lambda})} H \ket{\psi(\b{\lambda})}\) が効率的に評価できれば別にそれでも問題ない。)
そのような条件を満たす \(\ket{\psi(\b{\lambda})}\) の形には色々なものがあるようだ。ちょっと調べてみると、例えば量子化学や固体の計算、つまり電子系でよく使われるものとして、スレーター-ジャストロー (Slater-Jastrow) 型波動関数と呼ばれるものがあるらしい。これは \(D(\b{x})\) をスレーター行列式で表される波動関数として、それに Jastrow 因子と呼ばれる項 \(e^{J(\b{x})}\) をかけた、
\begin{align}
\psi_{SJ}(\b{x}) &:= e^{J(\b{x})} D(\b{x})\\
J(\b{x}) &= \sum_{i}\chi(\b{x}_i) + \sum_{i,j} u(\b{x}_i, \b{x_j})
\end{align}
という形の波動関数である。最適化されるパラメータ \(\lambda\) は、\(\chi\) や \(u\) と書いた関数の中に埋め込まれている。
\(u\) や \(\chi\) の形も研究されているらしく、なかなか奥が深そうな感じがする。
また、最近ではニューラルネットワークを使って \(\ket{\psi(\b{\lambda})}\) を構成するというのも流行ってるみたいだ。
変分量子モンテカルロ以外の量子モンテカルロ
長々と VMC について説明してしまったが、ここからはそれ以外の量子モンテカルロ法について見ていこう。"量子モンテカルロ法" とか、英語で "quantum monte carlo" と検索してみると、(VMCももちろん出てくるがそれ以外に)
- 経路積分量子モンテカルロ (path-integral QMC)
- 拡散モンテカルロ (Diffusion Monte-Carlo, DMC)
- グリーン関数量子モンテカルロ (Green's function QMC)
という言葉がよく出てくる。
初めて量子モンテカルロについて調べているときには、名前が違うのだから全てちょっとずつ異なる手法なのだろう、と思っていたが、調べていくと、
これら 3 つは全て本質的に同じである ようだった。どうも初学者にやさしくない。
多分彼らは分野によって違う言葉を使っている。スピン模型などをやっているひとは経路積分モンテカルロという言葉を使い、量子化学や固体の第一原理計算をやっている人は拡散モンテカルロという言葉を使っているみたいだ。 記事公開後に頂いたコメントによると、どうも分野によって、というよりは実装方法による違いで呼び分けられているようだった。違いは後述する。グリーン関数量子モンテカルロという言葉は、どうもこの言葉を使っている新しい文献が見つからないので、少し古い呼び方なのだと思う。
虚時間発展
さて、手法の具体的内容に入ろう。これらは全て、
虚時間発展
を行うためのアルゴリズムだ。虚時間発展とは、ハミルトニアン \(H\) と (正の) 実数 \(\beta\) について \(e^{-\beta H}\) を量子状態に作用させることである。量子状態の実時間発展は \(e^{-iHt}\) と表されるが、この時間 \(t\) を \(t=-i\beta\) とすることでこの演算子 \(e^{-\beta H}\) になることから「虚時間発展」と呼ばれる。
虚時間発展を使うと、量子系の基底状態を求めたり、有限温度での平衡状態 \(e^{-\beta H}\) を計算したりできる。ここでは基底状態を求める、という応用に的を絞って説明することにする。
実は適当な初期量子状態 \(\ket{\psi(0)}\) を虚時間発展させることで、その量子系の基底状態が得られるのだ。このことを示すには、次のようにする。まずハミルトニアン \(H\) の固有状態を \(\ket{E_i}\)、対応する固有エネルギーを \( 0\leq E_0\leq E_1 \leq \cdots \) としよう。基底状態のエネルギー \(E_0\) が0以上であると仮定したが、もし0以下になってしまうときは、適当に大きな正の数 \(a\) を取ってきて \(H \to H+a\) とすれば良い。すると初期状態は必ず
\begin{align}
\ket{\psi(0)} = \sum_{i} c_i(0)\ket{E_i}
\end{align}
と展開できる。この状態に \(e^{-\beta H}\) を作用させると、
\begin{align}
\ket{\psi(\beta)} &= e^{-\beta H}\left(\sum_{i} c_i(0)\ket{E_i}\right) \\
&= \sum_{i} c_i(0)e^{-\beta E_i}\ket{E_i}
\end{align}
となる。この式から、虚時間発展を行うとエネルギーの高い成分を指数関数的に消せることがわかる。実際 \(c_i(\beta) = c_i(0)e^{-\beta E_i}\) とすると
\begin{align}
\frac{c_i(\beta)}{c_0(\beta)} = e^{-\beta (E_i-E_0)}
\end{align}
なので、\(\beta\) を十分大きく取れば、\(\ket{\psi(\beta)}\) には近似的に \(\ket{E_0}\) だけが含まれるようになる。
虚時間発展を確率的に行う
この虚時間発展を確率的に行うのが上記3つの量子モンテカルロ法である。この節ではこれをどうやって実行しているのか説明する。
ハミルトニアン \(H\) は普通、適当な演算子の和として表されていることが多い。そのことを明示的に示すために、一般に
\begin{align}
H = \sum_{\alpha=1}^M H_\alpha
\end{align}
と書こう。作り出したい演算子は \(e^{-\beta H}\) であるが、\(H\) が指数的な大きさの行列であることから、これを直接計算するのは計算量的にしんどい。一方でハミルトニアンを構成している要素 \(H_\alpha\) であれば、\(e^{-\beta H_\alpha}\) を厳密に計算できることが多い。
そこで、これらの量子モンテカルロアルゴリズムでは鈴木-トロッター展開と呼ばれる近似公式
\begin{align}
e^{-\Delta\beta H} = e^{-\Delta\beta H_1}e^{-\Delta\beta H_2}\cdots e^{-\Delta\beta H_M} + O\left(\Delta\beta^2 \left(\sum_{\alpha}\|H_\alpha\|\right)^2\right)
\end{align}
を用いて、所望の演算 \(e^{-\beta H}\) を実現する。(誤差項は
[arXiv:1912.08854] を参考にした。) 上式において、\(\Delta\beta\) は小さな実数である。つまり、\(\Delta\beta = \beta/n\) として
\begin{align}
e^{-\beta H} \approx \left[e^{-\beta H_1/n}e^{-\beta H_2/n}\cdots e^{-\beta H_M/n}\right]^n
\end{align}
によって \(e^{-\beta H}\) を近似する。
さて、上記の量子モンテカルロ法では、上の近似的な居時間発展演算子を初期状態 \(\ket{\psi_0}\) に作用させた量
\begin{align}
\left[e^{-\beta H_1/n}e^{-\beta H_2/n}\cdots e^{-\beta H_M/n}\right]^n \ket{\psi_0}
\end{align}
を計算する際に、確率的過程を使う。アイデアは、各 \(e^{-\beta H_\alpha/n}\) の間に適当な正規直交基底 \(\ket{i}\) で作った恒等演算子 \(\sum_{i} \ket{i}\bra{i}\) を挿入するというものである。やってみると、
\begin{align}
&\left(\sum_{i} \ket{i}\bra{i}\right)e^{-\beta H_1/n}\left(\sum_{j} \ket{j}\bra{j}\right)e^{-\beta H_2/n}\left(\sum_{k} \ket{k}\bra{k}\right)\cdots \left(\sum_{l} \ket{l}\bra{l}\right)e^{-\beta H_M/n}\left(\sum_{m} \ket{m}\bra{m}\right)\ket{\psi_0}\\
&=\sum_{i,j,k,\cdots,l,m} \ket{i}\bra{i}e^{-\beta H_1/n}\ket{j}\bra{j}e^{-\beta H_2/n}\ket{k}\bra{k}\cdots \ket{l}\bra{l}e^{-\beta H_M/n}\ket{m}\braket{m}{\psi_0}
\end{align}
となる。(\(\cdots\) の中で \(e^{-\beta H_1/n}e^{-\beta H_2/n}\cdots e^{-\beta H_M/n}\) が \(n\) 回繰り返されていると思ってほしい。) ここで現れた \(\bra{i}e^{-\beta H_\alpha/n}\ket{j}\) という形の項を
\begin{align}
G_{ij}^{(\alpha)} := \bra{i}e^{-\beta H_\alpha/n}\ket{j}
\end{align}
とおく。すると上の式は
\begin{align}
\sum_{i,j,k,\cdots,l,m}\ket{i}G^{(1)}_{ij}G^{(2)}_{jk}\cdots G^{(M)}_{lm}\braket{m}{\psi_0}
\end{align}
となり、\(G^{(\alpha)}_{ij}\) という行列の掛け算を繰り返し実行するものだということがよく分かる。もちろん、この行列積を直接計算しようとすると指数的な計算量が必要となってしまう。そこで確率遷移行列とみなせる行列 \(P^{(\alpha)}_{ij}\) を次のように定義する。
\begin{align}
P^{(\alpha)}_{ij} := \frac{|G^{(\alpha)}_{ij}|}{\sum_{i}|G^{(\alpha)}_{ij}|}
\end{align}
逆に \(G^{(\alpha)}_{ij}\) は \(P^{(\alpha)}_{ij}\) を用いて
\begin{align}
G^{(\alpha)}_{ij} &= w_{ij} P^{(\alpha)}_{ij}\\
w_{ij} &= \frac{G^{(\alpha)}_{ij}}{|G^{(\alpha)}_{ij}|}\sum_{i}|G^{(\alpha)}_{ij}|
\end{align}
と表せる。すると上の式は、
\begin{align}
\sum_{i,j,k,\cdots,l,m} w_{ij}w_{kl}\cdots w_{lm} \ket{i}P^{(1)}_{ij}P^{(2)}_{jk}\cdots P^{(M)}_{lm}\braket{m}{\psi_0}
\end{align}
となる。さらに \(\braket{m}{\psi_0} = c_m\) と書こう。ここまでくると、(賢い方ならもっと前でわかっていたかもしれないが) この和がモンテカルロ的に評価できることがわかるだろう。次のようなアルゴリズムを使えば良い。まず、
- \(|c_m|\) に比例する確率 (具体的には \(p_m = |c_{m}|/(\sum_{m}|c_{m}|)\)) で \(m\) をサンプルする。
- \(P_{lm}^{(M)}\) に比例する確率で \(l\) をサンプルする。
- \(P_{l'l}^{(M-1)}\) に比例する確率で \(l'\) をサンプルする。
- \(\vdots\)
- \(P_{jk}^{(2)}\) に比例する確率で \(j\) をサンプルする。
- \(P_{ij}^{(1)}\) に比例する確率で \(i\) をサンプルする。
として、基底のインデックスの系列
\begin{align}
m \to l \to \cdots \to k \to j \to i
\end{align}
を得る。さらにこれを \(S\) 回繰り返して、多数の系列
\begin{align}
m_s \to l_s \to \cdots \to k_s \to j_s \to i_s~~(s=1,\cdots,S)
\end{align}
を得た後、
\begin{align}
\sum_{s=1}^S w_{i_sj_s}w_{k_sl_s}\cdots w_{l_sm_s} c_{m_s}' \ket{i_s}
\end{align}
(ただし \(c_{m}' = \frac{c_m}{|c_m|}\sum_{m}|c_m|\)) という量子状態を構成する。この和は、大数の法則により
\begin{align}
\sum_{s=1}^S w_{i_sj_s}w_{k_sl_s}\cdots w_{l_sm_s} c_{m_s}' \ket{i_s}\\
\to \sum_{i,j,k,\cdots,l,m} w_{ij}w_{kl}\cdots w_{lm} \ket{i}P^{(1)}_{ij}P^{(2)}_{jk}\cdots P^{(M)}_{lm}\braket{m}{\psi_0}
\end{align}
と収束していくので、これで虚時間発展をモンテカルロ的に行えたことになる。
各量子モンテカルロ法の名前に関して
各量子モンテカルロ法の名前の由来と、実装による呼び分けについて書いて終わりにしようと思う。(僕の想像も入っているので間違っていたら教えてほしい。)
グリーン関数量子モンテカルロ:
\(G^{(\alpha)}_{ij}\) がグリーン関数、あるいはプロパゲータと呼ばれているものだから。実装方法はわからなかった。(そもそもこの方法は最近はあまり使われてないようだ。)
経路積分モンテカルロ:
途中で出てきた
\begin{align}
\sum_{i,j,k,\cdots,l,m} \ket{i}\bra{i}e^{-\beta H_1/n}\ket{j}\bra{j}e^{-\beta H_2/n}\ket{k}\bra{k}\cdots \ket{l}\bra{l}e^{-\beta H_M/n}\ket{m}\braket{m}{\psi_0}
\end{align}
という式が、\(\ket{m}\) から始まって \(\ket{i}\) に至るまで、全ての取りうる経路について足し合わせたように見えるから。 全体の \(\beta\) をある値に固定して、上記の和をモンテカルロ法によって求めるという実装をするものを特にこの名前で呼ぶようだ。
拡散モンテカルロ:
電子のハミルトニアン
\begin{align}
H_{\text{elec}} = \sum_i \left(-\frac{\nabla_i^2}{2} + V(\b{r}_i)\right) + \sum_{ij} \frac{1}{r_{ij}}
\end{align}
に対して Trotter 展開を適用すると
\begin{align}
\exp(-\beta H_{\text{elec}}) \approx \exp\left(-\beta\sum_i\frac{\nabla_i^2}{2}\right) \exp\left[\sum_i V(\b{r}_i) + \sum_{ij} \frac{1}{r_{ij}}\right]
\end{align}
のようになる。このとき現れる \(\exp\left(-\beta\sum_i\frac{\nabla_i^2}{2}\right)\) は拡散方程式
\begin{align}
\frac{\partial}{\partial \beta} f(\b{r}_1, \cdots ,\b{r}_N, \beta) = -\frac{\nabla_i^2}{2} f(\b{r}_1, \cdots ,\b{r}_N, \beta)
\end{align}
のプロパゲータであり、つまり \(\exp\left(-\beta\sum_i\frac{\nabla_i^2}{2}\right)\) はある点 \((\b{r}_1, \cdots, \b{r}_N)\) から \(\beta\) の時間だけ拡散させる、という作用を表す演算子である。これが理由となって、特に電子系などの量子力学的粒子系を実空間で虚時間発展させる手法を、拡散モンテカルロと読んでいるようだ。ちなみに具体的な行列表現は
\begin{align}
\bra{\b{r}'}\exp\left(-\beta\sum_i\frac{\nabla_i^2}{2}\right)\ket{\b{r}} = (2\pi \beta)^{-3N/2}\exp\left(-\frac{|\b{r}-\b{r}'|^2}{2\beta}\right)
\end{align}
である。特に、全体の \(\beta\) を固定せず、逆に固定した \(\Delta\beta = \beta/n\) を使って \(\exp(-\Delta \beta H_i)\) をかけ続け、確率分布が平衡状態に達するまで時間発展を続けるという実装をするものをこの名前で呼ぶようだ。