変分モンテカルロ(VMC)をやってみる~準備~

Thijssenの計算物理学の12.2.2の部分の練習問題をやってみる。

www.maruzen-publishing.co.jp

変分モンテカルロ法(VMC)
VMCのための前準備
参考文献

変分モンテカルロ法(VMC)

量子力学の多くの問題においてSchrödinger方程式を厳密に解くことは難しく、何らかの形で近似的に解く必要がある。そこで変分による方法で基底状態について近似的に解くことを考えてみよう。

VMCのための前準備

量子力学での変分原理

Hamiltonianを $H$ とし、そのエネルギー固有値を $E_n$ としたときにそれに対応した固有関数を $\phi_n(X)$ （ $X$ は座標）と書くことにしよう*1。
ここで唐突だが、ある波動関数 $\Psi(X,a)$ （ $a$ はパラメータ）を考えてみよう。この波動関数をエネルギー固有関数でモード展開すると次のように展開できる：

$\Psi(X,a)=\displaystyle\sum_{n}\Psi_n\phi_n(X)$

ここで $\Psi_n$ は $\phi_n(X)$ で展開するときの展開係数であり、それは

$\Psi_n=\int dX\, \phi_n^{*}(X)\Psi(X,a)$

で与えられる。
量子力学においてエネルギーの期待値 $\langle E\rangle$ は、

$\langle E\rangle=\langle \Psi|H|\Psi\rangle=\int dX\,\Psi^{*}(X,a)H\Psi(X,a)$

で計算できるのであった。ここにモード展開された表式での波動関数を代入してあげると、

$\displaystyle \langle E\rangle =\int dX\,\sum_{n,m}\Psi_{m}^{*}\Psi_{n}\phi_m^{*}(X)H\phi_{n}(X)\\ =\displaystyle\sum_{n,m}\Psi_{m}^{*}\Psi_{n}\int dX\,\phi_m^{*}(X)E_{n}\phi_{n}(X)\\ =\displaystyle\sum_{n,m}E_{n}\Psi_{m}^{*}\Psi_{n}\delta_{n,m}\\ =\displaystyle\sum_{n}E_n|\Psi_n|^2\\ \geq\displaystyle\sum_{n}E_{0}|\Psi_n|^2\\ =E_0$

となる。ここで $\Psi(X,a)$ が規格化されているとして規格化条件 $\sum_{n}|\Psi_n|^2=1$ を用いた。
ここから波動関数の期待値 $\langle E\rangle$ は基底状態のエネルギー $E_0$ を下回ることはないことがわかる。なので $\Psi(X,a)$ のパラメータ $a$ を期待値 $\langle E\rangle$ が最小になるように調整すれば波動関数は基底状態の波動関数に近づくことになり、基底状態のエネルギーの近似値を得ることができる。

量子力学における変分原理

試行波動関数 $\Psi(X,a)$ （ $a$ は変分パラメータ）について、

$\delta\left(\int dX\,\Psi^*(X,a)H\Psi(X,a)\right)=0$

となるように変分パラメータを決めると基底状態の良い近似になる。

さて、ここまでは試行波動関数 $\Psi(X,a)$ が規格化されているものとして話を進めてきたが、一般に試行波動関数として複雑なものを採用した場合は規格化を行うことは面倒であるし、そもそも積分を解析的に行うことが可能かどうかすら怪しいであろう。そこでこの先では試行波動関数があらかじめには規格化されていないとし、試行波動関数 $|\Psi(a)\rangle$ を $(||\Psi(a)||\rangle)^{-1/2}|\Psi(a)\rangle$ と書くことにしよう。すると、エネルギー期待値は、

$\displaystyle \langle E\rangle=\frac{\langle\Psi(a)|H|\Psi(a)\rangle}{\langle\Psi(a)|\Psi(a)\rangle}$

と書き直される。よって変分原理は、

量子力学における変分原理

試行波動関数 $\Psi(X,a)$ （ $a$ は変分パラメータ）について、

$\delta\left(\displaystyle\frac{\int dX\,\Psi^*(X,a)H\Psi(X,a)}{\int dX\,\Psi^*(X,a)\Psi(X,a)}\right)=0$

となるように変分パラメータを決めると基底状態の良い近似になる。

のように書き直せばよい。
原理的にはこの変分原理に従ってエネルギーが最小になるように変分パラメータ $a$ を調整すればよいのだが、先に述べたように試行関数やHamiltonianが複雑になれば解析的に実行することは困難であり、積分を区分求積的に数値計算で行うにしても1次元系ならともかく3次元などになれば問題が複雑になればなるほど積分の次元が飛躍的に大きくなり次元の呪いに嵌ってしまう。例えば3次元の2粒子系（Helium原子の電子など）では積分に必要な次元は6次元になる。
そこで高次元でも数値計算を行うためにマルコフ連鎖モンテカルロ法を用いて変分における積分計算を行うことにする（だから変分"モンテカルロ"という。）。

マルコフ連鎖モンテカルロ法(MCMC)

前節での要請に従ってMCMCを導入する。

モンテカルロ法

先にも述べたように高次元の数値積分は難しい。このことを実感するために領域として長さ1の $n$ 次元超立方体を積分することを考えてみよう。仮に $[0,1]$ の間を1000分割して区分求積的に積分を行うことを考えると、積分領域は $10^{3n}$ 個の領域に分けられることになり、考えなければならない領域の数が指数関数的に増加してしまい高次元では計算不可能になってしまうことがわかるだろう。
この計算コストの指数的増加による次元の呪いから逃れるために積分領域からランダムに値を選んできてそれらの平均をとることで積分を求めるということを考えてみる。同様に $[0,1]^n$ の領域での積分を考えよう。ランダムに点を選んで持ってきて平均を取るというアイデアに従えば積分は、

$I=\displaystyle\frac{1}{N}\sum_{i=1}^{N}f(x_i)\hspace{10pt}(x_i\mbox{はランダムに選ぶ})$

と表される。この式自体は $x_i$ をランダムに選ぶことを除けば区分求積とほぼ同じ形なので $N$ が大きくなれば厳密な積分に近づくことは容易にわかる。そこでこの近づき方がどのように振る舞うかを確認してみよう。そのためにはこの式の分散の振る舞いを見ればよい。 $I$ の分散は次の式で表せる。

$\sigma^2=\displaystyle\left\langle\left(\frac{1}{N}\sum_{i=1}^{N}f(x_i)\right)^2\right\rangle-\left(\left\langle\frac{1}{N}\sum_{i=1}^{N}f(x_i)\right\rangle\right)^2\\$

この式において初項の $i\neq j$ の時は後ろの項と相殺して消え、 $i=j$ のときに初項から $f^2$ の期待値が $N$ 個でることを考えて計算すると、

$\sigma^2=\displaystyle\frac{1}{N}\left(\left\langle f^2\right\rangle-\left\langle f\right\rangle^2\right)$

となり、 $\sigma\sim1/\sqrt{N}$ のように次元によらない形で振る舞うことがわかる。
以上からモンテカルロ法を用いて積分計算をすれば次元の呪いを回避しながら数値計算できそうだと予想できるだろう。

マルコフ連鎖モンテカルロ法

さて、上の例ではランダムに点を選ぶときに点が選ばれる確率は一様にサンプリングされるとした。しかし、高次元の積分を行うときには一様にサンプリングすることは計算資源の無駄遣いになってしまう。というのも多くの被積分関数は領域によって積分への寄与が全く異なり、次元が上がれば積分への寄与が大きい適切なサンプリング点がサンプリングされる確率が下がってしまうからである。例えばGauss積分 $e^{-\sum_{i=1}^{n}x_i^2/2}$ を考え、サンプリングは $[-10000,10000]^n$ から行うとする。よく知られているようにGauss関数は $0$ に近いところでのみ大きな値を持ち原点から離れるにしたがって急速に $0$ に近づく。そこで有意な値を持つ範囲を $|x|>2$ であるとすると、 $n=1$ のときの有意なサンプリング領域は $2\times10^{-2}\%$ である。しかし、 $n=10$ では有意な領域は全積分領域のうち $4\times10^{-6}\%$ しかない。こうしたことから高次元で安易に一様サンプリングをしてしまうと大きな統計誤差を含んだ結果が出てしまう。つまり一見モンテカルロ法を用いて次元の呪いを回避できたように見えてサンプリングの問題として残ってしまっている。
このようなサンプリングの問題を踏まえると、サンプリング数を高次元でも膨大なものにしないために一様にサンプリングを行うのではなく、積分領域のうち有意な領域を重点的にサンプリングすればいいだろう。この有意な領域への重点的なサンプリングを行うために望んでいる確率分布を構成する方法の一つとしてマルコフ連鎖モンテカルロ法(MCMC)がある。

構成したい確率分布を $P(\{x\})=P(x_1,\dots,x_n)$ としよう。MCMCでは $x_1, \dots, x_n$ の値を各点が現れる滞在時間を $P(\{x\})$ に比例するように変化させるようにし、これによって十分長い時間で統計平均の良い近似を得ることができる。そこでこのようなMCMCの方法として今回用いるメトロポリス法によるMCMCを導入することにしよう*2。

メトロポリス法

MCMCの一つの方法としてメトロポリス法の手順を確認する。
先と同様に構成したい確率分布を $P(x)$ とする。メトロポリス法の手続きは次のようなものである。

メトロポリス法

1. $k$ 番目のサンプリング点を $x^{(k)}$ とする。実数 $\Delta x$ をランダムに選び、 $k+1$ 番目の候補を $x'=x^{(k)}+\Delta x$ として提案する。ただし、MCMCの条件を満たすために $\Delta x$ :と $-\Delta x$ が同じ確率で現れるように候補の選び方を提案する。
2. メトロポリステスト:
$[0,1]$ の間の一様乱数 $r$ を生成し、 $r< P(x')/P(x)$ であれば提案を受理して $x^{(k+1)}=x'$ とする。そうでなければ候補を棄却して $x^{(k+1)}=x^{(k)}$ とする(この条件は $\mathrm{min}(1, P(x')/P(x))$ で提案を受理するということと同じである。)。

1において述べたように $\Delta x$ と $\Delta -\Delta x$ が出現する確率が同じであればよいのでその提案確率は $[-c,c]$ の間の一様乱数でも $x^{k}$ を中心としたガウス乱数でも何でもよい。なお、記事でVMCをする際にはガウス乱数を用いて計算することにする。

VMC

このようにして得られるメトロポリス法によるMCMCを用いてどのように変分計算を行うかを考えてみよう。
量子力学のおける変分原理によれば、試行関数を最適化するためには、

$\langle E\rangle =\displaystyle\frac{\int dX\,\Psi^*(X,a)H\Psi(X,a)}{\int dX\,\Psi^*(X,a)\Psi(X,a)}$

を最小化すればよいのであった。唐突だがこの式を少し変形してみよう。形を確率分布を挟んだ積分らしく変形すると、 $\langle E\rangle$ は局所エネルギー $E_L$ を用いて

$\langle E\rangle =\displaystyle\frac{\int dX\,\Psi^*(X,a)\Psi(X,a)\displaystyle\frac{H\Psi(X,a)}{\Psi(X,a)}}{\int dX\,\Psi^*(X,a)\Psi(X,a)}\\ =\displaystyle\int dX\,\frac{|\Psi(X,a)|^2}{\int dX'\,|\Psi(X',a)|^2}\displaystyle\frac{H\Psi(X,a)}{\Psi(X,a)}\\ =\int dX\,\rho(X)E_L$

と書くことができる。ここで $\rho(X)$ と $E_L$ について

$\rho(X)= \displaystyle\int dX\,\frac{|\Psi(X,a)|^2}{\int dX'\,|\Psi(X',a)|^2}$

$E_L= \displaystyle\frac{H\Psi(X,a)}{\Psi(X,a)}$

とした。
このように書き換えることでエネルギーの期待値 $\langle E\rangle$ は確率分布が $\rho(X)$ である局所エネルギー $E_L$ の期待値と思うことができる。そこでMCMCを用いてサンプリング点を $\rho(X)$ から選び出して平均を取ってやれば大数の法則からエネルギーの期待値 $\langle E\rangle$ を求めることができるであろう。これがVMCにおける多次元積分の基本的な方針である。
最後にここまでの話を踏まえてVMCでのメトロポリス法を整理しておこう。

VMC

**REPAT**
始めの変分パラメータ $a$ を決める。ランダムな場所に $N$ 個のwalkerを初期配置として配置する。
**REPEAT**
**FOR 全walker**

1. $k$ 番目のサンプリング点を $X^{(k)}$ とする。walkerの移動後の座標として $k+1$ 番目の候補を提案分布に従って $X'=X^{(k)}+\Delta X$ として提案する。ただし、MCMCの条件を満たすために $\Delta X$ :と $-\Delta X$ が同じ確率で現れるように候補の選び方を提案する。
2. メトロポリステスト:
$[0,1]$ の間の一様乱数 $r$ を生成し、 $r< \rho(X')/\rho(X)=|\Psi(X')/\Psi(X)|^2$ であれば提案を受理して $X^{(k+1)}=X'$ とする。そうでなければ候補を棄却して $X^{(k+1)}=X^{(k)}$ とする。
3. サンプリング点での局所エネルギー $E_L$ を計算する。

**END FOR**
得られた局所エネルギーの平均を取ってエネルギーの期待値などを求める。
Newton法に従って変分パラメータ $a$ を更新する。
**END**

Newton法

最後に変分パラメータの調整の仕方を確認する。
エネルギーが最小になる $\alpha$ では $\langle E\rangle$ の $a$ による微分が $0$ になるので $dE/da$ のゼロ点を求めればよい。 $E$ 自体を数値微分したものを使ってもよいが二階導関数を数値微分の数値微分で求めることになるので大きな誤差が出てしまう。そこで一階導関数は解析的に計算することにしよう。 $E$ の表式に戻って微分を行うと、

$\displaystyle\frac{d\langle E\rangle}{da} =\displaystyle\frac{d}{da}\frac{\langle\Psi(a)|H|\Psi(a)\rangle}{\langle\Psi(a)|\Psi(a)\rangle}\\ =\displaystyle\frac{d\langle \Psi(a)|}{da}\frac{H|\Psi(a)\rangle}{\langle \Psi(a)|\Psi(a)\rangle}+\frac{d}{da}\frac{\langle\Psi(a)|H}{\langle\Psi(a)|\Psi(a)\rangle}\frac{d|\Psi(a)\rangle}{da}\\ -\displaystyle\frac{\langle\Psi(a)|H|\Psi(a)\rangle}{\left(\langle \Psi(a)|\Psi(a)\rangle\right)^2}\left(\frac{d\langle \Psi(a)|}{da}|\Psi(a)\rangle+\langle \Psi(a)|\frac{d|\Psi(a)\rangle}{da}\right)$

となる。ここで試行関数として実関数のみを用いることを仮定すると、計算結果として現れるのは実数であることに注意すると、

$\displaystyle\langle \Psi(a)|H\frac{d|\Psi(a)\rangle}{da}=\left(\langle \Psi(a)|H\frac{d|\Psi(a)\rangle}{da}\right)^{*}=\frac{d\langle \Psi(a)}{da}|H\Psi(a)\rangle$

$\displaystyle\langle \Psi(a)|\frac{d|\Psi(a)\rangle}{da}=\left(\langle \Psi(a)|\frac{d|\Psi(a)\rangle}{da}\right)^{*}=\frac{d\langle \Psi(a)}{da}|\Psi(a)\rangle$

と右側に落とした微分を左側に持ってくることができる。よって $dE/da$ は、

$\displaystyle\frac{d\langle E\rangle}{da}\\ =\displaystyle2\left(\frac{d\langle \Psi(a)|}{da}\frac{H|\Psi(a)\rangle}{\langle \Psi(a)|\Psi(a)\rangle}-\displaystyle\frac{\langle\Psi(a)|H|\Psi(a)\rangle}{\left(\langle \Psi(a)|\Psi(a)\rangle\right)^2}\frac{d\langle \Psi(a)|}{da}|\Psi(a)\rangle\right)\\ =\displaystyle 2\left(\int dX\frac{\Psi(X,a)H\Psi(X,a)}{\int dX'\Psi(X',a)^2}\frac{1}{\Psi(X,a)}\frac{d\Psi(X,a)}{da}\right.\\ \left.\displaystyle-\int dX\frac{\Psi(X,a)H\Psi(X,a)}{\int dX'\Psi(X',a)^2}\int dX\frac{\Psi(X,a)^2}{\int dX'\Psi(X',a)^2}\frac{1}{\Psi(X,a)}\frac{d\Psi(X,a)}{da}\right)\\ =\displaystyle 2\left(\int dX\frac{\Psi(X,a)H\Psi(X,a)}{\int dX'\Psi(X',a)^2}\frac{d\ln \Psi(X,a)}{da}\right.\\ \left.\displaystyle-\int dX\frac{\Psi(X,a)H\Psi(X,a)}{\int dX'\Psi(X',a)^2}\int dX\frac{\Psi(X,a)^2}{\int dX'\Psi(X',a)^2}\frac{d\ln \Psi(X,a)}{da}\right)\\ =\displaystyle 2\left(\left\langle E_L\frac{d\ln\Psi(a)}{da}\right\rangle-\left\langle E_L\right\rangle\left\langle\frac{d\ln \Psi(a)}{da}\right\rangle\right)$

と計算できる。これをMCMCを用いて数値計算し、あとはNewton法の更新式

$\displaystyle a_{\mathrm{new}}=a_{\mathrm{old}}-\left.\gamma\frac{d\langle E\rangle}{da}\right|_{a=a_{\mathrm{old}}}$

に従って変分パラメータ $a$ を更新すればよい。ただし $\gamma$ は $dE/da$ の数値微分である。

さて、これで準備は終わりなので以上を踏まえて具体的な系で計算をしてみよう。

参考文献

[1]. J. M. Thijssen, 「計算物理学」,丸善出版, 2012.
[2]. 花田政範, 松浦壮, 「ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門」, 講談社, 2020.

*1: $n$ が大きくなる順にエネルギー固有値が並べられているとしよう。

*2:もちろんMCMCになるような一般条件があるのだが、本題ではないのでこの記事では述べない。

MiyanTarumi’s blog

感想とか色々。数物系は書くか微妙。。。

変分モンテカルロ(VMC)をやってみる~準備~

変分モンテカルロ法(VMC)

VMCのための前準備

量子力学での変分原理

マルコフ連鎖モンテカルロ法(MCMC)

モンテカルロ法

マルコフ連鎖モンテカルロ法

メトロポリス法

VMC

Newton法

参考文献