アルゴリズム

変分モンテカルロ法

変分モンテカルロ法では, 試行波動関数を用意して, その試行波動関数が含むパラメータを 変分原理に従って最適化することで量子多体系の基底状態(または低励起エネルギー状態)の波動関数を 近似的に求めます。 試行波動関数に対する物理量の期待値を計算する部分で, マルコフ連鎖モンテカルロ法を利用し, 効率よく重み付きサンプリングを行います。

本パッケージでは, サンプリングに用いる完全系として電子の実空間配置 \(| x\rangle\) をとっています:

\[| x\rangle = \prod_{n=1}^{N_e/2} c_{r_{n\uparrow}}^{\dagger} \prod_{n=1}^{N_e/2} c_{r_{n\downarrow}}^{\dagger} |0 \rangle\]

ここで, \(r_{n\sigma}\) は \(n\) 番目の電子(スピン \(\sigma\))の位置, \(c_{r_{n\sigma}}^{\dagger}\) はその位置での 電子(スピン \(\sigma\))の生成演算子を表します。この基底を用いると, 演算子 \(A\) の期待値は

\[\langle A \rangle =\frac{\langle \psi| A| \psi \rangle}{\langle \psi | \psi \rangle} =\sum_x \frac{\langle \psi| A | x\rangle \langle x| \psi \rangle}{\langle \psi |\psi \rangle}\]

となるため, マルコフ連鎖の重みを

\[ \begin{align}\begin{aligned}\rho(x)=\frac{|\langle x| \psi \rangle|^2}{\langle \psi | \psi \rangle} \ge 0, \quad \sum_{x} \rho(x)=1\\と定義して,\end{aligned}\end{align} \]

\[\langle A \rangle =\sum_x \rho(x) \frac{\langle \psi| A | x\rangle }{\langle \psi |x \rangle}\]

と書き直した後、\(x\) に関する和をマルコフ連鎖モンテカルロ法により 評価しています。Local Green’s function \(G_{ij\sigma\sigma'}(x)\) は

\[G_{ij\sigma\sigma'}(x)=\frac{\langle \psi | c_{i\sigma}^{\dagger} c_{j\sigma'} | \psi \rangle}{\langle \psi | x \rangle}\]

と定義されますが, これも演算子 \(A\) を \(c_{i\sigma}^{\dagger} c_{j\sigma'}\) ととることで, 同じ方法により重み付きサンプリングを行うことができます。 なお, サンプリングに使用する乱数生成については, メルセンヌツイスター法を使用しています[Mutsuo2008 ]。

Bogoliubov表現

スピン系の計算において一体項( transfer), InterAll 形式での相互作用, 相関関数のインデックスの指定にはBogoliubov表現が使われています。 一般に、スピンの演算子は次のようにフェルミオンの生成 \(\cdot\) 消滅演算子 \(c_{i \sigma}^\dagger\), \(c_{i \sigma}\) によって書き換えることができます:

\[\begin{split}\begin{aligned} S_{i z} &= \sum_{\sigma = -S}^{S} \sigma c_{i \sigma}^\dagger c_{i \sigma} \\ S_{i}^+ &= \sum_{\sigma = -S}^{S-1} \sqrt{S(S+1) - \sigma(\sigma+1)} c_{i \sigma+1}^\dagger c_{i \sigma} \\ S_{i}^- &= \sum_{\sigma = -S}^{S-1} \sqrt{S(S+1) - \sigma(\sigma+1)} c_{i \sigma}^\dagger c_{i \sigma+1} \end{aligned}\end{split}\]

本パッケージでは、 \(S=1/2\) のスピン系のみ取り扱っており、上記の式で \(S=1/2\) と置いたものを用いています。

パフィアン-スレーター行列式の性質

この節では, パフィアン-スレーター行列式のもつ性質について簡単にまとめます。 次節 と 次々節 でパフィアン-スレーター行列式と単一スレーター行列式の間の関係を導出し、 最後 に \(f_{ij}\) の特異値分解の意味について説明します。

\(f_{ij}\) と \(\Phi_{in\sigma}\) の関係 (スピン反平行の場合)

多変数変分モンテカルロ法で試行波動関数の一体部分として用いられるパフィアン-スレーター行列式は

\[|\phi_{\rm Pf}\rangle=\Big(\sum_{i,j=1}^{N_{s}}f_{ij} c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\Big)^{N_{\rm e}/2}|0\rangle,\]

のように定義されます。ここで, \(N_{s}\) はサイト数, \(N_{e}\) は全電子数, \(f_{ij}\) は変分パラメータです。 簡単化のため, 以降 \(f_{ij}\) は実数と仮定します。また, 単一スレーター行列式として

\[\begin{split}\begin{aligned} |\phi_{\rm SL}\rangle&=\Big(\prod_{n=1}^{N_{e}/2}\psi_{n\uparrow}^{\dagger}\Big) \Big(\prod_{m=1}^{N_{e}/2}\psi_{m\downarrow}^{\dagger}\Big)|0\rangle, \\ \psi_{n\sigma}^{\dagger}&=\sum_{i=1}^{N_{s}}\Phi_{in\sigma}c^{\dagger}_{i\sigma}. \end{aligned}\end{split}\]

を定義します。ただし, \(\Phi\) は正規直交基底であり, クロネッカーのデルタ \(\delta_{nm}\) を用い

\[\sum_{i=1}^{N_{s}}\Phi_{in\sigma}\Phi_{im\sigma}=\delta_{nm},\]

で表されます。この直交性の関係から, 以下の関係式

\[\begin{split}\begin{aligned} {[\psi^{\dagger}_{n\sigma}, \psi_{m\sigma}]}_{+}&=\delta_{nm},\\ G_{ij\sigma}=\langle c_{i\sigma}^{\dagger}c_{j\sigma}\rangle &=\frac{\langle \phi_{\rm SL}| c_{i\sigma}^{\dagger}c_{j\sigma} | \phi_{\rm SL}\rangle}{\langle \phi_{\rm SL}|\phi_{\rm SL}\rangle } \\ &=\sum_{n} \Phi_{in\sigma} \Phi_{jn\sigma}. \end{aligned}\end{split}\]

が導かれます。

次に, \(|\phi_{\rm SL}\rangle\) を変形し, \(f_{ij}\) と \(\Phi_{in\sigma}\) の間に成り立つ関係式を示します。 \(\psi^{\dagger}_{n\sigma}\) の交換関係を用いると, \(|\phi_{\rm SL}\rangle\) は

\[\begin{aligned} |\phi_{\rm SL}\rangle \propto \prod_{n=1}^{N_{e}/2} \Big(\psi_{n\uparrow}^{\dagger}\psi_{\mu(n)\downarrow}^{\dagger}\Big)|0\rangle, \end{aligned}\]

と書き換えられます。ここで, \(\mu(n)\) は \(n= 1, 2, \cdots, N_{e}/2\) の置換を表します。 議論を簡単にするため, 同一のペア \(n=\mu(n)\) となる場合を考えましょう。 このとき, \(K_{n}^{\dagger}=\psi_{n\uparrow}^{\dagger}\psi_{n\downarrow}^{\dagger}\) として, \(K_{n}^{\dagger}K_{m}^{\dagger}=K_{m}^{\dagger}K_{n}^{\dagger}\) の関係を用いることで,

\[\begin{split}\begin{aligned} |\phi_{\rm SL}\rangle &\propto \prod_{n=1}^{N_{e}/2}\Big(\psi_{n\uparrow}^{\dagger}\psi_{n\downarrow}^{\dagger}\Big)|0\rangle =\prod_{n=1}^{N_{e}/2} K_{n}^{\dagger}|0\rangle \\ &\propto\Big(\sum_{n=1}^{\frac{N_{e}}{2}}K_{n}^{\dagger}\Big)^{\frac{N_{e}}{2}} |0\rangle =\Big(\sum_{i,j=1}^{N_{s}}\Big[\sum_{n=1}^{\frac{N_{e}}{2}}\Phi_{in\uparrow}\Phi_{jn\downarrow}\Big] c_{i\uparrow}^{\dagger}c_{j\downarrow}^{\dagger}\Big)^{N_e/2}|0\rangle, \end{aligned}\end{split}\]

の関係が得られます。これより \(f_{ij}\) は単一スレーター行列式の係数により

\[\begin{aligned} f_{ij}=\sum_{n=1}^{\frac{N_{e}}{2}}\Phi_{in\uparrow}\Phi_{jn\downarrow}. \end{aligned}\]

として表されることが分かります。なお, この形式は単一スレーター行列式で与えられる \(f_{ij}\) の表式の一つであり, 実際にはペアを組む自由度(どの \(\mu(n)\) を選ぶか)およびゲージの自由度 (すなわち \(\Phi_{in\sigma}\) の符号の自由度)に依存します。 この自由度の多さが \(f_{ij}\) の冗長性につながっています。

\(F_{IJ}\) と \(\Phi_{In}\) の関係 (スピン平行も含めた場合)

前節で考察したパフィアン-スレーター波動関数と単一スレーター波動関数の間の関係は、 同種スピンのペアリングも考えた場合に拡張することができます。 パフィアン-スレーター波動関数とスレーター波動関数をそれぞれ

\[\begin{split}\begin{aligned} |\phi_{\rm Pf}\rangle&=\Big(\sum_{I,J=1}^{2N_{s}}F_{IJ}c_{I}^{\dagger}c_{J}^{\dagger}\Big)^{N_{\rm e}/2}|0\rangle, \\ |\phi_{\rm SL}\rangle&=\Big(\prod_{n=1}^{N_{e}}\psi_{n}^{\dagger}\Big)|0\rangle,~~\psi_{n}^{\dagger}=\sum_{I=1}^{2N_{s}}\Phi_{In}c^{\dagger}_{I}.\end{aligned}\end{split}\]

と定義します。ここで \(I,J\) はスピン自由度も含めたサイトのインデックスです。 スピン反平行の場合とほぼ同様の議論を用いることで,

\[\begin{aligned} F_{IJ}=\sum_{n=1}^{\frac{N_{e}}{2}}\Big(\Phi_{I,2n-1}\Phi_{J,2n}-\Phi_{J,2n-1}\Phi_{I,2n}\Big). \end{aligned}\]

の関係を示すことができます。これはスピン反平行のペアリングにもそのまま適用できるので, mVMC ver1.0以降ではこの表式を使用しています。

\(f_{ij}\) の特異値分解

行列 \(F\), \(\Phi_{\uparrow}\), \(\Phi_{\downarrow}\), \(\Sigma\) を

\[\begin{split}\begin{aligned} &(F)_{ij}=f_{ij},~~~ (\Phi_{\uparrow})_{in}=\Phi_{in\uparrow},~~~ (\Phi_{\downarrow})_{in}=\Phi_{in\downarrow}, \\ &\Sigma={\rm diag}[\underbrace{1,\cdots,1}_{N_e/2},0,0,0],\end{aligned}\end{split}\]

として定義します。前節のように \(f_{ij}\) (すなわち:math:F)が単一スレーター行列と関係づけられて いるとき、 \(F\) の特異値分解は

\[\begin{aligned} F=\Phi_{\uparrow}\Sigma\Phi_{\downarrow}^{t}.\end{aligned}\]

となることを示すことができます。 この結果は、一般に \(F\) を特異値分解したとき、非ゼロの特異値が \(N_{e}/2\) 個存在し, かつ全ての \(F\) の非ゼロの特異値が \(1\) であった場合, \(f_{ij}\) が単一スレーター波動関数を 記述すること(つまり平均場近似解として記述できること)を表しています。 言い換えると, 特異値の非ゼロ成分の数とその値が, シングルスレータ行列式からパフィアンスレーター行列式がどのようにしてずれるのか, という点について定量的な基準を与えることを示しています。

Power-Lanczos法

このセクションでは, Power-Lanczos法での最適化パラメータ \(\alpha\) の決定方法について述べます。 また, ここではシングルステップのLanczos法を適用した後の物理量の計算についても説明します。

\(\alpha\) の決定

最初に, 変分モンテカルロ法のサンプリングに関して簡単に説明します。 物理量 \(\hat{A}\) は以下の手順で計算されます:

\[\begin{split}\begin{aligned} &\langle \hat{A}\rangle = \frac{\langle \phi| \hat{A}|\phi \rangle}{\langle \phi| \phi \rangle} = \sum_{x} \rho(x) F(x, {\hat{A}}),\\ & \rho(x)=\frac{|\langle \phi|x\rangle|^2}{\langle \phi | \phi \rangle}, ~~~~F(x, {\hat{A}}) = \frac{\langle x| \hat{A}|\phi \rangle}{\langle x| \phi \rangle}.\end{aligned}\end{split}\]

演算子の積 \(\hat{A}\hat{B}\) を計算する場合には, 以下の二通りの方法があります。

\[\begin{split}\begin{aligned} &\langle \hat{A} \hat{B}\rangle = \sum_{x} \rho(x) F(x, {\hat{A}\hat{B}}),\\ &\langle \hat{A} \hat{B}\rangle = \sum_{x} \rho(x) F^{\dagger}(x, {\hat{A})F(x, \hat{B}}).\end{aligned}\end{split}\]

後述するように, 後者の表記の方が数値的には安定します。 例えば, エネルギーの期待値の分散 \(\sigma^2=\langle (\hat{H}-\langle \hat{H}\rangle)^2\rangle\) を考えてみると, 分散は以下の2通りの方法で計算できます。

\[\begin{split}\begin{aligned} \sigma^2 &=\sum_{x} \rho(x) F(x, (\hat{H}-\langle \hat{H}\rangle)^2) = \sum_{x} \rho(x) F(x, \hat{H}^2) - \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 ,\\ \sigma^2 &=\sum_{x} \rho(x) F^{\dagger}(x, \hat{H}-\langle \hat{H}\rangle)F(x, \hat{H}-\langle \hat{H}\rangle) \nonumber \\ &= \sum_{x} \rho(x) F^{\dagger}(x, \hat{H}) F(x, \hat{H})- \left[ \sum_{x} \rho(x) F(x, \hat{H})\right]^2 \end{aligned}\end{split}\]

この定義から, 後者の方法では常に正の値となることが保証されているのに対して, 前者の方法では分散が正の値になることが必ずしも保証されないことが分かります。次に, シングルステップでのpower-Lanczos波動関数 \(|\phi\rangle =(1+\alpha \hat{H}) |\psi \rangle\) に対するエネルギーの期待値とその分散を考えます。エネルギーは以下の式で計算されます：

\[\begin{aligned} E_{LS}(\alpha) =\frac{\langle \phi| \hat{H} |\phi\rangle}{\langle \phi|\phi\rangle}=\frac{h_1 + \alpha(h_{2(20)} + h_{2(11)}) + \alpha^2 h_{3(12)}}{1 + 2\alpha h_1 + \alpha^2 h_{2(11)}}, \end{aligned}\]

ここで, \(h_1\), \(h_{2(11)},~h_{2(20)},\) and \(h_{3(12)}\) を以下のように定義しました：

\[\begin{split}\begin{aligned} &h_1 =\sum_{x} \rho(x) F^{\dagger}(x, \hat{H}),\\ &h_{2(11)}=\sum_{x} \rho(x) F^{\dagger}(x, \hat{H}) F(x, \hat{H}),\\ &h_{2(20)}=\sum_{x} \rho(x) F^{\dagger}(x, \hat{H}^2),\\ &h_{3(12)}=\sum_{x} \rho(x) F^{\dagger}(x, \hat{H})F(x, \hat{H}^2).\end{aligned}\end{split}\]

\(\frac{\partial E_{LS}(\alpha)}{\partial \alpha}=0\) の条件から \(\alpha\) の二次方程式が導出され, それを解くことで \(\alpha\) が決定されます。 分散に関しても同様の手法で計算することが可能です。

物理量の計算

最適化パラメータ \(\alpha\) を用いることで, 演算子 \(\hat{A}\) の期待値を以下の式から計算することが出来ます：

\[\begin{aligned} A_{LS}(\alpha) =\frac{\langle \phi| \hat{A} |\phi\rangle}{\langle \phi|\phi\rangle}=\frac{A_0 + \alpha(A_{1(10)} + A_{1(01)}) + \alpha^2 A_{2(11)}}{1 + 2\alpha h_1 + \alpha^2 h_{2(11)}}, \end{aligned}\]

ここで, \(A_0\), \(A_{1(10)},~A_{1(01)},\) and \(A_{2(11)}\) は

\[\begin{split}\begin{aligned} &A_0 =\sum_{x} \rho(x) F(x, \hat{A}),\\ &A_{1(10)}=\sum_{x} \rho(x) F^{\dagger}(x, \hat{H}) F(x, \hat{A}),\\ &A_{1(01)}=\sum_{x} \rho(x) F(x, \hat{A}\hat{H}),\\ &A_{2(11)}=\sum_{x} \rho(x) F^{\dagger}(x, \hat{H})F(x, \hat{A}\hat{H}).\end{aligned}\end{split}\]

として定義される変数を表します。プログラムでは, この表式に基づき一体グリーン関数および二体グリーン関数の計算を行っています。